Kinetic simulation of a 50 mTorr capacitively coupled argon discharge over a range of frequencies and comparison to experiments

Advanced chip design, for applications in areas such as nano-optics, artificial intelligence, and bio-sensing, presents challenges during the processing steps for fabricating sub-nanometer features on the silicon wafer. In the plasma etching steps, which are used for the removal of materials, this process requires optimization and control over etch rate, uniformity, and surface precision.¹ These critical parameters are strongly affected by various factors including plasma generation methods, gas-phase chemistry, and surface conditions. A traditional method of plasma generation involves applying radio frequency voltages, typically at 13.56 MHz, between sets of parallel electrodes, in a configuration known as a capacitively coupled plasma (CCP) discharge.^2–4 

Plasma discharges operated in the very-high-frequency (VHF) range (30–300 MHz), however, have been shown to offer numerous benefits. VHF discharges result in a high ion current, increasing the ion yield the substrate; and a reduced sheath thickness, aiding in ion directionality^5–7 and enabling micro- and nano-scale etching with a high aspect ratio.^8,9 A high ion yield allows the operation of the discharge at a lower pressure and a lower voltage amplitude, reducing the required power input and lessening damage to the substrate.

Due to high-frequency operation, the shape of the ion energy distribution function (IEDF) impinging on the wafer surface is narrow. On the other hand, electrons that are strongly heated by the sheath oscillations penetrate into bulk plasma and drive the gas-phase plasma chemistry via collisions with the background gas. These collisions lead to a modification in the electron energy distribution function (EEDF). Published experiments and simulations have confirmed that low-pressure discharges exhibit a bi-Maxwellian electron energy distribution, where a major population of electrons has low mean energy, but a small population of electrons has very high mean energy.^6,10–18 These fast electrons assist in ionization and sustain the discharge. Since the overall population of electrons has low mean energy, the plasma can be operated with a relatively low input power.

The nature of electron energy distribution function (EEDF) is determined by the mechanism through which electrons gain and lose energy. Typically, in CCP discharges, there are three modes¹² of electron heating: collisional or Ohmic heating, collisionless heating near the sheath, and secondary electron-induced heating. The third mechanism is dominant in high-voltage discharges.¹⁹ At low driving voltages, the dominant heating mechanism varies between collisionless heating at lower pressures and Ohmic heating at higher pressures. Two interpretations of collisionless heating are known. Collisionless heating can be the result of energy transfer between the oscillating sheath edge and the electrons near the sheath edge, also called stochastic heating,⁴ or due to compression of the electron cloud near the sheath edge, also called pressure heating.²⁰ Electrons lose their energy through inelastic collisions with the neutral gas particles and absorption at the boundaries.

CCP discharges are nonlinear due to the space–time dependence of the driving electric field and the nonuniform ion density in the sheath, which makes theoretical analysis challenging. Researchers have relied on experimental and computational techniques to study discharge behavior.^{6,10–13,16,17,21–23} Godyak and Piejak¹⁰ studied the transition of the electron energy distribution from a two-temperature distribution at lower pressures to a Druyvesteyn-like distribution at higher pressures (above 1 Torr) when the driving frequency was 13.56 MHz. They concluded that this transition was caused by the change in electron heating mode from collisionless heating at lower pressures to Ohmic heating at higher pressures.

The need to optimize etch parameters and increase the yield of ions has prompted the study of the frequency dependence of the behavior of these discharges.¹ Abdel-Fattah and Sugai^11,12 conducted experiments at 50–150 mTorr at various driving frequencies for helium and argon CCP discharges, and observed a transition from a Druyvesteyn-type distribution at low frequencies to a bi-Maxwellian distribution at high frequencies (greater than 50 MHz). The transition was associated with an increase in electron number density at the center of the discharge and a decrease in effective electron temperature. As the frequency was increased, the temperature of the bulk electrons in the center of the discharge decreased, whereas the temperature of the electrons in the tail of energy distribution increased.

In this paper, we present results of particle-in-cell (PIC) simulations for the conditions studied experimentally by Abdel-Fattah and Sugai.^11,12 They determined the electron energy distribution for a wide range of frequencies, and a parametric effect has not been widely studied. Colgan et al.⁶ performed experiments on an argon discharge at 250 mTorr for various frequencies, but the electron energy distributions were not studied. Godyak and Piejak¹⁰ studied electron energy distributions for high pressures up to 3 Torr, but only at 13.56 MHz driving frequency.

