Plasma Application Library Manual

Plasma Module
Application Library Manual

Plasma Module Application Library Manual
19982016 COMSOL
Protected by U.S. Patents listed on www.comsol.com/patents, and U.S. Patents 7,519,518; 7,596,474;
7,623,991; 8,219,373; 8,457,932; 8,954,302; 9,098,106; 9,146,652; and 9,323,503. Patents pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement (www.comsol.com/comsol-license-agreement) and may be used or copied only under the terms
of the license agreement.
COMSOL, the COMSOL logo, COMSOL Multiphysics, Capture the Concept, COMSOL Desktop,
LiveLink, and COMSOL Server are either registered trademarks or trademarks of COMSOL AB. All other
trademarks are the property of their respective owners, and COMSOL AB and its subsidiaries and products
are not affiliated with, endorsed by, sponsored by, or supported by those trademark owners. For a list of such
trademark owners, see www.comsol.com/trademarks.
Version: COMSOL 5.2a
Contact Information
Visit the Contact COMSOL page at www.comsol.com/contact to submit general
inquiries, contact Technical Support, or search for an address and phone number. You can
also visit the Worldwide Sales Offices page at www.comsol.com/contact/offices for
address and contact information.
If you need to contact Support, an online request form is located at the COMSOL Access
page at www.comsol.com/support/case. Other useful links include:
Support Center: www.comsol.com/support

Product Download: www.comsol.com/product-download
Product Updates: www.comsol.com/support/updates
COMSOL Blog: www.comsol.com/blogs
Discussion Forum: www.comsol.com/community
Events: www.comsol.com/events
COMSOL Video Gallery: www.comsol.com/video
Support Knowledge Base: www.comsol.com/support/knowledgebase
Part number: CM021404

Created in COMSOL Multiphysics 5.2a
Dielectric Barrier Discharge
This model is licensed under the COMSOL Software License Agreement 5.2a.
All trademarks are the property of their respective owners. See www.comsol.com/trademarks.
Introduction
This model simulates electrical breakdown in an atmospheric pressure gas. Because
electrical breakdown is a complicated process, a 1D model is considered. To highlight the
physics of the breakdown process, this example uses a simple argon chemistry which keeps
the number of species and reactions to a minimum.
Model Definition
The operating principle for a dielectric barrier discharge is as follows: there is a small gap
filled with a gas in between two dielectric plates. The gap between the two dielectric plates
is typically less than one millimeter. On one of the dielectric plates, a sinusoidal voltage is
applied. The other plate is electrically grounded. As the voltage applied to the top plate
increases, a stronger electric field forms in the gap between the plates. Any free electrons
in the gap1 are accelerated and if the electric field is strong enough they may acquire
enough energy to cause ionization. This can lead to a cascade effect where the number of
electrons in the gap increases exponentially on a nanosecond timescale. Electrons created
via electron impact ionization rush towards one of the dielectric plates, in the opposite
direction to the electric field. An equal number of ions are also generated during electron
impact ionization (electrons and ions must be created in equal pairs to preserve the overall
charge balance). The ions rush toward the opposite dielectric plate in the same direction
as the electric field. As a result, surface charge with opposite sign accumulates on both
dielectric plates. This causes the electric field to become shielded from the gas filled gap.
In fact, the electric field across the gap cannot exceed the breakdown electric field, which
is highly dependent on the gas. The breakdown electric field is also a function of the
surface properties of the dielectric material. Surface charge accumulation temporarily
terminates the discharge until the field reverses direction and the process repeats in the
opposing direction.
Modeling dielectric barrier discharges in more than one dimension is, of course possible,
but the results can be difficult to interpret due to the amount of competing physics in the
problem. In this simple model the problem is reduced to 1D by assuming the dielectric
gap is much smaller than the diameter of the plates. It also allows us to quickly gain some
insight into the characteristics of the discharge without excessive computation time.
1. There are typically around 1,000,000 m

-3 free electrons in air at sea level.
2 | DIELECTRIC BARRIER DISCHARGE

The geometry for a typical dielectric barrier discharge is shown in Figure 1. The dielectric
plates may be up to 15cm in diameter, the dielectric and gap thickness are typically less
than 1 millimeter.
V = V0sin(t)
Dielectric
Dielectric
Ground
Figure 1: Graphic illustration of a typical dielectric barrier discharge.
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 Theory for the Drift Diffusion Interface in the Plasma Module Users Guide.
( n ) + [ n ( E ) D n ] = R
t e e e e e e
n
( ) + [ n ( E ) D n ] + E e = R
t
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 computed from the
electron mobility using:
D e = e T e, = --- e, D = T e
5
3
The source coefficients in the above equations are determined by the plasma chemistry
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:

M
Re = x j kj Nn ne
j=1
where xj is the mole fraction of the target species for reaction j, kj is the rate coefficient for
reaction j (SI unit: m3/s), and Nn is the total neutral number density (SI unit: 1/m3). The
electron energy loss is obtained by summing the collisional energy loss over all reactions:
R = xj kj Nn ne j
j=1
where j is the energy loss from reaction j (SI unit: V). The rate coefficients can be
computed from cross section data by the following integral:

kk = 0 k ( )f ( ) d
where = (2q/me)1/2 (SI unit: C1/2/kg1/2), me is the electron mass (SI unit: kg), is
energy (SI unit: V), k is the collision cross section (SI unit: m2), and f is the electron
energy distribution function. In this case 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 Users Guide.
( w k ) + ( u )w k = j k + R k
t
The electrostatic field is computed using the following equation:
0 r V =
The space charge density is automatically computed based on the plasma chemistry
specified in the model using the formula:
N
= q Z k n k n e

k = 1
For detailed information about electrostatics see Theory for the Electrostatics Interface in
the Plasma Module Users Guide.

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:
n e = --- e, th n e
1
2 p ( p n ) (1)
p
and the electron energy flux:
n = --- e, th n
5
6 p p ( p n ) (2)
p
The second term on the right-hand side of Equation 1 is the gain of electrons due to
secondary emission effects, p being the secondary emission coefficient. The second term
in Equation 2 is the secondary emission energy flux, p being the mean energy of the
secondary electrons. 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:
n j k = M w R k + M w c k Z k ( E n ) [ Z k k ( E n ) > 0 ]
Surface charge accumulation is added to the dielectric surfaces that are adjacent to the gap
where the plasma forms by way of the following boundary condition:
n ( D1 D 2 ) = s
where s is the surface charge density, which is computed by solving the following
distributed ODE on the surfaces:
d s
= n Ji + n Je
dt
where n Ji is the normal component of the total ion current density at the wall, and n
Je is the normal component of the total electron current density at the wall.
The discharge is driven by a sinusoidal electric potential applied to the exterior boundary
of one of the dielectric plates:
V = V 0 sin ( t )
where the applied peak voltage, V0 is 750 V and the angular frequency, the RF frequency
being 50 kHz. The exterior boundary of the other dielectric plate is grounded.

PLASMA CHEMISTRY
Argon is an attractive gas to use in a benchmark problem since only a handful of reactions
and a few species need to be considered. The list of chemical reactions considered is as
follows:
TABLE 1: TABLE OF COLLISIONS AND REACTIONS MODELED
REACTION FORMULA TYPE (eV)

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 -
Initially a small number of seed electrons are present. These are necessary in order to
initiate the discharge on the first RF cycle. In addition to the 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 the ions reach the wall they are assumed to change back to neutral argon atoms and
donate their charge to the wall.
Results and Discussion

It is more convenient to analyze the results of this one dimensional problem by extruding
the solution into two dimensions. The extra dimension represents time. In COMSOL
Multiphysics this is accomplished by adding a Parametric Extrusion 1D data set. The
surface plot is convenient because you can immediately see how the variables of interest
evolve over time.
The mass fraction of excited argon atoms is plotted in Figure 2. The excited species have
a much longer lifetime in the gap than the electrons or ions. This is because the primary
mechanism for destruction of excited argon species is de-excitation upon contact with the
wall. The excited argon atoms can only reach the wall via diffusion whereas the electrons
and ions reach the wall very rapidly due to migration. It is also apparent from Figure 2 that

the discharge reaches a periodic steady state solution after only two RF cycles.
Start up
cycle
Figure 2: Mass fraction of excited argon.
The electric potential is plotted in Figure 3. The voltage is relatively uniform across the
discharge gap. This can be seen more clearly by examining the electric field in Figure 4.
There is a much stronger electric field in the dielectric materials than in the gap. This is
because the surface charge which accumulates on the dielectric surfaces tends to shield out
the electric field.

Figure 3: Electric potential (x-axis) vs. time (y-axis).
Figure 4: Electric field across the gap (x-axis) vs. time (y-axis).

start-up
cycle
Figure 5: Extruded plot of the electron density.
Figure 6: Extruded plot of the mean electron energy.

Implicit in the equations solved for the number of charged particles and electrostatic
potential is that the total electrical current is conserved. Mathematically, this means that:

J = 0
t
where J is the total plasma current density (SI unit: A/m2), and is the space charge
density (SI unit: C/m3). Since electrons and ions are created in equal pairs the time
derivative of the space charge density should be approximately zero. For this 1D model,
this means the total current density must be constant across the gap at any snapshot in
time. In Figure 7, the current density due to electrons (left) and ions (right) is plotted.
The current density is not symmetric because of the different secondary emission
coefficients used on the dielectric surfaces.
Figure 7: Plot of the electron current density (left) and the ion current density (right) in the
discharge, excluding the first RF cycle.
The total plasma current density is plotted in Figure 8. As expected, the total current
density is constant across the gap at any point in time.

Figure 8: Plot of the total plasma current density (sum of the electron and ion current density),
excluding the first RF cycle. Conservation of charge requires that the total current density be
constant across the gap at any point in time.
The total current at the grounded electrode is plotted in Figure 9. In the absence of the
plasma, the current would be a perfect cosine wave. However, the presence of the plasma
and flow of charged particles leads to a non-sinusoidal current waveform. The
instantaneously absorbed power in the plasma is plotted in Figure 10. Time averaging this
over 1 RF cycle yields the power absorbed by the plasma. The power is around 16.7 W on
one half cycle and 17.7 W on the other half cycle. The difference is because the secondary
emission coefficients are different on the upper and lower plates.

Start-.up cycle
Figure 9: Plot of the total discharge current vs. time.
Figure 10: Plot of power vs. time for the dielectric barrier discharge.

Application Library path: Plasma_Module/Direct_Current_Discharges/
argon_dbd_1d
Modeling Instructions
From the File menu, choose New.
NEW
In the New window, click Model Wizard.
MODEL WIZARD
1 In the Model Wizard window, click 1D.
2 In the Select Physics tree, select Plasma>Plasma (plas).
3 Click Add.
4 Click Study.
5 In the Select Study tree, select Preset Studies>Time Dependent.
6 Click Done.
GEOMETRY 1
You start by defining the geometry for the problem. This model has a simple 1D geometry
consisting of 3 domains. Two dielectric domains and a gap where the plasma forms.
Interval 1 (i1)
1 On the Geometry toolbar, click Interval.
2 In the Settings window for Interval, locate the Interval section.
3 From the Number of intervals list, choose Many.
4 In the Points text field, type -1e-4, 0, 2e-4, 3e-4.
5 Click Build All Objects.
DEFINITIONS
Variables 1
1 On the Home toolbar, click Variables and choose Local Variables.
2 In the Settings window for Variables, locate the Variables section.

3 In the table, enter the following settings:
Name Expression Unit Description

freq 50e3[Hz] Hz RF Frequency
omega 2*pi*freq Hz Angular frequency
mueN 1E25[1/(m*V*s)] 1/(Vms) Reduced electron
mobility
Vrf -750[V]*sin(omega*t) V Applied voltage
dplate 0.1[m] m Plate diameter
As 0.25*pi*dplate^2 m Plate area
Now add the charge conservation to the dielectric materials. This means that only a
charge conservation equation is solved in the dielectric material and all the plasma
components are solved for in the gap between the dielectrics.
PLASMA (PLAS)
Charge Conservation 1
1 On the Physics toolbar, click Domains and choose Charge Conservation.
2 Select Domains 1 and 3 only.
Load in the argon cross sections from file. They form the basis of the plasma chemistry
under investigation.
Cross Section Import 1

1 On the Physics toolbar, click Global and choose Cross Section Import.
2 In the Settings window for Cross Section Import, locate the Cross Section Import section.
3 Click Browse.
4 Browse to the applications Application Libraries folder and double-click the file
Ar_xsecs.txt.
5 In the Model Builder window, click Plasma (plas).

6 In the Settings window for Plasma, locate the Out-of-Plane Thickness section.
7 In the A text field, type As.
8 Locate the Plasma Properties section. Select the Use reduced electron transport properties
check box.
Reaction 1
1 On the Physics toolbar, click Domains and choose Reaction.
2 In the Settings window for Reaction, locate the Reaction Formula section.

3 In the Formula text field, type Ars+Ars=>e+Ar+Ar+.
4 Locate the Reaction Parameters section. In the kf text field, type 3.3734e8.
Reaction 2
3 In the Formula text field, type Ars+Ar=>Ar+Ar.
4 Locate the Reaction Parameters section. In the kf text field, type 1807.
Species: Ar
1 In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click Species: Ar.
2 In the Settings window for Species, locate the Species Formula section.
3 Select the From mass constraint check box.
4 Locate the General Parameters section. From the Preset species data list, choose Ar.
Species: Ars
1 In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click Species:
Ars.
2 In the Settings window for Species, locate the General Parameters section.
3 In the x0 text field, type 1e-11.
4 From the Preset species data list, choose Ar.
Now let the initial number density of Argon ions be the same as the initial number of
electrons. This forces the plasma to be initially charge neutral.
Species: Ar+
Ar+.
3 Select the Initial value from electroneutrality constraint check box.
Plasma Model 1
1 In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click Plasma
Model 1.
2 In the Settings window for Plasma Model, locate the Electron Density and Energy section.
3 In the eNn text field, type mueN.

4 Locate the Model Inputs section. In the T text field, type 400[K].
The initial number density of seed electrons is very small, only one million free electrons
per cubic meter. This corresponds to a near zero conductivity. So, the gap is truly acting
as an insulator initially.
Initial Values 1
1 In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click Initial
Values 1.
2 In the Settings window for Initial Values, locate the Initial Values section.
3 In the ne, 0 text field, type 1e6.
4 In the 0 text field, type 5.
Secondary emission of electrons is important when studying discharge curves from
DBDs. In this example you add a higher secondary emission coefficient on the left wall.
Surface Reaction 1
1 On the Physics toolbar, click Boundaries and choose Surface Reaction.
2 Select Boundary 2 only.
3 In the Settings window for Surface Reaction, locate the Reaction Formula section.
4 In the Formula text field, type Ar+=>Ar.
5 Locate the Secondary Emission Parameters section. In the i text field, type 0.01.
6 In the i text field, type 2.5.
Surface Reaction 2
5 Locate the Secondary Emission Parameters section. In the i text field, type 1E-6.
6 In the i text field, type 2.5.
Surface Reaction 3
2 Select Boundaries 2 and 3 only.
4 In the Formula text field, type Ars=>Ar.

Wall 1
1 On the Physics toolbar, click Boundaries and choose Wall.
Surface charge will begin to accumulate when the gas begins to break down. This will
cause the electric field to be shielded in the gap. This is the phenomena responsible for
terminating the discharge and also the reason why the breakdown voltage cannot be
exceeded across the gap. COMSOL automatically computes the amount of surface
charge accumulation when the feature is added to the model. The surface charge
accumulation is computed by integration the electron and ion fluxes to the wall.
Surface Charge Accumulation 1

1 On the Physics toolbar, click Boundaries and choose Surface Charge Accumulation.
Ground 1
1 On the Physics toolbar, click Boundaries and choose Ground.
Terminal 1
1 On the Physics toolbar, click Boundaries and choose Terminal.
3 In the Settings window for Terminal, locate the Terminal section.
4 In the Terminal name text field, type electrode.
5 In the V0 text field, type Vrf.
Now assign the relative permittivity to the dielectric material and the air gap where the
plasma forms.
MATERIALS
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose
Blank Material.
Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1)>Materials right-click Material
1 (mat1) and choose Rename.
2 In the Rename Material dialog box, type Dielectric 1 in the New label text field.
3 Click OK.

5 In the Settings window for Material, locate the Material Contents section.
Property Name Value Unit Property

group
Relative permittivity epsilonr 10 1 Basic
MESH 1
There must be sufficient mesh density to resolve the sharp gradients in the electron and
ion density in the gap. Therefore you specify that there are 200 elements across the width
of the gap.
1 In the Settings window for Mesh, locate the Mesh Settings section.
2 From the Element size list, choose Extremely fine.
3 Right-click Component 1 (comp1)>Mesh 1 and choose Edge.
Edge 1
In the Model Builder window, under Component 1 (comp1)>Mesh 1 right-click Edge 1 and
choose Distribution.
Distribution 1
1 In the Settings window for Distribution, locate the Domain Selection section.
2 Click Clear Selection.
3 Select Domain 2 only.
4 Locate the Distribution section. From the Distribution properties list, choose Predefined
distribution type.
5 In the Number of elements text field, type 200.
6 In the Element ratio text field, type 5.
7 From the Distribution method list, choose Geometric sequence.
8 Select the Symmetric distribution check box.
9 Click Build All.
STUDY 1
Step 1: Time Dependent

1 In the Model Builder window, expand the Study 1 node, then click Step 1: Time
Dependent.
2 In the Settings window for Time Dependent, locate the Study Settings section.

3 Click Range.
4 In the Range dialog box, choose Number of values from the Entry method list.
5 In the Stop text field, type 1e-4.
6 In the Number of values text field, type 201.
7 Click Replace.
8 On the Home toolbar, click Compute.
RESULTS
Parametric Extrusion 1D 1
1 On the Results toolbar, click More Data Sets and choose Parametric Extrusion 1D.
2 In the Settings window for Parametric Extrusion 1D, locate the Settings section.
3 Clear the Separate levels check box.
Now create a new parametric dataset which ignores the first startup RF cycle so the
current density can be visualized later
1 Right-click Parametric Extrusion 1D 1 and choose Duplicate.
2 In the Settings window for Parametric Extrusion 1D, locate the Data section.
3 From the Time selection list, choose From list. In the Times (s) list, choose 1E-5 through
1E-4.
2D Plot Group 4
1 On the Results toolbar, click 2D Plot Group.
2 In the Settings window for 2D Plot Group, type Excited Argon Mass Fraction in
the Label text field.
Excited Argon Mass Fraction

In the Model Builder window, under Results right-click Excited Argon Mass Fraction and
choose Surface.
Surface 1
1 In the Settings window for Surface, click Replace Expression in the upper-right corner of
the Expression section. From the menu, choose Component 1>Plasma (Heavy Species
Transport)>Mass fractions>plas.wArs - Mass fraction.
2 On the Excited Argon Mass Fraction toolbar, click Plot.
3 Click the Zoom Extents button on the Graphics toolbar.

Excited Argon Mass Fraction
Right-click Excited Argon Mass Fraction and choose Duplicate.
Excited Argon Mass Fraction 1

In the Settings window for 2D Plot Group, type Electric Potential in the Label text
field.
Surface 1
1 In the Model Builder window, expand the Results>Electric Potential node, then click
Surface 1.
the Expression section. From the menu, choose Component 1>Plasma (Electrostatics)>
Electric>V - Electric potential.
3 On the Electric Potential toolbar, click Plot.
Electric Potential
In the Model Builder window, under Results right-click Electric Potential and choose
Duplicate.
Electric Potential 1
In the Settings window for 2D Plot Group, type Electric Field in the Label text field.
Surface 1
1 In the Model Builder window, expand the Results>Electric Field node, then click Surface 1.
the Expression section. From the menu, choose Component 1>Plasma (Electrostatics)>
Electric>Electric field>plas.Ex - Electric field, x component.
3 On the Electric Field toolbar, click Plot.
Electric Field
In the Model Builder window, under Results right-click Electric Field and choose Duplicate.
Electric Field 1
In the Settings window for 2D Plot Group, type Electron Density in the Label text field.
Surface 1
1 In the Model Builder window, expand the Results>Electron Density node, then click
Surface 1.
the Expression section. From the menu, choose Component 1>Plasma (Drift Diffusion)>
Electron density>plas.ne - Electron density.

3 On the Electron Density toolbar, click Plot.
Electron Density
In the Model Builder window, under Results right-click Electron Density and choose
Duplicate.
Electron Density 1
In the Settings window for 2D Plot Group, type Mean Electron Energy in the Label text
field.
Surface 1
1 In the Model Builder window, expand the Results>Mean Electron Energy node, then click
Surface 1.
Electron energy density>plas.ebar - Mean electron energy.
3 On the Mean Electron Energy toolbar, click Plot.
Mean Electron Energy

In the Model Builder window, under Results right-click Mean Electron Energy and choose
Duplicate.
Mean Electron Energy 1

1 In the Settings window for 2D Plot Group, type Electron Current Density in the
Label text field.
2 Locate the Data section. From the Data set list, choose Parametric Extrusion 1D 2.
Surface 1
1 In the Model Builder window, expand the Results>Electron Current Density node, then
click Surface 1.
Current>Electron current density>plas.Jelx - Electron current density, x component.
3 On the Electron Current Density toolbar, click Plot.
Electron Current Density

In the Model Builder window, under Results right-click Electron Current Density and choose
Duplicate.

Electron Current Density 1
In the Settings window for 2D Plot Group, type Argon Ion Current Density in the
Label text field.
Surface 1
1 In the Model Builder window, expand the Results>Argon Ion Current Density node, then
click Surface 1.
Transport)>Species>Species wAr_1p>Ion current density>plas.Jix_wAr_1p - Ion current
density, x component.
3 On the Argon Ion Current Density toolbar, click Plot.
Argon Ion Current Density

In the Model Builder window, under Results right-click Argon Ion Current Density and choose
Duplicate.
Argon Ion Current Density 1

In the Settings window for 2D Plot Group, type Total Conduction Current Density
in the Label text field.
Surface 1
1 In the Model Builder window, expand the Results>Total Conduction Current Density node,
then click Surface 1.
2 In the Settings window for Surface, locate the Expression section.
3 In the Expression text field, type plas.Jix_wAr_1p+plas.Jelx.
4 On the Total Conduction Current Density toolbar, click Plot.
1D Plot Group 12
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Terminal Current in the Label text
field.
Global 1
1 On the Terminal Current toolbar, click Global.
2 In the Settings window for Global, click Replace Expression in the upper-right corner of
the y-axis data section. From the menu, choose Component 1>Plasma (Electrostatics)>
plas.I_electrode - Current, Terminal electrode.
3 On the Terminal Current toolbar, click Plot.

1D Plot Group 13
2 In the Settings window for 1D Plot Group, type Total Power Deposition in the Label
text field.
Global 1
1 On the Total Power Deposition toolbar, click Global.
the y-axis data section. From the menu, choose Component 1>Plasma>Power and
collisions>plas.Pcap_tot - Total capacitive power deposition.
3 On the Total Power Deposition toolbar, click Plot.
Global Evaluation 1
1 On the Results toolbar, click Global Evaluation.
Use the timeavg operator to compute the time averaged power deposition for cycles
2-10.
2 In the Settings window for Global Evaluation, locate the Data section.
3 From the Time selection list, choose Last.
4 Locate the Expressions section. In the table, enter the following settings:
Expression Unit Description

timeavg(1e-5,2e-5,plas.Pcap_tot) W Cycle 2
5 Click Evaluate.
As you can see, the average power deposited to the plasma remains the same for each
cycle after only 3 RF cycles. The power is around 16.7 W on one half cycle and 17.7 W
on the other half cycle. The difference is due to the fact that the secondary emission
coefficients are different on the upper and lower plates.

GEC CCP Reactor, Argon Chemistry, 1D
Introduction
The NIST (National Institute of Standards and Technology) GEC CCP reactor provides
a platform for studying capacitively coupled plasmas (CCP). Even the simplest plasma
models are quite involved, so a 1D example helps in understanding the physics without
excessive CPU time. The problem has no steady-state solution, although a periodic
steady-state solution is reached after a suitable number of RF cycles (usually >1000).
Model Definition
The operating principle of a capacitively coupled plasma is different to the inductive case.
In a CCP reactor, there are no induction currents. Instead, the plasma is sustained by
applying a sinusoidal electrostatic potential across a small gap (2.54 cm in this case) filled
with a low pressure (1 Torr in this case) gas. This case uses argon because the plasma
chemistry is very simple. The mechanism 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 be transient, with
the system allowed to reach a periodic steady state solution. This condition is reached after
a suitable number of RF cycles.
n
( ) + [ n e ( e E ) D e n e ] = R e
t e
n
( ) + [ n ( E ) D n ] + E e = R
t
D e = e T e, = --- e, D = T e
5
3
2 | G E C C C P R E A C T O R, A R G O N C H E M I S T R Y, 1 D
Re = x j kj Nn ne
j=1
R = xj kj Nn ne j
j=1

kk = 0 k ( )f ( ) d

( w ) + ( u )w k = j k + R k
t k
0 r V =
N
= q Z k n k n e
k = 1
BOUNDARY CONDITIONS
n e = --- e, th n e
1
2 p ( p n ) (1)
p
n = --- e, th n
5
6 p p ( p n ) (2)
p
and the fact that the electric field is directed towards the wall:
Surface charge accumulation is added to the dielectric surface by way of the following
boundary condition:
n ( D1 D 2 ) = s
d s
= n Ji + n Je
dt
The discharge is driven by a sinusoidal electric potential applied to the exterior boundary
of the dielectric plate:
V = V 0 sin ( t )
where the applied peak voltage, V0 is 300 V and the angular frequency, the RF frequency
being 13.56 MHz. The exterior boundary of the other dielectric plate is grounded. The
left boundary is grounded.
PLASMA CHEMISTRY
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:

In a CCP reactor the electron density and density of excited species is relatively low so
stepwise ionization is not as important as in high density discharges. In addition to
volumetric reactions, the following surface reactions are implemented:
1 Ars=>Ar 1
2 Ar+=>Ar 1
When a metastable argon atom makes contact with the wall, it will revert to the ground
state argon atom with some probability (the sticking coefficient).

Due to the fact that capacitively coupled discharges never reach steady state, you have to
run the problem long enough to obtain a quasi steady state, that is, until the time-averaged
properties do not change from cycle to cycle. Once reached, quantities of interest can be
plotted at different phases of the RF cycle.
The electric potential, electric field, electron density, and ion density are all quantities of
interest. Figure 1 plots the electron density at different phases of the RF cycle.
Figure 1: Plot of the electron density at different phases of the RF cycle.
The electrons, with their low mass can respond to the electric field very rapidly. This means
that the electrons can move back and forth in the x-direction as the electric field changes
sign. The electrostatic potential is the same at t=0 and t=2. At these times the electrostatic
potential is positive in the bulk plasma. The electrons are more mobile than the heavier
ions, so they respond to the applied electric field over the course of one RF cycle. The
pushing and pulling of electrons back and forth, coupled with the fact that the ions
essentially remain stationary leads to a separation of space charge in the plasma sheath.
Figure 2: Plot of electron temperature at different phases of the RF cycle.
The space charge separation produces an electric field as well as a flux of electrons. This is
how power is deposited into the plasma. Figure 2 plots the electron temperature in the
plasma and Figure 3 and Figure 4 plot the electric potential and electric field respectively.
Figure 5 plots the ion density at different phases of the RF cycle. Over 1RF cycle, there is
a net positive power deposition. At quasi steady state, the power deposition is balanced by
the power lost due to chemical reactions and energy loss to the wall.
Figure 3: Plot of the electric potential at different phases of the RF cycle.
Figure 4: Plot of the electric field at different phases of the RF cycle.
Figure 6 plots the power deposition at different phases of the RF cycle. The maximum
power deposition occurs when the applied voltage is 0. Power deposition is at a minimum
when the applied voltage is at a maximum (300V) or minimum (-300V).
Figure 5: Plot of the ion density at different phases of the RF cycle.
This model also demonstrates two of the shielding mechanisms that occur in a capacitively
coupled plasma. The electric field is shielded from the plasma bulk over several Debye
lengths. The Debye length is defined as:
0 Te
D = -----------
-
ne q
where Te is in electron volts, 0 is the permittivity of free space, ne is the electron density
and q is the electron charge. This is the distance over which the electric field is shielded
from the plasma. One can see from Figure 4 that the electric field in the plasma bulk does
not change over 1 RF cycle.
The other shielding mechanism is caused by the difference in transport time scales between
electrons and ions. Since the electrons are more mobile than the ions, they can accumulate
a substantial drift velocity directed toward the walls, even at the high operating
frequencies. This results in an accumulation of surface charge on the dielectric surface.
This accumulation of surface charge tends to further shield the electric field from the
plasma.
Figure 6: Plot of the power deposition into the plasma at different phases of the RF cycle.
Reference
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, 722733, 2005.
Application Library path: Plasma_Module/Capacitively_Coupled_Plasmas/

argon_gec_ccp_1d
NEW
MODEL WIZARD
3 Click Add.
4 Click Study.
6 Click Done.
GLOBAL DEFINITIONS
Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
Name Expression Value Description

V0 300[V] 300 V Applied voltage
freq 13.56[MHz] 1.356E7 Hz RF Frequency
omega 2*pi*freq 8.52E7 Hz Angular frequency
p0 1[torr] 133.3 Pa Gas pressure
T0 400[K] 400 K Gas temperature
mueN 4E24[1/(V*m*s)] 4E24 1/(Vms) Reduced electron
mobility
GEOMETRY 1
Interval 1 (i1)
3 In the Right endpoint text field, type 0.0254.
Interval 2 (i2)
3 In the Left endpoint text field, type 0.0254.
PLASMA (PLAS)

3 Click Browse.
Ar_xsecs.txt.

6 In the Settings window for Plasma, locate the Plasma Properties section.
7 Select the Use reduced electron transport properties check box.
Reaction 1
4 Locate the Reaction Parameters section. In the kf text field, type 3.734E8.
Reaction 2
Surface Reaction 1
Surface Reaction 2
Species: Ar
Species: Ar+
Ar+.
Plasma Model 1
Model 1.
2 In the Settings window for Plasma Model, locate the Model Inputs section.
3 In the T text field, type T0.
4 In the pA text field, type p0.
5 Locate the Electron Density and Energy section. In the eNn text field, type mueN.
Wall 1
The surface charge accumulation models the effect of charged particle build up on
dielectric surfaces. This shields the electric field from the plasma.

Ground 1
The terminal feature is used to drive the CCP with a sinusoidal voltage. In this case the
voltage is 300 volts and the pressure is 1 torr. If the ratio of the voltage to the pressure
is higher, numerical instabilities may occur. In these cases the mesh refinement may be
necessary.
Terminal 1
3 In the Terminal name text field, type plate.
5 From the Terminal type list, choose Voltage.
6 In the V0 text field, type V0*sin(omega*t).
MATERIALS
Blank Material.
Material 1 (mat1)

group
Relative permittivity epsilonr 4.7 1 Basic
4 Right-click Component 1 (comp1)>Materials>Material 1 (mat1) and choose Rename.

5 In the Rename Material dialog box, type Dielectric in the New label text field.
6 Click OK.
MESH 1
In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
Edge.
Edge 1
1 In the Settings window for Edge, locate the Domain Selection section.
2 From the Geometric entity level list, choose Domain.
Distribution 1
1 Right-click Component 1 (comp1)>Mesh 1>Edge 1 and choose Distribution.
2 In the Settings window for Distribution, locate the Distribution section.
3 From the Distribution properties list, choose Predefined distribution type.
Edge 2
1 In the Model Builder window, right-click Mesh 1 and choose Edge.
Distribution 1
7 Click Build All.
STUDY 1

