Simulations of standing wave effect, stop band effect,and skin effect in large-area very high frequency symmetric capacitive discharges

Jiankai LIU (刘建凯) , Yuru ZHANG (张钰如), Kai ZHAO (赵凯),Deqi WEN (温德奇) and Younian WANG (王友年)

1 Key Laboratory of Materials Modification by Laser, Ion, and Electron Beams (Ministry of Education),School of Physics, Dalian University of Technology, Dalian 116024, People’s Republic of China

2 Department of Electrical and Computer Engineering,Michigan State University,East Lansing,MI 48824,United States of America

Abstract In this paper, Maxwell equations are coupled with a radially localized global model and an analytical sheath model to investigate the electromagnetic effects under various frequencies and electron powers in large-area very high frequency symmetric capacitive argon discharges.Simulation results indicate that both the vacuum wavelength and the sheath width decrease with frequency, leading to the reduced surface wavelength.As a result, the standing wave effect becomes pronounced, causing the fact that the radial profiles of the electron density, radio frequency voltage,and sheath width shift from uniform over center-high to multiple-node.When the frequency is close to or higher than the series resonance frequency,the surface waves cannot propagate to the radial center because of the significant radial damping.Due to the lack of power deposition near the radial center, the electron density is nearly zero there, i.e.the stop band effect.As power increases, the higher electron density leads to the decrease of the skin depth.Therefore, the importance of the skin effect gradually exceeds that of the standing wave effect,giving rise to the transition from the center-high to edge-high electron density profiles.The method proposed in this work could help to predict the plasma distribution under different discharge conditions in a few minutes, which is of significant importance in optimizing the plasma processing.

Keywords: electromagnetic effects, hybrid model, plasma uniformity

1.Introduction

Capacitively coupled plasmas (CCPs) are widely used for plasma enhanced chemical vapor deposition and dry etching in the semiconductor industry [1, 2].In these plasma processes, large-area CCPs become popular because of the current trend of increasing wafer size.Besides, very high frequency (VHF) sources have attracted growing interest due to the higher ion flux and lower ion bombarding energy to the substrate [3-14].However, in large-area VHF CCPs, the electromagnetic (EM) effects negatively affect the plasma radial uniformity and thus limit the plasma processing[15-19].

In the past two decades, several theoretical models, i.e.uniform slab model [15, 16, 20] and transmission line model[21-27], have been adopted to investigate the propagation of surface waves and the EM effects.Lieberman et al [15]employed a uniform slab model and solved the Maxwell equations analytically for a geometrically symmetric CCP discharge.They first addressed the standing wave effect and the skin effect, and pointed out the conditions under which these EM effects could be neglected.Sansonnens et al [16]developed a similar uniform slab model with the same sheath width surrounded for a geometrically asymmetric CCP discharge.They found that only the first even (axially symmetric) and the first odd (axially antisymmetric) modes of surface waves, which were responsible for the standing wave effect and the telegraph effect, could propagate inside the chamber.After that,Lieberman et al[20]used a uniform slab model with different sheath widths for a geometrically asymmetric CCP discharge, and they presented a series of system resonances associated with the symmetric and antisymmetric mode surface waves under various frequencies.Moreover, a transmission line model coupled with a global discharge model was developed by Chabert et al [21-24] to study the standing wave effect and the skin effect self-consistently, and they concluded that the CCP discharge was mainly sustained by the capacitive (E) field at lower radio frequency(rf)voltages and by the inductive(H)field at higher rf voltages, i.e.the E-H mode transition.To consider the influence of nonlinearly excited harmonics on the standing wave effect, Lieberman et al [25] developed a nonlinear transmission line model with uniform plasma density and a single homogeneous sheath, revealing that the series resonance-enhanced harmonics induced the standing wave effect and made the center-high profile of the electron power density more pronounced.Subsequently, more versatile nonlinear transmission line models were developed by Wen et al[26, 27] to examine the series resonance effect and the corresponding wave effects.

