2D fluid-analytical simulation of electromagnetic effects in low pressure, high frequency electronegative capacitive discharges

A fast 2D fluid-analytical multifrequency CCP reactor model was developed using the finite elements simulation tool COMSOL in the MATLAB numerical computing environment. Since this model has been fully described in reference [10], we only give a brief description below. As shown in figure 1, the CCP model has an axisymmetric cylindrical geometry with center of symmetry at r  =  0 (z-axis). For this study, the plasma chamber is 5 cm in height and 25 cm in radius. There is a 15 cm radius showerhead gas inlet at the top and radial center, and a 3 cm width gas outlet at the bottom and radial edge. The wafer electrode has a radius of 15 cm, and is insulated from other conducting surfaces by a 2.5 cm wide quartz dielectric spacer. The bulk plasma region is surrounded by a sheath region with a nominal width of s₀  =  0.5 cm. The feedstock gas in the simulations was Cl₂, and we used the chlorine reaction rate set compiled by Thorsteinsson and Gudmundsson [21].

Zoom In Zoom Out Reset image size

The model treats each region of the reactor as a dielectric slab. The free-space magnetic permeability $\mu ={{\mu}_{0}}$ is assumed everywhere while $\epsilon =\kappa {{\epsilon}_{0}}$ depends on the relative dielectric constant κ of each region. The sheath relative dielectric constant ${{\kappa}_{\text{s}}}$ is initially set to one but is actually a function of the local electric field and plasma parameters as discussed in [10, 22]. The quartz dielectric spacer has $\kappa =4$. In the plasma region, the relative dielectric constant is

$$\begin{eqnarray}&&{{\kappa}_{\text{p}}}=1-\frac{\omega _{\text{pe}}^{2}}{\omega \left(\omega -j{{\nu}_{\text{m}}}\right)},\end{eqnarray} \tag{ 1 }$$

where ω is the applied radian rf frequency, ${{\nu}_{\text{m}}}$ is the electron-neutral momentum transfer collision frequency, and ${{\omega}_{\text{pe}}}={{({{n}_{\text{e}}}{{e}^{2}}/({{\epsilon}_{0}}{{m}_{\text{e}}}))}^{1/2}}$ is the electron plasma frequency with n_e and m_e the electron density and mass respectively. Note that ${{\kappa}_{\text{p}}}$ is complex with a dissipative imaginary component, so that our model treats the plasma region as a lossy dielectric.

The simulation consists of four basic parts: (1) an EM model which uses the time-independent Helmholtz equation to solve for the capacitive fields in the linearized frequency domain; (2) an ambipolar, quasineutral bulk plasma model which solves the time-dependent fluid equations for ion continuity and electron energy balance; (3) an analytical sheath model which solves for the sheath parameters (i.e. sheath voltage, sheath width and ${{\kappa}_{\text{s}}}$); (4) a gas flow model which solves for the steady-state composition, pressure, temperature and velocity of a reactive gas. A flow chart depicting the CCP model simulation is shown in figure 2. The total simulation time for Cl₂ is about 60–90 min on a moderate workstation with 2.7 GHZ CPU and 12 GB of memory. The finite elements mesh used in the simulation consisted of 5846 triangular elements and 3108 mesh points.

Zoom In Zoom Out Reset image size

Plasma reactors are typically run at fixed input power P_rf rather than fixed input current I_rf. To do this with our code, we used a method proposed in [8]. After every EM solve, we calculate the rf input power P_rf and compare it to the desired power P_set. Then, the applied rf current amplitude I_rf at the source is adjusted, the EM field solve is repeated and P_rf is calculated again. These steps are repeated until P_rf and P_set agree with each other within a specified relative tolerance of 4%. If there are multiple frequencies, then the above procedure is done for each of the driving frequencies.