Dependent.
3 Click Range.
4 In the Range dialog box, type 100/freq in the Stop text field.
5 In the Step text field, type 2/freq.
6 Click Replace.
8 Click Range.
10 In the Start text field, type 101/13.56E6.
11 In the Stop text field, type 102/13.56E6.
13 Click Add.
The model can solve faster by deactivating the convergence plots. Also, the default plots
are unnecessary since custom plots will be used during results processing.
14 In the Model Builder window, click Study 1.
15 In the Settings window for Study, locate the Study Settings section.
16 Clear the Generate default plots check box.
17 Clear the Generate convergence plots check box.
RESULTS
1D Plot Group 1
2 In the Settings window for 1D Plot Group, locate the Data section.
3 From the Time selection list, choose From list.
4 In the Times (s) list, choose 7.448E-6, 7.467E-6, 7.485E-6, 7.504E-6, and 7.522E-6.
Line Graph 1
1 On the 1D Plot Group 1 toolbar, click Line Graph.
3 On the 1D Plot Group 1 toolbar, click Plot.
4 In the Settings window for Line Graph, click to expand the Legends section.
5 Select the Show legends check box.
6 From the Legends list, choose Manual.
Legends
0
pi/2
pi
3pi/2
2pi
8 Click to expand the Coloring and style section. Locate the Coloring and Style section. Find
the Line markers subsection. From the Marker list, choose Cycle.
1D Plot Group 1
1 In the Model Builder window, under Results click 1D Plot Group 1.
2 In the Settings window for 1D Plot Group, type Electron Density in the Label text
field.
3 Right-click Results>Electron Density and choose Duplicate.
Electron Density 1
In the Settings window for 1D Plot Group, type Electron Temperature in the Label text
field.
Line Graph 1
1 In the Model Builder window, expand the Electron Density 1 node, then click Results>
Electron Temperature>Line Graph 1.
2 In the Settings window for Line Graph, click Replace Expression in the upper-right corner
of the y-axis data section. From the menu, choose Component 1>Plasma (Drift Diffusion)>
Electron energy density>plas.Te - Electron temperature.
3 On the Electron Temperature toolbar, click Plot.
Electron Temperature
In the Model Builder window, under Results right-click Electron Temperature and choose
Duplicate.
Electron Temperature 1
In the Settings window for 1D Plot Group, type Electric Potential in the Label text
field.
Line Graph 1
1 In the Model Builder window, expand the Electron Temperature 1 node, then click
Results>Electric Potential>Line Graph 1.
2 In the Settings window for Line Graph, locate the Selection section.
3 From the Selection list, choose All domains.
4 Click Replace Expression in the upper-right corner of the y-axis data section. From the
menu, choose Component 1>Plasma (Electrostatics)>Electric>V - Electric potential.
5 On the Electric Potential toolbar, click Plot.
Electric Potential
In the Model Builder window, under Results right-click Electric Potential and choose
Duplicate.
Electric Potential 1
In the Settings window for 1D Plot Group, type Electric Field in the Label text field.
Line Graph 1
1 In the Model Builder window, expand the Electric Potential 1 node, then click Results>
Electric Field>Line Graph 1.
of the y-axis data section. From the menu, choose Component 1>Plasma (Electrostatics)>
Electric>Electric field>plas.Ex - Electric field, x component.
3 Locate the Selection section. In the list, select 2.
4 Select the Active toggle button.
5 Click Remove from Selection.
7 Click to expand the Quality section. From the Resolution list, choose No refinement.
8 On the Electric Field toolbar, click Plot.
Electric Field
In the Model Builder window, under Results right-click Electric Field and choose Duplicate.
Electric Field 1
In the Settings window for 1D Plot Group, type Ion Density in the Label text field.
Line Graph 1
1 In the Model Builder window, expand the Electric Field 1 node, then click Results>Ion
Density>Line Graph 1.
of the y-axis data section. From the menu, choose Component 1>Plasma (Heavy Species
Transport)>Number densities>plas.n_wAr_1p - Number density.
3 On the Ion Density toolbar, click Plot.
Ion Density
In the Model Builder window, under Results right-click Ion Density and choose Duplicate.
Ion Density 1
In the Settings window for 1D Plot Group, type Power Deposition in the Label text field.
Line Graph 1
1 In the Model Builder window, expand the Ion Density 1 node, then click Results>Power
Deposition>Line Graph 1.
Power and collisions>plas.Pcap - Capacitive power deposition.
3 On the Power Deposition toolbar, click Plot.
GEC ICP Reactor, Argon Chemistry
Introduction
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.
5 turn coil, 13.56 MHz (blue)
Dielectric material (gray)
Plasma (pink)
Wafer and pedestal
Figure 1: GEC ICP reactor geometry consisting of a 5 turn copper coil, plasma volume,
dielectrics, and wafer with pedestal.
Note: 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 m3). High density plasma sources are popular because low pressure ion
bombardment can provide a greater degree of anisotropy on the surface of the wafer.
2 | G E C I C P R E A C T O R, A R G O N C H E M I S T R Y
n
( ) + [ n e ( e E ) D e n e ] = R e
t e

( n ) + [ n ( E ) D n ] + E e = R
t
D e = e T e, = --- e, D = T e
5
3
Re = x j kj Nn ne
j=1
R = xj kj Nn ne j
j=1

kk = 0 k ( )f ( ) d
( w k ) + ( u )w k = j k + R k
t
0 r V =
N
= q Z k n k n e

k = 1
For a nonmagnetized, nonpolarized plasma, the induction currents are computed in the
frequency domain using the following equation:
2 1 e
( j 0 )A + ( 0 A ) = J
The plasma conductivity needs to be specified as a material property, usually from the cold
plasma approximation:
2
ne q
= -------------------------------
m e ( e + j )
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.
BOUNDARY CONDITIONS
n e = --- e, th n e
1
(1)
2
n = --- e, th n
5
(2)
6
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:

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:
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.
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 v.s. time in GEC ICP reactor.
Initially the power dissipated is all dissipated in the coil (~500 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 1600 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:
2
= -----------

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
Sci. Technol, vol. 14, pp. 722733, 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.
Application Library path: Plasma_Module/Inductively_Coupled_Plasmas/

argon_gec_icp
NEW
MODEL WIZARD
1 In the Model Wizard window, click 2D Axisymmetric.
2 In the Select Physics tree, select Plasma>Inductively Coupled Plasma.
3 Click Add.
4 Click Study.
5 In the Select Study tree, select Preset Studies for Selected Multiphysics>
Frequency-Transient.
6 Click Done.
Start by importing the geometry of the GEC reference cell.
GEOMETRY 1
Import 1 (imp1)
1 On the Home toolbar, click Import.
2 In the Settings window for Import, locate the Import section.
3 Click Browse.
argon_gec_icp.mphbin.
5 Click Import.
Add some predefined selections for the geometric entities which will be referenced later
on.
DEFINITIONS
Explicit 1
1 On the Definitions toolbar, click Explicit.
2 In the Model Builder window, right-click Explicit 1 and choose Rename.
3 In the Rename Explicit dialog box, type Walls in the New label text field.
4 Click OK.
5 In the Settings window for Explicit, locate the Input Entities section.
6 From the Geometric entity level list, choose Boundary.
7 Select Boundaries 6, 8, 3538, 44, 45, and 5156 only.
Explicit 2
3 In the Rename Explicit dialog box, type Coils in the New label text field.
4 Click OK.
Explicit 3
3 In the Rename Explicit dialog box, type Coil Boundaries in the New label text field.
4 Click OK.
6 In the Settings window for Explicit, locate the Output Entities section.
7 From the Output entities list, choose Adjacent boundaries.
GLOBAL DEFINITIONS
Parameters

Psp 1500[W] 1500 W Power input
mueN 4E24[1/(m*V*s)] 4E24 1/(Vms) Reduced electron
mobility
p0 0.02[torr] 2.666 Pa Gas pressure
PLASMA (PLAS)
1 In the Model Builder window, under Component 1 (comp1) click Plasma (plas).

3 Click Browse.
Ar_xsecs.txt.

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
Reaction 2
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
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 timestep,
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+
Ar+.
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 timesteps.
Initial Values 1
Values 1.
3 In the ne, 0 text field, type 1E15.
The Magnetic Fields are computed everywhere except the wafer and the wafer pedestal.
MAGNETIC FIELDS (MF)

On the Physics toolbar, click Plasma (plas) and choose Magnetic Fields (mf).
1 In the Model Builder window, under Component 1 (comp1) click Magnetic Fields (mf).
The Coil feature is used to electrically excite the system. The coil operates with a fixed
total power of 1500 watts.
3 In the Model Builder window, click Magnetic Fields (mf).
Coil 1
1 On the Physics toolbar, click Domains and choose Coil.
2 In the Settings window for Coil, locate the Domain Selection section.
3 From the Selection list, choose Coils.
4 Locate the Coil section. From the Coil excitation list, choose Power.
5 Select the Coil group check box.
6 In the Pcoil text field, type Psp.
PLASMA (PLAS)
Plasma Model 1
Model 1.
Next, define the material properties. There is no need to define the material properties
in the plasma domain, as these are defined by the Plasma Conductivity Coupling feature.
MATERIALS
Blank Material.
Material 1 (mat1)
1 In the Settings window for Material, locate the Geometric Entity Selection section.
3 Locate the Material Contents section. In the table, enter the following settings:

group
Electrical conductivity sigma 6e7 S/m Basic
Relative permeability mur 1 1 Basic
Material 2 (mat2)
1 Right-click Materials and choose Blank Material.

group
group
Electrical conductivity sigma 0 S/m Basic
Material 3 (mat3)

group
PLASMA (PLAS)
Surface reactions must always be included in a plasma model since they describe how ionic,
excited and radical species interact with the wall.
Surface Reaction 1
4 Locate the Boundary Selection section. From the Selection list, choose Walls.
Surface Reaction 2
Now, add boundary conditions to describe how the electrons interact with the wall.
Wall 1
2 In the Settings window for Wall, locate the General Wall Settings section.
3 In the re text field, type 0.2.
Ground 1
2 In the Settings window for Ground, locate the Boundary Selection section.
3 From the Selection list, choose Walls.
MESH 1
Meshing is a critical step in any plasma model. Boundary layer meshing on the reactor walls
is nearly always necessary. This is needed to capture the separation of space charge between
the electrons and ions close to the wall. You also add a fine mesh in the coil domains since
the skin depth needs to be resolved.
1 In the Model Builder window, under Component 1 (comp1) click Mesh 1.

3 From the Element size list, choose Finer.
Edge 1
1 Right-click Component 1 (comp1)>Mesh 1 and choose More Operations>Edge.
2 Select Boundaries 6, 8, 44, 45, and 54 only.
Size 1
1 Right-click Component 1 (comp1)>Mesh 1>Edge 1 and choose Size.
2 In the Settings window for Size, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Entire geometry.
4 Locate the Element Size section. Click the Custom button.
5 Locate the Element Size Parameters section. Select the Maximum element size check box.
6 In the associated text field, type 1E-3.
Free Triangular 1
1 In the Model Builder window, right-click Mesh 1 and choose Free Triangular.
2 In the Settings window for Free Triangular, locate the Domain Selection section.
Size 1
1 Right-click Component 1 (comp1)>Mesh 1>Free Triangular 1 and choose Size.
2 In the Settings window for Size, locate the Element Size section.
3 From the Predefined list, choose Extra fine.
Boundary Layers 1
1 In the Model Builder window, right-click Mesh 1 and choose Boundary Layers.
2 In the Settings window for Boundary Layers, locate the Domain Selection section.
Boundary Layer Properties

1 In the Model Builder window, under Component 1 (comp1)>Mesh 1>Boundary Layers 1
click Boundary Layer Properties.
2 In the Settings window for Boundary Layer Properties, locate the Boundary Selection
section.
4 Locate the Boundary Layer Properties section. In the Number of boundary layers text field,
type 5.
5 In the Boundary layer stretching factor text field, type 1.4.
Mapped 1
1 In the Model Builder window, right-click Mesh 1 and choose Mapped.
2 In the Settings window for Mapped, locate the Domain Selection section.
Distribution 1
1 Right-click Component 1 (comp1)>Mesh 1>Mapped 1 and choose Distribution.
2 In the Settings window for Distribution, locate the Boundary Selection section.
3 From the Selection list, choose Coil Boundaries.
distribution type.
10 In the Settings window for Mesh, click Build All.
STUDY 1
Step 1: Frequency-Transient
1 In the Settings window for Frequency-Transient, locate the Study Settings section.
2 In the Times text field, type 0 10^{range(-8,5/20,-3)}.
3 In the Frequency text field, type 13.56E6.
RESULTS
Electron Density (plas)

Click the Zoom Extents button on the Graphics toolbar.
Electron Temperature (plas)

Electric Potential (plas)

Now add a global plot for the coil resistance. This is defined as the real part of the total
voltage drop over the coil divided by the applied current. The Coil feature creates
pre-defined expressions for the resistance.
1D Plot Group 6
2 In the Settings window for 1D Plot Group, type Coil Resistance in the Label text
field.
3 Locate the Plot Settings section. Select the x-axis label check box.
4 In the associated text field, type Time (s).
5 Select the y-axis label check box.
6 In the associated text field, type Coil resistance (Ohm).
Global 1
1 On the Coil Resistance toolbar, click Global.
2 In the Settings window for Global, locate the y-Axis Data section.

mf.RCoil_1 Coil resistance
4 Click the x-Axis Log Scale button on the Graphics toolbar.

5 On the Coil Resistance toolbar, click Plot.
Now verify that 1500 watts is being applied to the system.
1D Plot Group 7
2 In the Settings window for 1D Plot Group, type Coil Power in the Label text field.
4 In the associated text field, type Time (s).
6 In the associated text field, type Power (W).
Global 1
1 On the Coil Power toolbar, click Global.

mf.PCoil_1 W Coil power

5 On the Coil Power toolbar, click Plot.
2D Plot Group 8
2 In the Settings window for 2D Plot Group, type Ion Number Density in the Label text
field.
Surface 1
1 Right-click Ion Number Density and choose Surface.
3 On the Ion Number Density toolbar, click Plot.
Ion Number Density

In the Model Builder window, under Results right-click Ion Number Density and choose
Duplicate.
Ion Number Density 1

In the Settings window for 2D Plot Group, type High Frequency Electric Field in
Surface 1
1 In the Model Builder window, expand the Results>High Frequency Electric Field node, then
click Surface 1.
the Expression section. From the menu, choose Component 1>Magnetic Fields>Electric>
mf.normE - Electric field norm.
3 On the High Frequency Electric Field toolbar, click Plot.
Observe that the electric field is slightly shielded by the plasma. This is due to the skin
effect in the plasma. As the electron number density increases, the plasma tends to shield
itself from the electric field.
High Frequency Electric Field

In the Model Builder window, under Results right-click High Frequency Electric Field and
choose Duplicate.
High Frequency Electric Field 1

In the Settings window for 2D Plot Group, type Excited Argon Number Density in the
Label text field.
Surface 1
1 In the Model Builder window, expand the Results>Excited Argon Number Density node,
Transport)>Number densities>plas.n_wArs - Number density.
3 On the Excited Argon Number Density toolbar, click Plot.
Study 1/Solution 1 (sol1)
1 In the Model Builder window, expand the Results>Data Sets node.
2 Right-click Study 1/Solution 1 (sol1) and choose Duplicate.
Selection
1 On the Results toolbar, click Selection.
2 In the Settings window for Selection, locate the Geometric Entity Selection section.
2D Plot Group 11
2 In the Settings window for 2D Plot Group, type Power Deposition in the Label text
field.
Surface 1
1 Right-click Power Deposition and choose Surface.
2 In the Settings window for Surface, locate the Data section.
3 From the Data set list, choose Study 1/Solution 1 (2) (sol1).
4 Click Replace Expression in the upper-right corner of the Expression section. From the
menu, choose Component 1>Magnetic Fields>Heating and losses>mf.Qrh - Volumetric loss
density, electric.
The effect of the shielding of the electric field due to the skin depth of the plasma is also
apparent when plotting the power deposition.
Argon Boltzmann Analysis
Introduction
The electron energy distribution function (EEDF) plays an important role in the overall
behavior of discharges. Analytic forms of the EEDF exist such as Maxwellian or
Druyvesteyn, but in some cases they fail to fit the discharge physics. Factors such as
ionization degree (ne/N), mole fraction of excited species, and reduced angular frequency
of the electromagnetic field can all influence the shape of the EEDF. This can lead to
differences in the magnitude of electron impact rate coefficients which are in turn used in
a fluid model for the electron density and mean electron energy. This tutorial example
shows how these factors influence the EEDF and the rate coefficients for argon.
Model Definition
This model is based on the results obtained in Ref 1. In this model, the x-coordinate
represents electron energy, (eV), rather than physical space. First a set of electron impact
collision cross sections are read from a file. The first time the model is solved, the EEDF
is assumed to be Maxwellian:
3 2
f() = 1 exp ( ( 2 ) )
where
32 5 2 1
1 = ( 5 2 ) (3 2) , 2 = ( 5 2 ) ( 3 2 )
Subsequently, the Two-term Boltzmann equation is solved which computes the EEDF
rather than assuming a specific form. The Two-term Boltzmann equation is:
Wf D f = S

where f is the EEDF (eV3/2) and
ne
W = 3a ------- A 1
2
(1)
Nn
and
E 2 k b T 2 ne
D = --- ------- -------- + ------------- + 2a ------- ( A 2 + A 3 )
32
(2)
3 Nn m q Nn
2 | ARGON BOLTZMANN ANALYSIS

For definitions of the quantities in the equations Equation 1 and Equation 2, see the
chapter The Boltzmann Equation, Two-Term Approximation Interface in the Plasma
Module Users Guide.
At zero energy, the condition that energy flux is zero must hold:
f
n Wf D = 0

and as , f 0. In practice, this boundary condition is applied by setting f to zero for

a suitably large . The default maximum is 100 eV by default.

Solving the Boltzmann equation allows various conditions for ionization degree, mole
fraction of excited species and reduced angular frequencies to be tested.
Figure 1: Maxwellian EEDF for a range of different mean electron energies.
Figure 1 shows the distribution function assuming Maxwellian. Electron-electron

collisions tend to make the shape of the EEDF closer to a Maxwellian distribution. The

influence of different degrees of ionization is shown in Figure 2.
Figure 2: EEDF in argon for an applied reduced electric field of 10 Td for different degrees
of ionization.
As the ionization degree increases, the EEDF tends to become more Maxwellian (which
would be represented by a linear line on a log scale). The electron-electron collisions tend
to push the electrons at the tail to higher energies, which has an effect on the rate
coefficients for electron impact reactions with a high activation energy.
The effect on the rate coefficients can be seen in Figure 3. The ionization rate coefficient
is substantially higher at a given mean electron energy when the ionization degree is
higher.

Figure 3: Plot of ionization rate coefficient versus mean electron energy for different
ionization degrees.
If the mole fraction of electronically excited species is non zero then superelastic collisions
can influence the shape of the EEDF. Figure 4 plots the EEDF for different mole fractions
of electronically excited argon. The superelastic collisions leads to a high energy tail in the
distribution function. This results in an increase in the magnitude of the ionization rate
coefficients.

Figure 4: EEDF for argon at a reduced electric field of 10 Td for different mole fractions of
electronically excited argon.
Figure 5: Ionization rate coefficient for different mole fractions of electronically excited
argon.

Finally, the influence of a high frequency oscillating field on the EEDF is studied by
applying different reduced angular frequencies. The EEDF is plotted in Figure 6.
Figure 6: Plot of EEDF for different reduced angular frequencies for the same mean electron
energy.
Reference
Science and Technology, vol. 14, pp. 722733, 2005.
Application Library path: Plasma_Module/Two-Term_Boltzmann_Equation/

boltzmann_argon

NEW
MODEL WIZARD
2 In the Select Physics tree, select Plasma>Boltzmann Equation, Two-Term Approximation
(be).
3 Click Add.
4 Click Study.
5 In the Select Study tree, select Preset Studies>Mean Energies.
6 Click Done.
GLOBAL DEFINITIONS
Begin by defining beta, the ionization degree, as a global parameter. This makes it available
for selection as a parameter in a parametric sweep.
Parameters

beta 1e-6 1E-6 Ionization degree
GEOMETRY 1
Notice that the geometry and mesh are automatically generated for this physics interface
and they very rarely need to be changed.
B O L T Z M A N N E Q U A T I O N , TWO - TE R M A P P R O X I M A T I O N ( B E )
1 In the Model Builder window, under Component 1 (comp1) click Boltzmann Equation,
Two-Term Approximation (be).
2 In the Settings window for Boltzmann Equation, Two-Term Approximation, locate the
Cross Section Import section.
3 Click Browse.
Ar_xsecs.txt.

Zero Probability 1
1 On the Physics toolbar, click Boundaries and choose Zero Probability.
Zero Energy Flux 1

1 On the Physics toolbar, click Boundaries and choose Zero Energy Flux.
STUDY 1
Step 1: Mean Energies

1 In the Model Builder window, expand the Study 1 node, then click Step 1: Mean Energies.
2 In the Settings window for Mean Energies, locate the Study Settings section.
3 In the Mean energies text field, type range(2,0.5,10).
RESULTS
EEDF (be)
1 In the Model Builder window, click EEDF (be).
2 In the Settings window for 1D Plot Group, click to expand the Axis section.
4 On the EEDF (be) toolbar, click Plot.
Electron Energy Distribution Function Settings section.
3 From the Electron energy distribution function list, choose Boltzmann.
4 Locate the Boltzmann Properties section. From the Electron-electron collisions list,
choose On.
Now parameterize the ionization degree of the argon.
Boltzmann Model 1
1 In the Model Builder window, under Component 1 (comp1)>Boltzmann Equation,
Two-Term Approximation (be) click Boltzmann Model 1.
2 In the Settings window for Boltzmann Model, locate the Mole Fraction Settings section.

3 In the xAr text field, type 1.
4 In the xArs text field, type 0.
5 Locate the Boltzmann Settings section. In the text field, type beta.
6 In the ne text field, type 1E18.
Initial Values 1
Two-Term Approximation (be) click Initial Values 1.
3 In the E/N0 text field, type 10[Td].
ADD STUDY
1 On the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select Study tree, select Preset Studies>Reduced
Electric Fields.
4 Click Add Study in the window toolbar.
5 On the Home toolbar, click Add Study to close the Add Study window.
STUDY 2
Step 1: Reduced Electric Fields

1 In the Model Builder window, under Study 2 click Step 1: Reduced Electric Fields.
2 In the Settings window for Reduced Electric Fields, locate the Study Settings section.
3 In the Reduced electric fields text field, type 10[Td].
4 Click to expand the Study extensions section. Locate the Study Extensions section. Select
the Auxiliary sweep check box.
5 Click Add.
Parameter name Parameter value list Parameter unit

beta 1e-6 1e-5 1e-4 1e-3 5e-3 1e-2

RESULTS
EEDF (be) 1
1 In the Model Builder window, under Results click EEDF (be) 1.
2 In the Settings window for 1D Plot Group, click to expand the Axis section.
3 Select the Manual axis limits check box.
4 In the x minimum text field, type -2.
5 In the x maximum text field, type 30.
6 In the y minimum text field, type 10^(-9).
7 In the y maximum text field, type 1.
Line Graph 1
1 In the Model Builder window, expand the EEDF (be) 1 node, then click Line Graph 1.
Legends
1e-6
1e-5
1e-4
1e-3
5e-3
1e-2
5 On the EEDF (be) 1 toolbar, click Plot.

Now parameterize the mole fraction of excited argon atoms.
ADD STUDY
3 Find the Studies subsection. In the Select Study tree, select Preset Studies>Reduced
Electric Fields.

STUDY 3
Step 1: Reduced Electric Fields

1 In the Model Builder window, under Study 3 click Step 1: Reduced Electric Fields.
2 In the Settings window for Reduced Electric Fields, locate the Study Settings section.
3 In the Reduced electric fields text field, type 10[Td].
4 Locate the Study Extensions section. Select the Auxiliary sweep check box.
5 Click Add.

beta 0 1e-5 1e-4 1e-3
Boltzmann Model 1
2 In the Settings window for Boltzmann Model, locate the Boltzmann Settings section.
3 In the text field, type 0.
4 Locate the Mole Fraction Settings section. In the xAr text field, type 1-beta.
5 In the xArs text field, type beta.
STUDY 3
On the Home toolbar, click Compute.
RESULTS
EEDF (be) 2
1 In the Model Builder window, click EEDF (be) 2.
2 In the Settings window for 1D Plot Group, locate the Axis section.

Line Graph 1
1 In the Model Builder window, expand the EEDF (be) 2 node, then click Line Graph 1.
2 In the Settings window for Line Graph, locate the Legends section.
Legends
0
1e-5
1e-4
1e-3

Now use the parametric sweep feature to compute the rate coefficients over a wide
range of mean electron energies for different degrees of ionization.
Boltzmann Model 1
3 In the text field, type beta.
4 Locate the Mole Fraction Settings section. In the xAr text field, type 1.
ADD STUDY
3 Find the Studies subsection. In the Select Study tree, select Preset Studies>Mean Energies.
STUDY 4

1 In the Model Builder window, under Study 4 click Step 1: Mean Energies.

Parametric Sweep
1 On the Study toolbar, click Parametric Sweep.
2 In the Settings window for Parametric Sweep, locate the Study Settings section.
3 Click Add.

beta 1e-6 1e-5 1e-4 1e-3 1e-2

8 On the Study toolbar, click Compute.
RESULTS
1D Plot Group 10
3 From the Data set list, choose Study 4/Parametric Solutions 1 (sol5).
Global 1
1 On the 1D Plot Group 10 toolbar, click Global.
the y-axis data section. From the menu, choose Component 1>Boltzmann Equation,
Two-Term Approximation>Rate coefficients>be.k_4 - Rate coefficient.
3 Locate the x-Axis Data section. From the Axis source data list, choose Inner solutions.
4 Click to expand the Legends section. Find the Include subsection. Select the Expression
check box.
5 Clear the Description check box.
1D Plot Group 10

3 Select the y-axis log scale check box.
4 Click to expand the Legend section. From the Position list, choose Lower right.
7 In the associated text field, type Mean electron energy (eV).
9 In the associated text field, type Ionization rate coefficient (m3/
s).
10 Locate the Axis section. Select the Manual axis limits check box.
12 In the y maximum text field, type 10^(-14).
Next, use the parametric sweep feature to compute the rate coefficients over a wide
range of mean electron energies for different mole fractions of excited argon
Boltzmann Model 1
4 Locate the Mole Fraction Settings section. In the xAr text field, type 1-beta.
5 In the xArs text field, type beta.
ADD STUDY

STUDY 5

4 Click to expand the Study extensions section.
Parametric Sweep
3 Click Add.

beta 0 1e-5 1e-4 1e-3

RESULTS
1D Plot Group 10
In the Model Builder window, under Results right-click 1D Plot Group 10 and choose
Duplicate.
1D Plot Group 11
Finally, compute the EEDF for a mean electron energy of 2.15 eV for different reduced
angular frequencies.

Boltzmann Properties section.
3 From the Oscillating field list, choose On.
Boltzmann Model 1
3 In the /N text field, type beta.
4 Locate the Mole Fraction Settings section. In the xAr text field, type 1.
ADD STUDY
STUDY 6

3 In the Mean energies text field, type 2.15.
Parametric Sweep
3 Click Add.

beta 0 1e-13 2e-13

RESULTS
EEDF (be) 3
2 In the Model Builder window, click EEDF (be) 3.
5 In the x minimum text field, type 0.

Oxygen Boltzmann Analysis
Introduction
The Boltzmann equation can be solved to validate sets of electron impact collision cross
sections. In fact, sets of collision cross sections are traditionally inferred by solving a
two-term approximation to the Boltzmann equation and comparing the results to swarm
experiments. This model solves a two-term approximation to the Boltzmann equation and
compares the computed drift velocity and electron temperature to experimental data.
Model Definition
The model uses molecular oxygen as the gas. For molecular oxygen, the model considers
the following set of electron impact collisions:
TABLE 1: TABLE OF COLLISIONS

1 e+O2=>e+O2 Momentum 0
2 e+O2=>O+O Attachment 0
3 e+O2=>e+O2(rot) Excitation 0.02
4 e+O2=>e+O2(v=1) Excitation 0.19
10 e+O2=>e+O2(a1d) Excitation 0.977
11 e+O2(a1d)=>e+O2 Superelastic -0.977
12 e+O2=>e+O2(b1s) Excitation 1.627
13 e+O2(b1s)=>e+O2 Superelastic -1.627
14 e+O2=>e+O2(45) Excitation 4.5
15 e+O2(45)=>e+O2 Superelastic -4.5
16 e+O2=>e+O+O Dissociation 6.0
17 e+O2=>e+O+O(1D) Dissociation 8.4
18 e+O2=>e+O+O(1S) Dissociation 9.97
19 e+O2=>e+O2+ Ionization 12.06
2 | OXYGEN BOLTZMANN ANALYSIS

In a superelastic collision, the electrons gain energy from excited species. The mole
fraction of each species is given in the table below and is estimated from typical discharge
conditions that occur in a drift tube. The degree of ionization is 106.
TABLE 2: TABLE OF MOLE FRACTIONS OF EACH SPECIES
SPECIES MOLE FRACTION
O2 0.99997
O2(a1d) 1.5E-5
O2(b1s) 1E-5
O2(45) 5E-6
The Two-term Boltzmann equation is solved which computes the EEDF rather than
assuming a specific form. The two-term approximation of the Boltzmann equation is:
Wf D f = S

where f is the EEDF (eV3/2) and
ne
W = 3a ------- A 1
2
(1)
N n
and
E 2 k b T 2 ne
D = --- ------- -------- + ------------- + 2a ------- ( A 2 + A 3 )
32
(2)
3 Nn m q Nn
For definitions of the quantities in the equations Equation 1 and Equation 2, see the
chapter The Boltzmann Equation, Two-Term Approximation Interface in the Plasma
Module Users Guide.
At zero energy, the condition that energy flux is zero must hold:
f
n Wf D = 0

and as , f 0. In practice, this boundary condition is applied by setting f to zero for

a suitably large . The default maximum is 100 eV by default.
When the Boltzmann equation has been the drift velocity and electron temperature for
a given reduced electric field can be compared to experimental results. The drift velocity
is defined as:

f
w = --- -------
E - d
3 N n 0 -------
m
which is linear in the reduced electric field by has a rather complicated integral dependence
on the EEDF. The electron temperature is simply two thirds of the mean electron
energy which is defined as:

32
= 0 f d
The experimental data comes from Ref. 2.

The two-term approximation to the Boltzmann equation is solved by parametrically
varying the mean electron energy. This parametric analysis allows the EEDF to be plotted
for different mean electron energies. Figure 1 plots the EEDF for different values of the
mean electron energy. The mean electron energy has its lowest value indicated by the blue
line on the left of the plot, corresponding to a mean electron energy of 2 electron volts.
At such a low energy the population of electrons with energy above the ionization
threshold is very low. As the mean electron energy increases, the population of electrons
with higher energy increases. This makes ionization processes more likely. Ionization is
important as it is usually the primary mechanism for sustaining plasmas. Notice also that