Nonlinearities in capacitively coupled radio frequency (CCRF) discharges appear in the form of higher harmonic oscillations of the electric field due to production of energetic electron beams.^{16,17,21,22,24,25} In low-pressure discharges, interaction of the electron beams with moving sheaths leads to collisionless heating. Above a critical driving frequency, the electron beams are effectively confined in the plasma column, resulting in enhanced electron heating and a drastic increase in electron density.^16,21,22,24 The process of confinement is explained further in Sec. III C. At higher pressures, Wilczek et al.²⁴ found that there was no drastic increase in the electron density and that the variation of the electron density with driving frequency was closer to the theoretical estimate^2,26 of $n e ∝ f 2$⁠.

According to Wilczek et al.,²⁴ depending on the driving frequency, the electron and ion densities increase by one or two orders of magnitude as pressure increases from 1.3 Pa (about 9.7 mTorr) to 5 Pa (about 37 mTorr). The rapid increase in electron density restricts computational parameters like cell size and time step, which increases time taken for a single simulation. Due to this reason, many PIC simulations^{14–17,21,23–25} have been conducted at lower pressures (below 10 mTorr). Since higher pressure results in high ion yield for materials processing applications, the current study focuses on conducting simulations for a 50 mTorr argon plasma at various frequencies and comparing the results with published experimental data.^11,12 In addition to the frequency dependence of the discharge behavior, we also investigate the influence of nonlinearity on electron heating and electron density.

The paper is organized as follows. In Sec. II, we discuss the numerical procedures used to solve the governing equations (the Boltzmann equation for the species and the Poisson equation for the electric potential) along with the simulation parameters used for this study. In the first part of Sec. III, we discuss the electron energy probability functions obtained from the simulations at various driving frequencies. We calculate the electron densities and the effective electron temperatures, and compare the distributions and the aforementioned calculated plasma properties to results from published experimental data. In the second part of Sec. III, we use a simple homogeneous lumped-parameter model to determine the sheath thickness. In the last part of Sec. III, the power deposition by the field to the electrons is analyzed and all the features observed are discussed in detail. Finally, a summary of the results is provided.

The simulations are carried out using a modified version of the EDIPIC (electrostatic direct-implicit particle-in-cell) code²⁷ developed at the Princeton Plasma Physics Laboratory. The code simulates one dimension in physical space and three dimensions in velocity space (1D-3V). This code has been benchmarked for various cases and shown to be accurate in predicting the properties of a discharge over a wide range of pressures.^22,28–32

In order to carry out a parametric study of the influence of driving frequency, the conditions for the simulations in this study were chosen to match those of Abdel-Fattah and Sugai's published experiments.¹² The simulation parameters are summarized in Table I, and a schematic diagram of the computational domain is provided in Fig. 1. An argon plasma is operated using an 80 V (peak-to-peak) sinusoidal voltage [ $V R F = 40 sin ( 2 π f t )$] applied to the left electrode at driving frequencies between 13.56 and 73 MHz. The right electrode is grounded. The gap between the electrodes is 6.5 cm, and the neutral gas pressure is 50 mTorr (about 6.7 Pa). The associated neutral gas density is $1.61 × 10 21 m − 3$⁠, which is uniformly distributed throughout the simulation domain. The gas temperature is 300 K. Since the maximum magnitude of potential difference applied between the electrodes is 40 V, the boundary conditions for the walls are set to be 100% absorbing. Hence, the values of secondary emission coefficient and the effective reflection coefficient are taken to be 0.

For the calculations, the number of macroparticles per cell is chosen to be 200. The gap between the electrodes is divided into 3911 equal sized cells, $16.62 μ m$ wide. The cell size is chosen to satisfy the limit imposed by the Debye length (⁠ $λ D ≈ 300 μ m$⁠). An additional study with 1566 cells in the domain, each of width $41.5 μ m$⁠, was done to study the effect of grid resolution. Changing the cell size did not affect the results at low driving frequencies (like 13.56 MHz). At high driving frequencies, the larger gradients in electric field and current near the sheath were better resolved when the finer grid was used. Figure 2 shows the time average number density profiles for the 73 MHz driving frequency case for the two different grid sizes and it can be seen that the finer grid resolves the profile of density in the quasineutral region, whereas the coarse grid results in a lower density in the center.