In this section we show and discuss the results of single frequency 2D fluid-analytical simulations of the CCP reactor shown in figure 1 for a range of frequencies from 30 to 90 MHz. In each case, we assume a 500 sccm Cl₂ gas flow from the showerhead and set the pressure p at the outlet surface to 10 mTorr. For each frequency f, the input power P_rf is chosen to produce a central electron density of about $2.5\times {{10}^{16}}$ m⁻³. The 1D TL model is not self-consistent, but solves the EM equations for a given n_e and T_e. Since we chose an input n_e of $2.5\times {{10}^{16}}$ m⁻³ for the TL model, we did the same for the 2D fluid model. Since n_e depends primarily on P_e (power absorbed by the electrons), fixing n_e is essentially equivalent to fixing P_e. Both axially symmetric and antisymmetric radially-propagating wave modes exist in asymmetrically excited capacitively discharges [7]. The upper and lower axial sheath fields (E_z) are aligned (in phase) for the symmetric mode while they are opposed (180 degrees out of phase) for the antisymmetric mode.

At a frequency of f  =  30 MHz, an rf input power of ${{P}_{\text{rf}}}=1150$ W gives the base case electron density of about $2.5\times {{10}^{16}}$ m⁻³. In figure 3, we show the contour plots for (a) the electron density n_e, (b) the Cl⁻ density n₋, and (c) the electron temperature T_e. The central electronegativity $\alpha ={{n}_{-}}/{{n}_{\text{e}}}\approx 10$ for this chlorine discharge. Due to the fairly high α, n_e is fairly uniform above the wafer electrode. In a highly electronegative plasma with ${{n}_{-}}\gg {{n}_{\text{e}}}$, quasineutrality in the bulk is maintained mostly by the negative ions rather than the electrons. Since the less mobile negative ions have a much smaller temperature than the electrons, only a small ambipolar field with $|{{\phi}_{a}}(r)|\ll$ T_e is required to confine the negative ions to the plasma core. The electrons are in Boltzmann equilibrium with this small field so that ${{n}_{\text{e}}}(r)={{n}_{\text{e}0}}\exp \left({{\phi}_{a}}(r)/{{\text{T}}_{\text{e}}}\right)$, where ${{n}_{\text{e}0}}\equiv {{n}_{\text{e}}}(r=0)$. Note that we have chosen the zero of the potential to be at the radial center so that ${{\phi}_{a}}(0)=0$ and ${{\phi}_{a}}(r)<0$. Since $|{{\phi}_{a}}(r)|\ll$ T_e, ${{n}_{\text{e}}}(r)\approx {{n}_{\text{e}0}}$, and the electrons are approximately uniform in a highly electronegative plasma.

Zoom In Zoom Out Reset image size

As expected, T_e is highest near the wafer edge where the potential gradient and fields are high, leading to high radial field ohmic heating (electrostatic edge effect [8]). In figure 4, we show the results of the 30 MHz case for (a) the midplane electron density ${{n}_{\text{e}}}(r)$, (b) the midplane positive and negative ion densities n₊(r) and n₋(r), (c) the rf (black) and dc (grey) sheath voltages along the plasma-sheath boundary, and (d) the sheath width along the plasma-sheath boundary. The abscissa x in plots (c) and (d) is measured starting from the radial center of the lower wafer electrode and going around the plasma-sheath boundary until it reaches the radial center of the larger grounded electrode. So, from x  =  0 to 15 cm, the plasma-sheath boundary is along the smaller wafer electrode; from x  =  15 to 17.5 cm, it is along the dielectric spacer, and from x  =  17.5 to 52.6 cm, it is along the larger grounded electrode. The spikes prominent in the plots indicate the corners at the radial edges. The dielectric spacer region is also the source region as all the voltage drop between the powered and grounded electrodes falls across it. As expected at this low frequency, the sheath voltage and sheath width are much larger between the plasma and the smaller wafer electrode than between the plasma and the larger grounded electrode. For this configuration, the symmetric mode excitation would dominate, and we expect the the axial sheath fields (E_zs) of the wafer and grounded electrodes to be nearly in phase. The sheath voltage is fairly uniform along each electrode, indicating that at 30 MHz, the wavelengths ${{\lambda}_{\text{ps}}}$ of the symmetric mode surface waves are much larger than the system size. Figure 5(a) shows the imaginary (solid) and real (dotted) parts of the complex axial sheath field amplitude E_zs at the bottom (black) and top (grey) electrodes for the 30 MHz case. For both electrodes, the imaginary parts of E_zs point in the same axial direction and dominate over their respective real parts, indicating that at any point in time, the E_zs of the bottom and top electrodes are aligned and nearly in phase.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The wavenumber of radially propagating symmetric surface waves in the plasma of a highly asymmetric discharge is given by [5]