the EEDF is not linear on the log scale. This indicates that the EEDF is non-Maxwellian
under these conditions.
Figure 1: Plot of the EEDF for different values of the mean electron energy.
The reduced electron transport properties are plotted in Figure 2. The transport
properties have a much weak dependence on the EEDF compared to rate or Townsend
coefficients. The electron mobility and electron energy mobility decrease as the mean
electron energy increases. The electron diffusivity and electron energy diffusivity increase
as the mean electron energy increases. If the EEDF was Maxwellian then the following
relations would hold:
D e = e T e, = --- e, D = T e
5
3
In the case that the EEDF is non-Maxwellian, this relation does not necessarily hold true.
The drift velocity of the electrons at a given mean electron energy or reduced electric field
is an important parameter because it is easy to measure experimentally in drift tubes.
Traditionally, the electron drift velocity is tabulated v.s. the reduced electric field (E/N).
The simulated and experimental drift velocity is plotted in Figure 3 for a range of reduced
electric fields. The agreement between the two is good over a wide range of reduced
electric field, indicating that the cross section data is consistent with experimental
measurements.

Figure 2: Reduced transport properties v.s. mean electron energy.
Another quantity which is relatively easy to measure experimentally is the ratio of the
electron diffusivity and the mobility. The simulated and experimental ratio of the
diffusivity and mobility is plotted in Figure 4. Einsteins relation dictates that this ratio is
equal to the electron temperature when the EEDF is Maxwellian.

Figure 3: Computed and experimental drift velocity for oxygen.
Figure 4: Computed and experimental D/ for different reduced electric fields.

The fraction of power channeled into various reactions is shown in Figure 5. This is
important in design of plasma sources because often the desired reactive species is known
in advance. The plot gives an indication of what the target mean electron energy should
be in order to channel as much of the available power into a specific reaction. Of course,
the power channeled into ionization must be high enough to sustain the plasma.
Figure 5: Plot of the fraction of the total power channeled into each reaction v.s. mean electron
energy.
The Townsend coefficients are plotted in Figure 6. The Townsend coefficients offer an
alternative way of defining reaction rates. The reaction rate depends on the electron flux
rather than the electron density. Townsend coefficients should be used when modeling DC
discharges.

Figure 6: Townsend coefficients v.s. mean electron energy for oxygen.
References
Science and Technology, vol. 14, pp. 722733, 2005.
2. J. Dutton, A Survey of Electron Swarm Data, J. Phys. Chem. Ref. Data, vol. 4,
pp. 577866, 1975.
Application Library path: Plasma_Module/Two-Term_Boltzmann_Equation/

boltzmann_oxygen

NEW
MODEL WIZARD
2 In the Select Physics tree, select Plasma>Boltzmann Equation, Two-Term Approximation
(be).
3 Click Add.
4 Click Study.
5 In the Select Study tree, select Preset Studies>Mean Energies.
6 Click Done.
GEOMETRY 1
Notice that COMSOL generates the geometry and mesh for the Boltzmann Equation,
Two-Term Approximation interface. The geometry and mesh almost never need to be
modified.
Cross Section Import section.
3 Click Browse.
O2_xsecs.txt.
5 Locate the Electron Energy Distribution Function Settings section. From the Electron
energy distribution function list, choose Boltzmann.
6 Locate the Boltzmann Properties section. From the Electron-electron collisions list,
choose On.
Boltzmann Model 1
1 In the Model Builder window, expand the Boltzmann Equation, Two-Term Approximation
(be) node, then click Boltzmann Model 1.
2 In the Settings window for Boltzmann Model, locate the Mole Fraction Settings section.
3 In the xO2a1d text field, type 1.5E-5.
4 In the xO245 text field, type 5E-6.

Initial Values 1
Two-Term Approximation (be) click Initial Values 1.
3 In the 0 text field, type 0.25.
Zero Energy Flux 1

1 On the Physics toolbar, click Boundaries and choose Zero Energy Flux.
Zero Probability 1
1 On the Physics toolbar, click Boundaries and choose Zero Probability.
DEFINITIONS
Interpolation 1 (int1)
1 On the Home toolbar, click Functions and choose Local>Interpolation.
2 In the Settings window for Interpolation, locate the Definition section.
3 From the Data source list, choose File.
4 Click Browse.
O2_drift_velocity_expt.txt.
6 Click Import.
7 Locate the Units section. In the Arguments text field, type Td.
8 In the Function text field, type m/s.
Interpolation 2 (int2)
1 On the Home toolbar, click Functions and choose Local>Interpolation.
2 In the Settings window for Interpolation, locate the Definition section.
3 From the Data source list, choose File.
4 Click Browse.
O2_Te_expt.txt.
6 Click Import.

7 Locate the Interpolation and Extrapolation section. From the Extrapolation list, choose
Linear.
8 Locate the Units section. In the Arguments text field, type Td.
9 In the Function text field, type V.
STUDY 1

1 In the Model Builder window, expand the Study 1 node, then click Step 1: Mean Energies.
3 In the Mean energies text field, type range(0.25,0.25,10).
Solution 1 (sol1)
1 On the Study toolbar, click Show Default Solver.
2 In the Model Builder window, expand the Study 1>Solver Configurations node.
3 In the Model Builder window, expand the Solution 1 (sol1) node, then click Stationary
Solver 1.
4 In the Settings window for Stationary Solver, locate the General section.
5 From the Linearity list, choose Nonlinear.
RESULTS
EEDF (be)
1 In the Model Builder window, under Results click EEDF (be).
3 From the Parameter selection (freq) list, choose From list.
4 In the Parameter values (freq (V)) list, choose 2, 3, 4, 5, and 6.
5 Click to expand the Axis section. Select the Manual axis limits check box.
8 In the y minimum text field, type 1e-8.
10 On the EEDF (be) toolbar, click Plot.

1D Plot Group 4
2 In the Settings window for 1D Plot Group, locate the Plot Settings section.
3 Select the x-axis label check box.
4 In the associated text field, type Reduced electric field (Td).
6 In the associated text field, type Drift velocity (m/s).
7 Click to expand the Legend section. From the Position list, choose Upper left.
Global 1

be.w m/s Computed drift velocity
4 Locate the x-Axis Data section. From the Parameter list, choose Expression.
5 Locate the x-axis data section. Click be.EN - Reduced electric field in the upper-right
corner of the section. Locate the x-Axis Data section. From the Unit list, choose Td.
1D Plot Group 4
In the Model Builder window, under Results click 1D Plot Group 4.
Global 2

int1(be.EN) Experimental drift velocity
6 Click to expand the Coloring and style section. Locate the Coloring and Style section. Find
the Line style subsection. From the Line list, choose None.

7 Find the Line markers subsection. From the Marker list, choose Circle.
8 In the Number text field, type 25.
9 In the Width text field, type 3.
1D Plot Group 5
4 In the associated text field, type Reduced electric field (Td).
6 In the associated text field, type D/\mu (eV).
Global 1

be.DeN/be.muN Wb Computed D/\mu
1D Plot Group 5
2 In the Settings window for 1D Plot Group, locate the Legend section.
3 From the Position list, choose Upper left.
Global 2


int2(be.EN) Experimental D/\mu
6 Locate the Coloring and Style section. Find the Line style subsection. From the Line list,
choose None.
7 Find the Line markers subsection. From the Marker list, choose Diamond.
8 In the Number text field, type 25.
9 In the Width text field, type 3.
1D Plot Group 6
Global 1

be.xi_1 1 Fractional power loss

choose Cycle.
5 Click to expand the Legends section. Find the Include subsection. Select the Expression
check box.
1D Plot Group 7
5 Locate the Axis section. Select the Manual axis limits check box.
6 In the x minimum text field, type -0.2375.
8 In the y minimum text field, type 1e-24.
9 In the y maximum text field, type 4.052e-18.
Global 1

be.alpha_1 m^2 Townsend coefficient

choose Cycle.
5 Locate the Legends section. Find the Include subsection. Select the Expression check box.
7 Click the y-Axis Log Scale button on the Graphics toolbar.

Benchmark Model of a Capacitively Coupled

Plasma
Introduction
The underlying physics of a capacitively coupled plasma is rather complicated, even for
rather simple geometric configurations and plasma chemistries. This model benchmarks
the COMSOL Multiphysics Plasma Module Capacitively Coupled Plasma interface against
many different codes, the results of which are taken from Ref. 1.
Model Definition
The model geometry consists of a 1D gap of 0.067 m. A plasma forms in the gap provided
the driving voltage and fill pressure are high enough. The driving frequency in this model
is 13.56 MHz. Helium chemistry is used, as was the case in Ref. 1.
n
( ) + [ n e ( e E ) D e n e ] = R e
t e
n
( ) + [ n ( E ) D n ] + E e = R
t
D e = e T e, = --- e, D = T e
5
3
Re = x j kj Nn ne
j=1
2 | BENCHMARK MODEL OF A CAPACITIVELY COUPLED PLASMA

R = xj kj Nn ne j
j=1

kk = 0 k ( )f ( ) d

( w ) + ( u )w k = j k + R k
t k
0 r V =
N
= q Z k n k n e
k = 1
Using Ref. 1, the ion diffusivity is defined as:
1
---
3 2
D = 0,10901 ( a )

and the mobility, as:
1
---
2
E = 0,79788 ( a )
where a = qE m and is the charge exchange mean free path.
BOUNDARY CONDITIONS
condition for the electron flux
n e = ------------ --- e, th n e
1r 1
(1)
1 + r 2
and the electron energy flux
n = ------------ --- e, th n
1r 5
(2)
1 + r 6
In order to make the COMSOL implementation of the electron losses to the wall
consistent with the reference, the value of r must be set to 5/11. The second term on the
right-hand side of Equation 1 is the gain of electrons due to secondary emission effects,
p being the secondary emission coefficient. The second term in Equation 2 is the
secondary emission energy flux, p being the mean energy of the secondary electrons. For
the helium ions are lost to the wall due to surface reactions and the fact that the electric
field is directed towards the wall:
n j k = M w c k Z k ( E n ) [ Z k k ( E n ) > 0 ]
Surface charge accumulation is added to the dielectric surface by way of the following
boundary condition:
n ( D1 D 2 ) = s
d s
= n Ji + n Je
dt

The discharge is driven by a sinusoidal total current applied to the exterior boundary of
the dielectric plate:
I = I 0 cos ( t )
where is the angular frequency. This boundary condition is not possible to implement
directly - COMSOL requires that a voltage is specified on the driven electrode. In order
to be able to apply a voltage to the driven electrode, a differential algebraic equation is
added to the model. This differential algebraic equation then adjusts the voltage on the
electrode at each time step such that the above current constraint is satisfied. The equation
solved for the unknown electrode potential, Vin, is hence, in weak form:
D I cos ( t ) V in = 0
t 0
where the bar denotes test function. This is implemented as a bidirectional weak constraint
in COMSOL.
PLASMA CHEMISTRY
The reference paper suggests a simplistic plasma chemistry for helium consisting of only 3
reactions and 4 species:
TABLE 1: TABLE OF COLLISIONS AND REACTIONS MODELED.

1 e+He=>e+He Elastic 0
2 e+He=>e+Hes Excitation 19.5
3 e+Hes=>e+He+ Ionization 24.5
In addition to volumetric reactions, the following surface reactions are implemented:

TABLE 2: TABLE OF SURFACE REACTIONS.
1 Hes=>He 1
2 He+=>He 0
When a metastable helium atom makes contact with the wall, it reverts to the ground state
helium atom with some probability (the sticking coefficient).

The time averaged ion current density is plotted in Figure 1. The peak ion current density
occurs on the electrodes, at a value of around 0.2A/m2. At lower pressures the ion current
density profile is more smooth across the gap. At 300mtorr there is a pronounced
flattening off of the ion current density in the plasma sheath. These results agree well with
those presented in Figure 9 of Ref. 1.
Figure 1: Plot of the time averaged ion current as a function of the distance from the left
electrode.
The time averaged excitation and ionization rates are plotted in Figure 2 and Figure 3.

Figure 2: Plot of the time averaged ionization rate as a function of the distance from the left
electrode.
Figure 3: Plot of the time averaged excitation rate as a function of the distance from the left
electrode.

The results for the ionization and excitation rates agree well with Ref. 1 in both absolute
value and spatial distribution. The time averaged electron power deposition is plotted in
Figure 4. The COMSOL results are again, in good agreement with the reference paper.
Figure 4: Plot of the time averaged electron power deposition as a function of the distance from
the left electrode.
Figure 5 and Figure 6 plot the time averaged electron and ion density at different
operating pressures. The electron and ion density is the same in the plasma bulk but the
ion density is higher in the plasma sheath. This creates a net positive space charge density
in the plasma sheath which tends to hold electrons in the plasma and accelerate ions
towards the wall via the plasma potential.

Figure 5: Plot of the time averaged electron density at different operating pressures.
Figure 6: Plot of the time averaged ion density for different operating pressures.

Ref. 1 tabulates a wide range of lumped parameters for the discharge. These same
parameters for the COMSOL model are shown in Table 3. There is good agreement
between the COMSOL model and the values quoted in Ref. 1.
TABLE 3: RESULTS FROM THE DISCHARGE MODEL
3 2 2 2
P (MTORR) V rf ( V ) n e ( cm ) ( eV ) P tot ( mWcm ) j pos ( mAcm ) I rf ( mAcm )
30 392.5 1.38x108 14.7 2.92 0.019 1.00
100 220.0 6.19x108 7.3 4.0 0.0195 1.00
300 165.0 1.25x109 5.5 6.46 0.0175 1.00
Reference
1. M. Surendra, Radiofrequency discharge benchmark model comparison, Plasma
Sources Sci. Technol, vol. 4, pp 5673, 1995.
Application Library path: Plasma_Module/Capacitively_Coupled_Plasmas/

ccp_benchmark
NEW
MODEL WIZARD
3 Click Add.
4 Click Study.
6 Click Done.
7 In the Model Builder windows toolbar, click the Show button and select Advanced Physics
Options in the menu.

8 In the Model Builder windows toolbar, click the Show button and select Equation View
in the menu.
GEOMETRY 1
Interval 1 (i1)
Interval 2 (i2)
GLOBAL DEFINITIONS
Parameters

p0 0.3[torr] 40 Pa Gas pressure
freq 13.56[MHz] 1.356E7 Hz RF frequency
omega 2*pi*freq 8.52E7 Hz Angular frequency
J0 1[mA/cm^2] 10 A/m Current density
DEFINITIONS
Variables 1


I_sp J0*cos(omega*t)* A Current setpoint
1[m^2]
mue 3.33E24[1/(V*m*s)]/ m/(Vs) Electron mobility
plas.Nn
a e_const* Help variable
abs(plas.Ex)/
(plas.M_wHe/
N_A_const)
sigma 3e-15[cm^2] m Charge exchange cross
section
lambda 1/(plas.Nn*sigma) m Charge exchange mean
free path
arg2 a*lambda^3 Help variable
D_He_1p 0.10901*sqrt(arg2) Helium ion
diffusivity
v_He_1p 0.79788*sqrt(a* Helium ion drift
lambda) velocity
Integration 1 (intop1)
1 On the Definitions toolbar, click Component Couplings and choose Integration.
2 In the Settings window for Integration, locate the Source Selection section.
PLASMA (PLAS)
1 In the Model Builder window, expand the Component 1 (comp1)>Plasma (plas) node, then
click Plasma (plas).
2 In the Settings window for Plasma, locate the Transport Settings section.
3 From the Diffusion model list, choose Ficks law.
4 Locate the Out-of-Plane Thickness section. In the A text field, type 1[m^2].
5 In the Model Builder windows toolbar, click the Show button and select Stabilization in
the menu.
6 Click to expand the Stabilization section. Clear the Reaction source stabilization check
box.

7 Clear the Source stabilization check box.

3 Click Browse.
He_xsecs.txt.
Species: He
1 In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click Species: He.
4 Locate the General Parameters section. In the Df text field, type 0.8.
5 From the Preset species data list, choose He.
Species: Hes
Hes.
3 In the x0 text field, type 1E-6.
4 In the Df text field, type 0.8.
5 From the Preset species data list, choose He.
Species: He+
He+.
4 Locate the General Parameters section. From the Preset species data list, choose He.
5 In the Df text field, type D_He_1p.
Plasma Model 1
Model 1.

5 Locate the Electron Density and Energy section. In the e text field, type mue.
Surface Reaction 1
4 In the Formula text field, type Hes=>He.
Surface Reaction 2
4 In the Formula text field, type He+=>He.
5 Locate the Reaction Parameters section. In the f text field, type 0.
MATERIALS
Blank Material.
Material 1 (mat1)

group
PLASMA (PLAS)
Initial Values 1
Values 1.

5 In the V text field, type 1e-6[V/m]*x.
Wall 1
4 In the re text field, type 5/11.
Ground 1

Terminal 1
4 From the Terminal type list, choose Current.
5 In the I0 text field, type I_sp.
6 Click to expand the Constraint settings section. Locate the Constraint Settings section.
Select the Use weak constraints check box.
Species: He+
1 In the Model Builder window, expand the Species: He+ node, then click Equation View.
2 In the Settings window for Equation View, locate the Variables section.
Name Expression Unit Description Selection Details

plas.um_ v_He_1p/ m/(Vs) Mobility Domain 1
wHe_1p sqrt(Vx^2+
eps)

MESH 1
Edge.
Edge 1
Distribution 1
8 In the Model Builder window, right-click Mesh 1 and choose Edge.
Edge 2
choose Size.
Size 1
2 From the Predefined list, choose Extremely fine.
3 Click Build All.
GEOMETRY 1
MESH 1

STUDY 1

Dependent.
3 Click Range.
5 In the Start text field, type 0.
6 In the Stop text field, type 200/13.56e6.
8 Click Replace.
10 Click Range.
11 In the Range dialog box, type 200/13.56e6 in the Start text field.
14 Click Add.
17 Clear the Generate convergence plots check box.
Parametric Sweep
3 Click Add.

p0 0.03[torr] 0.1[torr] 0.3[torr]

RESULTS
In the Model Builder window, expand the Results node.
1D Plot Group 1
2 In the Settings window for 1D Plot Group, type Time Averaged Ion Current
Density in the Label text field.
3 Locate the Data section. From the Data set list, choose Study 1/Parametric Solutions 1
(sol2).
6 In the associated text field, type Distance (m).
8 In the associated text field, type Time averaged ion current density (A/m
2).
9 Click to expand the Title section. From the Title type list, choose None.
Line Graph 1
1 On the Time Averaged Ion Current Density toolbar, click Line Graph.
3 In the Settings window for Line Graph, locate the y-Axis Data section.
4 In the Expression text field, type timeavg(200/freq,201/freq,
pprint(plas.Jix_wHe_1p),1e-3).
5 Click to expand the Legends section. Select the Show legends check box.
Legends
0.03 torr
0.1 torr
0.3 torr
8 On the Time Averaged Ion Current Density toolbar, click Plot.


Time Averaged Ion Current Density
In the Model Builder window, under Results right-click Time Averaged Ion Current Density
and choose Duplicate.
Time Averaged Ion Current Density 1

1 In the Settings window for 1D Plot Group, type Time Averaged Ionization Rate in
2 Locate the Plot Settings section. In the x-axis label text field, type Distance (m).
3 In the y-axis label text field, type Time averaged ionization rate (1/(m3s)).
Line Graph 1
1 In the Model Builder window, expand the Results>Time Averaged Ionization Rate node,
then click Line Graph 1.
3 In the Expression text field, type timeavg(200/freq,201/freq,plas.r_3*
N_A_const,1e-3).
4 On the Time Averaged Ionization Rate toolbar, click Plot.

Time Averaged Ionization Rate

In the Model Builder window, under Results right-click Time Averaged Ionization Rate and
choose Duplicate.
Time Averaged Ionization Rate 1

1 In the Settings window for 1D Plot Group, type Time Averaged Excitation Rate in
3 In the y-axis label text field, type Time averaged excitation rate (1/(m3s)).
Line Graph 1
1 In the Model Builder window, expand the Results>Time Averaged Excitation Rate node,
3 In the Expression text field, type timeavg(200/freq,201/freq,plas.r_2*
N_A_const,1e-3).
4 On the Time Averaged Excitation Rate toolbar, click Plot.

Time Averaged Excitation Rate

In the Model Builder window, under Results right-click Time Averaged Excitation Rate and
choose Duplicate.
Time Averaged Excitation Rate 1

1 In the Settings window for 1D Plot Group, type Time Averaged Power Deposition
3 In the y-axis label text field, type Time averaged power deposition (W/m3).
Line Graph 1
1 In the Model Builder window, expand the Results>Time Averaged Power Deposition node,
3 In the Expression text field, type timeavg(200/freq,201/freq,
pprint(plas.Pcap),1e-3).
4 On the Time Averaged Power Deposition toolbar, click Plot.
Time Averaged Power Deposition

In the Model Builder window, under Results right-click Time Averaged Power Deposition and
choose Duplicate.
Time Averaged Power Deposition 1

1 In the Settings window for 1D Plot Group, type Time Averaged Electron Density
2 Locate the Plot Settings section. In the y-axis label text field, type Time averaged
electron density (1/m3).
Line Graph 1
1 In the Model Builder window, expand the Results>Time Averaged Electron Density node,
3 In the Expression text field, type timeavg(200/freq,201/freq,plas.ne,1e-3).
4 On the Time Averaged Electron Density toolbar, click Plot.

Time Averaged Electron Density
In the Model Builder window, under Results right-click Time Averaged Electron Density and
choose Duplicate.
Time Averaged Electron Density 1

1 In the Settings window for 1D Plot Group, type Time Averaged Ion Density in the
Label text field.
2 Locate the Plot Settings section. In the y-axis label text field, type Time averaged ion
density (1/m3).
Line Graph 1
1 In the Model Builder window, expand the Results>Time Averaged Ion Density node, then
click Line Graph 1.
3 In the Expression text field, type timeavg(200/freq,201/freq,plas.n_wHe_1p,
1e-3).
4 On the Time Averaged Ion Density toolbar, click Plot.
Cut Point 1D 1
1 On the Results toolbar, click More Data Sets and choose Cut Point 1D.
2 In the Settings window for Cut Point 1D, locate the Data section.
4 Locate the Point Data section. In the x text field, type 0.067/2.
1D Plot Group 7
2 In the Settings window for 1D Plot Group, type Electric Potential on Electrode
3 Locate the Data section. From the Data set list, choose Study 1/Parametric Solutions 1
(sol2).
4 From the Time selection list, choose Interpolated.
5 Click Range.
7 In the Start text field, type 200/13.56e6.

10 Click Replace.
Point Graph 1
1 On the Electric Potential on Electrode toolbar, click Point Graph.
3 In the Settings window for Point Graph, locate the y-Axis Data section.
4 In the Expression text field, type V.
5 On the Electric Potential on Electrode toolbar, click Plot.
Point Evaluation 1
1 On the Results toolbar, click Point Evaluation.
2 In the Settings window for Point Evaluation, locate the Data section.
5 From the Table columns list, choose Time.
6 From the Data set list, choose Cut Point 1D 1.
7 Click Replace Expression in the upper-right corner of the Expressions section. From the
menu, choose Component 1>Plasma (Drift Diffusion)>Electron density>plas.ne - Electron
density.

plas.ne 1/cm^3 Electron density
9 Click Evaluate.
Point Evaluation 2
1 Right-click Point Evaluation 1 and choose Duplicate.
2 In the Settings window for Point Evaluation, click Replace Expression in the upper-right
corner of the Expressions section. From the menu, choose Component 1>Plasma (Drift
Diffusion)>Electron energy density>plas.ebar - Mean electron energy.
3 Click Table 1 - Point Evaluation 1 (plas.ne).

Point Evaluation 3
1 On the Results toolbar, click Point Evaluation.
2 In the Settings window for Point Evaluation, locate the Data section.
menu, choose Component 1>Plasma (Heavy Species Transport)>Current>plas.nJt - Total
ion current density on wall.

plas.nJt mA/cm^2 Total ion current density on wall
9 Click Table 1 - Point Evaluation 1 (plas.ne).
Global Evaluation 1
1 On the Results toolbar, click Global Evaluation.
4 From the Time selection list, choose Interpolated.
5 Click Range.
7 In the Start text field, type 200/13.56e6.
10 Click Replace.
menu, choose Component 1>Plasma>Power and collisions>plas.Pcap_tot - Total capacitive
power deposition.


1e-4*plas.Pcap_tot mW
15 Locate the Data Series Operation section. From the Operation list, choose Average.
16 Click New Table.

Atmospheric Pressure Corona Discharge
Introduction
This model simulates a corona discharge occurring between two co-axially fashioned
conductors. The negative electric potential is applied to the inner conductor and the
exterior conductor is grounded. The discharge gas simulated is argon at atmospheric
pressure.
Model Definition
Figure 1 shows a cross section of the model. By considering a long and uniform coaxial
conductor configuration, the model can be viewed as axisymmetric and thus simplified to
a 1D problem.
The model presented in the following section is used to simulate the ionization of the
neutral gas (Ar) as well as the flux of charged particles (Ar+ and electrons) when the
negative electric potential is applied at the inner conductor (cathode). The high electric
field generated by the combination of high potential and small conductor curvature radius
(inner conductor, ri) causes electron drift and ionization of the neutral gas surrounding
the cathode. The resulting ions generate more electrons through secondary emission at
the cathode surface. These electrons are accelerated through a small region away from the
cathode, where they can acquire significant energy. This can lead to ionization which
creates new electron-ion pairs. The secondary ions migrate towards the cathode where
they eject more secondary electrons. This process is responsible for sustaining the
discharge.
The model is based on the fluid equations for electrons and ions as well as on Poissons
equation. Secondary electrons generated by ion bombardment of the cathode surface are
taken into account. The model uses a Scharfetter-Gummel upwind scheme to eliminate
numerical instabilities in the number density of the charged particles associated with the
finite element method. This is needed, in particular close to the cathode, where the ion
flux is particularly high.
2 | ATMOSPHERIC PRESSURE CORONA DISCHARGE

simulated
cathode (-Vin) domain (1D)
anode
ro
ri
Figure 1: Not-to-scale cross section of the co-axial configuration. The negative potential (-Vin)
is applied at the inner conductor (cathode) and the outer electrode is grounded (anode). The
shaded area represents the ionization region created by the positive space charge distribution
generated in the vicinity of the cathode.
DOMAIN EQUATIONS
The electron density and mean electron energy are computed by solving the drift-diffusion
equations for the electron density and mean electron energy. Convection of electrons due
to fluid motion is neglected. For more detailed information on electron transport see the
section Theory for the Drift Diffusion Interface in the Plasma Module Users Guide.
( n ) + [ n ( E ) D n ] = R
t e e e e e e
n
( ) + [ n ( E ) D n ] + E e = R
t
electron mobility using the relations
5
D e = e T e, = --- e, D = T e
3

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
Re = x j kj Nn ne
j=1
R = xj kj Nn ne j
j=1
Here j is the energy loss from reaction j (SI unit: V). The rate coefficients can be
computed from cross-section data as the integrals

kk = 0 k ( )f ( ) d
the energy (SI unit: V), k is the collision cross section (SI unit: m2), and f is the electron
energy distribution function. In this case, a Maxwellian EEDF is assumed.
species:
( w k ) + ( u )w k = j k + R k
t
For detailed information on the transport of the non-electron species see the section
Theory for the Heavy Species Transport Interface in the Plasma Module Users Guide.
The electrostatic field is computed using the equation
0 r V =
specified in the model using the formula

N
= q Z k n k n e
k = 1
Boundary Conditions
wall and gained due to secondary emission effects, resulting in the boundary condition
n e = --- e, th n e
1
2 p ( p n ) (1)
p
for the electron flux and
n = --- e, th n
5
6 p p ( p n ) (2)
p
for the electron energy flux. The second term on the right-hand side of Equation 1 is the
gain of electrons due to secondary emission effects, p being the secondary emission
coefficient. The second term in Equation 2 is the secondary emission energy flux, p being
the mean energy of the secondary electrons. For the heavy species, ions are lost to the wall
due to surface reactions and the fact that the electric field is directed towards the wall:
The discharge is driven by the electric potential applied to the inner conductor of the
coaxial geometry (at coordinate r = ri).
V = V 0 tanh --
t

where the function tanh(t/) is used to generate the voltage step function (1000 V). The
other boundary (at coordinate r = ro) is grounded.
Note that, during the simulation, the cathode is submitted to an intense flux of ions that
generate an important amount of secondary electrons which, in turn, increase the cathode
ion bombardment and so on. In order to prevent this avalanche from increasing
indefinitely, a RC circuit has been added in series with the system. In order to model the

circuit addition, the voltage at the inner conductor is modified using the differential
equation
V
V = V0 Ip Rb + Rb Cb
t
where Ip is defined as

I p = n J i + n J e + ( n D ) dS

t
and where n Ji is the normal component of the total ion current density at the wall,
n Je is the normal component of the total electron current density at the wall, and n D
the normal electrical displacement at the surface.
PLASMA CHEMISTRY
Argon is an attractive gas to use in a benchmark problem because only a handful of
reactions and a few species need to be considered. Table 1 lists the chemical reactions
considered.
Reaction Formula Type (eV) kf (m3/(smol))

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 Ar+Ars=>Ar+Ar Reaction - 1807
7 Ars+Ars=>Ars+Ar Reaction - 2.3107
Initially a small number of seed electrons are present. These are necessary in order to
initiate the discharge. In addition to the 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 the ion atoms reach the wall, they are assumed to change back to neutral argon
atoms and donate their charge to the wall. Note that the secondary emission coefficient is

set to 0.2 on the cathode boundary (at coordinate r = ri) and to 0 at the outer electrode
(at coordinate r = ro). The mean electron energy of the secondary electron is set to 4 eV.

When solving the problem, you can choose to plot the electron and ion density as the
solver runs. Doing so allows the behavior of the charged species to be observed as the
model is solving.
Looking at the evolution of the densities during the solving process shows that the gas is
initially weekly ionized (electrons and ions having relatively low densities compared to
neutral atoms). As the negative voltage is applied to the cathode, the highly mobile
electrons are accelerated towards the anode leaving behind a positively charged gas in the
cathode vicinity. With increasing negative potential, ion bombardment becomes more
persistent on the cathode surface, which generates more secondary electrons, ionizing
more neutral atoms, and engendering greater ion current. As the negative potential rises,
the population of charged particles becomes larger as a consequence of this avalanche
effect.
As the ion current becomes more significant, the RC circuit reduces the cathode negative
potential such that an equilibrium is reached between the generation of the charged
particles, preventing the transition of the plasma into an arcing regime.
Figure 2 shows the electron and ion densities at the end of the simulation (t = 0.1 s). Note
that the figure is plotted on a log-log scale. In the figure, one can see that the ion density
is approximately three orders of magnitude higher than the electron density in the vicinity
of the cathode. The tiny positive space charge distribution generated by the ion density
defines an ionization region that screens the cathode potential from the anode. This can
be observed by displaying the electric potential along the radius of the coaxial assembly;
see Figure 3.

Figure 2: Plot of electron (blue) and ion (green) number density at the end of the simulation
(t = 0.1 s).

Figure 3: Electric potential along the radius of the coaxial assembly at the end of the
simulation (t = 0.1 s).
The high voltage and important current density at the cathode enhance the electron
temperature near the electrode, thus boosting ionization of the neutral atoms in the
ionization zone. This can be seen in Figure 4, which displays the electron temperature at
the end of the simulation.

Figure 4: Electron temperature along the radius of the assembly at the end of the simulation
(t = 0.1 s).
To see the effect of the RC circuit on the system, plot the secondary electron flux as a
function of time at the cathode surface; see Figure 5. A close look at the figure shows a
rapid rise of the flux that tends to stabilize as the circuit sees higher currents. Figure 6 also
shows the potential at both electrodes as a function of time. Comparing figure Figure 5
with Figure 6 reveals the effect of the circuit on both potential and secondary emission at
the cathode surface. One can also observe a similar effect by plotting the electron current
density at the electrodes, Figure 6. Doing so shows a direct relation between the cathode
secondary emission (positive current flow), the anode recombination (negative flow), and
the RC circuit. Figure 7 displays a surface plot of the logarithm of the electron density at
the last time step.

