INTRODUCTION TO

Plasma Module
[Link] Page 2 Saturday, October 21, 2023 10:31 AM

Introduction to the Plasma Module

© 1998–2023 COMSOL
Protected by patents listed on [Link]/patents, and U.S. Patents 7,519,518; 7,596,474; 7,623,991; 8,219,373;
8,457,932; 9,098,106; 9,146,652; 9,323,503; 9,372,673; 9,454,625; 10,019,544; 10,650,177; and 10,776,541. Patents
pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement ([Link]/comsol-license-agreement) and may be used or copied only under the terms of the
license agreement.
COMSOL, the COMSOL logo, COMSOL Multiphysics, COMSOL Desktop, COMSOL Compiler, COMSOL Server,
and LiveLink 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 [Link]/
trademarks.
Version: COMSOL 6.2

Contact Information
Visit the Contact COMSOL page at [Link]/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 [Link]/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
[Link]/support/case. Other useful links include:

• Support Center: [Link]/support

• Product Download: [Link]/product-download
• Product Updates: [Link]/support/updates
• COMSOL Blog: [Link]/blogs
• Discussion Forum: [Link]/community
• Events: [Link]/events
• COMSOL Video Gallery: [Link]/video
• Support Knowledge Base: [Link]/support/knowledgebase

Part number. CM021402

[Link] Page 3 Saturday, October 21, 2023 10:31 AM

Contents

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Applications of Plasmas in Science & Technology . . . . . . . . . . . . . . 5
Plasma Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Physics Guide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Physics Interface List by Space Dimension and Preset Study Type15
Application Specific Interfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Argon Discharge in the GEC Reference Cell . . . . . . . . . . . . . . . 22

|3
[Link] Page 4 Saturday, October 21, 2023 10:31 AM

4|
[Link] Page 5 Saturday, October 21, 2023 10:31 AM

Introduction

The Plasma Module is tailor-made to model and simulate low-temperature plasma

sources and systems. Engineers and scientists use it to gain insight into the physics
of discharges and gauge the performance of existing or potential designs. The
module can perform analysis in all space dimensions — 1D, 2D, and 3D —.
Plasma systems are, by their very nature, complicated systems with a high degree
of nonlinearity. Small changes to the electrical input or plasma chemistry can result
in significant changes in the discharge characteristics.
Low-temperature plasmas represent the amalgamation of fluid mechanics,
reaction engineering, physical kinetics, heat transfer, mass transfer, and
electromagnetics. The Plasma Module is a specialized tool for modeling
nonequilibrium and equilibrium discharges, which occur in a wide range of
engineering disciplines.
The Plasma Module consists of a suite of physics interfaces with which you can
model arbitrary systems, a set of documented examples, and a manual. The
intended audience is researchers and engineers with a background in the physics
of low-temperature discharges.

Applications of Plasmas in Science & Technology

There are many different types of plasma that are typically of interest. The main
difference between the plasmas is that the mechanism of energy transfer between
the electrons and fields is different. In this section the following types of plasma
are discussed:
• Direct Current Discharges
• Inductively coupled plasmas
• Capacitively Coupled Plasmas
• Microwave Plasmas
• Plasma-enhanced Chemical Vapor Deposition
• Dielectric Barrier Discharges
• Equilibrium Discharges

D IRECT C URRENT D ISCHARGES

Direct current (DC) discharges can sustained through secondary electron
emission at the cathode due to ion bombardment. The electrons ejected from the
cathode are accelerated through the cathode fall region into the bulk of the
plasma. They may acquire enough energy to ionize the background gas, creating

|5
[Link] Page 6 Saturday, October 21, 2023 10:31 AM

a new electron-ion pair. The electron makes its way to the anode whereas the ion
will migrate to the cathode, where it may create a new secondary electron.
Several distinct regions typically exist in a DC discharge, depending on the specific
operating conditions and geometry. A typically DC discharge is shown below. The
cathode glow (often known as cathode fall) region is where most of the voltage
drop occurs. This is the region where the discharge current is built up due to the
high degree of ionization. The electron density and flux grow exponentially in this
region (Ref. 5). The negative glow and Faraday space regions connect the cathode
glow to the positive column. In this region the electron temperature drops rapidly
due to the fact that the electric field is very weak. In the positive column region,
the discharge is very uniform, and the electron-ion pairs generated are lost to the
walls in the radial direction due to the ambipolar field (Ref. 5).
Negative
Aston DS glow Faraday space Anode DS

Cathode
(-)
Anode
(+)
Cathode Positive column Anode
glow Cathode DS
(Crookes) glow

Diagram showing the most important regions in a typical DC discharge, after Ref. 5.

DC discharges can easily begin to arc, especially at high pressures. Arcing is

typically undesirable so, in practice, a ballast resistor or series RC circuit is placed
between the power supply and the cathode. As the discharge current at the
cathode begins to increase, the voltage applied at the cathode will begin to drop
due to the presence of the RC circuit. If the ballast resistance is high enough, the
DC discharge will remain smooth and uniform and not begin to arc. The pressure
range for a DC discharge is typically in the range of 10 Pa all the way up to

6|
[Link] Page 7 Saturday, October 21, 2023 10:31 AM

atmospheric pressure and the applied voltage is typically several hundred volts all
the way up to several kV.

Schematic of a positive column DC discharge including an external circuit. The ballast resistor and
blocking capacitor prevent the discharge from arcing.

I NDUCTIVELY COUPLED PLASMAS