Figure 3 shows the profiles of electron current at the beginning of the input voltage cycle (⁠ $ω t = 0$ according to the defined voltage) for the two grid sizes. The coarser grid produces a non-uniform profile of the total current near the sheath edge (below 0.325 cm) because the gradient in electron conduction current near the sheath edge is not resolved. A cell size of $16.62 μ m$ was used for the remaining results presented in this article.

The calculations were carried out on computing resources provided by the Purdue University Community Cluster Program.³⁴ Initially, a case was run for thousands of RF cycles to allow it to reach a periodic state. Convergence was confirmed by monitoring the time variation of electron and ion densities at the center of the domain after the initial run. In a periodic state, the value of electron and ion densities at the center of the domain should be constant, as shown in Fig. 4. After a periodic state was achieved, the simulation was run for an additional 100 RF cycles to obtain statistical data. The time steps were chosen to be smaller than period of plasma oscillations divided by $2 π$ (⁠ $1 / ω p ≈ 100 ps )$⁠.

Figure 5 shows the calculated electron energy probability functions (EEPF) at the center of the discharge gap for each case. (See Table I for the list of driving frequencies.) The data are normalized such that the EEPF, when multiplied by the square root of the energy and integrated over energy, gives the electron number density. Plotted in this form, the shape of the EEPF for each case exhibits three regions of roughly constant slope, one at very low electron energies, one in the mid-energy range, and finally one at high electron energies. Non-Maxwellian energy distributions in capacitively coupled discharges form due to nonlocal kinetic effects, a phenomenon that has been studied previously for low background gas pressure.^35–38 

As frequency increases, the shape of the distribution evolves, and the curve shifts upward because of an increase in electron number density. For the cases from 13.56 to 44 MHz, Fig. 5(a), the population of electrons in the low-energy region increases much faster than the population in the higher energy range. The production of low-energy electrons in ionization processes results in an increase in the population of the bulk electrons, while the rate of energy deposition into the electrons is low for low driving frequency cases, resulting in a disproportionate increase in the population of electrons in the tail of the EEPF (further analysis of energy deposition is provided in Sec. III C). For the cases from 50 to 73 MHz, Fig. 5(b), the population of both the bulk and the tail portions of the EEPF increases proportionally as the rate of energy deposition at the center becomes more efficient at high driving frequencies.

Figure 6 compares the distribution functions obtained from the present simulations to the corresponding experimental data.¹² For the purpose of comparison of simulation results with the experimental results, the data in Figs. 6 and 7 are plotted in the units of $cm − 3$⁠.

We see from Fig. 6 that the computed values for EEPF generally follow qualitatively with the experimental measurements in the range of 2–10 eV for all cases. The agreement is found to be in the mid-energy slope out of the three slopes observed in the numerical calculations of the EEPF. In comparing the numerical calculations to the measurements, it must be kept in mind that the probes used in the experiments only provide accurate data for a limited electron energy range.¹¹ The experimental data are not accurate for very low and very high energies. Furthermore, the experimental data in Ref. 12 were extrapolated to zero for the low energy values, which results in an unphysical local maximum at low energy for the experimental distribution. As a result, at low driving frequencies, the experimental distributions do not exhibit the first constant slope region at low electron energies, which are observed in the simulated distribution functions. At high electron energies (greater than about 12–15 eV), the experimental distributions fall steeply, whereas the computational results show a third region of constant slope. These experimental limitations also affect the normalization of the distribution, introducing a shift in the experimental EEPF and uncertainty in the electron number density derived from its integral.

Figure 7 shows the comparison of the electron density and effective electron temperature, calculated from the EEPF using Eqs. (2) and (3), between simulation and experiment.¹² The electron number density, Fig. 7(a), increases with frequency, and the predicted rate of increase roughly matches that of the experiments. The magnitude of the predicted electron number density, however, is significantly lower than the reported experimental values. As mentioned earlier, limitations in the experimental energy resolution introduce a shift in the experimental EEPF and uncertainty in the electron number density derived from its integral. From Fig. 6, we see that there is an additional area under the EEPF curve in the experimental data, which does not exist in the simulation results. This uncertainty may, in part, explain the level of discrepancy.