In addition, numerous self-consistent numerical models have also been developed to study the EM effects[17,28-50]and the modulation of the plasma radial uniformity[28, 30, 33, 34, 40, 43-48].Lee et al [17] combined a bulk fluid model with an analytical sheath model and the Maxwell equations,and they observed the stop band effect for the first time in the frequency range between the series resonance frequency and the electron plasma frequency, and this effect has been captured experimentally by Liu et al [19, 51] and Han et al [52] later.By a fast fluid-analytical model, Kawamura et al [28-30] presented that the higher harmonics contributed significantly to the center-high plasma profile, and this could be suppressed by adding an additional low frequency power or a dielectric layer.Recently,they proved that when the frequency was above or near an antisymmetric spatial resonance frequency,the symmetric and antisymmetric modes could coexist in a symmetric discharge,leading to the non-symmetric distribution of the plasma parameters about the midplane [31].Rauf and Bera et al [32-37] investigated the effects of external parameters, such as the source power,inter-electrode gap, plasma electronegativity and inhomogeneous magnetic field on the VHF CCPs by a two-dimensional fluid model.Moreover,Chen et al[38,39]developed a threedimensional fluid model and solved the Maxwell equations by the finite difference time domain technique for a large-area rectangle chamber, focusing on the influence of segmenting powered electrode on the plasma uniformity.Xu and Zhang et al [40-44] reported the phase shift effect on the plasma distribution, and recently Zhang et al studied the plasma characteristics in discharges sustained by multiple consecutive harmonics in the VHF regime [45].Yang et al[46-48] used the hybrid plasma equipment model to analyze the effects of different parameters on the plasma spatial profile and the ion energy distribution,and they proposed that the plasma uniformity could be improved by graded conductivity electrodes.Upadhyay and Sawada et al [49, 50] solved the fluid equations and the Maxwell equations under a highresolution computational mesh,and their results indicated the importance of higher harmonics on the center-high plasma profile, which has also been observed in [11, 13, 14, 25-27,29, 53, 54].

It is known that the EM effects,which originate from the propagation of surface waves, have a significant influence on the plasma spatial distribution.Vice versa, the propagation characteristics of surface waves are affected by the plasma properties,such as the plasma spatial profile,sheath width and so on [15, 16, 20].However, in the previous research mentioned above, the surface wave propagation was usually studied with the uniform slab model, in which the plasma properties were assumed to be spatially uniform.In the fluid simulations,the emphases were mainly put on the variation of the plasma properties with different discharge parameters, as well as the modulation of the plasma uniformity.The surface wave propagation characteristics were overlooked, and the simulation regions were limited due to the computational burden.

In order to have a deeper understanding of the VHF CCP discharges, a self-consistent model, which includes the interaction between the spatially non-uniform plasma parameters and the propagation of surface waves,is necessary.In this work,a radially localized global model is coupled with an analytical sheath model and the Maxwell equations to selfconsistently investigate the propagation characteristics of surface waves together with the plasma properties under various discharge parameters.Note that we assume a fixed non-uniform axial density profile, and the radial distribution of the plasma density is determined by applying the global model at different radial positions.Moreover,since the series resonance-enhanced harmonics cannot be excited in the geometrically symmetric discharges performed in this study,the Maxwell equations are solved in the frequency domain to reduce the calculation time [27].This strategy allows us to self-consistently predict the plasma uniformity under various discharge conditions in a few minutes.Besides,the transition from the standing wave effect to the stop band effect with frequency,as well as the expansion of the stop band region,is captured.To understand the evolution of the electron density distribution with discharge parameters, which derives from different EM effects, the self-consistently predicted electric field, and thus the power deposition are presented.Besides,the surface wavenumbers, as well as the surface wavelength,radial damping length and skin depth, are calculated for a deeper understanding of the underlying mechanisms behind the EM effects.

Figure 1.Schematic diagram of a geometrically symmetric cylindrical CCP reactor.

This paper is organized as follows.The radially localized EM global model is described in section 2.In section 3, we first investigate the appearance of the standing wave effect,and then we present the transition from the standing wave effect to the stop band effect with increasing frequency.After that, the influence of the electron power on the shifting from the standing wave effect to the skin effect is discussed.Finally, conclusions are given in section 4.

2.Radially localized EM global model

The schematic diagram of a geometrically symmetric cylindrical CCP reactor is shown in figure 1.A large chamber with a radius of R=50 cm and a gap of 2l=5 cm is adopted due to the current trend of increasing wafer size.The half width of the bulk plasma is = −d l s0withs0the sheath width.

2.1.EM equations

The EM equations solved in this study have been described in detail in the previous work[55],so only a brief description is given here.Assumed a transverse magnetic mode, the Maxwell equations in the frequency domain are expressed as

Here,φBis the angular component of the magnetic field,Eris the radial component of the electric field,Ezis the axial component of the electric field,ωis the driving angular frequency,ε0is the vacuum permittivity,μ0is the vacuum permeability, andκpis the plasma relative permittivity.

In the bulk region, the relative permittivity iswhereνmis the electron-neutral momentum transfer frequency andis the electron plasma frequency with the element charge e, the electron massme,and the electron densityne.The electron densitynein the bulk region is assumed to be spatially nonuniform along the axial direction

which is a good approximation at low pressures [1, 55].The radial profile of the electron density(i.e.ne(r, 0))atz=0 cm is determined by the following radially localized global model.

