Effect of higher driving frequencies in dual-frequency discharge on plasma generation in capacitive coupled plasmas: PIC-MCC simulation

The semiconductor industry has been advancing rapidly, driven by continuous miniaturization, three-dimensional (3D) integration, and increasing complexity in transistor design.^1,2 With these developments, the demand for more advanced plasma processing techniques has become crucial. As the industry pushes the limits of device scaling, traditional techniques are no longer sufficient to meet the precision and efficiency required for modern fabrication. Consequently, there is a need for new methodologies that can precisely control the interactions between various plasma species, such as ions, electrons, and radicals, within these increasingly intricate processes.³ 

In particular, the control of ions, electrons, and radicals in plasma plays a vital role in semiconductor fabrication, especially in processes like Plasma-Enhanced Chemical Vapor Deposition (PECVD).^4–6 PECVD has been widely used for thin-film deposition, where it is essential to control not only the plasma density but also the energy distribution of the ions and electrons, as well as the flux of radicals. Several discharge techniques have been proposed to address these needs, with many focusing on advanced methods for independently controlling these plasma components. To investigate the independent control of plasma properties, many studies have used various discharge techniques such as the dual-frequency (DF) excitation,^7–19 the electrical asymmetry effect (EAE),^20,21 the tailored voltage waveforms (TVWs)^22–25 discharge method, and the amplitude modulation method.^26–31 

Using dual- or multi-frequency excitation has opened up new possibilities for independently controlling different plasma properties. In general, high-frequency (f_H) sources are used to control plasma density and the resulting ion flux, while low-frequency (f_L) sources are employed to control the ion energy upon impact with the substrate.^7–19 Therefore, f_L is typically set below the ion plasma frequency. The ability to decouple the ion flux from the ion energy is particularly significant for applications in which uniform deposition and etching are crucial for maintaining device integrity.

In research focused on high-frequency excitation in single-frequency f_rf discharges (SFD), several studies, such as those by Sharma et al.,³² have explored the effects of increasing the driving frequency. Higher driving frequencies have been found to impact various plasma parameters, including the ion density and radical generation efficiency, which directly influence deposition rates and film properties.^33–35 However, the research on dual-frequency discharges (DFD), particularly where the lower frequency (f_L) is set above the ion plasma frequency, is still limited. The combination of a higher frequency f_H with f_L above the ion plasma frequency (f_L > f_pi) presents unique challenges, especially concerning the impact on electron heating effects. It is not yet fully understood how this setup influences electron dynamics and, consequently, the overall plasma characteristics.

Particle-In-Cell Monte Carlo Collision (PIC-MCC) simulations have become a valuable tool for addressing these uncertainties.^36–41 Unlike experimental methods, which often face difficulties in measuring ion and electron behaviors at fine scales, PIC-MCC simulations allow for detailed analysis of plasma dynamics, providing insights into how different parameters affect plasma behavior and interactions. The ability to model and visualize the motion of charged particles under various discharge conditions makes this approach indispensable for investigating the effects of dual-frequency excitations in plasma systems.

In light of these considerations, this study aims to investigate the electron heating mechanisms in a DFD configuration, specifically focusing on the effects of a higher driving frequency f_H and f_L = 13.56 MHz in capacitively coupled plasma (CCP) using the PIC-MCC model. We explore the role of electric field fluctuations in this setup and analyze how they contribute to the electron heating process.

The simulation was conducted using the PIC-MCC model code, PEGASUS.⁴¹ The PIC-MCC model^36–40 combines the particle-in-cell (PIC) method with the Monte Carlo collision (MCC) technique, treating ions and electrons as discrete particles to accurately compute the energy distribution of charged particles. This model is capable of simulating not only charged particles in electromagnetic fields but also low-temperature plasmas within plasma reactors. Moreover, to evaluate the density of the neutral species, direct Monte Carlo (DSMC) simulation was used. The simulations were conducted using PIC-MCC and DSMC methods (PEGASUS Software Inc.⁴¹) to analyze charged particles and neutral species, respectively. Both PIC-MCC and DSMC are computational techniques that derive velocity distribution functions by solving the Boltzmann equation. These simulations allow for the analysis of the spatial distribution of the electromagnetic field, as well as the movement of ions, electrons, and neutral species within the plasma. PIC-MCC was employed to analyze plasma density, temperature, ion/electron energy, ion sheath expansion, and the spatiotemporal behavior of the electric field in CCP. PIC-MCC computes the motion of ions and electrons in each cell by solving the equations of motion and Maxwell’s equations over time intervals. The effects of collisions are included in the calculations through the incorporation of the Monte Carlo Collision simulation. In addition, DSMC simulates the behavior of radicals within individual cells throughout the computational domain. Both the PIC-MCC and DSMC methods use superparticles, representing groups of several billion individual particles, in place of modeling each particle individually. For the DSMC simulation, all collisions are assumed to be elastic, with no external forces acting on the particles. While numerous studies have modeled the behavior of charged particles using PIC or PIC-MCC, few have simultaneously calculated the behavior of both charged particles and neutral species. In this work, the working gas is Ar and species and gas-phase reactions involved in argon plasma simulation are presented in Table I. The reactions considered include three types of collisions for electron–Ar interactions (elastic scattering, excitation, and ionization) and two types for Ar⁺–Ar interactions (backward scattering and isotropic scattering).^42,43 In addition, stepwise and pooling ionization related to metastable states (Ar*) are also taken into account.^44,45

