As the line width of semiconductor devices required by the semiconductor industry becomes finer, detailed processing techniques and equipment are required. This increases the importance of plasma processes with anisotropic etching characteristics [1]. Plasma is widely used in semiconductor processes such as plasma etching, sputtering, and plasma-enhanced chemical vapor deposition [2]. In semiconductor processes, a low process pressure, high ionization rate, high etching rate, and good etching uniformity are required; thus, equipment optimized for each process has been developed [3].

The semiconductor industry mainly uses capacitively coupled plasma (CCP) and inductively coupled plasma (ICP) sources. Although CCP is widely used, the usual drift-diffusion approximation cannot be performed efficiently at process pressures ≤ 20 mTorr; hence, development of a high-density plasma generating device is necessary. Thus, ICP devices, which can generate relatively high-density plasma at lower pressures (0.5–50 mTorr) than those in CCP, are widely used. In ICP, an antenna is installed at the side of a cylindrical device or on the top of a quartz window. In addition, a radio frequency (RF) power is applied to produce an electric field and maintain discharge. Moreover, an independent control of ion density and energy is achieved by applying an RF voltage to the electrode on which the wafer is placed. Several studies have shown the advantages of ICP equipment [4–6].

It is difficult to determine the optimal process conditions as most of the process parameters employed for the structure of the equipment used in semiconductor manufacturing are determined through experiments. Changes in plasma parameters, such as electron temperature, plasma density, and electron energy distribution, have a significant influence on the process results. Thus, for an accurate understanding of the plasma discharge phenomenon, an understanding of the physics of plasma must be preceded through simulation. This study examined physical phenomena that are difficult to understand experimentally through mathematical modeling of plasma discharge phenomena and found correlations between various plasma elements. It is expected that this study will help in determining the optimal process conditions when dealing with plasma equipment.

This study used the most efficient and universal fluid model among the most commonly used plasma modeling methods to calculate the quantity closest to the measurable physical quantity. The model used in this study is described in detail in Section 2. To justify the simulation model used in this study, Section 3.1 compares the model results and experimental data. The Boltzmann electron energy distribution function (EEDF), ion temperature, and ion mobility calculations were applied to the model to obtain stable convergence results in low-pressure discharge. A detailed description of these results is provided in Section 3.2. The plasma properties as a function of process pressure conditions (1–25 mTorr) are discussed in Section 3.3.

In this study, ICP discharge in semiconductor processing was simulated using the COMSOL Multiphysics software. The structure of an ICP model used in the production of a 300-mm wafer was described in [7]. Figure 1 shows the 2D axisymmetric structural model used in this study. An RF power of 13.56 MHz was applied through the antenna installed on the top of the cylindrical chamber (height = 17.2 cm and radius = 26 cm). The chamber and antenna coil were separated by an electric insulator. The coil had a central radius of 16 cm and a spacing of 1.1 cm and was located 20.1 cm above the wafer. The wafer (radius = 15 cm) was fixed with a chuck with a radius of 22.75 cm. Argon gas was used as the process gas and its characteristics were observed while changing the initial vacuum pressure from 1–25 mTorr. The reaction equations used in the simulation are listed in Table I. The inlet and exit ports of the gas were located at the top and bottom of the chamber. Except for the gas inlet and the dielectric material, the boundaries of the chamber were all grounded.

Table 1 . Set of gas reactions in Ar plasma..

Index	Reaction		Reference
1	e + Ar → Ar + e	-	[8]
2	e + Ar → Ars + e	11.55	[8]
3	e + Ar → Ars + e	13	[8]
4	e + Ar → Ars + e	14	[8]
5	e + Ar → Ars + e + e	15.7	[8]
6	Ars + Ars → Ar + Ars + e	-	[9]
7	Ar + Ars → Ar + Ar	-	[9]

Figure 1. Schematic of the cylindrically 2D symmetric inductively coupled plasma reactor geometry.

2.1. Plasma fluid model

The plasma fluid model was used to calculate the electromagnetic field and absorbed power by solving the transmission equation obtained from the Maxwell-Boltzmann equation. This model consists of a continuous equation, a motion equation, and a Poisson equation for each plasma species [10].

The Ampere's law is applied in the plasma region and the antenna, as described in the following equation:

where μ₀ is the permeability of free space, μ_r is the relative permeability of the material, and $\overrightarrow{A}$ is the magnetic vector potential. In the plasma region, the right hand side (RHS) is $\overrightarrow{J}={\overrightarrow{J}}_{ind}+{\overrightarrow{J}}_d$ and can be expressed as the sum of induced and displacement current densities. In the antenna area, the RHS is $\overrightarrow{J}={\overrightarrow{J}}_{ind}+{\overrightarrow{J}}_d+{\overrightarrow{J}}_A,$ and the external current density is added to the current density terms in the plasma region.

