Dipolar Microwave Plasma Source

This model presents a 2D axisymmetric dipolar microwave plasma source sustained through resonant heating of the electrons. This is known as electron cyclotron resonance (ECR), which occurs when a suitable high magnetic flux density is present along with the microwaves.

This is an advanced model that showcases many of the features that make COMSOL Multiphysics unique, including:

•

Infinite elements for the magnetostatic model.

•

Functional-based mesh adaptation to create a fine mesh on the ECR surface.

•

PMLs for the electromagnetic waves to represent infinite space.

•

Degrees of freedom for all 3 components of the high-frequency electric field despite the fact that the problem is geometrically axisymmetric.

•

Full anisotropic tensors for the plasma conductivity and charged particle transport properties.

•

Resonant power absorption in the ECR surface by the electrons.

•

Equation-based modeling using integrated quantities to fix the total absorbed power.

•

Solver sequencing to first compute the static magnetic field, then solve for all the plasma components.

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

Model Definition — Static Magnetic Field

For the static magnetic field, Ampère’s law governs the azimuthal component of the magnetic vector potential:

where the external current density,

only has an azimuthal component and is defined in the coil as:

where N is the number of turns in the coil I is the total current and A is the cross sectional area. To represent the fact that the coil is in free space, infinite elements are used far away from the coil, as shown in Figure 1. A stationary study type is used to compute the static magnetic field. This field is then fed into a self-consistent model for the plasma.

Figure 1: Basic concept for the plasma source. A stationary azimuthal current flows in the coil which generates a static magnetic field in the surroundings. Resonant heating of the electrons occurs on the contour of the critical magnetic flux density.

The plasma conductivity becomes a full tensor in the presence of a static magnetic field. At some critical magnetic field the electrons continually gain energy from both the electric and magnetic fields over one RF cycle. This leads to a resonance zone in the plasma where the incoming electromagnetic wave is completely absorbed over a very short distance. The critical magnetic field is only dependent on the angular frequency, the electron mass, and the charge:

At 2.45 GHz the critical magnetic flux density is 875 gauss or 0.0875 T. Therefore you can use functional-based mesh adaptation to ensure that the ECR surface is adequately meshed for the plasma model. The functional is somewhat arbitrary; it is chosen such that it is zero everywhere but becomes large at the resonant magnetic flux density. In this model, use the functional

(1)

where δ is a small number to prevent division by zero.

Model Definition — Microwave Plasma

In this example, you solve the following wave equation for the high-frequency component of the electric field in the frequency domain:

Here σ is the plasma conductivity, which is a full tensor and a function of the electron density, collision frequency, and the static magnetic flux density. Using the definitions

(2)

where q is the electron charge, me is the electron mass, νe is the electron-neutral collision frequency, and ω is the angular frequency. The inverse of the plasma conductivity is defined as

(3)

where ne is the electron number density. Using the inverse of the plasma conductivity is convenient because it can be written in a compact form. COMSOL Multiphysics automatically computes the tensor form of the plasma conductivity for you by inverting Equation 3. Because the plasma conductivity tensor is a full tensor, all three components of the electric field are computed despite the fact that the only excitation from the coaxial port occurs in the rz-plane. The nonlinearity in the plasma conductivity can be seen in Figure 2. The surface represents four of the components of the nondimensional plasma conductivity versus the r- and z-components of the magnetic flux density (indicated by the x-axis and y-axis, respectively) on a log scale. At the resonant flux density (0.0875 T) the plasma conductivity is more than 106 higher than the case where no static magnetic field is present.

In Ref. 1 the size of the resonance is smoothed over a distance which can be resolved by the mesh. It is argued that this has a physical basis corresponding to collision-less heating. In Ref. 2 the physical reasoning behind the broadening of the resonance zone is Doppler shifting of the electrons into resonance. The same smoothing used in Ref. 1 is available in the COMSOL software by selecting the check box in the Microwave Plasma interface properties. In this case, the collision frequency, νe in Equation 2 is replaced by an effective collision frequency:

(4)

where δ is chosen to be 20. This is very simple from an implementation point of view but does lead to unphysical power absorption away from the resonance zone. The approach taken in Ref. 2 leads to the ECR surface being broadened only at the resonance zone.

Figure 2: Plots of the four components of the plasma conductivity tensor.

Compute the electron number density and electron energy density by solving a pair of drift-diffusion equations:

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

The electron transport properties are, like the plasma conductivity, full tensors. The electron mobility in the direction of the magnetic field lines is up to 8 orders of magnitude higher than the cross-field electron mobility. As such, electrons are only transported along magnetic field lines. The inverse of the electron mobility can be written in compact form as

(5)

where μdc is the electron mobility in the absence of a magnetic field. The COMSOL software automatically inverts the matrix in Equation 5 for you. The source coefficients in the above equations are determined by the plasma chemistry and are written using rate coefficients. Suppose that there are M reactions that contribute to the growth or decay of electron density and P inelastic electron-neutral collisions. In general, P >> M. In the case of rate coefficients, 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 multiphysics interface. The rate coefficients can 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. In this model the distribution function is chosen to be Maxwellian:

where

is the mean electron energy:

and

where Γ is the gamma function. The heating term, E · Γe has two components, one due to electron motion in the ambipolar field in the rz-plane and one due to heating of the electrons by the microwaves. Heating due to the microwaves is handled in the same way as described in Ref. 1. The power transferred from the electromagnetic field to the electrons is normalized so that 10 W of total power is absorbed by the electrons. This is accomplished by multiplying the heating term by a factor, α, which is defined as:

(6)

The result of this is that exactly 10 W of total power is transferred from the electromagnetic field to the electrons. If you did not apply this renormalization of the absorbed power then there would be nothing to stop the plasma from simply self-extinguishing or absorbing an inordinate amount of power. This approach is perfectly valid because the microwave equations are linear. The only drawback from this method is that the S-parameters given on the coaxial port are not valid because the fields are decoupled from the plasma model. Furthermore, it is not possible to self-consistently compute the ratio of the absorbed and reflected power through the excitation port, a quantity that may be of interest.

For nonelectron species, the following equation is solved for the mass fraction of species k:

(7)

As with the electrons, the ion transport properties are functions of the static magnetic flux density. The magnetic force is included as it can generate a significant ion velocity in the azimuthal direction close to the antenna. The ion mobility is also a function of the ambipolar electric field and is specified as a lookup table. The ion diffusion velocity, vk, is related to the diffusive flux via

where

(8)

Equation 7 can be re-arranged to give an expression for the diffusion velocity as a function of the other variables:

where

COMSOL automatically inverts Equation 8 when defining the diffusion velocity for each of the ionic species. 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.

Plasma Chemistry