In this study, we investigate the effects of higher driving frequency in DFD on plasma parameters and higher harmonic generation in capacitively coupled plasmas, analyzing the behavior of ions, electrons, and radicals through a combination of PIC-MCC and DSMC methods.

Figure 1 illustrates a schematic of the parallel-plate CCP reactor used in this study. The reactor comprises a bottom electrode powered by RF and a grounded (GND) top electrode, with a gap of 50 mm between them. The z = 0 mm position was defined on the surface of the bottom RF-powered electrode. We sustained a discharge by applying 200 peak-to-peak voltage V_pp of dual-frequency discharge f_H + f_L and single-frequency discharge to a powered electrode. For dual-frequency discharges, the waveforms of V_pp are expressed by the following formula: $Vppt=Vrfcos2πfHt+cos2πfLt$⁠. On the other hand, for single-frequency discharges, the waveforms were given as $Vppt=Vrf⁡cos2πfrft$⁠. Here, f_L = 13.56 MHz was fixed although f_H and f_rf was varied from 27.12 MHz (=2 f_L) to 108.48 MHz (=8 f_L), which are integer multiples of f_L. The value of V_rf was fixed at 200 V, and the magnitudes of V_pp were adjusted to be equal for both the single-frequency discharge and the dual-frequency discharge. Furthermore, V_dc was set to 0 V to examine the impact of the discharge waveform itself on plasma parameters.

The left and right boundaries were defined as symmetrical boundaries. The calculation was conducted assuming an infinitely extended parallel plate. The coordinate system employed was a two-dimensional x–y (Cartesian) system. A typical value used for the ion-induced secondary electron emission coefficient (SEEC) is γ = 0.1. Gas pressure = 10 mTorr and a gas temperature of 300 K. The simulation region (50 mm electrode gap) is divided into 250 grids, with 200 superparticles per cell for all simulation sets. The time step is set to 5 × 10⁻¹¹ s.

Figure 2 shows the driving frequency dependences of (a) spatio-averaged electron density and (b) power density in single- and dual-frequency discharges. Here, f_L = 13.56 MHz in DFD, and the discharge voltage in both SFD and DFD was fixed at V_pp = 200  V. The values of the spatio-averaged electron densities among the electrodes in Fig. 2(a) were averaged over f_rf in SFD and f_L in DFD. In addition, in this calculation system, the Ar ion density and electron density are equal. The spatio-averaged electron density increases with increasing driving frequency in both single- and dual-frequency discharges. However, the electron density in SFD is larger than in DFD across all frequency ranges. For example, the electron density at f_rf = 27.12 MHz in SFD is 1.1 times larger than at f_H + f_L = 27.12 + 13.56 MHz in DFD. Furthermore, the electron density at f_rf = 108 MHz in SFD is 1.49 times larger than at f_H + f_L = 108.48 + 13.56 MHz in DFD. Next, the rf power density also increases with increasing driving frequency, similarly to the electron density, in both single- and dual-frequency discharges. In addition, the input power density in SFD is larger than in DFD across all driving frequency ranges. For example, the input power density at f_rf = 27.12 MHz in SFD is 1.51 times larger than at f_H + f_L = 27.12 + 13.56 MHz in DFD. Similarly, the input power density at f_rf = 108 MHz in SFD is 1.95 times larger than at f_H + f_L = 108.48 + 13.56 MHz in DFD. The increase in electron density with increasing driving frequency can be attributed to (1) enhanced electron acceleration (stochastic heating effect) due to the increase in driving frequency³⁴ and (2) the increase in input power density.³⁴  Figure 2(c) shows the dependence of the electron density per unit power density on the driving frequency. This result indicates that DFD achieves a higher electron density per unit input power density. Specifically, at f_rf = f_H = 27.12 MHz, the electron density is 1.36 times higher, and at f_rf = f_H = 108.48 MHz, the electron density is 1.31 times higher in DFD compared to SFD. This suggests that DFD generates plasma more efficiently with higher power efficiency. To investigate the cause of this, electric field fluctuations are analyzed. In this study, the operating gas pressure is low (10 mTorr), with the calculated electron mean free path at this pressure being ∼27 mm. The sheath width ranges from 3 to 7 mm, and the bulk plasma length is between 35 and 44 mm. As a result, electrons undergo at most one collision, making collisional heating in the bulk plasma negligible. Therefore, the efficient plasma generation in DFD is attributed to an alternative heating mechanism present within the discharge.