Inductively coupled plasmas (ICP) were first used in the 1960s as thermal plasmas
in coating equipment, Ref. 4. These devices operated at pressures on the order of
0.1 atm and produced gas temperatures on the order of 10,000 K. In the 1990s,
ICP became popular in the film processing industry as a way of fabricating large
semiconductor wafers. These plasmas operated in the low-pressure regime, from
0.002–1 Torr and, as a consequence the gas temperature remained close to room
temperature. Low-pressure ICP are attractive because they provide a relatively
uniform plasma density over a large volume. The plasma density is also high,
around 1018 1/m3, which results in a significant ion flux to the surface of the
wafer. Faraday shields are often added to reduce the effect of capacitive coupling
between the plasma and the driving coil. This reduces the effect of ion
bombardment in the dielectric window, which can lead to degradation and

|7
[Link] Page 8 Saturday, October 21, 2023 10:31 AM

contamination of the chamber. Some typical ICP configurations are shown in the
figure below.

RF

Dielectric
Substrate

RF bias

Dielectric RF

Substrate

RF bias

Sketch of some typical ICP configurations used in the semiconductor processing industry, after Ref. 5.

From an electrical point of view, inductively coupled plasmas act like transformers,
with the driving coil acting as the primary and the plasma acting as the secondary.
The current flowing through the coil induces a current in the plasma, causing the
electrons to flow in the opposite direction to the current flowing in the coil. This,
in turn, induces an opposing current back in the coil. The coil and plasma are
strongly coupled. Heating of the electrons only occurs within the skin depth of
the plasma. The higher the plasma density or the applied frequency, the smaller
the volume over which power is deposited from the fields to the electrons. The

8|
[Link] Page 9 Saturday, October 21, 2023 10:31 AM

design and operating frequency of the coil should be such that the plasma is
uniform over the reactor volume.

Inductively Coupled Plasma. The electron temperature is shown as a slice plot and the current density
induced in the plasma is shown as streamlines. The plasma is created and sustained via a rectangular
single spiral coil.

C APACITIVELY C OUPLED P LASMAS

Capacitively coupled plasmas I(t)
C
(CCP) are often used in the
semiconductor industry to
deposit dielectric films on
semiconductor wafers. A
typical CCP configuration is RF
shown to the right. By driving ~
the electrodes with a time
varying voltage, typically at Substrate
frequencies ranging from
100 kHz up to 100 MHz, the
discharge can sustain itself
without requiring secondary
electron emission (Ref. 4).
Typical operating conditions for capacitively coupled plasmas are driving voltages
from 100–1000 V, pressures of 2–200 Pa, electron densities in the region of
1015–1017 1/m3 and electron temperatures on the order of 2–20 eV. The

|9
[Link] Page 10 Saturday, October 21, 2023 10:31 AM

operating frequency is typically 13.56 MHz, which is one of the frequencies

reserved worldwide for the industrial, scientific, and medical uses.

M ICROWAVE P LASMAS
Microwave plasmas, or wave heated discharges, are sustained when electrons can
gain enough energy from an electromagnetic wave as it penetrates into the plasma.
The physics of a microwave plasma is quite different depending on whether the
TE mode (out-of-plane electric field) or the TM mode (in-plane electric field) is
propagating. In the 2D axisymmetric case, the TE mode means that only the
azimuthal component of the electric field is computed and the TM mode means
that the in-plane r and z components of the electric field are computed. In both
cases it is not possible for the electromagnetic wave to penetrate into regions of
the plasma where the electron density exceeds the critical electron density (around
7.6×1016 1/m3 for argon at 2.45 GHz). On this contour corresponding to the
critical electron density, the wave goes from being propagating to evanescent. This
means that all the power from the wave is absorbed in a very small region in space,
which can make the model numerically unstable. The pressure range for
microwave plasmas is very broad. For electron cyclotron resonance (ECR)
plasmas, the pressure can be on the order of 1 Pa or less. For non-ECR plasmas
the pressure typically ranges from 100 Pa up to atmospheric pressure. The power

10 |
[Link] Page 11 Saturday, October 21, 2023 10:31 AM

can range from a few watts all the way up to several kilowatts. Microwave plasmas
are popular due to the cheap availability of microwave power.

Electric field (surface) and critical electron density (white contour) in a cross-configuration microwave
plasma.

P LASMA - ENHANCED C HEMICAL V APOR D EPOSITION

Plasma-enhanced chemical vapor deposition (PECVD) is a discharge operating in
a reactive gas. It is used as an alternative to traditional chemical vapor deposition,
since it can generate the radical species required to initiate the deposition process
very efficiently, without the need for very high gas temperatures. In the case of
silicon deposition in a silane plasma, the electron impact reactions can dissociate
the silane into highly reactive silicon hydrides and hydrogen gas. Equally as
important as the gas phase reactions is the complicated set of chemical reactions

| 11
[Link] Page 12 Saturday, October 21, 2023 10:31 AM

that take place on the surface of the wafer. Typical PECVD reactors are inductively
coupled plasmas, but microwave discharges can also be used.

PECVD model of silicon deposition onto a wafer (surface) from a silane plasma. Silane is fed in from the
shower head at the top and broken into radical species due to the presence of the electrons, which are
sustained through induction heating.

D IELECTRIC B ARRIER D ISCHARGES

Dielectric barrier discharges (DBD) have many uses including UV and EUV
sources, ozone generation, and surface modification. 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. On one of the dielectric plates, a time varying
voltage is applied. The other plate is electrically grounded. As the voltage applied
increases, the electric field increases in the gap between the plates. Any free
electrons in the gap will be 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 time scale. Electrons created via electron impact ionization rush
toward 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

12 |
[Link] Page 13 Saturday, October 21, 2023 10:31 AM

gas filled gap. Surface charge accumulation temporarily terminates the discharge
until the field reverses direction and the process repeats in the opposing direction.

