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

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{\phi }}{\partial z}=-{\rm{j}}\omega {\varepsilon }_{0}{\mu }_{0}{\kappa }_{{\rm{p}}}{E}_{r},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{r}\displaystyle \frac{\partial \left(r{B}_{\phi }\right)}{\partial r}={\rm{j}}\omega {\varepsilon }_{0}{\mu }_{0}{\kappa }_{{\rm{p}}}{E}_{z},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{r}}{\partial z}-\displaystyle \frac{\partial {E}_{z}}{\partial r}=-{\rm{j}}\omega {B}_{\phi }.\end{eqnarray} \tag{ 3 }$$

Here, ${B}_{\phi }$ is the angular component of the magnetic field, ${E}_{r}$ is the radial component of the electric field, ${E}_{z}$ is the axial component of the electric field, $\omega$ is the driving angular frequency, ${\varepsilon }_{0}$ is the vacuum permittivity, ${\mu }_{0}$ is the vacuum permeability, and ${\kappa }_{{\rm{p}}}$ is the plasma relative permittivity.

In the bulk region, the relative permittivity is ${\kappa }_{{\rm{p}}}=1-{\omega }_{{\rm{p}}}^{2}/\omega \left(\omega -{\rm{j}}{\nu }_{{\rm{m}}}\right),$ where ${\nu }_{{\rm{m}}}$ is the electron-neutral momentum transfer frequency and ${\omega }_{{\rm{p}}}={\left({e}^{2}{n}_{{\rm{e}}}/{\varepsilon }_{0}{m}_{{\rm{e}}}\right)}^{1/2}$ is the electron plasma frequency with the element charge e, the electron mass ${m}_{{\rm{e}}},$ and the electron density ${n}_{{\rm{e}}}.$ The electron density ${n}_{{\rm{e}}}$ in the bulk region is assumed to be spatially non-uniform along the axial direction

$$\begin{eqnarray}&&{n}_{{\rm{e}}}(r,z)={n}_{{\rm{e}}}(r,0){\left[1-{(z/l)}^{2}\right]}^{1/2},\end{eqnarray} \tag{ 4 }$$

which is a good approximation at low pressures [1, 55]. The radial profile of the electron density (i.e. ${n}_{{\rm{e}}}(r,0)$) at $z=0$ cm is determined by the following radially localized global model.

In the sheath region, the relative permittivity varies along the radial direction ${\kappa }_{{\rm{p}}}={s}_{0}/s(r).$ This allows the model to use a simulation-convenient fixed sheath thickness ${s}_{0},$ and meanwhile to equivalently take the propagation of surface waves in an actual radial varying sheath $s(r)$ into account [17, 28]. Here, ${s}_{0}$ is determined by the position of ${u}_{{\rm{i}}}={u}_{{\rm{B}}},$ where ${u}_{{\rm{i}}}$ is the ion velocity and ${u}_{{\rm{B}}}=\sqrt{e{T}_{{\rm{e}}}/M}$ is the Bohm velocity with the electron temperature ${T}_{{\rm{e}}}$ and the ion mass M (equation (5.3.5) in [1]). Note that ${s}_{0}$ influences 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 width $s(r),$ which is determined by the following analytical sheath model.

From equations (1)–(3), a Helmholtz equation for ${B}_{\phi }$ is obtained

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial r}\left[\displaystyle \frac{1}{{\kappa }_{{\rm{p}}}r}\displaystyle \frac{\partial \left(r{B}_{\phi }\right)}{\partial r}\right]+\displaystyle \frac{\partial }{\partial z}\left(\displaystyle \frac{1}{{\kappa }_{{\rm{p}}}}\displaystyle \frac{\partial {B}_{\phi }}{\partial z}\right)+{k}_{0}^{2}{B}_{\phi }=0,\end{eqnarray} \tag{ 5 }$$

where ${k}_{0}=\omega /c$ is the vacuum wavenumber, with c the speed of light in vacuum. The boundary conditions for ${B}_{\phi }$ are as follows: (1) ${B}_{\phi }=0$ at the radial center. (2) ${{\bf{e}}}_{{\rm{n}}}\cdot {\rm{\nabla }}{B}_{\phi }=0$ at two electrodes with ${{\bf{e}}}_{{\rm{n}}}$ the unit normal vector. (3) ${B}_{\phi }={\mu }_{0}{I}_{{\rm{rf}}}/2\pi R$ at the radial edge, where ${I}_{{\rm{rf}}}$ is the source current. Substituting ${B}_{\phi }$ into equations (1) and (2), the EM fields can be determined, and the capacitive and inductive power densities in the bulk region are given by