In the calculations reported in this paper, both secondary emission and electron reflection at the electrodes were disregarded due to the relatively low pressure and voltage. Surface effects at the electrodes were found by Schulenberg et al.³⁹ to play an important role in contributing to the electron density at relatively high applied voltages (150–350 V). Those authors used a trial-and-error method to find that a reflection coefficient of 0.7, and a secondary emission coefficient of 0.07, provided a best match for their calculated electron densities to their measured values. To ensure that the discrepancy in electron density between the present simulations and Ref. 11 was not due to surface effects, a similar test for our low-voltage discharge cases was conducted. It was found that the electron density increased by only 26%, insufficient to explain the difference between the simulation and experimental results.

According to a simple theoretical analysis,^2,26 which ignores energy gained by electrons from collisionless heating, electron density is predicted to increase as $n e ∝ f 2$⁠. The present computations predict that the electron density increases as $n e ∝ f 2.05$ as shown by the dashed curve in Fig. 7(a), similar to the theoretical trend. For the present case, the gas pressure is not low enough for collisionless heating to be a dominant mode of power deposition. This can be further confirmed by calculating the ratio of power deposited by stochastic heating (P_stoc) to that by Ohmic heating (P_Ohm) using the process outlined in Abdel-Fattah and Sugai.¹² The values for the ratio $P stoc / P Ohm$ were found to be less than 1 for all the driving frequency cases studied here. Even for the highest driving frequency studied (73 MHz), $P stoc / P Ohm$ was found to be 0.959.

The trend of experimentally calculated electron density vs the driving frequency shows a change in the slope of the curve after a driving frequency of 27 MHz; this behavior is absent in the computational curve. It is to be noted that Wilczek et al.²⁴ studied the effect of background gas pressure on the dependence of bulk electron density vs the driving frequency and found that at pressures above 2 Pa (about 15 mTorr) the drastic increase in electron density after certain driving frequency disappeared, and the trend became closer to the theoretical approximation of $n e ∼ f 2$⁠, for which non-local kinetic effects are ignored.^2,26

Figure 8 shows the profiles of electron and ion densities for the lowest (13.56 MHz) and the highest (73 MHz) driving frequency cases. As the driving frequency increases from 13.56 to 73 MHz, the density at the center also increases from $5 × 10 14$ to $1.5 × 10 16 m − 3$⁠, respectively, as was previously observed, $5 × 10 8$ to $1.5 × 10 10 cm − 3$ in Fig. 7(a). Ionization increases at higher driving frequencies and causes the plasma to fill up the discharge gap. Consequently, the thickness of the combined sheath also decreases from 10 to 3.3 mm as the driving frequency increases from 13.56 to 73 MHz. The thickness of the sheath for all the driving frequencies is calculated from the simulation parameters a simple homogeneous model in Sec. III B.

In the VHF range, CCP discharges can be modeled by a simple lumped-parameter circuit,⁴ since the wavelength of the input voltage signal is much larger than all geometric scales involved in the discharge. In this homogeneous model, a uniform quasineutral plasma is sandwiched between two sheaths completely devoid of electrons. The ions are assumed to be stationary, with constant density throughout the discharge gap. In contrast, the electrons are assumed to respond instantaneously to the electric field and to oscillate at the driving frequency between the electrodes.

The two sheaths are represented by variable gaps of width $d 1$ and $d 2$⁠, as shown in Fig. 9(a). The plasma column of thickness $L − ( d 1 + d 2 )$⁠, where L is the gap between the electrodes, oscillates with the RF frequency, just touching the electrodes.

The equivalent circuit diagram for the discharge gap is shown in Fig. 9(b). The two sheaths are combined into a single equivalent vacuum capacitance $C sh$ with constant total size $d sh = d 1 + d 2$⁠. The plasma column is modeled by a combination of a resistance $R pl$ and a capacitance $C pl$ connected in parallel. A detailed explanation of this model is given in Ref. 4.

Figure 10 shows the comparison of the magnitude and phase of the impedance obtained with Eqs. (7) and (8) for all the values of driving frequencies. The magnitudes match very well for all frequencies, while the phases have a maximum deviation of about 7%. A sensitivity analysis was conducted using various ratios of electron to ion densities ranging from 1% to 30% and the corresponding sheath edge locations to determine the best fit for the left- and right-hand sides of Eq. (8), and it was found that defining the sheath edge with $n e / n i = 10 %$ gives the least deviation. Simple estimates for the combined sheath thickness and the electron density thus give accurate representation of the RF discharge as an electrical circuit.