Dielectric Barrier Discharge: Plot of the mass fraction of electronically excited argon atoms in a dielectric
barrier discharge. The x-axis represents space and the y-axis time. The applied frequency is 10 kHz.

Dielectric barrier discharges typically operate in the pressure range of 0.1–3 atm
with an applied voltage of 1–100 kV and a frequency of kHz up to MHz. The
dielectric barriers are usually made or quartz, glass, or a ceramic material.

E QUILIBRIUM D ISCHARGES
Equilibrium discharges (sometimes called thermal plasmas) have 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. The electron temperature is assumed
to be equal to the temperature of the background gas. Equilibrium discharges are
treated differently to the other examples above. The chemical composition of the
plasma is not computed, only the gas temperature and electric fields. These models
allow the thermal loads on surfaces to be computed in a more efficient and stable
way than using a nonequilibrium discharge model.

| 13
[Link] Page 14 Saturday, October 21, 2023 10:31 AM

Plasma Modeling

There are three methods typically employed to model plasmas. All methods
essentially involve solving an appropriate transport equation along with Maxwell’s
equations.

• The Kinetic approach solves for the distribution function for the ions and
electrons in a plasma by solving the Boltzmann equation. The distribution
function can also be computed using a Lagrangian approach, where ordinary
differential equations are solved for the electron position and velocity, as is
the case in particle-in-cell (PIC) codes. The advantage of this approach is
that it allows for the electron energy distribution function to take an
arbitrary shape in energy space and can also show effects that are not
captured by fluid models. Drawbacks include difficulties associated with
implementing arbitrary plasma chemistries and high computational cost.
• The Fluid approximation describes the plasma in terms of macroscopic
quantities by assuming a specific form of the distribution function and
taking velocity moments of the Boltzmann equation. This results in a set of
coupled partial differential equations that can be solved using the finite
element method. The Plasma Module uses the fluid approximation.
Advantages of the fluid approximation are that it is efficient to solve the
resulting equations and straightforward to define arbitrarily complex plasma
chemistries. It is also easy to couple the electron dynamics to the
electromagnetic fields, which are also computed using the finite element
method.
• Hybrid methods treat some of the components of a plasma as a fluid and
others kinetically. Usually, the electrons are treated kinetically by solving the
Boltzmann equation or using a Monte Carlo method, and the ions and other
radical species are computed using a fluid model. In this sense, hybrid
methods offer a compromise between the fidelity of kinetic models and the
speed and convenience of fluid models.
Modeling the interaction between the plasma and an external circuit is an
important part of understanding the overall characteristics of a discharge. This
module has tools to add circuit elements directly to a 1D, 2D, or 3D model, or
import an existing SPICE Netlist into the model. The plasma chemistry is

14 |
[Link] Page 15 Saturday, October 21, 2023 10:31 AM

specified either by loading in sets of collision cross sections from a file or by adding
reactions and species to the Model Builder.

Physics Guide
The complexity of plasma modeling lies in the fact that it combines elements of
reaction engineering, statistical physics, fluid mechanics, physical kinetics, heat
transfer, mass transfer, and electromagnetics. The net result is a true multiphysics
problem involving complicated coupling between the different physics. The
module is designed to simplify the process of setting up a self-consistent model of
a low-temperature plasma.
The physics interfaces include all the necessary tools to model plasma discharges,
beginning with a Boltzmann Equation, Two-Term Approximation solver that
computes the electron transport properties and source coefficients from a set of
electron impact collision cross sections. This interface makes it possible to
determine many of the interesting characteristics of a discharge by providing input
properties such as the electric field and the electron impact reactions that make up
the plasma chemistry, without solving a space-dependent problem.
For space-dependent models, the reactions and species which make up the plasma
chemistry are conveniently managed in the Model Builder. When the fluid velocity
and gas temperature are of interest, there are physics interfaces available for
laminar flow and heat transfer. There are several options available when coupling
the charged species transport to the electromagnetic fields.

Physics Interface List by Space Dimension and Preset Study Type

PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE

DIMENSION

AC/DC

Electrical Circuit cir Not space stationary; frequency

dependent domain; time dependent;
frequency domain;
eigenfrequency
Electrostatics1 es all dimensions stationary; time dependent;
stationary source sweep

| 15
[Link] Page 16 Saturday, October 21, 2023 10:31 AM

PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE

DIMENSION

Fluid Flow

Single-Phase Flow

Laminar Flow1 spf 3D, 2D, 2D stationary; time dependent

axisymmetric
Plasma

Boltzmann Equation, be 0D mean energies; reduced

Two-Term Approximation electric fields; time
dependent
Plasma, Time Periodic ptp 1D, 2D time periodic; time periodic
axisymmetric to time dependent
Plasma plas all dimensions time dependent; stationary

Inductively Coupled — 3D, 2D, 2D frequency-transient;

Plasma2,4 axisymmetric. frequency-stationary
Inductively Coupled Plasma — 2D, 2D Frequency-Time periodic;
with RF Bias2,4 axisymmetric. time periodic to time
dependent
Microwave Plasma3,4 — 3D, 2D, 2D frequency-transient
axisymmetric
Electric Discharges

Corona Discharge4 — all dimensions stationary

Electrical Breakdown ebd 3D, 2D, 2D time dependent

Detection axisymmetric
Equilibrium Discharges

Equilibrium Discharges, — 2D, 2D stationary; time dependent;

Out-of-Plane axisymmetric frequency-transient;
Currents2,4 frequency-stationary

16 |
[Link] Page 17 Saturday, October 21, 2023 10:31 AM

PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE

DIMENSION

