COMSOL 6.4 - GEC ICP Reactor, Argon Chemistry

The GEC cell was introduced by NIST (National Institute of Standards and Technology) in order to provide a standardized platform for experimental and modeling studies of discharges in different laboratories. The plasma is sustained via inductive heating. The Reference Cell operates as an inductively-coupled plasma in this model.

Figure 1: GEC ICP reactor geometry consisting of a 5 turn copper coil, plasma volume, dielectrics, and wafer with pedestal.

This application requires the Plasma Module and the AC/DC Module.

Model Definition

Electron transport is modeled by solving the continuity equation, the momentum equation under the drift-diffusion approximation, and the mean electron energy equation (for detailed information on electron transport, see Theory for the Drift Diffusion Interface in the Plasma Module User’s Guide)

The source coefficients in the above equations are determined by the plasma chemistry. The electron rate expression is defined as

where νe,j is the stoichiometric coefficient, and the reaction rate is defined as

where kjf is the forward rate constant and kjr is the reversed rate constant. Both the Electron Impact Reaction feature and Reaction feature can contribute to the electron rate expression. However, when using the Reaction feature it is important to note that the associated electron energy gain or loss is not included in the source term of the electron mean energy equation.

The rate constants can be computed from electron impact cross-section data

where γ = (2q/me)1/2 (SI unit: C1/2/kg1/2), me is the electron mass (SI unit: kg), ε is the electron energy (SI unit: V), σ is the electron impact collision cross section (SI unit: m2), and f is the electron energy distribution function.

When Townsend coefficients are used, the reaction rate is defined as

where αj/Nn is the reduced Townsend coefficient for reaction j (SI unit: m2) and Γe is the electron flux as defined above (SI unit: 1/(m2·s)). Townsend coefficients can increase the stability of the numerical scheme when the electron flux is field driven as is the case with DC discharges.

The total electron energy loss or gained is calculated by summing the collisional energy changes from all reactions defined with the Electron Impact Reaction feature as

where Δεj is the energy loss from reaction j (SI unit: V) and F is the Faraday constant (SI unit: C/mol). For excitation and ionization collisions Δεj corresponds to the energy of the excited state being excited/deexcited or ionized, for attachment Δεj is set to zero, and for elastic collisions

where me and mk are the electron and heavy species mass in kg, Te is the electron temperature in eV, and Tgas is the gas temperature in K.

For heavy species, the following equation is solved for the mass fraction of each species (for detailed information on the transport of the nonelectron species, see Theory for the Heavy Species Transport Interface in the Plasma Module User’s Guide):

The electrostatic field is computed using the following equation:

The space charge density ρ is automatically computed based on the plasma chemistry specified in the model using the formula

For detailed information about electrostatics see Theory for the Electrostatics Interface in the Plasma Module User’s Guide.

For a nonmagnetized, nonpolarized plasma, the induction currents are computed in the frequency domain using the following equation:

The electromagnetic wave “sees” a plasma defined by the plasma conductivity in the cold plasma approximation that is set in the multiphysics feature:

where ne is the electron density, q is the electron charge, me is the electron mass, νe is the collision frequency, and ω is the angular frequency. The Joule heating term that is responsible to heat the electrons is set in the multiphysics feature.

boundary conditions

Electrons are lost to the wall due to random motion within a few mean free paths of the wall and gained due to secondary emission effects, resulting in the following boundary condition for the electron flux:

and the electron energy flux:

For the heavy species, ions are lost to the wall due to surface reactions and the fact that the electric field is directed toward the wall:

The walls of the reactor are grounded.

plasma chemistry

Because the physics occurring in an inductively coupled plasma is rather complex, it is always best to start a modeling project with a simple chemical mechanism. Argon is one of the simplest mechanisms to implement at low pressures. The electronically excited states can be lumped into a single species, which results in a chemical mechanism consisting of only 3 species and 7 reactions (electron impact cross sections are obtained from Ref. 3):

Table 1: table of collisions and reactions modeled.
reaction	formula	type
1	e+Ar=>e+Ar	Elastic	0
2	e+Ar=>e+Ars	Excitation	11.5
3	e+Ars=>e+Ar	Superelastic	-11.5
4	e+Ar=>2e+Ar+	Ionization	15.8
5	e+Ars=>2e+Ar+	Ionization	4.24
6	Ars+Ars=>e+Ar+Ar+	Penning ionization	-
7	Ars+Ar=>Ar+Ar	Metastable quenching	-