The scaling of sheath thickness with respect to the driving frequency depends on the operating conditions, like driving source (voltage/current/power driven), and the electrode gap. Using a simple model comparing the ion drift flux to the ambipolar diffusion flux, Raizer et al.⁴ predicted the dependence of sheath thickness on driving frequency as $d sh ∝ f − 2 / 3$⁠. Vahedi et al.⁷ used an electron power balance to predict that the sheath thickness scales approximately inversely with the frequency and confirmed this using two-dimensional PIC simulations for an argon discharge driven by a constant voltage at 60 mTorr. Experimental studies have been done at various driving frequencies on relatively high-pressure (around 1 Torr) air discharges driven by a constant power,⁴⁰ and by a constant voltage amplitude,⁴¹ and have found the dependence of the sheath thickness calculated using the lumped-parameter model to be $d sh ∝ f − 0.43$ when driven by constant power, and $d sh ∝ f − 2 / 3$ when driven by constant voltage. In our case of a voltage-driven discharge, when the estimated combined sheath thickness is plotted against driving frequency (Fig. 11), we see its dependence as $d sh ∝ f − 0.66$⁠, which is very close to the estimates of $d sh ∝ f − 2 / 3$ predicted by the theory.⁴ 

Defining the sheath edge at the location where the electron density is 10% of the ion density gives a good estimate of the thickness of the sheath for the equivalent circuit model. Another way to define the sheath thickness is comparing the electric field at the sheath edge to the electric field at the electrode. In our study, we found that at the value of sheath thickness used for the calculations in Eq. (6), the ratio of the electric field at that location to the field at the electrode remains constant (about 0.4) irrespective of the driving frequency, as shown in Fig. 12.

While the lumped-parameter model gives a simple way to estimate sheath thickness, it does not take into account more complex phenomena happening at the sheath edge. The spatial variation of ion density and the production of electron beams in the sheath modulate the sheath edge such that it oscillates with higher harmonics in addition to the driving frequency.^{14,17,19,42,43} Investigating these nonlinear effects will aid in understanding the discharge dynamics.

The power deposition to the electrons was computed by time averaging $E · j e$ over the last 100 RF cycles after the simulations reached a steady state. Deposited power is transferred to the kinetic energy of electron motion, which in turn is randomized in part in collisions, leading to increased electron temperature. Negative values of the power deposition indicate that the electrons lose kinetic energy to the field.

Figure 13 shows the profiles of average power deposition $⟨ E · j e ⟩$ for different driving frequencies. Electrons get heated in vicinity of the sheath; this energy deposition increases as the frequency increases as shown in Fig. 13(a). Since electron density in the sheath is negligible compared to that in the bulk plasma, the electrons gain energy on average in a collisionless manner through interaction with the oscillating, highly nonlinear electric field rather than through Ohmic heating. In such a case, the energy gained depends on the velocity of oscillation of the plasma-sheath boundary. At high driving frequencies, the plasma-sheath boundary oscillates much faster compared to the low driving frequency cases. Hence, we have higher power deposition at high driving frequencies. Later in this section, we will see that there are higher harmonic oscillations present in the sheath, in addition to the fundamental driving frequency. These nonlinearities enhance power deposition by the field into the electrons for the high driving frequency cases.

As the driving frequency increases, the amplitude of electron oscillation decreases. Hence, at higher driving frequencies, the electrons lose their energy in a narrower region near the sheath edge, compared to the lower driving frequency cases, as seen in Fig. 13(a). In the center of the discharge, $⟨ E · j e ⟩$ is relatively small due to the lower values of electric field there. The power deposition in the center of the discharge gap increases as frequency increases, as seen in Fig. 13(b), where selected cases from the region denoted by the red box in Fig. 13(a) are plotted. This is due to the increase in center electron density as the frequency increases, giving rise to a larger electron current.

Careful consideration is required to assess the power deposition by the field to the electrons near the sheath edge. As we observe in Fig. 13(c), where selected cases from the region denoted by the blue box in Fig. 13(a) are plotted, just outside the sheath edge electrons lose energy on average to the field, an effect that increases in magnitude as the frequency increases. Closer to the center from this region of negative $⟨ E · j e ⟩$⁠, there is a local maximum in power deposition for the high driving frequency data shown by the green arrows in Fig. 13.