Equilibrium Discharges, — 2D, 2D stationary; time dependent;

In-Plane Currents2,4 axisymmetric frequency-transient;
frequency-stationary
Equilibrium — 3D stationary; time dependent;
Discharges2,4 frequency-transient;
frequency-stationary
Species Transport

Drift Diffusion dd all dimensions time dependent

Heavy Species hs all dimensions time dependent

Transport
Charge Transport ct all dimensions stationary

1
This physics interface is included with the core COMSOL package but has added functionality
for this module.
2 Requires the addition of the AC/DC Module.
3 Requires the addition of the RF Module.
4 This physics interface is a predefined multiphysics coupling that automatically adds all the

physics interfaces and coupling features required.

AC/DC I NTERFACES
The AC/DC Branch chapter describes the two physics interfaces available with
this module under the AC/DC branch of the Model Wizard. Use the
Electrostatics interface to compute the electrostatic field in the plasma caused by
separation of space charge between the electrons and ions. The Plasma Module
enhances the interface included with the basic COMSOL license.

F LUID F LOW I NTERFACES

The Fluid Flow Branch describes the Laminar Flow interface, which has a few
additional features available for this module compared to the basic license.

B OLTZMANN E QUATION , T WO -T ERM A PPROXIMATION I NTERFACE

The Boltzmann Equation, Two-Term Approximation interface computes the
electron energy distribution function (EEDF) from a set of collision cross sections
for some mean discharge conditions. The interface can be used as a preprocessing

| 17
[Link] Page 18 Saturday, October 21, 2023 10:31 AM

stage before solving a full space dependent model. The main purpose of this
interface is to compute electron source coefficients and transport properties.

D RIFT D IFFUSION I NTERFACE

The Drift Diffusion interface is used to compute the electron density and mean
electron energy for any type of plasma. A wide range of boundary conditions are
available to handle secondary emission, thermionic emission, and wall losses. This
interface rarely needs to be used by itself as it makes up part of the application
specific interfaces described later.

CHARGE TRANSPORT
The Charge Transport interface computes the density of charge carriers in a
background gas under the assumption that the transport is dominated by
migration. This is typically only used as part of the Corona Discharge multiphysics
interface (see below).

T HE H EAVY S PECIES I NTERFACE

The Heavy Species Transport interface solves a mass balance equation for all
nonelectron species. This includes charged, neutral, and electronically excited
species. The interface also allows you to add electron impact reactions, chemical
reactions, surface reactions, volumetric species, and surface species via the Model
Builder. This interface rarely needs to be used by itself as it makes up part of the
application specific interfaces described later.

ELECTRICAL BREAKDOWN DETECTION

The Electrical Breakdown Detection interfaces uses an approximate method to
determine if electrical breakdown will occur in a given design by integrating
Townsend growth coefficients along electric field lines.

Application Specific Interfaces

The Plasma Module provides four application specific multiphysics interfaces for
modeling the most common types of discharge.

P LASMA
The Plasma interface can be used to model positive columns, DBD discharges,
glow discharges, and corona discharges. The complicated coupling between the
electron transport, heavy species transport, and electrostatic field is handled
automatically by the software. Furthermore, the secondary emission flux from ion

18 |
[Link] Page 19 Saturday, October 21, 2023 10:31 AM

bombardment on an electrode is automatically computed and used in the

boundary condition for electrons.

P LASMA , T IME P ERIODIC

The Plasma, Time Periodic interface can be used to model capacitively coupled
plasmas. Instead of solving the problem in the time domain, the periodic
steady-state solution is computed. This avoids having to solve for tens or hundreds
of thousands of RF cycles, which is typically how long it takes before the plasma
reaches the periodic steady-state solution. This approach maintains all the
nonlinearity of the model while dramatically reducing computation time. The
physics interface accomplishes this by attaching an extra dimension to the
underlying mathematical equations representing one RF cycle, and enforcing
periodic boundary conditions in the aforementioned extra dimension.

I NDUCTIVELY C OUPLED P LASMA

The Inductively Coupled Plasma interface can be used to model discharges
sustained through induction currents. These discharges typically operate in the
MHz frequency range. Inductively coupled plasmas (ICP) are important in plasma
processing and plasma sources because the plasma density can be considerably
higher than in capacitively coupled discharges. Inductively coupled plasmas are
also attractive from the modeling perspective because they are relatively
straightforward to model, due to the fact that the induction currents can be solved
for in the frequency domain. This means that the RF cycle applied to the driving
coil does not need to be explicitly resolved when solving. As such, the quasi
steady-state solution is reached in relatively few time steps.

I NDUCTIVELY C OUPLED P LASMA WITH RF BIAS

The Inductively Coupled Plasma with Rf Bias interface can be used to model
discharges sustained through induction currents and that have a periodic RF
biased electrode. This interface uses the Plasma, Time Periodic interface to solve
the periodic steady state for the RF bias. The inductive currents are solved in the
frequency domain as in the Inductively Coupled Plasma interface. This type of
discharges are interesting because the plasma density and the ions flux at a surface
can be controlled somewhat independently.

M ICROWAVE P LASMA
The Microwave Plasma interface can be used to model discharges which are
sustained through heating of electrons due to electromagnetic waves. These
discharges typically operate in the GHz frequency range. Wave–heated discharges
usually fall into one of two categories: discharges with no external DC magnetic

| 19
[Link] Page 20 Saturday, October 21, 2023 10:31 AM

field and discharges with a high intensity static magnetic field. If a suitably high
DC magnetic field is present then electron cyclotron resonance (ECR) can occur
where electrons continually gain energy from the electric field over 1 RF period.
Modeling microwave plasmas involves solving equations for the electron density,
mean electron energy, heavy species, the electrostatic potential, and the high
frequency electric field. The high frequency electric field is computed in the
frequency domain and losses are introduced via a complex plasma conductivity.