$$\begin{eqnarray}&&{p}_{z}(r,z)=\displaystyle \frac{1}{2}\mathrm{Re}\left({\sigma }_{{\rm{p}}}\right){\left|{E}_{z}\right|}^{2},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{p}_{r}(r,z)=\displaystyle \frac{1}{2}\mathrm{Re}\left({\sigma }_{{\rm{p}}}\right){\left|{E}_{r}\right|}^{2},\end{eqnarray} \tag{ 7 }$$

where ${\sigma }_{{\rm{p}}}={\varepsilon }_{0}{\omega }_{{\rm{p}}}^{2}/\left({v}_{{\rm{m}}}-{\rm{j}}\omega \right)$ is the plasma conductivity.

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 velocity ${u}_{{\rm{B}}},$ and then they are lost at the electrodes

$$\begin{eqnarray}&&N{K}_{{\rm{iz}}}\displaystyle {\int }_{-d}^{d}{n}_{{\rm{e}}}(r,z){\rm{d}}z=2{u}_{{\rm{B}}}{h}_{l}{n}_{{\rm{e}}}(r,0).\end{eqnarray} \tag{ 8 }$$

Here $N={P}_{{\rm{g}}}/{k}_{{\rm{B}}}{T}_{0}$ is the background gas density, with the gas pressure ${P}_{{\rm{g}}},$ the gas temperature ${T}_{0}=293\,{\rm{K}},$ and the Boltzmann constant ${k}_{{\rm{B}}}.$ ${K}_{{\rm{iz}}}$ is the ionization rate coefficient, as shown in table 1. ${h}_{l}={\rm{0.86}}/\sqrt{3+d/{\lambda }_{{\rm{i}}}}$ is the axial edge to center density ratio, with the ion-neutral mean free path ${\lambda }_{{\rm{i}}}=1/330{P}_{{\rm{g}}}.$ Due to the fixed axial density profile (see equation (4)), the ratio of $\displaystyle {\int }_{-d}^{d}{n}_{{\rm{e}}}(r,z){\rm{d}}z$ to ${n}_{{\rm{e}}}(r,0)$ in equation (8) is a constant, so ${K}_{{\rm{iz}}}$ and thus ${T}_{{\rm{e}}}$ is independent on the radial position.

Table 1. Electron impact reaction rate coefficients in terms of the electron temperature ${T}_{{\rm{e}}}(\mathrm{eV})$ (section 3.5 in [1]).

Reaction		Collision rate coefficient (${{\rm{m}}}^{3}\,{{\rm{s}}}^{-1}$)
Elastic collision	${\rm{e}}+{\rm{Ar}}\to {\rm{e}}+{\rm{Ar}}$	${K}_{{\rm{elas}}}=2.336\times {10}^{-14}{T}_{{\rm{e}}}^{1.609}{{\rm{e}}}^{0.0618{\left(\mathrm{ln}{T}_{{\rm{e}}}\right)}^{2}-0.1171{\left(\mathrm{ln}{T}_{{\rm{e}}}\right)}^{3}}$
Ionization	${\rm{e}}+{\rm{Ar}}\to {\rm{2e}}+{{\rm{Ar}}}^{+}$	${K}_{{\rm{iz}}}=2.34\times {10}^{-14}{{T}_{{\rm{e}}}}^{0.59}{{\rm{e}}}^{-17.44/{T}_{{\rm{e}}}}$
Excitation	${\rm{e}}+{\rm{Ar}}\to {\rm{e}}+{\rm{Ar}}*$	${K}_{{\rm{exc}}}=2.48\times {10}^{-14}{T}_{{\rm{e}}}^{0.33}{{\rm{e}}}^{-12.78/{T}_{{\rm{e}}}}$

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

$$\begin{eqnarray}&&2{n}_{{\rm{e}}}(r,0){h}_{l}{u}_{{\rm{B}}}{\varepsilon }_{{\rm{e}}}e={P}_{{\rm{tot}}}(r).\end{eqnarray} \tag{ 9 }$$

