Inductively coupled plasma (ICP) discharge is widely used in semiconductor manufacturing processes [1–2]. Due to the advantage of increasing plasma uniformity and reducing charge-induced damage, pulsed ICP is used instead of continuous wave (CW) plasma [3–5]. As pulse power with a specific pulse frequency and duty cycle is applied to the antenna, the high ion energy is generated in the pulse-on state (high-state) and used for etching. In the off state (zero-state), the density of the electron and cations decrease as the plasma sheath disappears. In Ref. [6], Han et al. measured key plasma parameters, such as electron density, electron temperature, and electron energy distribution function (EEDF), in pulsed argon ICP used in industrial etching processes through different probes.

EEDF is required for the calculation of reaction coefficients for electron collision reactions, while ionization and excitation processes are greatly influenced by EEDF. Therefore, the plasma calculation results vary, based on which EEDF is considered. Maxwellian assumes that electrons are in thermal equilibrium and used under a high ionization degree. However, because this study has a low degree of ionization, Druyvesteyn may produce more appropriate results [7].

In plasma simulation based on a fluid model, the drift-diffusion approximation is commonly used in the momentum equations to reduce the computational resources. Some terms in the momentum equations are neglected by assuming that the electron motion is dominated by the collision and in a steady-state [8]. This approximation is suitable for use under certain conditions: collisional plasma, gas pressure over 10 mTorr, and very low ionization degree [7]. Because the process conditions in this study are suitable for using the approximate equation, the drift-diffusion approximation is applied to the momentum equation of electron density and electron energy density.

This study was carried out based on the information presented in the experimental paper [6], with the same process conditions used in the simulation model, assuming a two-dimensional (2D) axisymmetric structure. In the simulation, cases were appropriately partitioned for each probe used in the measurement, and all cases were performed with Druyvesteyn EEDF because the EEDF of ICP is often non- Maxwellian. To obtain results close to the experiment, heat transfer and gas flow effects are considered, as shown in Fig. 1.

Figure 1. Algorithm for computing ICP.

This paper is organized into 4 sections: 1) Introduction for the background and brief overview; 2) Model description, including the overall scheme, and explanations for the equations; 3) Results and discussion of comparing the simulation results with experimental measurements to validate the study; 4) Conclusion for the result analysis and discussion of future work.

2.1. Model description

For comparison, the model is constructed almost similar to the equipment used in the experimental paper, whose structural diagram is shown in Fig. 2. Although the experimental equipment is three-dimensional, a 2D axisymmetric model is assumed to reduce com- putational resources.

Figure 2. Structural diagram for the simulation: (a) the experiment equipment and (b) rotated 2D axisymmetric model used in simulation.

The pre-rotation structure in Fig. 2(b) is demonstrated in Fig. 3. There are six plasma chamber areas: antenna, air cap, ceramic window, wafer, chuck, and plasma region. The antenna is made of copper with a three-turn coil with a 3 mm radius and is surrounded by an air cap that divides the area from the outside. Pulsed power was applied with the following conditions: duty cycle 50 %, pulse frequency 1 kHz, driving frequency 2 MHz, and current 24 A for the peak power of 300 W. Below it is a ceramic window composed of Al₂O₃, which separates the antenna from the plasma region and allows the transmission of electromagnetic waves to the plasma region [9]. In the plasma region, Ar gas was injected into the inlet under the conditions of 100 sccm and 25 mTorr, with an outlet at the chamber bottom. Although it does not participate in the calculation, a chuck and a 300-mm wafer were constructed under the plasma region to aid visual understanding. Additionally, the collision processes in the Ar discharge are summarized in Table I.

Table 1 . Important collision processes in the argon discharge.

Reaction	Formula	Type	∆ε (eV)	Ref.
1	e+Ar→e+Ar	Elastic	0	[10]
2	e+Ar→e+Ars	Excitation	11.5	[10]
3	e+Ars→e+Ar	Superelastic	−11.5	[10]
4	e+Ar→2e+Ar+	Ionization	15.8	[10]
5	e+Ars→2e+Ar+	Ionization	4.24	[10]
6	Ars+Ars→e+Ar+Ar+	Penning ionization	-	[10]
7	Ars+Ar→Ar+Ar	Metastable quenching	-	[10]

Figure 3. Schematic of the ICP model.

This research is based on a fluid model and carried out with COMSOL Multiphysics 5.5.

2.2. Equations

When the AC current is driven to the antenna, a time-varying magnetic field is created in the rz direction. From this, the induced electric field is created in a vertical (θ) direction, transferring the energy that can cause collisions inside the chamber. Eqs. (1) to (3) enumerate the calculation formulas for these electromagnetic fields. Both the capacitive and inductive fields are considered when calculating the electric field, as indicated by Eq. (3).