In the sheath region,the relative permittivity varies along the radial directionκp=s0/s(r).This allows the model to use a simulation-convenient fixed sheath thicknesss0,and meanwhile to equivalently take the propagation of surface waves in an actual radial varying sheaths(r) into account[17, 28].Here,s0is determined by the position ofwhereuiis the ion velocity andis the Bohm velocity with the electron temperatureTeand the ion mass M(equation (5.3.5) in [1]).Note thats0influences the electron density and hence the power deposition near the bulk-sheath interface,and the propagation of surface waves is affected by the actual sheath widths(r) which is determined by the following analytical sheath model.

From equations (1)-(3), a Helmholtz equation forφBis obtained

wherek0=ω/cis the vacuum wavenumber,with c the speed of light in vacuum.The boundary conditions forφBare as follows: (1)Bφ=0 at the radial center.(2)en·∇Bφ=0at two electrodes with enthe unit normal vector.(3)at the radial edge, whereIrfis the source current.SubstitutingφBinto equations (1) and (2), the EM fields can be determined, and the capacitive and inductive power densities in the bulk region are given by

2.2.Radially localized global model

In the radially localized global model, the electron temperature is determined by the particle balance equation.We assume that the electrons reach the bulk-sheath interface with the Bohm velocityuB, and then they are lost at the electrodes

HereN=Pg/kBT0is the background gas density, with the gas pressurePg,the gas temperatureT0=293 K,and the Boltzmann constantk.BKizis the ionization rate coefficient,as shown in table 1.is the axial edge to center density ratio, with the ion-neutral mean free pathDue to the fixed axial density profile (see equation (4)), the ratio ofin equation(8)is a constant,soKizand thusTeis independent on the radial position.

Table 1.Electron impact reaction rate coefficients in terms of the electron temperature Te (e V) (section 3.5 in [1]).

In the model, the plasma radial transport is ignored and the local heating is assumed,which is similar to the models in[21, 23, 24].The electron energy balance equation, which is used to calculate the electron density, is given by

Here Ee= Ec+7.2Teis the total electron energy loss,whereT7.2erepresents the electron energy loss when the electrons cross the sheath to the electrode [20].Ec=15.76+12.14Kexc/Kiz+3meTeKelas/MKizis the collisional electron energy loss per electron-ion pair created,Kexcis the average energy loss-weighted excitation rate coefficient for average excitation energy 12.14 eV [1], andKelasis the elastic collision rate coefficient (see table 1).Ptot(r)=Pz(r)+Pr(r)+2Sstoc(r)+2Sohm,sh(r) is the total electron power area density.are the capacitive and inductive power area densities in the bulk region.Sstoc(r) andSohm,sh(r) are the stochastic heating and the Ohmic heating in each sheath, and they are calculated by the analytical sheath model below.

2.3.Analytical sheath model

The analytical collisionless sheath model is used to determine the sheath width and the electron sheath heating(section 11.2 in [1]).The maximum sheath width is given by

2.4.Radially localized EM global model overview

First, the EM fields are determined by solving equations(1)-(3) with the initialne(r,z) ands(r) ,and thenV(r),pz(r,z),andpr(r,z) are obtained.Subsequently,the analytical sheath model is implemented, the sheath widths(r) is updated, and the sheath heatingSohm,sh(r)andSstoc(r)are calculated.Then,the energy balance equation (9) is solved at each radial position forne(r, 0) ,and the axial distributionne(r,z) is determined according to equation (4).Finally, the electron temperatureTeis obtained by solving the particle balance equation (8).In each iteration, the source currentIrfis adjusted to ensure that the calculated total electron power,which is obtained by integratingPtot(r) over the whole electrode, equals to the fixed value.The iteration continues until the energy balance is achieved at each radial position.

3.Results and discussion

3.1.Standing wave effect

First, we study the appearance of the standing wave effect in the frequency range from 13.56 to 100 MHz, with the gas pressure of 50 mTorr and the electron power of 100 W.Figure 2(a) shows the radial profiles of the electron density(i.e.ne(r, 0)) atz=0 cm for various frequencies.It is clear thatne(r, 0) is quite uniform along the radial direction at 13.56 MHz, even near the radial edge, because the particle loss at the side wall and the plasma transport are not considered in the global model.As the frequency increases to 60 MHz,ne(r, 0)shows an obvious center-high profile,due to the standing wave effect.At 80 MHz,ne(r, 0) first decreases along the radial direction, with a density node appearing atr=37.3 cm,and then it increases from the node to the radial edge.As the frequency increases further to 100 MHz, the density node moves inward to =r26.0 cm, and the second density peak appears at =r40.0 cm, indicating that the standing wave effect becomes more pronounced.