C ORONA D ISCHARGE
The Corona Discharge interface uses a simplified charge transport model coupled
with electrostatics to provide an approximate method of computing the charge
density and the electrostatic field in corona discharges.

E QUILIBRIUM D ISCHARGES , O UT - OF -P LANE C URRENTS

The Equilibrium Discharges, Out-of-plane Currents multiphysics interface,
available in 2D and 2D axisymmetric, is used to study equilibrium discharges in a
magnetohydrodynamics (MHD) framework where the currents are out-of-plane.
This multiphysics interface adds three single physics interfaces: Magnetic Fields,
Heat Transfer in Fluids, and Laminar Flow, together with several multiphysics
coupling features. The multiphysics couplings add the MHD coupling between
the Magnetic Fields and the Laminar Flow interfaces. The multiphysics couplings
also add heating and cooling of the equilibrium plasma by enthalpy transport,
Joule heating and radiation loss.

E QUILIBRIUM D ISCHARGES , I N -P LANE C URRENTS

The Equilibrium Discharges, in-plane Currents multiphysics interface,
available in 2D and 2D axisymmetric, is used to study equilibrium discharges in a
magnetohydrodynamics (MHD) framework where the currents are in-plane. This
multiphysics interface adds three single physics interfaces: Magnetic and Electric
Fields, Heat Transfer in Fluids, and Laminar Flow, together with several
multiphysics coupling features. The multiphysics couplings add the MHD
coupling between the Magnetic and Electric Fields and the Laminar Flow
interfaces. The multiphysics couplings also add heating and cooling of the
equilibrium plasma by enthalpy transport, Joule heating and radiation loss as well
as special boundary conditions to model the ion and electron heating at the plasma
boundaries.

E QUILIBRIUM D ISCHARGES
The Equilibrium Discharges multiphysics interface, available in 3D, is used to
study equilibrium discharges in a magnetohydrodynamics (MHD) framework.

20 |
[Link] Page 21 Saturday, October 21, 2023 10:31 AM

This multiphysics interface adds three single physics interfaces: Magnetic and
Electric Fields, Heat Transfer in Fluids, and Laminar Flow, together with several
multiphysics coupling features. The multiphysics couplings add the MHD
coupling between the Magnetic and Electric Fields and the Laminar Flow
interfaces. The multiphysics couplings also add heating and cooling of the
equilibrium plasma by enthalpy transport, Joule heating, and radiation loss as well
as special boundary conditions to model the ion and electron heating at the plasma
boundaries.

L IMITATIONS OF THE P LASMA M ODULE

The Plasma Module cannot currently model magnetrons and space thrusters due
to the fact that the perpendicular versus cross field electron mobility and diffusivity
can be on the order of 108. This poses severe numerical challenges and inevitably
there is unphysical numerical diffusion across the magnetic field lines. As such, a
truly self-consistent model of a Hall thruster or magnetron is not currently
possible, although it is still possible to gain valuable insight into the physics. The
Plasma module also cannot model plasmas that are not collisional enough for the
fluid-type equations used to be valid. For reactors with characteristic dimensions
of 10 cm the lower possible pressure would be 20 mTorr.

References

1. J.P. Boeuf and others, “Physics and Modelling of a Negative Ion Source Prototype for
the ITER Neutral Beam Injection,” [Link]
PDFplus/2010/cn180/cn180_papers/itr_p1-[Link].

2. G.J.M. Hagelaar, and L.C. Pitchford, “Solving the Boltzmann Equation to Obtain
Electron Transport Coefficients and Rate Coefficients for Fluid Models,” Plasma Sources
Sci. Technol , vol. 14, pp. 722–733, 2005.
3. D.P. Lymberopoulos 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, pp. 473–494, 1995.

4. A. Bogaerts, E. Neyts, R. Gijbels, and J. van der Mullen, “Gas discharge plasmas and
their applications,” Spectrochimica Acta Part B , vol. 57, pp. 609–658, 2002.

5. M.A. Lieberman and A.J. Lichtenberg, “Principles of Plasma Discharges and Materials
Processing,” Wiley, 2005.

| 21
[Link] Page 22 Saturday, October 21, 2023 10:31 AM

Argon Discharge in the GEC Reference Cell

The GEC cell was introduced by NIST in order to provide a standardized platform
for experimental and modeling studies of discharges in different laboratories
(Ref. 1). The plasma is sustained via inductive heating. The Reference Cell
operates as an inductively-coupled plasma in this mode.
Inductively coupled discharges typically operate at low pressures (<10 Pa) and
high charge density (>1017 m-3). High density plasma sources are popular because
low pressure ion bombardment can provide a greater degree of anisotropy on the
surface of the wafer.

5 turn copper coil, 13.56 MHz (blue)

Dielectric material (gray)

Plasma domain (pink)

Wafer and pedestal

GEC ICP reactor geometry consisting of a 5-turn copper coil, plasma volume, dielectrics, and wafer with
pedestal.

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

22 |
[Link] Page 23 Saturday, October 21, 2023 10:31 AM

electronically excited states can be lumped into a single species, which results in a
chemical mechanism consisting of only 3 species and 7 reactions:
LIST OF VOLUMETRIC PLASMA CHEMICAL REACTIONS CONSIDERED
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 -

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, and 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:
LIST OF SURFACE REACTIONS CONSIDERED
REACTION FORMULA STICKING COEFFICIENT

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).
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, the coil is driven with a fixed power, 1500 W.