To investigate the cause of this electron heating mechanism, electric field fluctuations were analyzed. First, the analysis results for SFD at f_rf = 13.56 MHz are presented. Figure 3 shows (a) the spatiotemporal evolution of the electric field, (b) the time evolution of the electric field, and (c) the auto power spectrum of the electric field in the plasma bulk region (z = 25 mm) at SFD with f_rf = 13.56 MHz.

The spatiotemporal evolution of the electric field shows one period of 13.56 MHz (∼73.7 ns) on the vertical axis. The horizontal axis represents the distance between the electrodes. It can be observed that the electric field structure changes during each period of the driving frequency f_rf. The electric field intensity is larger from the plasma sheath edge toward the electrode. The time evolution of the electric field at the plasma bulk indicates that the time-averaged value is almost 0 V/m, but it oscillates around ±300 V/m. To evaluate these oscillation components, FFT analysis was performed (2 × 10⁵ sample points). From this result, it was found that not only the fundamental frequency f_rf = 13.56 MHz but also harmonics such as 2 f_rf, 5 f_rf and others exhibit sharp peaks. In addition, a broad peak around 0.21 GHz was observed. This peak corresponds to the electron plasma frequency (⁠$ωpe=2πfpe=nee2ε0me$⁠, where n_e is the electron density in the plasma bulk). It was found that the auto power spectrum decays sharply beyond this peak. The appearance of a peak at the electron plasma frequency in the electric field fluctuation spectrum has also been observed in previous studies,^34,35 suggesting that this is due to the resonance between electron motion and electric field fluctuations.^46,46

Focusing on the conditions with low driving frequency f_rf = f_H = 27.12 MHz and high driving frequency f_rf = f_H = 108.48 MHz, we explain the heating mechanisms in single- and dual-frequency discharges by discussing the spatiotemporal structure of the electric field and the electric field fluctuations and time-averaged electron heating in the plasma bulk and sheath regions. Figure 4 shows the spatiotemporal evolution of the electric field at (a) 27.12 MHz, (b) 27.12 + 13.56 MHz, (c) 108.48 MHz, and (d) 108.48 +13.56 MHz for one RF (13.56 MHz) period. The corresponding electric field in the plasma bulk and sheath of the discharge, as well as its fast Fourier transform (FFT), are presented in Fig. 5. Here, the sheath position refers to the location near the RF electrode where the electron power absorption is the highest. In both single- and dual-frequency discharges, the electric field fluctuation in the plasma bulk slightly increases as the driving frequency rises.

In plasma bulk, for SFD with f_rf = 27.12 MHz, E_z = 6.2 ± 172 V/m [as shown later in Fig. 5(a1)], and for DFD with f_H + f_L = 27.12 + 13.56 MHz, E_z = 43.5 ± 193 V/m [Fig. 5(c1)]. Similarly, for SFD with 108.48 MHz, E_z = −6.7 ± 273 V/m [Fig. 5(e1)], and for DFD with f_H + f_L = 108.48 + 13.56 MHz, E_z = 3.4 ± 252 V/m [Fig. 5(h1)], showing slightly higher electric field fluctuations in the plasma bulk as the driving frequency increases. In addition, in both single- and dual-frequency discharges, the sheath thickness decreases as the driving frequency increases. The sheath thickness is approximately 7–8 mm at f_rf = f_H = 27.12 MHz and about 2–3 mm at f_rf = f_H = 108.48 MHz. In the spatiotemporal structure of the electric field in single- and dual-frequency discharges, the major difference lies in the electric field structure between the sheath edge and the electrodes. In SFD, the electric field oscillates according to the f_rf frequency, but in DFD, the electric field fluctuation at f_H is modulated by f_L. For example, at t/T = 0.50, at the RF electrode, in SFD (27.12 MHz), the field is a few kV/m, whereas in DFD (27.12 + 13.56 MHz), the field strength exceeds 20 kV/m. Conversely, there are phases where the electric field strength in DFD is weaker than in SFD. This characteristic of DFD is that the electric field fluctuation at f_H is modulated by f_L, resulting in significant changes in the electric field structure from the sheath to the electrodes, which leads to sheath modulation. This can also be explained by the differences in the time-series and power spectra of the electric field fluctuations in the plasma sheath and bulk regions.