$$\begin{eqnarray}&&{{k}_{\text{s}}}\approx {{k}_{0}}\,{{\left(\frac{{{l}_{\text{eff}}}}{s}\right)}^{1/2}},\end{eqnarray} \tag{ 2 }$$

where ${{k}_{0}}\equiv \omega /c$ is the vacuum wavenumber, c is the speed of light, s is the sheath width of the dominant sheath, l_eff is an effective discharge width, including the skin effect, and we assume $s\ll {{l}_{\text{eff}}}$. If l is the actual discharge width and $\delta \equiv c/{{\omega}_{\text{pe}}}$ is the skin depth, then

$$\begin{eqnarray}&&{{l}_{\text{eff}}}\approx {\int}_{0}^{l}{{\text{e}}^{-z/\delta}}\text{d}z=\delta \left(1-{{\text{e}}^{-l/\delta}}\right).\end{eqnarray} \tag{ 3 }$$

Thus, at higher n_e, the skin depth δ and the effective plasma width l_eff decrease. At ${{n}_{\text{e}}}=2.5\times {{10}^{16}}$ m⁻³, ${{l}_{\text{eff}}}=2.75$ cm. From (2), if the frequency f increases or the ratio $s/{{l}_{\text{eff}}}$ decreases, then ${{\lambda}_{\text{ps}}}=2\pi /{{k}_{\text{s}}}$ also decreases, making wave effects more visible. An expected spatial resonance condition for the symmetric mode is

$$\begin{eqnarray}&&{{k}_{\text{s}}}{{R}_{\text{T}}}={{\chi}_{01}}\approx 2.405,\end{eqnarray} \tag{ 4 }$$

where ${{R}_{\text{T}}}=25$ cm is the discharge radius, and ${{\chi}_{01}}$ is the first zero of the zero-order Bessel function ${{J}_{0}}\left(\chi \right)$.

In figures 6 and 7, we show the same plots as in figures 3 and 4, but for a higher frequency of 45 MHz. At this higher frequency, a smaller drive power of 750 W gives approximately the same base case electron density with accompanying smaller sheath voltages and sheath widths. As with the 30 MHz case, the discharge is highly electronegative with $\alpha \approx 10$, and n_e is fairly uniform above the wafer electrode while ${{T}_{\text{e}}}$ is largest near the wafer electrode edge due to the electrostatic edge effect. The smaller wafer electrode dominates with sheath voltages and sheath widths much larger than that of the larger grounded electrode. The sheath voltage and sheath width are fairly uniform, but slightly center-peaked along the wafer electrode when compared to the 30 MHz case due to the reduced ${{\lambda}_{\text{ps}}}$ at the higher frequency.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At 60 MHz (as seen in figure 8), a smaller drive power of 500 W gives about the same base case electron density, but with smaller sheath voltages and sheath widths than at 45 MHz. Due to the reduction in ${{\lambda}_{\text{ps}}}$ with increasing f and decreasing sheath width s, the wave effects are larger than at 45 MHz, with the sheath voltage and sheath width noticeably center-peaked along the wafer electrode. There are also visible wave effects along the larger grounded electrode.

Zoom In Zoom Out Reset image size