Figure 2.Radial profiles of (a) the electron density (i.e. n e(r , 0)) at z =0 cm, (b) the rf voltage (i.e.Vr f(r )) between two electrodes, and(c) the sheath width s ( r) for different frequencies: 13.56, 60, 80, and 100 MHz, at 50 mTorr and 100 W.

To illustrate the standing wave effect, the dispersion relations of surface waves [15, 16] are solved under various frequencies by using the spatially averaged electron densityne,aveand the radially averaged sheath widthsave[29-31] (see appendix A).Since the total electron power is fixed as we mentioned above,ne,ave=8.4 ×109cm−3stays unchanged at four frequencies.At 13.56 MHz,the obtained radial wavenumber iskr=0.0050 −0.000042j withsave=1.15 cm.Note that Re (kr)R=0.25is much smaller than the frist zero of the zeroorder Bessel functionχ01= 2.40,indicating that the surface wavelength (i.e.λ=2π/Re(kr)=1257 cm) is much longer than the chamber radius R.This gives rise to the uniform rf voltage between two electrodesVrf(r) and the uniform sheath widths(r) (fgiures 2(b) and (c)), and consequently the uniform sheath heating (i.e.Psh(r) in fgiure 3(b)).As the frequency increases,Vrf(r) ands(r) show overall downward trend and significant radial non-uniformity, which is consistent with the radial proflies ofne(r, 0) ,and this will be discussed below.

The spatial profiles of the electric fieldE in the bulk region, as well as the radial distributions of the power area density, are presented in figure 3.The contour lines and the arrows in the first column indicate the amplitude and direction ofE, respectively.To reveal the electric field in the bulk region more clearly, we have excluded the sheath regions in all the figures ofE.In the second column of figure 3, the power area density in sheathsPsh(r) is calculated byPsh(r)=2Sstoc(r)+2Sohm,sh(r).Note that although the discharge gap is larger than the electron mean free path at 50 mTorr and 100 W, the stochastic heating plays an important role under this condition[21,56-60].For instance,when the frequency increases from 13.56 to 180 MHz, the stochastic heating increases from 57 to 84 W, while the Ohmic heating decreases from 43 to 16 W.In figure 3(a), ∣E∣is also radially uniform,because of the very long surface wavelength and the radially uniform electron density, as we mentioned above, and this results in the similar distribution ofPz(r)(figure 3(b)).In the axial direction, ∣E ∣shows a significant decrease from the bulk-sheath interface toz=0 cm,which is caused by the higher electron density and thus the higher∣κp∣at the axial center.Besides,Erand the correspondingPr(r)are very weak under this condition,because of the negligible skin effect at the relatively low electron density.As a result,Ptot(r)is mainly determined byPsh(r) andPz(r) ,and the radially uniformPtot(r) leads to the similar distribution ofne(r, 0) in figure 2(a).Note thatPsh(r)represents the capacitive power in the sheaths, indicating that the capacitive power rather than the inductive power dominates the discharge when the electron power is relatively low, i.e.100 W.

Figure 3.Spatial proflies of the electric field E in the bulk plasma(first column), as well as the radial proflies of the capacitive power area density in the bulk Pz ( r ) ,the inductive power area density in the bulk Pr (r ) ,the power area density in sheaths Ps h(r ) ,and the total power area density Pt ot(r )(second column)for different frequencies:13.56 MHz(a)and(b),60 MHz(c)and(d),80 MHz(e)and(f),and 100 MHz (g) and (h), at 50 mTorr and 100 W.