Figure 5 shows the time evolution and auto power spectrum of fluctuation components of E_z in plasma bulk (left side) and sheath regions (right side) in condition of SFD (a1), (a2) and (b1), (b2) f_rf = 13.56 MHz, (e1), (e2) and (f1), (f2) f_rf = 108.48 MHz and DFD (c1), (c2) and (d1), (d2) f_H + f_L = 27.12 + 13.56 MHz, (g1), (g2) and (h1), (h2) f_H + f_L = 108.48 + 13.56 MHz. First, in the power spectrum of electric field fluctuations, there is a peak near the electron plasma frequency in both the plasma sheath and bulk regions, and electric field fluctuations above this frequency decay significantly. This suggests that the electric field oscillates around the electron plasma frequency and that these oscillations are related to the motion of the electrons.

In the plasma bulk region, since the electric field fluctuation at the driving frequency is relatively small, the peak around the electron plasma frequency is dominant [as shown on the left side of Fig. 5(b1), (d1), (f1), and (h1)]. Therefore, the time evolutions of the electric field in the plasma bulk region [as shown on the left side of Fig. 5(a1), (c1), (e1), and (g1)] mainly oscillate around the electron plasma frequency. In addition, at lower frequency ranges, there are peaks corresponding to the driving frequency’s electric field fluctuation and harmonic components resulting from their nonlinear interactions. The auto power spectrum of the electric field in the plasma bulk region for SFD with f_rf = 13.56 MHz [as shown in Fig. 5(b1)] shows several peaks at f_rf, 3 f_rf, 5 f_rf, and so on. In contrast, for DFD in the bulk region, there are peaks not only at f_H and its harmonics but also sideband peaks at f_H + f_L. This shows that DFD has more nonlinear waveforms and multiple peaks.

Next, let us discuss the electric field fluctuations in the plasma sheath region. The electric field strength of the driving frequency is very strong, and its nonlinearity is also high, showing a filamentary structure [as shown on the left side of Fig. 5(b2), (d2), (f2), and (h2)]. Electric field transients show a filamentary structure indicative of an enhanced nonlinearity in the discharge. In SFD, harmonics of the fundamental frequency, such as f_rf, 2 f_rf, and 3 f_rf, appear [Fig. 5(b2) and (f2)]. In DFD, the nonlinear coupling waves of f_H (such as f_H, 2 f_H, 3 f_H, …) and f_L (such as mf_H + nf_L, where m = 0, 1, 2, 3,… and n = 1, 2, 3,…) appear [Fig. 5(d1) and (h1)]. Therefore, in the electric field fluctuation spectrum for DFD with f_H + f_L = 108.48 + 13.56 MHz, many sidebands at 108.48 MHz ± n × 13.56 MHz appear. This shows that DFD has more peaks in the electric field fluctuations than SFD. The electric field fluctuation near the electron plasma frequency exhibits a broad peak. This is suggested to be because the electron density in the plasma sheath region has a spatial gradient, and the electron plasma frequency in the sheath region has a width.

To compare the magnitude of the electric field fluctuations in SFD and DFD, the electric field fluctuation components up to the electron plasma frequency, where electrons are accelerated, were integrated from the auto power spectrum of the electric field. The integral value was then normalized by the electron plasma frequency to derive the electric field strength component per unit frequency interval (Normalized |E_z|²). Figure 6 shows the normalized total electric field component of SFD and DFD in (a) plasma bulk and (b) sheath regions. In both the plasma sheath and bulk regions, the normalized electric field fluctuation component in DFD is larger than in SFD at all frequency bands (except in the case of f_rf = f_H = 27.12 MHz in plasma bulk). Specifically, in the plasma sheath region, the normalized electric field component in DFD is significantly larger. At f_rf = f_H = 54.24 MHz, it is 1.42 times larger; at 81.36 MHz, 1.55 times; and at 108.48 MHz, 1.66 times larger. These results suggest that DFD has larger electric field fluctuation components than SFD, and this is due to the increased nonlinear coupling of mf_H + nf_L. To investigate how this larger electric field component affects the heating mechanism, we examine the power absorption density distribution per electron.