CCP discharges operated by both sinusoidal and non-sinusoidal waveforms have shown to produce fast electron beams during the sheath expansion phase.^{16,21,24,25,30,44,45} At low gas pressures, these electron beams travel uninhibited to the opposite sheath. If the opposite sheath is in the collapsed phase when the fast electron beams reach it, the electrons are lost to the opposite electrode and are not confined in the plasma, as was shown by Wilczek et al.²⁴ On the other hand, if the electron beams reach the opposite sheath during its expanded phase, they are reflected back into the plasma and cause additional ionization in the bulk. Researchers have observed that at low pressures, there is a critical driving frequency above which the electron beams are successfully confined in the plasma and cause a step increase in the electron density at the center of the plasma.

As mentioned before, Wilczek et al.²⁴ also studied the effect of background gas pressure on the dependence of bulk electron density vs the driving frequency and found that at pressures above 2 Pa (about 15 mTorr), the drastic increase in electron density disappeared, and the trend followed was closer to the theoretical approximation of $n e ∼ f 2$⁠, for which non-local kinetic effects are ignored.^2,26

Figures 14(a) and 14(b) show the different components of current for 13.56 and 73 MHz, respectively. Here, the currents are conditionally averaged based on the phase of the voltage input. They are plotted at four different phases of an RF cycle, specifically $ω t = 0 , 0.5 π , 1.0 π$⁠, and $1.5 π$⁠, based on the sinusoidal voltage waveform defined in Sec. II.

At 50 mTorr gas pressure, we see the effect of fast electron beams near the sheath edge as a local increase in electron conduction current (⁠ $j e$⁠), about 71% greater than the total current (⁠ $j t$⁠) value, as shown in Fig. 14(b) at $ω t = 0.5 π$ and $1.5 π$⁠. Since the total current (conduction plus displacement) should remain constant over the entire domain, the increase in $j e$ is compensated by an equivalent change in displacement current, $j d$⁠. We see from Fig. 14(a) that the local increase in $j e$ is absent at low driving frequencies and is only observed as the frequency increases. Since the higher gas pressure hinders uninhibited travel of the electron beams, they are unable to reach the opposite sheaths. Instead, local heating is observed up to short distance after the sheath edge, about 8–16 mm as seen in Fig. 13(c), before the electrons lose their energy in inelastic collisions.

As the fast electron beams exiting the sheath travel through the plasma, colliding with the neutral gas, they encounter a region of opposing electric field just after the sheath edge. Figure 15 shows the conditionally averaged (based on the phase of the voltage input) electric field near the left (powered) electrode at the same phases in the RF cycle as the plots of current in Fig. 14. The direction of travel of the plasma beams can be inferred from the sign of $j e$ at each phase. At $ω t = 0.5 π$⁠, the sheath is in the collapsed phase, top right plot of Fig. 14(b), and the increase in local conduction current is positive, which results from the fast electrons moving toward the electrode. As seen in Fig. 15(b), the electric field at this moment in the same region is negative, which results in an electric force on the electrons in the direction away from the electrode. Since the motion of beam electrons and the force on them are in opposite directions, the electron decelerates and loses energy to the field. Similarly, at $ω t = 1.5 π$⁠, the sheath is in the expanded phase, bottom right plot of Fig. 14(b), and the increase in local conduction current is negative, which is again opposite to the electric force on the electrons, with positive field in Fig. 15(b). Averaging this over the entire RF cycle results in a negative $⟨ E · j e ⟩$ just outside the sheath edge. This explanation is summarized in Table II. At 13.56 MHz driving frequency, we do not observe the production of fast electron beams near the sheath edge, since there is no local increase in the electron conduction current as seen in Fig. 14(a).