| 23
[Link] Page 24 Saturday, October 21, 2023 10:31 AM

Model Wizard

The first step to build a model is to open the COMSOL Desktop and then specify
the type of analysis you want to do — in this case, a frequency-transient inductively
coupled plasma analysis. The frequency-transient study type means that the high
frequency electromagnetic field is computed in the frequency domain and all other
variables are computed in the time domain.
Note: These instructions are for the user interface on Windows® but apply, with
minor differences, also to Linux® and macOS.
1 Open COMSOL Multiphysics. In the New window, click the Model Wizard
button .
2 In the Model Wizard window, click the 2D Axisymmetric button .
3 In the Select physics tree, select Plasma>Inductively Coupled Plasma .
4 Click the Add button.
5 Click the Study button .
6 In the tree, select Frequency-Transient .
7 Click Done .

I m p o r t i n g t h e G eo m e t r y

The next step is to create your geometry, which also can be imported from an
external program. COMSOL Multiphysics supports a multitude of CAD
programs and file formats. In this example, insert the prepared geometry sequence
from a file in the COMSOL Multiphysics file format (.mph). The file contains the
geometry of the GEC reference cell.
The location of the files used in this exercise varies based on your installation. For
example, if the installation is on your hard drive, the file path might be similar to
C:\Program Files\COMSOL\COMSOL62\applications\.
1 In the Geometry toolbar, click Insert Sequence .
2 Browse to the folder Plasma_Module\Inductively_Coupled_Plasmas under the
COMSOL installation directory and double-click on the file
argon_gec_icp_geom.mph.

24 |
[Link] Page 25 Saturday, October 21, 2023 10:31 AM

3 Click the Import button. The geometry should appear in the Graphics window
as shown below.

Definitions

1 In the Definitions toolbar, click Explicit .

2 In the Settings window, type Walls in the Label field.
3 Locate the Input Entities section.
4 From the Geometric entity level list, choose Boundary.
5 Go to the Home tab on the Model desktop toolbar, then choose
Windows>Selection List.
6 Select Boundaries 6, 8, 35–38, 44, 45, and 51–56 only (by holding shown shift
and clicking on the list), then click the Add to selection button at the top
of the Selection List settings window. The selected boundaries in the graphics
window will turn blue, indicating that the selection is confirmed.
7 Click back on the Model Builder tab.
8 In the Definitions toolbar, click Explicit .
9 In the Settings window, type Coils in the Label field.

| 25
[Link] Page 26 Saturday, October 21, 2023 10:31 AM

10Click on the Selection List tab and select Domains 6 and 8–11 only, then click
the Add to selection button at the top of the Selection List settings window.
11Click back on the Model Builder tab.
12In the Definitions toolbar, click Explicit .
13In the Settings window, type Coil Boundaries in the Label field.
14Click on the Selection List tab and select Domains 6 and 8–11 only, then click
the Add to selection button at the top of the Selection List settings window.
15Click back on the Model Builder tab. In the Explicit settings window, locate the
Output Entities section.
16From the Output entities list, choose Adjacent boundaries.
The settings for the selection should be as given below:

17Go to the Model desktop toolbar.

18In the Home toolbar, click Parameters and select Parameters 1 .
19In the Parameters settings window, locate the Parameters section.

26 |
[Link] Page 27 Saturday, October 21, 2023 10:31 AM

In the table, enter the following settings:

Plasma and Magnetic Fields

1 In the Model Builder window, under Component 1 click on the Plasma

node.
2 In the settings window, Locate the Plasma Properties section and select the Use
reduced electron transport properties check box.
3 Click on the Selection List window. Select Domain 3 only.
4 In the Model Builder window, right-click Plasma and choose Cross Section
Import.
5 In the Cross Section Import settings window, locate the Cross Section Import
section.
6 Click the Browse button.
7 Browse to the module Application Library folder and double-click the file
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.
1 In the Model Builder window, right-click Plasma and choose the domain setting
Heavy Species Transport>Reaction .
2 In the Reaction settings window, locate the Reaction Formula section.
3 In the Formula field, type Ars+Ars=>e+Ar+Ar+. Click off the settings window.
4 Locate the Kinetics Expressions section. In the kf field, type 3.734E8.
5 In the Model Builder window, right-click Plasma and choose the domain setting
Heavy Species Transport>Reaction .
6 In the Reaction settings window, locate the Reaction Formula section.

| 27
[Link] Page 28 Saturday, October 21, 2023 10:31 AM

7 In the Formula field, type Ars+Ar=>Ar+Ar. Click off the settings window.
8 Locate the Kinetics Expressions section. In the kf field, type 1807.
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.
1 In the Model Builder window, under Component 1>Plasma click Species: Ar
.
2 In the Species settings window, locate the Species Formula section.
3 Select the From mass constraint check box.
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 time step, 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.
1 In the Model Builder window, under Component 1>Plasma click Species: Ar+
.
2 In the Species settings window, locate the Species Formula section.
3 Select the Initial value from electroneutrality constraint check box.
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
time steps.
1 In the Model Builder window, under Component 1>Plasma click Initial Values
1 .
2 In the Initial Values settings window, locate the Initial Values section.
3 In the ne,0 field, type 1E15[1/m^3].
4 In the 0 field, type 5[V].
5 In the Model Builder window, under Component 1>Plasma click Plasma Model
1 .
6 In the Plasma Model settings window, locate the Model Inputs section.
7 In the T field, type T0.
8 In the pA field, type p0.
9 Locate the Electron Density and Energy section. In the eNn field, type mueN.

28 |
[Link] Page 29 Saturday, October 21, 2023 10:31 AM