Figure 5: Secondary emission flux along time at the cathode.

Figure 6: Potential at the electrodes as a function of time.

Figure 7: Electron current density at the electrodes as a function of time.

Figure 8: 2D electron density (log of ne) at the last time step.

corona_discharge_1d
NEW
MODEL WIZARD
3 Click Add.

4 Click Study.
6 Click Done.
To set up this model, you need access to some advanced settings.
PLASMA (PLAS)
Solve this model using the finite volume method and a Scharfetter-Gummel scheme. You
can switch from the finite element method to the finite volume method by suitable choices
in the Discretization section.
1 In the Model Builder windows toolbar, click the Show button and select Discretization in
the menu.
3 In the Settings window for Plasma, click to expand the Discretization section.
4 From the Formulation list, choose Finite volume (constant shape function).
GLOBAL DEFINITIONS
Parameters

V0 1000[V] 1000 V Voltage
GEOMETRY 1
Follow the steps below to create the model geometry: a simple 1D geometry consisting of
a single domain bounded by the cathode (left, inner conductor) and the anode (right,
outer conductor).
1 In the Model Builder window, under Component 1 (comp1) click Geometry 1.

2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose cm.
Interval 1 (i1)

4 In the Right endpoint text field, type 5.
PLASMA (PLAS)
Load the argon cross sections from file. They form the basis of the plasma chemistry under
investigation. To fix the species data, load the argon preset for each species.

3 Click Browse.
Ar_xsecs.txt.
When solving any type of reacting flow problem one species must always be chosen to
fulfill the mass constraint. This should be taken as the species with the largest mass
fraction.
Species: Ar
When solving a plasma problem the plasma must initially be charge neutral. COMSOL
initial electroneutrality constraint is satisfied. Now let the initial number density of
argon ions be the same as the initial number of electrons. This forces the plasma to be
initially charge neutral.
Species: Ar+
Ar+.
Species: Ar

Species: Ars
Ars.
Now add two more regular reactions that 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
3 In the Formula text field, type Ar+Ars=>Ar+Ar.
Reaction 2
3 In the Formula text field, type Ars+Ars=>Ars+Ar.
4 Locate the Reaction Parameters section. In the kf text field, type 2.3e7.
Surface reactions must always be included in a plasma model because they describe how
ionic, excited, and radical species interact with the wall.
Surface Reaction 1
4 Locate the Boundary Selection section. From the Selection list, choose All boundaries.
Secondary emission of electrons is important to sustain a DC discharge. In this example,
add a secondary emission coefficient on the left wall (cathode).
Surface Reaction 2

6 In the i text field, type 4.
Surface Reaction 3
8 Locate the Out-of-Plane Thickness section. In the dz text field, type 0.02[m].
Plasma Model 1
Model 1.
3 In the eNn text field, type 4e24.
Ground 1
Now, add an RC circuit in series with the system. The role of the circuit is to limit the
current at the cathode and avoid the arcing regime.
Metal Contact 1
1 On the Physics toolbar, click Boundaries and choose Metal Contact.
3 In the Settings window for Metal Contact, locate the Terminal section.
4 In the V0 text field, type -V0.
5 Locate the Quick Circuit Settings section. From the Quick circuit type list, choose Series
RC circuit.

Wall 1
2 In the Settings window for Wall, locate the Boundary Selection section.
3 From the Selection list, choose All boundaries.
MESH 1
Meshing is a critical step in any plasma model. A fine mesh is needed close to the electrodes
to capture the separation of space charge between the electrons and ions close to the wall.
1 In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
Edge.
Edge 1
Distribution 1
7 Click Build All.
STUDY 1
Compute an initial solution to generate the default plots and then set them up to show the
electron and ion density while the solver runs.
1 In the Model Builder window, right-click Study 1 and choose Get Initial Value for Step.
RESULTS

Set up an electron and ion density plot.

4 In the associated text field, type Density (1/m3).
Line Graph 1
1 In the Model Builder window, expand the Electron Density (plas) node, then click Line
Graph 1.
3 Select the Show legends check box.
6 Locate the Legends section. In the table, enter the following settings:
Legends
Electrons
7 Right-click Results>Electron Density (plas)>Line Graph 1 and choose Rename.

8 In the Rename Line Graph dialog box, type Electrons in the New label text field.
9 Click OK.
Electrons 1
1 Right-click Results>Electron Density (plas)>Electrons and choose Duplicate.
2 In the Model Builder window, under Results>Electron Density (plas) right-click Electrons
1 and choose Rename.
3 In the Rename Line Graph dialog box, type Ions in the New label text field.
4 Click OK.
6 In the Expression text field, type plas.n_wAr_1p.
Legends
Ions

3 Right-click Results>Electron Density (plas) and choose Rename.

4 In the Rename 1D Plot Group dialog box, type Electron and Ion Number Density in
the New label text field.
5 Click OK.
6 Click the y-Axis Log Scale button on the Graphics toolbar.
Now set up the last two default plots.

Line Graph 1
1 In the Model Builder window, expand the Electron Temperature (plas) node, then click
Line Graph 1.
2 In the Settings window for Line Graph, click to expand the Title section.
3 From the Title type list, choose None.

1 In the Model Builder window, under Results click Electric Potential (plas).
Line Graph 1
1 In the Model Builder window, expand the Electric Potential (plas) node, then click Line
Graph 1.
2 In the Settings window for Line Graph, locate the Title section.
3 From the Title type list, choose None.
4 Locate the y-Axis Data section. Select the Description check box.
5 In the associated text field, type Electron potential.
STUDY 1

1 In the Model Builder window, under Study 1 click Step 1: Time Dependent.
3 In the Times text field, type 0 .
4 Click Range.

6 In the Start text field, type -8.
7 In the Stop text field, type 0.
9 From the Function to apply to all values list, choose exp10.
10 Click Add.
11 In the Settings window for Time Dependent, click to expand the Results while solving
section.
12 Locate the Results While Solving section. Select the Plot check box.
RESULTS

Now plot the electric potential on the electrodes.
1D Plot Group 4
On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
Point Graph 1
1 On the 1D Plot Group 4 toolbar, click Point Graph.
2 In the Settings window for Point Graph, locate the Selection section.
4 Locate the y-Axis Data section. In the Expression text field, type V.
Legends
Cathode
Anode
8 Locate the y-Axis Data section. Select the Description check box.
9 In the associated text field, type Electric Potential.

1D Plot Group 4
3 From the Position list, choose Lower right.
4 Right-click Results>1D Plot Group 4 and choose Rename.
5 In the Rename 1D Plot Group dialog box, type Electric Potential on Electrodes
in the New label text field.
6 Click OK.
8 On the Electric Potential on Electrodes toolbar, click Plot.
Electric Potential on Electrodes

1 In the Model Builder window, collapse the Results>Electric Potential on Electrodes node.
Plot the electron current density on the electrodes as follows:
1D Plot Group 5
Point Graph 1
2 In the Settings window for Point Graph, locate the Selection section.
4 Locate the y-Axis Data section. In the Expression text field, type plas.nJe.
5 Locate the Legends section. Select the Show legends check box.
Legends
cathode
anode
8 Locate the Title section. From the Title type list, choose None.
1D Plot Group 5

4 Right-click Results>1D Plot Group 5 and choose Rename.
5 In the Rename 1D Plot Group dialog box, type Electron Current Density on
Electrodes in the New label text field.
6 Click OK.
8 On the Electron Current Density on Electrodes toolbar, click Plot.
Electron Current Density on Electrodes

1 In the Model Builder window, collapse the Results>Electron Current Density on Electrodes
node.
Next, plot the secondary electron emission at the cathode.
1D Plot Group 6
Point Graph 1
3 In the Settings window for Point Graph, locate the y-Axis Data section.
4 In the Expression text field, type plas.sflux.
5 Locate the Title section. From the Title type list, choose None.
1D Plot Group 6
1 In the Model Builder window, under Results right-click 1D Plot Group 6 and choose
Rename.
2 In the Rename 1D Plot Group dialog box, type Secondary Emission Flux in the New
label text field.
3 Click OK.
Secondary Emission Flux

1 In the Model Builder window, collapse the Results>Secondary Emission Flux node.
3 On the Secondary Emission Flux toolbar, click Plot.
Finally, create a nice plot for the model thumbnail.
Revolution 1D 1
On the Results toolbar, click More Data Sets and choose Revolution 1D.

Data Sets
In the Model Builder window, collapse the Results>Data Sets node.
RESULTS
2D Plot Group 7
On the Results toolbar, click 2D Plot Group.
Surface 1
1 In the Model Builder window, right-click 2D Plot Group 7 and choose Surface.
3 In the Expression text field, type log(plas.ne).
2D Plot Group 7
1 In the Model Builder window, under Results right-click 2D Plot Group 7 and choose
Rename.
2 In the Rename 2D Plot Group dialog box, type Log of Electron Density in the New
label text field.
3 Click OK.
Log of Electron Density

1 In the Model Builder window, collapse the Results>Log of Electron Density node.
2 On the Log of Electron Density toolbar, click Plot.

Dipolar Microwave Plasma Source
Introduction
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
unique, including:
Infinite elements for the magnetostatic model.

Functional-based mesh adaption 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.
Note: This application requires the Plasma Module, AC/DC Module, and RF Module.
Model DefinitionStatic Magnetic Field

For the static magnetic field, Ampres law governs the azimuthal component of the
magnetic vector potential:
1 1
r 0 ( A ) = J
where the external current density, J only has an azimuthal component and is defined in
the coil as:
NI
J = -------
A
2 | DIPOLAR MICROWAVE PLASMA SOURCE

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.
Antenna
Coil
Plasma
Infinite elements
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
charge:
m e
B cr = ------------
q
At 2.45 GHz the critical magnetic flux density is 875 gauss or 0.0875 T. Therefore you
can use functional-based mesh adaption 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
f = ---------------------------------------------- (1)
B 0.0875 +
where is a small number to prevent division by zero.
Model DefinitionMicrowave Plasma

In this example, you solve the following wave equation for the high-frequency component
of the electric field in the frequency domain:
0 ( E ) k 0 r --------- E = 0
1 2 j
0
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
q
= ------------------------------- (2)
m e ( e + j )
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
1 B z B y
1
qne = B z 1 B x (3)
B y B x 1
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 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
COMSOL by selecting the Doppler broadening 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:

e = e + ---- (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:
n
( ) + [ n e ( e E ) D e n e ] = R e
t e
( n ) + [ n ( E ) D n ] + E = R
t e

The electron diffusivity, energy mobility and energy diffusivity are calculated from the
electron mobility using
D e = e T e, = --- e, D = T e
5
3
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
1
-------- Bz By
dc
1 1
e = B z -------- B x (5)
dc
1
B y B x --------
dc
where dc is the electron mobility in the absence of a magnetic field. COMSOL

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 which 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
Re = x j kj Nn ne
j=1
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:
R = xj kj Nn ne j
j=1
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 may
be computed from cross section data by the following integral:

kk = 0 k ( )f ( ) d
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:
3 2
f() = 1 exp ( ( 2 ) )
where is the mean electron energy:
n
= -----
ne
and
32 5 2 1
1 = ( 5 2 ) (3 2) , 2 = ( 5 2 ) ( 3 2 )
where is the gamma function. The heating term, Ee 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:
10 W
= ----------------------------- (6)
Q ind dV
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 due to the fact that 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:

( w k ) + ( u )w k = j k + R k (7)
t
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 look up table. The ion diffusion velocity, vk, is
related to the diffusive flux via
j k = v k
where
v k = D m ln ( w ) + D m ln ( M ) + Z ( E + v k B ) (8)
Equation 7 can be re-arranged to give an expression for the diffusion velocity as a function
of the other variables:
1
v k = A [ D m ln ( w ) + D m ln ( M ) + ZE ]
where
1 ZB z ZB y
A = ZB z 1 ZB x
ZB y ZB x 1
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:
0 r V =
specified in the model using the formula
N
= q Z k n k n e
k = 1

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.

On surfaces, the following two reactions are considered:

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:
n e = --- e, th n e
1
2
n = --- e, th n
5
6

Losses at the wall for the heavy species is due to surface reactions and migration due to the
ambipolar field:
1 1
n j k = M w R k + M w c k Z k [ ( A E ) n ] [ ( Zk k ( A E) n) > 0]
The reactor walls are grounded.
Notes About the COMSOL Implementation

You solve this problem in two stages. First, compute the static magnetic field using
adaptive mesh refinement. Then, in a separate study step, you 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.
Results and DiscussionStatic 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.

ECR surface
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 adaption 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 adaption. The
mesh is very fine on the ECR surface and relatively coarse away from the resonance zone.
Results and DiscussionMicrowave Plasma Model

The electron density at the quasi steady state solution is plotted in Figure 5. The peak
electron density is around 51016 m3 and peaks radially outwards 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 tensors 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
4104 m2/(Vs) towards the coil edges and 104 m2/(Vs) radially outwards 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 towards the coil edges. This leads to
migration of electrons along the magnetic field lines when they are produced in the
ionization region, Figure 9.

Figure 11: Plot of the electron mobility tensors 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.41016 m3 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.
Application Library path: Plasma_Module/Wave-Heated_Discharges/

dipolar_ecr_source

NEW
MODEL WIZARD
2 In the Select Physics tree, select AC/DC>Magnetic Fields (mf).
3 Click Add.
4 Click Study.
5 In the Select Study tree, select Preset Studies>Stationary.
6 Click Done.
GLOBAL DEFINITIONS
Parameters

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)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type r0.
4 In the Height text field, type z0.
5 Locate the Position section. In the z text field, type -z0/2.
Rectangle 2 (r2)

3 In the Width text field, type 0.01.
4 In the Height text field, type 0.03.
5 Locate the Position section. In the z text field, type 0.04.
Rectangle 3 (r3)
5 Locate the Position section. In the r text field, type 0.006.
6 In the z text field, type 0.072.
Rectangle 4 (r4)
Rectangle 5 (r5)
3 In the Width text field, type r0-0.01.
Rectangle 6 (r6)
5 Locate the Position section. In the r text field, type r0.

Rectangle 7 (r7)
4 In the Height text field, type z0.
6 In the z text field, type -z0/2.
Rectangle 8 (r8)
6 In the z text field, type -0.02-z0/2.
Rectangle 9 (r9)
3 In the Width text field, type r0.
5 Locate the Position section. In the z text field, type -0.02-z0/2.
Rectangle 10 (r10)
Bzier Polygon 1 (b1)

1 On the Geometry toolbar, click Primitives and choose Bzier Polygon.
2 In the Settings window for Bzier Polygon, locate the Polygon Segments section.
3 Find the Added segments subsection. Click Add Linear.
4 Find the Control points subsection. In row 1, set r to 0.01 and z to 0.072.

5 In row 2, set r to 0.01 and z to 0.07.
DEFINITIONS
1 In the Model Builder window, expand the Component 1 (comp1)>Definitions node.
2 Right-click Definitions and choose View.
Axis
1 In the Model Builder window, expand the View 2 node, then click Axis.
2 In the Settings window for Axis, locate the Axis section.
3 In the r minimum text field, type -0.05.
4 In the r maximum text field, type 0.16.
5 In the z minimum text field, type -0.02.
6 In the z maximum text field, type 0.12.
7 Click Update.
View 2
1 In the Model Builder window, under Component 1 (comp1)>Definitions click View 2.
2 In the Settings window for View, locate the View section.
3 Select the Lock axis check box.
Explicit 1
5 Right-click Explicit 1 and choose Rename.
7 Click OK.
Infinite Element Domain 1 (ie1)

1 On the Definitions toolbar, click Infinite Element Domain.
2 In the Settings window for Infinite Element Domain, locate the Geometry section.
3 From the Type list, choose Cylindrical.
4 Select Domains 1, 5, and 811 only.

MATERIALS
Blank Material.
Material 1 (mat1)
1 Select Domains 1, 2, 4, 5, and 711 only.

group
Material 2 (mat2)

group
Material 3 (mat3)

group

DEFINITIONS
Integration 1a (intop1)
Coil 1
3 In the Settings window for Coil, locate the Coil section.
4 From the Conductor model list, choose Homogenized multi-turn.
5 Locate the Homogenized Multi-Turn Conductor section. In the N text field, type 5000.
6 Locate the Coil section. In the Icoil text field, type 14.
MESH 1
More Operations>Edge.
Edge 1
Select Boundary 19 only.
Size 1
4 Click the Custom button.
6 In the associated text field, type 0.0005.
Edge 2
1 In the Model Builder window, right-click Mesh 1 and choose More Operations>Edge.
2 Select Boundaries 6, 18, 20, and 22 only.
Size 1

Free Triangular 1
Size 1
Free Triangular 2
Size 1
Free Triangular 3
In the Settings window for Free Triangular, click Build All.

STUDY 1
Step 1: Stationary
1 In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Adaptive mesh refinement check box.
Solution 1 (sol1)
3 In the Model Builder window, expand the Solution 1 (sol1) node.
4 In the Model Builder window, expand the Study 1>Solver Configurations>Solution 1
(sol1)>Stationary Solver 1 node, then click Adaptive Mesh Refinement.
5 In the Settings window for Adaptive Mesh Refinement, locate the General section.
6 In the Maximum number of refinements text field, type 1.
7 Locate the Error Estimation section. From the Error estimate list, choose Functional.
8 In the Functional text field, type intop1(1/(abs(mf.normB-0.0875)+1e-4)).
9 Locate the Mesh Refinement section. From the Refinement method list, choose Mesh
initialization.
RESULTS
Reproduce the magnetic flux density plot in Figure 3 with the following steps.
1 In the Model Builder window, expand the Results node.

In the Model Builder window, expand the Results>Data Sets node, then click Study 1/Solution
1 (sol1).
Selection

4 Select Domains 2, 4, 6, and 7 only.
2D Plot Group 1
2 In the Settings window for 2D Plot Group, type Stationary Magnetic Flux
Density in the Label text field.
Contour 1
1 Right-click Stationary Magnetic Flux Density and choose Contour.
2 In the Settings window for Contour, locate the Expression section.
3 In the Expression text field, type log10(mf.normB+eps).
4 Locate the Levels section. In the Total levels text field, type 8.
5 Locate the Coloring and Style section. From the Contour type list, choose Filled.
6 From the Color table list, choose GrayScale.
7 Clear the Color legend check box.
8 Select the Reverse color table check box.
Stationary Magnetic Flux Density

In the Model Builder window, under Results right-click Stationary Magnetic Flux Density and
choose Contour.
Contour 2
1 In the Settings window for Contour, locate the Expression section.
2 In the Expression text field, type Aphi.
3 Locate the Coloring and Style section. From the Coloring list, choose Uniform.
4 From the Color list, choose Black.
5 Clear the Color legend check box.
Contour 3
1 Right-click Results>Stationary Magnetic Flux Density>Contour 2 and choose Duplicate.
2 In the Settings window for Contour, click Replace Expression in the upper-right corner
of the Expression section. From the menu, choose Component 1>Magnetic Fields>
Magnetic>mf.normB - Magnetic flux density norm.
3 Locate the Levels section. From the Entry method list, choose Levels.
4 Click Range.
6 In the Start text field, type 0.086.

7 In the Stop text field, type 0.089.
9 Click Replace.
10 On the Stationary Magnetic Flux Density toolbar, click Plot.
Next, visualize the refined mesh.
2D Plot Group 2
2 In the Settings window for 2D Plot Group, type Adpative Mesh in the Label text field.
3 Locate the Data section. From the Data set list, choose Study 1/Adaptive Mesh Refinement
1 (sol2).
Mesh 1
1 Right-click Adpative Mesh and choose Mesh.
2 In the Settings window for Mesh, locate the Color section.
3 From the Element color list, choose None.
4 From the Wireframe color list, choose Custom.
5 On Windows, click the colored bar underneath, or?if you are running the cross-platform
desktop?the Color button.
6 Click Define custom colors.
7 Set the RGB values to 192, 192, and 192, respectively.
8 Click Add to custom colors.
9 Click Show color palette only or OK on the cross-platform desktop.
10 On the Adpative Mesh toolbar, click Plot.
Compare the result with the plot in Figure 4.
ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the tree, select Plasma>Microwave Plasma.
4 Click Add to Component in the window toolbar.
PLASMA (PLAS)
On the Physics toolbar, click Magnetic Fields (mf) and choose Plasma (plas).
1 On the Home toolbar, click Add Physics to close the Add Physics window.

3 Select Domains 1, 2, 5, 6, and 811 only.
E LE C TRO MA G NE TI C WAVES , FRE QU E N C Y D OM A IN ( EM W )

On the Physics toolbar, click Plasma (plas) and choose Electromagnetic Waves, Frequency
Domain (emw).
1 In the Model Builder window, under Component 1 (comp1) click Electromagnetic Waves,
Frequency Domain (emw).
2 In the Settings window for Electromagnetic Waves, Frequency Domain, locate the
Domain Selection section.
PLASMA (PLAS)
On the Physics toolbar, click Electromagnetic Waves, Frequency Domain (emw) and choose
Plasma (plas).
3 Select the Full expression for diffusivity check box.
4 Select the Compute tensor ion transport properties check box.
5 Locate the Plasma Properties section. Select the Compute tensor electron transport
properties check box.
MULTIPHYSICS
1 In the Model Builder window, under Component 1 (comp1)>Multiphysics click Plasma
Conductivity Coupling 1 (pcc1).
2 In the Settings window for Plasma Conductivity Coupling, locate the Compute Tensor
Plasma Conductivity section.
3 Select the Compute tensor plasma conductivity check box.
PLASMA (PLAS)
2 In the Settings window for Plasma, locate the Domain Selection section.


3 Click Browse.
Ar_xsecs.txt.
Reaction 1
Reaction 2
Species: Ar
Species: Ars
Ars.
Species: Ar+
Ar+.

5 Click to expand the Mobility and diffusivity expressions section. Locate the Mobility and
Diffusivity Expressions section. From the Specification list, choose Specify mobility,
compute diffusivity.
6 Click to expand the Mobility specification section. Locate the Mobility Specification
section. From the Specify using list, choose Lookup table.
7 Click Load from File.
ion_mobility_data.txt.
Surface Reaction 1
Surface Reaction 2
Wall 1
Ground 1
PLASMA (PLAS)
On the Physics toolbar, click Electromagnetic Waves, Frequency Domain (emw) and choose
Plasma (plas).

Initial Values 1
Values 1.
Plasma Model 1
Model 1.
3 From the B list, choose Magnetic flux density (mf).
4 In the T text field, type 300.
5 In the pA text field, type 1.
6 Locate the Electron Density and Energy section. In the e text field, type 1E25[1/(m*V*
s)]/plas.Nn.

On the Physics toolbar, click Plasma (plas) and choose Electromagnetic Waves, Frequency
Domain (emw).
In the Model Builder window, under Component 1 (comp1) click Electromagnetic Waves,
Port 1
1 On the Physics toolbar, click Boundaries and choose Port.
3 In the Settings window for Port, locate the Port Properties section.
4 From the Type of port list, choose Coaxial.
5 From the Wave excitation at this port list, choose On.
6 Select the Specify deposited power check box.
7 In the Pdep text field, type Psp.
COMPONENT 1 (COMP1)
Mesh 3
On the Mesh toolbar, click Add Mesh.

MESH 3
In the Model Builder window, under Component 1 (comp1)>Meshes right-click Mesh 3 and
choose More Operations>Reference.
Reference 1
1 In the Settings window for Reference, locate the Reference section.
2 From the Mesh list, choose Mesh 2.
Refine 1
1 In the Model Builder window, right-click Mesh 3 and choose More Operations>Refine.
2 In the Settings window for Refine, locate the Domain Selection section.
5 Locate the Refine Options section. In the Number of refinements text field, type 2.
Boundary Layers 1
1 Right-click Mesh 3 and choose Boundary Layers.

1 In the Model Builder window, under Component 1 (comp1)>Meshes>Mesh 3>Boundary
Layers 1 click Boundary Layer Properties.
section.
type 4.
6 Click Build All.
ADD STUDY
3 Find the Studies subsection. In the Select Study tree, select Preset Studies for Multiphysics.

4 Find the Physics interfaces in study subsection. In the table, clear the Solve check box for
the Magnetic Fields (mf) interface.
5 Find the Studies subsection. In the Select Study tree, select Preset Studies for
Multiphysics>Frequency-Transient.
STUDY 2
3 In the Times text field, type 0 10^{range(-8,6/10,-2)}.
5 Locate the Physics and Variables Selection section. In the table, clear the Solve for check
box for the Magnetic Fields (mf) interface.
6 Click to expand the Mesh selection section. Locate the Mesh Selection section. In the
table, enter the following settings:
Geometry Mesh
Geometry 1 mesh3
7 Click to expand the Mesh selection section. Click to expand the Values of dependent
variables section. In the Settings window for Frequency-Transient, locate the Values of
Dependent Variables section.
8 Find the Values of variables not solved for subsection. From the Settings list, choose User
controlled.
9 From the Method list, choose Solution.
10 From the Study list, choose Study 1, Stationary.
11 From the Solution list, choose Adaptive Mesh Refinement 1 (sol2).
Solution 3 (sol3)
3 In the Model Builder window, expand the Solution 3 (sol3) node, then click Dependent
Variables 1.
4 In the Settings window for Dependent Variables, locate the General section.

5 From the Defined by study step list, choose User defined.
6 Locate the Values of Variables Not Solved For section. From the Method list, choose
Solution.
7 From the Solution list, choose Adaptive Mesh Refinement 1 (sol2).
RESULTS
Three new default plots show the electron density, electron temperature, and plasma
potential respectively; compare with those in Figure 5, Figure 6, and Figure 7.

Follow the instructions below to reproduce Figure 8 through Figure 15.
2D Plot Group 7
field.
3 Locate the Data section. From the Data set list, choose Study 2/Solution 3 (sol3).
Surface 1
the Expression section. From the menu, choose Component 1>Electromagnetic Waves,
Frequency Domain>Heating and losses>emw.Qrh - Resistive losses.
Compare with Figure 8.
Power Deposition
In the Model Builder window, under Results right-click Power Deposition and choose
Duplicate.
Power Deposition 1
In the Settings window for 2D Plot Group, type Electron Source in the Label text field.
Surface 1
1 In the Model Builder window, expand the Results>Electron Source node, then click Surface
1.
3 In the Expression text field, type plas.Re.

4 On the Electron Source toolbar, click Plot.
Electron Source
In the Model Builder window, under Results right-click Electron Source and choose
Duplicate.
Electron Source 1
In the Settings window for 2D Plot Group, type Electron Mobility, rr Component in
Surface 1
1 In the Model Builder window, expand the Results>Electron Mobility, rr Component node,
Electron transport properties>Electron mobility>plas.muerr - Electron mobility, rr
component.
3 On the Electron Mobility, rr Component toolbar, click Plot.
Electron Mobility, rr Component

In the Model Builder window, under Results right-click Electron Mobility, rr Component and
choose Duplicate.
Electron Mobility, rr Component 1

In the Settings window for 2D Plot Group, type Electron Mobility, zz Component in
Surface 1
1 In the Model Builder window, expand the Results>Electron Mobility, zz Component node,
Electron transport properties>Electron mobility>plas.muezz - Electron mobility, zz
component.
3 On the Electron Mobility, zz Component toolbar, click Plot.

2D Plot Group 11
2 In the Settings window for 2D Plot Group, type Conduction Current Density,
r-Component in the Label text field.
3 Locate the Data section. From the Data set list, choose Study 2/Solution 3 (sol3).
Surface 1
1 Right-click Conduction Current Density, r-Component and choose Surface.
3 In the Expression text field, type emw.Jr-emw.Jdr.
4 On the Conduction Current Density, r-Component toolbar, click Plot.
Conduction Current Density, r-Component

In the Model Builder window, under Results right-click Conduction Current Density,
r-Component and choose Duplicate.
Conduction Current Density, r-Component 1

In the Settings window for 2D Plot Group, type Conduction Current Density,
z-Component in the Label text field.
Surface 1
1 In the Model Builder window, expand the Results>Conduction Current Density,
z-Component node, then click Surface 1.
3 In the Expression text field, type emw.Jz-emw.Jdz.
4 On the Conduction Current Density, z-Component toolbar, click Plot.
Conduction Current Density, z-Component

In the Model Builder window, under Results right-click Conduction Current Density,
z-Component and choose Duplicate.
Conduction Current Density, z-Component 1

In the Settings window for 2D Plot Group, type Conduction Current Density,
phi-Component in the Label text field.

Surface 1
1 In the Model Builder window, expand the Results>Conduction Current Density,
phi-Component node, then click Surface 1.
3 In the Expression text field, type emw.Jphi-emw.Jdphi.
4 On the Conduction Current Density, phi-Component toolbar, click Plot.
Electron Mobility, rr Component

In the Model Builder window, under Results right-click Electron Mobility, rr Component and
choose Duplicate.
Electron Mobility, rr Component 1

In the Settings window for 2D Plot Group, type Mean Plasma Electrical
Conductivity in the Label text field.
Surface 1
1 In the Model Builder window, expand the Results>Mean Plasma Electrical Conductivity
node, then click Surface 1.
3 In the Expression text field, type (emw.sigmarr+emw.sigmaphiphi+emw.sigmazz)/3.
4 On the Mean Plasma Electrical Conductivity toolbar, click Plot.

Drift Diffusion Tutorial
Introduction
The foundation of the COMSOL Multiphysics Plasma Module is the Drift Diffusion
interface which describes the transport of electrons in an electric field. The Drift Diffusion
interface solves a pair of reaction/convection/diffusion equations, one for the electron
density and the other for the electron energy density. The mean electron energy is
computed by dividing the electron energy density by the electron number density.
Model Definition
This tutorial example computes the electron number density and mean electron energy in
a drift tube. Electrons are released due to thermionic emission on the left boundary with
an assumed mean electron energy. The electrons are then accelerated towards the right
boundary due to an imposed external electric field which is oriented in the opposite
direction from the electron drift velocity:
Electric field
Electron Wall
influx
Figure 1: In the drift tube the electrons enter the left boundary and are accelerated by the
electric field toward the wall.
MODEL EQUATIONS
A simple model is set up to test the Drift Diffusion interface. The equations solved are, for
the electron number density:
n
( ) + e = Re (1)
t e
where
e = ( e E )n e D e n e (2)
and, ne denotes the electron density (1/m3), Re is the electron rate expression (1/(m3s)),
e is the electron mobility which is either a scalar or tensor (m2/(Vs)), E is the electric
field (V/m), and De is the electron diffusivity, which is either a scalar or a tensor. The first
term on the right side of Equation 2 represents migration of electrons due to an electric
field. The seconds term on the right side of Equation 2 represents diffusion of electrons
2 | D R I F T D I F F U S I O N TU T O R I A L
from regions of high electron density to low electron density.
An equation for the electron energy density is solved in conjunction with Equation 1. The
equation is:
(n ) + + E = R
(3)
t e
where
= ( E )n D n
Here, n is the electron energy density (V/m3), R is the energy loss/gain due to inelastic
collisions (V/(m3s)), is the electron energy mobility (m2/(Vs)), E is the electric field
(V/m), and D is the electron energy diffusivity (m2/s). The subscript refers to electron
energy. The third term on the left side of Equation 3 represents heating of the electrons
due to an external electric field. Note that this term can either be a heat source or a heat
sink depending on whether the electrons are drifting in the same direction as the electric
field or not. For a Maxwellian electron energy distribution function, the following
relationships hold:
D e = e T e, = --- e, D = T e
5
3
where Te is the electron temperature. So, given the electron mobility, the other
transport properties required can be computed. The electron temperature is a function
of the mean electron energy, which is defined as:
n
= -----
ne
and then:
T e = ---
2
3
SOURCE COEFFICIENTS
The electron source term, Re is the sum of the electron impact reaction rates that make up
the plasma chemistry. The electron energy loss due to inelastic collisions, R is a function
of the electron impact reaction rates multiplied by the energy loss corresponding to each
reaction. Mathematically, the electron source is defined as:
M
Re = x j kj Nn ne
j=1
reaction j (m3/s), and Nn is the total neutral number density (1/m3). The energy loss is
defined as:
R = xj kj Nn ne j
j=1
where j is the energy loss from reaction j (V).
An influx of electron due to thermionic emission is specified on the left boundary. The
electrons that are emitted from the surface are then accelerated toward the wall by the
electric field. The acceleration leads to an increase in the mean electron energy and
subsequently ionization occurs. This creates new electrons, which are ultimately lost into
a wall opposite the emitting surface. The gain of electrons due to ionization is included in
the model by assuming the drift tube contains argon. This example assumes that the mole
fraction of electronically excited argon atoms and argon ions is very small. This means that
you do not solve additional equations for the mole fractions of these two species. The
following reactions are included in the model:
REACTION FORMULA TYPE (eV) A B E
1 e+Ar=>e+Ar Elastic 0 1.99E-014 0.93 0.41