The model considers argon plasma chemistry with the following set of collisions including elastic, excitation, direct ionization and stepwise ionization. Penning ionization and metastable quenching are also included in the model (electron 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	-

On surfaces, the following two reactions are considered:

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

Boundary Conditions

The above partial differential equations must be supplemented with a suitable set of boundary conditions. The coaxial port boundary condition is used to drive the electromagnetic waves. The port power is inconsequential due to the normalization scheme used on the absorbed power.

For the electrons, neglect reflection as well as secondary and thermal emission to get the following boundary condition on the electron flux:

and the electron energy flux:

Losses at the wall for the heavy species is due to surface reactions and migration due to the ambipolar field:

The reactor walls are grounded.

Results and Discussion — Static Magnetic Field Model

Figure 3 and Figure 4 present the results from the first study step. As expected, the azimuthal current in the coil generates a static magnetic field that has a “3”-shaped contour at a flux density of 0.0875 T. The magnetic field lines form a circular pattern around the coil, which is important to bear in mind when discuss the transport of the charged particles later.

Figure 3: Plot of the static magnetic flux density on a log scale (filled contour), magnetic field lines (thin lines), and the ECR surface at 0.0875 T (thick black line).

In Figure 4 the mesh, which has adapted based on the functional given in Equation 1 is shown. The mesh has clearly been significantly refined around the contour of the resonant magnetic flux density. This is required to accurately resolve the region where all the power deposition to the electrons occurs. Functional-based mesh adaptation is a feature that makes the finite element approach more attractive than the finite difference approach for ECR modeling.

Figure 4: Mesh generated after one refinement using functional-based mesh adaptation. The mesh is very fine on the ECR surface and relatively coarse away from the resonance zone.

Results and Discussion — Microwave Plasma Model

The electron density at the quasi steady state solution is plotted in Figure 5. The peak electron density is around 5·1016 m−3 and peaks radially outward from the center of the coil. The magnitude of the electron number density and its profile agree well with the results in Ref. 1.

Figure 5: Plot of the electron density at the quasi steady-state condition. The peak electron density is still below the critical electron density at the chosen operating frequency.

Despite the sharply peaked heating the electron temperature, plotted in Figure 6 does not show such peaks. Recall from Figure 3 that the magnetic field lines show the circular pattern away from the coil. This leads to strong energy transport along the field lines and very little transport across the magnetic field lines. Indeed the circular pattern along which the electron temperature is constant is consistent with the magnetic field lines. The peak electron temperature is around 3.8 eV and around 1.78 eV below the coil which is again, consistent with the results in Ref. 2.

Figure 6: Plot of the electron temperature, which peaks at around 3.8 eV.

Despite the fact that power is only deposited to the plasma on the ECR surface, the electron temperature, plotted in Figure 6, is not sharply peaked at the critical magnetic flux density. Recall from Figure 3 that the magnetic field lines show the circular pattern away from the coil. The high degree of anisotropy in the electron transport properties results in strong energy transport along the magnetic field lines and little transport across the magnetic field lines. Indeed, the circular pattern along which the electron temperature is constant is consistent with the magnetic field lines. The peak electron temperature is around 3.8 eV and around 1.78 eV below the coil, which is again consistent with the results in Ref. 2.

Figure 7: Plot of the plasma potential.

The electron density profile shows no signs of the resonance zone, which is clearly seen in Figure 8. The power deposition is very high, peaking at 35 W/cm3. All of the power deposition into the plasma from the electromagnetic field occurs in this resonance zone.

Figure 8: Plot of the power deposition into the plasma. Nearly all the power deposition occurs on the ECR surface.

The ionization source, plotted in Figure 9 is more highly localized around the coil. This corresponds to the region where the electron density and electron temperature are highest. Because the ionization rate scales linearly with the electron density and exponentially with the electron temperature this is to be expected.

Figure 9: Plot of the rate expression for electrons generated via ionization.

The plasma potential, plotted in Figure 7, peaks at around 16 V. The plasma potential is uniform throughout the plasma, even though the electron temperature shows large variations. The physical basis for the flat plasma potential is explained in Ref. 1.

Figure 10: Plot of the electron mobility tensor’s rr-component. The mobility varies by 8 orders of magnitude over the space of only a couple of centimeters.

The degree of anisotropy in the electron transport properties can be seen in Figure 10 and Figure 11. In Figure 10 the electron mobility varies by 8 orders of magnitude, it is 4·104 m2/(Vs) toward the coil edges and 10−4 m2/(Vs) radially outward from the coil center. In Figure 11 the opposite is true, the electron mobility is very high in the z direction at the center of the coil, and very small toward the coil edges. This leads to migration of electrons along the magnetic field lines when they are produced in the ionization region; see Figure 9.

Figure 11: Plot of the electron mobility tensor’s zz-component.

Figure 12: Unnormalized radial component of the microwave conduction current.

The conduction current due to the microwaves is plotted in Figure 12 – Figure 14. The largest component of the conduction current is actually in the azimuthal direction despite the coaxial port only propagating in the TM mode. Despite this, the heating (cooling) due to the dot product of the azimuthal components of the current and electric fields is small, due to the much lower value of the azimuthal component of the electric field.

Figure 13: Unnormalized axial component of the microwave conduction current.

Figure 14: Unnormalized azimuthal component of the microwave conduction current.

Finally, the trace of the plasma conductivity is plotted in Figure 15. The resonance zone is evident and the locally high electrical conductivity leads to the propagating electromagnetic waves to be absorbed.

It is worth mentioning that the electron density in this example model is below the critical plasma density everywhere (7.4·1016 m−3 at 2.45 GHz). If either the pressure or the power is increased, the power absorption can shift from the ECR surface to the contour where the plasma density is equal to the critical plasma density. On this contour the phase velocity approaches infinity whereas the group velocity approaches zero. The numerical instabilities caused by this are also smoothed out by adding an effective collision frequency to the actual collision frequency using Equation 4.

Figure 15: Plot of the trace of the plasma conductivity tensor.

Reference

1. G.J.M. Hagelaar, K. Makasheva, L. Garrigues, and J.-P. Boeuf, “Modelling of a dipolar microwave plasma sustained by electron cyclotron resonance”, J. Phys. D: Appl. Phys., vol. 42, p. 194019 (12pp), 2009.

2. R.L. Kinder and M.J. Kushner, “Consequences of mode structure on plasma properties in electron cyclotron resonance sources”, J. Vac. Sci. Technol. A, vol. 175, Sep/Oct 1999.

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

Notes About the COMSOL Implementation

This problem is solved in two stages. First, compute the static magnetic field using adaptive mesh refinement. Then, in a separate study step, solve for the electron density, electron energy density, mass fraction of argon ions, and mass fraction of electronically excited argon atoms, as well as the plasma potential and the 3 components of the high-frequency electric field. The magnetic flux density computed in the first study step is used to define the tensor plasma conductivity as well as electron and ion transport properties.

Plasma_Module/Wave-Heated_Discharges/dipolar_ecr_source

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

Click the

button in the toolbar.

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

Click .

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
r0	0.12	0.12	Plasma source radius
z0	0.24	0.24	Plasma source height
Psp	10[W]	10 W	Total power absorbed by the plasma setpoint

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type r0.

In the text field, type z0.

Locate the section. In the text field, type -z0/2.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.01.

In the text field, type 0.03.

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

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.004.

In the text field, type 0.048.

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

In the text field, type 0.072.

Rectangle 4 (r4)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.004.

In the text field, type 0.05.

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

Rectangle 5 (r5)

In the toolbar, click

In the window for , locate the section.

In the text field, type r0-0.01.

In the text field, type 0.02.

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

In the text field, type 0.12.

Rectangle 6 (r6)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.02.

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

In the text field, type 0.12.

Rectangle 7 (r7)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.02.

In the text field, type z0.

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

In the text field, type -z0/2.

Rectangle 8 (r8)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.02.

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

In the text field, type -0.02-z0/2.

Rectangle 9 (r9)

In the toolbar, click

In the window for , locate the section.

In the text field, type r0.

In the text field, type 0.02.

Locate the section. In the text field, type -0.02-z0/2.

Rectangle 10 (r10)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.01.

In the text field, type 0.02.

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

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

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

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

Click

Definitions

In the window, expand the node.

Right-click and choose .

Axis

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type -0.05.

In the text field, type 0.16.

In the text field, type -0.02.

In the text field, type 0.12.

Click

In the window, click .

In the window for , locate the section.

Select the check box.

Walls

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundaries 4, 6, 18–20, 22, and 26 only.

Right-click and choose .

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

Click .

Infinite Element Domain 1 (ie1)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Domains 1, 5, and 8–11 only.

Materials

Material 1 (mat1)

In the window, under right-click and choose .

Select Domains 1, 2, 4, 5, and 7–11 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 .

Select Domain 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	6e7	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

Material 3 (mat3)

Right-click and choose .

Select Domain 6 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	2	1	Basic