Due to recent advancements in semiconductor technology, the design rules for developing semiconductors is becoming more stringent [2,3]. Device manufacturing processes have become more important, especially when plasma discharges are used [4,5]. An inductively coupled plasma (ICP) has a high-density plasma source with a simple geometry [6–10]. In this work, a simulation method using two-dimensional (2D) symmetric geometry was studied. Although a simulation reflecting three-dimensional (3D) geometry can be used, 2D simulation methods have been widely adapted as one of the most effective methods, especially when they are combined with other physics.

Details of the various plasma simulation models can be found in the literature [11–16]. In this work, the simulation was based on a fluid model. The particle-in-cell (PIC) method, or the kinetic model, have been studied for plasma simulations; there are many results when combined with the fluid model [17–19]. These hybrid models are believed to be one of the most efficient methods. For low-pressure discharge modeling, it is well known that the electron energy distribution function (EEDF) and the ion inertial terms are important for discharge simulations. However, it is hard to obtain a hybrid simulation model for commercially available modeling packages when considering the ion inertial terms since it is not easy to get stable and convergent results.

The EEDF and the IEDF were considered in this simulation using COMSOL to improve the fluid model. However, the methods studied in this work are different from standard hybrid methods. Detailed explanations of the model descriptions are provided in Section 2. The results and discussion are presented in Section 3. Finally, the conclusion is described in Section 4.

The simulation model used in this study has a 2D axisymmetric geometry. The concept of the model used is shown in Fig. 1. The width of the chamber is 8.26 cm, the height is 8 cm, and the thickness of the window is 1 cm. The antenna features a circular ring shape and a square of 0.3 cm in the cross-section. There are two gas ports with a length of 1.26 cm. The length of the dielectric material around the gas port is 1.27 cm. The frequency and power applied to the coil were 13.56 MHz and 100 W, respectively, and the bias frequency and power were 3.4 MHz and 20 W, respectively. The blocking capacitor located in the bias power was set to 0.1 μF. The material of the window is aluminum oxide and the material of the antenna is copper. The pressure of Ar gas was assumed to be 30 mTorr.

Equations for the electron density and electron energy density are given in Equations (1–4).

where n_e is the electron density, Γe→ is the electron flux, μ_e is the electron mobility, D_e is the electron diffusivity, and R_e is the electron rate expression.

where n_ɛ is the electron energy density, Γe→ is the electron energy flux, E → is the electric field, μ_ɛ is the electron energy mobility, D_ɛ is the electron energy diffusivity, and S_en is the energy loss or gain as featured in Fig. 1. The reaction scheme used is presented in Table I.

The following equations are used to describe all of the heavy species, except the electron:

where ρ is the density of the mixture, w is the mass fraction of the species, R is the rate expression for the species, J → is the diffusivity flux vector, μ_m is the mixture averaged mobility for the species, V → is the multicomponent diffusion velocity for the species, D^m is the mixture averaged diffusion coefficient, M_n is the mean molar mass of the mixture, D^T is the thermal diffusion coefficient for the species, and z is the charge number for the species.

To find the electric potential, Poisson’s equation was solved:

For Poisson’s equation, e is the electric charge and ɛ₀ is the permittivity of the vacuum. In general, the EEDF is obtained through small cell units. In this work, the cell for calculating the distribution functions are larger than the conventional method. The EEDF and IEDF were used for the correction of the transport coefficients and the reactions. It is expected that if the size of the cell is sufficiently small, the result would be similar to the usual hybrid method. In particular, the ion mobility values extracted from the IEDF become reduced when they are in the vicinity of the sheath. This is supposed to complement the inertial terms that are ignored when describing the ion fluxes [21]. The area divided in this simulation is shown in Fig. 2. The IEDF of each region was calculated by using the particle model to change the ion mobility and diffusion coefficient. In addition, the ion temperature was estimated from the IEDF as demonstrated in Fig. 2. The algorithm for applying the EEDF and IEDF are presented in Fig. 3.

The following equation was applied to obtain the EEDF for each cell:

The variables W, D, A, γ, A₁, A₂, and A₃ are as follows:

where S is the energy loss due to the inelastic collisions, m_e is the electron mass, ɛ is the energy, σ_ɛ is the total elastic collision cross-section, q is the electron charge, T is the temperature of the background gas, k_b is the Boltzmann constant, and N_n is the background gas density.

The particle method was applied to obtain the kinetic characteristics and the IEDF of the ions.