2 e+Ar=>e+Ars Excitation 11.5 8.77E-015 0.62 18.16
3 e+Ar=>2e+Ar+ Ionization 15.8 2.15E-014 0.49 24.75
The rate constants for these reactions are given in Arrhenius form:
B
k j = ATe exp ( E T e )
The rate constants are computed from cross-section data for each reaction assuming a
Maxwellian electron energy distribution function.
The electron density is plotted in Figure 2. The peak electron density occurs close to the
far wall. The peak electron density is five times higher than the electron density on the left
wall due to new electrons being created through ionization.
Figure 2: Plot of the electron density in the drift tube.
Because there are no variations in the solution in the y-direction it can be more convenient
to create a 1D data set in the x-direction and plot various quantities v.s. the x-axis. Figure 3
plots the electron density as a function of the x-direction. The electron temperature is
plotted in Figure 4. On the left wall the electron temperature is fixed to 2 eV. The
temperature steadily increases over a narrow region. This is due to the fact that there is a
strong drift velocity in the opposite direction to the electric field. As the electron
temperature increases, so do the rate constants which are responsible for creating new
electrons. The increase in electron temperature also has a significant on the number of
inelastic collisions which occur in the tube. After the initial rise in electron temperature,
the electron temperature remains constant until close to the far wall. In this region the
joule heating caused by the electron drift velocity in the opposite direction to the electric
field is balanced by the energy loss due to inelastic collisions.
The highly nonlinear behavior in such a simple example showcases the fact that very
complicated dynamics occur in even the simplest of plasmas.
Figure 3: Plot of the electron temperature across the drift tube.
Figure 4: Cross section plot of electron density across the drift tube.
Figure 5: Cross plot of the electron temperature across the drift tube.
Reference
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. 722733, 2005.

drift_diffusion_tutorial
NEW
MODEL WIZARD
2 In the Select Physics tree, select Plasma>Drift Diffusion (dd).
3 Click Add.
4 Click Study.
6 Click Done.
GEOMETRY 1
Rectangle 1 (r1)
3 In the Width text field, type 5E-3.
4 In the Height text field, type 5E-4.
Now add a table of expressions for the external electric field, temperature, pressure, and
the electron impact reactions occurring in the drift tube. You also define the influx of
electrons due to thermionic emission. The influx of electrons and the external electric
field are both ramped up from an initial value of zero to aid convergence at early
timesteps.
DEFINITIONS
Variables 1

ramp tanh(1E7[1/s]*t) Ramp function
T 300[K] K Gas temperature
p 10[torr] Pa Gas pressure
Nn0 p/k_B_const/T 1/m Neutral number
density
k1 1.993039e-014* Elastic rate
dd.Te^0.93* coefficient
exp(-0.41/dd.Te)
k2 8.773932e-015* Excitation rate
exp(-18.16/dd.Te)
k3 2.153106e-014* Ionization rate
exp(-24.75/dd.Te)
r1 k1*dd.Nn*dd.ne Elastic collision
reaction rate
r2 k2*dd.Nn*dd.ne Electronic excitation
reaction rate
r3 k3*dd.Nn*dd.ne Ionization reaction
rate
de1 0[V] V Energy loss, elastic
collision
de2 11.56[V] V Energy loss,
electronic excitation
de3 15.8[V] V Energy loss,
ionization
Re r3 Electron production
rate
Sen -e_const*(r1*de1+r2* Collisional power
de2+r3*de3) loss
Ex -d(V,x) V/m Electric field
V -100[V]*ramp* V Electric potential
(0.005[m]-x)/
0.005[m]
Ain 0.5E-3^2[m^2] m Wall area
I 2E-6[A] A Thermal emission
current
influx I/Ain/e_const*ramp 1/(ms) Electron influx
mobility
DRIFT DIFFUSION (DD)

1 In the Model Builder window, under Component 1 (comp1) click Drift Diffusion (dd).
2 In the Settings window for Drift Diffusion, locate the Electron Properties section.
Drift Diffusion Model 1

1 In the Model Builder window, under Component 1 (comp1)>Drift Diffusion (dd) click Drift
Diffusion Model 1.
2 In the Settings window for Drift Diffusion Model, locate the Model Inputs section.
3 In the Nn text field, type Nn0.
4 In the V text field, type V.
5 In the Sen text field, type Sen.
Initial Values 1
1 In the Model Builder window, under Component 1 (comp1)>Drift Diffusion (dd) click Initial
Values 1.
Electron Production Rate 1
1 On the Physics toolbar, click Domains and choose Electron Production Rate.
3 In the Settings window for Electron Production Rate, locate the Electron Production
Rate section.
4 In the Re text field, type Re.
Wall 1
3 Select the Include migration effects check box.
Electron Density and Energy 1

1 On the Physics toolbar, click Boundaries and choose Electron Density and Energy.
3 In the Settings window for Electron Density and Energy, locate the Electron Density and
Energy section.
4 Select the Fix mean electron energy check box.
Wall 2
3 In the Settings window for Wall, locate the Electron Density Wall Settings section.
4 In the tn text field, type influx.
5 Locate the Electron Energy Wall Settings section. Clear the Use wall for electron energy
check box.
The above boundary conditions applied to boundary 1 means that the mean electron
energy will be fixed at 3 electron volts and there will be an influx of electrons resulting
in a net current influx of 2E-6 amps.
When you create the mesh you want to use a finer mesh close to the walls so that the
sharp gradients in the electron density and electron energy density are adequately
resolved. You accomplish this by using a graded mapped mesh.
MESH 1
Mapped.
Mapped 1
In the Model Builder window, under Component 1 (comp1)>Mesh 1 right-click Mapped 1 and
Distribution 1
Size
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 click Size.
4 Click Build All.
The mean electron energy and electron density can change on a sub nanosecond
timescale. The electron density and mean electron energy reach their steady state values
very quickly so it is only necessary to solve the problem for 1 microsecond.
STUDY 1

Dependent.
3 In the Times text field, type 0.
4 Click Range.
7 In the Stop text field, type -6.
10 Click Add.
RESULTS
Electron Density (dd)

1 In the Model Builder window, under Results click Electron Density (dd).
2 In the Settings window for 2D Plot Group, click to expand the Window settings section.
3 Click to expand the Color legend section. Locate the Color Legend section. From the
Position list, choose Bottom.
5 On the Electron Density (dd) toolbar, click Plot.
Electron Temperature (dd)

1 In the Model Builder window, under Results click Electron Temperature (dd).
2 In the Settings window for 2D Plot Group, locate the Color Legend section.
3 From the Position list, choose Bottom.
5 On the Electron Temperature (dd) toolbar, click Plot.
Cut Line 2D 1
1 On the Results toolbar, click Cut Line 2D.
3 In the Settings window for Cut Line 2D, locate the Line Data section.
4 In row Point 1, set y to 2.5e-4.
5 In row Point 2, set y to 2.5e-4 and x to 5-3.
6 Click Plot.
1D Plot Group 3
3 From the Data set list, choose Cut Line 2D 1.
Line Graph 1
of the y-axis data section. From the menu, choose Component 1>Drift Diffusion>Electron
density>dd.ne - Electron density.
1D Plot Group 4
3 From the Data set list, choose Cut Line 2D 1.
Line Graph 1
of the y-axis data section. From the menu, choose Component 1>Drift Diffusion>Electron
energy density>dd.Te - Electron temperature.
Electrodeless Lamp
Introduction
This model simulates an electrodeless lamp with argon/mercury chemistry. The low
excitation threshold for mercury atoms means that even though the mercury is present in
small concentrations, its interaction with electrons determines the overall discharge
characteristics. here is strong UV emission from the plasma at 185 nm and 253 nm
stemming from spontaneous decay of electronically excited mercury atoms. The UV
emission can stimulate phosphors coated on the surface of the bulb resulting in visible
light. From an electrical point of view, the lamp can be thought of as a transformer, where
the coil acts as the primary and the plasma acts as the secondary. If the efficiency of
discharge lamps could be increased by 1%, it would result in a saving of 109 kWh per year
worldwide.
Note: This application requires the Plasma Module and the AC/DC Module.
Model Definition
A schematic of the geometry used to solve the problem is given in Figure 1. A sinusoidal
current is applied to the copper coil (green) which creates a magnetic field in the ferrite
core (gray). When the plasma ignites, a magnetic circuit is created between the ferrite core
and the plasma. The free electrons in the plasma bulk are accelerated by the electric field.
This leads to creation of new electrons through ionization which sustains the plasma. In
2 | ELECTRODELESS LAMP
quasi steady-state, the creation of new electrons is balanced by the loss of electrons to the
wall.
475 mtorr Ar
25 mtorr Hg
Dielectric
material
Copper coil
Ferrite core
Figure 1: Diagram of electrodeless light source.
The presence of mercury leads to the formation of electronically excited mercury atoms.
Certain excited states emit a photon at a given wavelength with a certain emission
frequency. By solving for the number density of each of the excited species, you can
determine the amount of energy channeled into creating the excited mercury atoms. You
can then calculate the amount of energy emitted from the plasma as photons.
In order to simplify the analysis, the following assumptions are made:
The model is assumed to be axially symmetric.

The AC induction currents are solved in the frequency domain.
The electron energy distribution function (EEDF) is assumed to be Maxwellian.
Thermal quenching of excited atoms is not considered.
Energy losses in the ferrite core are not considered.
A trapping factor is used to specify an effective emission photon frequency for the
excited mercury atoms. These trapping factors are based on published data.
The electronically excited argon species are lumped into a single species.
PLASMA CHEMISTRY
The chemical mechanism comes from Ref. 1 and consists of 11 species and 96 reactions:

1 e+Ar=>e+Ar Momentum 0
2 e+Ar=>e+Ar* Excitation 11.56
+
3 e+Ar=>e+e+Ar Ionization 15.80
4 e+Ar*=>e+Ar Superelastic -11.56
+
5 e+Ar*=>e+e+Ar Ionization 4.24
6 e+Hg=>e+Hg Momentum 0
7 e+Hg=>e+Hg(63P0) Excitation 4.66
11 e+Hg=>e+Hg(73S1) Excitation 7.70
12 e+Hg=>e+Hg(63DJ) Excitation 8.85
+
13 e+Hg=>e+e+Hg Ionization 10.44
14 e+Hg(63P0)=>e+Hg(63P0) Momentum 0
15 e+Hg(63P0)=>e+Hg Superelastic -4.66
16 e+Hg(63P0)=>e+Hg(63P1) Excitation 0.21
19 e+Hg(63P0)=>e+Hg(73S1) Excitation 3.04
20 e+Hg(63P0)=>e+Hg(63DJ) Excitation 4.18
+
21 e+Hg(63P0)=>e+e+Hg Ionization 5.78
24 e+Hg(63P1)=>e+Hg(63P0) Superelastic -0.21
+
29 e+Hg(63P1)=>e+e+Hg Ionization 5.57

37 e+Hg(63P2)=>e+e+Hg+ Ionization 5.01
43 e+Hg(61P1)=>e+e+Hg+ Ionization 3.74
44 e+Hg(73S1)=>e+Hg(73S1) Momentum 0
45 e+Hg(73S1)=>e+Hg Superelastic -7.7
46 e+Hg(73S1)=>e+Hg(63P0) Superelastic -3.04
49 e+Hg(73S1)=>e+e+Hg+ Ionization 2.74
50 e+Hg(63DJ)=>e+Hg(63DJ) Momentum 0
51 e+Hg(63DJ)=>e+Hg Superelastic -8.85
52 e+Hg(63DJ)=>e+Hg(63P0) Superelastic -4.19
55 e+Hg(63DJ)=>e+e+Hg+ Ionization 1.59
56 Ar*+Ar*=>e+Ar+Ar+ Penning 0
57 Ar*+Hg=>e+Ar+Hg+ Penning 0
58 Ar*+Hg(63P0)=>e+Ar+Hg+ Penning 0
+
60 Ar*+Hg(63P2)=>e+Ar+Hg Penning 0

62 Ar*+Hg(73S1)=>e+Ar+Hg+ Penning 0
63 Ar*+Hg(63DJ)=>e+Ar+Hg+ Penning 0
64 Hg(63P2)+Hg(63P2)=>e+Hg+Hg+ Penning 0
66 Hg(63P2)+Hg(73S1)=>e+Hg+Hg+ Penning 0
67 Hg(63P2)+Hg(63DJ)=>e+Hg+Hg+ Penning 0
72 Hg(61P1)+Hg(73S1)=>e+Hg+Hg+ Penning 0
73 Hg(61P1)+Hg(63DJ)=>e+Hg+Hg+ Penning 0
74 Hg(73S1)+Hg(63P0)=>e+Hg+Hg+ Penning 0
78 Hg(73S1)+Hg(73S1)=>e+Hg+Hg+ Penning 0
79 Hg(73S1)+Hg(63DJ)=>e+Hg+Hg+ Penning 0
80 Hg(63DJ)+Hg(63P0)=>e+Hg+Hg+ Penning 0
84 Hg(63DJ)+Hg(73S1)=>e+Hg+Hg+ Penning 0
85 Hg(63DJ)+Hg(63DJ)=>e+Hg+Hg+ Penning 0
+ +
86 Ar +Hg=>Hg +Ar Charge exchange 0
+ +
87 Ar +Ar=>Ar+Ar Charge exchange 0
+ +
88 Hg +Hg=>Hg+Hg Charge exchange 0
89 Hg(63P1)=>Hg 253nm 0
90 Hg(61P1)=>Hg 185nm 0
91 Hg(73S1)=>Hg(63P0) 405nm 0

92 Hg(73S1)=>Hg(63P1) 436nm 0
93 Hg(73S1)=>Hg(63P2) 546nm 0
94 Hg(63DJ)=>Hg(63P0) 297nm 0
95 Hg(63DJ)=>Hg(63P1) - 0
96 Hg(63DJ)=>Hg(63P2) 365nm 0
The following surface reactions are considered:
TABLE 2: SURFACE REACTIONS.
REACTION FORMULA
1 Ars=>Ar
2 Ar+=>Ar
3 Hg1=>Hg
4 Hg2=>Hg
5 Hg3=>Hg
6 Hg4=>Hg
7 Hg5=>Hg
8 Hg6=>Hg
9 Hg+=>Hg
( n ) + [ n ( E ) D n ] = R
t e e e e e e
n
( ) + [ n ( E ) D n ] + E e = R
t
D e = e T e, = --- e, D = T e
5
3
Re = x j kj Nn ne
j=1
R = xj kj Nn ne j
j=1

kk = 0 k ( )f ( ) d
( w k ) + ( u )w k = j k + R k
t
0 r V =
N
= q Z k n k n e

k = 1
For a nonmagnetized and nonpolarized plasma, the induction currents are computed in
the frequency domain using the following equation:
2 1 e
( j 0 )A + ( 0 A ) = J
The plasma conductivity needs to be specified as a material property, usually from the cold
plasma approximation:
2
ne q
= -------------------------------
m e ( e + j )
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 lamp is operated by a fixed power of 80 watts. This means that the total power
dissipation in the system is 80 W. Some of the power is lost in the coil and the ferrite but
the bulk of the power is channeled into the plasma.
The results are presented below.
Figure 2: Surface plot of electron density inside the column.
The electron density is plotted in Figure 2. The electron density is high, as one would
expect in an inductively coupled plasma. The peak value of the electron density at the
driving frequency used in the model results in a peak plasma conductivity of around
180 S/m. The high value for the electron density and low excitation and ionization
threshold for mercury results in a very low electron temperature which is plotted in
Figure 3. The peak electron temperature is only 1.27 eV, which through Boltzmanns
relation results in a low plasma potential. The plasma potential is plotted in Figure 4 and
only peaks at 8 V.
Figure 3: Plot of the electron temperature.
Figure 4: Plot of plasma potential.
Figure 5: Plot of the resistive losses.
Figure 6: Plot of the mole fraction of ground state mercury.
Figure 7: Plot of the number density of mercury ions.
The resistive losses in the plasma are plotted in Figure 5. The plasma skin depth is a few
centimeters so there is no real shielding of the azimuthal electric field. The mole fraction
of the ground state mercury is plotted in Figure 6. The mole fraction is low in the core of
the plasma and higher on the walls. This is because the electrons consume the ground state
mercury in the plasma core, converting it to electronically excited states. The electronically
excited mercury atoms diffuse to the walls of the lamp where they de-excite back to the
ground state. This continuous consumption of ground state mercury in the plasma bulk
and release on the walls results in large gradients in mole fraction within the bulb.
There are two ion species present in the plasma, argon, and mercury. Despite the fact that
the number density of ground state argon is 25 times higher than mercury, the density of
mercury ions is several hundred times greater than the density of argon ions. This is
because the ionization energy for mercury is only 10.44 eV compared to 15.7 eV for
argon. Direct ionization of mercury is preferable to argon because the tail of the electron
energy distribution function drops dramatically at higher electron energies. Additionally,
any argon ions which encounter a ground state or electronically excited mercury atom
donate their charge because it is energetically favorable.
Figure 8: Plot of the mole fraction of Hg(63P1). Spontaneous decay of this species is responsible
for the generation of 253nm radiation.
The mole fraction of Hg(63P1) is plotted in Figure 8. These atoms spontaneously emit
photos at a frequency factor of 8106 s1. On the way to the walls of the lamp, the photons
continuously excite mercury atoms and then be released when spontaneous decay occurs.
This resonant absorption and re-absorption of the photons means that the frequency
factor appears to be much lower than it actually is. Since a self consistent model of the
radiation imprisonment of the photons is computationally impractical, a trapping factor is
used to approximate this effect. A trapping factor of 10 is used for the Hg(63P1) atoms
which means that the frequency factor is lowered by a factor of 10. In Figure 9 the mole
fraction of Hg(61P1) is plotted. A trapping factor of 1000 is used for the spontaneous
decay back to ground state mercury.
Figure 9: Plot of the mole fraction of Hg(61P1). Spontaneous decay of this species is responsible
for the generation of 185nm radiation.
Reference
1. K. Rajaraman, Radiation Transport in Low Pressure Plasmas: Lighting and
Semiconductor Etching Plasmas, Ph.D. thesis, Depart. of Physics, University of Illinois,
2005.

electrodeless_lamp
NEW
MODEL WIZARD
3 Click Add.
4 Click Study.
6 Click Done.
GEOMETRY 1

4 Find the Control points subsection. In row 2, set r to 0.015.
5 Find the Added segments subsection. Click Add Cubic.
6 Find the Control points subsection. In row 2, set z to 0.025.
9 Find the Added segments subsection. Click Add Quadratic.
11 In row 3, set r to 0 and z to 0.07.
13 Find the Control points subsection. In row 2, set z to 0.
Rectangle 1 (r1)
Rectangle 2 (r2)
Rectangle 3 (r3)
Chamfer 1 (cha1)
1 On the Geometry toolbar, click Chamfer.
2 Click the Zoom In button on the Graphics toolbar.
3 On the object r2, select Point 3 only.
4 In the Settings window for Chamfer, locate the Distance section.
5 In the Distance from vertex text field, type 1.5e-3.
Chamfer 2 (cha2)
1 On the Geometry toolbar, click Chamfer.
2 On the object r3, select Point 3 only.
3 In the Settings window for Chamfer, locate the Distance section.
4 In the Distance from vertex text field, type 2e-3.
Square 1 (sq1)
1 On the Geometry toolbar, click Primitives and choose Square.
2 In the Settings window for Square, locate the Size section.
3 In the Side length text field, type 0.001.
Array 1 (arr1)
1 On the Geometry toolbar, click Transforms and choose Array.
2 Select the object sq1 only.
3 In the Settings window for Array, locate the Size section.
4 From the Array type list, choose Linear.
5 In the Size text field, type 5.
6 Locate the Displacement section. In the z text field, type 5e-3.
DEFINITIONS
Variables 1

Tinit 350[K] K Gas temperature
pinit 101325*(500E-3/ Pa Initial total
760)[Pa] pressure
lamp_power 80[W] W Lamp power
tf1 10 Trapping factor
tf2 1000 Trapping factor
Explicit 1
2 Select Domains 59 only.
3 In the Settings window for Explicit, locate the Output Entities section.
4 From the Output entities list, choose Adjacent boundaries.
5 In the Label text field, type Coil boundaries.
Explicit 2
3 In the Settings window for Explicit, type Coil domains in the Label text field.
Explicit 3
4 Select Boundaries 8, 27, and 3538 only.
5 In the Label text field, type Boundary layers.
Explicit 4
3 In the Settings window for Explicit, type Discharge in the Label text field.
PLASMA (PLAS)

3 Click Browse.
Ar_xsecs.txt.

3 Click Browse.
Hg_xsecs.txt.
Reaction 1
3 In the Formula text field, type Ars+Ars=>Ar++Ar+e.
4 Locate the Reaction Parameters section. In the kf text field, type N_A_const*
1.00E-15[m^3/s].
Reaction 2
3 In the Formula text field, type Ars+Hg=>Hg++Ar+e.
9E-16[m^3/s].
58: Ars+Hg=>Hg++Ar+e
1 Right-click Reaction 2 and choose Duplicate.
3 In the Formula text field, type Ars+Hg1=>Hg++Ar+e.
59: Ars+Hg1=>Hg++Ar+e
1 Right-click Component 1 (comp1)>Plasma (plas)>58: Ars+Hg=>Hg++Ar+e and choose
Duplicate.
1 Right-click Component 1 (comp1)>Plasma (plas)>59: Ars+Hg1=>Hg++Ar+e and choose
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
3 In the Formula text field, type Hg3+Hg3=>Hg++Hg+e.
3.50E-16[m^3/s].
65: Hg3+Hg3=>Hg++Hg+e
Duplicate.
1 Right-click Component 1 (comp1)>Plasma (plas)>65: Hg3+Hg3=>Hg++Hg+e and choose
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
Duplicate.
3 In the Formula text field, type Ar++Hg=>Hg++Ar.
1.50E-17[m^3/s].
87: Ar++Hg=>Hg++Ar
Duplicate.
3 In the Formula text field, type Ar++Ar=>Ar++Ar.
4.60E-16[m^3/s].
88: Ar++Ar=>Ar++Ar
1 Right-click Component 1 (comp1)>Plasma (plas)>87: Ar++Hg=>Hg++Ar and choose
Duplicate.
3 In the Formula text field, type Hg++Hg=>Hg+Hg+.
1.00E-15[m^3/s].
89: Hg++Hg=>Hg+Hg+
1 Right-click Component 1 (comp1)>Plasma (plas)>88: Ar++Ar=>Ar++Ar and choose
Duplicate.
3 In the Formula text field, type Hg2=>Hg.
4 Locate the Reaction Parameters section. In the kf text field, type 8.00E6/tf1.
90: Hg2=>Hg
1 Right-click Component 1 (comp1)>Plasma (plas)>89: Hg++Hg=>Hg+Hg+ and choose
Duplicate.
4 Locate the Reaction Parameters section. In the kf text field, type 7.50E8/tf2.
91: Hg4=>Hg
1 Right-click Component 1 (comp1)>Plasma (plas)>90: Hg2=>Hg and choose Duplicate.
3 In the Formula text field, type Hg5=>Hg1.
92: Hg5=>Hg1
93: Hg5=>Hg2
1 Right-click Component 1 (comp1)>Plasma (plas)>92: Hg5=>Hg1 and choose Duplicate.
4 Locate the Reaction Parameters section. In the kf text field, type 2E7.
94: Hg5=>Hg3
95: Hg6=>Hg1
96: Hg6=>Hg2
Surface Reaction 1
4 Locate the Boundary Selection section. From the Selection list, choose Boundary layers.
2: Ars=>Ar
1 Right-click Surface Reaction 1 and choose Duplicate.
3: Ar+=>Ar
1 Right-click Component 1 (comp1)>Plasma (plas)>2: Ars=>Ar and choose Duplicate.
4: Hg1=>Hg
1 Right-click Component 1 (comp1)>Plasma (plas)>3: Ar+=>Ar and choose Duplicate.
5: Hg2=>Hg
6: Hg3=>Hg
7: Hg4=>Hg
8: Hg5=>Hg
9: Hg6=>Hg
3 In the Formula text field, type Hg+=>Hg.
Species: Hg
1 In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click Species: Hg.
3 In the Mw text field, type 0.2006.
4 In the text field, type 2.969[angstrom].
5 In the /kb text field, type 750.
6 In the x0 text field, type 0.05.
Species: Hg1
Hg1.
Species: Hg2
Hg2.
Species: Hg3
Hg3.
Species: Hg4
Hg4.
Species: Hg5
Hg5.
Species: Hg6
Hg6.
Species: Ar
Species: Ar+
Ar+.
3 In the n0 text field, type 1E16.
Species: Hg+
Hg+.
6 Locate the Species Formula section. Select the Initial value from electroneutrality
constraint check box.
Plasma Model 1
Model 1.
3 In the e text field, type 4E24[1/(m*V*s)]/plas.Nn.
4 Locate the Model Inputs section. In the T text field, type Tinit.
5 In the pA text field, type pinit.
Ground 1
3 From the Selection list, choose Boundary layers.
Wall 1
Initial Values 1
Values 1.
5 Click the Zoom Box button on the Graphics toolbar.
In the Model Builder window, under Component 1 (comp1) click Magnetic Fields (mf).
Coil 1
3 From the Selection list, choose Coil domains.
4 Locate the Coil section. Select the Coil group check box.
5 From the Coil excitation list, choose Power.
6 In the Pcoil text field, type lamp_power.
MATERIALS
Blank Material.
Material 1 (mat1)

group

5 In the Rename Material dialog box, type Coils in the New label text field.
6 Click OK.
Material 2 (mat2)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.

group
Relative permeability mur 1e3 1 Basic

6 In the Rename Material dialog box, type Ferrite in the New label text field.
7 Click OK.
Material 3 (mat3)

group

6 In the Rename Material dialog box, type Air in the New label text field.
7 Click OK.
Material 4 (mat4)

group

7 Click OK.
MESH 1
2 From the Element size list, choose Extra fine.
Edge 1
2 In the Settings window for Edge, locate the Boundary Selection section.
3 From the Selection list, choose Coil boundaries.
Distribution 1
3 From the Selection list, choose Coil boundaries.
distribution type.
Mapped 1
4 From the Selection list, choose Coil domains.
Edge 2
1 Right-click Mesh 1 and choose More Operations>Edge.
Size 1
5 In the associated text field, type 5e-4.
Edge 3
2 Select Boundaries 3638 only.
Size 1
5 In the associated text field, type 1e-3.
Boundary Layers 1
5 Click to expand the Transition section. Clear the Smooth transition to interior mesh check
box.

section.
4 Locate the Boundary Layer Properties section. In the Boundary layer stretching factor text
field, type 1.1.
5 In the Thickness adjustment factor text field, type 1.
6 From the Thickness of first layer list, choose Manual.
7 In the Thickness text field, type 2E-5.
Free Triangular 1
STUDY 1
2 In the Frequency text field, type 2[MHz].
4 Click Range.
7 In the Stop text field, type log10(2e-3).
10 Click Add.
RESULTS

1 In the Model Builder window, expand the Results>Data Sets node.
2 Right-click Study 1/Solution 1 (sol1) and choose Duplicate.
Selection
2D Plot Group 6
field.
Surface 1
2 In the Settings window for Surface, locate the Data section.
3 From the Data set list, choose Study 1/Solution 1 (2) (sol1).
4 Click Replace Expression in the upper-right corner of the Expression section. From the
menu, choose Component 1>Magnetic Fields>Heating and losses>mf.Qrh - Volumetric loss
density, electric.
2D Plot Group 7
2 In the Settings window for 2D Plot Group, type Ground State Mercury Mole
Fraction in the Label text field.
Surface 1
1 Right-click Ground State Mercury Mole Fraction and choose Surface.
Transport)>Mole fractions>plas.x_wHg - Mole fraction.
3 On the Ground State Mercury Mole Fraction toolbar, click Plot.
2D Plot Group 8
2 In the Settings window for 2D Plot Group, type Mercury Ion Number Density in the
Label text field.
Surface 1
1 Right-click Mercury Ion Number Density and choose Surface.
Transport)>Number densities>plas.n_wHg_1p - Number density.
3 On the Mercury Ion Number Density toolbar, click Plot.
2D Plot Group 9
2 In the Settings window for 2D Plot Group, type Mole Fraction of Excited Mercury
2 in the Label text field.
Surface 1
1 Right-click Mole Fraction of Excited Mercury 2 and choose Surface.
Transport)>Mole fractions>plas.x_wHg2 - Mole fraction.
3 On the Mole Fraction of Excited Mercury 2 toolbar, click Plot.
2D Plot Group 10
2 In the Settings window for 2D Plot Group, type Mole Fraction of Excited Mercury
4 in the Label text field.
Surface 1
1 Right-click Mole Fraction of Excited Mercury 4 and choose Surface.
Transport)>Mole fractions>plas.x_wHg4 - Mole fraction.
3 On the Mole Fraction of Excited Mercury 4 toolbar, click Plot.
ICP Torch
Introduction
Thermal plasmas have nowadays a large range of industrial applications including: cutting,
welding, spraying, waste destruction and surface treatment. Thermal plasmas are assumed
to be under partial to complete local thermodynamic equilibrium (LTE) conditions.
Under LTE, the plasma can be considered a conductive fluid mixture and therefore, be
modeled using the magnetohydrodynamics (MHD) equations. This model shows how to
use the Equilibrium Inductively Coupled Discharge interface to simulate the plasma
generated in an inductively coupled plasma torch.
Figure 1 displays the geometry of the to-be-modeled inductive plasma torch.
Figure 1: Geometry of an inductively coupled plasma torch. The torch is composed of three
concentric quartz tubes in which gas are injected from the bottom and exit from the top the
torch. In this model, a fixed power of 11 kW is transfered to the plasma by a three turn coil
operating at 3MHz.
Note: This application requires the Plasma Module and AC/DC Module.
Model Definition
This model is based on the work presented in Ref. 1 and uses the following assumptions:
1 The plasma torch is modeled by a fully axisymmetric configuration.
2 | I C P TO R C H
2 The coil consists of parallel current carrying rings with a circle cross section, 6 mm in
diameter. This implies neglecting the axial component of the coil current.
3 Steady state, laminar pure argon plasma flow at atmospheric pressure.
4 Optically thin plasma under local thermodynamic equilibrium (LTE) conditions.
5 Viscous dissipation and pressure work in the energy equation are neglected.
Figure 2 shows the geometry of the model.
Figure 2: Schematic of the ICP torch. Flow enters from the base (v1, v2 and v3) and leaves out
the top. The dimensions of the different part of the model are given in the Modeling
Instructions section.
In this model excitation is provided to a three turns coil at 3 MHz. The gas flowing in the
sheath tube (plasma confinement tube) is then ionized by Joule heating.
The model is solved using a frequency-transient study in combination with a single turn
coil feature which set a fixed power to the system (11 kW). By fixing the power, the current
and electric potential can vary in the coil as the plasma electrical conductivity builds up.
Steady-state is reached when the coil current stabilized to its nominal value.
In this model the three different gas stream velocities (v1 for the carrier tube, v2 for the
central tube and v3 for the sheath tube) are composed of pure argon. The temperature
dependent argons physical properties are loaded from the material library under
Equilibrium Discharge. Note that the temperature range of the physical properties span
3 | I C P TO R C H
from 500 K to 25 000K. Note also that a minimum electrical conductivity has been used
to initiate the plasma. The latter has been set to 1 S/m.