The electric field in the plasma was calculated using the Poisson equation:

where e is the electric charge, ε₀ is the permittivity of vacuum, n_i is the ion density, and n_e is the electron density.

The equation for the electron transport is given by:

where S_e is the net electron source term considering reactions related to the generation or loss of electrons.

where x_j is the mole fraction of the target species for reaction j, N_n is the total neutral number density, and k_j is the reaction coefficient calculated from the collision cross-section for the electron-impact reaction. In Eq. (1), the ${\overrightarrow{\Gamma }}_e$ the drift-diffusion equation is given by

where µ_e is the electron mobility and D_e is the electron diffusivity. The boundary conditions for the electron flux at the walls are

where $u_{e,th}=\sqrt{\frac{8k_B{T}_e}{m_e\pi }}$ and is defined as the thermal speed of the electrons.

The equation for the electron energy, similar to the electron continuity equation, is

where S_ε is a variable that includes energy loss/gain and absorption power affected by inductive discharge. As $n_{\epsilon }=\frac{3}{2}{n}_e{k}_B{T}_e$ and $T_e=\frac{2}{3}\overline{\epsilon },$ therefore, $n_{\epsilon }=\overline{\epsilon }{n}_e$ . It follows that

where ${\mu }_{\epsilon }$ is the electron energy mobility and $D_{\epsilon }$ is the electron energy diffusivity. The boundary conditions for the electron flux at the walls are

The diffusion model of neutrons and ions in the plasma domain uses a mixture-averaged model.

where ρ is the density of the mixture, ω is the mass fraction of the species, R is the rate expression for species, and $\overrightarrow{J}$ is the diffusivity flux vector.

where $\overrightarrow{V}$ is the multi-component diffusion velocity for species.

where D_m is the mixture-averaged coefficient, M_n is the mean molar mass of the mixture, and z is the charge number for species.

By specifying the diffusion coefficient and mobility of ions in the low-pressure plasma modeling, stable results can be obtained under the process condition of 1 mTorr. Ion mobility is defined as a function of the reduced electric field for an argon ion. The ion temperature was calculated using the local field approximation, as described in Eq. (13). The diffusivity was calculated using Einstein's relation given in Eq. (14).

where T_g is the background gas temperature, M_ion is the ion mass, M_n is the mean mass of the mixtur , k is the Boltzmann constant, N_A is the Avogadro’s constant, and ${\mu }_i$ is the ion mobility.

The EEDF was used to calculate the rate of the electron collision reaction. Hence, the choice of EEDF, which was used to calculate the plasma parameters, affected the process result. Figure 2 shows the comparison of the Maxwell, Druyvesteyn, and Boltzmann EEDF and this study confirmed the effect of these parameters on the plasma results. Using the Druyvesteyn distribution function, the inelastic collision of electrons and heavy species led to a decrease in the EEDF of high electron energy. When the Maxwell distribution was used, the number of particles with high energy slowly decreased, unlike the other functions. The Druyvesteyn distribution exhibited a maximum value and a distribution similar to the Boltzmann distribution. The Maxwell and Druyvesteyn istribution functions are described in Eqs. (15) and (16), respectively.

Figure 2. Electron probability distribution function of Maxwell, Druyvesteyn, and Boltzmann.

where ${\beta }_1=\text{Γ}{\left(\frac{5}{2}\right)}^{\frac{3}{2}}\text{Γ}{\left(\frac{3}{2}\right)}^{-\frac{5}{2}}, {\beta }_2=\text{Γ}\left(\frac{5}{2}\right)\text{Γ}{\left(\frac{3}{2}\right)}^{-1}.$

where ${\beta }_1=\text{Γ}{\left(\frac{5}{4}\right)}^{\frac{3}{2}}\text{Γ}{\left(\frac{3}{4}\right)}^{-\frac{5}{2}}, {\beta }_2=\text{Γ}\left(\frac{5}{4}\right)\text{Γ}{\left(\frac{3}{4}\right)}^{-1}$ , $\varphi$ is the mean electron energy, and ε is the electron energy.

The EEDF can be calculated from the Boltzmann equation, which can be solved through a drastic simplification.

where f is the electron distribution in six-dimensional phase space, v is the velocity coordinates, and C[f] is the rate of change in f due to collisions. The EEDF obtained by solving the Boltzmann equation in the two-term approximation was F ₀ [11].

The plasma parameters (rate coefficients and reduced transport properties) were calculated from the EEDF as follows:

where $\gamma$ is ${\left(\frac{2q}{m_e}\right)}^{1/2}$ , ε is the energy, $\overline{\epsilon }$ is the mean electron energy, ${\sigma }_j$ is the collision cross section, and ${\tilde{\sigma }}_m$ is the normalized total momentum collision cross section. Rate coefficients and transport properties were used in each source term of the continuity equation and flux calculation.

The following equation was applied to consider the gas velocity and pressure. The momentum conservation equation is given as:

and the continuity equation as:

where ρ is the gas density, $\overrightarrow{u}$ is the gas flow velocity, p is the pressure, and $\text{μ}$ is the dynamic viscosity. The following equation was applied to consider the thermodynamic properties.

3.1. Plasma parameters

This section is the result of the Ar plasma discharge calculated using the Druyvesteyn EEDF at 10 mTorr and absorbed power of 570 W, considering gas flow and heat transfer. A comparison between the simulation values shown in Fig. 3 and the experimental results shown in [7] indicates a difference in electron density; however, the tendency is consistent. Figures 4(a)–(d) show the electron density, electron temperature, plasma potential, and absorbed power, respectively. The electron density and plasma potential appeared to have the highest values at the center of the plasma-generating section. The electron temperature had a maximum value near the bottom of the antenna, where the power deposition was maximum.

Figure 3. Simulation results of plasma density for Ar 10 mTorr and a power of 570 W; (a) axial electron density profiles and (b) radial electron density profiles.

Figure 4. Plasma properties computed for Druyvesteyn EEDF at 10 mTorr and gas flow rate of 100 sccm, and absorbed power of 570 W; (a) electron density (1/m³), (b) electron temperature (eV), (c) plasma potential (V), and (d) power density (W/m³).

3.2. Comparison of plasma properties according to EEDF

Figure 5 shows the change in EEDF at 1 mTorr gas pressure, 100 sccm inflow gas flow rate, and 24 A coil current. When the Maxwell EEDF was used, it can be seen that the electron temperature, plasma potential, and absorbed power were lower than those obtained from the other EEDFs. As shown in Fig. 2, Maxwell EEDF has a large number of high energy electrons. Thus, the number of electrons involved in the ionization reaction increases, and the electron energy after the reaction decreases, resulting in a lower electron temperature. This was also the reason for the higher electron temperature obtained from the Druyvesteyn EEDF, which had a low energy maximum. When a coil current of 24 A was applied, the absorbed power in the Maxwell EEDF was 350 W, which was different from the 500 W obtained from the other EEDFs. EEDF is applied to the calculation of electron energy mobility, diffusivity, and reaction rate, and it was found that the effect of this on the simulation result.

Figure 5. Plasma properties distribution of electron density, electron temperature, plasma potential, power density with (a) Maxwell EEDF, (b) Druyvesteyn EEDF, and (c) Boltzmann EEDF.

3.3. Plasma properties under pressure

Figure 6 shows the simulation results for various pressures when z = 9 cm. For the Druyvesteyn EEDF, the gas flow rate used is 100 sccm, the pressure conditions are 1, 10, and 25 mTorr, and the absorption power is 500 W. As the pressure increased, the electron density increased while the electron temperature decreased, as shown in Fig. 6. When the pressure increased, the density of Ar colliding with electrons and the ionization reaction also increased. Thus, the electron density increased while the electron temperature decreased. As the pressure decreased, the density of Ar decreased and the time for the electrons to collide with the neutral particles increased. Thus, ion mobility increased, and particles quickly exited the chamber. By considering the temperature and mobility of the ions, it was possible to perform stable plasma analysis even at a pressure as low as 1 mTorr.

Figure 6. Radial plasma properties profiles (z = 9 cm). Distributions of (a) electron density, (b) electric potential, (c) Ar density, and (d) Ar+ density with pressures (1, 10, and 25 mTorr).

In this study, an equipment for Ar ICP for semiconductor manufacturing processes was modeled. The simulation results were analyzed based on the EEDF and pressure levels applied. Gas flow and heat transfer effects were applied to create conditions similar to those in experiments. To compare with the experimental values, the structure of an ICP equipment published in the literature was used. The comparison of plasma properties was analyzed using the Maxwell, Druyvesteyn, and Boltzmann EEDFs. The results of the model assuming the Maxwell EEDF were different from those when applying the other EEDFs. As the EEDF was used in the calculation of plasma parameters, it had a significant influence on the results. In addition, the electron density, electron temperature, and plasma potential were measured based on the applied pressure and their correlations were analyzed. As the pressure decreased, particle mobility increased, and the particles tended to quickly escape from the chamber. Specifying the temperature and mobility of ions resulted in a stable plasma density at a pressure of 1 mTorr. Including an analysis on the behavior of ions in the analysis of the phenomenon of low pressure plasma discharge was found important.

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