Stepwise ionization (reaction 5) can play an important role in sustaining low pressure argon discharges. Excited argon atoms are consumed via superelastic collisions with electrons, quenching with neutral argon atoms, ionization or Penning ionization where two metastable argon atoms react to form a neutral argon atom, an argon ion and an electron. In addition to volumetric reactions, the following surface reactions are implemented:

Table 2: table of surface reactions.
reaction	formula	Sticking Coefficient
1	Ars=>Ar	1
2	Ar+=>Ar	1

When a metastable argon atom makes contact with the wall, it reverts to the ground state argon atom with some probability (the sticking coefficient).

electrical excitation

From an electrical point of view, the GEC reactor behaves as a transformer. A current is applied to the driving coil (the primary) and this induces a current in the plasma (the secondary). The plasma then induces an opposing current back in the coil, increasing its resistance. The current flowing in the plasma depends on the current applied to the coil and the reaction kinetics. The total plasma current can vary from no current (plasma not sustained) to the same current as the primary which corresponds to perfect coupling between the coil and the plasma.

In this example a fixed power of 1500 W is applied to the coil. Some of this power is dissipated in the coil, some is deposited into the plasma.

Results and Discussion

The peak electron density occurs at the center of the reactor, underneath the RF coil. The electron density in this case is high enough to cause some shielding of the azimuthal electric field.

Figure 2: Plot of the electron density inside the GEC ICP reactor.

The electron “temperature” is highest directly underneath the coil, which is where the bulk of the power deposition occurs.

Figure 3: Plot of the electron “temperature” inside the GEC ICP reactor.

Figure 4: Plot of the electric potential inside the GEC ICP reactor.

From an electrical standpoint, the quantities of interest are total power deposition, coil resistance and inductance, and reactor efficiency. These “global” parameters are relatively easy to measure when the plasma is on or off, so such quantities provide an easy route for comparison with experimental data, without the need for expensive optical emission spectroscopy equipment or Langmuir probes.

The resistance of the coil increases by a little less than a factor of 4 when the plasma is on. When the plasma is on, there is a substantial opposing current induced back into the coil from the plasma. The electric potential applied across the coil needs to increase in order to maintain the same total current.

Figure 5: Plot of coil resistance versus time in GEC ICP reactor.

Initially the power dissipated is all dissipated in the coil (~1400 W). After about 1 microsecond, the plasma ignition begins and as the neutral gas atoms split into electrons and ions, the electrons begin to absorb more and more power. Over a period of 2 microseconds, the plasma goes from absorbing no power to absorbing around 1200 W.

Figure 6: Plot of total power versus time in GEC ICP reactor.

The ion density is exactly the same except from in a thin region close to the walls. In this region, the ion density dominates the electron density which leads to a positive potential in the plasma bulk with respect to the walls. The positive potential increases the flux of ions and reduces the flux of electrons to the wall.

Figure 7: Plot of number density of argon ions in the GEC ICP reactor.

Figure 8: Plot of the norm of the electric field due to the induction currents.

The number density of excited species is also greatest in the center of the reactor. Unlike the charged species, there is no rapid drop off in number density close to the walls. The physics of the excited species is relatively simple: they are formed in the center of the reactor by high energy electrons and are lost to the via either stepwise ionization or diffusion to the wall. Because the excited argon atoms are not susceptible to migration due to the electric field, they can exist in much higher quantities than ions. The peak number density of excited argon atoms represents a mass fraction of around 0.02.

Figure 9: Plot of the number density of excited argon atoms in the GEC ICP reactor.

The skin depth of the plasma is on the order of 1cm which prevents the electric field from penetrating into the core of the plasma. The skin depth is defined as:

where μ is the permeability, σ is the plasma conductivity, and ω is the angular frequency. This tells us that increasing the driving frequency does not necessarily couple more power into the plasma. As the frequency increases, the plasma tends to shield the region over which power is deposited into a thin layer close to the upper wall.

Figure 10: Plot of the power deposition into the plasma in the GEC ICP reactor. The region over which power is deposited to the plasma is governed by the plasma skin depth.

References