Here ${ {\mathcal E} }_{{\rm{e}}}={ {\mathcal E} }_{{\rm{c}}}+7.2{T}_{{\rm{e}}}$ is the total electron energy loss, where $7.2{T}_{{\rm{e}}}$ represents the electron energy loss when the electrons cross the sheath to the electrode [20]. ${ {\mathcal E} }_{{\rm{c}}}=15.76+12.14{K}_{{\rm{exc}}}/{K}_{{\rm{iz}}}+3{m}_{{\rm{e}}}{T}_{{\rm{e}}}{K}_{{\rm{elas}}}/M{K}_{{\rm{iz}}}$ is the collisional electron energy loss per electron-ion pair created, ${K}_{{\rm{exc}}}$ is the average energy loss-weighted excitation rate coefficient for average excitation energy 12.14 eV [1], and ${K}_{{\rm{elas}}}$ is the elastic collision rate coefficient (see table 1). ${P}_{{\rm{tot}}}(r)={P}_{z}(r)+{P}_{r}(r)\,+2{S}_{{\rm{stoc}}}(r)+2{S}_{{\rm{ohm}},{\rm{sh}}}(r)$ is the total electron power area density. ${P}_{z}(r)=\displaystyle {\int }_{-d}^{d}{p}_{z}(r,z){\rm{d}}z$ and ${P}_{r}(r)=\displaystyle {\int }_{-d}^{d}{p}_{r}(r,z){\rm{d}}z$ are the capacitive and inductive power area densities in the bulk region. ${S}_{{\rm{stoc}}}(r)$ and ${S}_{{\rm{ohm}},{\rm{sh}}}(r)$ are the stochastic heating and the Ohmic heating in each sheath, and they are calculated by the analytical sheath model below.

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

$$\begin{eqnarray}&&{s}_{\max }(r)={2}^{1/4}{\left(\displaystyle \frac{0.82{\varepsilon }_{0}}{e}\right)}^{1/2}\displaystyle \frac{\overline{V}{(r)}^{3/4}}{{{h}_{l}}^{1/2}{n}_{{\rm{e}}}{(r,0)}^{1/2}{{T}_{{\rm{e}}}}^{1/4}},\end{eqnarray} \tag{ 10 }$$

where $\overline{V}(r)=0.83V(r)$ is the DC sheath voltage and $V(r)=\left|\displaystyle {\int }_{-l}^{-d}{E}_{z}\,{\rm{d}}z\right|=\left|\displaystyle {\int }_{d}^{l}{E}_{z}\,{\rm{d}}z\right|$ is the rf voltage across each sheath. The time averaged sheath width is given by $s(r)={s}_{\max }(r)/1.23.$ Note that the minimum sheath width ${s}_{\min }=2.61{\lambda }_{{\rm{De}}}$ is added to the time averaged sheath width $s(r)$ to avoid that $s(r)$ is nearly zero at the radial positions where the sheath voltage drop $V(r)$ is very small, and ${\lambda }_{{\rm{De}}}={({\varepsilon }_{0}{T}_{{\rm{e}}}/e{h}_{l}{n}_{{\rm{e}}})}^{1/2}$ is the Debye length. Moreover, the stochastic heating and the Ohmic heating in each sheath are given by

$$\begin{eqnarray}&&{S}_{{\rm{stoc}}}(r)=0.45{\left(\displaystyle \frac{{m}_{{\rm{e}}}}{e}\right)}^{1/2}{\varepsilon }_{0}{\omega }^{2}{{T}_{{\rm{e}}}}^{1/2}V(r),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{S}_{\mathrm{ohm},\mathrm{sh}}(r)=0.407\displaystyle \frac{{m}_{{\rm{e}}}}{2\pi e}{\varepsilon }_{0}{\omega }^{2}{v}_{{\rm{m}}}{s}_{\max }(r)V(r).\end{eqnarray} \tag{ 12 }$$

