# Modified Fluid Simulation of an Inductively Coupled Plasma Discharge

Dong-gil Kim, YoungDo Jeong, Yejin Shon, Deuk-Chul Kwon, Jae-hong Jeon, and HeeHwan Choe

## Abstract

Recent developments in semiconductor technology have resulted in advancements in device
manufacturing processes, especially when plasma discharges are required. In this investigation,
an inductively coupled plasma (ICP) discharge with a radio frequency bias simulation
was explored. A 2D symmetric structure is presented and a modified fluid model is
applied using the electron distribution function (EEDF). This is achieved by using
a two-term Boltzmann approximation and an ion energy distribution function (IEDF)
with the particle-in-cell method. In addition, the heat transfer and gas flow effects
are considered to improve the simulation results. The application of multi-physics
components to the ICP discharge displays hybrid-like results using COMSOL[

**Keywords:**Plasma, ICP, EEDF, IEDF, Heat transfer, Flow

## I. Introduction

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.

## II. Model description

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

where *n** _{e}* is the electron density,

*μ*

*is the electron mobility,*

_{e}*D*

*is the electron diffusivity, and*

_{e}*R*

*is the electron rate expression.*

_{e}where *n** _{ɛ}* is the electron energy density,

*E⃗*is the electric field,

*μ*

*is the electron energy mobility,*

_{ɛ}*D*

*is the electron energy diffusivity, and*

_{ɛ}*S*

*is the energy loss or gain as featured in Fig. 1. The reaction scheme used is presented in*

_{en}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*

*is the mixture averaged diffusion coefficient,*

^{m}*M*

*is the mean molar mass of the mixture,*

_{n}*D*

*is the thermal diffusion coefficient for the species, and z is the charge number for the species.*

^{T}To find the electric potential, Poisson’s equation was solved:

For Poisson’s equation, *e* is the electric charge and *ɛ*_{0} 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_{1}, A_{2}, and A_{3} 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*

*is the background gas density.*

_{n}The particle method was applied to obtain the kinetic characteristics and the IEDF of the ions.

where mi is the ion mass and
*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,

*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*

*is the reaction rate.*

_{j}The gas temperature in the inlet was assumed to be 300 K, −*↛* • *q⃗* = 0 was used in the gas outlet, and
*↛* 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.

## III. Results and discussion

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.

## IV. Conclusions

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.

## Acknowledgements

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

## Article information

###### Articles from Applied Science and Convergence Technology are provided here courtesy of **Applied Science and Convergence Technology**

## References

- COMSOL 5.3a.Array. accessed Dec. 3, 2019
- Choe HH and Yoon NS. J Korean Phy Soc 2003;859:42.
- Makabe T and Petrović ZL. Plasma Electronics: Applications to Microelectronic Device Fabrication. New York: Taylor and Francis; 2006