Cite this Article
Wang Shuai, Long Hai-Feng, Bi Zhen-Hua, Jiang Wei, Xu Xiang, Wang You-Nian. One-dimensional hybrid simulation of the electrical asymmetry effect caused by the fourth-order harmonic in dual-frequency capacitively coupled plasma. Chinese Physics B, 2016, 25(11): 115202  
Permissions
One-dimensional hybrid simulation of the electrical asymmetry effect caused by the fourth-order harmonic in dual-frequency capacitively coupled plasma
Wang Shuai1, †, , Long Hai-Feng1, Bi Zhen-Hua2, Jiang Wei3, Xu Xiang4, Wang You-Nian4
Department of Physics, College of Sciences, Northeastern University, Shenyang 110891, China
Liaoning Key Lab of Optoelectronic Films & Materials, School of Physics and Materials Engineering, Dalian Nationalities University, Dalian 116600, China
School of Physics, Huazhong University of Science and Technology, Wuhan 430074, China
School of Physics and Optoelectronic Technology, Dalian University of Technology, Dalian 116024, China

 

† Corresponding author. E-mail: s_wang@mail.neu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11305032, 11305028, 11375163, and 11275039) and the Scientific Foundation of Ministry of Education of China (Grant No. N130405008).

Abstract
Abstract

A one-dimensional hybrid model was developed to study the electrical asymmetry effect (EAE) caused by the fourth-order harmonic in a dual-frequency capacitively coupled Ar plasma. The self-bias voltage caused by the fourth-order frequency changes periodically with the phase angle, and the cycle of self-bias with the phase angle is π/2, which is half of that in the second-order case. The influence of the phase angle between the fundamental and its fourth-order frequency on the ion density profiles and the ion energy distributions (IEDs) were studied. Both the ion density profile and the IEDs can be controlled by the phase angle, which provides a convenient way to adjust the sheath characters without changing the main discharge parameters.

PACS: 52.65.–y;52.80.Pi
Keyword:capacitively coupled plasmas;electrical asymmetry effect;Hybrid model;ion energy distributions
1. Introduction

Capacitively coupled plasmas (CCP) are of vital importance in materials processing such as etching and deposition in the semiconductor industry.[1,2] In the industrial processes, throughput is determined by the ion flux, etching and deposition processes are controlled by the ion energy, therefore, the independent control of the ion flux and ion energy is necessary for these applications.[3] However, separate control of ion flux and ion energy is difficult in the single frequency discharge. In order to overcome this problem, the dual-frequency capacitively coupled plasma (DF-CCP)[4] was proposed to control the ion flux and energy separately. The traditional DF-CCP[4–9] usually has two frequency sources which are applied to one or two electrodes,[4,6,7,9–17] wherein the low frequency source controls the ion energy and angular distributions, while the high frequency source controls the plasma density and ion flux. Although the traditional DF-CCP has been used for a long time, this mode has been investigated and confirmed that it might limit the separate control of the ion flux and energy recently.[18,19]

A novel idea based on a newly discovered electrical asymmetry effect (EAE) was introduced to achieve high-degree separate control of ion flux and energy in dual-frequency capacitively coupled plasmas.[20,21] The two sheaths close to the surface of electrode are electrically asymmetrical by applying a fundamental and a second harmonic to the discharge, even the discharge itself is geometrically symmetric. A DC self-bias develops both in geometrically symmetric and geometrically asymmetry discharges.[20] By changing the phase angle, one can modulate the ion energy distributions according to the linear relationship between DC self-bias voltage and the phase angle Θ. The applying of the second harmonic may lead to a best control of both the ion energy and the ion flux. Recently, the electrical asymmetry effect in DF-CCP discharge has been studied both numerically and experimentally[22,27] and has been analyzed in detail through the analytical model, fluid and particle in cell (PIC) simulations. In particular, Donkó et al.[22] studied EAE operated in the 13.56 MHz and 27.12 MHz in the geometric symmetry discharge. The result confirmed the theoretical prediction of Heil et al.[20] and verified that the separate control of ion flux and energy can be realized through the EAE in an ideal manner. Later, Schulze et al.[23] tested this result in the EAE experiment for the first time. Electron heating in electropositive capacitively coupled plasma has been studied by using experiments, simulations and models.[28–41] For example, Jiang et al.[42,43] developed a PIC/MC model and studied the heating mechanisms in detail. Lisovskiy et al. studied the three discharge modes and the mode transitions in low pressure nitrogen[44–46]. Schulze et al. investigated electrically asymmetric DF-CCP discharges in CF4 by particle-in-cell simulation at different pressures and values of Θ.[47]