As the frequency increases to 60 MHz, the radial wavenumber iskr=0.0459 −0.0011j withsave=0.20 cm.Note that Re (kr)R=2.30is close toχ01, indicating that the surface wavelength (i.e.λ= 137 cm) is comparable to the chamber size, and the standing wave effect becomes significant.As a result,Vrf(r) ands(r) are characterized by obvious centerhigh profiles in figures 2(b)and (c).However,∣E∣in the bulk region is slightly higher at the radial edge,especially near the bulk-sheath interface (figure 3(c)).This is because althoughVrf(r)shows a radial decrease under this condition,∣Ez∣in the bulk region is stronger at the radial edge due to the lower electron density and ∣κp∣ there (see equation (B5) in appendix B).Moreover, the radial increase of ∣Er∣also becomes more pronounced at 60 MHz than at 13.56 MHz due to the decrease of the surface wavelength,which is consistent with the analytical solution based on the uniform slab model[15, 16].Under the combined influence of∣Ez∣and∣Er∣ ,the maximum of ∣E∣appears at the radial edge of the bulk-sheath interface.In addition,Pz(r) exhibits a parabolic profile along the radial direction, which is similar toVrf(r) (see equation (B6) in appendix B), except for the peak due to the strongEz∣ ∣at the radial edge,as shown in figure 3(d).Besides,Psh(r) shows a center-high profile, which is caused by the similar distribution ofVrf(r) ands(r) (see equations (11) and(12)).AlthoughPr(r) becomes more significant than that at 13.56 MHz due to the stronger∣Er∣,Pr(r) is still negligible compared withPz(r)andPsh(r) ,because of the relatively low electron density.As a result, the maximum ofPtot(r) appears at the center, leading to the center-high profile ofne(r, 0)in figure 2(a).At 80 MHz, the radial wavenumber iskr=0.0711 −0.0026j withsave=0.16 cm,and the first node of surface wave appears atχ01/Re(kr)=33.8 cm,which is slightly different from the first node ofne(r, 0),Vrf(r) ands(r) in figure 2 (i.e.r=37.3 cm).This is because the dispersion relation is obtained by using the spatially averaged plasma parameters,which is similar to the approach in [29-31].Indeed, the electron density shows significant non-uniformity, which is clear from figure 2(a), so the predicted radial wavenumber by using the spatially averaged plasma parameters deviates slightly from the actual value.However, the evolution of radial wavenumbers and thus the surface wavelength with frequency agrees quite well with simulation results, and this can be used to qualitatively explain the EM effects under various discharge conditions.In figure 3(e), ∣E∣exhibits a peak near the node position of the bulk-sheath interface, which is again due to the combined effect of∣Ez∣and∣Er∣ ,as mentioned above.When the frequency increases further to 100 MHz, the radial wavenumber iskr=0.0998 −0.0052j withsave=0.14 cm,and the first node of surface wave takes place atχ01/Re (kr)=24.1 cm.Therefore, the first node ofne(r, 0),Vrf(r) ands(r) moves inward and the second peak appears in the chamber(figure 2).

3.2.Stop band effect

In this subsection, we focus on the transition from the standing wave effect to the stop band effect in the frequency range from 120 to 180 MHz at 50 mTorr and 100 W.In figure 4(a),ne(r, 0) shows a multiple-node profile at 120 MHz, i.e.ne(r, 0) has two nodes atr=19.1 cm andr=41.3 cm, respectively.The multiple-node density profile is again due to the significant standing wave effect,which has also been captured by experiments [19, 51, 52] and simulations [17, 38].As the frequency increases to 140 MHz, the stop band effect takes over the standing wave effect and dominates the radial profile of the electron density, which is similar to that obtained in [17, 19, 51, 52].As a result,ne(r, 0) is nearly zero in the region ofr＜4.5 cm, and it shows a radial increase towards the radial edge.As the frequency increases further to 160 MHz and 180 MHz, the stop band effect becomes more significant and the stop band region,wherene(r, 0)is nearly zero,expands tor＜19.5 cm andr＜26.0 cm, respectively.Note that the electron density distributions obtained in this work are different from[61,62],because of the low electron density, i.e.in the order of 1010cm−3.This value is calculated self-consistently by the EM global model under the electron power of 100 W.If a higher electron power is applied, the electron density increases, and thus the series resonance frequency may be several times higher than the driving frequency, as will be discussed below.Therefore,the stop band effect does not take place, and the discharge is dominated by the standing wave effect at VHF [61].Moreover, both the electron mean free path for argon under such discharge condition (i.e.about 0.6 cm) and the discharge gap (i.e.5 cm) are much smaller than the chamber radius (i.e.50 cm), so the plasma radial transport, which tends to decrease the stop band region,should be insignificant in these cases.In figures 4(b) and (c),Vrf(r) ands(r) show similar radial profiles tone(r, 0) at four frequencies.For instance,s(r)drops to the minimum value in the region wherene(r, 0) andVrf(r) are nearly zero.

