An inductively coupled plasma (ICP) discharge is used widely for wafer fabrication processes of microelectronic devices such as semiconductors or displays [1]. Further, an ICP has advantages on plasma-surface interactions such as etching or deposition with an elaborate scale in the manufacturing of micro-circuits because of the high density, low pressure, and good uniformity [2-8].

To increase the efficiency of ultrafine processes, it is necessary to gain an understanding of the physics of plasma through computational modeling, which is an appropriate approach to determine associations between plasma factors and process conditions that are difficult to detect in experiments. The ICP discharge simulations are conducted in a two-dimensional axisymmetric model because of the limited computing resources and time [9]; however, in this study, a three-dimensional asymmetric model is used to investigate the asymmetric effects in plasma properties caused by the geometrical structure, and to simulate the actual process conditions more closely.

In this study, a fluid-based model was used among various plasma analysis methods for the simulation and detailed explanations about this model are present in the literature [10,11]. The fluid model treats plasma particles in the chamber as a fluid, and thus, there are limitations because plasma phenomena occur owing to the movement of individual particles in low-pressure environments. However, fluid model has the benefits of saving resources by using the averaged statistical characteristics and assuming energy distribution of particles to simplify equations [12].

The power flowing through the antenna input generates a magnetic field, and the induced electric field transfers the force to the electrons in the plasma region [13]. The magnitude and distribution of the coil voltage vary based on the impedance of the antenna; further, they are closely related to the magnitude of the transmitted force, and in fact, the plasma characteristics change based on this factor. Among these impedances calculated by field, capacitance and inductance play a considerable role in calculating voltage. Electric energy―the source of capacitance―depends on the electromagnetic field, and therefore, it is important to solve both the inductive and capacitive electric fields to calculate magnetic and electrical physics accurately.

In this paper, a three-dimensional simulation was carried out to consider the asymmetric effect of the plasma characteristics along with the gas inlet and the antenna shape, which could not be considered in the two-dimensional axisymmetric structure. Fluid-based simulations were performed in COMSOL Multiphysics. This study focuses on the characteristics of the plasma depending on the coil impedance variation and the power induced in the chamber according to the method of analyzing the electromagnetic field

2.1. Model description

In this study, a three-dimensional asymmetric model was used, and the geometric scheme of this model is as shown in Fig. 1. There are four plasma chamber areas: plasma region, ceramic window, antenna, and air cap. The diameter of the plasma region is 35 cm, and each dielectric portion and outlet is located at a thickness of 2.5 cm at the bottom thereof. Four inlet ports, from which the reaction gas flows in, are located at 90° intervals on the side of the chamber. To separate the plasma region from the antenna, a ceramic window is placed on top of the plasma region at a height of 2 cm and they are composed of aluminum oxide. The antenna, which consists of copper, is wound three times with a cross section of 1 cm2 and applied 100 A at a frequency of 13.56 MHz. A shielding cap filled with air is set on the ceramic window to protect the coil; rectangular floating portions around the power supply part of the antenna are designed to avoid the effect of the upper part of the shielding cap that is set to ground.

Figure 1. Schematic of the ICP model.

Argon was used as the reaction gas under the condition of 20 mTorr; Table I summarizes the reactions used in the simulation process. This study was conducted with COMSOL Multiphysics using a fluid model.

Table 1 . Important collision processes in the argon discharge.

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

2.2. Equations

The equations for electron density and electron energy density, which follows a drift-diffusion approximation, are provided below.

where n_e is the electron density, ${\overrightarrow{\Gamma }}_e$ is the electron flux, (S-L)_e is the difference between source of electron and loss of electron, μ_e is the electron mobility, and 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, and S_en is the energy loss and gain due to the collisions, μ_ε is the electron energy mobility; and D_ε is the electron energy diffusivity.

The gas inlets and outlets need to satisfy the conditions of zero charge on the boundary, as indicated in Eq. (6). In addition, to apply the insulation effect, the normal components of the electron and electron energy flux must be zero at the boundary, as indicated in Eqs. (7) and (8).

The lower surface of the ceramic window in contact with the plasma region must meet the surface charge accumulation condition, which builds up charge in the normal direction, as

where $\overrightarrow{D}$ is the electric displacement field, ρ_s is the surface charge density and it is the solution of the Eq. (10), ${\overrightarrow{J}}_i$ is the total ion current density on the wall, and ${\overrightarrow{J}}_e$ is the total electron current density on the wall.