Surface reactions must always be included in a plasma model since they describe
how ionic, excited, and radical species interact with the wall.
1 In the Model Builder window, under Component 1 right-click Plasma and
choose the boundary condition Heavy Species Transport>Surface Reaction .
2 In the Surface Reaction settings window, locate the Reaction Formula section.
3 In the Formula field, type Ars=>Ar.
4 Locate the Boundary Selection section. From the Selection list, choose Walls.
5 In the Model Builder window, right-click Plasma and choose the boundary
condition Heavy Species Transport>Surface Reaction .
6 In the Surface Reaction settings window, locate the Reaction Formula section.
7 In the Formula field, type Ar+=>Ar.
8 Locate the Boundary Selection section. From the Selection list, choose Walls.
Now add boundary conditions to describe how the electrons interact with the wall
and specify that the walls are grounded.
1 In the Model Builder window, right-click Plasma and choose the boundary
condition Drift Diffusion>Wall .
2 In the Wall settings window, locate the General Wall Settings section.
3 In the re field, type 0.2.
4 Locate the Boundary Selection section. From the Selection list, choose Walls.
5 In the Physics toolbar, click Boundaries and choose Ground .
6 In the Ground settings window, locate the Boundary Selection section.
7 From the Selection list, choose Walls.
You need to compute the AC electric field both inside and outside the plasma. It
is not necessary to compute the high frequency fields in the wafer or wafer
pedestal, so start by modifying the selection for the Magnetic Fields interface.
1 In the Model Builder window, click Magnetic Fields.
2 Select Domains 3–6 and 8–12 only.
The Coil feature allows you to drive the system with a fixed total power. Some of
this power will be dissipated in the coil, the rest will be coupled into the plasma.
In this example 1500 W are applied to the system. This results in a plasma with a
high number density. You also need to specify the gas temperature and pressure.
1 In the Model Builder window, right-click Magnetic Fields and choose the
domain setting Coil .
2 In the Coil settings window, locate the Domain Selection section.
3 From the Selection list, choose Coils.

| 29
[Link] Page 30 Saturday, October 21, 2023 10:31 AM

4 Locate the Coil section. Select the Coil group check box.
In the Coil excitation, select Power and in the Pcoil field, type Psp.

Mate rials

Now assign material properties to the air, coil, and dielectric domains.
1 In the Model Builder window, under Component 1 right-click Materials and
choose Blank Material .
2 In the Material settings window, locate the Geometric Entity Selection section.
3 From the Selection list, choose Coils.
4 Locate the Material Contents section. In the table, enter the following settings:

5 In the Model Builder window, under Component 1 right-click Materials and

choose Blank Material .

30 |
[Link] Page 31 Saturday, October 21, 2023 10:31 AM

6 In the Material settings window, locate the Geometric Entity Selection section.
Select Domain 5 only. In the table, enter the following settings:

7 In the Model Builder window, under Component 1 right-click Materials and

choose Blank Material .

| 31
[Link] Page 32 Saturday, October 21, 2023 10:31 AM

8 In the Material settings window, locate the Geometric Entity Selection section.
Select Domains 4 and 12 only. In the table, enter the following settings:

Me sh

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 effect needs to be resolved.
1 In the Model Builder window, click Mesh 1 .
2 In the Mesh settings window, locate the Mesh Settings section.
3 From the Element size list, choose Finer.
4 Right-click Mesh 1 and choose More Operations>Edge .
5 Select Boundaries 6, 8, 44, 45, and 54 only.
6 Right-click Component 1>Mesh 1>Edge 1 and choose Size .
7 In the Size settings window, locate the Element Size section.
8 Click the Custom button.

32 |
[Link] Page 33 Saturday, October 21, 2023 10:31 AM

9 Locate the Element Size Parameters section. Select the Maximum element size
check box.
10In the associated text field, type 1E-3.
11In the Model Builder window, right-click Mesh 1 and choose Free
Triangular .
12In the Free Triangular settings window, locate the Domain Selection section.
13From the Geometric entity level list, choose Domain.
14Select Domain 3 only.
15Right-click Component 1>Mesh 1>Free Triangular 1 and choose Size .
16In the Size settings window, locate the Element Size section.
17From the Predefined list, choose Extra fine.
18In the Model Builder window, right-click Mesh 1 and choose Boundary
Layers .
19In the Boundary Layers settings window, locate the Domain Selection section.
20From the Geometric entity level list, choose Domain.
21Select Domain 3 only.
22In the Model Builder window, under Component 1>Mesh 1>Boundary Layers
1 click Boundary Layer Properties .
23In the Boundary Layer Properties settings window, locate the Boundary
Selection section.
24From the Selection list, choose Walls.
25Locate the Boundary Layer Properties section. In the Number of boundary
layers field, type 5.
26In the Boundary layer stretching factor field, type 1.4.
27In the Model Builder window, right-click Mesh 1 and choose Mapped .
28In the Mapped settings window, locate the Domain Selection section.
29From the Geometric entity level list, choose Domain.
30From the Selection list, choose Coils.
31Right-click Component 1>Mesh 1>Mapped 1 and choose Distribution .
32In the Distribution settings window, locate the Boundary Selection section.
33From the Selection list, choose Coil Boundaries.
34In the Distribution settings window, locate the Distribution section.
35From the Distribution type list, choose Predefined distribution type.
36In the Number of elements field, type 25.
37In the Element ratio field, type 20.

| 33
[Link] Page 34 Saturday, October 21, 2023 10:31 AM

38From the Growth formula list, choose Geometric sequence.

39Select the Symmetric distribution check box.
40In the Model Builder window, right-click Mesh 1 and choose Free
Triangular .
41Right-click Mesh 1 and choose Build All .