where mi is the ion mass and vi→ is the ion velocity. The probability of the collision P is calculated as follows:

where ν is the collision frequency and Δt is the time interval. The probability distribution function is assumed to be a drifting Maxwell distribution as follows:

where m_g is the mass, vg→ is the velocity, and u → is the mean or drift velocity of the background Ar.

In this simulation, heat transfer was considered to find the effect of the thermal properties for the plasma discharge. The heat sources used to obtain the thermodynamic properties were calculated from the enthalpy changes of the reactions. The following equation was used to obtain the thermodynamic properties:

where ρ is the gas density, C_p is the heat capacity, u → is the gas flow velocity, q → is the conductive heat flux; and Q is the heat source of the gas.

where H_j is the enthalpy of the reaction and r_j is the reaction rate.

The gas temperature in the inlet was assumed to be 300 K, - n → · q → = 0 was used in the gas outlet, and -n→·q→=-k1615πλTwallT (T-Twall) was used in the wall. Here, n → is normal unit vector, k is thermal conductivity, λ is mean free path, and T_wall is wall temperature.

For the gas flow properties, the following equation was applied:

where ρ is the gas density, u → is the gas flow velocity, p is the pressure, and μ is the dynamic viscosity. A boundary condition was applied where the gas flow rate was Q_in = 100 sccm, the gas outlet pressure was adjusted to 30 mTorr, and u → = 0 at the wall.

Figure 4 illustrates the consideration of the heat transfer effect. The gas temperature is highest near the center of the plasma region and decreases along the edge of the chamber. A neutral density appears along the lowest distribution near the center of the plasma region. This increases towards the edge of the chamber and the maximum value appears at the gas inlet. The temperature and density profiles around the center of the chamber are noticeably different. This is because the enthalpy change is most active at this point. The neutral density was almost the same, except for the gas inlet and outlet. The heat transfer effect changed the maximum electron density as demonstrated in Fig. 4.

Figure 5 shows the Maxwellian EEDF and Boltzmann EEDF values for cell numbers 1–7. As the electron energy increases, the difference between the Maxwellian EEDF and Boltzmann EEDF increases, which is noticeable at more than 10 eV. For the Boltzmann EEDF, the EEDF value decreases toward the upper (cell 7) or lower (cell 1) direction with respect to the z axis of the cell. These results indicate that the EEDF should be considered as a function of the positions to obtain more accurate results.

Figure 6 displays the relative count according to the particle kinetic energy at the bottom electrode. The result using the empirical ion temperature model is different from that of the IEDF by the PIC method. The empirical ion temperature model with the Maxwellian EEDF has a maximum value near 10 eV and a broad peak near 20–30 eV. The model using 28-cell EEDF and IEDF featured more clear peaks around 10 eV and 47 eV.

Figure 7 presents the ion temperature and mobility. They have different values for every 28 cells. In particular, the ion mobility shows the highest value in the middle region of the chamber, and the value decreases toward the outside edge of the chamber except for the gas outlet. The decreased ion mobility near the sheath region is supposed to have a similar effect when considering the ion inertial term. This is because it reduced the drift part of the ion flux [21]. The ion temperature shows the usual profiles.

Figure 8 shows the plasma characteristics for the Maxwellian EEDF and the empirical ion energy in comparison to the EEDF and the IEDF in this article. Profiles and values of the electron density were slightly different due to the differences in the transport and reaction coefficients because of the EEDF. Also, the electron temperature profiles and the uniformity were changed. In terms of the electric potential, relative differences appeared near the wall due to the change in the ion transport coefficients. The ion inertial terms are more important near the sheath region since the plasma characteristics change near the sheath. Although the sizes of the cells used in this model are larger than the traditional hybrid model, this model reflected the meaning of the distribution functions.

Simulation for the plasma characteristics in the ICP discharge with a radio frequency bias was investigated using COMSOL. In this model, the two-term Boltzmann approximation EEDF and the IEDF, with the PIC method, were used for very large cell sizes. Although rough cells were used, the results confirmed that this method reflected the EEDF and the ion inertial terms. Moreover, the heat transfer and gas flow were considered to improve the background Ar gas characteristics. This method shows a possibility of improving the commercial simulation tools. This is based on the fluid model with a numerical stability and a relatively small amount of computational resources.

In the near future, a comparable model to the traditional one will be developed, and the results will be published soon.

This work was supported by the National Research Council of Science & Technology (NST) grant by the Korea government (MSIP) (No. CAP-19-02-NFRI).