COMSOL 6.4 - Model of an Argon/Chlorine Inductively Coupled Plasma Reactor with RF Bias

Model of an Argon/Chlorine Inductively Coupled Plasma Reactor with RF Bias

This tutorial uses the interface to model an argon/chlorine inductively coupled plasma reactor with RF bias (also known as ICP/CCP). The interface is used to study discharges that are sustained by induction currents and have an additional time-periodic electric excitation. This interface adds the and interfaces. The magnetic field is solved in the frequency domain and the plasma solves for a periodic solution. The multiphysics couplings couple the time-averaged plasma conductivity to the interface and couple the resulting electron heating due to the induction currents back to the interface.

It is also shown how to prepare a model with a mixture of different elements (in this case Ar and Cl2) in which one of the species can dissociate by electron impact (Cl2 dissociates into Cl) and where negative ions exist (the dissociative electron attachment of Cl2 creates Cl-).

A simplified plasma chemistry is used to discuss the main aspects of such discharges. It is important to keep in mind that a benchmark is not attempted and the idea is to provide a base case that can be used to develop more complex chemistries. In fact, it might be necessary to modify the data used and to add more reactions to achieve experimental verification.

Model Definition

The plasma model is solved self-consistently with the , , and interfaces. The plasma transport equations coupled with Poisson’s equation are solved in a time-periodic fashion, allowing for the description of the different quantities along the period. In this way, electrons exhibit periodic modulation caused by the RF bias. The interface is solved in the frequency domain. The Joule heating term obtained from the interface is thus an averaged quantity and is given to the electron mean-energy equation as a constant value along the period. The plasma conductivity used in the magnetic field equations is computed from period-averaged quantities obtained using the interface.

The and interfaces are solved in a stationary form and consequently all period evolution is neglected. The fluid and thermodynamic properties computed in the plasma model and passed on to the and interfaces are period-averaged quantities.

In particular, this model assumes that the heavy species only exist in the . This means that ions and neutrals do not have a period modulation and only experience a time-averaged electric field.

Electric Excitation

There are two independent forms of electric excitation present in this model

•

The inductively coupled power controlled by the coil

•

The capacitively coupled power controlled by the RF bias

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

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

Here, <ne> is the period averaged electron density, q is the electron charge, me is the electron mass, <νe> is the period averaged collision frequency, and ω is the angular frequency. The Joule heating term, which represents the heating of the electrons, is set in the multiphysics feature.

The biased electrode is driven by a sinusoidal excitation at constant voltage amplitude

(1)

The DC bias voltage Vdc,b is computed so that the net current over a period is zero.

Plasma 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 with the exception of the biased electrode.

Plasma Chemistry

Negative ions are created in certain molecular gaseous discharges (like chlorine, oxygen, hydrogen, fluorocarbons, and so on) and these discharges tend to have complex plasma chemistries with many ions, dissociative products, and excited states. Here a simple plasma chemistry is used and no benchmark is attempted. In fact, it might be necessary to modify the data used and add more reactions to achieve experimental verification. Nevertheless, this plasma chemistry allows to show the main aspects of an electronegative discharge. The plasma chemistry should be adapt and improved for specific applications as needed. One aspect that emerges from analyzing model results of the present model is that the dissociation degree is high. Knowing this it might be needed to add more reactions involving Cl.

The plasma chemistry used here includes electron impact reactions from Ref. 1and Ref. 2, and some reactions are estimated from similar reactions based on data from Ref. 3. A good discussion of plasma chemistry for electronegative gases can be found in Ref. 4. In particular, in the section “A Data Set for Oxygen,” page 270 of Ref. 4 a set of reactions for oxygen discharges is discussed.

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 three species and seven reactions presented in Table 1 (electron impact cross sections are obtained from Ref. 1).

Table 1: Argon reactions.
Reaction	Formula	Type
1	e+Ar=>e+Ar	Elastic	-
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	-

Chlorine has a much richer reaction set that includes vibrational and rotational excitations, excitation of several electronic excited states, electron impact dissociation, dissociative attachment, and many others. Electron impact reactions with Cl2 are from Ref. 2 except for e+Cl=>Cl+e+e and e+Cl-=>Cl+e+e, which are estimated from similar reactions from Ref. 3. Electron impact reactions are neglected except for e+Cl=>Cl+e+e. For simplicity, rotational, vibrational, and electronic excited states are not treated explicitly but energy losses are considered. The chlorine electron impact reactions used in this model are presented in Table 2.