The following equations explain the mixture averaged method for calculating heavy species, which satisfies mass conservation [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, Z is the charge number for the species, and μ_m is the mixture averaged mobility for species.

While formulating the electric field, the magnetic field (MF) physics considers the inductive field by considering magnetic vector potential as shown below.

where $\overrightarrow{H}$ is the magnetic field, $\overrightarrow{J}$ is the current density, $\overrightarrow{B}$ is the magnetic flux density, and $\overrightarrow{A}$ is the magnetic vector potential.

Meanwhile, the electric field calculated in magnetic and electric field physics (MEF), which considers both of inductive field and capacitive field by using the electric potential additionally, is given by

where V is the electric potential. The equations for computing the magnetic flux density and magnetic field are same as those for MF physics.

To prove the physics setting that can calculate the inductance and impedance of the antenna accurately, three cases were considered. Case (1) was simulated with a coil analysis node in MF physics. Cases (2) and (3) were simulated with MEF physics; however, the difference is that coil analysis was used only in Case (2).

The expression of coil inductance and capacitance used in Cases (1) and (2) is given by

and the voltage equation according to the impedance change is defined based on the phasor form of Ohm’s law as

where L_coil is the inductance of the coil, W_m is the total magnetic energy, and I_coil is the current that imposed at coil; C_coil is the capacitance of the coil, W_e is the total electric energy, and V_coil is the voltage of the coil; and Z_coil is $R_{\mathrm{coil}}+j({\omega L}_{\mathrm{coil}}-\frac{1}{{\mathrm{\omega c}}_{\mathrm{coil}}})$ , where $\omega$ is the angular frequency.

Meanwhile the impedance is calculated from a field definition in Case (3) by using Poynting’s theorem, as

where w_e is the electric energy density, w_m is the magnetic energy density, R is the resistance, I_i is the input current, $\overrightarrow{H}$ is the magnetic field, $\overrightarrow{D}$ is the electric flux density, s_i is the input surface, cross section of the coaxial line, $\overrightarrow{S}$ is the Poynting vector, and X is the reactance.

3.1. Case comparison focused on the electromagnetic field analysis

First, the resistance, inductance, and capacitance of the coil calculated by the electromagnetic field and the coil voltage because of the impedance change for all cases are listed in Table II. Further, the power absorbed by electrons because of inductive discharge is listed for when the current applied to the antenna is 100 A. The absorbed power contributes to computing electron energy density (n_ε ).

Table 2 . Comparison of coil impedance, voltage, and absorbed power between all cases.

	Case (1)	Case (2)	Case (3)
Rcoil [Ω]	0.39	0.42	0.44
Lcoil [H]	3.6 × 10-7	4.2 × 10-7	4.1 × 10-7
Ccoil [F]	−7.1 × 10-15	−3.83 × 10-15	−5.7 × 10-12
Vcoil [V]	39.5 + j6273	42.1 + j7162	44.38 + j6667
Pabs [W]	1977	2103	2195

In coil geometry analysis, the voltage out-of-plane to the coil was the source of the current density. This current density was used to calculate the reactance of the coil. Therefore, in this coil geometry analysis, the calculation method was suitable for simple static materials; however, the plasma was very complex for applying this method.

Figure 2 shows the real and imaginary parts of the antenna voltage in Cases (2) and (3). It is common for the power injection port to show the highest numerical value (red) and the output part connected to the ground (blue) with the lowest value distribution as shown in Case (3). However, there is an absence of settings for the ground connection to the antenna output part in the coil analysis node. Thus, Case (2) shows a nonzero potential in the output section, which causes a numerical error. In addition, the coil geometry analysis study step, in which the current direction of the antenna is solved, produces optimized results when the coil has a constant cross-section without shape change or bottleneck phenomena. Therefore, the current calculation at the coil bending part could be incorrect.

Figure 2. Coil voltage distribution for Cases (2) and (3).

Figure 3(a) shows the magnetic flux density of the magnetic field generated by the power applied to the coil. Eq. (14) indicates that the magnetic field can be expected to be generated in the rotating direction of the current flow. Further, Fig. 3(b) represents the electric field induced by the magnetic field transmitted to the plasma region. The power generated from the coil is transferred to electrons in the plasma via a time-varying magnetic field, which leads to an azimuthal direction induced electric field that generates high-density plasma.

Figure 3. Magnetic flux density and induced electric field for Case (3).

3.2. Case comparison focused on the plasma characteristics

Figure 4 illustrates normalized plasma characteristics with varying transmitted power caused by coil voltage variations. As shown in Fig. 4, there are slightly differences in the maximum value of electron density, electron temperature, and electric potential between all cases. Since Case (3) received the most power in the plasma region, it has the highest electron density and temperature compared to the other cases. Further, because the time to reach the plasma wall changes based on the speed difference caused by the differences in weight between electrons and ions, the electrons collide into the wall faster than the ions. Thus, the electrons are lost faster than ions, and the inside of the plasma region has a positive charge because it repeats the light emission process of generating electrons and ions. Negative potential near the wall accelerates ions and slows electrons because of the repulsion. The same number of electrons and ions reach the wall, and thus, the wall has a negative charge until the electron current density and ion current density are balanced, which is called the sheath region. Therefore, the electric potential has a negative value at the side of the chamber.

Figure 4. Comparison of electron density, electron temperature, and electric potential between all cases.

Figures 5 - 8 show the major plasma characteristics by sliced on xy plane at 9 cm where the inlets are located, 5 cm as the center of the chamber, and 1 cm located symmetrically with the inlets relative to the center of the chamber. Since Case (3) has the highest absorbed power, as shown in Table II, exhibits the highest values for collisional power loss and electron density. In addition, due to the structural asymmetry, the distribution of variables is skewed to one side that cannot be observed in the two-dimensional axisymmetric structure.

Figure 5. Comparison of electron density on sliced xy plane.

Figure 6. Comparison of electron temperature on sliced xy plane.

Figure 7. Comparison of electric potential on sliced xy plane.

Figure 8. Comparison of collisional power loss on sliced xy plane.

This study observed the appropriate physics settings required for calculating the correct coil voltage by determining the optimal coil impedance using COMSOL Multiphysics. The simulation is conducted using a three-dimensional asymmetric ICP discharge model for semiconductors manufacturing process, and this work is based on fluid models for plasma. Owing to the difference in the calculation of the electric field, magnetic field physics and magnetic and electric field physics show large differences in the total electric energy that leads to variances in impedance. This shows that accurate electric field calculations are essential because they are closely related to the power induced to the plasma chamber. In addition, due to the absence of a setting to apply ground to the output of the antenna and the diversified geometric characteristics of the antenna structure, the use of coil geometry analysis node results in incorrect current calculations which leads to a numerical error. Therefore, accurate impedance calculations computed from the electromagnetic field based on Poynting’s theorem increase the simulation accuracy.

To obtain more suitable plasma properties, future studies will include heat transfer physics that considers changes in enthalpy caused by collisions and flow physics to consider the change in pressure and velocity fields along the flow of fluid.

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