1. G.J.M. Hagelaar and L.C. Pitchford, “Solving the Boltzmann Equation to Obtain Electron Transport Coefficients and Rate Coefficients for Fluid Models,” Plasma Sources Sci. Technol., vol. 14, pp. 722–733, 2005.

2. D.P. Lymberopoulos and D.J. Economou, “Two-Dimensional Self-Consistent Radio Frequency Plasma Simulations Relevant to the Gaseous Electronics Conference RF Reference Cell,” J. Res. Natl. Inst. Stand. Technol., vol. 100, pp. 473–494, 1995.

3. Phelps database, www.lxcat.net, retrieved 2017.

Plasma_Module/Inductively_Coupled_Plasmas/argon_gec_icp

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Geometry 1

Insert the prepared geometry sequence from file. You can read the instructions for creating the geometry in the appendix.

In the toolbar, click and choose .

Browse to the model’s Application Libraries folder and double-click the file argon_gec_icp_geom_sequence.mph.

Add some predefined selections for the geometric entities which will be referenced later on.

Definitions

Walls

In the toolbar, click

In the window for , type Walls in the text field.

Locate the section. From the list, choose .

Select Boundaries 6, 8, 35–38, 44, 45, and 51–56 only.

Coils

In the toolbar, click

In the window for , type Coils in the text field.

Select Domains 6 and 8–11 only.

Coil Boundaries

In the toolbar, click

In the window for , type Coil Boundaries in the text field.

Select Domains 6 and 8–11 only.

Locate the section. From the list, choose .

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
Psp	1500[W]	1500 W	Power input
mueN	4E24[1/(mVs)]	4E24 1/(V·m·s)	Reduced electron mobility
T0	300[K]	300 K	Gas temperature
p0	0.02[torr]	2.6664 Pa	Gas pressure

Plasma (plas)

In the window, under click .

Select Domain 3 only.

Cross Section Import 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file Ar_xsecs.txt.

Click

In the window, click .

In the window for , locate the section.

Select the checkbox.

Now you add two more regular reactions which describe how electronically excited argon atoms are consumed on the volumetric level. The rate coefficients for these reactions are taken from the literature.

Reaction 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type Ars+Ars=>e+Ar+Ar+.

Locate the section. In the kf text field, type 3.734E8.

Reaction 2

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type Ars+Ar=>Ar+Ar.

Locate the section. In the kf text field, type 1807.

When solving any type of reacting flow problem there always needs to be one species which is selected to fulfill the mass constraint. This should be taken as the species with the largest mass fraction.

Species: Ar

In the window, click .

In the window for , locate the section.

Select the checkbox.

When solving a plasma problem the plasma must be initially charge neutral. COMSOL automatically computes the initial concentration of a selected ionic species such that the initial electroneutrality constraint is satisfied. Once the simulation begins to time step, the plasma need not be charge neutral. In fact, the separation of space charge between the ions and electrons close to the wall is a critical component in sustaining the discharge.

Species: Ar+

In the window, click .

In the window for , locate the section.

Select the checkbox.

Initial conditions for the electron number density and mean electron energy are critical for any plasma model. If the initial electron density is too low then the plasma may not be able to sustain itself and may self extinguish. If the initial electron density is too high then convergence problems may occur during initial time steps.

Initial Values 1

In the window, click .

In the window for , locate the section.

In the ne,0 text field, type 1E15.

In the ε0 text field, type 5.

The are computed everywhere except the wafer and the wafer pedestal.

Magnetic Fields (mf)

In the window, under click .

Select Domains 3–6 and 8–12 only.

The feature is used to electrically excite the system. The coil operates with a fixed total power of 1500 watts.

Domain Coil 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Select the checkbox.

In the Pcoil text field, type Psp.

Plasma (plas)

Plasma Model 1

In the window, under > click .

In the window for , locate the section.

In the T text field, type T0.

In the pA text field, type p0.

Locate the section. In the μeNn text field, type mueN.

Next, define the material properties. There is no need to define the material properties in the plasma domain, as these are defined by the feature.

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	6e7	S/m	Basic
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

Material 2 (mat2)

Right-click and choose .

Select Domain 5 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	0	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

Material 3 (mat3)

Right-click and choose .

Select Domains 4 and 12 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	0	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	4.2	1	Basic