Table 2: Chlorine Electron impact reactions.
Reaction	Formula	Type
1	e+Cl2=>Cl+Cl-	Dissociative attachment	-
2	e+Cl2=>e+Cl2	Elastic	-
3	e+Cl2=>e+Cl2	Vibrational excitation	0.069
4	e+Cl2=>e+Cl2	Vibrational excitation	0.139
5-9	e+Cl2=>e+Cl+Cl	Dissociative excitation	3.36-7.02
10	e+Cl2=>e+Cl2	Excitation	10.54
11	e+Cl2=>e+Cl2	Excitation	10.7
12	e+Cl2=>e+Cl-+Cl+	Excitation	11
13	e+Cl2=>2e+Cl2+	Ionization	11.49
14	e+Cl2=>2e+Cl+Cl+	Ionization	11.49
15	e+Cl=>2e+Cl+	Ionization	14.25
16	e+Cl-=>2e+Cl	Ionization	14.25

Table 3 presents the heavy species reactions involving ions. All these reactions rates are estimated from similar reactions from Ref. 3.

Table 3: Heavy species reactions involving ions.
Reaction	Formula	Type
1	Cl++Cl-=>2Cl	Mutual recombination
2	Cl2++Cl-=>3Cl	Mutual recombination
3	Cl-+Ar+=>Cl+Ar	Mutual recombination

In addition to volumetric reactions, Table 4 lists the surface reactions implemented.

Table 4: Surface reactions.
Reaction	Formula	Sticking coefficient	Secondary emission coefficient	Mean energy of secondary electrons (V)
1	Ars=>Ar	1	0.07	5.8
2	Ar+=>Ar	1	0.07	5.8
3	Cl=>0.5Cl2	0.01	0	0
4	Cl2+=>Cl2	1	0.07	5.8
5	Cl-=>Cl	1	0	0
6	Cl+=>Cl	1	0.07	5.8

Boundary conditions for heavy species are introduced in the model by using surface reactions. If no surface reactions that lead to the loss of a given species at a surface are introduced in the model, that species will not have losses by transport. This can lead to the unbounded growth of a given species and a steady-state solution might not be possible.

Atomic recombination (reaction 3 in Table 4) at a surface is an important aspect of plasma discharges with molecular species since it influences the dissociation degree in the discharge. The sticking coefficient for atomic recombination is a function of the surface type and temperature.

Electronegative Plasmas

Electronegative plasmas are plasmas that contain negative ions. Negative ions are mainly created by electron dissociative attachment (for example, e+Cl2=>Cl+Cl-). This reaction tends to be very effective at low electron energies and can reduce the electrons in a discharge to the point that an ion–ion discharge is obtained. The transport and volume creation/destruction mechanisms tend to be more complex than electropositive plasmas in many respects. Here, only a few are mentioned with emphasis on the numerical difficulties that they introduce. More information can be found in Ref. 4, section 10.3 and references therein.

In electronegative discharges negative ions are well confined by the ambipolar electric field and losses by transport are very small. This means that to achieve a steady state volume losses need to be included for negative ions. The mechanisms by which negative ions are lost depend on the gas mixture and pressure and they are: mutual recombination with positive ions (for example, Cl-+Cl+=>2C or Cl-+Ar+=>Cl+Ar), detachment in collisions with excited or neutral atoms or molecules (for example, Cl-+Cl=>Cl2+e or Cl-+Cl2=>Cl+Cl2+e), and electron-impact detachment (for example, e+Cl-=>Cl+2e).

In electronegative discharges it is often possible to identify two spatial regions using the electronegativity (ratio of the negative ion density to the electron density): (i) one in the core of the discharge (the electronegative core) with high electronegativity where the dominant charge species are positive and negative ions; (ii) and the other close to the boundaries (electropositive edges) where the dominant charged species are electron and positive ions. In the transition between these two regions the negative ion density drops abruptly causing a chock-like phenomena. This transition needs to be well resolved spatially. If not, oscillations can be seen in the negative ion density and the model might not converge. Some strategies to deal with this are:

•

Increase the negative ion diffusion coefficient or decreasing its mobility. When using the option the mobility is computed using Einstein relation which by default uses the gas temperature. By specifying a higher ion temperature results in a smaller mobility.

•

Enable in the section (the stabilization sections are visible when is selected in ). This option adds artificial diffusion to all ions and helps smoothing the sharp transition of the negative ion density between the electropositive edge and the electronegative core, and also increase the density of the negative ions in the electropositive edge effectively increasing its losses by transport. This option should be used very carefully since completely wrong results can be obtained if too much diffusion is used (the tuning parameter for ions should not be larger than 0.1). A useful strategy is to start with a large (for example, 0.5) and ramp it down using a .