In previous studies, a common choice for the two frequencies was fundamental frequency f with its second harmonic 2 f. To achieve a separate control of the ion energy and the ion flux, a large disparity of the two frequencies is needed. Therefore, a higher-order harmonic would be needed in order to achieve both the EAE and the separate control of the ion flux and the energy. In this paper, we choose the fundamental frequency f and its four-order harmonic 4 f. Firstly, an analytic discussion is presented based on the physical model by Heil et al.[20] By using the four-order frequency, the self-bias varies with the phase angle more rapidly, the period abridges from π to π/2. In this case, more data are needed to study the relationship between the phase angle and the self-bias, therefore, in this paper, we study the EAE of Ar plasmas using a one-dimensional fluid/MC hybrid model which can provide an accurate simulation result more efficiently. The concept of a hybrid model[48–51] was proposed in the 1990s. The model consists of two parts which are a fluid model and a Monte Carlo model respectively. The study of the ion energy distributions (IEDs) is quite necessary since the ion flux plays an important role during the surface etching procedure. In particular, we focus on the IEDs striking the electrode, and make a comparison of DC self-bias and IEDs with different phase angles between the two frequencies.

The main contents of this paper are as follows. The analytic physics model is presented in Section 2. The one-dimensional fluid/MC hybrid model for the discharge, the numerical results, and the discussion are given in Section 3, and a brief summary is in Section 4.

2. The analytical calculation of the self-bias

The applied RF voltage is

where Θ is the phase angle between the fundamental and its fourth harmonic, ω = 2πf. The total voltage drop cross the discharge can be written as

By normalizing the voltages to V0 and defining Vdc = V0 η, ϕ = ωt, equation (1) becomes

where ϕ̃ represents the oscillating part, and the total voltage ϕ becomes

then we take the derivative with φ. In order to obtain the extreme value of Eq. (4), the result should be set to 0, that is

As equation (5) is too complicated to solve analytically, it has not been solved directly. Instead, we show the relationship between φm and Θ satisfying Eq. (5) in Fig. 1. It can be seen that all the solutions of φm can be approximately set as

By substituting the approximate value of φm into Eq. (3), the oscillating part of the voltage drop can be written as

The pictures of the eight functions are plotted in Fig. 2, from which we can construct the maximum voltage ϕ̃m1 and the minimum voltage ϕ̃m2 of the oscillating part. The maximum voltage is

and the minimum voltage is

The self-bias could be calculated by η = −(ϕ̃m1 + ϕ̃m2), where we use Eq. (22) in Ref. [20] and make the assumption that the symmetry parameter ɛ equals to 1. The expression of self-bias varies in eight sections from −π to π, which is

Equation (17) is plotted in Fig. 3, which shows that the self-bias η varies linearly with the phase angle Θ within each of the eight sections, as the Θ increases, the self-bias η changes periodically with a period π/2. This shows that a self-bias η can also be induced by changing the phase angle of the fourth-order harmonic 4 f. Compared with the second harmonic, the extreme value of η is relatively smaller and the period changes from π to π/2.

Fig. 1. Solutions of Eq. (5), the RF phases corresponding to the max voltage drop as a function of the phase angle between the first and the fourth harmonics.
Fig. 2. The picture of Eqs. (7)–(14) which shows how the extreme voltages are constructed from the above eight functions.
Fig. 3. The value of η calculated from Eq. (8), using the assumption that the symmetry parameter equals to 1.
3. Numerical simulation of the EAE
3.1. Hybrid model

In this section, the hybrid model containing a fluid model and a Monte Carlo model is introduced. The fluid model mainly calculates the physical properties of the bulk plasma in the transport process and the Monte Carlo method is used to compute the ion energy distributions (IEDs) after the collision occurs. This method can not only improve the computational efficiency, but also help to obtain accurate results.

The schematic diagram of the capacitively coupled plasma reactor is shown in Fig. 4, in which two electrodes separate from each other with a distance of 2 cm. A sinusoidal RF voltage is used on the upper electrode at x = 2 cm, and the self-bias on the lower electrode at x = 0 cm is calculated to balance the electron and ion flux.

Fig. 4. Schematic diagram of the DC/RF combined driven capacitively coupled plasma reactor.