First, the EM fields are determined by solving equations (1)–(3) with the initial ${n}_{{\rm{e}}}(r,z)$ and $s(r),$ and then $V(r),$ ${p}_{z}(r,z),$ and ${p}_{r}(r,z)$ are obtained. Subsequently, the analytical sheath model is implemented, the sheath width $s(r)$ is updated, and the sheath heating ${S}_{\mathrm{ohm},\mathrm{sh}}(r)$ and ${S}_{{\rm{stoc}}}(r)$ are calculated. Then, the energy balance equation (9) is solved at each radial position for ${n}_{{\rm{e}}}(r,0),$ and the axial distribution ${n}_{{\rm{e}}}(r,z)$ is determined according to equation (4). Finally, the electron temperature ${T}_{{\rm{e}}}$ is obtained by solving the particle balance equation (8). In each iteration, the source current ${I}_{{\rm{rf}}}$ is adjusted to ensure that the calculated total electron power, which is obtained by integrating ${P}_{{\rm{tot}}}(r)$ over the whole electrode, equals to the fixed value. The iteration continues until the energy balance is achieved at each radial position.

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. ${n}_{{\rm{e}}}(r,0)$) at $z=0$ cm for various frequencies. It is clear that ${n}_{{\rm{e}}}(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, ${n}_{{\rm{e}}}(r,0)$ shows an obvious center-high profile, due to the standing wave effect. At 80 MHz, ${n}_{{\rm{e}}}(r,0)$ first decreases along the radial direction, with a density node appearing at $r=37.3\,{\rm{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 $r={\rm{26.0}}\,{\rm{cm}},$ and the second density peak appears at $r={\rm{40.0}}\,{\rm{cm}},$ indicating that the standing wave effect becomes more pronounced.

Zoom In Zoom Out Reset image size

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 density ${n}_{{\rm{e}},\mathrm{ave}}$ and the radially averaged sheath width ${s}_{{\rm{ave}}}$ [29–31] (see appendix A). Since the total electron power is fixed as we mentioned above, ${n}_{{\rm{e}},\text{ave}}=8.4\times {10}^{9}\,{{\rm{cm}}}^{-3}$ stays unchanged at four frequencies. At 13.56 MHz, the obtained radial wavenumber is ${k}_{r}=0.0050-0.000042{\rm{j}}$ with ${s}_{{\rm{ave}}}=1.15\,{\rm{cm}}.$ Note that $\mathrm{Re}({k}_{r})R=0.25$ is much smaller than the first zero of the zero-order Bessel function ${\chi }_{01}=2.40,$ indicating that the surface wavelength (i.e. $\lambda =2\pi /\mathrm{Re}({k}_{r})=1257\,{\rm{cm}}$) is much longer than the chamber radius R. This gives rise to the uniform rf voltage between two electrodes ${V}_{{\rm{rf}}}(r)$ and the uniform sheath width $s(r)$ (figures 2(b) and (c)), and consequently the uniform sheath heating (i.e. ${P}_{{\rm{sh}}}(r)$ in figure 3(b)). As the frequency increases, ${V}_{{\rm{rf}}}(r)$ and $s(r)$ show overall downward trend and significant radial non-uniformity, which is consistent with the radial profiles of ${n}_{{\rm{e}}}(r,0),$ and this will be discussed below.

Zoom In Zoom Out Reset image size

The spatial profiles of the electric field ${\bf{E}}$ 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 of ${\bf{E}},$ respectively. To reveal the electric field in the bulk region more clearly, we have excluded the sheath regions in all the figures of ${\bf{E}}.$ In the second column of figure 3, the power area density in sheaths ${P}_{{\rm{sh}}}(r)$ is calculated by ${P}_{{\rm{sh}}}(r)=2{S}_{{\rm{stoc}}}(r)+2{S}_{\mathrm{ohm},\mathrm{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), $| {\bf{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 of ${P}_{z}(r)$ (figure 3(b)). In the axial direction, $| {\bf{E}}|$ shows a significant decrease from the bulk-sheath interface to $z=0\,{\rm{cm}},$ which is caused by the higher electron density and thus the higher $| {\kappa }_{{\rm{p}}}|$ at the axial center. Besides, ${E}_{r}$ and the corresponding ${P}_{r}(r)$ are very weak under this condition, because of the negligible skin effect at the relatively low electron density. As a result, ${P}_{{\rm{tot}}}(r)$ is mainly determined by ${P}_{{\rm{sh}}}(r)$ and ${P}_{z}(r),$ and the radially uniform ${P}_{{\rm{tot}}}(r)$ leads to the similar distribution of ${n}_{{\rm{e}}}(r,0)$ in figure 2(a). Note that ${P}_{{\rm{sh}}}(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.

As the frequency increases to 60 MHz, the radial wavenumber is ${k}_{r}=0.0459-0.0011{\rm{j}}$ with ${s}_{{\rm{ave}}}=0.20\,{\rm{cm}}.$ Note that $\mathrm{Re}({k}_{r})R=2.30$ is close to ${\chi }_{01},$ indicating that the surface wavelength (i.e. $\lambda =137\,{\rm{cm}}$) is comparable to the chamber size, and the standing wave effect becomes significant. As a result, ${V}_{{\rm{rf}}}(r)$ and $s(r)$ are characterized by obvious center-high profiles in figures 2(b) and (c). However, $| {\bf{E}}|$ in the bulk region is slightly higher at the radial edge, especially near the bulk-sheath interface (figure 3(c)). This is because although ${V}_{{\rm{rf}}}(r)$ shows a radial decrease under this condition, $| {E}_{z}|$ in the bulk region is stronger at the radial edge due to the lower electron density and $| {\kappa }_{{\rm{p}}}|$ there (see equation (B5) in appendix B). Moreover, the radial increase of $| {E}_{r}|$ 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 $| {E}_{z}|$ and $| {E}_{r}| ,$ the maximum of $| {\bf{E}}|$ appears at the radial edge of the bulk-sheath interface. In addition, ${P}_{z}(r)$ exhibits a parabolic profile along the radial direction, which is similar to ${V}_{{\rm{rf}}}(r)$ (see equation (B6) in appendix B), except for the peak due to the strong $| {E}_{z}|$ at the radial edge, as shown in figure 3(d). Besides, ${P}_{{\rm{sh}}}(r)$ shows a center-high profile, which is caused by the similar distribution of ${V}_{{\rm{rf}}}(r)$ and $s(r)$ (see equations (11) and (12)). Although ${P}_{r}(r)$ becomes more significant than that at 13.56 MHz due to the stronger $| {E}_{r}| ,$ ${P}_{r}(r)$ is still negligible compared with ${P}_{z}(r)$ and ${P}_{{\rm{sh}}}(r),$ because of the relatively low electron density. As a result, the maximum of ${P}_{{\rm{tot}}}(r)$ appears at the center, leading to the center-high profile of ${n}_{{\rm{e}}}(r,0)$ in figure 2(a). At 80 MHz, the radial wavenumber is ${k}_{r}=0.0711-0.0026{\rm{j}}$ with ${s}_{{\rm{ave}}}=0.16\,{\rm{cm}},$ and the first node of surface wave appears at ${\chi }_{01}/\mathrm{Re}({k}_{r})=33.8\,{\rm{cm}},$ which is slightly different from the first node of ${n}_{{\rm{e}}}(r,0),$ ${V}_{{\rm{rf}}}(r)$ and $s(r)$ in figure 2 (i.e. $r=37.3\,{\rm{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), $| {\bf{E}}|$ exhibits a peak near the node position of the bulk-sheath interface, which is again due to the combined effect of $| {E}_{z}|$ and $| {E}_{r}| ,$ as mentioned above. When the frequency increases further to 100 MHz, the radial wavenumber is ${k}_{r}=0.0998-0.0052{\rm{j}}$ with ${s}_{{\rm{ave}}}=0.14\,{\rm{cm}},$ and the first node of surface wave takes place at ${\chi }_{01}/\mathrm{Re}({k}_{r})=24.1\,{\rm{cm}}.$ Therefore, the first node of ${n}_{{\rm{e}}}(r,0),$ ${V}_{{\rm{rf}}}(r)$ and $s(r)$ moves inward and the second peak appears in the chamber (figure 2).

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), ${n}_{{\rm{e}}}(r,0)$ shows a multiple-node profile at 120 MHz, i.e. ${n}_{{\rm{e}}}(r,0)$ has two nodes at $r=19.1\,{\rm{cm}}$ and $r=41.3\,{\rm{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, ${n}_{{\rm{e}}}(r,0)$ is nearly zero in the region of $r\lt 4.5\,{\rm{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, where ${n}_{{\rm{e}}}(r,0)$ is nearly zero, expands to $r\lt 19.5\,{\rm{cm}}$ and $r\lt 26.0\,{\rm{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 ${10}^{10}\,{{\rm{cm}}}^{-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), ${V}_{{\rm{rf}}}(r)$ and $s(r)$ show similar radial profiles to ${n}_{{\rm{e}}}(r,0)$ at four frequencies. For instance, $s(r)$ drops to the minimum value in the region where ${n}_{{\rm{e}}}(r,0)$ and ${V}_{{\rm{rf}}}(r)$ are nearly zero.

Zoom In Zoom Out Reset image size

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 ${\delta }_{r}$ (see appendix A). At 120 MHz, the radial wavenumber is ${k}_{r}=0.1378-0.0106{\rm{j}}$ with ${s}_{{\rm{ave}}}=0.13\,{\rm{cm}}.$ The radial damping length is ${\delta }_{r}=94.3\,{\rm{cm}},$ which is almost twice the chamber radius R, indicating that there is no significant radial damping. Since $\mathrm{Re}({k}_{r})R=6.89$ is larger than the second zero of the zero-order Bessel function ${\chi }_{02}=5.52,$ ${n}_{{\rm{e}}}(r,0),$ ${V}_{{\rm{rf}}}(r)$ and $s(r)$ exhibit two nodes at 120 MHz (figure 4), and so is the power area density, i.e. ${P}_{z}(r),$ ${P}_{{\rm{sh}}}(r)$ and ${P}_{{\rm{tot}}}(r)$ (figure 5(b)). Besides, it is clear that $| {\bf{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 is ${k}_{r}=0.1911-0.0223{\rm{j}}$ with ${s}_{{\rm{ave}}}=0.12\,{\rm{cm}}.$ Since the radial damping length ${\delta }_{r}=44.8\,{\rm{cm}}$ becomes smaller than R, the surface waves have significant damping. As a result, ${n}_{{\rm{e}}}(r,0),$ ${V}_{{\rm{rf}}}(r),$ and $s(r)$ decrease rapidly from the radial edge to the center (figure 4), and ${V}_{{\rm{rf}}}(r)$ and $| {\bf{E}}|$ are nearly zero in the region of $r\lt 4.5\,{\rm{cm}}.$ As the frequency increases to 160 MHz and 180 MHz, the radial wavenumbers are ${k}_{r}=0.2835-0.0512{\rm{j}}$ with ${s}_{{\rm{ave}}}=0.11\,{\rm{cm}},$ and ${k}_{r}=0.4337-0.0865{\rm{j}}$ with ${s}_{{\rm{ave}}}=0.10\,{\rm{cm}},$ and the radial damping lengths decrease strikingly to ${\delta }_{r}=19.5\,{\rm{cm}}$ at 160 MHz and ${\delta }_{r}=11.6\,{\rm{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).

Zoom In Zoom Out Reset image size

Figure 6 shows the variation of wavenumber with frequency in the frequency range of 0–300 MHz at the electron density of $8.4\times {10}^{9}\,{{\rm{cm}}}^{-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 averaging ${s}_{{\rm{ave}}}$ at 120, 140, 160, and 180 MHz (figure 4(c)), and the averaged value only differs from ${s}_{{\rm{ave}}}$ at 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 frequency ${f}_{{\rm{res}}}={\omega }_{{\rm{p}}}/2\pi \sqrt{1+d/s}={\rm{171MHz}},$ i.e. in the frequency range $f\lt 120\,{\rm{MHz}},$ $\mathrm{Re}({k}_{r})R$ increases almost linearly, and ${\rm{Im}}({k}_{r})R$ is much smaller than 1, which is consistent with the approximate expression of the radial wavenumber at low density and low frequency ${k}_{r}\approx {k}_{0}{(1+d/s)}^{1/2}$ [15]. As a result, the surface wave could propagate from the radial edge to the center almost without damping. When the frequency is close to ${f}_{{\rm{res}}},$ ${\rm{Im}}({k}_{r})R$ increases 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.

Zoom In Zoom Out Reset image size

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 density increases from $8.4\times {10}^{9}\,{{\rm{cm}}}^{-3}$ at 100 W to $6.7\,\times {10}^{10}\,{{\rm{cm}}}^{-3}$ at 800 W. Moreover, ${n}_{{\rm{e}}}(r,0)$ has a density node at $r=37.3\,{\rm{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, ${n}_{{\rm{e}}}(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, and ${n}_{{\rm{e}}}(r,0)$ exhibits an edge-high profile, as shown in figure 7(a). In figures 7(b) and (c), ${V}_{{\rm{rf}}}(r)$ and $s(r)$ show significant non-uniformity, and the node positions of ${V}_{{\rm{rf}}}(r)$ and $s(r)$ are almost unchanged at different electron powers, because $\mathrm{Re}({k}_{r})R$ only varies slightly from 3.56 at 100 W to 3.70 at 800 W.

Zoom In Zoom Out Reset image size

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 ${\delta }_{{\rm{p}}}$ (see appendix A). At 100 W, the axial wavenumber in the bulk ${k}_{z{\rm{p}}}=0.0242-0.1790{\rm{j}}$ is obtained with ${n}_{{\rm{e}},\text{ave}}=8.4\times {10}^{9}\,{{\rm{cm}}}^{-3}$ and ${s}_{{\rm{ave}}}=0.16\,{\rm{cm}}.$ Since the skin depth ${\delta }_{{\rm{p}}}=5.59\,{\rm{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 influence, and leads to the radial non-uniformity of ${n}_{{\rm{e}}}(r,0),$ ${V}_{{\rm{rf}}}(r),$ and $s(r)$ in figure 7, as we discussed above. At 200 W, the electron density increases to ${n}_{{\rm{e}},\mathrm{ave}}=1.7\times {10}^{10}\,{{\rm{cm}}}^{-3},$ and ${s}_{{\rm{ave}}}=0.14\,{\rm{cm}}$ has a slight decrease. The axial wavenumber is ${k}_{z{\rm{p}}}=0.0364-0.2439{\rm{j}}$ and the skin depth decreases to ${\delta }_{{\rm{p}}}=4.10\,{\rm{cm}}.$ Correspondingly, $\left|{E}_{r}\right|$ contributes more to $| {\bf{E}}|$ (figure 8(a)) than that at 100 W (figure 3(e)), especially near the bulk-sheath interface and near the node position. Therefore, ${P}_{r}(r)$ in figure 8(b) has a slight increase, and the radial decrease of ${P}_{{\rm{tot}}}(r)$ near the radial center becomes less obvious than that shown in figure 3(f), which leads to the similar variation of ${n}_{{\rm{e}}}(r,0)$ in figure 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.

Zoom In Zoom Out Reset image size

As the electron power increases to 400 W, the electron density increases to ${n}_{{\rm{e}},\mathrm{ave}}=3.4\times {10}^{10}\,{{\rm{cm}}}^{-3},$ and the sheath width has a further decrease, i.e. ${s}_{{\rm{ave}}}=0.12\,{\rm{cm}}.$ Under this condition, ${k}_{z{\rm{p}}}=0.0528-0.3384{\rm{j}}$ and the skin depth decreases to ${\delta }_{{\rm{p}}}=2.96\,{\rm{cm}},$ which is close to the half width of the bulk region, indicating that the skin effect becomes significant. As a result, $\left|{E}_{r}\right|$ is the main component of $| {\bf{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, ${P}_{r}(r)$ shows a significant increase and leads to the flat distribution of ${P}_{{\rm{tot}}}(r)$ and ${n}_{{\rm{e}}}(r,0)$ near the radial center (figures 7(a) and 8(d)). At 800 W, ${k}_{z{\rm{p}}}=0.0757-0.4736{\rm{j}}$ is obtained with ${n}_{{\rm{e}},\mathrm{ave}}=6.7\times {10}^{10}\,{{\rm{cm}}}^{-3}$ and ${s}_{{\rm{ave}}}=0.10{\rm{cm}},$ and the skin depth ${\delta }_{{\rm{p}}}=2.11\,{\rm{cm}}$ becomes smaller than the half width of the bulk region. Therefore, $\left|{E}_{r}\right|$ is responsible for the distribution of $| {\bf{E}}|$ in the most region, as shown in figure 8(e). Moreover, the maximum value of ${P}_{r}(r)$ at 34.7 cm becomes larger than the sum of ${P}_{{\rm{sh}}}(r)$ and ${P}_{z}(r)$ at the center (figure 8(f)), leading to the edge-high profile of ${n}_{{\rm{e}}}(r,0).$

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 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].

Zoom In Zoom Out Reset image size