Inflow and Outflow

The boundary condition fixes the mass fraction or mole fraction of specified species. It is used in this model to set the mole fraction of the inlet mixture at a boundary. The fluid velocity is computed and set at the same boundary by the feature of the interface. At the boundary condition the model fraction of Cl is set to a small number so that the mole fraction of Ar/Cl2 is always respected at the boundary. When dissociation is high, the mole fraction of Ar can be different than specified because the mass fraction of Ar is obtained from a mass constraint and, in fact, no explicit constraint is applied at the boundary.

The boundary condition removes transport by diffusion at the boundary but it still keeps transport by convection and migration in the electric field. The boundary condition only sets boundary conditions for heavy species. To model a net outflow out of the system the boundary condition of the Laminar flow interface should be used.

When solving for plasmas with chemistries that contain more than one element (for example, Ar and Cl2) with a stationary solver the mass fraction of each element is not conserved if no constraint is used. This problem is similar in nature to the one found when solving for Navier–Stokes equations in steady state without fixing the pressure somewhere. As a strategy to overcome this problem, one can fix the mole fraction of a given species using the boundary condition even if no fluid flow exists in the system.

In the present model, it is assumed that the charged species (electrons and ions) are lost by transport at the inflow and outflow boundaries. This is modeled by applying the surface reactions for ions and the feature for electrons. This could represent a metal grid boundary that allows the neutral gas to flow through it, while the plasma “sees” a metal boundary. It is also assumed that Ars and Cl can react at the inflow and outflow boundaries. To adjust the probability of a reaction to occur at a surface, the in the feature can be adjusted. As an example, to model a 50% probability of Ars to react at the outflow grid, the should be set to 0.5.

Solution Strategy

The recommended solution strategy to model an inductively coupled plasma with RF bias is as follows. First solve the problem with inductively coupled plasma only. To speed up this study, set the in the to 1. Since there is no RF bias applied yet it is not necessary to have resolution along the period. Afterward, add a second study that uses the solution of the previous study as initial conditions. Increase the to a value in the range 30–50 and use the to increase the amplitude of the applied voltage.

For pure CCP problems, it is advised to set the to in the feature because it makes the numerical problem much easier. For ICP/CCP models, it is advised to set to .

At the biased electrode, the electrons can exhibit strong period modulation and can attain as high energies as in a pure CCP reactor. For these reasons, the space and period resolutions need to be as high as in a pure CCP reactor model.

Results and Discussion

Figure 1 through Figure 3 show the period-averaged spatial distributions of the electron density, electron temperature, and electric potential for 250 W of coil power and Vrf = 100 V. The electron density is relatively low (almost 1016 m−3) and has a flat profile in the core of the discharge, where the dominant charged species are the ions of chlorine (not shown) with number densities of the order of 1017 m−3. This makes the discharge have high electronegativity in most of the reactor. At the edges, the negative ion density decreases fast and the plasma becomes electropositive.

The dominant power coupled to the system is inductive and is responsible for maintaining a high plasma density. The RF bias influences the plasma density significantly but the electron temperature changes from having a maximum just below the coil for Vrf = 0 V to a maximum at the biased electrode for higher voltages.

From the period-averaged potential it is possible to see that a DC bias of about −75 V is developed at the biased electrode. This is also possible to see from the period evolution of the applied voltage Figure 4. This figure also shows the current at the electrode, which strongly deviates from a sinusoidal and has a strong contribution from the electron current when the sheath collapses.

Figure 5 and Figure 6 show the magnitude of the fluid velocity and the gas temperature, respectively. The flow pattern, as expected, has its maximum at the inlet and the trajectory strongly connects to the outlet with the region at the symmetry axis having low fluid velocities. The background gas temperature attains 1000 K with elastic collisions of electrons with the background gas being the main heating mechanism.

Figure 1: Period averaged electron density for 250 W of coil power and Vrf = 100 V.

Figure 2: Period averaged electron temperature for 250 W of coil power and Vrf = 100 V.

Figure 3: Period averaged electric potential for 250 W of coil power and Vrf = 100 V.

Figure 4: Current and voltage at the biased electrode for 250 W of coil power and Vrf = 100 V.

Figure 5: Magnitude of the fluid velocity for 250 W of coil power and Vrf = 100 V.

Figure 6: Gas temperature for 250 W of coil power and Vrf = 100 V.

References

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

2. SIGLO database, www.lxcat.net, retrieved 2022.