In figure 9, we present results as in figure 8 at a higher frequency of 72.5 MHz. For this case, a smaller drive power of 394 W gives the same base case electron density. In contrast to the previous lower frequency cases, the smaller wafer electrode sheath no longer dominates, but instead, the peak sheath voltages and sheath widths at the radial center are nearly equal for the smaller and larger electrodes. Figure 5(b) shows the imaginary (solid) and real (dotted) parts of the complex axial sheath field amplitude E_zs at the bottom (black) and top (grey) electrodes for the 72.5 MHz case. For the bottom electrode the imaginary part of E_zs dominates and points in the negative axial direction, while for the top electrode, the real part of E_zs dominates and points in the positive axial direction. Thus, at any point in time, the axial sheath fields of the bottom and top electrodes point in opposite directions and are nearly 90 degrees out of phase, indicating the presence of a strong antisymmetric component. Also, the voltage drop across the dielectric spacer drops sharply from about 220 V at 60 MHz to about 29 V at 72.5 MHz, indicating that the antisymmetric mode is nearly spatially resonant. Note that even near an antisymmetric mode resonance, we don't have a pure antisymmetric mode, but a linear combination of both symmetric and antisymmetric modes.

Zoom In Zoom Out Reset image size

The wavenumber of an antisymmetric mode surface wave is given by [7]

$$\begin{eqnarray}&&{{k}_{a}}\approx {{\left(k_{0}^{2}-\frac{2}{{{\kappa}_{\text{p}}}\,s{{l}_{\text{eff}}}}\right)}^{1/2}}.\end{eqnarray} \tag{ 5 }$$

The first spatial resonance for the anti-symmetric mode would be at

$$\begin{eqnarray}&&{{k}_{a}}R={{\chi}_{11}}\approx 3.83,\end{eqnarray} \tag{ 6 }$$

where ${{\chi}_{11}}$ is the first zero of the first-order Bessel function ${{J}_{1}}\left(\chi \right)$. For the 72.5 MHz case, using $s\approx 0.7$ mm and ${{l}_{\text{eff}}}=2.75$ cm in (6), we obtain ${{k}_{a}}\,R\approx 4.1$.

In figure 10, we present results at 90 MHz near the first symmetric resonance, where a smaller drive power of 350 W gives approximately the same base case electron density. In this case, the wafer electrode sheath is again dominant, with significantly larger average sheath voltages than the grounded electrode, indicating a return to a primarily symmetric mode. However, a pure symmetric mode has no net current entering the exterior region, so there must still be a significant antisymmetric component, and the symmetric mode resonance condition (4) is not completely satisfied.

Zoom In Zoom Out Reset image size

At 30 MHz, the midplane electron and ion densities are fairly uniform across the wafer with the electrons slightly peaked and the ions slightly depressed at the radial center. In highly electronegative discharges with high ion densities, the recombination rate starts to exceed the ionization rate at the radial center, causing the ion densities to become depressed in the center; this in turn, increases the ion density gradients and diffusive ion fluxes towards the center until a steady-state is achieved. As discussed above, the electrons in a highly electronegative plasma are fairly uniform with a slight maximum at the center where ${{\phi}_{a}}(r)$ goes to zero. At 45 MHz, the electrons and ions behave much as in the 30 MHz case, but the slight center high of n_e is less noticeable at 45 MHz. For the same n_e and P_e, the sheath widths and voltages are smaller at 45 MHz compared to 30 MHz. So, the bulk ambipolar field required to confine the negative ions to the core is smaller for the 45 MHz case, making the center electron peak less noticeable. At 60 MHz, the behavior of the midplane ion densities are similar to the lower frequency cases, but now the electrons are slightly peaked at the radial edge of the wafer electrode rather than at the radial center. As the frequency increases to 60 MHz, the system approaches the antisymmetric mode resonance near 72.5 MHz. In this case, not only the bottom wafer electrode, but also the top electrode exhibits a large radial potential gradient (high E_r field) near the radial wafer edge at r  =  15 cm. Thus, the bulk Ohmic heating and T_e(r) become increasingly peaked at this radius, leading to an increased ionization rate and an ${{n}_{\text{e}}}(r)$ peak near the radial wafer edge. This enhancement of the edge effect reaches its maximum at the antisymmetric mode resonance near 72.5 MHz. Once the frequency increases to 90 MHz, the system approaches the symmetric mode resonance, leading to increased fields and electron T_e(r) heating at the radial center, which flattens the ${{n}_{\text{e}}}(r)$ profile.