The fluid model consists of continuity and momentum equations of electrons and ions, the energy balance equation of electrons, and Poisson’s equation. In the discharge process, the collisions between ions and surrounding neutral gases are frequent, therefore, the ion temperature is basically consistent with the neutral gas, and it is unnecessary to solve the energy equation of ions. The drift-diffusion approximation is valid for describing the electron flux and ion flux in the fluid model since the pressure in discharges is not too low. In this mode, the ionization occurs mainly in the region of the bulk plasma, and the ionization does not play an essential role. Therefore, we ignore the secondary electron emission, and the electron flux and the electron energy flux to the electrodes are assumed as follows:

where uth is the electron thermal velocity and Λ is the electron reflection coefficient. Besides, the ion density and velocity gradients are equal to 0 at the electrodes, that is, ∂ni/∂x = 0 and ∂ui/∂x = 0.[52]

Furthermore, we can describe the behaviors of ions and electrons and the electric potential by the following equations:

Here, the letters e and i represent electron and positive argon ions, ni and Γi are the density and flux of argon ions, respectively; Si and Se are source terms; μi and Di are the mobility and diffusion coefficients of ions; ne, Γe, and Te are the density, flux, and temperature of electrons, respectively; We represents the energy loss due to inelastic collisions between electrons and neutrals;[53] V is the electric potential; E = −∂V/∂x is the electric field; and ɛ0 is the vacuum permittivity. The initial distributions of the plasma density and initial temperature are considered to be homogeneous with n0 = 109 cm−3 and Te = 2 eV. The MC simulation begins on the condition that the system achieves a steady state, typically after 1000 RF cycles.

After obtaining the plasma density, the electric field, and other macroscopic quantities in the fluid module, we can trace the trajectories of test particles in the discharge process. The grid step of space in the MC model is Δx = 0.5 × 10−4 m, the time step is one thousandth of the RF period which equals to Δt = 0.25 × 10−10 s, and the velocity is normalized by the factor of Δx/Δt = 0.5 × 10−4/0.25 × 10−10 m · s−1. The motion equation is calculated by Newton’s equation and the collision process is simulated by the Monte Carlo method. The energy and angular distributions of each particle will be recorded when the test particle reaches the designated position. According to the statistical result of the distributions by using a large number of test particles, we can accurately calculate the energy and angle distributions of the electrons, ions, and groups anywhere we want to know. In this paper, for the argon plasma, we assume that the neutral species have a uniform distribution in space with a Maxwellian velocity distribution, and the temperature of neutral gas is the same as that in the fluid model.

For the argon plasma, the collisions of electrons are elastic scattering (e + Ar → e + Ar), excitation (e + Ar → e+Ar*), and ionization (e + Ar → 2e + Ar+). For the ion–neutral collisions, the elastic scattering (Ar+ + Ar → Ar+ +Ar) and the charge exchange collision (Ar+ + Ar → Ar+ Ar+) are included in the present model.

Besides, we use the null-collision method to reduce the computational time. The following equation shows that the collision probabilities for all the particles are independent of the energy and position

where νmax = max(ng(x))max (σT(ɛ)ν), ng is the density of the target species, and σT is the total cross section for all types of collisions.[54–56] Since the electron energy distribution will affect the transport process in the fluid model, we need to bring the results calculated in the MC model back to the fluid model after a period of time to ensure the accuracy of the calculation.

3.2. Numerical results

The initial value of the plasma potential is set to 0 V, the discharge pressure is set to P = 50 mTorr, and the voltage applied to the lower electrode is

where V0 = 200 V, f = 10 MHz, VDC is the self-bias which is calculated by the hybrid model in order to balance the electron and ion flux toward the upper and the lower electrodes.

Figure 5 shows the temporally averaged ion density profiles with different phase angles. The density profiles are asymmetric due to the effect of the self-bias. When the phase angle changes from 10° to 50°, the self-bias on the lower electrode increases from negative to positive, leading to the decrease of the sheath length near the lower electrode and the increase of the density near the electrode. Therefore, the sheath structure and the ion density near the electrode can be adjusted by changing the phase angle to fit different requirements.

Figure 6 shows the influence of the phase angle Θ on the self-bias VDC. The self-bias changes periodically with the phase angle, the period is π/2, VDC increases linearly with the phase angle Θ within the first half period. By further increasing the phase angle within the second half period, VDC decreases linearly with Θ. Comparing with the self-bias induced by the second-order harmonic, which is presented in Fig. 7 with the same V0, VDC induced by the fourth-order harmonic is relatively small, the max value in Fig. 6 is 30 V less than the maximum 108 V in Fig. 7. This means that the influence on the self-bias caused by a higher-order harmonic is relatively smaller. However, the higher-order harmonic is a good choice to obtain both better separate control of the ion flux and the ion energy, and precise control of the self-bias by adjusting the phase angle Θ.