Nonlinearities play an important role in the generation of fast electrons, electron cooling near the sheath edge, the eventual local heating after the sheath edge, and the broad positive power deposition plateau around the center of the discharge gap. Nonlinear dynamics in the discharge show up in the form of various high-frequency oscillations of properties like the electric field and electron current. Figures 16(a) and 16(b) show the time variation of the electric field in one RF cycle of the 73 MHz case at the location of minimum $⟨ E · j e ⟩$ (around 4 mm) and the electric field at that location in the frequency domain. Figures 16(c) and 16(d) show the same for the electron current, and Fig. 16(e) shows the variation of $E · j e$ over one RF cycle at the same location. The vertical axis in Figs. 16(b) and 16(d) is normalized such that the energy of all the oscillating components adds up to one. The phase difference between the electric field and the electron current is clearly visible in Figs. 16(a) and 16(c). From Figs. 16(b) and 16(d), we see that along with oscillation at the fundamental frequency (74.2% of total wave energy), the second (16.6%), third (4.2%), fourth (1.7%), and the sixth (1.7%) harmonics are also present in the electric field, while the electron current oscillates at the second (70.9% of total wave energy), third (13%), and the sixth (0.4%) harmonics, along with the fundamental frequency (12.6%).

In general, as we move from the sheath edge to the center of the plasma, these nonlinearities interact to either enhance or reduce the different components of oscillation. Figures 17(a) and 17(b) show that at both low and high driving frequencies, the second harmonic oscillation of the electric field is much more significant than the fundamental frequency at the center of the discharge. We have only observed oscillation in the electric field at a few harmonics higher than the fundamental frequency for our 50 mTorr case. The high number of collisions in the discharge gap suppresses very high-frequency oscillation. Research has demonstrated that low-pressure discharges (less than 10 mTorr) display harmonics at very high frequencies (larger than 10 or 20 times the fundamental frequency) at the center of the gap.^15–17,25

As an example of the role these harmonics play at higher driving frequencies, the electric field and the current were low-pass filtered so that only the fundamental frequency was retained and $⟨ E · j e ⟩$ was recalculated from the filtered values. Figure 18 shows the comparison of $⟨ E · j e ⟩$ at the location of minimum $⟨ E · j e ⟩$ and at the center of the discharge gap using unfiltered and filtered field and current. We see from Fig. 18(a) that the electrons lose energy on average only when the driving frequency is 50 MHz and higher. At frequencies below 50 MHz, filtering the current and the field does not have much effect on $⟨ E · j e ⟩$⁠, whereas above 50 MHz, including oscillations only at the fundamental frequency reduces the magnitude of electron cooling.

Nonlinearities contribute to local electron heating in the center of the discharge gap as well. From Fig. 18(b), we see that at higher driving frequencies, the magnitude of electron heating at the center is lower when the filtered values of current and electric field were used.

The first term on the right-hand-side of Eq. (11) corresponds to the field due to temporal changes in the electron current and is attributed to electron inertia. The third term is also an inertial term (the convective acceleration), which is attributed to the spatial changes in electron density. The second term accounts for collisions, and the fourth term corresponds to diffusion.

This is a simplified expression obtained under a number of assumptions. It is valid only in the quasineutral region of the plasma where electron conduction current is dominant over the displacement current. It does not account for the oscillation of the sheath. The electron density profile is assumed to be a step-like function, which increases from zero in the sheath to the plasma density at the sheath edge. Ionization in the sheath is neglected by assuming sheath width much smaller than the bulk width, which eliminates the source term in the equation for the conservation of mass of electrons. When this simple analytical expression is plotted using current and electron density from the simulation and compared to the electric field obtained directly from PIC (Fig. 19), they match very well in the bulk of the plasma and near the sheath edge, but both of the curves do not show signs of field reversal as was found in collisionless discharges. This was true for all the driving frequencies studied.

Using Eq. (11) as the representation of the electric field in the discharge, we can analyze the contribution of the individual terms in the equation to the average power deposition. Figure 20 shows zoomed in profiles of the individual terms of Eq. (11) multiplied by the electron current, which correspond to power deposition (a) from electron inertia due to instantaneous changes in the electron current (which we will call inertia term 1), (b) owing to the effect of collisions (collision term), (c) from inertia due to changes in the electron density profile (inertia term 2), and (d) owing to the effect of diffusion of electrons (diffusion term). The red and orange curves correspond to the same phase in the RF cycle as the red and orange curves in Fig. 15(b).