The model nonlinearly couples azimuthally symmetric waves and series resonance harmonics in a CCP consisting of two circular conducting electrodes of radius R separated by a gap width l, as shown in figure 11. All the rf quantities are assumed to be functions of radius r and time t. The discharge is driven at r  =  R by a source having resistance R_T and voltage ${{V}_{\text{rf}}}={{V}_{\text{rf}0}}\cos \omega t$ through a large blocking capacitor C_B. The discharge is assumed to be highly asymmetric with the sheath at the bottom wafer electrode much larger than that of the upper grounded electrode so that the grounded electrode sheath can be neglected. The sheath width s is also assumed to be much smaller than the bulk plasma width d so that $d\approx l=\text{const}.$. In addition to the rf voltage, a dc voltage appears across C_B, placing a dc bias voltage $-{{V}_{\text{B}}}$ across the electrodes. The model is not self-consistent since it solves for the discharge currents, voltages and sheath widths as a function of r and t assuming a specified uniform and constant n_e and T_e.

Zoom In Zoom Out Reset image size

The nonlinearity of the charge-voltage relation generates harmonics which may be enhanced by series resonance. The series resonance between the plasma inductance and the sheath capacitance occurs at a frequency [23–26]

$$\begin{eqnarray}&&{{\omega}_{\text{SR}}}={{\omega}_{\text{pe}}}\,{{\left(\frac{s}{{{l}_{\text{eff}}}}\right)}^{1/2}}.\end{eqnarray} \tag{ 7 }$$

For a typical rf drive frequency ω, ${{\omega}_{\text{SR}}}\gg \omega$. However since the sheath width s oscillates with time, ${{\omega}_{\text{SR}}}$ sweeps across harmonics ${{\omega}_{\text{N}}}\equiv N\omega$ of the driving frequency so that

$$\begin{eqnarray}&&{{\omega}_{\text{N}}}=N\omega ={{\omega}_{\text{SR}}}={{\omega}_{\text{pe}}}\,{{\left(\frac{s}{{{l}_{\text{eff}}}}\right)}^{1/2}}.\end{eqnarray} \tag{ 8 }$$

For a given n_e (which determines ${{\omega}_{\text{pe}}}$), (8) defines the sheath widths that excite the various series-resonance enhanced harmonics of the driving frequency.

From (2), the wavenumber ${{k}_{\text{p}}}=2\pi /{{\lambda}_{\text{p}}}$ of radially-propagating surface waves in the plasma is given by

$$\begin{eqnarray}&&{{k}_{\text{p}}}=\frac{{{\omega}_{\text{SW}}}}{c}{{\left(\frac{{{l}_{\text{eff}}}}{s}\right)}^{1/2}},\end{eqnarray} \tag{ 9 }$$

where ${{\omega}_{\text{SW}}}$ is the surface wave frequency. For a voltage driven discharge (source resistance R_T much smaller than discharge impedance), we expect to observe spatial resonances at

$$\begin{eqnarray}&&{{k}_{\text{p}}}={{\chi}_{0m}}/R,\end{eqnarray} \tag{ 10 }$$

where ${{\chi}_{0m}}$ is the mth zero of the zero-order Bessel function ${{J}_{0}}\left(\chi \right)$. We expect the strongest spatial excitations to occur at the first resonance m  =  1. Equating (9) with (10), we obtain

$$\begin{eqnarray}&&{{\omega}_{\text{SW}}}=\frac{c{{\chi}_{01}}}{R}{{\left(\frac{s}{{{l}_{\text{eff}}}}\right)}^{1/2}}.\end{eqnarray} \tag{ 11 }$$

As the sheath width s oscillates with time, ${{\omega}_{\text{SW}}}$ sweeps across harmonics ${{\omega}_{M}}=M\omega$ of the driving frequency, leading to