Fig. 5. The temporally averaged ion density profiles. The phase angles between the first and the fourth harmonics are Θ = 10°,30°,50°.
Fig. 6. The self-bias VDC on the electrode with different phase angles Θ caused by the fourth-order harmonic with V0 = 400 V.
Fig. 7. The self-bias VDC on the electrode with different phase angles Θ caused by the second-order harmonic with V0 = 400 V.

A comparison of the self-bias between the analytic method result and the hybrid model result is presented in Fig. 8. The self-bias calculated from both models changes linearly with the phase angle Θ, and the self-bias from the hybrid model has a phase shift about 10° comparing with that from the analytical model. The maximum self-bias from the analytical model is a little lower than that from the hybrid model, because we keep the symmetry parameter in the analytical model as a constant. However, the changing cycle and the slope of self-bias with the phase angle from these two models are in good agreement.

Fig. 8. The comparison of the self-bias calculated by the analytic method and the hybrid model with V0 = 400 V.

Figure 9 shows the ion energy distributions at the electrode with different phase angles Θ. Like the self-bias VDC, the ion energy distributions also vary periodically with Θ, and within each half of the period, the maximum energy changes linearly with the phase angle. When the phase angle Θ is adjusted from 50° to 100°, the maximum ion energy increases from 150 eV to 180 eV, which proves that the IEDs can be adjusted by changing the phase angle between the first- and fourth-order harmonics.

Fig. 9. The ion energy distributions on the lower electrode with phase angle Θ changing from 0 to 90°.
4. Conclusion

In the present paper, an analytical model and a one-dimensional hybrid model have been used to investigate the electrically asymmetric effect in a capacitively coupled argon discharge driven by the basic 10 MHz and its fourth-order frequency 40 MHz. The self-bias which balances the electron and ion flux toward the two electrodes changes with different phase angles between the basic and the fourth-order frequency. Comparing with the second-order frequency 20 MHz, the cycle of self-bias with the phase angle changing from π to π/2 (half of the 10 MHz–20 MHz case), the self-bias voltage varying linearly within each of the π/4 regions, and the influence on the self-bias caused by the fourth-order frequency is relatively small. However, the density profile and the IEDs on the electrode can also be controlled by the phase angle Θ, the max energy, and the high energy peaks in IEDs change simultaneously with the self-bias. This means that the profile of the density and the ion energy on the electrode can also be controlled by fourth-order frequency, which provides an opportunity to adjust the discharge profile in an easy and convenient way.

