3D ICP Reactor, Argon Chemistry

Although not recommended, 3-dimensional plasma modeling is possible to do in COMSOL. A square coil is placed on top of a dielectric window and it electrically excited at 13.56 MHz. A plasma is formed in the chamber beneath the dielectric window, which contains argon gas at low pressure (20 mtorr). The gas flows into the process chamber from two 2 inch ports and the gas is extracted through a single 4 inch port. The plasma is sustained via electromagnetic induction where power is transferred from the electromagnetic fields to the electrons.

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

Model Definition

Inductively coupled discharges typically operate at low pressures (<10 Pa) and high charge density (>1017 m−3). High density plasma sources are popular because low pressure ion bombardment can provide a greater degree of anisotropy on the surface of the wafer.

Domain equations

The Inductively Coupled Plasma (ICP) reactor interface solves a system of coupled partial differential equations for:

•

The electron density.

•

The mean electron energy.

•

The mass fraction of each of the non-electron species.

•

The electrostatic potential.

•

The electric field due to induction currents.

The electron density and mean electron energy are computed by solving a pair of drift-diffusion equations for the electron density and mean electron energy. Convection of electrons due to fluid motion is neglected. For detailed information on electron transport see the Plasma Module User’s Guide.

The electron source Re and the energy loss due to inelastic collisions Rε are defined later. In this model, the electron diffusivity, energy mobility and energy diffusivity are calculated from the electron mobility, μe using:

The source coefficients in the above equations are determined by the plasma chemistry and are written using either rate coefficients. Suppose that there are M reactions which contribute to the growth or decay of electron density and P electron-neutral collisions. In general P >> M. The electron source term is given by

where xj is the mole fraction of the target species for reaction j, kj is the rate coefficient for reaction j (m3/s), and Nn is the total neutral number density (1/m3). The electron energy loss is obtained by summing the collisional energy loss over all reactions:

where Δεj is the energy loss from reaction j (V). The electron source and inelastic energy loss are automatically computed by the interface. The rate coefficients may be computed from cross section data by the following integral:

where γ = (2q/me)1/2 (C1/2/kg1/2), me is the electron mass (kg), ε is energy (V), σk is the collision cross section (m2) and f is the electron energy distribution function (EEDF). In this model a Maxwellian EEDF is assumed.

For non-electron species, the following equation is solved for the mass fraction of each species. For detailed information on the transport of the non-electron 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 plasma conductivity is automatically defined using the cold plasma approximation:

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.

plasma chemistry

Since 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-section 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 impedance. The current flowing in the plasma depends on the current applied to the coil and the reaction kinetics. The total plasma current may 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. Since this is a 3D model it is not practical to resolve the skin depth in the driving coil. Instead an impedance boundary condition is used which only computes the coil current density in the tangential direction of the coil surface. This means the mesh does not have to resolve the skin depth in the coil, which at 13.56MHz for copper is only 17 microns.

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 1: Plot of electron temperature (slice), magnetic vector potential (streamlines) and electric field norm (streamline color).

The electron “temperature” is highest directly underneath the coil, which is where the bulk of the power deposition occurs. The streamlines in Figure 1 show the norm of the electric field. The streamlines aren’t perfectly circular since the coil is asymmetric.

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. Lymberopolous 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, p. 473, 1995.

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

Plasma_Module/Inductively_Coupled_Plasmas/argon_3d_icp

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

In the tree, select .

button in the toolbar.

In the dialog box, in the tree, select the check box for the node .

Click .

Import 1 (imp1)

In the toolbar, click

In the window for , locate the section.

Click .

Browse to the model’s Application Libraries folder and double-click the file argon_3d_icp.mphbin.

Click .

Definitions

Coil boundaries

In the toolbar, click

Click the

button in the toolbar.

Select Domains 4–6, 8, and 9 only.

Right-click and choose .

In the dialog box, type Coil boundaries in the text field.

Click .

In the window for , locate the section.

From the list, choose .

Current density

In the toolbar, click

Select Domain 7 only.

In the window for , locate the section.

From the list, choose .

Right-click and choose .

In the dialog box, type Current density in the text field.

Click .

Plasma

In the toolbar, click

Select Domain 1 only.

Right-click and choose .

In the dialog box, type Plasma in the text field.

Click .

Inlets

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundaries 1 and 4 only.

Right-click and choose .

In the dialog box, type Inlets in the text field.

Click .

Outlet

In the toolbar, click

Select Domain 1 only.

In the window for , locate the section.

From the list, choose .

Select Boundary 100 only.

Right-click and choose .

In the dialog box, type Outlet in the text field.

Click .

Click the

button in the toolbar.

ICP domains

In the toolbar, click

Select Domains 1–3 only.

Right-click and choose .

In the dialog box, type ICP domains in the text field.

Click .

Plasma walls

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundaries 2, 3, 5–9, 12, 17–23, 70, 73–75, and 96–99 only.

Right-click and choose .

In the dialog box, type Plasma walls in the text field.

Click .

In the window for , locate the section.

Clear the check box.

Click the

button in the toolbar.

Variables 1

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Mw	0.04[kg/mol]	kg/mol	Molecular weight
T	300[K]	K	Model temperature
p0	0.02[Torr]	Pa	Outlet pressure
rho_stp	101325[Pa]*Mw/R_const/273.15[K]	kg/m³	Density at standard pressure and temperature
Qin	0.5E-2[m^3/s]	m³/s	Volumetric flowrate
Js0	1E3[A/m]tanh(1E7t[1/s])	A/m	Surface current density

Plasma (plas)

In the window, under click .

In the window for , locate the section.

Click

Select Domain 1 only.

Magnetic Fields (mf)

In the window, under click .

Select Domains 1–3 only.

Laminar Flow (spf)

Since the density variation is not small, the flow cannot be regarded as incompressible. Therefore set the flow to be compressible.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. Click

Select Domain 1 only.

Define the pressure reference level in the interface properties.

Locate the section. In the pref text field, type 0.

Plasma (plas)

In the window, under click .

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.

Find the subsection. Select the check box.

Species: Ar

In the window, click .

In the window for , locate the section.

Select the check box.

Species: Ar+

In the window, click .

In the window for , locate the section.

Select the check box.

Species: Ars

In the window, click .

In the window for , locate the section.

In the x0 text field, type 1E-4.

Initial Values 1

In the window, click .

In the window for , locate the section.

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

In the ε0 text field, type 4.

Plasma Model 1

In the window, click .

In the window for , locate the section.

In the μe text field, type 4E24[1/(m*V*s)]/plas.Nn.

Locate the section. From the u list, choose .

In the T text field, type T.

From the pA list, choose .

Magnetic Fields (mf)

In the window, under click .

Impedance Boundary Condition 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Surface Current Density 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Specify the Js0 vector as


Js0	x
0	y
0	z

Plasma (plas)

In the window, under click .

Ground 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Wall 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Surface Reaction 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type Ar+=>Ar.

Surface Reaction 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type Ars=>Ar.

Species: Ars

In the window, click .

Outflow 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Species: Ar+

In the window, click .

Outflow 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Laminar Flow (spf)

In the window, under click .

In the window for , locate the section.

In the pref text field, type p0.

Fluid Properties 1

In the window, under click .

In the window for , locate the section.

From the ρ list, choose .

From the μ list, choose .

Inlet 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 .

Locate the section. Click the button.

In the V0 text field, type Qin.

Outlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Materials

Material 1 (mat1)

In the window, under right-click and choose .

Select Domains 1–3 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
Electrical 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 2 (mat2)

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
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Electrical conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	6e7	S/m	Basic