where $\overrightarrow{H}$ is the magnetic field, $\overrightarrow{J}$ is the current density, $\overrightarrow{B}$ is the magnetic flux density, while $\overrightarrow{A}$ and V indicate the magnetic vector and electric potentials, respectively.

The electrons that obtain energy from the induced electric field generate plasma through the collisions, as shown in Table I. The electron density is defined by the drift-diffusion approximation, as shown in Eq. (5).

The electron energy density is defined by Eq. (6) and follows the drift-diffusion approximation.

where n_e is the electron density, ${\overrightarrow{\Gamma }}_e$ is the electron flux, ${\left(S-L\right)}_e$ is the difference between the electron source and electron loss, ${\mu }_e$ is the electron mobility, D_e is the electron diffusivity, k is the Boltzmann constant, n_ε is the electron energy density, ${\overrightarrow{\Gamma }}_{\epsilon }$ is the electron energy flux, S_en is the energy loss and gain due to the collisions, ${\mu }_{\epsilon }$ is the electron energy mobility, and $D_{\epsilon }$ is the electron energy diffusivity.

The generalized electron energy distribution function is shown in Eq. (9), and the definition of the variable used in the equation is shown in Eqs. (10) and (11): Depending on whether the number of ‘g’ is 1 or 2, it turns out to be Maxwellian EEDF and Druyvesteyn EEDF, respectively.

where $\varphi$ is the mean electron energy, ε is the electron energy, and $\Gamma$ is the upper incomplete gamma function.

Heavy species transport properties are calculated according to the mixture-averaged diffusion model, as given in Eqs. (12)–(14).

where ρ is the density of the mixture, w is the mass fraction of the species, $\overrightarrow{u}$ is the mass-averaged fluid velocity vector, $\overrightarrow{j}$ is the diffusive flux vector, R is the rate expression for species, $\overrightarrow{V}$ is the multicomponent diffusion velocity for species, D_m is the mixture-averaged diffusion coefficient, M_n is the mean molar mass of the mixture, T is the gas temperature, $D_m^T$ is the thermal diffusion coefficient, z is the charge number for the ${\mu }_m$ ixture averaged mobility for species.

The velocity and pressure fields are calculated based on the Navier–Stokes equations, where Eqs. (15) and (16) represent the momentum and mass conservation, respectively. Eq. (16) shows that the motion of the fluid follows the continuity equation.

where ρ is the gas density, $\overrightarrow{u}$ is the gas flow velocity, p is the pressure, $\overrightarrow{I}$ is the identity matrix, and $\mu$ is the dynamic viscosity. The above equations are for the domain. For the boundary condition of a plasma wall, the fluid velocity is zero as the wall is stationary.

Considering heat transfer effects means that the convection of energy, thermal conduction, radiation, and a heat source are calculated as shown in Eq. (17) to control the internal energy over time and finally conserve energy.

where C_P is the heat capacity, $\overrightarrow{q}$ is the heat source for gas, and k is the thermal conductivity. The temperature at the inlet and wall is constant at 300 K, while the gradients of temperature at the outlet in the normal direction are zero.

3.1. Comparison of EEDF and the physics effect in CW plasma

In this section, simulations were performed under continuous wave conditions to determine the effects of EEDF and additional physics. Figure 4 illustrates the distribution of normalized plasma parameters, representing the maximum values for each by the numbers.

Figure 4. Plasma parameters classified according to EEDF and physical effects in CW plasma conditions.

First, based on the difference in EEDF, it can be observed that the results considering Druyvesteyn EEDF have relatively lower electron density values, higher electron temperatures, and electric potential values than those of the Maxwellian EEDF. The Maxwellian EEDF exhibits a high maximum electron density by overestimating the ionization rate, and Druyvesteyn assumes that elastic collisions occur more frequently than inelastic collisions, which explains the above differences. When comparing the simulation results, the degree of ionization of Maxwellian EEDF was approximately 4.7 times higher than that of Druyvesteyn EEDF. Looking at the reaction rate of Druyvesteryn EEDF, elastic collisions showed 496 times higher than the excitation reaction and 781 times higher than the ionization reaction.

Next, the changes due to the effects of additional physics are compared. In Fig. 5, there are fluctuations in pressure and temperature, which are fixed as constant when not considering the effects of gas flow and heat transfer. As temperature increases, the electron density decreases according to the ideal gas law. Because the velocity and pressure fields are calculated by applying the gas flow effect, the pressure has the highest value at the inflow part, and the movement of fluid from inlet to outlet is considered as shown in a streamline.

Figure 5. Influences of additional physics effects.

3.2. Pulsed ICP with additional physical effects