Figure 7 shows the radial profiles (z = 25–50 mm, from plasma bulk to RF electrode) of normalized electron power absorption in the case of SFD and DFD with f_rf = f_H = (a) 27.12 MHz, (b) 54.12 MHz, (c) 81.36 MHz, and (d) 108.48 MHz. Figure 7(e) indicates the driving frequency dependence of spatio-integrated normalized electron power absorption, which is the integrated value of normalized electron power absorption from plasma bulk to RF electrode. Normalized electron power absorption refers to the amount of power absorbed by an electron per unit volume, which is an indicator of how much power electrons receive from the electric field. These results show that at all frequency bands, the normalized electron power absorption in the plasma sheath region is larger for DFD than for SFD. In addition, the spatio-integrated normalized electron power absorption is ∼1.7–3 times larger in DFD at all frequency bands. The larger normalized electron power absorption indicates that energy is being efficiently received from the electric field. In other words, these results show that in DFD, the electric field fluctuations are larger, and electrons can receive power more efficiently. Therefore, the power efficiency is higher, and electric field fluctuations are a key factor in accelerating electrons. The results suggest that DFD has higher efficiency.

In the sheath region, it is difficult to precisely derive the EEDF information. Therefore, in this study, we focused on the electric field fluctuations and clarified the plasma generation mechanism in DFD by comparing it with SFD. Although this study was a simulation with constant applied voltage, simulations under constant input power density conditions remain as future work. The control of plasma generation and sheath modulation in the sheath region is crucial for controlling the composition ratio of radical, ion, and neutral species during one RF of the f_L period. These results suggest that plasma generation can be more precisely controlled.

We investigated the electron heating mechanism in dual-frequency discharge (DFD) with higher driving frequencies f_H = 27.12 MHz (=2 f_L) to 108.48 MHz (=8 f_L) and f_L = 13.56 MHz, focusing on electric field fluctuations and comparing it with single-frequency discharge (SFD). In typical dual-frequency discharges, f_L is selected near the ion plasma frequency (e.g., 2 MHz) and is used to control ion behavior. However, in this study, we investigated how setting f_L to a higher value of 13.56 MHz, exceeding the ion plasma frequency, affects the electron heating process.

These results revealed that DFD achieves higher electron density per unit input power density compared to SFD. To investigate the underlying cause, we analyzed the spatiotemporal structure of the electric field and the auto-power spectrum of electric field fluctuations. In DFD, more electric field fluctuation peaks were observed, including not only harmonics of f_H but also nonlinear coupling waves of f_H and f_L. These fluctuations were particularly prominent in the sheath region, and the electric field fluctuation components in the frequency range up to the electron plasma frequency were larger in DFD than in SFD. Furthermore, the normalized electron power absorption densities were also higher in DFD, suggesting that the large electric field fluctuations significantly contribute to electron heating (stochastic heating). A higher normalized electron power absorption indicates more efficient energy transfer from the electric field to electrons. These findings suggest that DFD with a higher frequency f_L achieves greater electric field fluctuations, enabling electrons to efficiently gain energy from the electric field, leading to higher plasma generation efficiency compared to SFD.

Since the electron motion in the sheath region does not reach thermal equilibrium, the EEDF cannot be derived and is typically used only as an evaluation metric for the plasma bulk. This study highlights that electric field fluctuation analysis is a valuable approach for investigating electron heating mechanisms, particularly stochastic heating, in the plasma sheath region, which significantly contributes to plasma generation.

This work was partly supported by JSPS KAKENHI (Grant No. J.P.23K03368 and J.P.24H00205) and research grants from the Amano Institute of Technology and the Kakihara Science and Technology Foundation.

The authors have no conflicts to disclose.

Hiroshi Otomo: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Jian-Syun Lai: Data curation (equal); Investigation (equal); Methodology (equal); Visualization (equal). Kunihiro Kamataki: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Yuma Yamamoto: Investigation (equal); Methodology (equal). Masaharu Shiratani: Conceptualization (equal); Funding acquisition (equal); Writing – review & editing (equal).

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