Figure 3, and Figure 4 respectively shows the plasma temperature distribution, and
velocity magnitude of the argon plasma after 0.1 s. Figure 5 shows the electrical
conductivity of the plasma at the same time (0.1 s). Note that, for this figure, the electrical
conductivity of the other constituents of the model has been set to 0 for sake of
visualization.
Figure 6 displays the magnetic flux norm at steady state (0.1 s). Note that the electrical
conductivity of the plasma screens the magnetic flux as a consequence of the skin effect.
Figure 7 shows the coil current as a function of the simulation time. The steady state is
reached when the current stabilizes, i.e. around t=0.1 s.
4 | I C P TO R C H
Figure 3: Surface plot of the LTE plasma temperature.
Figure 4: Plot of the velocity magnitude.
5 | I C P TO R C H
Figure 5: Plot of the plasma electrical conductivity.
Figure 6: Norm of the magnetic flux. Note the effect of the resistivity on the penetration of the
field (skin effect).
6 | I C P TO R C H
Figure 7: Coil current as a function of time for a fixed excitation power. Note the stabilization
of the current density as the system reach the steady state.
Reference
1. S. Xue, P. Proulx and M.I. Boulos, Extended-field electromagnetic model for
inductively coupled plasma, J. Phys. D. 34, 1897, 2001
Application Library path: Plasma_Module/Equilibrium_Discharges/icp_torch
NEW
7 | I C P TO R C H
MODEL WIZARD
2 In the Select Physics tree, select Plasma>Equilibrium Discharges>Equilibrium Inductively
Coupled Plasma.
3 Click Add.
4 Click Study.
6 Click Done.
ROOT
Select the mm units.
GEOMETRY 1
1 In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose mm.
GLOBAL DEFINITIONS
Parameters
8 | I C P TO R C H

T0 300[K] 300 K Ambient temperature
Pext 11[kW] 1.1E4 W Coil excitation
power
freq 3[MHz] 3E6 Hz Coil excitation
frequency
r_3 125[mm] 0.125 m Axial length:
Computational domain
L_3 200[mm] 0.2 m Height:
Computational domain
and sheath tube
d_1 2[mm] 0.002 m Thickness: Carrier
tube
L_0 50[mm] 0.05 m Height: Carrier tube
and central tube
r_1 3.7[mm] 0.0037 m Inner radius:
Carrier tube
d_2 2.2[mm] 0.0022 m Thickness: Central
tube
r_2 18.8[mm] 0.0188 m Inner radius:
Central tube
d_3 3.5[mm] 0.0035 m Thickness: Sheath
tube
r_0 25[mm] 0.025 m Inner radius: Sheath
tube
d_c 6[mm] 0.006 m Diameter: Coils
r_c 33[mm] 0.033 m Axial length: Center
of the coils
L_1 63[mm] 0.063 m Height: Center of
the lower coil
L_2 121[mm] 0.121 m Height: Center of
the upper coil
Q_1 1[l/min] 1.667E-5 m/s Gas stream: Carrier
tube
Q_2 3[l/min] 5E-5 m/s Gas stream: Central
tube
Q_3 31[l/min] 5.167E-4 m/s Gas stream: Sheath
tube
M 0.04[kg/mole] 0.04 kg/mol Molar mass: Argon
9 | I C P TO R C H
mv_stp 22.4[l/mole] 0.0224 m/mol Molar volume at stp
mdot1 M*Q_1/mv_stp 2.976E-5 kg/s Mass flow rate:
Carrier tube
Central tube
Sheath tube
rho_stp 1.91[kg/m^3] 1.91 kg/m Density of argon at
stp
A1 pi*(r_1)^2 4.301E-5 m Cross section:
Carrier gas stream
A2 pi*(r_2^2-(r_1+ 0.001008 m Cross section:
d_1)^2) Central gas stream
A3 pi*(r_0^2-(r_2+ 5.781E-4 m Cross section:
d_2)^2) Sheath gas stream
v1 mdot1/rho_stp/A1 0.3623 m/s Velocity: Carrier
gas stream
v2 mdot2/rho_stp/A2 0.04636 m/s Velocity: Central
gas stream
v3 mdot3/rho_stp/A3 0.8356 m/s Velocity: Sheath gas
stream
Define the computational domain.
GEOMETRY 1
Rectangle 1 (r1)
3 In the Width text field, type r_3.
4 In the Height text field, type L_3.
Define the carrier tube.
Rectangle 2 (r2)
3 In the Width text field, type d_1.
10 | I C P TO R C H
5 Locate the Position section. In the r text field, type r_1.
Define the central tube.
Rectangle 3 (r3)
Define the sheat tube.
Rectangle 4 (r4)
Define the coils.
Circle 1 (c1)
1 On the Geometry toolbar, click Primitives and choose Circle.
2 In the Settings window for Circle, locate the Size and Shape section.
3 In the Radius text field, type d_c/2.
4 Locate the Position section. In the r text field, type r_c.
5 In the z text field, type L_1.
Circle 2 (c2)
5 In the z text field, type (L_1+L_2)/2.
Circle 3 (c3)
11 | I C P TO R C H
5 In the z text field, type L_2.
Define a ramp function that will be used to bring the initial gas stream to their steady
state values.
DEFINITIONS
Ramp 1 (rm1)
1 On the Home toolbar, click Functions and choose Global>Ramp.
Set the slope of the ramp function such that the steady state values are reached over a
short interval of time, e.g. one excitation cycle. The steady state values are reached when
the ramp function is equal to unity. Check the cutoff box and make sure the value of the
cutoof is set to 1 in order to limit the value of the ramp function to one.
2 In the Settings window for Ramp, locate the Parameters section.
3 In the Slope text field, type freq.
4 Select the Cutoff check box.
Define the different domain type for easy selection.
Explicit 1
3 In the Rename Explicit dialog box, type Air in the New label text field.
4 Click OK.
Explicit 2
3 In the Rename Explicit dialog box, type Plasma in the New label text field.
4 Click OK.
Explicit 3
12 | I C P TO R C H
3 In the Rename Explicit dialog box, type Quartz in the New label text field.
4 Click OK.
Explicit 4
3 In the Rename Explicit dialog box, type Coils in the New label text field.
4 Click OK.
Add the different materials used in the model using the material library.
ADD MATERIAL
1 On the Home toolbar, click Add Material to open the Add Material window.
2 Go to the Add Material window.
3 In the tree, select Built-In>Air.
MATERIALS
Air (mat1)
1 In the Model Builder window, under Component 1 (comp1)>Materials click Air (mat1).
3 From the Selection list, choose Air.
ADD MATERIAL
2 In the tree, select AC/DC>Copper.
MATERIALS
Copper (mat2)
1 In the Model Builder window, under Component 1 (comp1)>Materials click Copper (mat2).
13 | I C P TO R C H
ADD MATERIAL
2 In the tree, select AC/DC>Quartz.
MATERIALS
Quartz (mat3)
1 In the Model Builder window, under Component 1 (comp1)>Materials click Quartz (mat3).
3 From the Selection list, choose Quartz.
ADD MATERIAL
2 In the tree, select Equilibrium Discharge>Argon.
MATERIALS
Argon (mat4)
1 In the Model Builder window, under Component 1 (comp1)>Materials click Argon (mat4).
3 From the Selection list, choose Plasma.
4 On the Home toolbar, click Add Material to close the Add Material window.
Adjust the selection and features of each physics composing the model.
The magnetic field interface is used over the whole computational domain. The
Single-Turn Coil feature is used here to transfer the excitation power to the plasma.

the menu.
2 In the Model Builder window, under Component 1 (comp1) click Magnetic Fields (mf).
3 In the Settings window for Magnetic Fields, click to expand the Discretization section.
4 From the Magnetic vector potential list, choose Linear.
14 | I C P TO R C H
Coil 1
4 Locate the Coil section. Select the Coil group check box.
6 In the Pcoil text field, type Pext.
The heat transfer in the air is neglected in this model.
H E A T TR A N S F E R I N F L U I D S ( H T )
1 In the Model Builder window, under Component 1 (comp1) click Heat Transfer in Fluids
(ht).
Solid 1
1 On the Physics toolbar, click Domains and choose Solid.
Initial Values 1
1 In the Model Builder window, under Component 1 (comp1)>Heat Transfer in Fluids (ht)
click Initial Values 1.
2 In the Settings window for Initial Values, type T0 in the T text field.
Add a heat transfer in solids feature for the solid part of the heat transfer model (tubes
and coils).
Temperature 1
1 On the Physics toolbar, click Boundaries and choose Temperature.
3 In the Settings window for Temperature, locate the Temperature section.
4 In the T0 text field, type T0.
The single phase flow is only applied to the plasma.
LAMINAR FLOW (SPF)

Since the density variation is not small, the flow can not be regarded as incompressible.
Therefore set the flow to be compressible. Add some isotropic diffusion which is initially
bery high then ramps down to zero after the plasma has ignited.
15 | I C P TO R C H
1 In the Model Builder window, under Component 1 (comp1) click Laminar Flow (spf).
2 In the Settings window for Laminar Flow, locate the Physical Model section.
3 From the Compressibility list, choose Compressible flow (Ma<0.3).
4 Locate the Domain Selection section. From the Selection list, choose Plasma.
5 Click to expand the Equation section. From the Equation form list, choose Stationary.
the menu.
7 Click to expand the Consistent stabilization section. Locate the Consistent Stabilization
section. Find the Navier-Stokes equations subsection. Clear the Crosswind diffusion check
box.
8 Click to expand the Inconsistent stabilization section. Locate the Inconsistent
Stabilization section. Select the Isotropic diffusion check box.
9 In the id text field, type 10-10*tanh(1000*(t-0.08)).
Inlet 1
1 On the Physics toolbar, click Boundaries and choose Inlet.
Add the inlets with their proper velocities. Note the presence of the ramp function. The
later is used to ramp the velocities from their initial values to the steady state.
3 In the Settings window for Inlet, locate the Velocity section.
4 In the U0 text field, type v1*rm1(t[1/s]).
Inlet 2
Inlet 3
Outlet 1
1 On the Physics toolbar, click Boundaries and choose Outlet.
16 | I C P TO R C H
3 In the Settings window for Outlet, locate the Pressure Conditions section.
4 Clear the Suppress backflow check box.
The Lorentz force can be neglected in this model. Neglecting this force also makes the
model much easier to solve.
MULTIPHYSICS
In the Model Builder window, under Component 1 (comp1)>Multiphysics right-click Lorentz
Force 1 (lf1) and choose Disable.
MESH 1
Edit Physics-Induced Sequence.
Size 1
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 click Size 1.
Size 2
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 click Size 2.
Boundary Layers 2
2 Right-click Boundary Layers 2 and choose Move Up.

1 Click the Select Box button on the Graphics toolbar.
3 Select Boundaries 2132 only.
17 | I C P TO R C H
4 In the Settings window for Boundary Layer Properties, locate the Boundary Layer
Properties section.
5 In the Number of boundary layers text field, type 4.
6 From the Thickness of first layer list, choose Manual.
7 In the Thickness text field, type 8[um].
8 Click Build All.
STUDY 1
1 In the Model Builder window, expand the Study 1 node, then click Step 1:
3 In the Times text field, type range(0,0.05,1)*0.1.
4 In the Frequency text field, type freq.
Change some solver settings.
Solution 1 (sol1)
3 Right-click Time-Dependent Solver 1 and choose Segregated.
4 Right-click Study 1>Solver Configurations>Solution 1 (sol1)>Time-Dependent Solver 1>
Segregated 1 and choose Lower Limit.
5 In the Settings window for Lower Limit, locate the Lower Limit section.
6 In the Lower limits (field variables) text field, type comp1.T 300.
7 Right-click Segregated 1 and choose Segregated Step.
8 In the Settings window for Segregated Step, locate the General section.
9 Under Variables, click Add.
10 In the Add dialog box, In the Variables list, choose Pressure (comp1.p) and Velocity field
(comp1.u).
11 Click OK.
12 In the Settings window for Segregated Step, click to expand the Method and termination
section.
18 | I C P TO R C H
13 Locate the Method and Termination section. From the Jacobian update list, choose On
every iteration.
14 In the Model Builder window, under Study 1>Solver Configurations>Solution 1 (sol1)>
Time-Dependent Solver 1>Segregated 1 click Segregated Step.
15 In the Settings window for Segregated Step, locate the General section.
16 In the Variables list, choose Pressure (comp1.p) and Velocity field (comp1.u).
17 Under Variables, click Delete.
18 Locate the Method and Termination section. From the Jacobian update list, choose On
every iteration.
19 In the Model Builder window, under Study 1>Solver Configurations>Solution 1 (sol1) click
Time-Dependent Solver 1.
20 In the Settings window for Time-Dependent Solver, click to expand the Time stepping
section.
21 Locate the Time Stepping section. From the Steps taken by solver list, choose
Intermediate.
22 Select the Nonlinear controller check box.
Create some relevant figures.
The temperature.
RESULTS
2D Plot Group 1
2 In the Model Builder window, right-click 2D Plot Group 1 and choose Rename.
3 In the Rename 2D Plot Group dialog box, type Temperature in the New label text field.
4 Click OK.
Surface 1
1 Right-click Temperature and choose Surface.
3 In the Expression text field, type T.
4 On the Temperature toolbar, click Plot.
Duplicate the figure to display the fluid velocity magnitude.
19 | I C P TO R C H
Temperature
In the Model Builder window, under Results right-click Temperature and choose Duplicate.
Temperature 1
1 In the Model Builder window, under Results right-click Temperature 1 and choose
Rename.
2 In the Rename 2D Plot Group dialog box, type Velocity in the New label text field.
3 Click OK.
Surface 1
1 In the Model Builder window, expand the Results>Velocity node, then click Surface 1.
3 In the Expression text field, type spf.U.
4 On the Velocity toolbar, click Plot.
Duplicate the figure to display the electrical conductivity.
Velocity
In the Model Builder window, under Results right-click Velocity and choose Duplicate.
Velocity 1
1 In the Model Builder window, under Results right-click Velocity 1 and choose Rename.
2 In the Rename 2D Plot Group dialog box, type Electrical conductivity in the New
label text field.
3 Click OK.
Surface 1
1 In the Model Builder window, expand the Results>Electrical conductivity node, then click
Surface 1.
3 In the Expression text field, type mf.sigmarr*(dom==1).
4 On the Electrical conductivity toolbar, click Plot.
Duplicate the figure to display the norm of the magnetic flux. Note the effect of the
plasma conductivity on the skin depth.
Electrical conductivity
In the Model Builder window, under Results right-click Electrical conductivity and choose
Duplicate.
20 | I C P TO R C H
Electrical conductivity 1
1 In the Model Builder window, under Results right-click Electrical conductivity 1 and
choose Rename.
2 In the Rename 2D Plot Group dialog box, type Magnetic flux in the New label text field.
3 Click OK.
Surface 1
1 In the Model Builder window, expand the Results>Magnetic flux node, then click Surface
1.
3 In the Expression text field, type mf.normB.
4 On the Magnetic flux toolbar, click Plot.
Display the coil current as a function of time. Note the time it takes to get the steady
state (constant current in the coils).
1D Plot Group 5
3 In the Rename 1D Plot Group dialog box, type Coil current in the New label text field.
4 Click OK.
Global 1
1 On the Coil current toolbar, click Global.

mf.ICoil_1 A Coil current
4 On the Coil current toolbar, click Plot.

Create a nice 3D plot for the model thumbnail.
Create first a revolution data set.
Revolution 2D 1
1 On the Results toolbar, click More Data Sets and choose Revolution 2D.
2 In the Settings window for Revolution 2D, click to expand the Revolution layers section.
3 Locate the Revolution Layers section. In the Start angle text field, type -90.
21 | I C P TO R C H
4 In the Revolution angle text field, type 225.
Then create the 3D plot.
3D Plot Group 6
3 In the Rename 3D Plot Group dialog box, type Temperature 3D in the New label text
field.
4 Click OK.
Volume 1
1 Right-click Temperature 3D and choose Volume.
2 In the Settings window for Volume, locate the Expression section.
3 In the Expression text field, type T.
4 On the Temperature 3D toolbar, click Plot.
5 Click the Go to Default 3D View button on the Graphics toolbar.
Set the figure as a model thumbnail by clicking on the root folder in the model builder
than expand the model thumbnail section and click on set model thumbnail.
22 | I C P TO R C H
In-Plane Microwave Plasma
Introduction
Wave-heated discharges can be very simple, such as when a plane wave is guided into a
reactor using a waveguide, or very complicated, such as in the case with ECR (electron
cyclotron resonance) reactors. In this simple example, a wave is launched down a
waveguide where it intersects a flowing gas at low pressure, resulting in formation of an
argon plasma. Microwave plasmas typically have high number density without requiring
significant power absorption. The plasma potential is also quite low compared to capacitive
or DC discharges. Therefore, microwave plasmas share many of the characteristics of
inductive discharges.
Electromagnetic wave in,

2.45 GHz
Gas flow in, 1 Torr
Figure 1: Diagram of geometry modeled. A TE or TM mode wave enters from the top port and
intersects the gas flow leading to the formation of a plasma.
Note: The model requires the Plasma Module and RF Module.
Model Definition
n
( ) + [ n e ( e E ) D e n e ] = R e
t e
2 | IN-PLANE MICROWAVE PLASMA

( n ) + [ n ( E ) D n ] + E = R
t e
D e = e T e, = --- e, D = T e
5
3
Re = x j kj Nn ne
j=1
R = xj kj Nn ne j
j=1

kk = 0 k ( )f ( ) d
( w k ) + ( u )w k = j k + R k
t

0 r V =
N
= q Z k n k n e

k = 1
In a microwave reactor the high frequency electric field is computed in the frequency
domain using the following equation:
( r E ) k 0 r --------- E = 0
1 2 j
0
The relationship between the plasma current density and the electric field becomes more
complicated in the presence of a DC magnetic field. The following equation defines this
relationship:
1
J = E
Here, is the plasma conductivity tensor, which is a function of the electron density,
collision frequency, and magnetic flux density. Using the definitions:
q
= -------------------------------, = n e q
m e ( e + j )
where q is the electron charge, me is the electron mass, ne is the collision frequency, and
is the angular frequency of the electromagnetic field. In this example, the inverse of the
plasma conductivity is diagonal because there is no external DC magnetic field:
1
--- 0 0

1
= 0 1 --- 0

1
0 0 ---


The gas flow is modeled assuming a constant velocity in the x-direction.
PLASMA CHEMISTRY
The chemical mechanism for the plasma consists of only 3 species and 7 reactions:

electron. Reaction number 7 is responsible for heating of the gas. The 11.5 eV of energy
which was consumed in creating the electronically excited argon atom is returns to the gas
as thermal energy when the excited metastable quenches. In addition to volumetric
reactions, the following surface reactions are implemented:
1 Ars=>Ar 1
2 Ar+=>Ar 1
The plasma is sustained through absorption of electromagnetic waves. The Port boundary
condition is used to excite the plasma. A total absorbed power of 30 W is fed into the port.
In a second study, the electrical excitation is changed to the TM mode, where the electric
field has only an in-plane component. The total absorbed power is the same as the TE
mode case.

The electron density is plotted in Figure 2 and peaks slightly downstream of the crossing
point. The electron density is also slightly asymmetric in the y-plane due to the fact that
the electromagnetic waves are absorbed asymmetrically. The electron temperature is
plotted in Figure 3. The electron temperature is relatively low everywhere, in part due to
the high operating pressure (1 Torr). The electron temperature peaks directly
underneath the waveguide where the wave is absorbed. The norm of the electric field can
be seen in Figure 5. The electric field is high inside the waveguide and there are no losses.
Once the wave is exposed to the plasma, the energy is absorbed by the electrons, raising
the electron temperature enough to generate new electrons through ionization. The
ionization rate is high enough to sustain the plasma. The contour of the critical plasma
density is also plotted in Figure 5. The electromagnetic wave cannot penetrate into regions
exceeding the critical plasma density. Since the electron temperature is relatively low,
one would expect the plasma potential to be low. The plasma potential is plotted in
Figure 4 and is only around 10 volts.
In the TE mode, electrons do not experience any change in the high-frequency electric
field during the microwave time scale. This means that the phase coherence between the
electrons and electromagnetic waves is only destroyed through collisions with the
background gas. The loss of phase coherence between the electrons and high-frequency
fields is what results in energy gain for the electrons. Therefore, the momentum collision
frequency is simply given by:
m = e
where e is the collision frequency between the electrons and neutrals.
When switching to the TM mode, the electron density, electron temperature and plasma
potential are quite similar to the TE mode case. This can be seen in Figure 7, Figure 8 and
Figure 9. The electric field is very different however, Figure 9. The electric field cannot
penetrate past the contour of critical electron density, and has its greatest magnitude in this
location. The power deposition, Figure 10 is also highly localized to the contour of critical
electron density. The TM mode causes in-plane motion of the electrons on the microwave
time scale, so in regions where the high-frequency electric field is significant (the contour
where the electron density is equal to the critical density), the time-averaged electric field
experienced by the electrons may be non-zero. This destroys the phase coherence between
the electrons and the fields, causing the electrons to gain energy. This is an example of a
non-local kinetic effect, which is difficult to approximate with a fluid model. However,
since this effect is similar to collisions with a background gas, the non-local effects can be

approximated by adding an effective collision frequency to the momentum collision
frequency:
m = e + eff
where eff is the effective collision frequency to account for non-local effects. In this
example, since the Doppler broadening parameter is set to 20, this corresponds to an
effective collision frequency of:

eff = ------
20

Figure 2: Plot of the electron density.
Figure 3: Plot of the electron temperature in the reactor.

Figure 4: Plot of the plasma potential in the reactor.\
Figure 5: Plot of the electric field norm. The white contour represents the critical plasma

density, where the electron density is equal to 7.6E16[1/m3].
Figure 6: Plot of the power deposition into the plasma. The white contour represents the critical
plasma density, where the electron density is equal to 7.6E16[1/m3].

Figure 7: Electron density for the TM mode case.

Figure 8: Plot of the electron temperature for the TM mode case.
Figure 9: Close up of the high frequency electric field norm for the TM mode case.
Figure 10: Close up of the power deposition into the plasma for the TM mode case.

Reference
1. M.A. Lieberman and A.J. Lichtenberg, Principles of Plasma Discharges and Materials
Processing, John Wiley & Sons, 2005.
Application Library path: Plasma_Module/Wave-Heated_Discharges/

inplane_microwave_plasma
NEW
MODEL WIZARD
2 In the Select Physics tree, select Plasma>Microwave Plasma.
3 Click Add.
4 Click Study.
6 Click Done.
GLOBAL DEFINITIONS
Parameters

P0 30[W] 30 W Absorbed power

GEOMETRY 1
Rectangle 1 (r1)
Rectangle 2 (r2)
5 Locate the Position section. In the x text field, type -0.1.
6 In the y text field, type 0.1.
Rectangle 3 (r3)
5 Locate the Position section. In the y text field, type 0.15.
DEFINITIONS
Explicit 1

7 Click OK.
PLASMA (PLAS)

3 Click Browse.
Ar_xsecs.txt.
Reaction 1
Reaction 2
Species: Ar
Species: Ar+
Ar+.

Surface Reaction 1
Surface Reaction 2
7 Select the Convection check box.
check box.
Plasma Model 1
Model 1.
3 Specify the u vector as
10 x
0 y
4 In the T text field, type 350.

5 In the pA text field, type 1[torr].
6 Locate the Electron Density and Energy section. In the eNn text field, type 4E24.
Initial Values 1
Values 1.

Wall 1
Ground 1
Electron Outlet 1
1 On the Physics toolbar, click Boundaries and choose Electron Outlet.
Species: Ars
In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click Species: Ars.
Outflow 1
1 On the Physics toolbar, click Attributes and choose Outflow.
2 In the Settings window for Outflow, locate the Boundary Selection section.
Species: Ar+
In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click Species: Ar+.
Outflow 1

Components section.
3 From the Electric field components solved for list, choose Out-of-plane vector.

4 In the Model Builder window, click Electromagnetic Waves, Frequency Domain (emw).
Port 1
1 On the Physics toolbar, click Boundaries and choose Port.
3 In the Settings window for Port, locate the Port Properties section.
4 From the Type of port list, choose Rectangular.
5 From the Wave excitation at this port list, choose On.
6 Select the Specify deposited power check box.
7 In the Pdep text field, type P0.
MATERIALS
Blank Material.
Material 1 (mat1)

group
Electrical conductivity epsilonr 5 1 Basic
Relative permittivity sigma 0 S/m Basic
MESH 1
Boundary Layers.
Boundary Layers 1


section.
Size
3 From the Calibrate for list, choose Plasma.
4 From the Predefined list, choose Fine.
5 Click Build All.
STUDY 1
3 Click Range.
9 Click Add.
11 In the Frequency text field, type 2.45[GHz].
RESULTS
Electric Field (emw)

In the Model Builder window, under Results right-click Electric Field (emw) and choose
Contour.
Contour 1
1 In the Settings window for Contour, locate the Levels section.
2 From the Entry method list, choose Levels.

3 In the Levels text field, type 7.6E16.
4 Locate the Coloring and Style section. Clear the Color legend check box.
5 From the Coloring list, choose Uniform.
6 From the Color list, choose White.
7 On the Electric Field (emw) toolbar, click Plot.
Electric Field (emw)

Right-click Electric Field (emw) and choose Duplicate.
Electric Field (emw) 1

1 In the Model Builder window, expand the Electric Field (emw) 1 node.
2 Right-click Results>Electric Field (emw) 1 and choose Rename.
3 In the Rename 2D Plot Group dialog box, type Resistive Heating in the New label text
field.
4 Click OK.
Surface 1
the Expression section. From the menu, choose Component 1>Electromagnetic Waves,
Frequency Domain>Heating and losses>emw.Qrh - Resistive losses.
2 On the Resistive Heating toolbar, click Plot.
Now change to the in-plane electric field, which makes the problem much more difficult
to solve. This is because all the power will be absorbed on the contour of critical electron
density. Setting a Doppler broadening parameter of 20 smooths out the region over
which power is deposited to help with convergence.
MULTIPHYSICS
1 In the Model Builder window, under Component 1 (comp1)>Multiphysics click Plasma
Conductivity Coupling 1 (pcc1).
2 In the Settings window for Plasma Conductivity Coupling, locate the Compute Tensor
Plasma Conductivity section.
3 Select the Compute tensor plasma conductivity check box.


Components section.
3 From the Electric field components solved for list, choose In-plane vector.
Port 1
1 In the Model Builder window, under Component 1 (comp1)>Electromagnetic Waves,
Frequency Domain (emw) click Port 1.
2 In the Settings window for Port, locate the Port Mode Settings section.
3 From the Mode type list, choose Transverse magnetic (TM).
ADD STUDY
3 Find the Studies subsection. In the Select Study tree, select Preset Studies for
Multiphysics>Frequency-Transient.
STUDY 2
2 In the Model Builder window, under Study 2 click Step 1: Frequency-Transient.
5 Click Range.
11 Click Add.
13 In the Frequency text field, type 2.45[GHz].

RESULTS
Electric Field (emw) 1

1 In the Model Builder window, under Results click Electric Field (emw) 1.
3 From the Data set list, choose Study 2/Solution 2 (sol2).
4 On the Electric Field (emw) 1 toolbar, click Plot.
5 Click the Zoom In button on the Graphics toolbar.
Resistive Heating
In the Model Builder window, under Results right-click Resistive Heating and choose
Duplicate.
Resistive Heating 1
3 On the Resistive Heating 1 toolbar, click Plot.
1D Plot Group 11
2 In the Settings window for 1D Plot Group, type Port Power in the Label text field.
Global 1
1 On the Port Power toolbar, click Global.
the y-axis data section. From the menu, choose Component 1>Electromagnetic Waves,
Frequency Domain>Global>emw.Pin - Port input power.
3 On the Port Power toolbar, click Plot.
4 Click to expand the Legends section. From the Legends list, choose Manual.
Legends
TE Mode
Port Power
1 In the Model Builder window, under Results click Port Power.
2 In the Settings window for 1D Plot Group, click to expand the Legend section.

Global 1
In the Model Builder window, under Results>Port Power right-click Global 1 and choose
Duplicate.
Global 2
1 In the Settings window for Global, locate the Data section.
Legends
TM Mode
5 On the Port Power toolbar, click Plot.

Ion Energy Distribution Function
Introduction
One of the most useful quantities of interest after solving a self-consistent plasma model is
the ion energy distribution function (IEDF). The magnitude and shape of the IEDF
depends on many of the discharge parameters; pressure, plasma potential, sheath width
and so forth. At very low pressures the plasma sheath is said to be collisionless, meaning
that the ion energy is not retarded by collisions with the background gas. At higher
pressures the ions collide with the background gas molecules in the sheath and their energy
at the moment of impact with a surface is reduced.
Note: This application requires the Plasma Module, Particle Tracing Module, and AC/
DC Module.
Model Definition
The equations of motion for ions in an electric field and background gas are
d ( mv ) = qZE
dt
where m is the ion mass (SI unit: kg), v is the particle velocity (SI unit: m/s), q is unit
charge (SI unit: C), Z is the ion charge number (SI unit: dimensionless).
When an ion undergoes a collision with the background gas, its velocity vector changes.
The probability of a collision event occurring depends in the ion-neutral collision
frequency, (SI unit: 1/s) which is defined as:
= Nd vp vg
where Nd is the background number density (SI unit: 1/m3), is the ion-neutral charge
exchange collision cross section (SI unit: m2) and vp is the particle velocity and vg is the
velocity of the background gas atoms or molecules. In this example the collision cross
section is assumed to be constant, 6x10-19m2, as given in Ref. 2.
The collision probability defined as:
P = 1 exp ( t )
If P is greater than a random number between 0 and 1 then the particle velocity is
reinitialized to the following expression:
2 | ION ENERGY DISTRIBUTION FUNCTION

v p' = ( m p v p + m g v g m g v g v p R ) ( m p + m g )
where vp is the precollision particle velocity, mg is the mass of the background gas atoms
or molecules, R is a uniformly distributed random unit vector, and vg is the velocity of the
background gas atoms or molecules, which is sampled from a Maxwellian distribution
function:
2
1 ( vg u )
f ( v g ) = --------------------------------------
- exp -------------------------------
-
2k B T m g
32 ( 2k B T ) m g
where T is the temperature of the background gas.
The particles are released 5 mm away from the wafer surface and 20 mm radially inwards
from the center of the reactor. There are 22,500 ions modeled. They all start at the same
point in space and have an initial Maxwellian distribution function, which is for each
velocity direction
2
m - exp ----------
mv i
f ( vi ) = -------------- -
2kT 2kT
where T is the initial temperature of the ions, in this case 400K. This results in a
distribution function of
nsdim
f(v) = f ( vi )
i=1
When the ions strike the wafer their velocity is frozen for all subsequent timesteps. This
allows the velocity and energy distribution function to be recovered once all the ions have
made contact with the wall.
The angle at which the ions strike the wafer is also of interest. This can be recovered by
plotting a histogram of the inverse tangent of the radial and axial particle velocities.