The stop band effect derives from the radial damping of surface waves when they propagate from the radial edge to the center, which can be described by the radial damping lengthδr(see appendix A).At 120 MHz, the radial wavenumber iskr=0.1378 −0.0106j withsave=0.13 cm.The radial damping length isδr= 94.3 cm,which is almost twice the chamber radius R, indicating that there is no significant radial damping.Since Re(kr)R=6.89is larger than the second zero of the zero-order Bessel functionχ02= 5.52,ne(r, 0),Vrf(r) ands(r) exhibit two nodes at 120 MHz(figure 4), and so is the power area density, i.e.Pz(r),Psh(r)andPtot(r) (figure 5(b)).Besides, it is clear that ∣E ∣exhibits two peaks near the node positions at the bulk-sheath interface(figure 5(a)), as we discussed above.At 140 MHz, the radial wavenumber iskr=0.1911 −0.0223j withsave=0.12 cm.Since the radial damping lengthδr= 44.8 cm becomes smaller than R, the surface waves have significant damping.As a result,ne(r, 0),Vrf(r) ,ands(r)decrease rapidly from the radial edge to the center (figure 4), andVrf(r) and ∣E ∣are nearly zero in the region of ＜r4.5 cm.As the frequency increases to 160 MHz and 180 MHz,the radial wavenumbers arekr=0.2835 −0.0512j withsave=0.11 cm,andkr=0.4337 −0.0865j withsave=0.10 cm,and the radial damping lengths decrease strikingly toδr= 19.5 cm at 160 MHz andδr= 11.6 cm at 180 MHz,respectively.Due to the enhanced radial damping, the stop band region expands,as shown in figures 4(b), 5(e) and (g).

Figure 4.Radial profiles of (a) the electron density (i.e. n e(r , 0)) at z =0 cm, (b) the rf voltage (i.e.Vr f(r )) between two electrodes, and(c) the sheath width s ( r) for different frequencies: 120, 140, 160, and 180 MHz, at 50 mTorr and 100 W.

Figure 5.Spatial profiles of the electric feild E in the bulk plasma(first column), as well as the radial proflies of the capacitive power area density in the bulk Pz (r ) ,the inductive power area density in the bulk Pr (r ) ,the power area density in sheaths Ps h(r ) ,and the total power area density Pt ot(r )(second column)for different frequencies:120 MHz(a)and(b),140 MHz(c)and(d),160 MHz(e)and(f),and 180 MHz (g) and (h), at 50 mTorr and 100 W.

Figure 6 shows the variation of wavenumber with frequency in the frequency range of 0-300 MHz at the electron density of 8.4 ×109cm−3, the sheath width of 0.11 cm, and the electron temperature of 2.1 eV.Note that the spatially averaged electron density and the electron temperature are almost unchanged in this frequency range.Besides,the above sheath width is obtained by averagingsaveat 120, 140, 160,and 180 MHz (figure 4(c)), and the averaged value only differs fromsaveat each frequency by no more than 20%.Therefore, the dispersion relation curves in figure 6 can qualitatively describe the propagation characteristics of surface waves in this frequency range.As we can see, when the frequency is much lower than the series resonance frequencyi.e.in the frequency range ＜f120 MHz,k RRer( ) increases almost linearly, andis much smaller than 1, which is consistent with the approximate expression of the radial wavenumber at low density and low frequency[15].As a result, the surface wave could propagate from the radial edge to the center almost without damping.When the frequency is close tofres,Im(kr)Rincreases rapidly with frequency.Therefore, the radial damping becomes significant when the frequency is higher than 140 MHz, and this leads to the pronounced stop band effect.

Figure 7.Radial profiles of (a) the electron density (i.e. n e(r , 0)) at z =0 cm, (b) the rf voltage (i.e.Vr f(r )) between two electrodes, and(c) the sheath width s ( r) for different electron powers: 100, 200, 400, and 800 W, at 50 mTorr and 80 MHz.

3.3.Skin effect

The skin effect has been studied in a wide frequency range,i.e.40.7-180 MHz, in previous literatures [15, 17, 18, 32, 41, 51].In this subsection, 80 MHz is adopted to study the transition from the standing wave effect to the skin effect as the electron power increases from 100 to 800 W at 50 mTorr.In figure 7(a),it is obvious that the electron density increases significantly with electron power.For instance, the spatially averaged electron 1010cm−3at 800 W.Moreover,ne(r, 0) has a density node at density increases from8.4 ×109cm−3at 100 W to 6.7×r=37.3 cm and shows significant non-uniformity at 100 W,as we discussed above.At 200 W,the radial position of the density node almost keeps unchanged,but the decreasing trend near the radial center becomes less obvious.As the electron power increases to 400 W, the skin effect becomes significant due to the high electron density.Under the combined influence of the standing wave effect and the skin effect,ne(r, 0)near the radial center becomes uniform, and the density node becomes less obvious.When the electron power increases further to 800 W,the skin effect dominates the discharge,andne(r, 0)exhibits an edge-high profile, as shown in figure 7(a).In figures 7(b) and(c),Vrf(r) ands(r) show signifciant non-uniformity, and the node positions ofVrf(r) ands(r) are almost unchanged at different electron powers, because Re(kr)Ronly varies slightly from 3.56 at 100 W to 3.70 at 800 W.

