GEC CCP Reactor, Argon Chemistry

The NIST (National Institute of Standards and Technology) Gaseous Electronics Conference (GEC) Capacitively Coupled Plasma (CCP) reactor provides a standardized platform for studying capacitively coupled plasmas. Even the simplest plasma models are quite involved, so a 2D example helps in understanding the physics without excessive CPU time. The periodic steady state solution of an Argon discharge is computed and good agreement is obtained when comparing with measurements and simulations in the literature Ref. 1.

Model Definition

In this example, the GEC CCP reactor is simulated using the COMSOL Multiphysics Plasma Module. The simulations are for an Argon plasma sustained at a pressure of 100 mTorr by a periodic electric excitation of 13.56 MHz. The model is 2–dimensional and describes the space and time–periodic evolution of several macroscopic properties of the discharge. The reactor is electric asymmetric with a powered electrode of approximately 10 cm diameter and a gap between electrodes of 2.45 cm.

The mechanisms of power deposition into a CCP reactor is highly nonlinear and occurs at multiple different frequencies. Therefore, the electrostatic potential cannot be solved for in the frequency domain; the model must describe the periodic evolution of the charged particles to capture the nonlinear power absorption behavior.

electrical excitation

The driven electrode has a fixed power which computes the self DC bias. This corresponds to the following expression and set of constraints on the electric potential:

(1)

(2)

(3)

The constraint in Equation 2 is used to compute the self DC bias, Vdc,b. The constraint in Equation 3 is used to compute the RF potential, Va such that a fixed amount of power is deposited into the plasma.

plasma chemistry

Argon has 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 2 species and 3 reactions presented in Table 1. As in Ref. 1, the mass fraction of the metastable state is not computed. Electron impact cross-sections are obtained from Ref. 2.

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
4	e+Ar=>2e+Ar+	Ionization	15.8

In addition to volumetric reactions, the following surface reaction is implemented:

Table 2: table of surface reactions.
reaction	formula	Sticking Coefficient	Secondary emission coefficient	Mean energy of secondary electrons (V)
1	Ar+=>Ar	0	0.07	5.8

The ion use his internal energy to extract one electron from the wall with a probability of 0.07 and a mean energy of 5.8 V. The sticking coefficient is zero meaning that losses to the wall are assumed to be due to migration only.

Results and Discussion

The results presented in this section are for 1 W of power absorbed by the plasma. The voltage and current at the power electrode are presented in Figure 1. The amplitude of the applied voltage excitation and the DC self bias are numerically found to be approximately 100 V and -78 V, respectively. A DC self bias is present because the reactor is electrically asymmetric. This voltage ensures that the period averaged conduction current through the powered electrode is zero. Note also that the current collected at the electrode is not sinusoidal, meaning that there is power absorbed at harmonics higher than the fundamental.

Figure 2 to Figure 6 present different period–averaged plasma quantities. Figure 2 shows that the plasma is more dense within the small gap and attains the maximum density off–axis, as in the simulations of Ref. 1 and the measurements from Overzet et al (taken from Ref. 1). As expected, the period averaged potential presents an off–axis maximum related with the maximum of the plasma density, and a DC component at the powered electrode imposed by the calculated DC self–bias.

The ionization rate presented in Figure 6 also exhibits the off–axis behavior but the maximum does not coincide with the maximum electron density as reported in Ref. 1.

Due to the reactor electric asymmetry, the smaller electrode (the powered electrode in this simulation) has a larger sheath with more intense electric fields. The consequences are higher electron temperatures and power deposition near the powered electrode.


	The period averaged electron temperature in the COMSOL Multiphysics Plasma Module is computed from (4) where the brackets represent the period averaged. To obtain electron temperatures similar to the ones in Ref. 1the average should be computed as (5).

Figure 7 and Figure 8 shows profiles of the electron density along the axial and radial coordinate. These results are in good agreement with the data of figure 7 in Ref. 1 where simulations and measurements are compared.

Figure 1: Plot of the V-I characteristics of the discharge. Note the significant self DC bias due to the asymmetry of the discharge.

Figure 2: Plot of the period averaged electron density.

Figure 3: Plot of the period averaged electron temperature.

Figure 4: Plot of the period averaged electric potential.

Figure 5: Plot of the period averaged power deposition to the electrons.

Figure 6: Plot of the period averaged ionization rate (note the units are non-SI in this plot, for easier comparison with Ref. 1).

Figure 7: Plot of the on–axis electron density. Here, x=0 corresponds to the driven electrode.

Figure 8: Plot of the radial electron density profile half way across the discharge gap.

References

1. J.P.Boeuf and L.C. Pitchford, “Two-dimensional model of a capacitively coupled rf discharge and comparisons with experiments in the Gaseous Electronics Conference reference reactor,” Physical Review E, vol. 51, no .2, pp. 1376–1390, 1995.

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

Plasma_Module/Capacitively_Coupled_Plasmas/argon_gec_ccp

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, The interface will be used to compute the periodic steady state solution for the GEC reference cell, using the study.

click

In the tree, select .

Click .

Click

In the tree, select .

Click