The electric potential for an argon plasma at an operating pressure of 20mtorr is plotted
in Figure 1. The initial starting position for the ions is also shown:
Initial ion position
Figure 1: Plot of the plasma potential used to compute the IEDF.
The IEDF is plotted in Figure 2. Most of the ions have kinetic energy between around 20
and 22 eV. At the initial starting coordinate, the plasma potential is 20.98V , so it is
expected that the kinetic energy of the ions at the wall is of similar order. As the pressure
decreases the collision frequency between ions and neutrals is reduced, so the IEDF should
shift to the right, towards the maximum value given by the plasma potential. As the
pressure increases the sheath becomes more collisional which inhibits the ions from
reaching higher energies. Of course, changing the pressure results in a change in the

discharge characteristics which may alter the shape and magnitude of the IEDF in a
nonlinear way.
Figure 2: Plot of the ion energy distribution function (IEDF) in an inductively coupled
plasma.
The angle at which the ions are striking the wafer surface are plotted in Figure 3. Due to
the fact that the ions are released relatively close to the wafer surface, the range of angles
is small, typically between -20 and 20 degrees. The plot is not quite symmetric due to the
presence of a small outward component of the ambipolar electric field.

Figure 3: Plot of the angle at which the ions strike the surface of the wafer.
The ion angular energy distribution function is plotted in Figure 4. The angle is not quite
symmetric about zero due to the profile of the electric field at the release point.

Figure 4: Plot of the ion angular energy distribution function.
Notes About the COMSOL Implementation

This model is most conveniently solved by opening an existing model in the Plasma
Module Application Library, then computing the ion trajectories.
References
1. O.V. Vozniy, G.Y. Yeom, and A. Yu. Kropotov, Plasma Potential Influence on Ion
Energy Distribution Function in ICP Source, PSE, vol 5, no 1-2, pp. 2833, http://
www.pse.scpt.org.ua/en/jornal/1-2_07/3.pdf.
2. M. Surendra, Radiofrequency Discharge Benchmark Model Comparison, Plasma

Sources Sci. Technol, vol. 4, pp 5673, 1995.
3. A. V. Phelps, The application of scattering cross sections to ion flux models in

discharge sheaths, J. Appl. Phys. vol. 76, pp. 747-753, 1994.

ion_energy_distribution_function
Start by opening the model of the inductively coupled GEC reference cell from the Plasma
Module Application Library.
From the File menu, choose Open.
Browse to the applications Application Libraries folder and double-click the file
argon_gec_icp.mph.
COMPONENT 1 (COMP1)
Now add a Charged Particle Tracing interface to compute the ion energy distribution
function.
ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the tree, select AC/DC>Particle Tracing>Charged Particle Tracing (cpt).
5 On the Home toolbar, click Add Physics to close the Add Physics window.
ADD STUDY
3 Find the Studies subsection. In the Select Study tree, select Preset Studies.
4 Find the Physics interfaces in study subsection. In the table, clear the Solve check box for
Plasma (plas) and Magnetic Fields (mf).
5 Find the Studies subsection. In the Select Study tree, select Preset Studies>Time
Dependent.

STUDY 2

On the Home toolbar, click Add Study to close the Add Study window.
C H A R G E D PA R T I C L E TR A C I N G ( C P T )
1 In the Model Builder window, expand the Component 1 (comp1) node, then click Charged
Particle Tracing (cpt).
2 In the Settings window for Charged Particle Tracing, locate the Domain Selection
section.
GLOBAL DEFINITIONS
Parameters

Mw 0.04[kg/mol] 0.04 kg/mol Ion molecular weight
mi Mw/N_A_const 6.642E-26 kg Ion mass
DEFINITIONS
Analytic 1 (an1)
1 On the Home toolbar, click Functions and choose Global>Analytic.
Enter the analytic approximation for momentum cross section for elastic scattering
between Ar+ ions and neutral Ar atoms, which which depends on the kinetic energy of
the particles.
2 In the Settings window for Analytic, type Qm in the Function name text field.
3 Locate the Definition section. In the Expression text field, type 1.15e-18*x^(-0.1)*
(1+0.015/x)^0.6.
4 Locate the Units section. In the Arguments text field, type eV.
5 In the Function text field, type m^2.

Analytic 2 (an2)
1 On the Home toolbar, click Functions and choose Global>Analytic.
Enter the analytic approximation for isotropic elastic collision between Ar+ ions and
neutral Ar atoms from, which which depends on the kinetic energy of the particles.
2 In the Settings window for Analytic, type Qi in the Function name text field.
3 Locate the Definition section. In the Expression text field, type 2e-19/(x^(0.5)*(1+
x))+3e-19*x/(1+x/3)^(2.3).
4 Locate the Units section. In the Arguments text field, type eV.
5 In the Function text field, type m^2.
C H A R G E D PA R T I C L E TR A C I N G ( C P T )
Particle Properties 1
1 In the Model Builder window, under Component 1 (comp1)>Charged Particle Tracing (cpt)
click Particle Properties 1.
2 In the Settings window for Particle Properties, locate the Particle Mass section.
3 In the mp text field, type mi.
4 Locate the Charge Number section. In the Z text field, type 1.
Now add an Electric Force feature. The electric potential comes from the solved plasma
model. Specify that piecewise polynomial recovery should be used when computing the
electric force. This results in a more accurate reconstruction of the electric field.
Electric Force 1
1 On the Physics toolbar, click Domains and choose Electric Force.
3 In the Settings window for Electric Force, locate the Electric Force section.
4 From the Specify force using list, choose Electric potential.
5 From the V list, choose Electric potential (plas/pes1).
6 Locate the Advanced Settings section. Select the Use piecewise polynomial recovery on field
check box.
Now you add the collisional force between the ions and background gas. The collision
frequency is a function of the neutral number density (plas.Nn), the elastic and charge
exchange cross section and the particle velocity (cpt.V).
Collisions 1
1 On the Physics toolbar, click Domains and choose Collisions.

3 In the Settings window for Collisions, locate the Fluid Properties section.
4 In the Nd text field, type plas.Nn.
5 In the Mg text field, type Mw.
6 In the T text field, type plas.T.
Elastic 1
1 On the Physics toolbar, click Attributes and choose Elastic.
2 In the Settings window for Elastic, locate the Collision Frequency section.
3 In the text field, type Qi(cpt.Ep).
Collisions 1
In the Model Builder window, under Component 1 (comp1)>Charged Particle Tracing (cpt)
click Collisions 1.
Resonant Charge Exchange 1

1 On the Physics toolbar, click Attributes and choose Resonant Charge Exchange.
2 In the Settings window for Resonant Charge Exchange, locate the Collision Frequency
section.
3 In the text field, type (Qm(cpt.Ep)-Qi(cpt.Ep))/2.
Release from Grid 1

1 On the Physics toolbar, click Global and choose Release from Grid.
2 In the Settings window for Release from Grid, locate the Initial Coordinates section.
3 In the qr, 0 text field, type 0.02.
4 In the qz, 0 text field, type 0.005.
5 Locate the Initial Velocity section. From the Initial velocity list, choose Maxwellian.
6 In the Nv text field, type 150.
7 In the T0 text field, type 400.
STUDY 2

1 In the Model Builder window, under Study 2 click Step 1: Time Dependent.
3 In the Times text field, type range(0,5.0e-6/30,5.0e-6).

4 Click to expand the Values of dependent variables section. Locate the Values of Dependent
Variables section. Find the Values of variables not solved for subsection. From the Settings
list, choose User controlled.
5 From the Method list, choose Solution.
6 From the Study list, choose Study 1, Frequency-Transient.
RESULTS
1D Plot Group 13
2 In the Settings window for 1D Plot Group, type Ion Energy Distribution
Function in the Label text field.
3 Locate the Data section. From the Data set list, choose Particle 1.
Histogram 1
1 On the Ion Energy Distribution Function toolbar, click Histogram.
2 In the Settings window for Histogram, click Replace Expression in the upper-right corner
of the Expression section. From the menu, choose Component 1>Charged Particle
Tracing>Velocity and energy>cpt.Ep - Particle kinetic energy.
3 Locate the Expression section. From the Unit list, choose eV.
4 Locate the Bins section. From the Entry method list, choose Limits.
5 In the Limits text field, type range(0,25/500,25).
6 On the Ion Energy Distribution Function toolbar, click Plot.
1D Plot Group 14
2 In the Settings window for 1D Plot Group, type Ion Angular Distribution
Histogram 1
1 On the Ion Angular Distribution Function toolbar, click Histogram.
2 In the Settings window for Histogram, locate the Expression section.
3 In the Expression text field, type atan(cpt.vr/cpt.vz).

4 Locate the Bins section. In the Number text field, type 500.
5 Locate the Expression section. From the Unit list, choose .
6 On the Ion Angular Distribution Function toolbar, click Plot.
Ion Angular Distribution Function

1 In the Model Builder window, under Results click Ion Angular Distribution Function.
4 In the associated text field, type Ion angle of incidence (deg).
5 On the Ion Angular Distribution Function toolbar, click Plot.
2D Plot Group 15
2 In the Settings window for 2D Plot Group, type Ion Angular Energy Distribution
4 From the Time (s) list, choose 5E-6.
Histogram 1
1 On the Ion Angular Energy Distribution Function toolbar, click More Plots and choose
Histogram.
2 In the Settings window for Histogram, locate the x-Expression section.
3 In the Expression text field, type atan(cpt.vr/cpt.vz).
4 From the Unit list, choose .
5 Select the Description check box.
6 In the associated text field, type Angle of incidence (deg).
7 Locate the y-Expression section. In the Expression text field, type cpt.Ep.
8 From the Unit list, choose eV.
9 Locate the Bins section. Find the x bins subsection. From the Entry method list, choose
Limits.
10 In the Limits text field, type range(-10,20/79,10).
11 Find the y bins subsection. From the Entry method list, choose Limits.
12 In the Limits text field, type range(0,22/79,22).
Ion Angular Energy Distribution Function

1 In the Model Builder window, under Results click Ion Angular Energy Distribution Function.

3 Clear the Plot data set edges check box.
4 On the Ion Angular Energy Distribution Function toolbar, click Plot.

DC Glow Discharge, 1D
Introduction
DC glow discharges in the low pressure regime have long been used for gas lasers and
fluorescent lamps. DC discharges are attractive to study because the solution is time
independent. This model shows how to use the Plasma interface to set up an analysis of a
positive column. The discharge is sustained by emission of secondary electrons at the
cathode.
Model Definition
The DC discharge consists of two electrodes, one powered (the anode) and one grounded
(the cathode):
Cathode Anode
Figure 1: Schematic of the DC discharge. The voltage applied across the electrodes leads to
formation of a plasma.
DOMAIN EQUATIONS

( n ) + [ n e ( e E ) D e n e ] = R e
t e
n
( ) + [ n ( E ) D n ] + E e = R
t
where:
e = ( e E )n e D e n e
The electron diffusivity, energy mobility and energy diffusivity are computed from the
2 | DC GLOW DISCHARGE, 1D
D e = e T e, = --- e, D = T e
5
3
using rate coefficients. Suppose that there are M reactions which contribute to the growth
Re = x j kj Nn ne
j=1
reaction j (m3/s), and Nn is the total neutral number density (1/m3). For DC discharges
it is better practice to use Townsend coefficients instead of rate coefficients to define
reaction rates. Townsend coefficients provide a better description of what happens in the
cathode fall region Ref 1. When Townsend coefficients are used, the electron source term
is given by:
Re = xj j Nn e
j=1
where j is the Townsend coefficient for reaction j (m2) and e is the electron flux as
defined above (1/(m2s)). 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 electron energy loss is obtained by summing the collisional energy loss over all
reactions:
R = xj kj Nn ne j
j=1
where j is the energy loss from reaction j (V). The rate coefficients may be computed
from cross section data by the following integral:

kk = 0 k ( )f ( ) d
the collision cross section (m2) and f is the electron energy distribution function. In this
case a Maxwellian EEDF is assumed. When Townsend coefficients are used, the electron
energy loss is taken as:
R = xj j Nn e j
j=1
( w k ) + ( u )w k = j k + R k
t
0 r V =
N
= q Z k n k n e

k = 1
BOUNDARY CONDITIONS
Unlike RF discharges, the mechanism for sustaining the discharge is emission of secondary
electrons from the cathode. An electron is emitted from the cathode surface with a
specified probability when struck by an ion. These electrons are then accelerated by the
strong electric field close to the cathode where they acquire enough energy to initiate
ionization. The net result is a rapid increase in the electron density close to the cathode in
a region often known as the cathode fall or Crookes dark space.
n e = --- e, th n e
1
2 p ( p n ) (1)
p
n = --- e, th n
5
6 p p ( p n ) (2)
p
PLASMA CHEMISTRY

In this discharge, the electron density and density of excited species is relatively low so
1 Ars=>Ar 1
2 Ar+=>Ar 1
The electric potential, electron density, and mean electron energy are all quantities of
interest. Most of the variation in each of these quantities occurs along the axial length of
the column. Figure 2 plots the electron density in the column. The electron density peaks
in the region between the cathode fall and positive column. This region is sometimes
referred to as Faraday dark space. The electron density obtained in this 1D model is
different than that obtained in the 2D model due to the fact that the diffusive loss of
electrons to the outer walls and the accumulation of surface charge on the walls is not
modeled.
Figure 2: Plot of electron density inside the column.
In Figure 3 the electric potential is plotted along the axial length of the column. Notice
that the potential profile is markedly different from the linear drop in potential which
results in the absence of the plasma. The strong electric field in the cathode region can lead
to high energy ion bombardment of the cathode. Heating of the cathode surface occurs
which may in turn lead to thermal electron emission where additional electrons are emitted
from the cathode surface.
Figure 3: Plot of the electron temperature along the axial length of the positive column.
Figure 4: Plot of electron temperature along the axial length of the positive column.
Figure 5: Plot of the mass fraction of excited argon atoms.
Figure 6: Plot of the number density of argon ions.
The plasma current due to electrons, ions and their sum is plotted in Figure 7. As
expected, the ion current is highest at the cathode and increases sharply in the cathode fall
region. The ion bombardment of the cathode results in an electron current released from
the electrode. The electron current increases sharply in the cathode fall region because the
high electron temperature results in production of new electrons which then contribute to
the total electron current. Once the electrons pass the cathode fall region the electron
current density further increases due to further production of new electrons through
electron impact ionization with the background gas.
Figure 7: Plot of the electron current density (blue), the ion current density (green) and the
total current density (red) in the positive column.
Reference

positive_column_1d
NEW
MODEL WIZARD
3 Click Add.
4 Click Study.
6 Click Done.
the menu.
GEOMETRY 1
The geometry interval is defined to be consistent with the 2D version of the model, which
is available in the Application Library.
Interval 1 (i1)
DEFINITIONS
Variables 1

ramp tanh(1E7[1/s]*t) Ramp function for the
applied voltage
mueN 1E25[1/(m*V*s)] 1/(Vms) Electron mobility
V0 100[V]*ramp V Applied voltage
Wf 5 Surface work function
p0 0.1[torr] Pa Gas pressure
PLASMA (PLAS)

3 Click Browse.
Ar_xsecs.txt.

Because you will examine the electron, ion, and net currents flowing in the plasma, raise
the element order to 2. The current density is computed from space derivatives of the
charge carrying degrees of freedom, so using 2nd order shape functions gives a more
accurate value for the current density.
8 Click to expand the Discretization section. From the Formulation list, choose Finite
element, log formulation (quadratic shape function).
Plasma Model 1
Model 1.
Now change the way the source coefficients for electronic excitation and ionization are
specified. By default, COMSOL Multiphysics computes rate coefficients based on the
cross section data you supplied. For DC discharges, Townsend coefficients provide a
more accurate description of the cathode fall region so they should be used. The
Townsend coefficients are typically computed using the Boltzmann Equation, Two-Term
Approximation interface.
2: e+Ar=>e+Ars
1 In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click 2: e+Ar=>
e+Ars.
2 In the Settings window for Electron Impact Reaction, locate the Collision section.
3 From the Specify reaction using list, choose Use lookup table.
4 Locate the Reaction Parameters section. From the Rate constant form list, choose
Townsend coefficient.
town2.txt.
4: e+Ar=>2e+Ar+
2e+Ar+.
town4.txt.
Reaction 1
Reaction 2
Species: Ar
When solving a reacting flow problem there always needs to be one species which is
selected to fullfill the mass constraint. This should be taken as the species with the
largest mass fraction.
Species: Ars
Ars.
electroneutrality constraint is satisfied.
Species: Ar+
Ar+.
Now add a surface reaction which describes the neutralization of argon ions on the
electrode. Secondary emission of electrons is required to sustain the discharge, so enter
the emission coefficient and an estimate of the mean energy of the secondary electrons
based on the ionization energy threshold and the work function of the surface.
Surface Reaction 1
Make the secondary emission coefficient 0.35 and set the mean energy of the secondary
electrons to be the ionization energy (given by the expression plas.de_4) minus twice
the work function of the electrode.
6 In the i text field, type plas.de_4-2*Wf.
Surface Reaction 2
Surface Reaction 3
Wall 1
Ground 1
Metal Contact 1
3 In the Settings window for Metal Contact, locate the Terminal section.
4 In the V0 text field, type V0.
MESH 1
Edge.
Edge 1
Distribution 1
7 Click Build All.
STUDY 1

3 Click Range.
9 Click Add.
RESULTS

1 In the Model Builder window, under Results click Electron Density (plas).
7 In the associated text field, type Electron density (1/m3).

1 In the Model Builder window, under Results click Electron Temperature (plas).
7 In the associated text field, type Electron temperature (eV).

7 In the associated text field, type Potential (V).
1D Plot Group 4
3 Locate the Data section. From the Time selection list, choose Last.
7 In the associated text field, type Mass fraction of excited Argon (1).
Line Graph 1
1 On the Excited Argon Mass Fraction toolbar, click Line Graph.
1D Plot Group 5
2 In the Settings window for 1D Plot Group, type Argon Ion Number Density in the
Label text field.
7 In the associated text field, type Argon ion number density (1/m3).
Line Graph 1
1 On the Argon Ion Number Density toolbar, click Line Graph.
4 On the Argon Ion Number Density toolbar, click Plot.
Argon Ion Number Density

In the Model Builder window, under Results right-click Argon Ion Number Density and
choose Duplicate.
Argon Ion Number Density 1

1 In the Settings window for 1D Plot Group, type Current Density in the Label text
field.
3 Locate the Plot Settings section. In the y-axis label text field, type Current density
(A/m2).
Line Graph 1
1 In the Model Builder window, expand the Results>Current Density node, then click Line
Graph 1.
Current>Electron current density>plas.Jelx - Electron current density, x component.
3 Click to expand the Quality section. From the Resolution list, choose No refinement.
4 From the Recover list, choose Within domains.
7 On the Current Density toolbar, click Plot.
Legends
Electron current density
Line Graph 2
1 Right-click Results>Current Density>Line Graph 1 and choose Duplicate.
Transport)>Species>Species wAr_1p>Ion current density>plas.Jix_wAr_1p - Ion current
density, x component.
Legends
Ion current density
Line Graph 1
Right-click Line Graph 1 and choose Duplicate.
Line Graph 3
2 In the Expression text field, type plas.Jix_wAr_1p+plas.Jelx.
Legends
Total current density
DC Glow Discharge
Introduction
DC glow discharges in the low pressure regime have long been used for gas lasers and
fluorescent lamps. DC discharges are attractive to study because the solution is time
independent. This model shows how to use the Plasma interface to set up an analysis of a
positive column. The discharge is sustained by emission of secondary electrons at the
cathode.
Model Definition
The DC discharge consists of two electrodes, one powered (the anode) and one grounded
(the cathode). The positive column is coupled to an external circuit:
1000
Cathode Plasma Anode V

1 pF
Figure 1: Schematic of the DC discharge and external circuit.
DOMAIN EQUATIONS
n
( ) + [ n e ( e E ) D e n e ] = R e
t e
n
( ) + [ n ( E ) D n ] + E e = R
t
where:
e = ( e E )n e D e n e
2 | DC GLOW DISCHARGE
The electron diffusivity, energy mobility and energy diffusivity are computed from the
D e = e T e, = --- e, D = T e
5
3
using rate coefficients. Suppose that there are M reactions which contribute to the growth
Re = x j kj Nn ne
j=1
reaction j (m3/s), and Nn is the total neutral number density (1/m3). For DC discharges
it is better practice to use Townsend coefficients instead of rate coefficients to define
reaction rates. Townsend coefficients provide a better description of what happens in the
cathode fall region Ref 1. When Townsend coefficients are used, the electron source term
is given by:
Re = xj j Nn e
j=1
where j is the Townsend coefficient for reaction j (m2) and e is the electron flux as
defined above (1/(m2s)). 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 electron energy loss is obtained by summing the collisional energy loss over all
reactions:
R = xj kj Nn ne j
j=1
where j is the energy loss from reaction j (V). The rate coefficients may be computed
from cross section data by the following integral:

kk = 0 k ( )f ( ) d
the collision cross section (m2) and f is the electron energy distribution function. In this
case a Maxwellian EEDF is assumed. When Townsend coefficients are used, the electron
energy loss is taken as:
R = xj j Nn e j
j=1
( w k ) + ( u )w k = j k + R k
t
0 r V =
N
= q Z k n k n e

k = 1
BOUNDARY CONDITIONS
Unlike RF discharges, the mechanism for sustaining the discharge is emission of secondary
electrons from the cathode. An electron is emitted from the cathode surface with a
specified probability when struck by an ion. These electrons are then accelerated by the
strong electric field close to the cathode where they acquire enough energy to initiate
ionization. The net result is a rapid increase in the electron density close to the cathode in
a region often known as the cathode fall or Crookes dark space.
n e = --- e, th n e
1
2 p ( p n ) (1)
p
n = --- e, th n
5
6 p p ( p n ) (2)
p
PLASMA CHEMISTRY

In this discharge, the electron density and density of excited species is relatively low so
1 Ars=>Ar 1
2 Ar+=>Ar 1

The electric potential, electron density and mean electron energy are all quantities of
interest. Most of the variation in each of these quantities occurs along the axial length of
the column. Figure 2 plots the electron density in the column. The electron density peaks
in the region between the cathode fall and positive column. This region is sometimes
referred to as Faraday dark space. The electron density also decreases rapidly in the radial
direction. The is caused by diffusive loss of electrons to the outer walls of the column
where they accumulate a surface charge. The build up of negative charge leads to a positive
potential in the center of the column with respect to the walls.
In Figure 4 the electric potential is plotted along the axial length of the column. Notice
that the potential profile is markedly different from the linear drop in potential which
results in the absence of the plasma. The strong electric field in the cathode region can lead
to high energy ion bombardment of the cathode. Heating of the cathode surface occurs
which may in turn lead to thermal electron emission where additional electrons are emitted
from the cathode surface.
Figure 3: Plot of electron temperature along the axial length of the positive column.
Figure 4: Plot of the electric potential along the axial length of the positive column.
Figure 5: Plot of electric potential along the axial length of the positive column.
Figure 6: Plot of the electron temperature along the axial length of the positive column.
Figure 7: Plot of the electron density along the axial length of the positive column.
Figure 8: Plot of the number density of excited argon atoms.
Figure 9: Plot of the number density of argon ions.
Reference

positive_column_2d
NEW
MODEL WIZARD
3 Click Add.
4 Click Study.
6 Click Done.
GEOMETRY 1
Rectangle 1 (r1)
Rectangle 2 (r2)
4 In the Height text field, type 6e-3.
Rectangle 3 (r3)
4 In the Height text field, type 6e-3.
Compose 1 (co1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Compose.
2 Select the object r1 only.
3 In the Settings window for Compose, locate the Compose section.
4 In the Set formula text field, type r1-r2-r3.
Mesh Control Edges 1 (mce1)

1 On the Geometry toolbar, click Virtual Operations and choose Mesh Control Edges.
2 On the object fin, select Boundary 7 only.
3 On the Geometry toolbar, click Build All.
DEFINITIONS
Variables 1

mobility
V0 125[V] V Applied voltage
Wf 5 Work function
p0 0.5[torr] Pa Gas pressure
Explicit 1
3 In the Rename Explicit dialog box, type Cathode in the New label text field.
4 Click OK.
Explicit 2
3 In the Rename Explicit dialog box, type Anode in the New label text field.
4 Click OK.
Explicit 3
4 Click OK.
Explicit 4
3 In the Rename Explicit dialog box, type All Walls in the New label text field.
4 Click OK.
Explicit 5
3 In the Rename Explicit dialog box, type Non Cathode Walls in the New label text field.
4 Click OK.
7 Select Boundaries 2, 6, 8, 9, 11, and 12 only.
PLASMA (PLAS)

3 Click Browse.
Ar_xsecs.txt.

Plasma Model 1
Model 1.
You change the way the source coefficients for electronic excitation and ionization are
specified. By default, COMSOL computes rate coefficients based on the cross section
data you supplied. For DC discharges, Townsend coefficients provide a more accurate
description of the cathode fall region so they should be used. The Townsend coefficients
are typically computed using the Boltzmann Equation, Two-Term Approximation interface.
2: e+Ar=>e+Ars
e+Ars.
town2.txt.
4: e+Ar=>2e+Ar+
2e+Ar+.
town4.txt.
Reaction 1
Reaction 2
When solving a reacting flow problem there always needs to be one species which is
selected to fullfill the mass constraint. This should be taken as the species with the
largest mass fraction.
Species: Ar
Species: Ars
Ars.
electroneutrality constraint is satisfied.
Species: Ar+
Ar+.
Wall 1
3 From the Selection list, choose All Walls.
Ground 1
3 From the Selection list, choose Cathode.
Metal Contact 1
2 In the Settings window for Metal Contact, locate the Boundary Selection section.
3 From the Selection list, choose Anode.
4 Locate the Terminal section. In the V0 text field, type V0.
5 Locate the Quick Circuit Settings section. From the Quick circuit type list, choose Series
RC circuit.
Dielectric Contact 1
1 On the Physics toolbar, click Boundaries and choose Dielectric Contact.
2 In the Settings window for Dielectric Contact, locate the Boundary Selection section.
Now you add a surface reaction which describes the neutralization of Argon ions on the
electrode. Secondary emission of electrons is required to sustain the discharge, so you
enter the emission coefficient and an estimate of the mean energy of the secondary
electrons based on the ionization energy threshold and the work function of the surface.
Surface Reaction 1
4 Locate the Boundary Selection section. From the Selection list, choose Cathode.
Make the secondary emission coefficient 0.25 and set the mean energy of the secondary
electrons to be the ionization energy (given by the expression plas.de_4) minus twice
the work function of the electrode.
6 In the i text field, type plas.de_4-2*Wf.
Surface Reaction 2
4 Locate the Boundary Selection section. From the Selection list, choose Non Cathode Walls.
Surface Reaction 3
4 Locate the Boundary Selection section. From the Selection list, choose All Walls.
MESH 1
More Operations>Edge.
Edge 1
Select Boundaries 3, 5, 6, 8, 10, 11, and 14 only.
Size 1
Free Triangular 1
In the Model Builder window, under Component 1 (comp1)>Mesh 1 right-click Free
Triangular 1 and choose Size.
Size 1

section.
type 4.
5 Click Build All.
STUDY 1

Dependent.
4 Click Range.
10 Click Add.
RESULTS

In the Model Builder window, expand the Data Sets node, then click Study 1/Solution 1
(sol1).
Selection
Cut Line 2D 1
1 On the Results toolbar, click Cut Line 2D.
2 In the Settings window for Cut Line 2D, locate the Line Data section.
3 In row Point 1, set z to 0.016.
4 In row Point 2, set r to 0 and z to 0.384.
Mirror 2D 1
On the Results toolbar, click More Data Sets and choose Mirror 2D.

1 In the Model Builder window, under Results click Electron Density (plas).
3 From the Data set list, choose Mirror 2D 1.
4 On the Electron Density (plas) toolbar, click Plot.

1 In the Model Builder window, under Results click Electron Temperature (plas).
4 On the Electron Temperature (plas) toolbar, click Plot.

4 On the Electric Potential (plas) toolbar, click Plot.
1D Plot Group 4
2 In the Settings window for 1D Plot Group, type Electric Potential on Axis in the
Label text field.
5 In the associated text field, type Distance (x).
7 In the associated text field, type Electric Potential (V).
8 Locate the Data section. From the Data set list, choose Cut Line 2D 1.
Line Graph 1
1 On the Electric Potential on Axis toolbar, click Line Graph.
of the y-axis data section. From the menu, choose Component 1>Plasma (Electrostatics)>
Electric>V - Electric potential.
Electric Potential on Axis

1 In the Model Builder window, under Results click Electric Potential on Axis.
2 On the Electric Potential on Axis toolbar, click Plot.
1D Plot Group 5
2 In the Settings window for 1D Plot Group, type Electron Temperature on Axis in
7 In the associated text field, type Electron Temperature (V).
Line Graph 1
1 On the Electron Temperature on Axis toolbar, click Line Graph.
Electron energy density>plas.Te - Electron temperature.
3 On the Electron Temperature on Axis toolbar, click Plot.
1D Plot Group 6
2 In the Settings window for 1D Plot Group, type Electron Density on Axis in the
Label text field.
7 In the associated text field, type Electron Density (1/m^3).
Line Graph 1
1 On the Electron Density on Axis toolbar, click Line Graph.
Electron density>plas.ne - Electron density.
3 On the Electron Density on Axis toolbar, click Plot.
1D Plot Group 7
2 In the Settings window for 1D Plot Group, type Excited Argon Number Density on
Axis in the Label text field.
7 In the associated text field, type Excited Argon Density (1/m^3).
Line Graph 1
1 On the Excited Argon Number Density on Axis toolbar, click Line Graph.
Transport)>Number densities>plas.n_wArs - Number density.
3 On the Excited Argon Number Density on Axis toolbar, click Plot.
1D Plot Group 8
2 In the Settings window for 1D Plot Group, type Argon Ion Number Density on Axis
7 In the associated text field, type Argon +1 Density (1/m^3).
Line Graph 1
1 On the Argon Ion Number Density on Axis toolbar, click Line Graph.
3 On the Argon Ion Number Density on Axis toolbar, click Plot.
Surface Chemistry Tutorial
Introduction
Surface chemistry is often the most important and most overlooked aspect of reacting flow
modeling. Surface rate expressions can be hard to find or not even exist at all. Often it is
preferable to use sticking coefficients to describe surface reactions because they can be
estimated intuitively.
Model Definition
The tutorial model simulates outgassing from a wafer during a chemical vapor deposition
(CVD) process. Careful attention is paid to the overall mass balance in the system and the
difference between the mass averaged velocity and diffusion velocity is explored. The
geometry and operating principle for the model is shown in Figure 1. Initially the closed
container is full of 99.9% silicon hydride and 0.1% hydrogen (measured by molar content).
SiH4 converted to Si and H2 H2 outgassing
Wafer (in red)
Figure 1: Geometry and basic operating principle for the surface chemistry tutorial model.
A surface reaction begins to occur on the wafer which consumes the silicon hydride,
releases hydrogen into the domain, alters the composition of the absorbed species on the
2 | S U R F A C E C H E M I S T R Y TU T O R I A L
wafer surface and deposits bulk silicon. The following reactions are considered on the
surface of the wafer:
TABLE 0-1: SURFACE REACTIONS CONSIDERED
REACTION STICKING COEFFICIENT
SiH4+2Si(s)=>Si(b)+2SiH(s)+H2 1E-4
SiH(s)=>Si(s)+0.5H2 1E-4
The (s) here denotes surface species which means that the species only exists on the
surfaces where the reaction is occurring. To indicate that a species is bulk, append (b) to
the end of the species. Since bulk species cannot participate in surface reactions, they must
only be products and not reactants in a surface reaction. The net result of these two
competing reactions is SiH4=>Si(b)+2H2, which means it is expected that the silane is
replaced by hydrogen inside the reactor, and layers of silicon deposited on the wafer
surface.
MODEL EQUATIONS - SURFACE REACTIONS AND SURFACE SPECIES