Reference
1Lieberman M ALichtenberg A J1994Principles of Plasma Discharges and Material ProcessingNew YorkWeley
2Kim H CLee J KShon J W 2003 Phys. Plasmas 10 4545
3Lieberman M ALichtenberg A J2005Principles of Plasma Discharges and Materials Processing2nd edn.New JerseyWiley Interscience
4Goto H HLowe H DOhmi T 1992 J. Vac. Sci. Technol. 10 3048
5Xu D SZou SXin YSu X DWang X S 2014 Chin. Phys. 23 065201
6Wakayama GNanbu K 2003 IEEE Trans. Plasma Sci. 31 638
7Boyle P CEllingboe A RTurner M M 2004 Plasma Sources Sci. Technol. 13 493
8Turner M MChabert P 2006 Phys. Rev. Lett. 96 205001
9Donkó ZPetrovic L Z 2007 J. Phys. Conf. Ser. 86 012011
10Schulze JDonkó ZHeil B GLuggenhölscher DMussenbrock TBrinkmann R PCzarnetzki U 2008 J. Phys. D: Appl. Phys. 41 105214
11Boyle P CEllingboe A RTurner M M 2004 J. Phys. D: Appl. Phys. 37 697
12Denda TMiyoshi YKomukai YGoto TPetrovic Z LMakabe T 2004 J. Appl. Phys. 95 870
13Gans TSchulze JO’Connell DCzarnetzki UFaulkner REllingboe A RTurner M M 2006 Appl. Phys. Lett. 89 261502
14Schulze JGans TO’Connell DCzarnetzki UEllingboe A RTurner M M 2007 J. Phys. D: Appl. Phys. 40 7008
15Schulze JDonk’o ZLuggenölscher DCzarnetzki U2008The 19th European Sectional Conference on Atomic and Molecular Physics of Ionized GasesJuly 15–19, 2008Granada, Spain206
16Semmler EAwakowicz Pvon Keudell A 2007 Plasma Sources Sci. Technol. 16 839
17Salabas ABrinkmann R P 2005 Plasma Sources Sci. Technol. 14 S53
18Lee J KManuilenko O VBabaeva N YKim H CShon J W 2005 Plasma Sources Sci. Technol. 14 89
19Kawamura ELieberman M ALichtenberg A J 2006 Phys. Plasmas 13 053506
20Heil B GCzarnetzki UBrinkmann R PMussenbrock T 2008 J. Phys. D: Appl. Phys. 41 165202
21Heil B GSchulze JMussenbrock TBrinkmann R PCzarnetzki U 2008 IEEE Trans. Plasma Sci. 36 1404
22Donkó ZSchulze JHeil B GCzarnetzki U 2009 J. Phys. D: Appl. Phys. 42 025205
23Schulze JSchungel ECzarnetzki U 2009 J. Phys. D: Appl. Phys. 42 092005
24Schulze JSchungel ECzarnetzki UDonkó Z 2009 J. Appl. Phys. 106 063307
25Donkó ZSchulze JCzarnetzki ULuggenholscher D 2009 Appl. Phys. Lett. 94 131501
26Schulze JDonkó ZLuggenholscher DCzarnetzki U 2009 Plasma Sources Sci. Technol. 18 034011
27Schulze JSchungel EDonkó ZCzarnetzki U 2010 Plasma Sources Sci. Technol. 19 045028
28Lieberman M A 1988 IEEE Trans. Plasma Sci. 16 638
29Lieberman M AGodyak V A 1998 IEEE Trans. Plasma Sci. 26 955
30Surendra MGraves D B 1991 Phys. Rev. Lett. 66 1469
31Turner M M 1995 Phys. Rev. Lett. 75 1312
32Gozadinos GTurner M MVender D 2001 Phys. Rev. Lett. 87 135004
33Kaganovich I D 2002 Phys. Rev. Lett. 89 265006
34Kaganovich I DPolomarov O VTheodosioue C E 2006 IEEE Trans. Plasma Sci. 34 696
35Gans TSchulz-von der Gathen VDöbele H F 2004 Europhys. Lett. 66 232
36Mussenbrock TBrinkmann R P 2006 Appl. Phys. Lett. 88 151503
37Mussenbrock TBrinkmann R PLieberman M ALichtenberg A JKawamura E 2008 Phys. Rev. Lett. 101 085004
38Czarnetzki UMussenbrock TBrinkmann R P 2006 Phys. Plasmas 13 123503
39Ziegler DMussenbrock TBrinkmann R P 2008 Plasma Sources Sci. Technol. 17 045011
40Schulze JHeil B GLuggenhölscher DBrinkmann R PCzarnetzki U 2008 J. Phys. D: Appl. Phys. 41 195212
41Turner M M 2009 J. Phys. D: Appl. Phys. 42 194008
42Jiang WXu XDai Z LWang Y N 2008 Phys. Plasmas 15 033502
43Zhang Q ZJiang WZhao S XWang Y N 2010 J. Vac. Sci. Technol. 28 287
44Lisovskiy V AYegorenkov V D 1994 J. Phys. D: Appl. Phys. 27 2340
45Lisovskiy V AKharchenko N DYegorenkov V D 2008 J. Phys. D: Appl. Phys. 41 125207
46Lisovskiy V AKharchenko N DYegorenkov V D 2010 J. Phys. D: Appl. Phys. 43 425202
47Schulze JDerzsi ADonk Z 2011 Plasma Sources Sci. Technol. 20 045008
48Sommerer T JKushner M J 1992 J. Appl. Phys. 71 1654
49Ventzek P L GSommerer T JHoekstra R JKushner M J 1993 Appl. Phys. Lett. 63 605
50Kratzer MBrinkmann R PSabischa WSchmidt H J 2001 J. Appl. Phys. 90 2169
51Sato NTagashira H 1991 IEEE Trans. Plasma Sci. 19 102
52Nitschke T EGraves D B 1994 J. Appl. Phys. 76 5646
53Bukowski J DGraves D B 1996 J. Appl. Phys. 80 2614
54Vahedi VSurendra M 1995 Comput. Phys. Commun. 87 179
55Vahedi VBirdsall C KLieberman M A 1993 Phys. Fluids 5 2719
56Wang H YJiang WSun PKong L B 2014 Chin. Phys. 23 035204