Study

The model is now ready to be solved. You just need to enter the list of output
times and the frequency for the induction currents.
1 In the Model Builder window, expand the Study 1 node, then click Step 1:
Frequency-Transient .
2 In the Frequency-Transient settings window, locate the Study Settings section.
3 In the Times field, type 0 (the number zero followed by a space)
4 Click the Range button .
5 Go to the Range dialog box.
6 From the Entry method list, choose Number of values.
7 In the Start field, type -8.
8 In the Stop field, type -3.
9 In the Number of values field, type 21.
10From the Function to apply to all values list, choose exp10. The text field should
now contain 0 10^{range(-8,5/20,-3)}.
11Click the Add button.
12In the Frequency field, type 13.56E6.
13In the Model Builder window, right-click Study 1 and choose Compute .

Results

There are three default plots generated, one for the electron number density, one
for the electron temperature and one for the electrostatic potential. After the
model has solved, the default plot is of the electron density. The peak electron
density occurs at the center of the reactor, underneath the RF coil. The electron

34 |
[Link] Page 35 Saturday, October 21, 2023 10:31 AM

density in this case is high enough to cause some shielding of the azimuthal electric
field.

1 To view the electron temperature, click on the Electron Temperature (plas) plot
group .
The electron “temperature” is highest directly underneath the coil, which is where
the bulk of the power deposition occurs.

| 35
[Link] Page 36 Saturday, October 21, 2023 10:31 AM

2 To view the electric potential, click on the Electric Potential (plas) plot
group .

36 |
[Link] Page 37 Saturday, October 21, 2023 10:31 AM

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 current.
3 In the Model Builder window, right-click Results and choose 1D Plot
Group .
4 Right-click Results>1D Plot Group 6 and choose Global .
5 In the Global settings window, locate the y-Axis Data section. Click Replace
Expression . Select Magnetic Fields>Coil Parameters>Coil Resistance
(mf.RCoil_1) from the list.
6 Click the x-Axis Log Scale button in the Graphics toolbar.
7 Click the Plot button .

Now add a plot for the total power dissipated in the system. This is one half of the
real part of the total voltage drop over the coil multiplied by the applied current.
8 In the Model Builder window, right-click Results and choose 1D Plot
Group .
9 Right-click Results>1D Plot Group 5 and choose Global.
10In the Global settings window, locate the y-Axis Data section. Click Replace
Expression . Select Magnetic Fields>Coil Parameters>Coil Power
(mf.PCoil_1) from the list.
11Click the x-Axis Log Scale button in the Graphics toolbar.

| 37
[Link] Page 38 Saturday, October 21, 2023 10:31 AM

12Click the Plot button .

Now you add some additional two dimensional plots.

1 In the Model Builder window, right-click Results and choose 2D Plot
Group .
2 Right-click 2D Plot Group 6 and choose Surface .
3 In the Surface settings window, click Replace Expression in the
upper-right corner of the Expression section. From the menu, choose Plasma
(Heavy Species Transport)>Number densities>Number density
(plas.n_wAr_1p).
4 Click the Plot button .
5 Click the Zoom Extents button in the Graphics toolbar.

38 |
[Link] Page 39 Saturday, October 21, 2023 10:31 AM

A quick way of creating additional plots is to use the Duplicate option. Now create
a plot of the norm of the high frequency electric field.
6 In the Model Builder window, right-click 2D Plot Group 6 and choose
Duplicate .
7 In the Model Builder window, expand the 2D Plot Group 7 node, then click
Surface 1 .
8 In the Surface settings window, click Replace Expression in the
upper-right corner of the Expression section. From the menu, choose Magnetic
Fields>Electric>Electric field norm ([Link]).
9 Click the Plot button .
10Click the Zoom Extents button in the Graphics toolbar.

| 39
[Link] Page 40 Saturday, October 21, 2023 10:31 AM

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. Now create a plot of the number
density of electronically excited argon atoms.
11In the Model Builder window, right-click 2D Plot Group 7 and choose
Duplicate .
12In the Model Builder window, expand the 2D Plot Group 8 node, then click
Surface 1 .
13In the Surface settings window, click Replace Expression in the
upper-right corner of the Expression section. From the menu, choose Plasma
(Heavy Species Transport)>Number densities>Number density (plas.n_wArs).
14Click the Plot button .
15Click the Zoom Extents button in the Graphics toolbar.

40 |
[Link] Page 41 Saturday, October 21, 2023 10:31 AM

Finally, create a new dataset only active on the plasma domain in order to make it
easier to visualize the power deposition into the plasma.
1 In the Model Builder window, expand the Results>Datasets node.
2 Right-click Solution 1 and choose Duplicate .
3 Right-click Results>Datasets>Study 1/Solution 1 (2) and choose Selection.
4 In the Selection settings window, locate the Geometric Entity Selection section.
5 From the Geometric entity level list, choose Domain.
6 Select Domain 3 only.
7 Deactivate Propagate to lower dimensions.
8 In the Model Builder window, right-click Results and choose 2D Plot
Group .
9 Right-click 2D Plot Group 9 and choose Surface .
10In the Surface settings window, locate the Data section.
11From the Dataset list, choose Study 1/Solution 1(2).
12Click Replace Expression in the upper-right corner of the Expression
section. From the menu, choose Magnetic Fields>Heating and
losses>Volumetric loss density ([Link]).
13Click the Plot button .

| 41
[Link] Page 42 Saturday, October 21, 2023 10:31 AM

14Click the Zoom Extents button in the Graphics toolbar.

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.

Bibliography

1. D.P. Lymberopoulos 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,
pp. 473–494, 1995.

42 |