Applied Science and Convergence Technology

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

Additional article information


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[ 1].

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 [610]. 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 [1116]. 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 [1719]. 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.

Figure F1
(Color online) Schematic diagram of the simulated ICP.

Equations for the electron density and electron energy density are given in Equations (14).

n e t + · Γ e = R e Γ e = - n e ( μ e · E ) - D e n e

where ne is the electron density, Γ e is the electron flux, μe is the electron mobility, De is the electron diffusivity, and Re is the electron rate expression.

n ɛ t + · Γ ɛ + E · Γ ɛ = S e n Γ ɛ = - n ɛ ( μ ɛ · E ) - D ɛ n ɛ

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 Sen 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:

ρ w t + ρ ( u · w ) = · J + R J = ρ w V V = D m In ( w ) + D m In ( M n ) + D T ρ w In ( T ) - z μ m E

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, Dm is the mixture averaged diffusion coefficient, Mn is the mean molar mass of the mixture, DT 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:

2 V = - e ɛ 0 ( n i - n e )

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.

Figure F2
(Color online) Cells for EEDF and IEDF applications.
Figure F3
Algorithm for the adapting EEDF and IEDF.

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

ɛ ( W f - D f ɛ ) = S

The variables W, D, A, γ, A1, A2, and A3 are as follows:

W = - γ ɛ 2 σ ɛ - 3 a ( n e N n ) A 1 D = γ 3 ( E N n ) 2 ( ɛ σ m ) + γ k b T q ɛ 2 σ ɛ + 2 a ( n e N n ) ( A 2 + ɛ 3 / 2 A 3 ) a = q 2 γ 24 π ɛ 0 2 ln ( 12 π ( 2 ɛ 0 q ɛ ¯ / 3 ) 3 / 2 q 3 n e 1 / 2 ) ɛ 2 σ ɛ + 2 a ( n e N n ) ( A 2 + ɛ 3 / 2 A 3 ) γ = ( 2 q m e ) 1 / 2 [ C 1 / 2 / k g 1 / 2 ] A 1 = 0 ɛ u 1 / 2 f ( u ) d u , A 2 = 0 ɛ u 3 / 2 f ( u ) d u , A 3 ɛ f ( u ) d u

where S is the energy loss due to the inelastic collisions, me 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, kb is the Boltzmann constant, and Nn is the background gas density.

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

( m i v l ) t = - e V

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

p = 1 - exp ( - ν Δ t ) , ν = j N ν j

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:

f ( v g ) = 1 ( 2 π k b T / m g ) 3 / 2 exp [ m g | v g - u | 2 2 π k b T ]

where mg is the mass, v g 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:

ρ C p u · T + · q = Q

where ρ is the gas density, Cp 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.

Q = j Q j = j - ( H j r j )

where Hj is the enthalpy of the reaction and rj is the reaction rate.

The gas temperature in the inlet was assumed to be 300 K, −q⃗ = 0 was used in the gas outlet, and - n · q = - k 16 15 π λ T w a l l T ( T - T w a l l ) was used in the wall. Here, is normal unit vector, k is thermal conductivity, λ is mean free path, and Twall is wall temperature.

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

ρ ( u · ) u = · [ - p I + u ( u + ( u ) T ) - ( 2 / 3 ) u ( · u ) I ] · ( ρ u ) = 0

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 Qin = 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 F4
(Color online) Effects of heat transfer.

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 F5
(Color online) Axial scans of the EEDF in comparison to the Maxwellian EEDF.

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 F6
(Color online) IEDFs at the electrode.

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 F7
(Color online) The (a) ion mobility m2/(V · s) and (b) ion temperature[eV] using the PIC method.

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.

Figure F8
(Color online) Plasma characteristics obtained by (a) Maxwellian EEDF and the empirical ion temperature, and (b) 28-cell EEDF and IEDF.

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.


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

Applied Science and Convergence Technology.Nov 30, 2019; 28(6): 221-225.
Published online 2019-11-30. doi:  10.5757/ASCT.2019.28.6.221
aSchool of Electronics and Information Engineering, Korea Aerospace University, Gyeonggu-do 10540, Republic of Korea
bPlasma Technology Research Center, National Fusion Research Institute, Jeollabuk-do 54004, Republic of Korea
Received November 18, 2019; Accepted December 3, 2019.
Articles from Applied Science and Convergence Technology are provided here courtesy of Applied Science and Convergence Technology


  • 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

Figure 1

(Color online) Schematic diagram of the simulated ICP.

Figure 2

(Color online) Cells for EEDF and IEDF applications.

Figure 3

Algorithm for the adapting EEDF and IEDF.

Figure 4

(Color online) Effects of heat transfer.

Figure 5

(Color online) Axial scans of the EEDF in comparison to the Maxwellian EEDF.

Figure 6

(Color online) IEDFs at the electrode.

Figure 7

(Color online) The (a) ion mobility m2/(V · s) and (b) ion temperature[eV] using the PIC method.

Figure 8

(Color online) Plasma characteristics obtained by (a) Maxwellian EEDF and the empirical ion temperature, and (b) 28-cell EEDF and IEDF.

Table I

Ar reactions used in this simulation.

No. Reaction Reference
1 e + Ar → e + Ar [ 20]
2 e + Ar → e + Ars [ 20]
3 e + Ars → e + Ar [ 20]
4 e + Ar → 2e + Ar+ [ 20]
5 e + Ars → 2e + Ar+ [ 20]
6 Ars + Ars → e + Ar + Ar+ [ 20]
7 Ars + Ar → Ar + Ar [ 20]