This section compares the simulation results and plasma parameters measured using the double-sided planar Langmuir probe from Ref. [6]. The effects of gas flow and heat transfer were merged to create conditions similar to the actual process. The analysis of the results was performed in one pulse period, and the plasma parameters are observed at the center of the chamber(z = 81.5 mm), as shown in Fig. 6.

Figure 6. Plasma parameters observed at the chamber center (z = 81.5 mm).

Overall, for all three parameters, the trends of increasing or decreasing with the pulse-on and - off states are similar, but there exist slight numerical differences. This is because the simulation was conducted with Druyvesteyn EEDF, rather than the Boltzmann EEDF. The Boltzmann EEDF can calculate the drift velocity, a significant parameter in plasma reactions, and achieve maximum consistency in the kinetic and fluid description of electron dynamics [7]. Despite these advantages owing to the fatal drawback of the Boltzmann EEDF being computationally expensive, Druyvesteyn EEDF was used in this simulation to save computational resources.

Fundamental plasma parameters according to pulse-on and -off states were compared under the same conditions. The characteristics of the pulsed plasma in the off state are clearly illustrated in Fig. 7. The plasma sheath, which has a relatively large number of positive ions and neutrons than electrons, disappears as the power turned off. Therefore, a very low numerical value can be observed as the electrons are cooled. In addition, the density of electrons and cations decrease, which is a great feature of the pulse-off state. To more effectively observe the density decrease, a comparison according to the r direction at the center of the plasma region is shown in Fig. 8.

Figure 7. Comparison of pulse on-off state.

Figure 8. Comparing electron density in pulse on-off state.

3.3. Pulsed ICP without additional physical effects

Section 3.3 is prepared for comparison with the Hairpin resonator probe results in Ref [6], which is affected by the presence of a plasma sheath near the tip [6]. As shown in Fig. 4, not considering the heat transfer and gas flow effects results in a wider sheath. Since the actual measurements were affected by the sheath, additional effects were not considered. The time-varying electron density in one pulse cycle is presented in Fig. 9. Figure 9(a) is a graph obtained by comparing the normalized simulation result with the experimental value, while the solid and dotted lines represent the simulated and measured values with the hairpin probe, respectively. As mentioned in section 3.2, the numerical value is not exactly suitable because the Boltzmann EEDF is not considered, but the trend along the r direction is similar. Figure 9(b) is a graph showing only the simulation results. One of the noticeable results in the experiment, it can be seen through simulation that the lowest profile does not appear at the start of discharge.

Figure 9. Time-varying electron density in pulse duration. (a) Normalized electron density comparison of measured and simulated values and (b) computed electron density in the simulation.

Additionally, conditions for obtaining valid plasma simulation results were studied. In the plasma region, the electron density along the z-direction at r = 0 mm was updated for every pulse cycle. The 1^st pulse period infers invalid plasma parameter values, as shown in Fig. 10. To obtain accurate simulation results, the values after the 2^nd cycle have to be observed.

Figure 10. Electron density according to the pulse cycle.

This study presented a simulation based on an experimental paper that measures the plasma parameters of pulsed ICP in a 300-mm wafer size chamber used in a semiconductor etching process. To reduce computational resources, a 2D axisymmetric structure was assumed and performed based on a fluid model. Because the electron energy distribution functions in ICPs are frequently Non-Maxwellian, Druyvesteyn EEDF was considered. The calculation accuracy is improved in calculating the electromagnetic field by considering not only the induced electric field but also the capacitive one. The gas flow effects account for changes in pressure and velocity fields, and the heat transfer physics for considering changes in enthalpy due to collision were merged. Appropriate physics settings were applied to each of the Langmuir and hairpin probes used in the experiment, and the simulation was validated through a comparison of the results. As the Boltzmann EEDF was not considered, which can obtain relatively more accurate electron transport properties than the Druyvesteyn EEDF, there were some numerical differences, but the profiles of the plasma parameters were similar. In addition, the characteristics of the pulsed plasma in the off state were examined. As the power is turned off, the sheath disappears, electrons cool down, and densities of electrons and cations decrease. In pulsed ICP simulation, it is necessary to observe after the 2^nd cycle to obtain accurate results.

To obtain more suitable plasma properties, the Boltzmann EEDF calculated using the two-term Boltzmann equation should be considered. Moreover, by applying the RF bias for capacitively coupled ICP, uniformity can be improved by controlling the ion collision energy at the substrate independently of the plasma density. Moreover, the etch rate and profile could be predicted through the ion energy and angular distribution functions according to the RF bias.

This work was supported by the National Research Council of Science & Technology (NST) grant by the Korean government (MSIP) (No. CAP-20-02-NFRI). And we would like to express our gratitude and appreciation to the Altsoft for supporting this work.