3. A. Friedman, Plasma Chemistry, Cambridge University Press, 2008

4. M.A. Lieberman and A.J. Lichtenberg, Principles of Plasma Discharges and Materials Processing, John Wiley & Sons, 2005.

Plasma_Module/Inductively_Coupled_Plasmas/icp_ccp_argon_chlorine

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

Create the geometry of a generic ICP reactor.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 30.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 30.

In the text field, type 3.

Locate the section. In the text field, type 20.

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 3.

In the text field, type 23.

Array 1 (arr1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

In the text field, type 4.

Locate the section. In the text field, type 4.5.

Rectangle 4 (r4)

In the toolbar, click

In the window for , locate the section.

In the text field, type 15.

In the text field, type 5.

Rectangle 5 (r5)

In the toolbar, click

In the window for , locate the section.

In the text field, type 6.

In the text field, type 5.

Locate the section. In the text field, type 15.

Rectangle 6 (r6)

In the toolbar, click

In the window for , locate the section.

In the text field, type 15.

Locate the section. In the text field, type 5.

Rectangle 7 (r7)

In the toolbar, click

In the window for , locate the section.

In the text field, type 10.

In the text field, type 6.

Locate the section. In the text field, type 20.

In the text field, type 14.

Line Segment 1 (ls1)

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. In the text field, type 20.

In the text field, type 17.

Locate the section. In the text field, type 20.

In the text field, type 18.

Mesh Control Edges 1 (mce1)

In the toolbar, click

and choose .

On the object , select Boundaries 6 and 29 only.

In the toolbar, click

Click the

button in the toolbar.

Geometry 1

In the window, collapse the > node.

Add parameters for use in the model, such as the mole fractions of Cl2 and Ar as well as the amplitude of the RF bias Vrf.

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
Isopar	0.1	0.1	Ions isotropic diffusion
Picp	250[W]	250 W	Coil power
Qf	50	50	Mass flow in SCCCM
xCl2	0.9	0.9	Inflow mole fraction of Cl2
xAr	1-xCl2	0.1	Inflow mole fraction of Ar
Vrf	0[V]	0 V	Voltage amplitude of rf bias

Create explicit selections to simplify the model creation.

Definitions (comp1)

Coil

In the toolbar, click

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

Select Domains 5–7 and 9 only.

Coil Boundaries

Right-click and choose .

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

Locate the section. From the list, choose .

Walls

In the toolbar, click

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

Select Domain 2 only.

Locate the section. From the list, choose .

Add materials to define the properties of the coil, dielectric window, and air surrounding the coil.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the tree, select > .

Click the button in the window toolbar.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Materials

Air (mat1)

In the window for , locate the section.

Click

Select Domain 4 only.

Copper (mat2)

In the window, click .

In the window for , locate the section.

From the list, choose .

Glass (quartz) (mat3)

In the window, click .

Select Domain 3 only.

Definitions (comp1)

In the window, collapse the > node.

Materials

In the window, collapse the > node.

Click the

button in the toolbar.

In the dialog, select > in the tree.

In the tree, select the checkbox for the node > .

Click .

Select to add a correction term to the heavy species flux.

This model needs stabilization because the density of the negative ions can drop sharply and attain very small values at the reactor edges. The reaction source stabilizations add artificial creation to prevent the densities to reach very small values, and isotropic diffusion prevents the appearance of instabilities when the negative ion density drops sharply.

Set the in the extra dimension to 1 because in the first study only the ICP power is added and resolution along the period is not necessary.

Plasma, Time Periodic (ptp)

In the window, under click .

In the window for , click to expand the section.

Select the checkbox.

Locate the section. Find the subsection. Select the checkbox.

Locate the section. Select the checkbox.

Locate the section. From the list, choose .

Locate the section. In the N text field, type 1.

From the list, choose .

Click to expand the section. Select the checkbox.

In the δid,i text field, type Isopar.

Locate the section. Click

Select Domain 2 only.

The Plasma Chemistry Import feature

The next steps have instructions to use the feature to import a file that automatically creates the argon-chlorine plasma chemistry.

The following is set or created automatically:

Species properties using

Reaction group features for argon and chlorine

Surface reaction group features

The documentation accompanying the feature contains more information about the file structure and what can be set automatically.

Plasma Chemistry 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_Cl2_plasma_chemistry.txt.

Click

In the following, initial mole fractions and initial number density for ions are specified. The mass fraction of Ar is found from a mass constraint and the initial density of Cl2+ is found by requiring electroneutrality.

Species: Cl2

In the window, expand the > > node, then click .

In the window for , locate the section.

In the x0 text field, type xCl2.

Species: Cl

In the window, click .

In the window for , locate the section.

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

Species: Cl-

In the window, click .

In the window for , locate the section.

In the n0 text field, type 1e10.

Species: Cl+

In the window, click .

In the window for , locate the section.

In the n0 text field, type 1E10[1/m^3].

Species: Cl2+

In the window, click .

In the window for , locate the section.

Select the checkbox.

Species: Ar

In the window, click .

In the window for , locate the section.

Select the checkbox.

Species: Ar+

In the window, click .

In the window for , locate the section.

In the n0 text field, type 1E10[1/m^3].

Set the field for species with internal energy so that the reaction enthalpy is correctly computed.

Click to expand the section. In the Δh text field, type 15.80.

Species: Ars

In the window, click .

In the window for , locate the section.

In the Δh text field, type 11.5.

Species: Cl2+

In the window, click .

In the window for , locate the section.

In the Δh text field, type 11.49.

Species: Cl+

In the window, click .

In the window for , locate the section.

In the Δh text field, type 14.25.

Group - Species

In the window, collapse the > > node.

The surface reactions used in the model were created automatically but it is still necessary to specify at which boundary they are going to exist.

Surface Reactions - Walls

In the window, click .

Select Boundaries 6, 27, 33–36, and 38–40 only.

Surface Reactions - Wafer

In the window, click .

Select Boundary 4 only.

In the following, select the fluid velocity, temperature and pressure to come from the and interfaces.

Set the electron mobility. All other electron transport properties are computed from the given electron mobility.

Plasma Model 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the μeNn text field, type 0.7e24[1/V/m/s]*(ptp.Te/1[V])^-0.6*exp(1.3[V]/ptp.Te).

Initial Values 1

In the window, click .

In the window for , locate the section.

In the ne,0 text field, type 1E15[1/m^3].

In the ε0 text field, type 2[V].

Wall 1

In the toolbar, click

and choose .

The node sets boundary conditions for the electron transport equations.

In the window for , locate the section.

From the list, choose .

Ground 1

In the toolbar, click

and choose .

Select Boundaries 6, 33–36, 39, and 40 only.

The node is used here to fix the mole fraction of Ar, Cl2, and Cl at the boundary that represents the inlet.

The mole fraction of xCl is fixed at a low value to prevent it to attain large values due to dissociation of Cl2.

Inflow 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click

In the table, enter the following settings:


Species names	Mole fraction (1)
Ar	xAr

Click

In the table, enter the following settings:


Species names	Mole fraction (1)
Cl2	xCl2

Click

In the table, enter the following settings:


Species names	Mole fraction (1)
Cl	1e-4

Select Boundary 35 only.

Outflow 1

In the toolbar, click

and choose .

Select Boundary 39 only.

Dielectric Contact 1

In the toolbar, click

and choose .

Select Boundaries 27 and 38 only.

Add a node to set a RF bias.

Metal Contact 1

In the toolbar, click

and choose .

Select Boundary 4 only.

In the window for , locate the section.

From the list, choose .

Locate the section. In the Va text field, type Vrf.

Locate the section. Select the checkbox.

Plasma, Time Periodic (ptp)

In the window, collapse the > node.

Select the domains where the magnetic fields are to be solved and add the coil that is responsible for exciting the plasma.

Magnetic Fields (mf)

In the window, under click .

Select Domains 2–7 and 9 only.

Coil 1

In the toolbar, click

and choose .

In the window for , locate the section.

Select the checkbox.

From the list, choose .

In the Pcoil text field, type Picp.

Locate the section. From the list, choose .

Magnetic Fields (mf)

In the window, collapse the > node.

Select the domain where the fluid flow is to be solved and select the fluid properties to come from the interface and the fluid temperature to come from the interface.

Also, set an inlet and outlet for the flow in the system, and set the operation pressure.

Laminar Flow (spf)

In the window, under click .

In the window for , locate the section.

Click

Select Domain 2 only.

Locate the section. In the pref text field, type 0.015[torr].

Inlet 1

In the toolbar, click

and choose .

Select Boundary 35 only.

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the Qsccm text field, type Qf.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 39 only.

Select the domain where the interface is to be solved and select the fluid properties to come from the interface and the fluid velocity and pressure to come from the interface.

Laminar Flow (spf)

In the window, collapse the > node.

Heat Transfer in Fluids (ht)

In the window, under click .

In the window for , locate the section.

Click