The skin effect results from the shielding of surface waves by plasma at high electron density,and its importance could be evaluated by the skin depthδp(see appendix A).At 100 W,the axial wavenumber in the bulkkzp=0.0242 −0.1790j is obtained withne,ave=8.4 ×109cm−3andsave=0.16 cm.Since the skin depthδp= 5.59 cm is larger than the half width of the bulk region, the skin effect is negligible under this condition.Therefore, the standing wave effect has a dominant infulence, and leads to the radial non-uniformity ofne(r, 0),Vrf(r) ,ands(r) in figure 7, as we discussed above.At 200 W,the electron density increases tone,ave=1.7 ×1010cm−3, andsave=0.14 cm has a slight decrease.The axial wavenumber iskzp=0.0364 −0.2439j and the skin depth decreases toδp= 4.10 cm.Correspondingly,∣Er∣contributes more to ∣E ∣(figure 8(a))than that at 100 W(figure 3(e)),especially near the bulk-sheath interface and near the node position.Therefore,Pr(r)in figure 8(b)has a slight increase,and the radial decrease ofPtot(r) near the radial center becomes less obvious than that shown in figure 3(f), which leads to the similar variation ofne(r, 0) in fgiure 7(a).However, since the skin depth is still larger than the half width of the bulk region,the standing wave effect dominates the discharge under this condition.

As the electron power increases to 400 W, the electron density increases tone,ave=3.4 ×1010cm−3,and the sheath width has a further decrease, i.e.save=0.12 cm.Under this condition,kzp=0.0528 −0.3384j and the skin depth decreases toδp= 2.96 cm,which is close to the half width of the bulk region, indicating that the skin effect becomes significant.As a result,∣Er∣is the main component of ∣E∣in a large region near the bulk-sheath interface and near the node position, which is clear from the horizontal arrows in figure 8(c).Consequently,Pr(r) shows a signifciant increase and leads to the falt distribution ofPtot(r) andne(r, 0) near the radial center(fgiures 7(a) and 8(d)).At 800 W,kzp=0.0757 −0.4736j is obtained withne,ave=6.7 ×1010cm−3andsave=0.10cm,and the skin depthδp= 2.11 cm becomes smaller than the half width of the bulk region.Therefore,∣Er∣is responsible for the distribution of ∣E∣in the most region, as shown in figure 8(e).Moreover, the maximum value ofPr(r) at 34.7 cm becomes la rger than the sum ofPsh(r)andPz(r)at the center(fgiure 8(f)),leading to the edge-high proflie ofne(r, 0).

3.4.Comparison with experiments

To validate our model, the simulation results have been compared with experimental measurements.The cylindrical CCP reactor in the experiment has an inner diameter of 40 cm,containing two stainless steel parallel disk electrodes, 30 cm in diameter and separated by 4 cm.A highly purified(99.999%) argon gas is evenly injected into the discharge region through the showerhead-like top electrode.A signal generated by a function generator (Tektronix AFG31252) is amplified by a broadband power amplifier (AR, Model 1000A225) and then delivered through a Π-type matching network to the bottom electrode.The voltage waveform on the powered (bottom) electrode is measured by a high-frequency, high-voltage probe (Tektronix P5100A, 500 MHz bandwidth)and acquired via a digital oscilloscope(Tektronix MSO56).

Figure 8.Spatial profiles of the electric feild E in the bulk plasma(first column), as well as the radial proflies of the capacitive power area density in the bulk Pz (r ) ,the inductive power area density in the bulk Pr (r ) ,the power area density in sheaths Ps h(r ) ,and the total power area density Pt ot(r ) (second column) for different electron powers:200 W(a)and(b),400 W(c)and(d),and 800 W(e)and(f),at 50 mTorr and 80 MHz.

Figure 9 shows the measured and calculated electron density under 13.56, 30 and 100 MHz, with pressure of 3 Pa.In experiments, the source power is fixed at 40 W, and the measured rf voltages are 150, 50 and 30 V for the three frequencies, and these values are input into simulations as the voltages at the radial center.It is clear from figure 9(a)that the measured electron density is quite uniform at 13.56 and 30 MHz, except for a slightly higher value at the radial edge.At 100 MHz, the electron density exhibits a center-high profile, due to the standing wave effect.The evolution of the electron density with frequency in the simulation(figure 9(b))agrees well with the experimental measures,but there exists a small discrepancy in the absolute value,which may be caused by the simplified chamber geometry and chemistry set in the model.Besides, the radial decrease of the electron density at 100 MHz in experiments is more obvious than in simulations.This may be caused by the existence of the nonlinearly excited harmonics due to the geometrically asymmetry in experiments, which enhances the electron power near the center [25, 27, 50].