$$\begin{eqnarray}&&{{\omega}_{M}}=M\omega ={{\omega}_{\text{SW}}}=\frac{c{{\chi}_{01}}}{R}{{\left(\frac{s}{{{l}_{\text{eff}}}}\right)}^{1/2}}.\end{eqnarray} \tag{ 12 }$$

For a given n_e (which determines l_eff), (12), defines the sheath widths that excite the various m  =  1 spatial resonance enhanced harmonics of the driving frequency.

In summary, the oscillating sheath traverses a number of values of the series resonance corresponding to harmonics of the driving frequency. These interactions lead to harmonic powers well above those that would exist in the absence of the resonances. We look for the interaction of these harmonic resonances with the spatial wave resonances that can occur at the driving frequency and at higher harmonics of the driving frequency.

The 1D nonlinear TL model is used to estimate the amplitude and phases of the harmonics that would have been generated in the 2D fluid-analytical simulations of section 2. Since the TL model assumes that the sheath of the bottom powered electrode is much larger than that of the upper grounded electrode, we are able to apply it to all the cases discussed in section 2 with the exception of the 72.5 MHz case, where there was no dominant electrode. In accordance to the configuration in section 2, we set l  =  5 cm, R  =  15 cm, p  =  10 mTorr, ${{n}_{\text{e}}}=2.5\times {{10}^{16}}$ m⁻³, ${{T}_{\text{e}}}=2.5$ V, and the electron-neutral scattering rate coefficient ${{K}_{\text{el}}}=2.24\times {{10}^{-13}}$ m³ s⁻¹ for the 1D nonlinear TL simulations. The last quantity (K_el) is obtained from [21] for a Cl₂ discharge, assuming a gas temperature of 300 K and an electron temperature of 2.5 V. In order to take the skin effect into account, we calculate l_eff from the input l and n_e, and use it in place of l in the 1D TL calculations. As with typical etching CCP reactors, we assume that the source resistance R_T is small (i.e. ${{R}_{\text{T}}}=0.5\Omega$) so that the discharge is essentially voltage-driven, and that the wafer electrode is coated with an insulating dielectric.

From the series and spatial resonance conditions, (8) and (12), we see that is important to match the sheath widths rather than the sheath voltages. So, for each applied frequency f, the magnitude of the voltage source V_rf0 in the 1D nonlinear TL simulation is chosen to keep the maximum sheath width s_m at the source (r  =  R) to be the same as that of the corresponding 2D fluid-analytical simulation in section 2. In the TL simulations, ${{s}_{\text{m}}}(R)$ is the maximum value of the time-varying sheath width s(R, t). For the 2D fluid-analytical simulations, ${{s}_{\text{m}}}(R)$ is related to the time-averaged sheath width s(R) by ${{s}_{\text{m}}}(R)=1.23s(R)$ [10, 22]. For the four cases with applied frequency f  =  30, 45, 60, 90 MHz, the corresponding ${{s}_{\text{m}}}(R)$ and V_rf0 values are ${{s}_{\text{m}}}(R)\approx 6.8$, 3.7, 1.9, 0.55 mm, and ${{V}_{\text{rf}0}}=5000$, 1400, 350, 50 V respectively.

Figure 12 shows the TL simulation results of the electron power of the N harmonics P_eN relative to the electron power of the fundamental frequency P_e1 for (a) 5000V@30MHz, (b) 1400V@45MHz,  (c)  350V@60MHz,  and  (d) 50V@90MHz. For the 5000V@30MHz case, due to the high voltage, the nonlinearities generate small harmonic powers over a broad spectrum [1]. However, at this relatively low frequency, the fundamental dominates over the higher harmonics, as expected. For both the 1400V@45MHz and 350V@60MHz cases, we see significant higher harmonics that may be near spatial resonances which lead to center-high nonuniformities. For 50V@90 MHz, the fundamental dominates over the harmonics again since at 90 MHz, the fundamental is near the first symmetric spatial resonance, leading to center-high nonuniformities, as seen in the 2D fluid-analytical simulation at 90 MHz in section 2.

Zoom In Zoom Out Reset image size