Add some parameters for the geometric dimensions, the frequency and the input power.

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
L	2.54[cm]	0.0254 m	Discharge gap
R1	5.38[cm]	0.0538 m	Inner radius
R2	10.16[cm]	0.1016 m	Outer radius
Hd	10.16[cm]	0.1016 m	Chamber height
dThick	3[mm]	0.003 m	Dielectric thickness
f0	13.56[MHz]	1.356E7 Hz	Frequency
P0	1[W]	1 W	Power

Geometry 1

The geometry includes a point for the dielectric break, and three additional points to help create a mapped mesh later on. Since the interface generates a tremendous number of degrees of freedom, creating a geometry which can be efficiently meshed will help reduce computation time and memory requirements when solving.

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type R1.

In the text field, type L.

Click

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type R2-R1.

In the text field, type Hd.

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

In the text field, type -Hd/2+L/2.

Click

Click the

button in the toolbar.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type R1-dThick.

Click

Point 2 (pt2)

Right-click and choose .

In the window for , locate the section.

In the text field, type L.

Click

Point 3 (pt3)

In the window, under right-click and choose .

In the window for , locate the section.

In the text field, type R2.

Click

Point 4 (pt4)

In the window, under right-click and choose .

In the window for , locate the section.

In the text field, type R2.

Click

Line Segment 1 (ls1)

In the toolbar, click

and choose .

On the object , select Point 1 only.

It might be easier to select the correct by using the window. To open this window, in the toolbar click and choose . (If you are running the cross-platform desktop, you find in the main menu.)

In the window for , locate the section.

Find the subsection. Click to select the

toggle button.

On the object , select Point 1 only.

Add some mesh control edges to facilitate meshing without affecting the setup of the physics or postprocessing.

Mesh Control Edges 1 (mce1)

In the toolbar, click

and choose .

In the window for , locate the section.

Clear the check box.

On the object , select Boundaries 4 and 9 only.

In the toolbar, click

Now that the geometry is defined, add some selections to facilitate the setup of the physics.

Definitions (comp1)

Walls

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select the check box.

Select Boundaries 2–12 only.

In the text field, type Walls.

Driven Electrode

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundary 2 only.

In the text field, type Driven Electrode.

Dielectric Contact

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundaries 4 and 6 only.

In the text field, type Dielectric Contact.

Grounded Walls

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundaries 3, 5, and 7–12 only.

In the text field, type Grounded Walls.

For the physics, we will solve for the heavy species only in the base geometry, and use 30 elements in the extra dimension which represents one RF period.

Plasma, Time Periodic (ptp)

In the window, under click .

In the window for , locate the section.

In the Pxd text field, type 1/f0.

In the N text field, type 30.

From the list, choose .

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_reduced.txt.

Species: Ar

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. Select the check box.

Species: Ars

Disable the metastable species, since this is not included in the reference paper.

In the window, right-click and choose .

For the Argon ions, use the built in lookup table for the mobility, and the local field approximation for the ion temperature.

Species: Ar+

In the window, click .

In the window for , locate the section.

Select the check box.

Locate the section. From the list, choose .

Click to collapse the section. Click to expand the section. From the list, choose .

From the list, choose .

Click to expand 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.

Locate the section. In the γf text field, type 0.

Enter the temperature, pressure and electron mobility as defined in the reference.

Plasma Model 1

In the window, click .

In the window for , locate the section.

In the T text field, type 300[K].

In the pA text field, type 0.1[torr].

Locate the section. From the list, choose .

In the μe text field, type 3E5[cm^2/(V*s)]/0.1.

Wall 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the re text field, type 5/11.

It is much more stable to drive the electrode with a fixed power rather than a voltage. Activate computation of the DC self bias since that is included in the reference.

Metal Contact 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Prf text field, type P0.

In the fp text field, type f0.

Locate the section. Select the check box.

Dielectric Contact 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Ground 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

To minimize the number of degrees of freedom, take some time to create a mesh which is fine where the discharge is expected to show large gradients, and use a much coarser mesh elsewhere.

Mesh 1

Mapped 1

In the toolbar, click

In the window for , click to expand the section.

Clear the check box.

Distribution 1

Right-click and choose .

Select Boundaries 4 and 5 only.

In the window for , locate the section.

In the text field, type 4.

Distribution 2

In the window, right-click and choose .

Select Boundaries 2 and 3 only.

In the window for , locate the section.

From the list, choose .

In the text field, type 20.

In the text field, type 5.

Distribution 3

Right-click and choose .

Select Boundaries 7 and 9 only.

In the window for , locate the section.

From the list, choose .

In the text field, type 15.

In the text field, type 10.

Select the check box.

Distribution 4

Right-click and choose .

Select Boundaries 6, 8, 10, and 12 only.

In the window for , locate the section.

From the list, choose .

In the text field, type 15.

In the text field, type 8.

Select the check box.

Distribution 5

Right-click and choose .

Select Boundaries 1, 11, 13, and 14 only.

In the window for , locate the section.

From the list, choose .

In the text field, type 25.

In the text field, type 3.

Select the check box.

Click

The model will take about 20 minutes to solve, depending on hardware. Using , it is possible to visualize the evolution of the electron density at each nonlinear iteration. To do this, first get the initial value.

Time Periodic Study

In the window, click .