Figure 9.(a)Measured and(b) calculated electron density for argon discharges under different frequencies of 13.56, 30 and 100 MHz,with pressure of 3 Pa.

4.Conclusions

In this work, Maxwell equations are coupled with a radially localized global model and an analytical sheath model to study the EM effects and the underlying mechanisms for argon discharges in a geometrically symmetric cylindrical CCP reactor.The electron density is assumed to have a fixed non-uniform axial distribution, and the radial profile is determined by applying the radially localized global model,which allows us to improve the simulation efficiency and meanwhile take the non-uniform density profile into account.

First, we investigated the standing wave effect in the frequency range from 13.56 to 100 MHz at 100 W and 50 mTorr.At low frequencies, both the vacuum wavelength and the sheath width are very large,leading to the fact that the surface wavelength is much larger than the chamber radius.Therefore, no significant standing wave effect is observed,and the electron density, the rf voltage, and the sheath width are radially uniform.As the frequency increases, both the vacuum wavelength and the sheath width show a remarkable decrease, and so is the surface wavelength, which becomes comparable to or even smaller than the chamber radius.As a result, the standing wave effect becomes pronounced and gives rise to the center-high or even multiple-node profiles.

Besides, the stop band effect has been illustrated in the frequency range from 120 to 180 MHz at 100 W and 50 mTorr.Under the conditions we studied, the frequency only has a slight influence on the spatially averaged electron density and the sheath width, so the series resonance frequency almost keeps constant in this frequency range.When the frequency is close to or higher than the series resonance frequency, a significant radial damping of surface waves takes place,resulting in that they cannot propagate to the radial center.Therefore,the stop band effect dominates the discharge and results in a region near the radial center where the electron density is almost zero.

Finally, the skin effect has been demonstrated in the electron power range from 100 to 800 W at 80 MHz and 50 mTorr.At low electron powers, the electron density is relatively low,and the skin depth is much larger than the half width of the bulk plasma,so the skin effect is not obvious.As the electron power increases, the increased electron density leads to the reduced skin depth.Thus,the skin effect becomes pronounced,and the inductive power area density at the radial edge has a remarkable increase.Under the combined influence of the skin effect and the standing wave effect, the density profile becomes radially uniform at 400 W.As the power increases further, the skin effect dominates the discharge and gives rise to the edge-high profile of the electron density.

Acknowledgments

This work was supported by National Natural Science Foundation of China (NSFC) (Nos.11935005, 11875101,12020101005,12005035),the Natural Science Foundation of Liaoning Province (No.2020-MS-114), the China Postdoctoral Science Foundation (No.2020M670741), and the Fundamental Research Funds for the Central Universities(No.DUT20LAB201).The author Jiankai Liu also appreciates the financial support from the China Scholarship Council.

Appendix A.

The dispersion relations given by the uniform slab model are solved with the spatially averaged plasma parameters to understand the propagation characteristics of surface waves in detail.Since the higher symmetric modes rapidly damp out near the radial edge, and the antisymmetric modes do not exist in this work, only the dispersion relations for the first symmetric mode of surface waves are solved [15, 16, 31]

wherekr,kzpandkzsare the radial wavenumber, the axial wavenumber in the bulk and the axial wavenumber in the sheath, respectively.dave=l−saveis the radially averaged half width of the bulk plasma.

The above equations are calculated using the Matlab solver fsolve,and some important parameters,which describe the propagation of surface waves, are determined.For instance, Re(kr)Rcould help to explain the radial variation of the fields caused by the finite wavelengths.The radial damping lengthδr=1 /Im(kr)describes the radial decay of the fields from the radial edge to the center.The skin depthδp=1 /∣I m(kzp)∣illustrates the axial decay of the fields when they penetrate into the bulk plasma.

Appendix B.

To demonstrate the radial profiles ofpz(r,z) andEzat the bulk-sheath interface, two positions within the sheath region(ra,−d−) and(rb,−d−) ,and two positions within the bulk region(ra,−d+) and(rb,−d+) ,are adopted.

At the bulk-sheath interface, the axial electrical displacement is continuous

Using the expression forκpin the bulk and in the sheath regions, we have

and the ratio of∣E z(r,−d+)∣atr=raandr=rbis given by

BecauseEzin the sheath region almost keeps constant along the axial direction, we find

Substituting the expressions forωp,s(r) ,and∣Ez(r,−d−)∣into equation (B3), we obtain

Using equations (6) with (B5) and the expression forσp, the ratio ofpz(r,−d+) atr=raandr=rbis achieved

It is clear from equation (B6) thatpz(r,−d+) is proportional toVrf(r)1/2.

ORCID iDs