Past the sheath edge, the contributions of inertia term 2 and the collision term are negligible. During both halves of the cycle, irrespective of the potential at the powered electrode, inertia term 1 contributes negative values to the average power deposition, as seen in Fig. 20(a). In a perfectly sinusoidal system, averaging the term $∂ j ∂ t · j$ over an RF cycle should equate to zero. In our case of a highly nonlinear system, the contribution from inertia term 1 to the average power deposition remains negative. While the diffusion term also significantly contributes to $E · j e$ as seen in Fig. 20(d), it oscillates between positive and negative values during different phases of the cycle. Therefore, electron inertia, as a result of fast electron beams disrupting the linearity of the discharge, is responsible for electron cooling near the sheath edge.

Thus, simulations and analysis show that the negative heating of electrons near the sheath edge at the conditions considered here is not associated with field reversal. Instead, it is the result of electron inertia in rapidly oscillating fields combined with the discharge nonlinearity.

A computational study using the PIC methods was carried out of CCP discharges in argon at 50 mTorr and driving frequencies from 13.56 to 73 MHz. The energy distributions at the center of the plasma were complex and non-Maxwellian for all driving frequencies. When compared to published experimental results, the distributions matched reasonably well in a mid-energy range, 2–10 eV. Outside that range, the experiment data may be unreliable due to the resolution limitations of the probe used in the experiments. The numerical calculations matched the experimental trend of increasing electron density and decreasing electron temperature at the center of the discharge.

In the context of a lumped-parameter circuit model, the calculated impedance of the discharge matched well with the same parameter calculated through the voltage and current from the simulation. The estimates of the sheath thickness used to obtained the impedance values at all the frequencies followed the power law trend of $d sh ∝ f − 0.66$⁠. The ratio of electron to ion number densities at the best estimates of the sheath edge location was 10%, and the ratio of the electric field at the estimated sheath edge to the maximum electric field revealed a constant value of 0.4. These two ratios can be used to approximate the sheath thickness for the equivalent circuit model.

Interesting features were observed in profiles of averaged power deposition to the electrons $⟨ E · j e ⟩$⁠. Power deposition increases with increasing frequency in the sheath due to enhanced heating of electrons at high driving frequencies. Electron heating at the center of the discharge increases with increasing driving frequency due to the increased conductivity of the plasma at higher frequencies, aided by higher harmonic oscillation of the electric field. Near the sheath edge, power deposition first decreases, becomes negative, and then increases above zero, before finally reducing to the value at the center. The negative $⟨ E · j e ⟩$ past the sheath edge is due to the production of fast electron beams in the sheath, which travel in a direction against the electrical force on the electrons giving rise to electron cooling there. After this region of cooling, the electron beams aid in enhanced heating before losing their energy in inelastic collisions. The production of fast electron beams is due to the nonlinear dynamics of the sheath. The fast oscillations in the sheath, and consequently the field, energize the electrons on a timescale much smaller than the time required for the electrons to react, adding to the electron heating in the sheath, increasing electron cooling near the sheath edge, and enhancing the electron heating in the bulk of the plasma.

A simple theoretical approach derived from the momentum equation was used to analyze the contributions from various plasma phenomena to the behavior of electron heating. The analysis showed that electron inertia plays a significant role, resulting, in combination with harmonic generation, in an overall average negative power deposition near the sheath edge.

This work was conducted as a collaborative research project at the Princeton Collaborative Research Facility (PCRF), which is supported by the U.S. Department of Energy (DOE) under Contract No. DE AC02-09CH11466. The work by Purdue University co-authors (S.S., J.P., and S.M.) was supported by the U.S. DOE Fusion Energy Sciences Office under the grant DE-SC0021076 (Dr. Nirmol Podder, Program Manager). We express our gratitude to Purdue University's Rosen Center for Advanced Computing (RCAC) for computer resources to run the simulations.

The authors have no conflicts to disclose.

Saurabh Simha: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Visualization (lead); Writing – original draft (lead). Sarveshwar Sharma: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal). Alexander V. Khrabrov: Conceptualization (supporting); Data curation (supporting); Investigation (supporting); Methodology (equal); Software (equal); Validation (equal). Igor D. Kaganovich: Funding acquisition (lead); Software (lead); Supervision (supporting); Writing – review & editing (supporting). Jonathan Poggie: Conceptualization (supporting); Data curation (supporting); Formal analysis (equal); Investigation (equal); Project administration (supporting); Resources (supporting); Supervision (equal); Writing – review & editing (equal). Sergey Macheret: Conceptualization (equal); Data curation (supporting); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Project administration (supporting); Supervision (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.