The surface reaction rate for reaction i is given by:
ck
f
q i = k f, i ki
k=1
where the rate constant, kf is given by:
i
k f, i = -------------------- ------------------
1 1 8RT
- --- ------------
1 i 2 ( ) m 4 M
tot n
and ck is the concentration of species k which may be a volumetric or surface species, m is

the reaction order minus 1, T is the surface temperature, R is the gas constant and Mn is
the mean molecular weight of the gas mixture. i is the dimensionless sticking coefficient.
For the surface species the following equations are solved:
N
d
dt
= qi i
i=1
where is the surface site concentration (mol/m2), qi is the reaction rate for reaction i
(mol/m2), i is the change in site occupancy number for reaction i (dimensionless).
dZ k R
= ----------
dt tot
where tot is the total surface site concentration (mol/m2), Zk is the site fraction
(dimensionless) and R is the surface rate expression (mol/m2).
For the bulk surface species, the following equation is solved for the deposition height:
dh RM w
= --------------
dt
where h is the total growth height (m), M is the molecular weight (kg/mol) and is the
density of the bulk species (kg/m3).
DOMAIN EQUATIONS
Inside the domain, the Navier-Stokes equations are solved for the fluid velocity. The mass
fraction of hydrogen is computed by solving:
( w ) + ( u )w = j
t
where w is the mass fraction of hydrogen. For detailed information on the transport of the
non-electron species see Theory for the Heavy Species Transport Interface in the Plasma
Module Users Guide. The mass fraction of silane is not directly computed. Its value
comes from the fact that the sum of the mass fractions must equal one. The gas
temperature is computed by solving the energy equation.
BOUNDARY CONDITIONS
Fluid Flow
The following boundary conditions are used. The mass averaged velocity is constrained
using:
Mf
u = ------- n

where Mf is the inward (or outward in this example) mass flux which is defined, from the
surface chemistry as:
Kg

Mf = Mk sk
k=1
where sk is the surface rate expression for each species which comes from summing the
surface reaction rates multiplied by their stoichiometric coefficients over all surface
reactions:
I

sk = ki qi
i=1
Hydrogen Mass Fraction

The flux of hydrogen at the surface of the wafer comes from the surface reactions:
n j = Mk sk
Energy Equation
For the energy equation, the following boundary condition is used:
n T = qi hi
i=1
where hi is the molar enthalpy change due to reaction i.

The y-component of the mass averaged velocity is plotted in Figure 2. The mass averaged
velocity is negative at the surface of the wafer. This means that overall, mass is leaving the
system. This is to be expected since silane is being consumed, which has a molecular weight
of 0.032 kg/mol and replacing it with hydrogen, which has a lower molecular weight of
0.002 kg/mol.
Figure 2: Plot of the mass averaged velocity field after 5 seconds.
The surface reactions begin to consume the silane which results in concentration gradients
within the reactor. The outward mass flux at the wafer surface leads to a mass averaged
velocity everywhere inside the reactor. The combination of concentration gradients and
convection due to the mass averaged velocity tends to draw silane to wards the wafer. The
first surface reaction is exothermic while the second surface reaction is endothermic. The
amount of heat released by the exothermic reaction dominates, so the temperature begins
to increase. The temperature is plotted in Figure 3 and is highest on the wafer surface.
Figure 3: Plot of the change in gas temperature after 5 seconds. The higher temperature is
observed at the wafer surface due to the heat released in the surface reactions.
From the net reaction SiH4=>Si(b)+2H2 there should be a molar inflow of hydrogen at
twice the rate that silane is leaving. Although this condition is not applied explicitly, it is
implicit from the equations solved in the Heavy Species transport interface. Despite the
fact that there is a negative mass averaged velocity, there is a positive diffusion velocity for
the hydrogen at the wafer surface. The y-component of the diffusion velocity for hydrogen
is plotted in Figure 4.
Figure 4: Plot of the y-component of the diffusion velocity for hydrogen after 5 seconds. Even
though the mass averaged velocity is negative, the diffusion velocity is positive indicating an
inflow of hydrogen and a net outflow of total mass.
Figure 5 plots the integrated pressure in the reactor divided by the initial average pressure
in the reactor. Once the problem reaches steady state the average pressure in the reactor
has increased by a factor of 2. This is expected since every mole of silane is being replaced
by two moles of hydrogen. The silane in the reactor is completely consumed after around
200 seconds.
Figure 5: Plot of the ratio of the average pressure at time t divided by the average initial
pressure.
In Figure 6 the total mass initially present in the reactor divided by the total mass in the
reactor is plotted. The total mass in the reactor drops by a factor of 8, which is expected
since the silane with molecular mass 0.032 kg/mol is replaced with two moles of hydrogen
which has a molecular weight of 20.002=0.004 kg/mol.
Figure 6: Plot of the total mass in the reactor at t = 0 divided by the total mass at t = t.
Another way of verifying the correctness of the model is to compare the mass lost inside
the reactor to the mass accumulated on the wafer surface. The two quantities should be
equal. The total gain in mass from surface and bulk species is compared to the total loss in
mass from the reactor in Figure 7. The two curves agree very well indicating that the total
mass in the entire system is conserved.
Figure 7: Plot of the total mass lost in the domain (blue) and the total mass gained due to
Silicon deposition on the wafer surface (green).
The ultimate goal of most CVD models is to determine the total growth height and
growth rate on the surface of the wafer. The total growth height is plotted in Figure 8 and
saturates at about 159 Angstroms. The total growth rate saturates because all the silane in
the reactor is consumed after around 200 seconds. The accumulated growth height is also
very uniform across the surface of the wafer.
Figure 8: Plot of the growth height (surface and z-axis) vs. the wafer arc length (x-axis) and
time (y-axis). The final height of deposited silicon is 158 .
Reference
1. R.J. Kee, M.E. Coltrin, and P. Glarborg, Chemically Reacting Flow Theory and
Practice, John Wiley & Sons, 2003.
Application Library path: Plasma_Module/Chemical_Vapor_Deposition/

surface_chemistry_tutorial
NEW
MODEL WIZARD
2 In the Select Physics tree, select Plasma>Heavy Species Transport (hs).
3 Click Add.
4 In the Select Physics tree, select Fluid Flow>Single-Phase Flow>Laminar Flow (spf).
5 Click Add.
6 In the Select Physics tree, select Heat Transfer>Heat Transfer in Fluids (ht).
7 Click Add.
8 Click Study.
9 In the Select Study tree, select Preset Studies for Selected Physics Interfaces>Time
Dependent.
10 Click Done.
the menu.
the menu.
GEOMETRY 1
Square 1 (sq1)
1 On the Geometry toolbar, click Primitives and choose Square.
2 In the Settings window for Square, locate the Size section.
3 In the Side length text field, type 0.1.
Polygon 1 (pol1)
1 On the Geometry toolbar, click Primitives and choose Polygon.
2 In the Settings window for Polygon, locate the Coordinates section.
3 In the x text field, type 0.03 0.07.
4 In the y text field, type 0 0.
DEFINITIONS
Variables 1

mass_domain at(0, Mass change in the
intop1(hs.rho))-int domain
op1(hs.rho)
mass_bulk intop2(h_Si_bulk* Total mass of bulk
2329[kg/m^3]) species
mass_Si_surf intop2(Zk_Si_surf* Total mass of Si
gamman_Si_surf/ surface species
hs.sigmak_Si_surf*
0.028[kg/mol])
mass_SiH_surf intop2(Zk_SiH_surf* Total mass of SiH
gamman_SiH_surf/ surface species
hs.sigmak_SiH_surf*
0.029[kg/mol])
dmass_Si_surf mass_Si_surf-at(0, Total mass minus
mass_Si_surf) initial mass of Si
surface species
dmass_SiH_surf mass_SiH_surf-at(0, Total mass minus
mass_SiH_surf) initial mass of SiH
surface species
mass_surf dmass_Si_surf+ Total mass minus
dmass_SiH_surf initial mass of all
surface species
2 In the Settings window for Integration, locate the Source Selection section.
H E AV Y S P E C I E S TR A N S P O R T ( H S )
1 In the Model Builder window, under Component 1 (comp1) click Heavy Species Transport
(hs).
2 In the Settings window for Heavy Species Transport, locate the Transport Settings
section.
3 Select the Calculate thermodynamic properties check box.
4 Clear the Migration in electric field check box.
Surface Reaction 1
3 In the Formula text field, type SiH4+2Si(s)=>Si(b)+2SiH(s)+H2.
5 Locate the Reaction Parameters section. In the f text field, type 1E-4.
Surface Reaction 2
3 In the Formula text field, type SiH(s)=>Si(s)+0.5H2.
5 Locate the Reaction Parameters section. In the f text field, type 1E-4.
Species: SiH4
1 In the Model Builder window, under Component 1 (comp1)>Heavy Species Transport (hs)
click Species: SiH4.
4 Locate the General Parameters section. From the Preset species data list, choose SiH4.
Species: H2
click Species: H2.
3 From the Preset species data list, choose H2.
Species: Si(s)
click Species: Si(s).
2 In the Settings window for Surface Species, locate the General Parameters section.
3 From the Preset species data list, choose Si.
4 Locate the Surface Species Parameters section. In the Zk, 0 text field, type 0.995.
Species: Si(b)
click Species: Si(b).
3 From the Preset species data list, choose Si.
Species: SiH(s)
click Species: SiH(s).
3 From the Preset species data list, choose SiH.
4 Locate the Surface Species Parameters section. In the Zk, 0 text field, type 0.005.
Convection and Diffusion

click Convection and Diffusion.
2 In the Settings window for Convection and Diffusion, locate the Model Inputs section.
3 From the u list, choose Velocity field (spf).
4 From the T list, choose Temperature (ht).
5 From the pA list, choose Absolute pressure (spf).
LAMINAR FLOW (SPF)

On the Physics toolbar, click Heavy Species Transport (hs) and choose Laminar Flow (spf).
Fluid Properties 1
1 In the Model Builder window, expand the Component 1 (comp1)>Laminar Flow (spf) node,
then click Fluid Properties 1.
2 In the Settings window for Fluid Properties, locate the Fluid Properties section.
3 From the list, choose Density (hs/cdm1).
4 From the list, choose Dynamic viscosity (hs/cdm1).
Therefore set the flow to be compressible.
5 In the Model Builder window, click Laminar Flow (spf).
box.
9 Click to expand the Discretization section. From the Discretization of fluids list, choose
P2+P1.
Define the pressure reference level in the interface properties.
10 Locate the Physical Model section. Find the Reference values subsection. In the pref text
field, type 13.3.
Inlet 1
3 In the Settings window for Inlet, locate the Boundary Condition section.
4 From the list, choose Mass flow.
5 Locate the Mass Flow section. From the Mass flow type list, choose Pointwise mass flux.
6 From the Mf list, choose Inward mass flux (hs).
Reaction forces must apply to the velocities only in order to avoid undesired couplings
to the surface reactions, so use a unidirectional constraint.
8 Click to expand the Constraint settings section. Locate the Constraint Settings section.
From the Apply reaction terms on list, choose Individual dependent variables.
Fluid 1
1 In the Model Builder window, expand the Component 1 (comp1)>Heat Transfer in Fluids
(ht) node, then click Fluid 1.
2 In the Settings window for Fluid, locate the Model Inputs section.
5 Locate the Heat Conduction, Fluid section. From the k list, choose Thermal conductivity
(hs/cdm1).
6 Locate the Thermodynamics, Fluid section. From the Fluid type list, choose Ideal gas.
7 From the Gas constant type list, choose Mean molar mass.
8 From the Mn list, choose Mean molar mass (hs/cdm1).
9 From the Cp list, choose Mass-averaged mixture specific heat (hs/cdm1).
10 In the Model Builder window, click Heat Transfer in Fluids (ht).
Temperature 1
Initial Values 1
2 In the Settings window for Initial Values, type 300 in the T text field.
Boundary Heat Source 1

1 On the Physics toolbar, click Boundaries and choose Boundary Heat Source.
3 In the Settings window for Boundary Heat Source, locate the Boundary Heat Source
section.
4 From the Qb list, choose Total surface heat source of reaction (hs).
MESH 1
Size.
Size 1
1 In the Settings window for Size, locate the Geometric Entity Selection section.
2 From the Geometric entity level list, choose Point.
3 Select Points 3 and 4 only.
4 Locate the Element Size section. From the Predefined list, choose Extremely fine.
Edge 1
Size 1
Free Triangular 1
STUDY 1

Dependent.
3 In the Times text field, type range(0,5,300).
Solution 1 (sol1)
4 In the Model Builder window, under Study 1>Solver Configurations>Solution 1 (sol1) click
Time-Dependent Solver 1.
5 In the Settings window for Time-Dependent Solver, click to expand the Absolute
tolerance section.
6 Locate the Absolute Tolerance section. From the Global method list, choose Unscaled.
7 Click to expand the Time stepping section. Locate the Time Stepping section. In the Initial
step text field, type 0.001.
RESULTS
Velocity (spf)
1 In the Model Builder window, under Results click Velocity (spf).
3 From the Time (s) list, choose 5.
Surface
1 In the Model Builder window, expand the Velocity (spf) node, then click Surface.
the Expression section. From the menu, choose Component 1>Laminar Flow>Velocity and
pressure>Velocity field>v - Velocity field, y component.
3 On the Velocity (spf) toolbar, click Plot.
Temperature (ht)
1 In the Model Builder window, under Results click Temperature (ht).
3 From the Time (s) list, choose 5.
Surface
1 In the Model Builder window, expand the Temperature (ht) node, then click Surface.
3 In the Expression text field, type T-300[K].
4 On the Temperature (ht) toolbar, click Plot.
2D Plot Group 6
2 In the Settings window for 2D Plot Group, type Diffusion Velocity, Hydrogen in
3 Locate the Data section. From the Time (s) list, choose 5.
Surface 1
1 Right-click Diffusion Velocity, Hydrogen and choose Surface.
the Expression section. From the menu, choose Component 1>Heavy Species Transport>
Species>Species wH2>Diffusion velocity>hs.Vdy_wH2 - Diffusion velocity, y component.
3 On the Diffusion Velocity, Hydrogen toolbar, click Plot.
1D Plot Group 7
2 In the Settings window for 1D Plot Group, type Average Pressure Divided by
Average Initial Pressure in the Label text field.
Global 1
1 On the Average Pressure Divided by Average Initial Pressure toolbar, click Global.
2 In the Settings window for Global, click to expand the Legends section.
3 Find the Include subsection. Select the Expression check box.
5 Locate the y-Axis Data section. In the table, enter the following settings:

intop1(spf.pA)/at(0, 1 Average pressure divided by
intop1(spf.pA)) average initial pressure
6 On the Average Pressure Divided by Average Initial Pressure toolbar, click Plot.
1D Plot Group 8
2 In the Settings window for 1D Plot Group, type Density Ratio in the Label text field.
3 Locate the Legend section. From the Position list, choose Lower right.
Global 1
1 On the Density Ratio toolbar, click Global.
2 In the Settings window for Global, locate the Legends section.

at(0,intop1(hs.rho))/ 1 Initial mass in the reactor
intop1(hs.rho) divided by the mass in the
reactor
6 On the Density Ratio toolbar, click Plot.
1D Plot Group 9
2 In the Settings window for 1D Plot Group, type Total Reactor Mass in the Label text
field.
3 Locate the Legend section. From the Position list, choose Lower right.
Global 1
1 On the Total Reactor Mass toolbar, click Global.
2 In the Settings window for Global, locate the Legends section.

mass_domain kg/m Total mass in the domain
mass_surf+mass_bulk Total mass of surface plus bulk
species
6 On the Total Reactor Mass toolbar, click Plot.
Edge 2D 1
1 On the Results toolbar, click More Data Sets and choose Edge 2D.
On the Results toolbar, click More Data Sets and choose Parametric Extrusion 1D.
2D Plot Group 10
2 In the Settings window for 2D Plot Group, type Bulk Silicon Growth Height in the
Label text field.
3 Locate the Data section. From the Data set list, choose Parametric Extrusion 1D 1.
Surface 1
1 Right-click Bulk Silicon Growth Height and choose Surface.
the Expression section. From the menu, choose Component 1>Heavy Species Transport>
Etching and deposition>Si bulk>h_Si_bulk - Accumulated growth height.
3 Locate the Expression section. From the Unit list, choose .
4 Right-click Results>Bulk Silicon Growth Height>Surface 1 and choose Height Expression.
Height Expression 1
On the Bulk Silicon Growth Height toolbar, click Plot.
Thermal Plasma
Introduction
Low pressure discharges are characterized by the fact that the electron temperature is
much higher than the neutral gas temperature. As the gas pressure increases, the number
of collisions between the electrons and neutrals increases. At high enough pressures the
electron temperature becomes equal to the gas temperature. At this point the plasma is in
local thermodynamic equilibrium and a much simpler MHD model can be used to model
the plasma.
This model simulates a plasma at medium pressure (2 torr) where the gas temperature
cannot be assumed to be constant but the plasma is still not in local thermodynamic
equilibrium. In Figure 1 the electron (blue) and gas (black) temperatures are plotted as a
function of pressure. At low pressures the two temperatures are decoupled but as the
pressure increases the temperatures tend towards the same limit. There are no axes on the
plot since the exact temperature and pressure depends strongly on the gas in question.
Molecular gases tend to heat more easily than atomic gases for to the following reasons:
Continual dissociation and recombination of the molecule. An electron impact reaction

may dissociate a molecule into its atomic components. The energy lost by the electron
when this reaction occurs is given back to the gas in the form of thermal energy during
recombination.
Vibrational excitation and relaxation. The threshold energy for vibrational excitation is
much smaller than the threshold energy for dissociation. The continuous vibrational
excitation and relaxation of the molecules can cause the gas temperature to increase.
Note: This application requires the Plasma Module and AC/DC Module.
2 | THERMAL PLASMA
Figure 1: Plot of electron (blue) and gas (black) temperature v.s. pressure. At higher pressures
the two temperatures become equal.
Model Definition
n
( ) + [ n e ( e E ) D e n e ] = R e
t e
n
( ) + [ n ( E ) D n ] + E e = R
t
3 | THERMAL PLASMA
D e = e T e, = --- e, D = T e
5
3
Re = x j kj Nn ne
j=1
R = xj kj Nn ne j
j=1

kk = 0 k ( )f ( ) d
( w k ) + ( u )w k = j k + R k
t
0 r V =
4 | THERMAL PLASMA
N
= q Z k n k n e
k = 1
PLASMA CHEMISTRY

electron. Reaction number 1, elastic collisions with electrons is primarily responsible for
heating of the gas. In addition to volumetric reactions, the following surface reactions are
implemented:
1 Ars=>Ar 1
2 Ar+=>Ar 1
5 | THERMAL PLASMA
The reactor geometry is simply a cylindrical glass tube with a 4 turn coil wrapped around
it. Gas flows in from the bottom and exits out of the top. The gas is heated through elastic
and inelastic collisions. The inelastic collisions are responsible for the bulk of the gas
heating. A fixed power of 700 W is applied to the coil.
Coil
Inlet
Figure 2: Schematic of the ICP reactor. Flow enters from the base and leaves out the top.
6 | THERMAL PLASMA
7 | THERMAL PLASMA
Figure 4: Plot of the electron temperature inside the plasma source.
Figure 5: Plot of the plasma potential inside the plasma source.
8 | THERMAL PLASMA
Figure 6: Plot of the velocity field.
Figure 7: Revolved plot of the temperature inside the plasma source.
9 | THERMAL PLASMA
Figure 8: Plot of the mass fraction of ground state argon.
Figure 9: Plot of the mass fraction of electronically excited argon atoms.
10 | THERMAL PLASMA
Reference

thermal_plasma
NEW
MODEL WIZARD
3 Click Add.
4 In the Select Physics tree, select Fluid Flow>Single-Phase Flow>Laminar Flow (spf).
5 Click Add.
6 In the Select Physics tree, select Heat Transfer>Heat Transfer in Fluids (ht).
7 Click Add.
8 Click Study.
10 Click Done.
the menu.
11 | THERMAL PLASMA
GEOMETRY 1
Import 1 (imp1)
1 On the Home toolbar, click Import.
2 In the Settings window for Import, locate the Import section.
3 Click Browse.
thermal_plasma.mphbin.
5 Click Import.
GLOBAL DEFINITIONS
Parameters

p0 1[torr] 133.3 Pa Initial and outlet pressure
DEFINITIONS
Explicit 1
4 In the Rename Explicit dialog box, type Plasma in the New label text field.
5 Click OK.
Explicit 2
4 Select Boundaries 9, 10, 34, and 35 only.
7 Click OK.
12 | THERMAL PLASMA
Explicit 3
6 In the Rename Explicit dialog box, type Outlet in the New label text field.
7 Click OK.
Explicit 4
4 Select Boundaries 12, 13, 1517, 19, 21, 22, 2426, and 2832 only.
6 In the Rename Explicit dialog box, type Coil Walls in the New label text field.
7 Click OK.
Start by importing the cross sections for argon and by activating the convection and
thermodynamic property evaluation.
PLASMA (PLAS)
2 From the Selection list, choose Plasma.

3 Click Browse.
Ar_xsecs.txt.

8 Select the Calculate thermodynamic properties check box.
13 | THERMAL PLASMA
check box.
Reaction 1
Reaction 2
Ground 1
Surface Reaction 1
2 In the Settings window for Surface Reaction, locate the Boundary Selection section.
4 Locate the Reaction Formula section. In the Formula text field, type Ar+=>Ar.
Surface Reaction 2
2 In the Settings window for Surface Reaction, locate the Boundary Selection section.
4 Locate the Reaction Formula section. In the Formula text field, type Ars=>Ar.
Wall 1
14 | THERMAL PLASMA
Species: Ars
In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click Species: Ars.
Outflow 1
3 From the Selection list, choose Outlet.
Species: Ar+
In the Model Builder window, under Component 1 (comp1)>Plasma (plas) click Species: Ar+.
Outflow 1
3 From the Selection list, choose Outlet.
Species: Ar
Species: Ars
Ars.
5 Click to expand the Species thermodynamic parameters section. Locate the Species
Thermodynamic Parameters section. In the h text field, type 11.5.
The thermodynamic properties for the electronically excited Argon atoms can be the
same as for the ground state species plus the threshold energy for the electron impact
reaction. In this case this corresponds to an energy of 11.5eV. This is added in the text
field Additional enthalpy contribution.
Species: Ar+
Ar+.
15 | THERMAL PLASMA
5 Click to expand the Species thermodynamic parameters section. Locate the Species
Thermodynamic Parameters section. In the h text field, type 15.8.
The thermodynamic properties for the Argon ions can be the same as for the ground
state species plus the threshold energy for ionization. In this case this corresponds to an
energy of 15.8eV. This is added in the text field Additional enthalpy contribution.
You can set the gas temperature and pressure in the plasma model to the computed gas
pressure and temperature from other physics interfaces. The velocity field is also set to
the velocity field computed from the Laminar Flow interface.
Plasma Model 1
Model 1.
4 From the T list, choose Temperature (ht).
6 Locate the Electron Density and Energy section. In the eNn text field, type 4E24.
Initial Values 1
Values 1.

In the Model Builder window, under Component 1 (comp1) click Magnetic Fields (mf).
Coil 1
3 In the Settings window for Coil, locate the Coil section.
16 | THERMAL PLASMA
5 Select the Coil group check box.
6 In the Pcoil text field, type 700[W].
LAMINAR FLOW (SPF)

Therefore set the flow to be weakly compressible.
1 In the Model Builder window, under Component 1 (comp1) click Laminar Flow (spf).
box.
Define the pressure reference level to be zero. This means that absolute pressure values
are used throughout the model.
5 Locate the Physical Model section. Find the Reference values subsection. In the pref text
field, type 0.
7 In the Settings window for Laminar Flow, click to expand the Consistent stabilization
section.
8 Locate the Consistent Stabilization section. Find the Navier-Stokes equations subsection.
Clear the Crosswind diffusion check box.
9 Locate the Domain Selection section. Click Clear Selection.
11 Click to expand the Equation section. From the Equation form list, choose Stationary.
Fluid Properties 1
1 In the Model Builder window, expand the Laminar Flow (spf) node, then click Fluid
Properties 1.
2 In the Settings window for Fluid Properties, locate the Fluid Properties section.
3 From the list, choose Density (plas/pes1).
4 From the list, choose Dynamic viscosity (plas/pes1).
17 | THERMAL PLASMA
Initial Values 1
1 In the Model Builder window, under Component 1 (comp1)>Laminar Flow (spf) click Initial
Values 1.
3 In the p text field, type p0.
Inlet 1
3 In the Settings window for Inlet, locate the Boundary Condition section.
4 From the list, choose Mass flow.
5 Locate the Mass Flow section. From the Mass flow type list, choose Standard flow rate
(SCCM).
6 From the Mn list, choose Mean molar mass (plas/pes1).
7 In the Qsccm text field, type 100*tanh(1E5*t[1/s]).
Outlet 1
1 On the Physics toolbar, click Boundaries and choose Outlet.
3 In the Settings window for Outlet, locate the Pressure Conditions section.
4 In the p0 text field, type p0.
1 In the Model Builder window, under Component 1 (comp1) click Heat Transfer in Fluids
(ht).
2 In the Settings window for Heat Transfer in Fluids, locate the Domain Selection section.
5 Click to expand the Equation section. From the Equation form list, choose Time
dependent.
Fluid 1
1 In the Model Builder window, expand the Heat Transfer in Fluids (ht) node, then click Fluid
1.
2 In the Settings window for Fluid, locate the Model Inputs section.
18 | THERMAL PLASMA
5 Locate the Heat Conduction, Fluid section. From the k list, choose Thermal conductivity
(plas/pes1).
6 Locate the Thermodynamics, Fluid section. From the list, choose Density (plas/pes1).
7 From the Cp list, choose Mass-averaged mixture specific heat (plas/pes1).
8 From the list, choose User defined.
Initial Values 1
2 In the Settings window for Initial Values, type 300 in the T text field.
3 In the Model Builder window, click Heat Transfer in Fluids (ht).
Heat Source 1
1 On the Physics toolbar, click Domains and choose Heat Source.
3 In the Settings window for Heat Source, locate the Heat Source section.
4 From the Q0 list, choose Heat source for gas (plas/pes1).
Temperature 1
Outflow 1
1 On the Physics toolbar, click Boundaries and choose Outflow.
MATERIALS
Blank Material.
Material 1 (mat1)
19 | THERMAL PLASMA

group

6 Click OK.
Material 2 (mat2)

group

6 In the Rename Material dialog box, type Air in the New label text field.
7 Click OK.
Material 3 (mat3)
20 | THERMAL PLASMA

group
Electrical conductivity sigma 6E7 S/m Basic

6 In the Rename Material dialog box, type Copper coil in the New label text field.
7 Click OK.
A boundary layer mesh is used on the reactor walls so that the region of space charge
separation between the ions and electrons can be resolved.
MESH 1
2 From the Element size list, choose Extra fine.
Edge 1
Size 1
4 Click to expand the Element size parameters section. Locate the Element Size section.
Click the Custom button.
Free Triangular 1
21 | THERMAL PLASMA
Size 1
Boundary Layers 1

section.
type 5.
Mapped 1
Distribution 1
1 Right-click Component 1 (comp1)>Mesh 1>Mapped 1 and choose Distribution.
3 From the Selection list, choose Coil Walls.
distribution type.
22 | THERMAL PLASMA
Free Triangular 2
STUDY 1
3 Click Range.
9 Click Add.
RESULTS

2D Plot Group 11
2 In the Settings window for 2D Plot Group, type Argon Mass Fraction in the Label
text field.
Surface 1
1 Right-click Argon Mass Fraction and choose Surface.
Transport)>Mass fractions>plas.wAr - Mass fraction.
23 | THERMAL PLASMA
3 On the Argon Mass Fraction toolbar, click Plot.
2D Plot Group 12
Surface 1
1 Right-click Excited Argon Mass Fraction and choose Surface.
24 | THERMAL PLASMA

Plasma Application Library Manual

Uploaded by

Copyright:

Available Formats

Plasma Application Library Manual

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Plasma Application Library Manual

Uploaded by

Copyright:

Available Formats

Plasma Module

Application Library Manual

Support Center: www.comsol.com/support

Part number: CM021404

Dielectric Barrier Discharge

1. There are typically around 1,000,000 m

2 | DIELECTRIC BARRIER DISCHARGE

Figure 1: Graphic illustration of a typical dielectric barrier discharge.

3 | DIELECTRIC BARRIER DISCHARGE

The electrostatic field is computed using the following equation:

4 | DIELECTRIC BARRIER DISCHARGE

and the electron energy flux:

5 | DIELECTRIC BARRIER DISCHARGE

REACTION FORMULA TYPE (eV)

REACTION FORMULA STICKING COEFFICIENT

Results and Discussion

6 | DIELECTRIC BARRIER DISCHARGE

Figure 2: Mass fraction of excited argon.

7 | DIELECTRIC BARRIER DISCHARGE

8 | DIELECTRIC BARRIER DISCHARGE

Figure 5: Extruded plot of the electron density.

Figure 6: Extruded plot of the mean electron energy.

9 | DIELECTRIC BARRIER DISCHARGE

10 | DIELECTRIC BARRIER DISCHARGE

11 | DIELECTRIC BARRIER DISCHARGE

Figure 9: Plot of the total discharge current vs. time.

12 | DIELECTRIC BARRIER DISCHARGE

13 | DIELECTRIC BARRIER DISCHARGE

Name Expression Unit Description

Cross Section Import 1

5 In the Model Builder window, click Plasma (plas).

14 | DIELECTRIC BARRIER DISCHARGE

15 | DIELECTRIC BARRIER DISCHARGE

16 | DIELECTRIC BARRIER DISCHARGE

Surface Charge Accumulation 1

17 | DIELECTRIC BARRIER DISCHARGE

Property Name Value Unit Property

Step 1: Time Dependent

18 | DIELECTRIC BARRIER DISCHARGE

Excited Argon Mass Fraction

19 | DIELECTRIC BARRIER DISCHARGE

Excited Argon Mass Fraction 1

20 | DIELECTRIC BARRIER DISCHARGE

Mean Electron Energy

Mean Electron Energy 1

Electron Current Density

21 | DIELECTRIC BARRIER DISCHARGE

Argon Ion Current Density

Argon Ion Current Density 1

22 | DIELECTRIC BARRIER DISCHARGE

Expression Unit Description

23 | DIELECTRIC BARRIER DISCHARGE

GEC CCP Reactor, Argon Chemistry, 1D

The electrostatic field is computed using the following equation:

and the electron energy flux:

REACTION FORMULA TYPE (eV)

REACTION FORMULA STICKING COEFFICIENT

Results and Discussion

Figure 4: Plot of the electric field at different phases of the RF cycle.