 

 

LIBRARY
MIchIgan State
UnIversIIy

 

 

 

PLACE II RETURN BOX to remove this mutton. your mood.
' To AVOID FINES mum on or baton duo duo.

DATE DUE DATE DH; DATE DUE

 

 

 

 

JUN . « m?
E3 7:; '5‘ ’3“,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU lsAnNﬂnnIIIvo Action/Emil OppoﬂunIIy Inﬂation
WWI

 

MODELING THE ELECTROMAGNETIC FIELD AND THE PLASMA
EXCITATION IN A MODERATE PRESSURE MICROWAVE CAVITY PLAMSA
SOURCE

By

Wen-yi Tan

A DISSERTATION

Submitted to
Michigan State University
in partial fulﬁllment of the requirements

for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering

1 994

ABSTRACT

MODELING THE ELECTROMAGNETIC FIELD AND THE PLASMA
EXCITATION IN A MODERATE PRESSURE MICROWAVE CAVITY
PLAMSA SOURCE

By

Wen-yi Tan

Microwave plasma sources are used for the generation of plasma discharges in a
number of applications including diamond thin ﬁlm deposition. A numerical model has
been developed for these microwave plasma sources which includes a electromagnetic
ﬁeld model and a ﬂuid plasma model. The microwave plasma source modeled was a
cylindrical, single mode excited cavity with an input power probe for coupling the
microwave energy into the cavity. The time-varying electromagnetic ﬁelds inside the
resonant cavity, b0th inside and outside the discharge region, are obtained using the
electromagnetic ﬁeld model which incorporates a ﬁnite-difference time-domain (FDTD)
method to solve Maxwell's equations. The microwave electric ﬁeld interactions with the
plasma discharge are described using a ﬁnite difference solution of the electron
momentum transport equation. The characteristics of the discharge are simulated by the
ﬂuid plasma model which solves the electron and ion continuity equations, electron
energy balance equation, and the Poisson equation. A ﬁnal self-consistent solution is

obtained by iteratively solving the electromagnetic ﬁeld model and the ﬂuid plasma

model.

In particular, a TMOIn mode, 17.78 cm i.d. microwave cavity plasma reactor used
for diamond thin ﬁlm deposition is simulated using this numerical model. The spatial
eleCtric ﬁeld patterns, power absorption patterns, and quality factor of the cavity loaded
with a hydrogen discharge are investigated in the moderate pressure range (1 Torr - 100
Torr). The physical behavior of the hydrogen discharge, such as plasma density, electron
temperature, and plasma potential, are also simulated and analyzed for various input
conditions. The simulated results are compared with experimental data. Calculations of
the electromagnetic ﬁelds using the FDTD techniques are also done for a compact electron
cyclotron resonance (ECR) plasma source and for an argon coaxially loaded microwave
cavity plasma source. Additionally, microwave cavities loaded with lossy loads of various
conductivities were simulated and compared to exact analytical solution. Agreement was

found between the numerical FDTD solution and the analytical solution.

To my parents

iv

ACKNOWLEDEGMENTS

The author wishes to thank Dr. Timothy A. Grotjohn for his encouragement,
guidance and support thoughtout this research. Thanks are also given to guidance
committee members, Dr. Jes Asmussen, Dr. Donnie Reinhard, and Dr. David Yen for
taking time to serve on the examining committee.

This research was supported in part by grants from National Science Foundation.

TABLE OF CONTENTS

LIST OF TABLES ............................................................................................................. x
LIST OF FIGURES .......................................................................................................... xi
Chapter 1 Introduction ......................................................................................... 1
1.1 Introduction ......................................................................................................... 1

1.2 Objective ............................................................................................................ 4

1.3 Outline ................................................................................................................ 5
Chapter 2 Electromagnetic Field Model .............................................................. 8
2.1 Introduction ........................................................................................................ 8
2.2 Review of electromagnetic ﬁeld model ........................................................... 10
2.2.1 Frequency domain ................................................................................... 11

2.2.2 Time domain ........................................................................................... 13

2.3 FDTD model implementation in cylindrical coordinates ................................ 20
2.4 Conductivity models for FDTD model implementation .................................. 25
2.4.1 Conductivity model for non-dispersive lossy materials .......................... 26

2.4.2 Plasma conductivity model ..................................................................... 27

vi

2.4.2.1 Time harmonic domain solution .............................................. 28

2.4.2.2 Time domain solution .............................................................. 32
2.5 Summary .......................................................................................................... 33
Chapter 3 Numerical Models of Plasma Discharges ......................................... 34
3.1 Introduction ...................................................................................................... 34
3.2 Review of plasma discharge models ................................................................ 34
3.2.1 Particle model ......................................................................................... 35
3.2.2 Continuum (ﬂuid) model ........................................................................ 37
3.3 Model development ......................................................................................... 40
3.3.1 Distribution function ............................................................................... 40
3.3.2 Boltzmann transport equation ................................................................. 43
3.3.3 Moment equations ................................................................................... 44
3.3.4 Drift-Diffusion model ............................................................................ 47
3.3.5 Fluid model ............................................................................................ 49
3.4 Numerical solution implementation ................................................................. 52
3.4. l Discretization .......................................................................................... 53
3.4.2 Numerical solution .................................................................................. 58
3.5 Summary .......................................................................................................... 62
Chapter 4 Transient State Simulation of Electromagnetic Fields Inside
Empty and Loaded Microwave Resonant Cavities ........................... 63
4.1 Introduction ...................................................................................................... 63
4.2 Simulation of electromagnetic ﬁelds inside empty cavities ............................ 64

vii

4.3 Natural frequency response for loaded lossy cavity ........................................ 71
4.4 Simulation results of an argon plasma loaded cavity ....................................... 85
4.5 Implementation and results for a compact ECR ion source ............................. 89
4.6 Summary .......................................................................................................... 94
Chapter 5 Steady State Simulation of Electromagnetic Fields Inside

Plasma Loaded Microwave Resonant Cavities ............................... 100
5. 1 Introduction .................................................................................................... 100
5.2 Non-reﬂecting boundary conditions .............................................................. 103
5.3 Implementation for a coaxial waveguide ....................................................... 107
5.4 Implementation for a Microwave Cavity Plasma Reactor ............................. 112
5.5 Simulation results for a cavity reactor loaded with a H2 discharge .............. 124
5.6 Conclusions .................................................................................................... 135

Chapter 6 Self-consistent Simulation of a Microwave

Cavity Plasma Reactor .................................................................... 138
6. 1 Introduction .................................................................................................... 138
6.2 Model description .......................................................................................... 138
6.3 Reactor description and Q value measurement .............................................. 142
6.4 Implementation for H2 discharge simulation ................................................ 150
6.5 Simulation results for ﬂuid plasma model ..................................................... 156
6.6 Self-Consistent results for microwave cavity plasma reactor simulation ...... 159
6.7 Summary ........................................................................................................ 174

viii

LIST OF TABLES

Table 1: H2 reaction rates. ............................................................................................ 155

Figure 1-1:
Figure 2-1:
Figure 2-2:
Figure 2-3:

Figure 3-1:

Figure 3-2:

Figure 4-1:

Figure 4-2:

Figure 4-3:

Figure 4-4:

Figure 4-5:

Figure 4-6:

LIST OF FIGURES

Overview diagram. ....................................................................................... 6
Schematic description of microwave discharge generation. ........................ 9
A FDTD unit cell in three-dimensional cylindrical coordinates. ............... 16
Leap-frog scheme. ...................................................................................... 18
The distribution function f can be described by

three moments: n, <v>, and <§>. ............................................................... 42
Staggered—mesh arrangement. .................................................................... 54

TE011’ TE111, Wow and TM012 mOdc
excitation techniques in a cylindrical cavity. ............................................. 65

Electric ﬁeld energy, {7 j 80E2d V, versus time for the TED” mode. The
V

empty cylindrical cavity was 15.25 cm in diameter and 6.5 cm high. ....... 68
Radial dependence of IE¢I at four different 2 locations

for the TECH mode. ................................................................................... 69
Grid structure for FDTD simulation in cylindrical coordinate. ................. 70

2
Electric ﬁeld energy, ‘1]! 80E d V, versus time for the TMOH mode.

V
The cylindrical empty cavity was 17.8 cm in diameter
and 7.2 cm high. The resonant frequency was 2.45 GHz. .......................... 72
Total stored energy, U, versus time for the TMOU mode. ......................... 73

xi

Figure 4-7:

Figure 4-8:

Figure 4-9:

Figure 4-10:

Figure 4-11:

Figure 4—12:

Figure 4-13:

Figure 4-14

Figure 4- 15:

Figure 4—16:

Figure 4-17:

Figure 4-18:
Figure 4-19:

Figure 4-20:

Figure 4-21:

IE,I (arbitrary units) versus 2 at four different r locations
for the TMOH mode. .................................................................................. 74

IE,| (arbitrary units) versus 2 at four different r locations for the

TM012 mode. The cylindrical empty cavity was 17.8 cm in diameter

and 14.41 cm high. The resonant frequency was 2.45 GHz. ...................... 75
|E,l (arbitrary units) versus 4) at three different r locations

for the TB] 11 mode. The cylindrical empty cavity was 8.9 cm

in diameter and 6.9 cm high. The resonant frequency was 2.45 GHz. ....... 76
Conﬁguration of a coaxially loaded cylindrical cavity. ............................. 79
Fine grid structure for TED 11 mode lossy cavity simulation. ..................... 80
Space averaged IEI2 (arbitrary units), versus time, for the TED“ mode.
Results are shown for three different load conductivities after

the source is turned off ............................................................................... 81

Total stored energy, U, versus time for the TED“ mode
for three different load conductivities. ....................................................... 82

Radial dependence of IE¢| (arbitrary units) for the TECH mode for three
different conductivity loads in a 15.25 cm i.d. cylindrical cavity. ............. 83

S-plane plot for the complex frequency as a function of
load conductivity for the TEO 11 mode. ...................................................... 84

The geometry of a 17.8 cm i.d. cylindrical cavity with a 12 mm i.d.
quartz tube loaded with an argon plasma .................................................. 86

Electron density, ne, and effective collision frequency, veﬁ;

versus pressure for the argon plasma shown in Figure 4-16. ..................... 88
Loaded cavity quality factor, Q, versus pressure. ...................................... 90
Absorbed power density versus pressure. ................................................... 91
Simulation structure for an compact ECR ion source. ............................... 92

. l 2 . .
Electric ﬁeld energy, ‘7] EOE d V, versus time for the structure shown in
V

Figure 4-20 with no plasma load present.The excitation of the
source was turned off at time = 0.6 nsec ..................................................... 93

xii

Figure 4—22:
Figure 4-23:

Figure 4—24:

Figure 4-25:

Figure 5- 1:

Figure 5-2:
Figure 5-3:

Figure 5-4:

Figure 5-5:

Figure 5-6:

Figure 5-7:
Figure 5-8:

Figure 5-9:

Figure 5—10:

Figure 5-1 1:

Figure 5-12:

Figure 5-13:

lErl (arbitrary units) versus r and z. ............................................................ 95
IEZI (arbitrary units) versus r and z. ............................................................ 96

IE,| (arbitrary units) versus r and z for the compact ion source
operated with a helium plasma discharge. ................................................. 97

IEZI (arbitrary units) versus r and z for the compact ion source
operated with a helium plasma discharge. ................................................. 98

2
Steady-state response of the electric ﬁeld energy, JV] 80E (1 V,
V

versus time for the TED 11 mode. .............................................................. 101
Grid termination structure. ....................................................................... 106
Simulation structure for the cylindrical coaxial waveguide. ................... 108

Non-reﬂecting boundary simulation results:

(above) spatial average E, vs. time at two different cross section

(solid: 2:30, dash: 2:0). (below) spatial average power ﬂow vs. time

(solid: z=10, dash: 2:20). ......................................................................... 110

Spatial variation (r vs. 2) of the IE,I component for the
TEM mode in a cylindrical coaxial waveguide. ...................................... 111

Conﬁguration of a microwave cavity plasma reactor used for
diamond ﬁlm deposition. The cavity diameter is 17.78 cm. .................... 113

Simpliﬁed simulation structure for a microwave cavity plasma reactor. . 114
Electric ﬁeld energy vs. time for TMO12 mode cavity simulation. .......... 117

IE,I ﬁeld distribution for a TMO 12 mode microwave cavity
plasma reactor. ......................................................................................... 118

lEzl ﬁeld distribution for a TMO 12 mode microwave cavity
plasma reactor. ......................................................................................... 119

Field patterns of TM012 mode cylindrical cavity. .................................... 120

Q vs. cavity height for three different conductivity
loaded TM012 mode cavity simulation. ................................................... 122

Reﬂected power vs. cavity height for TM012 mode cavity simulation. 123

Figure 5-14:
Figure 5-15:
Figure 5-16:

Figure 5-17:

Figure 5-18:

Figure 5-19:

Figure 5-20:

Figure 5-21:

Figure 5-22:

Figure 6—1:
Figure 6-2:

Figure 6-3:

Figure 6-4:

Figure 6—5:

Figure 6-6:

Figure 6-7:

Figure 6-8:

Simulation structure cross-section for microwave cavity plasma reactor. 126
Q value vs. electron density at a pressure of 20 Torr. .............................. 128
Q value vs. electron density at a pressure of 70 Torr. .............................. 129

lE,l distribution in the r - 2 plane at 20 Torr for a plasma

density=3 x1012cm -3. ......................................................................... 131
lEzl distribution in the r - zzplane at 20 Torr for

a plasma density = 3 x 101 cm '3. ........................................................... 132
IErl distribution in the r - 2 plane at 70 Torr

for a plasma density = 8 x 1013 cm '3. ..................................................... 133
lEzl distribution in the r - 2 plane at 70 Torr

for a plasma density = 8 x 1013 cm '3. ..................................................... 134
Power density distribution in the r - 2 plane at 20 Torr

for a plasma density = 3 x 1012 cm '3. ..................................................... 136
Power density distribution in the r - 2 plane at 70 Torr

for a plasma density = 8 x 1013 cm ‘ . .................................................. 137
Flow chart for microwave cavity plasma reactor simulation. .................. 141
Microcoax probe and holes in the cavity walls. ....................................... 143
Electric ﬁeld strength vs. cavity height for microwave

cavity plasma reactor during diamond deposition process. ..................... 148
Field patterns for mom mode. ............................................................... 149
Simulation structure cross-section for microwave cavity

plasma reactor with quartz tube. .............................................................. 151
Grid structure for ﬂuid plasma model. ..................................................... 153

Simulation results of microwave cavity plasma reactor using

ﬂuid plasma model with uniform absorbed power density.

(a) Electron density (b) Plasma potential

(0) Electron temperature (c) Neutral temperature (assigned values). ...... 157

Electron temperature and density variations
during ﬂuid plasma simulation. ............................................................... 158

xiv

Figure 6—9: E, ﬁeld distribution. ................................................................................. 161

Figure 6-10: Ez ﬁeld distribution. ................................................................................. 162
Figure 6-1 1: Conductivity distribution. ........................................................................ 163
Figure 6-12: Absorbed power patterns. ........................................................................ 164
Figure 6-13: B, ﬁeld pattern inside the quartz disk region. . .......................................... 165
Figure 6-14: Ez ﬁeld pattern inside the quartz disk region. .......................................... 166

Figure 6-15: Simulation results of a microwave cavity plasma reactor
for pressure = 50 Torr and input microwave power = 1500 watt.
(a) Electron density (b) Plasma potential
(c) Electron temperature ((1) Neutral temperature (assigned values). ...... 168

Figure 6-16: Spatial variation of ionization rate and recombination rate
for pressure = 50 Torr and input microwave power = 1500 watt. ........... 169

Figure 6-17: Spatial distribution of various energy loss
for pressure = 50 Torr and input microwave power = 1500 watt. ........... 170

Figure 6-18: Absorbed power vs. number of iterations. ............................................... 171

Figure 6-19: Neutral temperature and plasma volume
vs. incident power for various pressures. ................................................. 172

Figure 6-20: Electron temperature vs. absorbed power for three different pressures. .. 173

Figure 6-21: Q value vs. absorbed power for three different pressures. ....................... 175

XV

Chapter 1

Introduction

1.1 Introduction

Microwave excited plasmas have demonstrated excellent potential in high-
pressure (>100 Torr), moderate-pressure (1 Torr - 100 Torr), and low-pressure (0.1 mTorr
- 1 Torr) plasma processing applications because of their ability to create high densities of
excited and charged species without high plasma sheath potentials and contamination
from electrodes. One way to create a microwave excited plasma is utilizing a single mode
microwave resonance cavity which couples microwave energy into the discharge with or
without a static multipolar magnetic ﬁeld [1]. Previous studies have shown that
microwave plasma reactor sources can efﬁciently produce microwave discharges at high
pressures for thermal processing applications [2], at moderate pressure (10 -100 Torr) for
diamond ﬁlm deposition [3], and at low pressure for silicon and III-V etching applications
[4115]-

The further development of microwave plasma processing technology requires
that the physical behavior of discharges inside microwave plasma reactors be better

understood and characterized. One of the major issues regarding this ﬁeld is the plasma

2

heating mechanism, and in particular, the electromagnetic excitation mechanism of
discharges inside microwave plasma sources. Generally, the discharges inside these
sources are generated and sustained by the microwave ﬁelds either through collisional
(Joule) heating or through electron cyclotron resonance (ECR) heating if an appropriate
static magnet ﬁeld is present. Inside the discharges, since the mass of electrons is much
less than the mass of ions, electrons are usually accelerated to higher energies than the
ions by microwave ﬁelds. Therefore, the electromagnetic ﬁeld energy is mainly imparted
to the electron gas, and the discharge heating mechanism is mainly determined by the

behavior of electrons.

Under moderate pressure conditions, the ECR heating mechanism becomes
insigniﬁcant since the mean free path of electrons, i.e., the distance between two
consecutive collisions, is much smaller than that in low pressure conditions. The electrons
can not resonate in the microwave electric ﬁeld long enough to gain energy between
collisions. Therefore, under moderate pressure conditions, the plasma is usually
unmagnetized, and created mainly by Joule heating. During the Joule heating process, the
electrons are accelerated by the input microwave electric ﬁeld, and a net transfer of energy
to the electrons occurs as the electrons undergo elastic collisions with heavier particles
such as ions and neutrals. In particular, the direction of electron motion changes during
elastic collisions, causing a net transfer of energy from the oscillating microwave ﬁelds to
the electron. This energetic electron gas can further transfer energy to heavier gas particles
by inelastic collision processes, such as gas ionization and dissociation. Thus, the input

gas is partially ionized and a discharge inside a plasma reactor is created.

3

Once ionization and dissociation occurs, the energy loss mechanism of discharges
depends on a number of factors including a key factor which is pressure. When the
pressure is low, the collisions between particles are infrequent, so electrons, ions, and
dissociated atoms will diffuse to the boundary wall where wall recombination takes place.
Hence, the discharge spreads out and ﬁlls the whole chamber, and the main energy loss is
due to diffusion processes. As the pressure is increased the rate of volume recombination
increases. Electrons, ions, and free radicals recombine quickly and convert the ionization,
dissociation and excitation energy into thermal energy. The discharge then contracts away
from the wall boundary and becomes thermally inhomogenious. Heat conduction then
becomes an important energy loss mechanism of the discharges [2H3].

A key consideration for the operation of microwave plasmas is the impedance
matching between the microwave generator and the applicator which couples the
microwave power into the plasma. An impedance mismatch causes a partial reﬂection of
microwave power by the applicator; hence not all the power delivered by the generator
can be used to sustain the plasma. This behavior essentially arises from the fact that the
impedance of the applicator, and the resonance characteristics of a cavity, strongly depend
on the actual plasma state, e.g., density, temperature and electrical conductivity.
Therefore, in order to understand the coupling efﬁciency from the generator to the
applicator, the characteristics of the resonant cavity, and the interactions between
discharge and resonant microwave ﬁeld, have to be studied in a systematic manner.

In order to understand these behaviors of the plasma inside a cavity reactor source,
microwave ﬁelds, gas transport, energy transport, as well as, their relations should be

investigated carefully. The coupling among these is complex and analytic models are hard

4

to apply to handle these problems. Thus, a systematic numerical model, which can self-
consistently solve the plasma behavior as well as the exciting microwave ﬁelds, is
essential to provide insight, interpret experiments, and help in the design, development,

and control of microwave plasma reactors.

1.2 Objective

The overall objective of this study is to develop a self-consistent numerical model
of a single mode microwave cavity plasma reactor source. This model will provide the
steady state electromagnetic ﬁeld solutions inside cylindrical resonant cavity structures
loaded with H2 discharges, as well as, be used to investigate the excitation mechanism and
characteristics of these discharges across a pressure range (100 mTorr to 100 Torr) and
across an input power range. The overall signiﬁcance of this study is to build up a self-
consistent model of microwave cavity plasma excitation.

The overall objective may be divided into several speciﬁc objectives which
include:

1. To develop a time domain numerical model to simulate the steady state
electromagnetic ﬁelds inside a cylindrical cavity with a speciﬁc resonant mode or input
coupling structure. The cavity can be either empty or loaded with a plasma.

2. To develop a plasma model which can describe the physical behavior and
investigate the characteristics of H2 discharges inside a microwave cavity plasma reactor
source.

3. To couple these two models together to build up a self-consistent microwave

cavity plasma reactor model. The discharge excitation mechanisms and the microwave

5
power absorbed by the discharge will be investigated in detail by this model.

4. To verify that the models developed agree with experimental diagnostic

measurements.

1.3 Outline

This dissertation is organized as shown in Figure 1-1. Chapter 2 presents the ﬁnite
difference time-domain (FDTD) model which is used throughout this study for
electromagnetic ﬁeld simulation. The plasma conductivity models used for determining
the microwave power absorbed by the discharges are developed in both the frequency
domain and the time domain. Chapter 3 reviews the numerical models for plasma
discharge simulation and presents the ﬂuid plasma model used in this study for microwave
discharge simulation. The physical background and numerical techniques used for model
development are discussed. Chapter 4 provides the simulation results of transient state
response of empty and loaded resonant cavities. The FDTD model is applied to investigate
the characteristics of resonant cavities, such as the natural frequency response, cavity
quality factor Q, and the electric ﬁeld patterns. The resonant cavities which are
investigated in this chapter include empty cavities operating in various resonant modes,
cavities loaded with lossy materials and argon discharges, and a compact ECR ion source.
Chapter 5 presents the simulation results of the electromagnetic behavior of a steady-state
microwave cavity plasma reactor used for diamond ﬁlm deposition. The method used to
model the input power coupling is included in the FDTD model. Plasma discharges used
here are assumed spatially uniform and with a constant temperature. Electric ﬁeld and

power absorption patterns under various input conductions are presented. Chapter 6

 

 

C N r \

 

 

Electromagnetic Field Model Plasma Model
Mgggl Dgyelgpmgn; Model Develo ment
(Chapter 2) Absorbed (Chapter 3)
Power

0 FDTD model F_. 0 Fluid equations
0 Plasma conductivity model " $923351)”; Sﬂfaﬁons
-- Frequency domain ’

__ Time domain -- Energy equation

 

 

 

 

 

 

 

 

(electron)
Electron -- Poisson equation
Tnin lin Temperature 0Nm '1 thd
_ta_§_e_t_So_ut_tL l u errca me 0
(Chapter 4) -- Newton’s method
Plasma
° -- Iterative method
0 Empty cavity Density
0 Loaded cavity \ J
-- Lossy material
-- Ar discharge l

 

o ECR compact ion source

 

 

Self-Consistent Solution

 

 

Chapter 6
W 0 Microwave cavity plasma
(Chapter 5) reactor loaded with H2 gas

. -- Electromagnetic behavior
0 Input power coupling

0 Uniform H2 discharges

k j

-- Discharge characteristics
-- Q value measurement
-- Model veriﬁcation

 

 

 

 

 

 

 

 

Figure 1-1: Overview diagram.

7

couples the FDTD electromagnetic ﬁeld model with the ﬂuid plasma model to provide a
self-consistent microwave cavity plasma reactor model. The characteristics of discharges,
such as densities and electron temperature, are investigated under various input
conditions. The simulated cavity quality factor Q is compared with experimental results.
Chapter 7 presents the conclusions and some speculations on future research and

numerical model development for microwave cavity plasma reactor simulation.

Chapter 2
Electromagnetic Field Model

2.1 Introduction

Microwave ﬁelds serve as the energy source to excite and sustain the discharge
inside cavity plasma reactors. The microwave energy is imparted into the cavity as shown
in Figure 2-1 through an input electromagnetic coupling structure, such as a power
coupling probe, and conﬁned by the cavity walls which serve as the electromagnetic
boundaries. When the exciting wave frequency and the cavity size match a certain
resonant condition, the electromagnetic ﬁeld will oscillate back and forth inside the cavity.
These oscillating electromagnetic ﬁelds usually have speciﬁc ﬁeld distribution patterns, or
so called resonant modes, so the electromagnetic energy can focus on some local regions
inside the cavity and breakdown the gases. These ionized gas species form the plasma
discharges which are usually conﬁned in a container, such as a quartz chamber. For an
empty lossless cylindrical cavity, the electric ﬁeld pattern, or the resonant cavity mode is
well characterized by the size and resonant frequency. When a plasma is present, the
plasma acts as a volume of lossy material loaded inside the cavity. It changes the resonant

electric ﬁeld patterns and the characteristics of the resonant cavity because microwave

Electromagnetic wave

111

Input power I H l

coupling m

Cavity walls

/

 

 

 

 

Quartz chamber

\\

 

Plasma

 

 

 

 

 

 

 

Figure 2-1: Schematic description of microwave discharge generation.

10

power is absorbed by the plasma. Moreover, since the plasma discharges are generated by
the resonant electromagnetic ﬁelds, changes of the resonant behavior of the cavity also
affect the characteristics of the discharges.

Since the resonant electromagnetic ﬁelds play a vital role for microwave discharge
characterization, a mathematical model is required to describe the electromagnetic ﬁeld
behavior inside the plasma reactor system. In this chapter, the history of electromagnetic
ﬁeld models, both frequency domain or time domain modules, used for solving the
resonant electromagnetic ﬁelds inside cavities are reviewed. The technique and
background of the ﬁnite-difference time-domain method, which is used in this study to
model the electromagnetic ﬁelds, are examined in detail. The method to implement the
FDTD model in cylindrical coordinate is next described. Finally, the techniques used to
solve electromagnetic ﬁeld interactions with lossy materials and plasma discharges by
coupling the FDTD model with appropriate conductivity models and by simultaneously
solving the Maxwell’s equations with the electron momentum transport equations are

discussed.

2.2 Review of electromagnetic ﬁeld model

One way to model the electromagnetic ﬁelds in a system is directly solving the
Maxwell equations, which govern all well-behaved electromagnetic quantities. The

Maxwell equations are

11

BI!

VXE = — E
3E
VXH — 83—" +1 (21)
VoE = E
E
V-H = 0

where E is the electric ﬁeld strength (Volts per meter), H is the magnetic ﬁeld strength
(Amperes per meter), J is the electric current density (Amperes per square meter), p is the
electric charge density (Coulombs per cubic meter), 8 is the electric permittivity of the
medium (Coulombs per Volt-meter), and it is the magnetic permeability of the medium
(Volt-seconds per Ampere-meter). This differential form of Maxwell’s equations can be

solved either in the frequency domain or in the time domain.

2.2.] Frequency domain

To solve the electromagnetic ﬁelds for a resonant cavity in the frequency domain,
the electromagnetic quantities such as E and H are often assumed to be time-harmonic,
i.e., their time dependence can be described by a periodic sinusoid and can be included in
Maxwell’s equations as a factor of the form e’m. Then, the Maxwell equations are
converted to vector Helmholtz equations, and the E-ﬁeld and the H-ﬁeld are analytically
solved in cylindrical coordinates in terms of a set of Bessel function type solutions [49].
The problem becomes an eigenvalue problem and the electromagnetic ﬁeld solution inside
the cavity can be expressed in terms of the eigenmode solutions such as the transverse
magnetic mode (TM mode) and transverse electric mode (TE mode) solutions[49]. For a

loaded cavity, i.e., different mediums ﬁlled inside the cavity, the electromagnetic ﬁelds

12

can be expanded in terms of a series of eigenmode solutions in each medium region, and
appropriate boundary conditions are applied to ensure the electromagnetic continuity
conditions are satisﬁed at the boundaries. A characteristic equation then be formed to
determine the unknowns of the ﬁeld solutions and can be solved either by graphical
techniques or by numerically searching for the eigenvalues using a complex root ﬁnding
method. This so called mode-matching or mode-expansion technique of solving time
harmonic electromagnetic ﬁeld solutions has been well established for studying the
electromagnetic ﬁeld behaviors inside coaxially loaded cylindrical waveguides and
cylindrical resonant cavities loaded with lossy materials [8].

In terms of microwave cavity plasma reactors, the same techniques can be applied
to investigate the microwave and plasma behavior by treating the discharges as lossy
materials. J. R. Roger [6] used this technique to solve the electromagnetic behavior inside
a cavity containing a quartz-tube-conﬁned Ar plasma by treating the plasma as a lossy,
homogeneous, isotrOpic, rod-shaped material. M. L. Passow, et a1. [7] used the same
electromagnetic analysis to determine the electron density, conductivity and collision
frequency of a tube-contained air discharge. Modeling the inhomogeneous properties of
plasmas using the mode-matching technique has been investigated by S. Offermann
[9][10]. He coupled the Maxwell’s equations together with an energy balance equation to
solve the local dielectric constants of a high pressure Hg plasma. In order to reduce the
complexity of this problem, only TMOmO modes were considered, which reduced the ﬁeld
conﬁguration to one dimension only.

Generally, there are two major difﬁculties in applying the frequency domain

electromagnetic ﬁeld solutions and mode-matching methods to solve the plasma-loaded

13

cavity problem. First, the plasma discharges inside microwave cavity reactors are
inhomogeneous, anisotropic, and dispersive (the electric properties dependent on
excitation frequency) in nature. This behavior makes the mode-matching method hard to
apply to solve the plasma-loaded cavity problem. Another problem results from the
geometry of the lossy load (discharge) or resonant cavity. If the load or cavity are
geometrically complex, the number of boundary conditions needed for ﬁeld solutions

increases leading to intractable computer solutions.

2.2.2 Time domain

An alternative way to solve the electromagnetic ﬁeld inside the cavity is to directly
solve the time-dependent Maxwell’s equations by using a ﬁnite-difference time-domain
(FDTD) method. The FDTD technique was ﬁrst proposed by K. S. Yee in 1966 to solve
the interactions of electromagnetic waves with perfectly conducting material [11]. In his
model, the time-dependent Maxwell’s equations in a rectangular coordinate system were

expanded to a set of six scalar equations:

29H,r [[35y 352)

6H, _ [[352 35,) 23
a “aa'a ‘”
8H, [[35, 35,] 24
"a? " a 'a? ‘3: ‘"

IO

14

as, {a}!z arty] II 25
'57 "as; ‘5: "g H
BE 13H 3H J

y _. _ I___2 -_y
'37 ‘13". 8x] 8 _ (2'6)
at:z 1 aHy 3ij J2
79': ‘42: "a3: '5 (2'7)

These scalar equations were further discretized in both the time and space domain

get a set of ﬁnite difference equations:

 

 

 

 

 

 

n+l/2 . . '1 l _ n-1/2 . . I 1 AI
Hx (I,j+§,k+§)—Hx (t,1+2,k+2)+ (i 4.1/(+1)
u ” 2’ 2
..1 n..l n.. 1 n.. 1 (2-8)
X E;(I,j+§,k+1)-Ey(l,j+-2-,k)-Ez(l,j+1,k+§)—Ez(l,j,k+§)
A2 A)’
H"+1/2(i+l,j,k+l)= H"_1/2(i+-1-,j,k+l)+ At
2 2 2 y 2 2 “(H1 .k+1)
2”’ 2
. .1 n.. 1 .1. ".1. (2-9)
x E?(I+l,j,k+§)—Ez(l,j,k+§)—E:(I+-i,j,k+1)-Ex(l+§,j,k)
Ax A2
n+1/2 . I. 1 _ n—l/2. 1.1 At
Hz (I+§,j+§,k)—Hz (1+2,j+2,k)+ (54.14.10
“ 2” 2’
.1. ".1. . .1 n..1 (2-10)
X E:(I+§,j+l,k)-Ey(l+§,j,k)—E;(I+1,j+§,k)-Ez(l,j+§,k)

 

 

Ay Ax

15

Ef"(t+1,j,k) = E:(i+l,j,k)+ A‘

 

 

 

 

 

 

 

 

 

 

 

 

 

2 2 e(i+%,j,k)
f n+1/2 . 1 . 1 ) n+1/2(. l . 1 ) \
Hz (1+2,j+2,k -Hz 1+2,j—§,k (2.11)
x Ay
H’””2(t+l,j,k+.1-)-H’””2(i+1,131:—1)
y 2 2 y 2 2 —./n(i+l .k)
K Az " 2’], /
£""(t,j+l,k) = E"(i,j+l,k)+ A‘
y 2 y 2 (. . 1 )
e l,j+-,k
2
( n+1/2(.. 1 1)_ n+l/2(.. 1 _1) )
Hz l,j+2,k+2 HJr I,j+2,k 2 (2.12)
x Az
H"+1/2(i+l,j+-1-,k)—H"+1/2(i—l,j+-l—,k)
z 2 2 z 2 2 _Jn(i 411:)
\ Ax y ’1 2’ J
£’”‘(i.j,k+1) = E"(i.j,k+-‘-)+ A’
z 2 z 2 £(i 'k-t-l)
,1, 2
( n+1/2(. _1_ . 1)_ n+1/2(._1 . 1) )
H, ”2””“2 ”2 ‘ 2”"“2 ' (2.13)
Ax
x
H"+1/2(i,j+l,k+l)-H"+U2(i,j—l,k+1)
1‘ , 2 2 x 2 2 4,10 .k+1)
K Ay ’ ’j’ 2)

 

 

where (i, j, k) denotes a grid point in space and n denotes the time step, Ax, Ay, and A2 are
the space grid unit in x, y, and 2 directions, and At is the time step. These ﬁeld components
in a unit cell are assigned in an interleaved manner similar to that shown in Figure 22.
With the system of equations, (2.8) to (2.13), the new value of a ﬁeld vector component at

any grid point (e.g. Ex"), depends only on its previous value (e. g. E X“) and the value of

16

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2 ..f
‘ +l’y/
(E y) ,.
f E (Lt-“2, j; k-r-l)
tamer/2, lit/5F) ,,/
(Hx),» /./'/ rdO
t2!” . ..... Hli¢112,.jf.kw1ﬁ) /(y)
E,( j, k+1/2) (Pry) f,
./""f
( ) fr Hz(i+l,|i+1/, )
./"
..| ’ r
(i, j, k) E (HI/2,13 k)
(ﬁx) (X)
2
(y)
6
> '0‘)

Figure 2-2: A FDTD unit cell in three-dimensional cylindrical coordinates.

17

the components of the other ﬁeld vector at the adjacent point, which are perpendicular to it
(e.g. Hymn/2 and Hle/Z). Therefore, with appropriate initial and boundary conditions,
the electromagnetic ﬁelds on the grid points of the space mesh can be evaluated in time
domain at alternate half-time steps (leapfrog method). By the stagger mesh arrangement in
space and by the leapfrog time marching method, as shown in Figure 2-3, the numerical
solutions of the FDTD model can achieve second order accuracy both in time and space
[55].

This method has been used extensively to investigate the electromagnetic wave
scattering phenomenon with dielectric or conducting objects. R. Holland, et. al.,
developed a three-dimensional computer code based on the FDTD method, to solve the
electromagnetic pulse ﬁeld scattering and surface currents for an aircraft [12]. Steady-
state, sinusoidal electromagnetic scattering and penetration problems were investigated by
A. Taﬂove, et. al., using the FDTD method for arbitrary metal or dielectric structures
[13][14]. The electromagnetic power deposition in a human body has also been simulated
by the FDTD method [17][18]. The model of the human body is an anatomically based
model, which is an inhomogeneous, human shaped structure and composed of different
types of tissues. Each tissue is assigned a set of speciﬁc dielectric constant and
conductivity values, which can be used to determine the electromagnetic power absorbed
by the tissue [19].

By using the Fourier Transform technique, the results of the FDTD method can be
used to calculate the parameters in the frequency domain. This technique has been applied
to analyze the characteristics of microsuips, waveguides, as well as, resonant cavities

[201-[26]. X. Zhang et. al. [20][21] use the FDTD method with Gaussian pulse excitation,

18

o -- E ﬁeld

0 ---~ H ﬁeld
Time

n+1 O
A

n+1/2 o ---------- :-\-O

 

 

n-l/2 O O

 

 

’ Space
i-l (i-1/2) i (i+1/2) i+l

Figure 2-3: Leap-frog scheme.

19

which provides a wide band of frequency, to investigate the dispersive characteristics and
discontinuity of microstrips. In terms of resonant cavities, D. H. Choi, et al., used the
FDTD method and discrete Fourier transformation (DFI') techniques to solve three
dimensional eigenvalue problems of inhomogeneous rectangular resonators [24]. The
steady-state FDTD solution was treated as a time-harmonic function, from which the
eigenvalues (resonate frequencies) were extracted by a discrete Fourier Transform

technique given by
s m = 25" (10,1'0, k0) exp (—j21tsnf) (2.14)
n

where S is the frequency (f) response spectrum, F is the FDTD solution at grid point (10, jO,
k0) of the 11‘h iteration, and s is the stability factor. The same method was also used by A.
Navarro, et. al., [25] to calculate the resonant frequencies of TEO and TMO mode
cylindrical cavities with a axially symmetric dielectric load at the base. The FDTD
method was also applied to simulate the microwave interactions with materials inside a
resonant cavity. M. F. Iskander, et al., [26] used the FDTD method to simulate the
sintering process of ceramics inside a rectangular, single mode microwave cavity. They
provided a model which utilized a realistic wave excitation arrangement to improve the
calculation of the cavity Q and the shift in resonant frequency. The steady state spatial
electric ﬁeld distributions were solved for both empty and loaded resonant cavities. The
power absorption for different materials loaded inside the resonant cavity under different
operating conditions was also calculated and analyzed.

There are several advantages of using the FDTD method to simulate the

20

electromagnetic ﬁeld behavior for a plasma loaded resonant cavity. First, it is easy to
implement for complicated loads (like plasmas), because spatially varying dielectric
parameters can be assigned using grid points. Second, complex simulation structure, such
as the discharge or cavity reactor geometry, can be solved by the FDTD method by using
appropriate grid structures. Moreover, most of the plasma numerical models, both particle
and ﬂuid models, can be easy coupled to the FDTD model because both are in the time
domain. Therefore, the FDTD method has great potential to solve the electromagnetic

ﬁelds inside a discharge loaded cavity.

2.3 FDTD model implementation in cylindrical coordinates

As mentioned above, the ﬁnite-difference time-domain (FDTD) method is a good
choice for solving the electromagnetic ﬁelds inside a plasma loaded cavity. The FDTD
method used in this dissertation is adopted from Yee’s algorithm [11] and transferred into
three-dimensional cylindrical coordinates. In cylindrical coordinates, (2.1) can be

expanded as six scalar equations given by

BHr ltagq, 13132]
u

a -.-.- 5.2. —;§$ (2.15)
31¢ _—. l[a_EZ_a_E'] (216)
a: [.1 Br 82 '

3H BE
2 _ l l r_la
— [r375 r57"E¢’] ‘2'”)

21

fr _ l[l?:12_a_H¢]-:’ (218)
at - 8 r34) 32 8 '
95¢ =l[31'_a_HZJ_J_¢ (219)
a: e 82 ar e '
8E 8H J

z - .1. 13 _l_' __2
'67 - el:r-a7(rH¢) ra¢ ] 8 ° (220)

These six differential equations can be discretized in both time and space domain by the
central ﬁnite difference method. The resulting ﬁnite difference equations in cylindrical

coordinate are:

 

 

 

 

 

 

+1/2 . 1 _Hn- 1/2(l.,1+. l I)+ At
H" ("1 2””2)‘ r 2"”2 11(1 .+1k+1)
a} 2: 2
l ( 151%,“)
E" +% ,k+1 ,
"j H (2.21)
x
E:(i,j+l, k+- :)- E:(i,j,k+- 2)
_ ’1(A¢) _
+1/2 . _l_ . 1) 1/2(i+1 1)+ 1At
H; (z+2.1.k+2 H; +2.j.k+

11:

mg)

(”i
-E:('+l’j’k+l 2A) E:("2j’k+ 2) (2.22)

 

E:(t+§,j.k+1)-E’,‘(i+1-,j.k)
- Az ..

 

 

 

22

 

 

 

 

 

 

 

 

+r/2._1_.1)_ —1/2(.l._1_) 2At
H: (1+2, +2,k —H: 1+2,]+2,k + (i+.l_ 4.1-[(+1)
’1 2’] 2’
_ £"(i+1j+1 k)-E"(i+lj k)-
rt+1"rix r 2’ ’ ’ 2”
ri2+l—ri2 A¢
X l 1
n . . . .
ri+lE¢(l+l,j+§,k)-riE;(l,j-l-§,k)
2 2
L ’t+1"’1
a:+1(~2m)=E:-(~1m+ )
€(l+§,j,k
(+1/2.l.1 +1/2.1.l
”3 (”5.15.042 (”24791
'1+1/2A¢’
x Hn+l/2.l. 1 Hn+l/2.l. 1
q: I+§,j,k+§ — ¢ l+§,],k-§
Az
Jn. 1 'k
K _r l+§,], )
n+1 . . l _ n . . I At
Ed: (l,j+§,k)—E¢(l,j+§,k)+ .. 1
8(l,j+2,k)

 

 

 

 

l 1 +1/2 1
+§,k+§)-H: (a +§rk
A2
1 1 +1/2 . 1 .
2 ski-”3 ("5’1"

+1/2 . .
x H: (l+-,j+

1
2’

k)

 

 

’i+r/2"1—1/2

. . l
"13("H2’kl

 

(2.23)

(2.24)

(2.25)

23

E""(i,j,k+l)=E:(i,j,k+1)+ A’ 1
€(i,j,k+§)

 

Z 2 2

 

 

 

( n+l/2 . I . 1 n+l/2 . 1 . 1 j
ri+l/2H¢ (r+§,j,k+§)—ri_l/2H¢ (l-§,],k+§)
2 2
rt+1/2"1—1/2 (2.26)
+1/2 . . 1 1 +1/2 .. 1 l
X 1, +-ak+-)-’Hn (I, "",k+-)
ri+1/2_ri-1/2x ’ (J 2 2 ’ j 2 2
,2 Ar
1+1/2 i-1/2
.. 1
\ —J:("J’k+2) /

 

 

where n is the time index, i, j, k are space indexes, At is the time step, Ar, At]; and A2 are
the dimensions of the unit cell in r, t1) and z direction, and rindex is position in r direction.
Similar to rectangular coordinates, the six ﬁeld locations in cylindrical coordinates are
assigned to be interleaved in space to satisfy the continuity of tangential ﬁeld components,
as shown in Figure 2-2. If the dielectric parameters, such as permittivity e and
permeability rt, are space dependent, appropriate values are assigned to lattice points for
each ﬁeld component. Therefore, inhomogenous, anisotropic, or arbitrary shape materials
can be easily handled by the FDTD model.

If the simulation domain is (1) symmetry, such as the TMO mode resonant cavity
problem, the three dimensional FDTD model can be reduced to a two dimensional model
by neglecting all the variation in the 4) direction, and the index (i, j, k) in (2.21) to (2.26)
reduce to (i, k) only. The discretized two dimensional Maxwell’s equations by assuming 4)

symmetry can be rewritten as:

24

n+l/2 . n-
”r (bk-+5) = Hr l/2(j,k+l)
2
X_ At l)[5l(i,k+l) -E;(r,k)] (2.27)

',k _ AZ
“(' +2

+1/2 . 1 1 n-
”; (I+§,k+§)=l‘l¢1/2(i+%k+l)+_ At

E:(z+1,t+;)-.»(. 1.2) we
2 ’ 2

—k)+ 2At
2’ (+1 k
l _

x[ ri+lE;(i+1,k) —riEg(i,k)] (2‘29)

1+1—

 

"1 . 1 n
E,+ (z+§,k)=Er(i+%,k)+ At
e(i+

”n+1/2(i+ll)_n+l/2 1 2
,, +-2-,k+i 11¢ (1+ “7) (.30)

 

En+1(i,k) =En(ik)+ £(At

k)
( H"+l/2(ik+ 14]) ”n+1/2( 1 \
“k‘§)

A2 (2.31)

H"+1/2(i+l )_ n+1/2, 1
2 +2" HZ l-§,k)

rt+1/2"t—1/2

 

 

44(tnj

25

+1 . l n . 1 At
E: (l’k+§)=Ez("k+§)+——(_ I)
8 l,k+-2-

( n+l/2 . 1 1 n+1/2 . 1 I \
’1+1/2H¢ ('+§’k+§)‘rt—1/2H¢ ("§’k+§ (2.32)

2 2
’t+1/2"t—r/2

[ —J:(i,k+%) J

 

 

 

The time step in (2.21) to (2.26) and (2.27) to (2.32) is determined by the cell size

and must satisfy the stability condition [13]:

 

5min . .
At 5 for three drrnensrons (2.33)

vme

and

At S m" for two dimensions (2.34)

vmax J5

 

where 8m." is the smallest space grid size, and vmaJlr is the maximum phase velocity within

the model. vm is typically selected as the speed of light, c.

2.4 Conductivity models for FDTD model implementation

For the resonant cavity simulation, the current density (J, 1,1,, Jz) in equations
(2.24) to (2.26) represents the current density which is induced by the applied
electromagnetic ﬁelds. It plays an important role in the FDTD model because it describes
the electromagnetic power absorbed by the loaded material or cavity wall. For an empty or

perfect dielectric (loseless) medium loaded cavity with perfect conducting walls, the

26

current density is zero and there is no power loss in the cavity. When losses occur due to
the presence of a lossy medium, such as a plasma load, the power density P absorbed by

the medium is:

P(r, t) = J(r,r) 0E(r, t) . (2.35)

The current density J can be determined by using an appropriate conductivity
model for lossy materials or plasmas. The plasma conductivity model can be established
in the frequency (time—harmonic) domain as a ﬁxed conductivity or in the time domain by
solving the momentum transport equation of electrons, where J = env. Additionally, the
current density J also can be calculated by averaging over the particle velocities in a

particle type plasma simulation.

2.4.1 Conductivity model for non-dispersive lossy materials

When electromagnetic ﬁelds interact with a non-dispersive (the dielectric
parameters are independent of frequency) lossy material, the current density (J ,, J¢, J2) in

(2.24) to (2.26) can be evaluated using the conductivity of the load, i.e.

J,,,(i.j.k) = o<i.j.k)E,,(i.j.k), m=r,¢,z (2.36)

where i, j, and k denote the grid indices in the r, (p, and z directions, respectively. In this
case, like 6 and 1.1, o can be assigned to each lattice point to calculate the current density.
By substituting (2.36) into (2.24) to (2.26), the equation in center ﬁnite difference form for

E ,, for example, can be rewritten as

27

Ef”(i+%,j,k) = CA,(i+%,j,k)Ef(i+%,j,k)+ Al’
£(i+ §,j,k)

fH:+l/2(i+%,j+%,k)—Hn+l/2(i+-1' '—-1-k) \

 

 

 

 

 

 

 

’ 2” 2’ (2.37)
x ’i+1/2A¢
H n+l/2(i+l . k+l)_Hn+l/2(i+_1_. k-l)
4* 2’1’ 2 e 2”’ 2
\- Az )
where
. l . . l . -l
1 Orl+§,j,k AI Or l+§,j,k A!
CA,(i+-,j,k) = 1— x 1+ (2.38)
2 . l . . l .
2811 + 2’1’ k) Ze,(r + 5,11)

Similarly, E,» and E2 can be rewritten in this manner.

2.4.2 Plasma conductivity model

Plasma discharges can also be viewed as lossy materials, and microwave power
absorbed by the plasma can be described using a plasma conductivity model. The
conductivity of a plasma can be derived by solving the macroscopic motion of the
electrons under the inﬂuence of applied forces and particle collision processes. The
microscopic electron motion can be described by the electron momentum transpon
equation, i.e., the Langevin equation [50] under the assumption of the cold plasma:

me%v(r,t) = —e(E(r, t) + v(r, t) xB (r, t)) —meveff(r, t) v(r, t) (2.39)

where v is the average electron velocity, E is the applied electric ﬁeld, B is the applied

28

magnetic ﬁeld, q is the electron charge, me is the electron mass, and veﬂ is the effective
collision frequency which represents the elastic collision frequency between the electron
and neutrals for momentum transfer. This equation assumes that the electrons are free to
move in a stationary uniform background of ions and neutral which provide a viscous
damping force to the motion of the electrons. It is suitable. for partially ionized and elastic
collision dominate discharges such as those of interest in this study.

From (2.39) the current density can be determined as

J (r, t) = —ene (r, t) v (r, t) (2.40)

where ne is the electron density. By equation (2.39) and (2.40), the characteristics of the
plasma (tie and veﬂ) can be coupled with the J quantity in FDTD model to investigate the
electromagnetic excitation of discharges inside microwave cavity plasma reactor.
Equation (2.39) can be solved and modeled in either the time harmonic domain or the time

domain.

2.4.2.1 Time harmonic domain solution

In the time harmonic domain, the applied electric ﬁeld is considered oscillating
harmonically with time, E = Eoeimt. For unmagnetized plasmas, the steady state

solution for v must have this time dependence given from (2.39)
jwmev = — eE — mvveff (2.41)

therefore,

29

 

 

 

 

 

_ _ v
v = m (v at: 'w) = :1 2 eff 2”j'§—'2']' (2‘42)
e eff 1 e veff+ (n veff+ 0)
From (2.40), the current density is
J = CE = —enev (2°43)
and the conductivity of an unmagnitized plasma is
ezn v
o: m‘[ 2 eff 2—j 2‘” 2]. (2.44)
e veff+w veff+w
Also, the dielectric parameter of an unmagnitized plasma can be written as
e = 808,
2 2
(o v (o (2.45)
£,=1_j£= 1___L_ +j_€[. P
(”80 v2 + m2 (n v2 + (02
eff eff
where (0,, is the plasma frequency, which is deﬁned by [50]
2
nee
(up = . (2.46)
"1:80

The time-average electromagnetic power absorbed by the plasma can be expressed

as[1]

<P>,,,,(r) = §Rewm . (o(r)E(r))*) (2.47)

30

where P has the units of power density (Watt/m3). Substituting equations (2.43) and (2.44)
into equation (2.47), the time-average-absorbed power density of an unmagnitized plasma

is given by

 

n e2 ve
(”abs”) = 2:" [ 2 ff ZJIEUHZ (2.48)
e veff+w

where lE(r)| is the magnitude of the electric ﬁeld as a function of r. From above, the power
absorbed by the discharge is proportional to the real part of the plasma conductivity, or the
imaginary part of the plasma dielectric constant.

For magnetized plasmas, the conductivity is in a complex tensor form, which
indicates the anisotropic properties of a magnetized plasma. By assuming the magnetic
ﬁeld B exits only in the z direction, i.e., B = B02, the relation of the current density to the

applied electric ﬁeld can be expressed as [50]:

 

r '- r 1 r -
J)r ol —0x 0 EJr
J, = ox <3l 0 Ey (2.49)
J O O 0,, E2

 

 

 

 

 

where

31

n eez veff + jw

 

 

 

OJ- - m . 2 2
e (Veff-l-jCD) +0)“
2
n e w
Ox = [:1 . cez 2 (2.50)
e (veg-um) +0)“
nee2 1
011 = -
me Veff-l-jO)

a)“ is deﬁned as the electron cyclotron frequency, which is equal to eBO/ m 9.
Similarly, by substituting equations (2.49) and (2.50) into equation (2.47), the

time-averaged power density absorbed by the magnetized plasma becomes:
1 2 2
(P)abs (r) = 5(Re (oi) |Ei (r)| + Re (on) '511 (r)| ) (2.51)

where the symbol J. means perpendicular to the static magnetic ﬁeld (B2), and // means
parallel to the static magnetic ﬁeld (B2). If the time-varying electric ﬁeld is perpendicular
to the static magnetic ﬁeld only, the time-averaged power density absorbed by the plasma

becomes:

 

 

2v
(P)a,,_,.(r) = "‘8 e” 1 + 1 )lEml2 (2.52)

4me (Vesz'i' (to—(ow)2 (Vesz+ (w+coce)2)
Note that when to >> veﬁr and race ~ to, (2.52) represents the electron cyclotron heating
mechanism.

The above conductivity expressions can be coupled with the FDTD method by
assigning the calculated conductivity value at each grid point. This method can

approximately provide the power absorption solutions of the plasmas. However, because

32

plasma discharges are dispersive in natural, i.e., its macroscopic dielectric parameters are
dependent on its excitation frequency, these conductivity models expressed in frequency
domain are highly dependent on the real exciting frequency. This will cause difﬁculties
while applying them in time-domain simulations (such as the FDTD model) where the
time variation of the electric ﬁelds may not be known before the simulation begins [16]. In
this study, in order to not lose generality, the Maxwell’s equations are solved self-
consistently in the time domain with the electron momentum transport equation (2.39).

The details are provided in the next section.

2.4.2.2 Time domain solution

To model the plasma conductivity due to an applied electric ﬁeld in the time
domain, 8 and u in the discretized Maxwell’s equations (2.21) to (2.26) are set equal to £0
and #0 (dielectric parameters in vacuum or air), and the current density in equations (2.24)
to (2.26) is evaluated using the electron momentum transport equation (2.40). To get the
time domain solution of equation (2.40), the ﬁnite difference method is used to discretize
this ﬁrst order differential equation in time. Neglecting the static magnetic ﬁeld, (2.40)

becomes

n .. .. n—l.. Arn-r..
vr,¢,z(z,j,k) = (1.0—At-veff(z,1,k)) -v,,¢,z(z,),k)—%n: “wow/c) (2.53)

and the current density is

J”,_,(i.j, k) = —qn,(i.j. k) vf,,,,(i.j, k) (2.54)

r,

33

where n denotes the time iteration count, At is the time step, and i, j, k denote the grid
location. Equations (2.53) and (2.54) are solved simultaneously with equations (2.21) to
(2.26). The absorbed power density can be solved directly by coupling (2.53) with (2.35).

From the above equations, the space-dependent current density J is determined by
the space-dependent discharge characteristics such as electron density ne and effective
collision frequency veﬁr. Moreover, the effective collision frequency which describes the
momentum transfer due to electron-neutral elastic collisions is a function of electron
energy distributions (electron temperature). For a complete microwave plasma excitation
solution, the electron density and electron temperature must be solved self-consistently by
using appropriate plasma models, such as ﬂuid (continuum) or particle plasma models.
These models used for solving the characteristics of plasma discharges are discussed in

the next chapter.

2.5 Summary

The electromagnetic ﬁeld model for this study was developed in this chapter. The
model developed solves Maxwell’s equations in two and three dimensional cylindrical
coordinates using a FDTD technique. The microwave induced current density which
determines the power absorption inside the cavity was described by using conductivity
models. The methods developed to describe the conductivity of discharges included the
solution of the electron momentum transport equation either in the frequency domain or in
the time domain. This electromagnetic ﬁeld model will be veriﬁed by the transient state

simulation results of empty and loaded cavities in Chapter 4.

Chapter 3

Numerical Models of Plasma Discharges

3.1 Introduction

In this chapter, numerical models which describe the behavior of plasma
discharges are presented and developed. First, the models used for plasma simulation in
the literature are reviewed and the physical background of discharge transport behavior is
discussed. Next, the ﬂuid description of the plasma based on the moments of Boltzmann
transport equation is derived. Finally, the numerical methods used for solving the ﬂuid

model description are introduced.

3.2 Review of plasma discharge models

Generally, for plasma-aided processing simulation, there are two major approaches
to solve the transport behavior of discharges. One is the particle approach, which is done
using the particle simulation technique which treats the plasma as a combination of
particles (electrons, ions, and neutrals) and calculates the trajectories of many particles in
time to obtain the macroscopic quantities of the plasma. The other one is the continuum
(or ﬂuid) approach, which neats the plasma as a ﬂuid and solves the moment equations of

the Boltzmann transport equation (BTE) which include the continuity equation,

34

35

momentum transport equation and energy transport equation of each plasma species by
appropriate numerical techniques. The ﬂuid models and particle models are sometimes
combined together as a hybrid ﬂuid/particle model to provide self-consistent solutions.

These models are reviewed in this section.

3.2.1 Particle model

The particle models simulate the trajectories of particles inside the discharge as
they move through a system under the inﬂuence of the applied ﬁelds and random
scattering forces. If the number of simulated trajectories is large enough, the average
results provide a good approximation of the particle transport behavior inside the
discharges. There are two types of particle models used in plasma processing simulation.
The ﬁrst is a single particle Monte Carlo method, which precisely simulates a single
particle trajectory and its collision events. It is usually used to solve the transport
parameters such as mobility, diffusion coefﬁcient and various collision rates of the plasma
species. M. J. Kushner [31] and B. E. Thompson, et. al. [32] studied the secondary
electrons phenomena and ion bombardment properties respectively in RF discharges by
the Monte Carlo techniques. A similar technique was developed by J. I. Ulacia and J. P.
Mchttie [33], to perform a two-dimensional plasma etching simulation. The motion of
electrons in applied magnetic ﬁelds has also been investigated using Monte Carlo methods
by G. R. Govinda Raju [34] and T. E. Sheridan, et. al. [35].

The major problem of the single particle Monte Carlo methods for plasma
simulation is the inclusion of the local electric ﬁeld induced by charge densities variations.

To overcome the above problem, a many-particle simulation method, particle-in-cell

36

method (PIC), is introduced in combination with the Monte Carlo collision method [47].
Originally, the PIC method was used in hot fully ionized collisionless bulk plasmas
(especially for fusion plasmas) [36]. The deﬁning characteristic of the PIC simulations is
the method of calculating the force acting on each particle. The simulation region is
divided into a number of cells and the resulting grid is used in the solution of a ﬁeld
equation (e.g. Poisson’s equation) from which the force on each particle can be
determined. Each simulation ‘particle’, also called a superparticle, is actually a group of
charges with a speciﬁc charge density. The history of the development of the PIC method
and the details of melding of PIC with Monte Carlo collision models was reviewed by C.
K. Birdsall [36].

Recently, the PIC method merged with the Monte Carlo collision method has been
used to model bounded, partially ionized, collisional discharges, which are used
intensively in industrial applications. D. Vender and R. Boswell [37], M. Suerndra and D.
B. Graves [38], H. W. Trombley et. al. [39], and R. W. Boswell and D. Vender [40] have
used the PIC method together with the Monte Carlo method to investigate RF discharge
behaviors for various operating conditions. T. A. Grotjohn [41] used a PIC/Monte Carlo
particle model coupled with FDTD solutions of Maxwell’s equations as mentioned in
Chapter 2 to model a compact ECR ion source in the 2d3v domain. In this magnetized
discharge simulation, the electron and ion trajectories were determined by the Lorentz
force equation and elastic/inelastic collision events.

In many cases, the particle models are the most accurate techniques available for
analyzing and simulating the plasma transport behaviors. Moreover, it is well suited for

low pressure, ECR, and non-equilibrium discharge simulation. However, its time-

37

consuming nature when implemented in computer programs largely reduces its efﬁciency.
This is especially unattractive for large simulation domains (>> particle mean free path)
and three dimensional (3d3v) simulations. Therefore, in this study, the microwave

discharges are modeled based on ﬂuid techniques.

3.2.2 Continuum (ﬂuid) model

In continuum plasma models, the transport behavior of the plasma species are
described by moment equation solutions of the Boltzmann transport equation, including
the continuity equation, momentum transport equation, and energy transport equation.
These moment equations are often combined with Poisson equation to provide self-
consistent solutions. Under the assumption of a steady state thermal equilibrium existing
in the system, the ﬂuid model might be simpliﬁed to solve only the continuity equation
and Poisson equation and becomes the drift-diffusion model. The details of these models
will be discussed in next section.

Fluid models have been used in the simulation of capacitively coupled DC and RF
discharges. The ﬁrst global self-consistent continuum DC and RF discharge models were
build up by D. Graves and K. F. Jensen [27]. They solved the electron and ion continuity
equations, the Poisson’s equation, and an energy conservation equation for electrons.
Assuming collisionally dominated particle motion, the electron and ion momentum
transport equations were reduced to expressions for particle ﬂux incorporating mobilities
and diffusivities. A ﬁnite-element method with mesh points was used for solving the
equations. D. Graves later coupled an electron excitation term into this model to

qualitatively examine the time and space dependence of rates of electron impact excitation

38

[28]. J. P. Boeuf added negative ions into the continuity equations, and analyzed the
inﬂuence of frequency and gas composition (electropositive v.s. electronegative gas) on
RF glow discharge properties [29]. A local equilibrium approximation was assumed, so
the transport parameters were dependent on the local electric ﬁeld. An implicit ﬁnite
difference scheme was used for the numerical technique. The same ﬂuid model and
assumptions were extended to two-dimensional simulations for N2 and SF6 RF glow
discharges by J. H. Tsai and C. Wu [30]. They used a more accurate ﬂux-corrected
transport method and a reconstructed fast-Fourier—transform technique to solve the
continuity equations and Poisson’s equation respectively.

The transport parameters and collision rates used in ﬂuid plasma model are
sometimes determined by Monte Carlo techniques. M. S. Barnes, et. al. [42] used this
method to solve ion and electron continuity equations for RF glow discharge. In his later
work, Barnes further coupled the Monte Carlo method with three moments of the
Boltzmann transport equation for steady state RF discharge simulation [43]. The
nonequilibrium characteristics of electron transport was taken into account by N. Sato and
H. Tagashira [44] in their ﬂux-corrected transport ﬂuid model. In their study, the Monte
Carlo simulation showed the electron transport properties were not in equilibrium with the
local electrical ﬁeld for SiH4/H2 mixture RF discharges.

For microwave ECR plasmas, due to the low operating pressure (0.1 mTorr to few
mTorr), ﬂuid models may no longer be valid. However, a direct application of a particle
model (like PIC/MC) cost a lot of time for realistic ECR reactors, which normally require
simulation in two (or three) dimensions. R. K. Porteous and D. B. Graves [45] developed a

hybrid electron ﬂuid-particle ion model for magnetically conﬁned low pressure discharges

39

with cylindrical axisymmetry. Ions are treated as particles, which can u'avel in two
dimensions under the combined inﬂuence of applied static magnetic ﬁelds and the self-
consistent electric ﬁelds solved by Poisson’s equation. Non-uniform grids were used
ranging from 0.5 Debye lengths to 300 Debye lengths. Electrons were treated as a ﬂuid,
which were strongly magnetized and had only axial directional motion. A Maxwell-
Boltzmann distribution was assumed and both the continuity and energy equations were
solved for electrons. Power was directly deposited into electrons with an assumed proﬁle.
This work provide the possibility of treating a two-dimensional system which had a very
small Debye length to system chamber length scale ratio.

Y. Weng and M. J. Kushner [46] used a hybrid Monte Carlo- ﬂuid model to
investigate the electron energy distributions in ECR discharges. The Monte Carlo
simulation of the electron swarm was capable of resolving the electron-electron collisions
and provided the details of electron energy distributions. This MC simulation was
iteratively combined with a ﬂuid model, which solved the electron and ion continuity
equations, to calculate ambipolar electric ﬁelds. The electric ﬁeld cycled back to the MC
simulation and was included in the equations of motion for electrons. The microwave
electric ﬁelds were assumed to be a plane wave with uniform amplitude as a function of
radial position. The results of the ﬂuid model for the ambipolar potential, though applied
in low pressure, successfully agreed with experimental results.

Two major difﬁculties which arise while applying a continuum model for plasma
simulation are: (1) If the pressure is low in the system, the collisions of particles are
infrequent. The mean free path for the particles (especially electrons) becomes large, and

in some case larger than the chamber dimensions. It is not suitable to apply ﬂuid models

40

when the mean free paths are not the smallest characteristic length in the system. (2) For
discharges with non-equilibrium or large density (or energy) gradient regions, e.g. sheath
or ECR regions, the transport parameters and rates used in the moments of Boltzmann’s
equation are difﬁcult to determined. The ﬂuid model is inappropriate to apply under these

conditions.

3.3 Model development

The ﬂuid plasma model used in this study is established in this section. The
physical background for model development, such as the concept of distribution function
and Boltzmann transport equations, is discussed. The moments of Boltzmann transport
equations are derived in detail. Finally, the moment equations used in the ﬂuid plasma

model are introduced.

3.3.1 Distribution function

Plasmas consist of several categories of particles, including the electrons, ions, and
neutral species. In order to give a complete speciﬁcation of the properties of this system of
particles, it is necessary to know the phase space location, i.e., the position r and velocity
v, of each particle. The number of particles in an elemental volume of the phase space at
any instant of time can be statistically speciﬁed by a distribution function f (r, v, t) . Then,
the density of particles n (r, t) at a given position inside the system can be determined by

integrating the distribution function over all velocity volume as

n(r,t) = [f(r, v, r) dv. (3.1)

41

Then, the average quantity ¢ of the entire collection of particles is given by

1
n (r, I)

 

(4» = [4» (r. v. t)f(r. v. z) dv . (3.2)

For example, the average velocity (v (r, t)) , and energy (i (r, t)) of all particles are

 

(v(r,t)) = "(1”)Jv(r,t)f(r, v, t) (IV (3.3)
(§(r t)) = —-—1——Ilmv2f(r v t)dv (3 4)
’ n(r,t) 2 ’ ’ '

respectively. Thus, the distribution function f can be described by the macroscopic
quantities (also called moments) deﬁned by (3.1), (3.3) and (3.4) as shown in Figure 3-1.
The area under the distribution function f is described by the density n, the displacement
off described by the average velocity <v>, and the width off is described by the average
energy <§>.

If the system is in thermodynamic equilibrium and not subject to any external

force, the distribution function f becomes the Maxwellian distribution

m )3/2 —-mv2/2kBT

f(r. v. r) = ntr.:)(m 7. e (3.5)
B

which is isotropic in velocity space. Sometimes, it is convenient to express f as a energy

distribution

42

 

 

 

 

0 <v> v

Figure 3-1: The distribution function f can be described by three moments: n, <v>, and
<§>.

43

LIE e—g/ZkBT

 

f (i) = (3.6)
J?! (kBT) 3/2

where f (E) has been normalized so that

gm.) at. = 1 . (3.7)

Using the Maxwellian distribution function, the average energy of all particles in

thermodynamic equilibrium can be determined by (3.4) as

(i) = ngT . (3.8)

3.3.2 Boltzmann transport equation

When an external force is applied to the system of particles, the transport
behaviors of the particles inside the discharge can be characterized by the Boltzmann
transport equation (BTE). The Boltzmann transport equation, which describes the motion
of the distribution function due to an external force in coordinate and velocity space, can

be written as

3f . E. _(§£)
37+v Vf+m Vf— at can (3.9)

where V, is a gradient in coordinate (x, y, 2) space, Vv is a gradient in velocity space, F
represents external forces, and m is particle mass. The right hand side of equation (3.9)

includes the randomly-timed scattering events that the particles experience, such as elastic

44

collisions, ionization collisions, and recombination collisions.

Since equation (3.9) does not have a direct, closed form solution, an alternative
way to solve this equation is to make a series of simplifying assumptions to obtain the
moments of the Boltzmann transport equation. The basic assumption is that the
distributions of each plasma species can be described adequately by three macroscopic
quantities, namely, particle density n, drift velocity (v), and temperature T. The use of the
moment equations, or ﬂuid equations, considers the plasma to act as a ﬂuid, rather than as
individual particles. The condition for validity of this approximation is that the distance
between particles be small with respect to the interparticle forces, i.e., r2230 » 1 , where AD
is the Debye length, and the mean free path << scale of change of macroscopic quantities.
This condition is well satisﬁed in the plasma discharges investigated for this study. In the
following section, the ﬂuid models used in this study are constructed by obtaining the

moment equations for the Boltzmann transport equation.

3.3.3 Moment equations

The moment equations are obtained from equation (3.9) by multiplying it by
various functions of velocity, <1>(v) , and integrating over velocity space. The result is
changing from the BTE (one equation) which is a function of (x, y, z, vx, vy, vz, t) to a
series of equations which are just functions of (x, y, z, t). Various (I) (v) functions include
1, v, and vv, etc. which give rise to the zero-order, ﬁrst-order, and second-order, etc.
moment equations.

Hence, multiplying (3.9) by (I) (v) and integrating over velocity space, the general

moment equation may be written as

45

a ne ad) _ a
5Euler») +V0(n(v¢))—;E- (57> — [502mm] (3.10)

coll

where F is replaced by (:19, where e is the particle charge and E is the electric ﬁeld.

Thus the BTE for the distribution function f is replaced by a set of equations
containing averaged quantities. Each moment equation introduces the next higher-order
velocity moment due to the second term in equation (3.10). The moment equations are
then an inﬁnite set of equations unless some additional assumptions are used to break the
chain of equations and restrict the variables to a manageable number. These additional
assumptions will be considered below.

The zero-order moment equation is obtained from equation (3.10) by putting (I) =1

[0 get
an _ an
37+V.(n<v>) - [§]c011. (3.11)

The zero-order moment equation (3.11) is the particle continuity equation. According to
(3.11) the change of particle density plus the divergence of particles equals the change in
density due to collisions.

The velocity v of the particles inside discharges can be expressed as the sum of the
drift velocity and thermal velocity c, i.e., v = (v) + c, and in most case, 6 » (v). By
putting (D (v) = mv into equation (3.10) and assuming the particle pressure is given by
P = rim (cz)/ 2 = nkBT (thermodynamic equilibrium condition), the ﬁrst-order moment

equation can be derived to be

46

col!

6 1 e _ a
37(v)+(v)-V(v)+'-;-,-IV(nkBT)—;IE _ [576)] (3.12)

where the collision term can be written as

8 1 a (v) a
— = — — —— — .1
[at<v>]coll mu [3! (mn<v>)]coll n [8100] col! (3 3)
which represents the rate of drift velocity (momentum) change due to elastic and inelastic
collisions.
Similarly, with the same assumption for particle pressure, the second-order

moment equations can be obtained by substitute (I) (v) = émvv into equation (3.10),

3?, + V0 ( (v) W) = en(v) o E — V0 ((v)nkBT) - VOQ + [gr/lo“ (3.14)
where W is the average energy density (W = énm<v>2 + gnkBT), and Q is the heat ﬂux
(Q = —3.VT, 3. is the heat conductivity). Equation (3.14) contains on the left-hand side
the rate of change of the average energy density plus the outﬂow of average energy density
W. On the right-hand side, the ﬁrst term is the power supplied by the electric ﬁeld, the
second term is the work performed by the particle pressure, the third term is the
divergence of the heat ﬂow Q, and the last term is the rate of change of particle energy
density due to collisions.

The zero, ﬁrst, and second moment equations shown in (3.11), (3.12) and (3.14)
are also called the continuity equation, momentum conservation equation, and energy
conservation equation, respectively. These moment equations sometimes are simpliﬁed to

construct the drift-diffusion model as shown in the next section.

47

3.3.4 Drift-Diffusion model

The Drift-Diffusion model consists of solving the Poisson equation and the
continuity equation (3.11). This model makes a simpliﬁed approximation for the
momentum equation (3.12) and neglects the energy conservation equation (3.14). The
drift-diffusion model assumes that all the particles in the discharges are at thermal
equilibrium and the particle temperature gradient is zero.

Assuming the particle considered here is an electron, e = —q, the ﬁrst step to
simplify the moment conservation equation is using a simple relaxation model to replace

the collision term in (3.12), i.e.

[%(v)]cou = —'< v)vm (3'15)

where V", is deﬁned as the collision frequency for momentum transfer. This assumption is
valid because the major momentum loss of electrons in partially ionized discharges results
primarily from the elastic collision with neutrals. Moreover, assuming the collision term

dominates over the inertial term, the momentum conservation equation becomes

(v) = -——————. (3.16)

11,. = .1— (3.17)

and

48

kBT
”a = 7“":

then the expression for the electron velocity can written as

V
<v,> = -u,,E-D,.7"

and the electron ﬂux becomes

1,, = n(vn) = —nunE—DnVn

(3.18)

(3.19)

(3.20)

This representation of electron transport is referred to as the ‘Drift-Diffusion’

model and can be extended similarly for ions in the discharges. The equations in the Drift-

Diffusion model for electrons and ions are the electron continuity equation, ion continuity

equation, and Poisson equation, which is used to provide a self-consistent space charge

ﬁeld. The equations are written as

(3.21)

(3.22)

(3.23)

where ‘I' is the plasma potential, It, and n,- are the electron and ion density, and the electron

and ion current density 1,, and J,- can be expressed as

49

Je = —neueE—Devne (3.24)
and
J. = nipiE—DiVn‘. . (3.25)

1

3.3.5 Fluid model

The Drift-Diffusion model discussed above is valid only under the assumptions
that the discharge is in thermodynamic equilibrium condition, i.e. Tam,” = Tim
=anmls, and has a spatially uniform temperature distribution. However, this is not true
for microwave discharges. For the microwave discharges investigated in this study,
elecu’ons absorb microwave energy more efﬁciently and hence have a higher temperature
than other plasma species, such as ions and neutrals. Due to most of the microwave energy
being imparted into electron gas, and the electron temperature plays an important role in
determining the transport parameters and reaction rates. Therefore, the energy
conservation equation of electron has to be included in the ﬂuid model in order to provide
a self-consistent solution. The ﬂuid model used in this model then includes the ion
continuity equation (3.23), the ion current equation (3.25), the electron continuity
equation (3.22), the electron current equation (3.24) and a variation of the electron energy
conservation equation (3.14).

The energy conservation equation of electron used in the ﬂuid model for this study
is based on equation (3.14). By assuming that the energy distribution of electrons is

ngTe » %m(v)2 , (3.14) can be simpliﬁed to be [51]

50

() = wants-:4...

where qe is the electron heat ﬂux, which includes both heat convection and conduction

terms
5 5
q, = gksTJe- EkBDene VTe. (3.27)

Note that the thermal conduction coefﬁcient A in (3.14) is approximated to be gDene

[52].

For partially ionized discharges, where the neutral densities are much greater than
the electron and ion densities, the right~hand side of the continuity equations of electrons
and ions, (3.22) and (3.23), are primarily determined by the inelastic collision events of

electrons, such as ionization (electron-neutral interaction) and recombination (electron-ion

interaction). These can be written as

call

a a . . . .
1;] [.41 3.2.3.7.).
(:0 j j

where Rio" is the ionization rate (sec'l), or, is the electron-ion recombination coefﬁcient
(m3 sec'l), and index j stands for different ion and neutral species. Equation (3.28) also
implies that the electron attachment mechanism is not signiﬁcant in the discharge of
interest. This is true for the argon gases portion of this study.

The electron energy loss due to collisions can also be more completely speciﬁed as

the energy loss due to several electron-neutral and electron-ion interactions, including

51

electron-neutral elastic collisions, ionizations, excitations, dissociations, and electron-ion
recombinations. Thus the collision term on the right-hand side of equation (3.26) can be

written as

[314,” = «(24.2... gem.” 122.22...)

(3.29)

2me 3 3 . . 3 . .
- 27( iksTe ‘- §ks7i)"’m"e - iksTeZdr'éne
J J

n

where Rio", Rm, and Rdis together with Siam em, and 3dr: represent the rate and threshold
energy for ionization, excitation, and dissociation. vm is the electron-neutral momentum
transfer frequency. Moreover, the collision rate R for neutral ionization, excitation and
dissociation can be expressed as

R nk

ion, ext, dis = n ion, m, dis
where n" is the neutral density, and k is the rate coefﬁcient with unit, m3/sec, which is
usually a function of electron temperature [50].

In summary, the ﬂuid model for plasma discharges involving both electrons and
ions solves self-consistently the Poison equation, the continuity equation for electron, the

continuity equation for ions, and the energy conservation equation for electrons. Since in

this study we are interested in steady state phenomenon, these equations can be written as

V”? 1:01, — 72,.) (3.30)

V.Je Znennlgiorrzalrnline (3°31)
1° j

52

V011.

2’12”" 1' kino—Zdnlne

V.qe = -—qJe OE

—n n "(2” [0,, ,0" +23, ext “#253.“in
—;—nze(§khT 4k7“)v'me—n ~23kTZajn121e

where

k
a
II

—neueE—DeVne

and

=5kBTJe— (SkBDn )VTe.

3.4 Numerical solution implementation

(3.32)

(3.33)

(3.34)

(3.35)

(3.36)

In this section, the implementation of numerical solutions for the models discussed

in the last section will be presented. To solve the ﬂuid plasma model numerically involves

discretization of the equations at a set of grid points. Setting up the discretized equations

for the grid points constructs a nonlinear matrix problem which typically needs to be

solved iteratively. If the iterations result in a convergent solution, the unknowns on the

grid points are solved. Several numerical methods have been used in the literature to solve

53

the ﬂuid equations including the ﬁnite difference and ﬁnite element methods. For this
study, the ﬁnite difference method has been used to solve these equations. The details of

the numerical solution are presented below.

3.4.1 Discretization

The numerical solution of the ﬂuid equations requires the discretization of the
equations (3.30) to (3.33) for an appropriate grid structure superimposed on the simulated
discharge geometry. For cylindrical microwave plasma sources, the discretization of the
equations should be done in cylindrical coordinate. Since the resonant electromagnetic
modes of the microwave plasma sources used for diamond thin ﬁlm deposition are
normally TMorn mode, the electric ﬁeld spatial distributions are (1) symmetric in nature and
the 13¢, component is zero. Therefore the discharge behavior can be assumed to be c)
symmetric and the discretization of the equations reduced to a two dimensional problem.
Thus, the simulation region remains in the r-z plane only.

In order to obtain a symmetric and stable solution from the simulation, a gird
structure arrangement is developed according to the staggered-mesh arrangement [43][59]
as shown in Figure 3—2. The area within the dotted square represents a simulation cell, (i,
k). The indices along the radial and axial axes are represented by i and k, respectively. The
positions of plasma densities, potential and electron temperature are deﬁned at the center
of each cell which is represented by the dark circle. The r and 2 components of the current
density and electric ﬁeld are deﬁned at the positions marked with an empty circle and gray
circle, respectively. The quantity values at undeﬁned positions can be interpolated from

their value-deﬁned neighbors. For example, the electron density at the boundary between

54

(i-l, k+1) (i, k+1) (i+1, k+l)

O O O O O

   

(i+l, k)

 

 

 

Figure 3-2: Staggered-mesh arrangement.

55

cell (i, k) and cell (i, k+1) can be determined by averaging the electron density of positions
(i, k) and (i, k+l).

The discretization of the Poisson equation is done using a ‘ﬁve-point’ ﬁnite
difference approximation derived from truncated Taylor series. The details of the

derivation are shown in Appendix. It gives in cylindrical coordinates

 

 

 

- 2 ‘ ri+l/2 .

2 2 - Ar. -‘I’(r+1,k)
’i+1/2"i—1/2 ‘

. 2 ‘ 0-1/2 .
+ 2 2 .Ar- -‘l’(t—1,k)

’t+1/2”t-r/2 "1

+ 2 .‘I’(i,k+1)

(Azk+Azk_l) Azk

2 .trl (i, k- 1) (3.37)

 

I (Azk+Azk_l) Azk_1

__ 1 2 -. ’z+1/2+’z-1/2
2 2 Ar. A,_
’i+1/2”z-1/2 ‘ "1

2 1 1 .
+ (Azk+Azk_1) (2T2); ”Lawn W”)

.. _%(n‘.(i, k) —n,(i, k))

 

 

 

 

For electron continuity equation, (3.31) can be written in discretized form as

2 . l
‘2 2 "’i+1/2'Jer(‘+§’k)
('t+1/2”i-1/2)
: 2 MI 1 l
2 2 1/2 "
L j ' 2 (3.38)

’i+r/2 ’i— 1/2

(Azk+2Azk_ 1) 'Je( ("k+2)'1¢2("k'2))

= ne (1', k) [nn (1', k) kion (i, k) —a, (i, k) ni (i, k) ]

 

56

The electron ﬂux in (3.34) can also be discretized as

 

 

 

 

. l - 0.. §W(i+1,k)—W(t,k) . .
Jer(l+2,k) — [—Z—ri+ 2( Ari )] ne(t+l,k)
(3.39)
+ &+“_¢(‘¥(i+1,k)—‘l’(i,k)) .n (ik)
Ari 2 Ari ‘3 ’
and
D . .
1,2(i,k+l)=[_—’.+”—’(‘P("k“)“W(""))]-ne(i+1,k)
2 2" 2 Azk (340)
De lie ‘P(i,k+1)—‘i’(i,k) . '
+[—k+-2—( Azk ) -ne(z,k).

The ion continuity equation is similar to that of the electron continuity equation
and is simply the substitute of index i for e and lie for -ui in (3.38), (3.39), and (3.40).

Finally, the electron energy conservation equation can be discretized using the
same methods discussed above. Note that the heating term we 0 E in equation (3.33)
primarily represents the microwave power absorbed by the plasma, which is determined
by the microwave electric ﬁeld and the microwave induced current density using the
electromagnetic ﬁeld model developed in Chapter 2. Therefore, for simplicity, this heating
term can be viewed as an input parameter, Pabs ,for the energy conservation equation, and

it is not necessary that it be discretized. Equation (3.33) can be discretized as:

57
4 2

A . 1k
2 2 "rust/2’4” “’2’
(ri+1/2_ri-1/2)
, 2

‘ . 1
2 2 ' ri— 1/2 ' qer(l_§’ k)
l' )

t+1/2"t—1/2

+ 2 -( (w)- (+1))
(Azk+Azk_l) 482’ 2 qez’ 2

= ne(i, IozH’u, k) +Pabs(i, k)
,-

(3.41)

 

where 211’ (i, k) represents all the collisional energy loss in (3.33). The heat ﬂux terms

1
are discretized as

 

. 1 _5 (Te(i+1,k)+Te(i,/()) ( 1 )
qer(l+2’k) '- EkB 2 J” l+'2',k

 

 

 

 

 

5 D (ne(i+1,k) +ne(i,k)) Te(i+ 1,k) —Te(i’ k) (3.42)
—§ 8 e 2 Ari
and
. 1 __5 (Te(i,k+1)+Te(i,k)) (.1)
qez(l,k+§) - ikg 2 Jez t,k+i 4
5k D (ne(i,k+1) +ne(i,k)) Tam/(+1) -Te(i,k) (3. 3)
_§ 8 e 2 Azk

Equation (3.37) to (3.43) are derived for the case that r is not equal to zero, i.e. not
at the center of the cylinder. The discretized equations for r equal to zero and their

derivations are also given in Appendix.

58
3.4.2 Numerical solution
The solution of the discretized ﬂuid plasma equations on a grid of N nodes
involves a total of 4N unknowns for the ﬂuid model since at each point the unknowns are
‘l’, ne, ni, and Te. The system of 4N dicretized equations including the Poisson equation,
the continuity equation for electrons and ions, and the energy conservation equation for

electrons are

”F“, (‘1’, ne, ”i3 T) ‘

F (‘P T ) Fn‘ (W, "e, n,, Te) 0 (3 44)
, n , n., = = ,

e ‘ e F"; (‘1’, ne, n,» T)

fr, (\P, ne, ni, Te) .

 

 

where F t}. denotes the system of the N poisson equations, F M and F n) denote the system of
2N continuity equations for electrons and ions, and FTC. denotes the system of N energy
conservation equations for electrons. Two principle methods of solving the system of non-
linear equations (3.44) are the fully coupled method or Newton method and the decoupled
method or iterative method [60] [61].

The Newton’s method linearizes the system of partial differential equations and,
by starting from an initial guess, the solution of the nonlinear equations is obtained by

iterating the matrix equation

59

 

 

 

 

r “k
317., an, an, an,
3‘17 3: Tn: 57¢ , -
- - F ‘l’k k nk 71‘
3F" 8F" 8F" 3F" N? ‘P ’ e’ i’ e
_._. _v_' k k k
a? ane ani aTe Ant F": ‘P’ e’ "PTI; (3 45)
, — _ k Ir. k °

3F". 8F": 8F". 3F": A", p4? HI.» "1., 71f)
W T T 57‘— AT

me n. e .. e. FT(‘Pk’ "I; "f, 71;)
arr. aFT‘ 3F... BFT‘ - ‘
as! arts Bn‘. 8T3

 

 

where k denotes the iteration count. The correction vector for the k-th iteration is given by

A‘Pk = ‘Pk+ l _ \Pk
An: = nelH-l—n k
Anik = nik + l _ n‘k (3'46)
ATek = Tek + l _ Tek

The iterative method decouples the equations such that each system of equations
can be treated independently for each iteration cycle. Each system of equations are solved
sequentially and variables are updated after each solution. When all the updates are less
than a set of criterion for an iteration, the convergence is reached. This iterative method

can be written as

 

[ant 1:, ..g, :79)”. ",4 r ,. .7.) (347)

8‘1, 9 "C, ni’

 

3F ‘Pk+1,nk’nk,Tk
[ 111‘}, an "e "t e) Ank=—Fn(‘l’ k+1 ’ , "UT/2) (3.48)

e

 

an.

I

k+1 nk+l k
aF (W 9 , "’7‘:
[ n‘. nen1e)] An.k = _Fn.(‘Pk+l,n/C+l,nﬁ 7*) (349)

 

8F ‘Pk+1n k+l nk+1Tk
[ 13(9) ATek=—FT(‘l’an k+1 "@1713 (3.50)

87’

C

where (3.47) to (3.50) are solved sequentially. At each stage only one system of N
equations is being linearized and solved by the Newton method, so the matrix being
solved has N rows by N columns regardless of the number of coupled equations being
solved. Whereas, in the full Newton method, all of the coupling between variables are
taken into account and the matrix size will be 4N x 4N.

Because of the tight—coupling, the full Newton’s method requires a smaller number
of iterations for convergence, however, since the matrix size is larger, it takes more time to
solve each iteration. Moreover, the Newton method needs larger computer memory to
store the elements in the matrix. The iterative method usually takes considerably more
iterations to reach the convergence while the time for each iterations is often less than the
full Newton method because of the smaller matrix size. However, due to loosely coupled
feature, the iterative method is less stable than the full Newton method, especially for high
density plasma simulation.

In this study, both Newton’s method and iterative method are used to solve the

discretized equations. The Newton’s method is used to solve Poisson’s equation, electron

61

continuity equation, and ion continuity equations. The matrix size is 3N x 3N and the
plasma potential, electron density and ion density are solved in a tightly-coupled manner.
These solutions are iteratively coupled with the discretized electron energy balance
equation to solve the electron temperature and this temperature is then put back to the
continuity equation and Poisson equation to modify the solution of densities and potential.
The final solution is achieved by iteratively solving the continuity/Poisson equations and
electron energy equation.

Since the electron energy balance equation and charge continuity equations are not
tightly coupled, a damping technique has to be used in order to prevent unstable numerical
problems such as oscillation or divergence from occurring. One damping method applied
here is to use only a portion of the updated data as the input to the next coupled equation.

This method can be written as

n+1

Q = Q" + LQ" +]— Q") x damping factor (351)

where Q represents a quantity to be solved (such as ne or Te) and n is the time index. (3.51)
indicates that the new quantity which is used as the input for next coupled equation in the
iterative solution is the sum of the previous quantity and a portion of the difference
between the previous quantity and updated quantity. The damping factor is usually less
than one to secure solution stability. If the damping factor is equal to one, that means no
damping is applied for the solution. The value of the damping factor depends on the
current situation of the solution. If damping factor is too high, the solution might diverge
very fast, however, if too low, it will take more loops to get a converged solution. A

automatic damping factor algorithm was implemented by checking the value of Ftp, Fne,

62

Pm and FTe deﬁned in (3.24). When an exact solution is reached, the value of Fry. Fne, Fm
and FTe should equal or be close to zero. If the current values of Fry, Fm, Pm and FTe are
higher than the previous ones, that means the solutions are starting to diverge and a
smaller value should be used to secure stability. Otherwise, the damping factor can
become larger to speedup the convergence. This damping technique is also used for the
Newton method while calculating the potential and densities to ensure the stability of the
solution.

Both the Newton method and iterative method produce a linear matrix equation of
the form Ax = B. This matrix equation can be solved by direct or iterative techniques [55].
In this study, this linear matrix equation is solved by the LU-decomposition method. Since
matrix A is a banded, sparse matrix, after a careful arrangement of matrix A, the size of A
can be reduced from 3N x 3N to (6n,+5)x3N, where nr is the number of grids in the r
direction. The reduced matrix A can save substantial memory space and increase the speed

of calculation [60].

3.5 Summary

The ﬂuid plasma model of this study was developed in this chapter. The model
consists of the Poisson equation, electron and ion continuity equations, and electron
energy equation. The numerical methods used to solve these equations were also
discussed. This ﬂuid model can be used to investigate the characteristics of discharges and
it can be coupled with the FDTD model to self-consistently simulate the electromagnetic
excitation of discharges inside a microwave cavity plasma reactor. These simulation

results will be present in Chapter 6.

Chapter 4

Transient State Simulation of Electromagnetic Fields
Inside Empty and Loaded Microwave Resonant Cavities

4.1 Introduction

In Chapter 2, we discussed the fundamentals of the FDTD method used for solving
Maxwell’s equations in order to get the electromagnetic ﬁeld solutions inside cylindrical
cavities. Before simulating the complete microwave cavity plasma reactors operating in
steady state with complex plasma loads, a simple-structure cylindrical resonant cavity will
be simulated. Transient simulations of empty and loaded cylindrical resonant cavities will
be done in order to establish the accuracy of the FDTD model for microwave cavity
plasma reactors.

Several simulation results based on the FDTD model are provided in this chapter
to study the electromagnetic mode variation and energy transfer inside cylindrical
microwave resonant cavities. The electromagnetic resonant mode solutions of cylindrical,
lossless, empty cavities, operating in either TE or TM modes, are determined in the time
domain and are compared with theoretical results to verify the accuracy of the FDTD
model. Two types of lossy loads in the resonant cavity are considered here. The ﬁrst type

is a non-dispersive load where the current density needed to solve Maxwell’s equations is

63

(,4

determined from the conductivity using J = (SE. The second type is a dispersive load
where the current density is determined by simultaneously solving the momentum
transport equation in the time domain. For the non-dispersive load type, FDTD solutions
were directly compared and found to agree with known analytical solutions. For the
dispersive load type, FDTD solutions were compared to experimental measured argon
discharge data. The simulations model the interactions between cavity characteristics,
such as natural frequency and quality factor, and the properties of the plasma in
microwave resonating cavities. Finally, the FDTD simulation results of a compact ECR

ion source loaded with a helium discharge will be presented.

4.2 Simulation of electromagnetic ﬁelds inside empty cavities

In order to demonstrate the application and validity of the above FDTD
formulation in the microwave plasma cavity, empty resonant cavity solutions are ﬁrst
described. Consider a cylindrical cavity geometry with perfectly conducting walls as
shown in Figure 4-1. The boundary condition for the electric ﬁelds in such a cavity is that
only the normal component of electric ﬁelds exists on the perfectly conducting wall and
the tangential electric ﬁelds at the wall are zero, i.e., E2 = 0 and E4, = 0 when r = a, and E,
andE¢=0whenz=Oand L.

The way chosen to numerically incorporate an electromagnetic ﬁeld excitation
source inside the cavity for these empty cavity simulation is based on the speciﬁc cavity
mode of resonance. The resonant mode depends on the diameter and height of the
cylindrical cavity. The natural frequency (or so called eigenfrequencies) of a lossless,

empty cavity are given by:

65

APerfectly conducting walls©

xcitation circle

 

 

 

 

 

 

 

 

 

 

Excitation line
a ’ (Er)
TEou mode TElll mode
Excitation surface
Er Er
(+) (+)
z
TMon mode TM012 mode

Figure 4-1: TEOH, 'I'Em, TMOH' and TM012 mode excitation techniques in a
cylindrical cavity.

66

 

 

(it: = shelf—.3112?
’ 2
2.3.2112») (2.3-‘1’

where TE and TM are the transverse electric mode and transverse magnetic mode

(4.1)

 

 

respectively, npq are integers for TB and TM mode speciﬁcation, A!!!) is the p‘h root of the
Bessel function Jn(x), and Amp, is the plh root of the ﬁrst derivative of the Bessel function
1,100 (J ’n(X)).

One way to implement a source is selecting several lattice points as source points
and assigning magnitudes of an electric ﬁeld component at these points based on
theoretical cavity ﬁeld solutions [24]. These source points are driven a few cycles at a
frequency close to the resonant frequency, then turned off. The electromagnetic ﬁelds then
resonant inside the cavity at a frequency which depends on the cavity size as given by
equation (4.1). This technique of exciting a mode inside the cavity and then turning off the
source gives the natural frequency response.

For example, consider the excitation of the TECH mode. Assume a cavity is 15.25
cm in diameter and 6.5 cm height, the resonant frequency for the TED“ mode can be
calculated using equation (4.1) as 3.33 GHz. The source points are chosen along a circle in
the middle of the cavity height as shown in Figure 4-1, where the Eq, component has its
maximum amplitude based on the theoretical resonant cavity ﬁeld distribution. A uniform
amplitude of E¢ is assigned along this circle and forced to oscillate at this frequency for a
few cycles. The electromagnetic ﬁelds propagate from these points and resonant inside the

cavity. The electromagnetic ﬁeld energy U stored in the cavity is deﬁned as

67

U(t) = 0.51 (8(r) (E(r,t))2+p.(r) (H(r,t))2)dV (4.2)
v

where the integration is performed over the volume of the cavity, V. The electromagnetic
ﬁeld energy starts to increase as shown in Figure 4-2 until the source is turned off (0.5
nanosecond), at which time the stored energy remains constant because no energy loss
occurs on the perfectly conducting walls. From the plots of the electric ﬁeld variation
versus time, the oscillation frequency after the source is turned off at 0.5 nsec can be
obtained as the natural resonant frequency (3.33 GHz) as shown in Figure 4-2. The
electric ﬁeld has E¢ components only, and its spatial variation in the r-z plane matches the
theoretical empty cavity results as shown in Figure 4-3. The theoretical empty cavity
mode solution can be obtained by solving the time-harmonic Maxwell’s equation in
cylindrical coordinates as mentioned in Chapter 2. For TED“ mode, the electric ﬁeld

solution is given by

I

A01 r r r - N
E = E¢ = 8710(A015)Sln(z) (4.3)

where 1’01 is the ﬁrst zero of the ﬁrst derivative of the Bessel function JO, and B is the
amplitude.

The number of grids used is 31 by 36 by 31 in the r, q), and z directions,
respectively. The grid spacing in each direction is uniform. The grid structure is shown in
Figure 4-4. The unit time step is chosen as 0.5 psec to satisfy the ﬁnite difference stability
condition as discussed in Section 2.3.

Other examples of mode excitation are the TM), 1, W012 and TE”, modes. For

68

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

120
m 111(1111111111111
26A 80-
8.?
:5
5V4o- -
20~
0 UeSUourUseUflfU U UU U U U U U UU UU
o. 0.5 1 1.5 2 2.5 3

Time (ns)

2
Figure 4-2: Electric ﬁeld energy, "VUEOE dV, versus time for the TED” mode. The
V

empty cylindrical cavity was 15.25 cm in diameter and 6.5 cm high.

69

 

50 I I I F I I T

45- 0006.6"- 0 Simulation -
z=3.03 cmO

40- ,0 O ..Theory ‘

35. C's-s. -

.0
. z=1.95 cm 0 “O O"
30- O- 0.. o

J. . ‘0..O-d.o.' . ..
25" O.- 0"0 o 0'0... 0'2. 1
0.0“ 2:52 cm 0", e, .

.o. .

IE¢| (Arbitrary Units)
0

o " z=0.87 cm '02 P,

10- 8
@

 

 

 

Figure 4-3: Radial dependence of IE¢| at four different 2 locations for the TED“ mode.

70

(r - (p) Plane:

 

(r - z) Plane:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4-4: Grid structure for FDTD simulation in cylindrical coordinate.

71

TMm,n modes, instead of using a circle as a source, a cylindrical surface can be used as an
excitation source as shown in Figure 4-1. This surface has a height equal to the cavity
height, and a radius where E, has a maximum for these modes. E, is varied sinusoidally
along the excitation surface based on theoretical resonant solutions. Using the same
procedures mentioned above, the simulation results of the resonant electric ﬁeld variations
with time and space for the TMOH mode are shown in Figure 4-5 to Figure 4-7. In
particular, the total stored energy as shown in Figure 4-6 is constant versus time after the
excitation source is turned off at 1.25 nsec. Figure 4-8 shows the electrical ﬁeld variation
found using the FDTD technique and ideal empty cavity solutions versus r and z for the
TM012 mode. Figure 4-9 shows the simulated E, variation versus (1) for the TEm mode.
The TMOH, TMmz, and TE,“ modes solved by the FDTD method are all well matched
with the theoretical solutions.

Due to the perfectly conducting wall assumptions, there is no energy loss for the
empty cavity model, and the total energy inside the cavity remains constant once the
source is turned off as shown earlier in Figure 4-6. Steady state solutions which are
obtained by the source being always in an “on”‘condition can not be obtained for these
ideal empty cavity simulations because no loss exists. The simulation technique which
gives the steady-state response (the source is always on) for loaded cavities will be shown

in the next chapter.

4.3 Natural frequency response for loaded lossy cavity

One application of this FDTD model is to investigate the power absorption inside a

cavity when a lossy non-dispersive load is present, i.e., the dielectric property of the load

72

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

140 r r r r r
111111111“
>, 100’ i
is” a
,5 E
D 80 -
E o
.93. 2.
Li" ‘5
.0 "a"; 60- -
51 es
.2
m 40- -
20- -
«1011141 U U 1 U U U.
0 l l l l l
0 0.5 1 1.5 2 2.5 3

Time (ns)

Figure 4-5: Electric ﬁeld energy, ‘1’] eoEde, versus time for the TMOH mode. The
V

cylindrical empty cavity was 17.8 cm in diameter and 7.2 cm high. The
resonant frequency was 2.45 GI-lz.

73

 

 

 

1.2 l I I I l
1)-
it;
‘3 8
LL}
1: D
2 E‘
9 E! 0.6-
m 23
g E
e— V 0.4»
0.2-
<-Source off
00 0.5 1 1.5 2 2.5
Time (ns)

Figure 4-6:

Total stored energy, U, versus time for the TM01, mode.

 

74

 

_ 9. .o -o .0. 0 Simulation

0 o
.0" r= 7.36 cm “Q
Q. 'o. Theory

30*

O" .I. ........ _-... 'Q .
25L .-"o O O O O o ‘
.0. ~30 O..'. O.
. "0 r= 4.29 cm 0".
20- -‘ 33' b. ‘

”6.6.6.6..“
..-'o O O o‘-..

O ,' .‘O O". .
15- ,0 ..o r: 2.76 cm 0... on. ‘

IErl (Arbitrary Units)
0

.Q-O'O'O‘O'o.. .
.00 0.0.0. 0.

,g'p' 0,.0’ r= 1.23 cm "0-
o'
o

 

 

 

..L
N
ml-
.b
01
O)
‘1
m

2 (cm)

Figure 4-7: IE,| (arbitrary units) versus 2 at four different r locations for the TMOH
mode.

IE,I (Arbitrary Units)

75

 

35

30*

20

15

10

25~ .5'0 0‘2. Q. .95 =55

. o. .. , .. .
r ,- b r= 4.29 cm Q -
~ :0 o. _o

’ '- "o e -'

I T I I I I I

0 Simulation
Theory

3'0 O, '. .‘ :0 .0 Oh. .. ..

of r: 2.76 cm '5 ('21-. {F5 ‘5' "2&9

O on. ..i O
, .. 0.00., - .- .o-o-o, -.
o .0 6 - G 202‘. 59: o- ' 0 0-. 0_

r: 1.23 cm "Ga. 63;. 5‘9 .0" "o. '

 

 

 

Figure 4—8:

16

lErl (arbitrary units) versus 2 at four different r locations for the TM012
mode. The cylindrical empty cavity was 17.8 cm in diameter and 14.41 cm

high. The resonant frequency was 2.45 GHz.

76

 

100” I I I .0. ﬁ I I
530.61??? 0 Simulation
gym. Q": .0. O. O, Th 0': ED
' ' 0’ ' o eo . o'
8 _ O-._o-. r=2.15 gm ry A
. . . 0' .0. .2 of :0

a, o .0' 6 '25.

lErl (Arbitrary Units)
‘1‘
O»
O

o :9 '53; o
o 5’ e

—L
I
Ila-5"“...
.‘
H

 

 

l l l .6: l
0 so 100 150 200 250 300 350

 

¢ (degree)

Figure 4—9: ”5.4 (arbitrary units) versus (1) at three different r locations for the TB“,

mode. The cylindrical empty cavity was 8.9 cm in diameter and 6.9 cm
high. The resonant frequency was 2.45 GHz.

77

is independent of the 'esonant frequency. The natural frequency describes the transient
response of the cavity to a pulse or to having a CW source shut off. For a cavity with a
lossy load, the natural frequency is no longer a real number, but a complex number which

can be written as
6) = (0’ +jw” (4.4)

where the real part (0’ of the natural frequency 6) describes the frequency of field
oscillation while the imaginary part 0)” of the natural frequency 6) is the exponential
decay coefﬁcient due to power absorption.

The FDTD model described in Chapter 2 can be used to analyze the natural
frequencies of a lossy loaded cavity, by assigning a set of parameters including a, u and 0,
at each lattice point inside the cavity. For a natural frequency response simulation, after
the source is turned off, the electromagnetic ﬁelds inside the cavity will decay at a natural
oscillating frequency, and the stored ﬁeld energy will decrease due to energy dissipated
into the lossy material. The relations can be expressed as follows:

If ' I
e-O) tgwt

E(t) E

O 2 ,, (4.5)
um = er' 0’ ’

where E0 is the initial electric ﬁeld amplitude and U0 is the initial stored energy at the time
the source is turned off. By studying the electric ﬁeld oscillation frequency after a source
is turned off and the rate of the stored energy decay, the real part and imaginary part of the

natural frequency 6) can be obtained.

78

As a speciﬁc example, consider a 15.24 cm diameter cavity loaded with a lossy,
homogeneous, non-dispersive load. The cavity is the same size as the TED” mode empty
cavity demonstrated in the previous section. The lossy load is 2.54 cm in diameter and it is
loaded coaxially with the same height as the cavity, as shown in Figure 4-10. To study the
effect of the load conductivity on the electromagnetic ﬁelds inside the load region, the
permittivity e in equations (2.15) to (2.20) is set equal to 3.080, and o is varied from zero
to several hundred. Near the boundary of lossy dielectrics, where the most power is
absorbed, ﬁne grids are used to improve the accuracy of results. This ﬁne grid structure is
shown in Figure 4-11. The plots of electric ﬁeld strength and total stored energy decay are
shown in Figure 4-12 and Figure 4-13. It is clear from these plots that the decay rate and
oscillation frequency varies with different conductivities. For example, at o = 0.1 mho/m
the result was 0)’ = 20.68 x 109 rad/sec and (1)” = 0.09 x 109 rad/sec, at o = 1.0 mho/m
the result was (0' = 21.13 x 109 rad/sec and 0)” = 0.45 x 109 rad/sec, and at o = 10 mho/m
the result was (0’ = 21.59 x 109 rad/sec and to" = 0.16 x 109 rad/sec. The radial variation
of 15¢, for different conductivity loads is shown in Figure 4-14. For the case of the larger
conductivity the electric ﬁeld decays rapidly into the 2.5 cm diameter load.

A s-plane frequency chart can be used to further demonstrate the effect of
conductivity variations on the electromagnetic ﬁelds. As shown in Figure 4-15, the
horizontal axis represents the imaginary part of the conjugate of the complex natural
frequency, i.e. —a)”, and the vertical axis is the real part of the conjugate of the complex
natural frequency, i.e. (0’. For 6 = O, the cavity is lossless and hence co” = 0. When 6
increases from zero, 0)” starts to increase, which is due to power dissipated into the lossy

material. After 6 is larger than 1.0, (1)” starts to decrease, because the cavity load

79

TED“ mode

B=15.24 cm
A=2.54 cm

L=6.5 cm

 

 

 

‘_ lossy load

 

 

 

 

 

 

 

Figure 4-10: Conﬁguration of a coaxially loaded cylindrical cavity.

80

(r - (p) Plane:

   

Boundary

/

(r - z) Plane:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

r=0

 

Lossy medium 4

 

Figure 4-11: Fine grid structure for TEOH mode lossy cavity simulation.

81

 

 

 

 

 

 

 

 

 

 

 

 

100 r t .
0.1 (mho/m)
N l- -1
E1 50
8 5 1 1.5 2 2 5
100 t 1 ,
1.0 (mho/m)
N .
g 50* a
8 5 1 1.5 2 2 5
100 t w t
10.0 (mho/m)
N a
Q 50
8.5 1 1.5 2 2.5
Time (ns)

Figure 4-12: Space averaged lEl2 (arbitrary units), versus time, for the TED“ mode.

Results are shown for three different load conductivities after the source is
turned off.

82

 

      
 

 

 

 

 

6 r t .
-—0.1(mho/m)
5- , ~---~1.0(mho/m)-
:-~\
1’ 'i:- -------- 10.0 (mho/m)
:3 4_ I \A‘.
g g ’ \\.-\
LL} :1 l \.
'20:! 2" l \\.\
3- \ -
9 E ii ‘\_.\
U) z: \_
g E l I“.
[— V 2‘ IX ‘
5 -‘\-~
‘3‘." \"\
1- “"\-~\ -
0 in Source offl . .
0 0.5 1 1.5 2 2.5

Time (ns)

Figure 4-13: Total stored energy, U, versus time for the T1301] mode for three different

load conductivities.

83

 

 

 

 

200- -—0.3(mho/m) r
-- - - 1.0 (mho/m) '
----- 20.0 (mho/m)
150 ~ a
”a“ "
'E
D
E: ..
..E 100'" .z' ...... \.\ d
'8 z". \ . \.
s /'/. \.\
_e .'.,- \\
m ’ ’\
_ , .
50 ' .x' _. \.\ "
’0/ .\
I' \.
/ \.
1' ; \. 2
/' .' \ '.
0 I ---- .1. . . <-- Boyndary l l l l l .
O 1 2 3 4 5 6 7 8
r (cm)

Figure 4-14 Radial dependence of IE¢| (arbitrary units) for the TEo“ mode for three
different conductivity loads in a 15.25 cm i.d. cylindrical cavity.

84

 

 

 

 

21.8 . 1 ' '
o . XXX
X FDTD Simulation result x
>.‘...-"".(5=200.O
. ~ ....A l ' I ‘
21 6 na yttc resut >.< 0:200
A . o=5.0
§ 21.4” x q
a: ' =
E . o 2.0
°é> 21.2~ -
... x
3: 5. 0:1.0
* ..
g 21~ "- i
0
a: o=0.5
‘x-..
o=0.3
20.8r ......... x .......... 0:01 6:0 01
.................. x......x.....)()(
20. ‘ ‘ ' '
-85 -0 4 -0.3 -0.2 -0.1 0

Im((o*) x 109 (rad/sec)

Figure 4-15: S-plane plot for the complex frequency as a function of load conductivity
for the TED“ mode.

85

becomes a conducting-like material and the electromagnetic ﬁelds only penetrate the load
a small amount (i.e. a skin depth). Thus, the resonant frequency changes (increases) as the
TED“ empty cavity mode (0 = 0) varies to a coaxial type mode when 0 is large.

For the case of a cavity loaded by a lossy, non-dispersive coaxial load, analytical
solutions of (0’ and to” exist. The FDTD technique solutions are compared to an
analytical solution in Figure 4-15. In particular the FDTD result is compared with
Maming’s data [8], which is solved by an analytical method in the frequency domain. In
the analytical method, the electric and magnetic ﬁelds are solved in cylindrical coordinates
in terms of a set of Bessel function type solutions. The boundary conditions are applied to
form a characteristic equation which is used to determine (0’ and to” of the various
modes. This method provides analytical solutions for non-dispersive lossy loads. The
analytical solution and FDTD solution show agreement. However, if a more complex
plasma with dispersive, inhomogeneous properties is used, the problem becomes very
difﬁcult to solve using analytical solutions. Moreover, if the geometry of the load is more
complicated, it is also hard to solve for an analytical solution. Hence, the analytical
solutions are only used in this work to verify the FDTD solution on simple geometry, non-

dispersive loads.

4.4 Simulation results of an argon plasma loaded cavity

A FDTD calculation has been done to compare with the experiment results of
Rogers [6]. In Rogers’ study, a 17.78 cm i.d., microwave cavity excited at 2.45 (3112 was
coaxially loaded by an argon plasma in a 12 mm i.d. quartz tube with the plasma

extending from the top to the bottom of the cavity as shown in Figure 4-16. The radial

86

 

Cavity walls

Quartz tube

 

. rgon discharge

 

 

 

 

Figure 4-16: The geometry of a 17.8 cm i.d. cylindrical cavity with a 12 mm i.d. quartz
tube loaded with an argon plasma.

87

component of electric ﬁeld strength was measured by micro-coax electrical probes
inserted through the cavity wall versus 2 and 0. The measured ﬁeld distribution was
compared with theoretical results to verify the mode sustained in the cavity with the
discharge present was approximately a TM012 mode. The discharges inside the quartz tube
had diameters of 2-4 mm depending on the pressure which was varied from 40 to 200 Torr.
Rogers measured the cavity resonant length, average discharge diameter, total power
absorption, cavity quality factor and calibrated electrical ﬁeld strength at a reference point
versus variations in pressure. This measured data was then used to extract the electron
density and effective collision frequency versus pressure using (2.44) as shown in Figure
4-17.

Treating these values of electron density, effective collision frequency, and
discharge diameter as input data, the cavity Q and the power density of the discharge can
be simulated using the FDTD model. The deﬁnition of the quality factor, Q, is the time
average energy stored in the cavity divided by time average energy dissipated in the cavity
per cycle. Therefore, by calculating the time average energy stored and the time average
power dissipated in the cavity for the ﬁrst cycle after the source is turned off, a quality
factor can be determined. In the simulation of the discharge loaded cavity, the electron
density and the effective collision frequency are assumed uniform across the discharge
volume. Additionally, the discharge volume is evaluated as a cylindrical rod with the same
height as the cavity. The quartz tube was assumed lossless with a relative permittivity of e,
= 3.78 and a tube wall thickness equal to 1.0 mm. The electromagnetic source excitation
used is the same as the TM012 mode empty cavity shown earlier in Figure 4-1. Current

densities are calculated in time domain by equations (2.53) and (2.54) and coupled with a

88

13

 

 

 

 

 

 

 

 

3 x 10 r 1 1 I
2.5 - -
A O
c? o
E 2 - 0 -
8
c” 1.5 - O .
1 O J l l l
o 50 100 150 200 250
10
1.5 x 10 I I I I
l8 1 P O _
o
>
O
0.5 l I J L
0 50 100 150 200 250

Pressure (Torr)

Figure 4-17: Electron density, ne, and effective collision frequency, vejf, versus pressure
for the argon plasma shown in Figure 4-16 (from Rogers’ data [6]).

89

FDTD solution of Maxwell’s equations to solve the electromagnetic ﬁelds inside and
outside the discharge volume. Then, the power absorbed by the discharge is determined by
equation (2.35), and the absorbed power density is evaluated by dividing the power
absorbed in the discharge by the discharge volume. The simulation data of loaded cavity Q
and absorbed power density is compared to experimental data in Figure 4- 18 and Figure 4-
19. The FDTD calculations show good agreement with the experimental data, which
demonstrates the applicability of the FDTD technique to model plasma loaded microwave

cavities.

4.5 Implementation and results for a compact ECR ion source

The FDTD model also has been used to investigate the electric ﬁeld distribution
inside a compact ion source [56]. The simpliﬁed ion source geometry for the FDTD
simulation is shown in Figure 4-20. Basically, this ion source is 5.4 cm in diameter and
12.4 cm in height, with a 3.0 cm 1d, 3.5 cm high quartz disk which is used to conﬁne the
plasma discharges. The electromagnetic ﬁeld boundaries for the structure in Figure 4-20
are the sliding short, center conductor, and the outer shell. Again, all the boundaries are
assumed perfect conducting, therefore only the normal component of the electric ﬁeld
exists on these boundaries. The amplitude of an electric ﬁeld component assigned on the
source points is based on the TEM mode wave solution for a cylindrical coaxial
waveguide. After an initial period of a few microwave oscillations (2.45 GHz), the source
term is turned off and the natural response for an empty ion source is observed as in Figure
4-21. The ﬁeld excitation plane chosen for electromagnetic ﬁeld generation is shown in

Figure 4-20. The resonant frequency is calculated from Figure 4-21 to be approximately

90

 

 

 

 

100 . , , r
0 Simulation
x
90* at Experiment ~
0
80- -
K
0’ at:
70" o * -
O
O
60- -
50~ -
0 50 100 150 200

Figure 4-18: Loaded cavity quality factor, Q, versus pressure.

Pressure (Torr)

250

91

 

 

 

 

60 . , , a
0 Simulation
50 ~ It Experiment -
g at
Q 40' "
7:“ a
N
E
3» 30~ -
.3 O
é) %
a; 20— 0 ~
3 X
o
a.
10- 9 r
0 l l l I
0 50 100 150 200

Pressure (Torr)

Figure 4-19: Absorbed power density versus pressure.

250

92

Center conductor
Sliding short

Field excitation plane

Cavity walls

Perm nent magnets

 

 

Ion discharges

 

 

 

 

 

 

 

 

 

 

 

 

F...

Figure 4-20: Simulation structure for an compact ECR ion source.

93

 

5000'

Electric Field Energy
(Relatrve Unrts)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

LUSUour’cecuiffU U U U U U U U

0 . . . . . .
0 100 200 300 400 500 600 700 800 900 1000

 

Time (x 3ps)

2
Figure 4-21: Electric ﬁeld energy,-‘17J‘EOE dV, versus time for the structure shown in
V

Figure 4-20 with no plasma load present. The excitation of the source was
turned off at time = 0.6 nsec.

94

2.6 GHz. The electric ﬁeld magnitude inside the resonating source is shown in Figure 4-22
and Figure 4-23.

The electromagnetic ﬁelds inside the resonant source when a helium plasma load
is present have also been simulated. The conductivity of the magnetized plasma is
determined based on equations (2.49) and (2.50), where the static magnetic ﬁelds
produced by the permanent magnets were simulated by TA. Grotjohn [41]. The plasma
density and electron temperature of the helium discharge are measured by plasma
diagnostic techniques in order to determine the effective collision frequency and
subsequently the conductivity of helium discharges [56]. The calculated conductivity for
the helium discharge has a peak value equal of 178 mho-m around the ECR region where
the microwave electric ﬁeld is perpendicular to the permanent magnetic ﬁeld. The
simulation results of the electric ﬁeld distribution are shown in Figure 4-24 and Figure 4-
25. The electric ﬁelds maintain the general structure shown for the empty source, although
the inﬂuence of the plasma load can be seen. Also, the cavity was simulated as having a
high quality factor, Q value, (several hundred to 1000) for the lower density helium
discharges. For a more accurate simulation, a particle plasma model would need to be
used to simulate the characteristics of this low pressure discharge and ECR heating

mechanism.

4.6 Summary

In this chapter, simulation results for transient state solutions of lossless and lossy
microwave cavities, i.e., the sources are turned off after several microwave excitation

cycles, were presented. These solutions provide information for studying the

95

 

 

 

2000\

1500s
1000~
500~

Ea: @283 E.

lErl (arbitrary units) versus r and 2.

Figure 4-22:

96

he»

                                                
                

4”,

  

_
2:: .
f r. — — ,
aswagzwg
1,: e. .
2,3,3. ..
Z, ”’ 4f‘u‘
€z€$~\“.

     

    

 

222,2 .

 
  

’
g m
its
’41(
Z

 

2000 \

N \ x >
m m w o w
w m 5
ﬁe: SEeS 1m.

lEzl (arbitrary units) versus r and 2.

Figure 4-23:

97

 

 

 

 

1200\

1000‘

m m m m
Ea: Sees 1m.

IBTI (arbitrary units) versus r and z for the compact ion source operated with
a helium plasma discharge.

Figure 4-24:

1000‘

800x

6004

400x

lEzl (Arbitrary Units)

200\

 

98

 

 

 

- ‘ ‘
‘
Elam
"1.1-2:52: :12?” 3
I .

Figure 4-25: lEzl (arbitrary units) versus r and z for the compact ion source operated

with a helium plasma discharge.

99

electromagnetic ﬁeld distributions of empty cavities and the natural frequency responses
of lossy cavities. The application of transient simulations to two different plasma sources

was detailed.

Chapter 5

Steady State Simulation of Electromagnetic Fields
Inside Plasma Loaded Microwave Resonant Cavities

5.1 Introduction

In Chapter 4, all the electromagnetic ﬁeld simulations of empty or lossy-loaded
cavities were based on transient state solutions, i.e., the source is turned off after several
microwave cycles. An alternative simulation approach is the steady state method where
the source is always on. This approach can not be used inside an ideal empty cavity
- because no losses exist. However, if a loss is present inside the cavity, it is possible to get
a steady state excitation when the input power (from the source points) is equal to the
power absorbed by the lossy material. For steady state solutions, the source points which
were described in Chapter 4 are kept oscillating all the time at a frequency equal to or
close to the natural frequency of the lossy loaded cavity. The stored energy increases and
achieves a steady state when the total stored energy is a constant versus time. This occurs
when the rate of power absorption matches the power input. One example is given in
Figure 5-1. The excitation mode, the cavity and the cavity load size are the same as that
used in natural frequency simulations discussed in Section 4.3 (see Figure 4-10). For the
data shown in Figure 5-1, the lossy material had 0 = 1.0 and e = 3.080.

100

101

 

50l- ..

 

40

I
,2:-
:3
:3-

30

I
1

Electric Field Energy
(Relatrve Unrts)

 

 

 

 

 

 

 

 

 

 

 

tel -

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Time (ns)

. 2 .
Figure 5-1: Steady-state response of the electric ﬁeld energy, £71808 dV , versus nme,
V

fOI' the TED“ mode.

102

Another consideration for plasma-loaded resonant cavity simulation is the method
of exciting the ﬁelds. For a microwave resonant cavity, the electromagnetic ﬁelds in the
cavity are excited by an external circuit or waveguide, by means of coaxial-line probes,
loops and small apertures. Moreover, the input power tuning (or matching) methods are
important to adjust the cavity input impedance and control the microwave power reﬂected
from the cavity.

In this chapter, the steady state electromagnetic excitation of discharge loaded
microwave plasma reactors is numerically modeled including the input power coupling
structure. The inclusion of the input coupling structure, which must be terminated in the
simulation at some position, creates an open boundary in the simulation. For this open
boundary, a non-reflecting (absorbing) boundary condition is used for truncation of the
grid points in the FDTD model. The methods to excite input electromagnetic wave for
simulation and to evaluate net power coupling into the plasma-loaded cavity are
investigated. The simulation results of microwave discharge reactor sources used for
diamond ﬁlm deposition will be provided. For simplicity, the discharges used in this
chapter are assumed to be a uniform, cold plasmas by using the simple conductivity model
discussed earlier in Chapter 2. The self-consistent plasma model which is coupled with the
steady state electromagnetic model will be present later in Chapter 6. The techniques used
to experimentally determined the electromagnetic mode and cavity factor Q are

introduced here and also in Chapter 6.

103

5.2 Non-reflecting boundary conditions

The simpliﬁed geometry for plasma cavity plasma reactor simulation in this study
has been shown earlier in Chapter 2 (Figure 2-1). The electromagnetic waves are coupling
into the cavity via an input coaxial probe. The plasma conﬁned in the quartz chamber can
be treated as a lossy material. The cavity walls serve as electromagnetic boundaries where
only perpendicular electric ﬁeld components can exist based on perfectly conducting wall
assumptions. All the electromagnetic waves propagate to this boundary will be totally
reﬂected under this assumption and no ﬁeld can penetrate through these boundaries.
Therefore, the grids used for FDTD simulation can be truncated at these region because no
wave or ﬁeld can exist outside these boundaries. However, at the input end of the coaxial
probe is an open boundary, i.e., the domain in which the ﬁeld has to be computed is
unbounded. The electromagnetic waves, either incident or reﬂected, do exist outside this
boundary and the number of grids needed for simulation becomes unlimited. Therefore, at
this open boundary a truncation method must be used for limiting the domain in which the
ﬁeld is computed. Moreover, this method has to prevent any artiﬁcial reﬂection of the
outgoing wave. Boundary conditions of this type are called absorbing boundary
conditions, or non-reﬂecting boundary conditions. This absorbing boundary condition has
been investigated extensively in solving electromagnetic wave scattering problems
[13][14][15]. One technique which is widely used is provided by Mur [15].

From the FDTD cell shown in Figure 2-2, it can be found that all components of
the electric ﬁeld vector E applied to a particular point on the boundary of the mesh are
tangential to this boundary while the relevant components of the magnetic ﬁeld vector H

are normal to it. For grid points on the mesh boundary, the E ﬁeld components which are

104

tangential to the boundary plane cannot be evaluated by the ﬁnite-difference techniques
since this would require H ﬁeld components that are outside the mesh. Therefore, the non-
reﬂecting boundary conditions for Maxwell’s equations on the mesh require the boundary
conditions for the electric ﬁeld components tangential to the boundary surface. Mur
proposed that each of the E ﬁeld components satisfy the three-dimensional sealer wave
equation independently

(32 32 02 -232

where W is the E ﬁeld component and c is the wave propagation speed, which is equal to
1/ J53 Assuming that the mesh is locate in the region 2 S Ito, and given the boundary
conditions for the plane 2 = k0, then a space-time plane-wave constituent traveling in the
direction of increasing 2, with inverse velocity components 5,, sy, 32 such that

2 2 2 -2 .
sx+sy+sz = c ,can be wrrttenas

W = Re(w(t+s{r+syy+(c-2—si-siy/22D (5.2)

. -2 2 2 1/2 . . . . .
wrth Re(c — 3x — 3,) S 0 , where w rs a function of trme. For this outgorng wave, the

ﬁrst-order boundary condition

a -1 2 2)1/28) _
(5+c (1— (csx) - (csy) 5 W|z=ko — 0 (5.3)

would, for ﬁxed values of 3x and sy, determine a W on the outer surface that is consistent

with an outgoing wave, i.e., it is absorbed. Assuming

105

(l— (csx)2— (csy)2)1/2 = 1+0( (csx)2+ (csy)2), (5.4)

the ﬁrst order approximation can be obtained as

(7+ (lg—JW '2: =0. (5.5)

For the grid termination structure shown in Figure 5-2, which is in three dimensional
cylindrical coordinates, the ﬁnite-difference approximation of (5.5) can be derived using

centered differences in both the space and time increments as

Ef+1(i,j,k k=0) E:(i’j’kO—1)+(f—Tr-Az)/ (ﬁt-+282) (56
.)
x(E:+l(i,j,kO-l)—E,(i.j,ko))
and
”"01.k01= E; (nu/<0 —1)+(j‘—L;-Az)/ (77th) 7
(5.)

x(E"+](i1j,,k0—1)-E;(i,j.k0))

where k0 is the grid terminating point in the z direction. Note that only two E ﬁeld
components needed to be evaluated at the boundary since they are the only components

tangential to the grid truncation plane, i.e., z = k0 plane.

106

 

_ k0 (Grid termination)

 

 

 

 

 

 

 

 

Figure 5-2: Grid termination structure.

107

5.3 Implementation for a coaxial waveguide

The non-reﬂecting boundary condition described in (5.6) and (5.7) was
implemented in a coaxial waveguide in order to verify its validity and accuracy. The
reason for using a cylindrical waveguide as the simulation structure is because the
microwave input power coupling probe discussed in this study can be view as a cylindrical
coaxial probe where the TEM mode electromagnetic waves propagate inside. The
simulation structure for the cylindrical coaxial waveguide is shown in Figure 5-3. The
coaxial waveguide (9 cm in length) consists of an inner cylindrical conductor and an outer
conducting shell, and the electromagnetic waves propagate between these two conductors.
The diameter for the outer shell is 9 cm, and for the inner conductor is 3 cm. The
electromagnetic wave is excited at the z = 0 plane by assigning the time-varying electric
ﬁeld component at the grid points on this plane based on theoretical TEM wave solutions
in a coaxial structure. The electric ﬁeld for the TEM wave in a cylindrical coaxial
waveguide has an E, component only, which is constant in the z and 4) directions and
proportional to a l/r variation in the r direction. The electromagnetic wave will then
propagate to the end of the coaxial waveguide which is at the z = 9 cm plane, and be
terminated by the boundary condition without any artiﬁcial wave reﬂection. One way to
check the direction of wave propagation, as well as, the power ﬂow is using Poynting’s

theorem:

W = §(ExH) . ds (5.3)
S

where W is the total power ﬂow through a surface, S, which is the cross sectional plane of

108

 

z = 9 cm "" "' '—""'""""""""—" _— ' Non-reﬂecting
boundary
condition

2 = 5.9 cm —————————— F - — — a

I-
3
U
8
e. J 4’
O
0
2'3

z: 2.8 cm ————— — -=-- — —————

'5 z
z=0 cm _, Field excitation

 

 

 

 

plane
l——> 1'

Figure 5-3: Simulation structure for the cylindrical coaxial waveguide.

109

the coaxial waveguide, such as the z = 2.8 cm and z = 5.9 cm planes shown in Figure 5-3.
The direction of the surface S is toward the +2 direction for the power ﬂow to the end of
the waveguide.

The simulation grid structure is the same as shown in Figure 4-4, and the total
grids are 30 x 36 x 30 in r, d) and 2 directions respectively. The inner and outer conductors
are assumed to be perfect conductors. The source points on the z = 0 plane oscillate at 2.45
GHz and the time step used here is one picosecond. The simulation results are obtained
after several microwave cycles and are shown in Figure 5—4 and Figure 5-5. In Figure 5-4,
the spatial average B, ﬁeld at the excitation plane (2 = 0) kept the same shape and
amplitude as the one on the absorption plane (2 = 9 cm) all the time. These waveforms
also kept a constant phase difference between the z = 0 and z = 9 cm plane due to the wave
propagation distance of 9 cm. This shows that the waveform generated on the excitation
plane propagated though the Open end of the waveguide without any interference or
distortion resulting from unwanted reﬂecting waves. Figure 5-4 also shows that the
electromagnetic power ﬂow though the waveguide is always positive, and the amplitude
of the power ﬂow is constant. This proves the power ﬂow is toward the end of the
waveguide and no other wave sources exist inside the waveguide and at the open
boundary. Figure 5-5 shows the time-average spatial variation of the r component of the
electric ﬁeld. The 1:“.T is constant in the z direction and follows a 1/r decay in the r direction,
which is the same behavior as the theoretical TEM wave in a cylindrical coaxial structure

[53].

llO

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

   

 

 

 

 

1000 1500 2000 2500 3000 3500 4000

Time (x 10'12 sec)

 

900 T ‘l T I ﬂ I I

 

o
. n' ' ’ O n
. . .
aoo - . .- .. -. . .~ .. .. « ~. .. . . .- .- . -. . 1
-a . . .. -. o . n .- a. .. .. a .g . A ,. o
u :' ,1 .. I. ‘c :. .1 .0 o. -‘ .. ,n '- .0 ._ 't ._ .-
°~ . , . -- '- . . ,- a. -. .. .. 0. .0 .- - -. .4
4! n 't I I l ' . I I ' 9
to ". .' ... 1. '- ‘. a' .n " 'I .. 'I ‘s a ’. t‘ .1 -
.n ' I' . I ' " u' 'e '. a. ' ' 'u 'a .' o' ' e
a , . . ~ ' v - . . u I
700- .. - , .2 . , . . - . t: .2 j. :2 .f . _ _ . - . . ,1 .: .' .-t
' ' a . ' . . -. ‘u ' . n . . a . c . .- ., o . . . - . . .
° ' g . - . ' - ' - . , - ‘ . . .. , . - - ‘. , . 1~ , .
o
:2 .: .. -- -- T' .'- -: -- :. '- -‘ -' “. ,; -.:- 7
. . ' . . . . ' o . - . . ' . .: . . . ' : - . '. _ . . . . . , I .: ~
600 " " " " -1 < " -- ‘. -- .~ .1 ' . . 't .- '. .
.— _ . . - . ‘ ' ' . . . - . - - . . ~ . ' ' . ' , c-r
. . _ a - n u ‘ . . ' _ . - . u . ' . ‘ ' _ .
. 2 u . . . . . : n r - ' . . . . a . . - _ . . . . : ' . ~ - : 2 . ‘ _
. u . D ‘ . ‘ ' . ' . - u ' . u ' . . ' ‘ ‘ ’ I 3 . . ' . ‘ ‘ . a ‘ .
. a . u . ~ . . u . . ' . . . o ‘ . . . _ 1 .
, . n o . _ - u . ' ‘ . ‘ _ , . . .
' ' . . . . . . ' . o . _ a . . . u . o . . o
' . c . ' o ’ I . v . u u a . ’ - ' o a . ~ 6 . o
. u u I ‘ . . . I - . . . . l
. . - . . . . u
' - - - - . r ' . ; . t . - ~ . . - - . . - . - . _ ; ; -
' t. n c ' e - . r - - a u - r . . -
500b ' ' ‘ .1 . ‘ . 1 ‘ d e a ' : » ‘ . a ' : . - a ' - I ‘ ' I s : ' I n v.
- ' ’ ‘ . e n ' . ’ 1 ‘ u . ‘ u ; 1 , ‘ D l , n
b I I
D ' n b r r ,
r b ’

400- ‘

 

Power ﬂow

i l
u ’ ~ ' ' ' I : u i ' t b
u e
p s p a I - I _ 4 l- b . . - t ~
- ' ‘ ‘ 1 ' . . r . 5 - ‘ ' a ' ‘ s - ' . . 1 . I - ' , . ' - l- . . .
. 1 0 ‘ - q ' a . ‘ ‘ . . I ‘
300- : .~ -- , - -: 'r - - : .~ :. . r -- -. - ' .. -- - r- . ‘tﬁ
. . - .1 t . . _ . ,1 .
, c n u ' v . . I . ' . I ° e ‘ ' . . . . o ' ‘ . u '
~ . - a ' _ . . ‘ . ' e ' ’ D I . ' '
n - : . . . , - . . 1 . . u . . . . - . . . . ~ . . u . . . . . . _ _
- I . . o . ~ ' I
c ‘ _ . , _ . o . o . . s . _ . . ' . _ . . . . . . , . . . . , . .
c _ . . ‘ . . . ~ . o I - ‘ u a . o ' . , . u o g - t
‘ . c I - . . - o n c . o .
. _ ' , _ , . . . I . o - n ‘ n . ' . . , . . o u u -
. , . . . . . . . t
_ : . . . . , u . 2 - _ . _ . . ' . _ . . Z . _ . _ . . . 1 o 1 , . u
- ' ’ ' ' . ’ . I o ' ' n ‘ . o ‘ I ' I t . ' I 1 c ' ' ‘ I -
u _ . u . . n . - ' . _ . . . . . . . , . o . . .
u o c ' . . . . a . . c z . - : . . - l _ . . o o . . n : . . . . , . .
‘ ' . . . - o ' o - I t ‘ I a - ‘ ' . a ‘ . . o o ‘ I -
. n 1 . u . . ' . o ' - ' . . i ' ‘

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

100:J
o ...‘

-1090 500 1000 1500 2000 2500 3000 3500 400C

Time (x 10"12 sec)

 

 

 

Figure 5-4: Non-reﬂecting boundary simulation results: (above) spatial average E, vs.

time at two different cross section (solid: z=30, dash: z=0). (below) spatial
average power ﬂow vs. time (solid: z=10, dash: z=20).

111

25\
20~
15x
10~

5~

lErl (arbitrary unit)

 

    

Non-reﬂecting 20
boundary 15

z (x 0.3 cm)

  
 

 

15
10

r (x 0.3 cm)

Field excitation plane

Figure 5-5: Spatial variation (r vs. 2) of the IE,I component for the TEM mode in a
cylindrical coaxial waveguide.

112

5.4 Implementation for a Microwave Cavity Plasma Reactor

The FDTD model is implemented for a microwave cavity plasma reactor such as
used for diamond thin ﬁlm deposition as shown schematically in Figure 5-6. This
cylindrical microwave cavity plasma reactor has a symmetrical construction which can be
excited in single electromagnetic modes. As shown, the reactor consists of a cylindrical
sidewall which forms the outer conducting shell of the cavity reactor which has an inside
diameter of 17.78 cm. The sliding short, the baseplate and the cavity sidewalls, form the
cylindrical excitation cavity. Microwave energy at 2.45 GHz is coupled into the cavity
through the coaxial input probe. The speciﬁc resonant mode is determined by the
geometrical size of the cylindrical cavity and it can be adjusted by the movable sliding
short. Generally, the electromagnetic resonant modes of this type of microwave cavity
plasma reactor for diamond thin ﬁlm deposition are Tan modes.

The simpliﬁed simulation structure of this reactor used in this section is shown in
Figure 5-7. The cavity is assumed to be a simple cylindrical cavity with a movable sliding
short on the top of the cavity. The short is used to adjust the height of the cavity (Ls) for
cavity resonant condition control[1]. A coaxial probe structure, including the inner and
outer conductor, are taken into account in the FDTD simulation region. The cavity
sidewalls, bottom, sliding short and coaxial probe (both the inner and outer conductor) are
assumed to be made of perfectly conducting material, so only the normal components of
the electric ﬁelds exist on these surface. These surfaces make up the electromagnetic
boundaries for the FDTD simulation model, and the only unbounded region (open
boundary) is at the input end of the coaxial probe. The non-reﬂecting boundary condition

discussed in the previous section is implemented at this boundary for grid truncation. The

113

  
 

Sliding Short

 

 

 

 

     

 

 

 

 

 

 

 

 

 

 

 

 

Cavity Walls 3
Outer Conductor
Input Coupling Probe
Quartz Disk
Plasmlal Discharge
Baseplate
Qua z Tube Substrate Screen

Figure 5-6: Conﬁguration of a microwave cavity plasma reactor used for diamond ﬁlm
deposition. The cavity diameter is 17.78 cm.

114

Non-reﬂecting Boundary

 

 

 

Sliding Short Wave Excitation
\ \ Plane (TEM Wave)
_ _ 9 - _
7T?
L.
: Plane for Power L Perfect
Flow Calculation Conductor
5 /
5 2
AL!

 

 

 

 

Plasma Discharges
(Lossy Material)

Figure 5-7: Simpliﬁed simulation structure for a microwave cavity plasma reactor.

115

plasma discharge region, at this point, is treated as a uniform lossy material with a speciﬁc
assigned conductivity and a constant dielectric constant (6:80).

The major difference between this electromagnetic ﬁeld model and the model
discussed in Chapter 4 for cylindrical resonant cavity simulation is, instead of selecting
several points inside the cavity region based on theoretical resonant cavity mode solutions
for electric ﬁeld generation, the electric ﬁeld here is provided by the coaxial probe which
is outside the cavity region. Hence, this method allows the simulation of the real
electromagnetic coupling mechanism for discharge excitation in the microwave cavity
plasma reactor. The technique used to excite the electromagnetic ﬁeld in this numerical
model is to select grid points on a cross sectional plane of the coaxial probe as source
points and assign the time-varying electric ﬁeld component at these points based on
theoretical TEM wave solutions in a coaxial.structure. The electromagnetic wave then
propagates down into the cavity region where power is absorbed by the discharge. Any
reﬂected electromagnetic wave will propagate to the end of the open boundary of the input
power probe and be terminated by the non-reﬂecting boundary condition without
producing any artiﬁcial wave reﬂection. In order to prevent the reﬂection of waves from

the source points, the electric ﬁelds on the source points are assigned as [13]

+ 1 . . . . .
E: (1,},k1) = E:(l,j,k1) +Csrn(21tft) (5.9)
where C is the amplitude of the electric ﬁeld source, k1 is the grid points where the source
plane is located, and f is the excitation frequency.
For the microwave cavity plasma reactor simulation discussed in this section, the

inner diameter of the cavity is 17.78 cm, and the excitation frequency is 2.45 GHz. The

116

discharge (or lossy material) is located at the center of the bottom of the cavity and is
assumed 8.85 cm in diameter and 2.0 cm in height. The conductivity and the cavity height
are treated as input parameters to investigate the steady state electric ﬁeld patterns and
resonant characteristics of the loaded cylindrical cavity. The number of grids used in this
simulation is 30 x 36 x 45 in the r, d) and z direction respectively, and the time step is 0.5
psec.

The simulation results of a 14.41 cm high cavity loaded with a 0.5 mho/m lossy
material are shown in Figures 5-8 to 5-10. This height (14.41 cm) is basically the height
for the theoretical TMmz mode in an empty cylindrical cavity with the excitation
frequency and diameter mentioned previously. This empty cavity mode is described by
equation (4.1). Figure 5-8 shows the electric ﬁeld energy stored in the cavity versus time.
After several microwave cycles, the power transfer to the cavity is equal to the power
absorbed by the lossy materials. When this balance occurs, the amplitude of the electric
ﬁeld energy is constant and the system has reached its steady state. The spatial electric
ﬁeld distribution in the r-z plane are shown in Figure 5-9 and Figure 5- 10. The simulation
results are in agreement with the theoretical TMmz mode pattern as shown in Figure 5-11.

Especially for the E, component in the z direction, which follows a sin(ZTm-) variation

3

proﬁle as expected for the theoretical TMmz mode variation in the z direction. Inside the
discharge region, the electric ﬁeld strength decays due to microwave power dissipated in
the lossy materials (or discharges). Moreover, the calculated Br and E2 components are
symmetric in the (1) direction, and the E¢ components remain zero all the time during the
simulation. Therefore, the three dimensional FDTD model used here can be further

simpliﬁed to be a two dimensional model, that is, the simulation domain are considered in

117

2% v v V V 1 v v v f

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1800-

1600'
£9 A 1400.
o 2
LS 1: 12001
“c D
7, g 1000. A
LI; 3 800. v
.: é)
a t. 600. l
2
1.1-J 4001

200’

0

 

0 100 200 300 400 500 600 700 800 9001000

Time (x 10 psec)

Figure 5-8: Electric ﬁeld energy vs. time for W012 mode cavity simulation.

118

ltrary uni

Sliding Short

t8)
—‘-‘ ”(900
88§888
[Ill/J

Cavity Walls

IEZI (arb

,8

 

 

 

 

 

 

 

 

 

 

 

0) \ = "0'0 o t
20 » ‘\ o’o':v%“v:,o:t: «9:. -
O O O O. 0.0” I” 0’ o‘
/ / ll ’0 -
- go‘ 02o? 293;.“ ¢::::‘
9 ’1 :I’v’" \x‘a‘ ‘ ‘ “
15 ., ,;/ ” ””’ “‘ “‘ “ 1 0
, o'f’3’“( “ ‘\“‘
In ut Cou lin " r“ , ’ xo'o'o‘. ‘. "‘ﬁ ‘ 8
P P g 335...... 0‘
Pro '0‘.“.“.‘ 6
«.355»
.O
‘ 2 r (cm)

 

Discharge Region

Figure 5-9: lE,l ﬁeld distribution for a TMmz mode microwave cavity plasma reactor.

119

 

 

 

    
 

 

 

 

300\

250~
:1? 200v 4A
8
E 150‘ Sliding Short
.5
5 100\ .. - Cavity Walls
m o’;:‘:“““‘:“:‘33:'9‘
— 50 V \ oﬁégﬁe:::\“‘€\‘l‘}}z

20

. 15 10
Input Coupling
Probe

Discharge Region

Figure 5- 10: IEZI ﬁeld distribution for a TM012 mode microwave cavity plasma reactor.

120

_, E - Field

-——-> H-Field

 

 

 

 

 

,al, in“)

 

 

Figure 5-11: Field patterns of TM012 mode cylindrical cavity.

121

the r and z direction only. This model can be built based on solving the two-dimensional
FDTD equations developed in Chapter 2 ((2.27) to (2.32)).

By using the two-dimensional FDTD model, the quality factor Q of the cavity,
which is determined from the stored electromagnetic energy and the power absorbed by
the lossy materials as discussed in Chapter 4, can beanalyzed for various discharge
conductivities and cavity heights. The simulation results of Q vs. cavity height for three
different conductivity materials are Shown in Figure 5-12. For a load with a given
conductivity, the cavity Q changes as the cavity length varies. When the Q value drops,
that means power is absorbed more efﬁciently by the lossy load. Figure 5-12 shows that
the lowest Q value occurred at a cavity height of about 14 cm, which is close to the mom
resonant empty cavity height of 14.41 cm. The lowest Q value doesn’t exactly occur at the
theoretical TM012 mode empty cavity height (14.41 cm) due to the loaded lossy material
changing the resonant behavior of the cavity. This phenomenon has also been shown
earlier in Figure 4-15, where the resonant frequency of a cylindrical cavity changes with
the conductivity of material loaded inside the cavity. When the cavity height increases
beyond 14 cm, the Q value increases, which means less power is coupled into the load and
the more power is reﬂected from the cavity. The increase of reﬂected power is shown in
Figure 5-13. The reﬂecting power is determined by using Poyting’s vector in (5.8). The
cross section plane used for determining the power reﬂected from the cavity is the open
boundary surface and the direction of the surface is toward the +2 direction for the
reﬂected power ﬂow out of the cavity. When the height approaches 21 cm, the Q value
drops again because the cavity is reaching another resonant condition, namely, TM013

mode, where the theoretical resonant cavity height is about 21 cm. Figure 5-12 also shows

122

 

 

250 . . . . i
-- -- 0.1mho/m '
- - 0.5 mho/m .'
200 ' — 1.0 mho/m ‘
150 I "' ‘I‘
O’ .
100 I- a 4
K.
50 P \ 1

 

 

 

10 11 12 13 14 15 16 17 18 19 20 21
Cavity Height (cm)

Figure 5-12: Q vs. cavity height for three different conductivity loaded TMmz mode
cavity simulation.

123

 

0.12 - - . . _ . . . - f
0.11 i ...... ......... N“'\.\.‘ \\ ff’wﬂﬂga ....... ~\ \
.‘2 \\ \ //'/
g 0'1 . l I/ / .
l l . l
l x
g 0.09 t it \ ’1’, I
3 l \ I: I
:3 0.08 h ‘1‘ \ :3 I
g \ l; I
E 0.07 l- \\\(;{\ I! _‘. _ 0] th/m q
H ----~ 1.0 mho/m
3.2. 0.06 . l l ‘
° \ I
m \I
0.05 . ' ‘

 

 

 

0'0410 1'1 1‘2 1‘3 1'4 1‘5 1‘6 f7 1‘8 1'9 2'0 21

Cavity Height (cm)

Figure 5-13: Reﬂected power vs. cavity height for TMmz mode cavity simulation.

124

the resonant behaviors for different conductivities. When 0 = 0.1 mho/m, less power is
absorbed by the loaded material, i.e., less perturbation for the resonant cavity, and the Q is
higher and cavity height for resonant condition is about equal to that of the TM012 mode
empty cavity. When 0 increases, more power is absorbed by the loaded material.
Therefore, the Q values drop, and the cavity height is starting to shift away from the empty

cavity height in order to achieve the resonant condition for the lossy loaded cavity.

5.5 Simulation results for a cavity reactor loaded with a H2 discharge

Another application made of this FDTD model was to simulate the power
absorption and electric ﬁeld distribution in a microwave cavity plasma reactor which was
investigated by Zhang, et. al [3]. In Zhang’s study, a 17.78 cm i.d. microwave cavity
excited at 2.45 GHz was loaded by a hydrogen dominated plasma in a 9.25 cm i.d., 4.35
cm high quartz disk. The radial component of electric ﬁeld strength at the cavity outer
wall was measured by a micro-coax electric probe inserted through the cavity wall versus
2 and q). The measured ﬁeld distribution was compared with theoretical results to verify
that the mode sustained in the cavity with the discharge present was approximately a
TMOU mode. The discharge inside the quartz disk had volumes as estimated visually of 30
cm3 to 100 cm3 depending on the pressure which was varied from 20 Torr to 70 Torr.
Zhang also found the quality factor, Q, of the cavity by measuring the absolute microwave
electric ﬁeld amplitude and the electromagnetic ﬁeld mode distribution to determine the
stored electromagnetic energy according to (4.2). The Stored energy, U, and the input (or
total absorbed) microwave power, P, then determined the Q value as Q=coU/P. In

particular, Q values determined this way were about 60 and they varied little over the

125

pressure range studied (20 - 70 Torr). Additional details of the experimental Q value
derivation and coax-probe calibration method will be discussed in next chapter.

The simulation structure of this reactor is shown in Figure 5-14, which is similar to
the structure used in the previous section. However, in this section, the baseplate included
in the numerical model is a more realistic reactor structure. The substrate holder is
assumed to be aligned with the bottom surface of the baseplate, and to have approximately
the same diameter as the quartz disk. Therefore, the sliding short, the baseplate, the
substrate holder, the coaxial probe (both the inner and outer conductor), and the cavity
sidewalls form the boundary of the FDTD simulation region. The boundary conditions for
the electric ﬁelds on these surfaces and the open boundary condin'on are the same as in last
section. A discharge is excited and sustained inside the disk-shaped quartz dome region by
the input microwave power. Moreover, the quartz disk is included in this model and
assumed lossless with a relative permittivity of e, = 3.78.

The simulation of electromagnetic ﬁelds requires knowledge of both the electron
density, he, and effective collision frequency, ve/f, in order to solve the electron
momentum transport equation and induced current density ((2.53) and (2.54)). A plasma
discharge description which includes several idealizing assumptions has been adopted for
this section. The discharge is assumed to have a cylindrical volume shape with a size
found from visual observation. The density and effective collision frequency within this
volume are assumed uniform. The idealizing assumptions ignore such behaviors as gas
ﬂow, non-constant plasma temperature proﬁle and plasma density variations. Hence, care
must be used in interpreting the simulation results in the plasma discharge on a local or

small scale. However, the macroscopic quantities such as electromagnetic ﬁeld mode, Q,

126

Non-reﬂecting Boundary
Wave Excitation

'lane (TEM Wave)

 
       

   
   

     
  
 

v1
9
d

.9

   
   
 
 

v
9..

V

A

    
  
  

V

"
ﬁ§

l

.0

  
  
  
 
 

      
  
  
 
 

      
   
 

      
  
 

   

Plane for Power
Flow Calculation

i?!
I...

   

v
9

  
  

  
   
   
 
  
   

O

.o

V
O

     
  

O

H

.3?

v5
6.0:.
c»

6

   
   
   
   
   
  
  
   
  

v
0

O
o

v
Q

   

O
o

,‘

v

%’

   

vv
'9
.‘O

     
 
 

. O

  
    
 

 
    

  
       

  
   
   
  
   
  

 

   

    
    
 

        
  
   
   
  
  
  

     
  

 

   
  
   
  

   
  
 
 
 
  

    
   

S .

Quartz Dlsk (8 r - 3.75) Perfect
\ Conductor
‘p\ \\\\\\\\\\\\\§
\
‘ \
~ \
‘ \
~ \

1. ~ ‘

'vvvvvvvvv 'VVV'IVVVV'

ﬁﬁéﬁéﬁéﬁﬁg iﬁﬁﬁéﬁﬁﬁQQQ

9909999990 eodboeooo-

ﬁﬁﬁﬁﬂﬁﬁﬁﬁ’ ’Qﬂﬁﬁﬁﬁﬁﬁﬁﬁ

9990009090 »oooeeeoo¢¢

0999999999 ,909999990. 2

ﬁﬂﬁﬂﬁﬁﬁﬁﬁ’ iﬁvvvvvvvvﬁ

9909900000 rooeeooeeo

$§WV~VV¢V§ iJVVVVVVVVQ

.A.A.A.A.A.A.A.A.A.A\ A’A’A’A’A.A.A.A.A.A'

Plasma Discharge | I r Substrate Holder

Figure 5-14: Simulation structure cross-section for microwave cavity plasma reactor.

127

and absorbed power are expected to be simulated accurately by using the measured values
for plasma quantities. A more accurate plasma discharge model in the next chapter will
remove this restriction.

From the measured gas temperature [54] and the corresponding gas pressure, the
effective collision frequencies in hydrogen discharges at different pressures can be
calculated approximately as v eff = 1.44 x 1012 x {1: (P is the pressure in Torr and T is the
temperature in Kelvin) according to [50]. Treating the discharge volume, collision
frequency and electron density (assumed uniform) as input data, the microwave
absorption properties can be evaluated by (2.35), and the electric ﬁeld distributions, power
deposition distributions and the quality factor Q of the cavity are simulated using the
FDTD model. By calculating the time average energy stored and the time-average power
dissipated in the cavity using equations (2.35) and (4.2) after the simulation reaches steady
state, a quality factor can be determined. Moreover, at steady state, the power dissipated in
the cavity is equal to the net power transfer to the cavity, which can be calculated by
Poynting’s theorem given by (5.8).

For the hydrogen discharges in this study, the electron density has not been
measured experimentally yet with high accuracy so it is treated as a variable parameter.
The cavity quality factor Q is simulated versus electron density under different pressure
conditions using the corresponding measured gas temperature and discharge volume. The
results are shown in Figure 5-15 and Figure 5-16. These simulations are done by using the
two-dimensional FDTD model since ¢ symmetry exists for the TMon mode cavity. The Q
values calculated by using either the power dissipation model (2.35) or the power ﬂow

model (5.8) are consistent with each other. For these simulations, when the electron

128

 

75 ...”..- . ...“--.
: Power Flow Model

 

7O

............. : Power Dissipation Model

I

65

60-

50*

45L

35*

 

 

 

30
10

Plasma Density (cm'3)

Figure 5-15: Q value vs. electron density at a pressure of 20 Torr.

 

10

129

 

 

 

 

 

so - . - . ~ - - ., . -
70 ~ : Power Flow Model
I ------------- : Power Dissipation Model
60*
O’ 50-
40*
30)-
20 A A A A A A A A l A A A A A
10" 10‘2 10

Plasma Density (cm'3)

Figure 5-16: Q value vs. electron density at a pressure of 70 Torr.

130

density increases, the conductivity of the discharge increases. At lower electron densities
the increase in conductivity results in increased total power absorbed. As the electron
density increases further, the Q value reaches a minimum and then increases. The increase
of Q occurs because the plasma has a high conductivity, and hence power absorption
occurs primarily at the plasma boundaries and not throughout the entire plasma volume as
the microwave skin depth decreases. The skip depth is deﬁned as the depth at which the
magnitude of the penetrating electric ﬁelds decrease to l/e (about 36.9%) of their value at
the boundary surface [53]. Following this deﬁnition, the skin depth can be obtained from
the simulation data by examining the magnitude of elecnic ﬁeld inside the discharge and
at discharge boundary. For example the simulated skin depth for transverse waves is about
0.7 cm at a plasma density of 3 x 1012 cm’3. Theoretically, the skin depth 8 for transverse

waves penetrating into plasmas is given by [62]
8 = c/(ope _ (5.10)

where (ope is the plasma frequency deﬁned in (2.46). Equation (5.10) is valid under the
assumption (0 « (ope. For plasma density of 3 x 1012 cm'3, rope is 9.75 x 1010, then the skin
depth is calculated to be 0.3 cm, which is on the order of the simulation result. By
matching the calculated Q values with the experimental data (~ 60), the electron densities
can be predicted around 3 x 1012 cm'3at 20 Torr and about 8 x 10 12 cm'3 at 70 Torr. At a
Q value of about 60, the higher plasma density value is selected based on initial
experimental measurements[54].

The spatial electric ﬁeld distributions in the r-z plane are shown in Figure 5-17 to

Figure 5-20 for the discharge maintained at 20 Torr with an electron density equal to 3.0 x

131

Quartz Dome

Sli ing Short

o I
o a. o o , z,
.. ......vazo,»
a C3.
aoao
co

2 o
2 3.

o
9’30;

¢oo 00 ,3: O

oo o

 

300\

200\
100\

35. bases 3.

o
ooso
...“...2............".. s . .
of , \\ ‘ ".

‘ O I

/ \
aeratev‘OOO

Baseplate

Discharge Region

o a
3

9,

ooo a
2’

¢

« I

a
’69

i

I
I”...

Input Coupling 3
Probe

3 x1012

the r - 2 plane at 20 Torr for a plasma density

rstribution in

115,! d'

Figure 5-17

3

cm

132

§

  
     

Quartz Dome

 
  
  
 
  

E 200 Sliding Short

3

E‘ 150

,5 1 Cavity Walls
5.

—N 50

a;

—L
No

  
  

10 10

Input Coupling 3
Probe

   
 
 

6

4

z (cm) Baseplate

0 0
Discharge Region

Figure 5- 18: IEZI distribution in the r - 2 plane at 20 Torr for a plasma density = 3 x 1012

cm ’3.

133

Sli ing Short

0
0’4
o
0"
Ox
//
//

Quartz Dome
:j.
0"
O

 

300\

200~
100~

Discharge Region

2 (cm)

Probe

IB,I 'stribution in the r - 2 plane at 70 Torr for a plasma density = 8 x 1013

Figure 5-19

cm '3.

134

250

 
  
  
  
  
 
 
  

200 Quartz Dome

Sliding Short

{9‘

Cavity Walls

IEZI (arbitrary units)
29 £33

0

0.12 0.1

     
  
 
  
  
    

0.1

Input Coupling”-O8
Probe 0.06

0.08

  

004 0-04 Baseplate

2 (cm)
0.02 0'02 r (cm)

0 0

Discharge Region

Figure 5-20: IEZI distribution in the r - 2 plane at 70 Torr for a plasma density = 8 x 1013

cm '3.

135

1012 cm'3, and for the discharge maintained at 70 Torr with an electron density equal to
8.0 x 1013 cm'3 respectively. The calculated E, and 152 components are symmetric in the (1)
direction, and the 13,, components are zero during the simulation. The simulation is in
approximate agreement with the theoretical TMOH mode wave pattern. In particular, from
these ﬁgures, the distribution of E, and E2 ﬁelds along the z and r directions are in
approximate agreement with the theoretical electric ﬁeld distribution for the TMOH mode
which was the experimentally measured mode[3]. Due to the presence of the discharge,
the electric ﬁelds decay further into the plasma region as microwave energy is absorbed by
the plasma. Examples of absorbed power distribution for the substrate location indicated
in Figure 5-14 are shown in Figure 5-21 and Figure 5-22 for 20 Torr and 70 Torr

conditions respectively.

5.6 Conclusions

The FDTD numerical model has been applied to solve the steady state
electromagnetic ﬁelds inside a plasma loaded microwave cylindrical resonant cavity in
this chapter. The simulated structure included the input coupling power probe, as well as,
the geometry of the base plate, substrate location, and quartz disk structure. The
simulation results of a hydrogen discharge loaded plasma reactor, as compared to
experimental data, showed good agreement for the electric ﬁeld distributions. The FDTD
electromagnetic model described in this paper will next be coupled with a more complete
plasma continuum model to obtain a more accurate simulation of microwave plasma
discharges, especially the spatial variations in the plasma. This self-consistent

electromagnetic ﬁeld and plasma model will be discussed in the next chapter.

136

      

Base Plate

r (cm)

   

o . o o o
" Xx»; «a. x?

    

    
      

 

   

 

 

o
3"...) 4.3/7”
o0... ooooo oo 00%
o
o

    

2

       

 
  

3x,

             
     

Discharge Region

    

   

o o o o: o o
o’owpnbuooHJHMonooMoH'. "Won“ bore
o» o o» a I».
a o o o o o a o
32:. I to .4

         

  

a
x, o 3:

         

 

Quartz Dome

 

 

150\
100x
50

0
12

m.
2

3.8: Embﬁav bacon— 330m

Power density distribution in the r - 2 plane at 20 Torr for a plasma density

0
0

Figure 5-21

3 x 1012 cm '3.

137

     

Base Plate

 

o
o o 0
.... 2...”... 0

    

 

 

Discharge Region

Quartz Dome

 

% m. m. m new
432

$55 saga—EV bacon 530m

i

Power density distribution in the r - 2 plane at 70 Torr for a plasma dens

= 8 x 1013 cm '3.

Figure 5-22

Chapter 6

Self-consistent Simulation of a Microwave Cavity
Plasma Reactor

6.1 Introduction

In this chapter, a self-consistent model is presented which couples the FDTD
electromagnetic ﬁeld model and the ﬂuid plasma model together in order to investigate the
electromagnetic excitation of discharges inside a microwave cavity plasma reactor. The
techniques used to coupling these two model are discussed in detail. The simulated plasma
reactor is loaded with a H2 discharge which is used for diamond thin ﬁlm deposition.

Simulation results are provided and compared with diagnostic results.

6.2 Model description

The numerical model used in this chapter for microwave cavity plasma reactor
simulation combines the FDTD model described in Chapter 2 and the ﬂuid plasma model
described in Chapter 3. The FDTD model which solves the electromagnetic ﬁelds inside
the discharges can provide output information such as the microwave power absorption of
the discharges by using an appropriate plasma conductivity model. This absorbed

microwave power is the input to the ﬂuid plasma model, which is the heating term of the

138

139

electron energy balance equation (3.33). The discharge characteristics such as the plasma
density and electron temperature are determined in Steady state under this power
absorption condition. The discharge characteristics information is then sent to the FDTD
model to modify the plasma conductivity and calculate a new discharge power absorption.
Therefore, in this iterative manner, the power absorbed by the discharges will converge to
a stable value, and the electromagnetic ﬁelds inside the reactor source and the plasma
characteristics can be solved self-consistently. The plasma conductivity model which is
the link between the FDTD model and the ﬂuid plasma model determines the relationship
between the microwave ﬁeld and induced microwave current inside the discharges. This
relation can be determined by solving the electron momentum transport equation
(Langevin equation) equations (2.39) and (2.40) in the time domain as stated in Chapter 2.

The ﬂuid plasma model describes the characteristics of discharges by solving the
steady state electron continuity equation, ion continuity equation, electron energy balance
equation, and Poisson equation (3.30)-(3.36). Finite difference techniques are used to
discretize these equations in cylindrical coordinates, and the resulting discretized
equations were shown in (3.37)-(3.43). To solve these non-linear type, discretized
equations, both direct method (Newton’s method) and iterative method techniques are
used (see section 3.4.2). In particular, the Poisson equation, electron continuity equation
and ion continuity equation are tightly coupled and solved by Newton’s method. The
unknowns solved at each grid point are the electron density, ion density and plasma
potential, and the electron and ion ﬂux (steady state DC ﬂux) can be determined by (3.34)
and (3.35) using the solution of these unknowns. Then, the electron density and ﬂux are

coupled into the discretized electron energy balance equation (3.41) to solve the electron

140

temperature. Another input for the discretized electron energy balance equation is the
absorbed microwave power density Pubs which comes from the FDTD electromagnetic
ﬁeld model. The calculated electron temperature is then feedback to the continuity
equations and Poisson equation to modify the reaction rates and parameters and update the
electron density, ion density, and plasma potential. The ﬁnal stable solution of density and
temperature is achieved by iteratively solving the continuity/Poisson equations and energy
balance equation.

The ﬂow chart of the microwave cavity plasma simulation is shown in Figure 6-1.
The initial plasma density and potential data are determined by solving the Poisson
equation and electron and ion continuity equations with some given constant reaction
rates, such as the ionization rate and recombination rate. This information is coupled with
initial absorbed microwave power, which might be assumed constant, to be the input data
of the energy balance equation of electrons. Then, as previously discussed, the steady state
electron density and temperature are determined by repeatedly solving the ﬂuid plasma
equations. These two quantities are then used to evaluate the microwave induced
conductivity of the discharges in the FDTD model. Next, new power absorption data is
obtained after a few microwave cycles of FDTD simulation, and it is coupled back to the
ﬂuid plasma model to update the plasma density and electron temperature. These
processes form a closed loop and are repeated until the steady state absorbed microwave
power converges to a stable value.

Since the solution of electron temperature and plasma densities are loosely
coupled, the damping technique discussed in Chapter 3 is used during the simulation in

order to prevent divergence phenomenon from occurring. The damping factors for

141

‘ "er Je’ P absorbed (initial) )

Solve electron energy balance equation

 

 

’1

 

—' to determine electron temperature (steady state)

 

 

 

l Te(damped)

 

Solve Poisson equation and continuity equations
to determine particle densities, ﬂux and plasma

potential (steady state) by Newton’s method

 

 

 

   

No

 

  
   

"er Te
converge?

 

ne, Je (damped)

he, Veﬁr

 

Solve Maxwell’s equation by FDTD model
to determine the power absorption
(time domain: three microwave cycles)

 

 

 

 

 

P absorbed (damPCd)

 

Figure 6-1: Flow chart for microwave cavity plasma reactor simulation.

142

electron temperature and density is normally 0.01, and might reduce to 10'5 if severe
oscillation or divergence occurs. The damping factor of the absorbed power density, Pubs,
which is the input to the ﬂuid plasma model, is normally 0.1 to 0.01. For solving the
Poisson equation and continuity equations, since they are tightly coupled by the Newton
method, the damping factors are normally set to 1, and might be reduced to 0.01 if

divergence occurs.

6.3 Reactor description and Q value measurement

The basic structure of the microwave cavity plasma reactor for diamond thin ﬁlm
deposition investigated in this chapter is similar to that in Chapter 5, which is shown in
Figure 5-6. The inner diameter of the cavity is still 17.78 cm, but the height of the cavity,
L, is usually adjust to 20.4 cm for TM013 mode excitation with 2.45 GHz input frequency.
The plasma discharge is conﬁned by a quartz disk which is 14.1 cm in diameter and 10.3
cm in high. The substrate holder and wafer size is 10.4 cm diameter, and the location of
the substrate holder can be adjusted in the z direction.

The electromagnetic resonant mode and cavity quality factor Q for this cavity
reactor are experimentally investigated by inserting a micro-coaxial probe into small holes
drilled though the resonant cavity walls. The center conductor of the probe is inserted
2mm beyond the inside surface of the cavity side walls, and the other end is connected to a
power meter during the experiment, as shown in Figure 6-2. This probe samples a small
amount of electromagnetic power near the cavity inner side walls which can be measured
by the power meter. Since the probe power reading is proportional to the square of the rrns

electric ﬁeld strength normal to the inside cavity, the cavity resonant mode and the cavity

143

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

SMA Connector 17.3 mm
Power ‘
meter H
:l
. Center
Microcoax ‘__.l Conductor
Cable 13 mm
Collar
1
Holes > [:1
E
/ I:
1:]
Cavity Walls [:1

 

~\

 

+N

 

 

 

 

 

 

 

Figure 6-2: Microcoax probe and holes in the cavity walls.

144

stored energy can be determined using this techniques.

In order to establish the relationship between the power reading of the probe and
the strength of electric ﬁeld, a calibration procedure has to be performed. In this
procedure, the calibration of the probe is done by using a 17.78 cm i.d. empty cylindrical
brass cavity which is excited by a 2.45 GHz microwave frequency. This is the same brass
cavity used on the diamond deposition system to take electric ﬁeld measurements. The
bottom of the cavity was covered with a brass plate in order to form a microwave resonant
cavity. The cavity height is adjusted to 20.4 cm for TMUZ mode excitation. The TM 112
mode was selected Since it gave the highest empty cavity quality factor (about 5000). This
quality factor was determined by the resonant frequency f0 and the bandwidth Af of the
reﬂected power measured from a frequency sweet of the resonant cavity, i.e., Q = f0/Af.
By inserting the micro-coax probe into the hole at z = 17.3 cm above the cavity bottom
surface on the cavity side walls, the probe power (Pp) was read from the connected power
meter. In order to established the relation between the electric ﬁeld strength E near the
cavity inside cavity walls and the probe power reading Pp, the relationship between E at
the probe sampling position on the cavity walls and the absorbed power Pa has to be ﬁrst
derived. The power absorbed by the empty cavity, Pa’ is deﬁned by Pa=Pi-P,, where Pi is
the incident power imparted to the cavity (approximately 100 mWatts) in this case and Pr
is the reﬂected power from the cavity [48].

Using perturbation theory for the lossy conducting surfaces and assuming that the
power absorbed by the empty cavity is dissipated on the cavity walls, the following power

balance equation is established [49],

145

Pa = R§|H|2ds (6.11)

where R is the surface resistance. For brass at 2.45 GHz excitation frequency,

R = 5.01 x 10‘7 f0 = 0.0248. (6.12)

For the W112 mode, the components of the H ﬁeld are:

l r . 21tz
Hr = —A;Jl(klla)s1n(¢)cos(-T) (6.13)
H — AL‘J' (i. I) ( (LIZ) (614)
(D — — a l “a cos (t))cos L .

and H z = 0 , where A is the amplitude constant, and 2.11 is ﬁrst root of the Bessel function
J1. For a = 8.85 cm (cavity radius), L = 20.4 cm, and 2.11 = 3.83, equation (6.11) can be

evaluated by (6.12) to (6.14) giving

2nL 21w
Pa = R]“11¢(r:a)|2r(d¢)dz+2RH|H¢(z=0)|2r(d¢)dr
O 0 0 0

21m (6.15)
+ 218“]ng = 0) |2r (do) dr

0 0
= 03993.42

This evaluation was performed using numerical integration. For a cylindrical cavity, the
only component the electric ﬁeld normal to the cavity side walls is the B, ﬁeld, which can

be written for the TM] 12 mode as

146
1 21: A11 r . 21tz
Er = —A(jw——E)(T)(7)JI(KHE)COS ((11) 8111(T). (6.16)

For r = a, q) = 0, and z = 17.3 cm, (6.16) can be calculated to be:

2
= 1.0344 x107A . (6.17)

r=a

 

I542

From (6.15) and (6.17), the relationship between the electric ﬁeld strength E near the

cavity inside walls and absorbed power is

= 2.59 x 107100. (6.18)

r=0

 

lErl2

Since the ratio of Pa to Pp, which is equal to 3243.24, is a constant, the relationship
between the electric ﬁeld strength E near the cavity inside walls and the probe power is

obtained as

152 = IEJZ = 8.4><1010xPp (6.19)

 

r=a

where the unit of E is Volt/m and the unit of PI) is Watt.

After the probe calibration is completed, the electric ﬁeld strength E near the inner
Side walls of the microwave cavity plasma reactor can be determined from the power
reading of the micro-coax probe. Probe experiments are done by measuring the electric
ﬁeld Strength at various 2 positions for the microwave cavity plasma reactor during
diamond thin ﬁlm deposition. The microwave excitation frequency is 2.45 GHz, and the

cavity height is adjusted to a position near 20.4 cm for a minimum reﬂected power

147

condition. The result is Shown in Figure 6-3 where the resonant condition for the
microwave cavity during deposition process matches closely to the ideal TMOB mode.
The ﬁeld pattern for the TM013 mode is shown in Figure 6-4. Additionally, the electric
ﬁeld was probed circumferential around the cavity. The ﬁeld Strength showed no variation
along the (1) direction as expected for the TMOB mode.

The Q factor of the plasma discharge loaded cavity reactor can be estimated by
calculating the electromagnetic energy stored in the cavity and the power dissipated in the
cavity. The energy stored in the cavity is determined under the assumption that the
electromagnetic ﬁelds are only slightly altered by the presence of the discharge and that
the electromagnetic ﬁelds in the discharge are not much different from those when the

discharge is not present. The electric ﬁeld for the theoretical TMOB cavity mode is given

 

by:

E_A3\'0131rjxr.3rtz 620

'_j(t)ea_L_l°1c-rsm—L_ (°’

E¢ = 0 (6.21)
2.2

Ez = A2):—12]o()‘01g)005(§'g‘z) (6.22)
1 a

where 9.01 is the ﬁrst root of Bessel function JO, and A can be determined in (6.20) and
(6.22) with the measured probe power Pp. Thus, the energy stored in the cavity can be

determined by:

148

 

 

 

 

2.5 I I I I I I I I I I

E :1: Experiment Theory

2 2~ x a

>
<1-

0 it: at : it

5 1.5 - . ~
5

5.0 '.

G ‘.

0 all . .

b '- .-

m _ 3* _-‘ . ..
2 1 .- ..

.2

LL 3" :1 3‘

.o ‘2

2 0.5 - '2. '

. ' at
0 .' 1 J 1 AA 1 1 1 J 1 1 1
0 2 4 6 8 10 12 14 16 18 20
Cavity Height (cm)

Figure 6-3: Electric ﬁeld strength vs. cavity height for microwave cavity plasma
reactor during diamond deposition process.

149

-——-> H-Field

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 6-4: Field patterns for TM013 mode.

150

U = we = ej|E|2dv = ej|Erlde+ej|Ez|2dv (6.23)
and the Q factor is given by,

_ wU .
Q _ 7 (6.24)

a

where P,l is the power absorbed with the plasma discharge present.

6.4 Implementation for H2 discharge simulation

The simulation structure for the FDTD model and the plasma ﬂuid model are
shown in Figure 675. The simulation structure for the FDTD model is basically the same
as that discussed in Chapter 5, but the quartz tube which is used to adjust the height of the
substrate holder is included in this chapter. The assignments of the perfect conductor
boundary condition and non-reﬂecting boundary condition are the same as in Chapter 5.
The location of substrate holder is again assumed as aligned with the upper surface of the
base plate. The boundaries of the ﬂuid plasma simulation are the quartz disk and the
substrate holder, where the diameter of the substrate holder is the same as that of the

quartz disk. The boundary conditions of the ﬂuid plasma model are:

an

e
r

atr = 0: = — = 0 (6.25)

QJIQJ
" -E
Q.)
$1

an
- I?

and

151

tion

'lane (TEM Wave)

Conductor

m
.n
e
P

 

Wave Exci

  

3????333
.9099...
.9990099.
vooeoeoe
0069009.
ooeoeooe
.ooeoeoo.
9909999:
.eooeoeo.
v0.9.9.0
oooeooo.
999999..
99900.0.
I??33333
99999$99

 

 

     
      

      

L

nonououou nouononouononouon.

23232
90950
ﬁﬁﬁﬁﬁ
HPQSPA

   

  

N on-reﬂectin -

Quartz Disk (er=3

 

 

, Boundary
75)
I r

Simulation structure cross-section for microwave cavity plasma reactor

with quartz tube.

Figure 6-5:

at the quartz wall and substrate: ne = n‘. = ‘1’ = 0
(6.26)
T? = Tn

At r = 0, which is the center of the cylindrical reactor, the boundary condition (6.25) is
based on assuming all the unsolved quantities (ne, ni, Te, and ‘1’) have an r symmetry when
r approaches to zero. The boundary conditions (6.26) at the quartz wall and substrate
assume wall recombination occurs when electrons and ions reach the walls and assume
that thermal equilibrium occurs at the walls. The potential is assumed zero on the substrate
and on the quartz walls. This assumption is good for the conducting substrate. However on
the insulating walls a ﬂoating potential is more appropriate. It was observed, though, that
the plasma is conﬁned mostly to the region just above the substrate with only weak
interaction with the quartz disk walls. Hence, the selection of the quartz boundary
condition of ‘1’ = 0 did not alter the simulation results signiﬁcantly.

The grid structure used for the ﬂuid plasma model is shown in Figure 6-6. The
total number of grids used in the plasma simulation is 11 x 11 in the r and z direction
respectively. In the z direction, a uniform grid spacing of 0.57 cm is used. In the r
direction, the grids are constructed in such a way that all the unit cell areas (or volumes,
Aerr21t) are equal. Therefore, the largest grid spacing in the r direction is 2.22 cm, and
the smallest is 0.36 cm. Since the grid structure and grid space for the FDTD model are
different from those of the ﬂuid plasma model, a linear interpolation technique was used to
interpolate the absorbed power to each grid point of the plasma model, and the plasma
density and electron temperature to each grid point of the FDTD model.

For the partially ionized and partially dissociated H2 plasma discharge simulation

developed in this study, the major particle interaction processes are the electron-H2 gas

153

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 6-6: Grid structure for ﬂuid plasma model.

154

inelastic collisions, electron-H2 gas elastic collision and electron-hydrogen ion
recombination. The electron-H2 inelastic collisions include the H2 gas ionization,
excitation and dissociation processes. The rate coefﬁcient for these collision processes are

expressed using the Arrhenius relationship [51] as

k' n = Aioncxp (—€ion/KBTe)

to
kext = Aextcxp (‘Eext/KBTe) (6.27)

kdis = Adiscxp (-8dis/KBTe)

where km, ken, and kdt's are the collision rate coefﬁcient with unit of m3/sec introduced in
Chapter 3; 830", sex, and edis are the threshold energy for H2 gas ionization, excitation and
dissociation; and A ion» Am and Adi: are the pre-exponential factors which are obtained by
approximating the rate constant data at low Te in Janev at a1.[57] to these relationships.
For simplicity, only the reactions with higher rate coefﬁcients are considered in this study
for H2 gas. The types of inelastic collisions and their corresponding rate parameters are
summarized in Table 1. The collision frequency for electron-H2 gas momentum transfer

ven is relatively independent of the electron temperature, so it can be written as [50]

ven (”2) = 1.44 x1012 x Pressugi (Torr)
TnL K)

where Tn is the neutral temperature which can be represented by the translational

 

(6.28)

temperature of H2 gas. Since ven is relatively independent of the electron temperature, the
effective collision frequency Veﬂ‘IS set equal to ven. It should be noted that the only neutral
Species considered in the ion and electron simulations was the H2 species. Hence these

simulations assume a low hydrogen dissociation percentage.

155

Table 1: H2 Reaction Rates

 

 

pre-
exponential
. . Threshold factor or
Reaction Expressron
energy rate
coefﬁcient
(m3/Sec)
Ionization e + H2 _, e + H; + 3 15.4 eV 1.0x10'14
Excitation 12.0 eV 6.5;th:15

*
e+H2—)H2+e

 

. . . -14
Dlssoc1atlon e + H2 _) e + H + H 10.0 eV 1.0x10

 

Recombination 0 eV 1.0x10'l4
e + ion —> neutral

 

 

 

 

 

 

156

6.5 Simulation results for ﬂuid plasma model

The simulation results for the microwave cavity plasma reactor loaded with a H2
discharge are present in this section. Before coupling the FDTD electromagnetic model
with the ﬂuid plasma model to get self-consistent power absorption solutions, the validity
and stability of the ﬂuid plasma model will examined ﬁrst by providing a constant
absorbed power, 500 Watt, as input data. The neutral temperature proﬁle is assumed
neutral as 2000 0K at the discharge center and 1000 0K at the quartz disk boundary. The
2000 0K neutral temperature and the volume of this 2000 0K region was selected based on
experimental measurements of the hydrogen translational temperature and the plasma
volume [54] as shown in Figure 6-7(d). Assuming the pressure is 50 Torr and the input
absorbed power is 500 Watt, the electron density, ion density, and electron temperature are
determined by solving the Poisson equation, the electron/ion continuity equations, and the
electron energy balance equation. The results are shown in Figure 6-7. Part (a) and part (b)
are the electron density and plasma potential. The electron density is higher at the center
and drops gradually toward the boundary to match the boundary condition (ne = ni = 0).
The plasma potential is quite uniform due to the quasi-neutral condition in the bulk plasma
region. Part (c) shows the electron temperature proﬁle. Basically, it is uniform over the
discharge region. The Stability of the numerical solution of the ﬂuid plasma model is
Shown in Figure 6-8. The electron temperature, which is very sensitive to the plasma
density oscillates severely during the beginning iterations since the electron energy
equation and the continuity equations are loosely coupled. By using a damping factor
equal to 10'5 for electron temperature during the beginning iterations, the electron

temperature converges to a stable value after about 30 iterations.

157

 

Tc (eV)

Figure 6-7:

 

Simulation results of microwave cavity plasma reactor using ﬂuid plasma
model with uniform absorbed power density. (a) Electron density (b)

Plasma potential (c) Electron temperature ((1) Neutral temperature
(assigned values).

158

80000 . - . . - - . . -
70000 1
600001

50000.
40000. I 1
30000. H
20000. .
10000. 4*-

 

T, (Kelvin)

 

 

 

 

~10000)
-20000.

-30000 - - - - - - - - -
10 20 30 40 50 60 70 80 90 100

Iteration

 

 

 

 

 

3.7e+187 - r - - - -
3.68e+18
3.66e+18
3.64e+18 .
3.62e+18

3.6e+18
3.58e+l8 , .
3.56e-1-18 t .
3.54e+18 U
3.52e+18 -

. 1 - .. - - - - - - -
”6+8 10 20 30 40 50 60 70 80 90 100

Iteration

 

 

ne (m'3)

 

 

 

 

 

Figure 6-8: Electron temperature and density variations during ﬂuid plasma Simulation.

159

This example shows the application of the ﬂuid model for microwave discharge
Simulation with a constant, uniform absorbed microwave power input. The resulting data
of this simulation can be used as the initial conditions (shown in Figure 6-1) for the self-

consistent microwave cavity plasma reactor Simulation.

6.6 Self-Consistent results for microwave cavity plasma reactor
simulation

The microwave cavity plasma reactor loaded with H2 gas for diamond ﬁlm
deposition is simulated and investigated by a self-consistent numerical model which has
been described in the previous sections. The input parameters for this model includes the
pressure and input microwave power. A set of empirical equations developed by G. King
[58] are used to establish the neutral temperature proﬁle for this Simulation. These
empirical equations are based on the statistically-designed experimental data obtained by
discharge diagnostics in a parameter Space including pressure and input microwave power.
The empirical equations used to predict the translational temperature of H2 gas and the

discharge volume are [58]:

Translational Temperature (0K) = 228.6 + 374.3 x Incident Power (kW) + 16.5 x

Pressure (Torr) i 94.2 (6.29)

Plasma Volume (cm3) = 449.7 + 116.2 x Incident Power (kW) -1 8.1 x Pressure

(Torr) + 57.1 x [Incident Power (kW)]2 + 0.25 X [Pressure (Torr)]2 - 5.4 X Pres-

sure (Torr) X Incident Power (kW) i 15.4 (6.30)

The plasma volume serves as a boundary for neutral temperature (Tn). Inside the plasma

160

volume region, the translational temperature given by (6.29) is the neutral temperature.
Outside the volume, the neutral temperature is equal to the temperature assigned on
boundaries (Tn = 1000 OK). Between these two regions, a linear temperature change
proﬁle is used to prevent an abrupt change of Tn, This is done by assigning the boundary
temperature (Tn = 1000 0K) to the region several grids, namely three grids, away from the
plasma volume in both the r and z directions. The temperature assigned on the grids
between these two region are linearly interpolated from the neutral temperature given by
(6.29) and boundary temperature ("1‘n = 1000 oK).

The simulation results for the electric ﬁeld distribution are shown in Figure 6-9
and Figure 6-10 for 1'5.r and E2 ﬁeld components respectively. The input condition for this
simulation is 50 Torr in pressure and 1500 Watts input microwave power. The electric
ﬁeld patterns basically follow the mom mode electric ﬁelds distribution. The affect of
the electric ﬁeld caused by the presence of discharge can be observed. It reduces the
amplitude of the electric ﬁeld in the plasma discharge region since electromagnetic power
is absorbed in the discharge region. The discontinuity of the electric ﬁeld resulting from
the difference of dielectric constant between the quartz and air can also be seen. The
conductivity of the plasma can be determined from (3.42) using the solved density and
effective collision frequency. Figure 611 shows the maximum conductivity of the
discharge is around 0.4 mho-m. Figure 6-12 shows the power absorption pattern and
indicates that the most microwave power is absorbed in the center region and closed to the
substrate. The electric ﬁeld patterns inside the quartz disk (plasma simulation region) are
shown in Figure 6-13 and Figure 6-14. It can be observed that the Ez ﬁeld doesn’t just

decay into the discharge but also is concentrated in the region above the substrate. This is

161

Quartz Dome

Sliding Short

x104

  

Baseplate

  

  

      

  

  

e
b
O
h

Quartz Tube

      

Substrate

Plasma Discharge

E, ﬁeld distribution.

Figure 6-9

162

 

 

 

m o
l /// 8 )
D. 44%, m ..w
.../x.,////// C T
V// //./I. l/l/l/ (
,, no 4 o o o
B ,,,,..,.|o¢o/oo/o r
o o ‘/ /// 9”” o 6
..u
m 4
m
Z // 2
m 8,4,,
zoo/o.
o, o ,.
a 1.2...
., ooze/084,4
... ‘G/x/xo/o/
«‘6,
5
O
h 0
S 1
g
..m
.1
$
\t/
0 m
2
C
(
4 m. z
0
1 H
X
\ N x > m;
6 4 2 0 0 0
3 C
m m

Substrate

Plasma Discharge

Figure 6-10: Ez ﬁeld distribution.

163

      

         

 

 

   

 

 

/
/ 1.6 ./
z x z 1/ x I/ x a
4,42,42,33

 
    

   
  

  
 

 

 

86:8

          
          

m
. . z 6 9 v” /
O oJIv v'xxI/x/x/
. 5 a .v v! x, xx”
D 00;. x/o/Zte‘xxxxzz
.. [zoo/4444444422,... nil/zoo,” o I
z 444444444? ....I loo/Io
434 2,424? «trivia?
m eeaeeeeeeeeiﬁn-Irea?
u xxx/0848406444,???“ .n§l%oxxozx/xz
././... .. .. .../.2. ... ‘| x x
x ‘ a /
lo/zxz/oh/xx/xoﬁl/zz/xxxzxzx/zérw. @xﬂz“ III/’0 xxx/xx

       

  

Substrate

harge

Plasma Disc

Conductivity distribution.

0

Figure 6-11

164

Quartz Dome

x107

        

   

                    

 

    

      

       

 

            

 
 

   

  

m
C
(
r
4
III/I //
x z I: ’4/4/x
44444442227??? 2,,
4444442222¢ /
xxx/44440.23... /
ZZZ/2,5,3:
o/¢//¢///,¢. I
4x0 xxx/H/II/x/z/z/xz/xzz 1409/0 xxx. I’M/M767 f. N?“ hi‘x/z/x/x/z/x/x/x/I/z/x/x /
z. . .....n.4.4.4///./// \ xxx/xxx”
. . .1./4040.000/4/ \‘xx/I/zz
eel/,o/o/z/oo/o/Ioeo/zioo/o/
xxx/xxx/z/xzx/xxo/x/o/z/xxz/ local/x
o/x/x/zx/xo/zzz/x/zo/z/z. l .4
xxxxzxrzzxzxxxzzzxzxzxxzxzxé.l‘.zzz
oeoeooeooeooooo. .44
3243434342 I‘ll/I
o/z/xI/z/x/x/x/zzxxzzx/z/x/x/x. I‘xxx
oooooozoooo/Ioo/Xfliz
Zeooeeooaoooeoa II o
eel/342442424311]
loo/43434335., II"
oozeoeeoeeéeooé<1ihﬁﬂ
Zoooeeze/Zeaoeeoé .1.”
40444444044444444/454/
[ooze/oo/ooxoo/oox/oo/él
xoéxéxxéxoo/oo/éxééxeé
zooooooo/oo/o/X/o/oo 4
ooezoezoeox/ooo/oozoo/x
zoo/oooo’oozoozoooo/
zooooxoo/oo/oooo/oo
ooo/oo/ooooooo/oo/o
xxxxxzzzx/x/z/z/x/zzxzz/x/z/x/z/ 0
[o/oo/oo/aoxoooo/ 1
oooxoozooooooo/oo
x44/4014/44/44/o
14049404/44/044
ooo/oéxéézoézéz
/////////////I
4.2. 242,442,, 2,,

    

               
 
        

Substrate

  
  

Plasma Discharge

Absorbed power patterns.

Figure 6-12

165

 

Substrate

Figure 6-13: B, ﬁeld pattern inside the quartz disk region.

166

volt/m

 

 

Substrate

Figure 6-14: Ez ﬁeld pattern inside the quartz disk region.

167

because the discharge has a ﬁnite size and the plasma density is not uniform. Therefore,
the electric ﬁeld might propagate into the quartz disk region from other directions and
concentrate at the region above the substrate without being blocked by the discharge.

The simulation result of the characteristics of a H2 discharges at 50 Torr and 1500
Watts input power is shown in Figure 6-15. The neutral temperature proﬁle determined by
the empirical equations is shown at part ((1). The plasma density is concentrated more at
the center. The electron temperature is also highest at the center of the discharge. The
spatial variation of ionization and recombination rate are shown in Figure 6-16. The
ionization rate is higher near the substrate due to the electron temperature being higher at
that region. The recombination rate proﬁle basically follows the electron density variation.
Figure 6-17 shows spatial variations of various energy loss mechanism of electrons,
including ionization, excitation, dissociation, recombination, and elastic collisions. It
Shows the energy loss due to elastic collision is the major energy loss for electrons in H2
discharge under moderate pressure. The convergence behavior of the absorbed power is
Shown in Figure 6-18. When the power absorbed by the plasma converges to a Stable
value, the steady state solution of this numerical model has been obtained.

The characteristics of H2 discharges and cavity quality factor Q are analde for
various pressure and absorbed power by this numerical model. The neutral temperature
and plasma volume determined by (6.29) and (6.30) are the input data for simulations and
their variations for different incident power and pressure are shown Figure 6-19. The
simulation results of average electron temperature vs. absorbed microwave power are
shown in Figure 6-20. These average values are calculated by integrating the temperature

inside the plasma volume and divided by the plasma volume. The average electron

168

(a) Electron Temperature (b) Plasma Potential

CV

 

Figure 6-15: Simulation results of a microwave cavity plasma reactor for pressure = 50
Torr and input microwave power = 1500 watt. (a) Electron density (b)
Plasma potential (c) Electron temperature ((1) Neutral temperature
(assigned values).

169

Ionization

 

 

Recombination

x1021

 

 

=50

Figure 6-16: Spatial variation of ionization rate and recombination rate for pressure
Torr and input microwave power = 1500 watt.

170

Ionization Excitation

x 104 x 105
10

000'!
CON-b

 

6 Dissociation

x10

CON-b

   

x 107 Total Loss

CON-h

   

Figure 6-17: Spatial distribution of various energy loss for pressure = 50 Torr and input
microwave power = 1500 watt. The unit of the energy loss is Watt/m3.

171

 

2900 - . . . . - . .
2800 - r
2700 -
2600 r
2500 .

2400 1

Absorbed Power (Watt)

2300 -

2200 r

2100f ‘\\~\-‘.h .

2000 - - - . 4 - - -
0 20 40 60 80 100 120 140 160 180

 

 

 

Iteration

Figure 6-18: Absorbed power vs. number of iterations.

172

 

 

 

 

 

 

 

 

2000_ *2 40 Torr .
o: 50 Torr X
A x x: 60 Torr
>4
Q, 1800' o -
G
H o
x
1600~ a
x
1 I l l l 1
4(19300 1600 1700 1800 1900 2000 2101
Watt
6? I I I I I
[3, 150- alt: 40 Torr x d
I: o: 50 Torr
E x x: 60 Torr
.E.’
O
:3 100- ~
g o
.52 o x x
a.
m l L l 1 l ‘1
1500 1600 1700 1800 1900 2000 2101
Watt

Figure 6-19: Neutral temperature and plasma volume vs. incident power for various
pressures.

173

 

 

 

 

1.6 r 1 a r .
X
15* *1 40 Torr a
o: 50 Torr

9 1'4_ x. 60 Torr _
<1)
V
G)
g... *
g 1.3- o x ‘
g arr
g 1.2- ~
2
8

1.1 . -
13 o x
2
LL] 0 x

1~ .
0.9- -
1200 1400 1600 1800 2000 2200
Absorbed Power (Watt)

Figure 6-20: Electron temperature vs. absorbed power for three different pressures.

2400

174

temperature is about 1.5 eV. The electron temperature increase with higher absorbed
microwave power. Also, the electron temperature drops with increases of pressure. This
can be understood since higher pressure means higher collision frequency for electron
with neutrals. Therefore, electrons more easily transfer their energy to neutral particles
and reduce their kinetic energy. The results for Q vs. absorbed microwave power are
shown in Figure 6-21 respectively. It shows when the more power absorbed by the
discharge the lower the Q factor. The pressure has only a little effect on the Q factor. At
the 50 Torr pressure and 1500 microwave power condition, the Q factor of the microwave
cavity plasma reactor with a H2 discharge has been experimental determined by using the
technique developed earlier Section 6.3. The Q value determined experimentally is about
100 which is close to the Simulation result (Q = 107). This close agreement is an

additional indication of the validity of the model.

6.7 Summary

The self-consistent simulation results for a microwave cavity plasma reactor were
given in this chapter by coupling the FDTD model and the ﬂuid plasma model developed
in previous chapters. The electric ﬁeld patterns and power absorption patterns in a TMOB
mode microwave cavity plasma reactor operating with a hydrogen discharge were
simulated. The simulated electric ﬁeld distributions showed a good agreement with the
results of electric ﬁeld measurements done using a micro—coax probe technique. The
characteristics of the H2 plasma where also investigated and analyzed for various
microwave powers and pressures. The quality factor Q for the discharge loaded cavity was

determined by numerical simulation and by experimental measurement The two results

175

 

 

 

 

150 I T I
140- at: 40 Torr ‘
130- x x. 50 Torr _
o: 60 Torr
120- -
110- -
X
0’ x o
100- -
O
90- ~
80- -
* X
701- ..
60- a
5 l l L
800 1000 1500 2000

Figure 6-21:

Power absorbed (watt)

Q value vs. absorbed power for three different pressures.

2500

176

showed good agreement which serves as a veriﬁcation of the accuracy of the numerical

model developed in this study.

Chapter 7

Conclusions

A self-consistent numerical model has been developed to simulate the
electromagnetic excitation of discharges inside microwave cavity plasma reactors. This
software includes a electromagnetic ﬁeld model and a ﬂuid plasma model. The microwave
cavity plasma reactors simulated by this numerical software were loaded with a H2
discharge and used for diamond thin ﬁlm deposition. The complex reactor geometry and
input power coupling probe Structure were included in the simulation. The input
parameters for the microwave cavity plasma reactor simulation included the pressure and
input microwave power. The operating pressure investigated in this study was the
moderate pressure range (1 Torr to 100 Torr). The input microwave power was from
several hundred watts to several kilowatts. The electromagnetic behavior of the discharge
loaded resonant cavity, such as the electric ﬁeld distributions, power absorption patterns,
and cavity quality factor Q, were studied and analyzed for various input condition. The
characteristics of H2 discharges, including the plasma density, electron temperature, and
plasma potential, were also investigated and studied. The electromagnetic mode and
cavity Q value of a H2 discharge loaded reactor were simulated and shown to be in

agreement with measured experimental values.

177

178

The electromagnetic ﬁeld model developed in this research is based on using the
ﬁnite-difference time-domain (FDTD) numerical method to solve the time-dependent
Maxwell’s equations. The major advantage of this method is that it can be easily
implemented to solve the ﬁeld in complex geometry Structures with complex loaded
materials, such as plasmas. In order to investigate the microwave power absorbed by the
lossy materials or discharges, appropriate conductivity models were identiﬁed and used to
solve the microwave ﬁeld induced current density. For plasma discharges, this induced
current density was determined by solving the electron momentum transport equation.

Chapter 4 presented the transient simulation results of microwave resonant cavities
using the three-dimensional FDTD simulation method for empty and loaded cylindrical
cavities. Comparisons to exact analytical models using simple geometry cavities were
used to verify the accuracy of the electromagnetic ﬁeld model. The results of empty cavity
simulations agreed well with the theoretical resonant mode. The simulated natural
frequency of a lossy loaded cavity also match the analytical solutions. A compact ECR ion
source and a Ar loaded cavity were also simulated by using a magnetized discharge
conductivity model for the ECR source and by simultaneously solving the electron
momentum transport equation in the time domain, respectively The steady state
simulation method was introduced by including the input power coupling structure into
the electromagnetic model as shown in Chapter 5. By implementing a non-reﬂecting
boundary condition, the open boundary problem at the end of the input coupling probe
was solved. The Q value and the power reﬂected from the microwave cavity plasma
reactor were analyzed for various cavity height and loaded conductivity. This simulated

reﬂected power characteristic behavior is similar to that found in experimental systems.

179

A ﬂuid plasma model was developed in Chapter 3. It solved the electron and ion
continuity equations, the electron energy balance equation, and the Poisson equation. The
input to this plasma model was the microwave power density absorbed by the discharge,
which was calculated by the electromagnetic ﬁeld model. The discharge characteristics
simulated by the plasma model were feedback to the electromagnetic ﬁeld model to
update the electromagnetic ﬁeld and power absorption. The ﬁnal self-consistent results
were solved iteratively to a converging solution by these two models. This self-consistent
solution was done using two dimensional solution of the FDTD and ﬂuid plasma models.

The electromagnetic ﬁeld solution technique present in this study allows the
calculation of the power absorption proﬁle in microwave discharges which has often been
a difﬁculty in plasma simulations. The understanding obtained from the electromagnetic
solution allows the reactor design structure to be analyzed for improvement and
optimization of such quantities as uniformity of microwave power absorption.
Additionally, the simulation and understanding of the electromagnetic ﬁelds is expected to
be key to microwave reactor control.

The techniques developed in this study could in the future be applied to a wide
range of plasma sources including the various high density sources, e.g. inductively
coupled sources, helicon source and ECR sources. Other future work extending this study
includes self-consistently solving the temperatures and densities of various neutral
species, expanding the two dimensional ﬂuid plasma model to three dimensions, and
developing a numerical model to simulate the side-feed input coupling structure for TED”

mode cavities.

APPENDIX

APPENDIX

The derivation of discretized Poisson equation (3.37) is provided here in details.

The Poisson equation can be written as
VOE = E , (A-l)

In 2-D cylindrical coordinates, it can be expanded to

BE

Z

a _ 12
-a-;(rEr) +ra—z — r8. (A-Z)

Integrating from r,-_1/2 to ri+ 1,2 and Zk- 1/2 to 2 “1,2, (A.1) becomes

. 1 . 1
’i+1/2Er(‘+§rk)° (Zk+l/2‘Zk—1/2)"ri—l/zEr('-§rk)' (Zk+1/2'Zk—l/2)

2 2
r0 +ro
I ””2 "1”](521‘3“ilEzi‘tk-ill -<A-3>

2 2
_ p(i,k) ri+l/2+ri-l/2 ( )
" e 2 Zk+l/2—Zk—l/2

 

 

Since

180

181

Azk+AZk—l
Zk+l/2—zk-l/2 = ‘—_2—— , (A4)

for uniform grid spacing in the z direction and
P (13k) = e(n,-(i. k) —n,.(i. k)) . (A5)

(A3) can be rewritten as

 

2 . 1 .1
( 2 2 )'(ri+1/2Er('+§rk)'ri—l/2Er(”§rk))

’1+l/2“’i—l/2

((M1)

= Ema, k) —n,(i. k))

 

Moreover, using

V‘i’ = E, (A.7)

E, and E2 are discretized to be

 

 

. 1 _ ‘I’(i+1,k)-‘I’(i,k)
Er(r+-2-,k) — Ar(i) (A8)
and

. _l_ _ ll’(i,k+1)—‘1’(i,k)
Ez(1,k+2) — Az (k) . (A.9)

Substitute (A8) and (A9) into (A.6) the discretized Poisson equation is

182

 

 

 

r.
42 22 ‘-l+l/2°\P(I+I,k)
Ar.
(rt-l-l/Z—ri—l/Z) ‘
r.
+- 2 2 2 --"“2.\11(i—1,k)
Ar.
(ri+1/2_ri-l/2) “1
+ 2 .‘i’(i,k+l)
(Azk+Azk_,) A2,,
2 .‘I’(i,k—1) (A.10)

 

+ (Azk+Az,,_1) Az,_,

_[, 2 [ .[ri+l/2+ri-l/2]
2 2 Ar. Ar.
(’i+1/2"’i-1/2) ‘ "l
2 l l .
+ . —+ ~‘l’ r,k
(Azk+AZk-l) (All. AZk—l)) ( )

= —g (n,(i, k) -ne(i. k))

 

 

 

 

as shown in (3.37). The other discretized equations, such as electron continuity equation
(3.38) and electron energy balance equation (3.41), can be derived in a Similar manner.
In case of the point r = 0 (i = 1), then (A2) is integrated from 0 to rm instead of

from 0-1/2 to 'i+1/2- With (A.4) and (A5), (A.6) becomes

 

3 2 1 l
’3/2Er15’*)+ (Azk+Azk_,) '(Ezlltkwl‘Ezilrk‘ill

= gln,(1.k)-n,(l.k))

(A.ll)

Then, substituting (A8) and (A9) into (A. 11), the discretized Poisson equation for r = 0 is

183

2 l 2 1
—-—-‘P(2,k)+—-—-‘l’(l,k)
'3/2 “I ’3/2 ml
2 .‘I’(l,k+1)+ 2 _‘I’(l,k—1)

 

+ (AZk+AZk-1) AZ]: (Azk+AZk-l) AZk—l

2 1 l
— — — -‘l’ l,k
(AZk+AZk_1)(Azk+Azk_1) ( )

(A.12)

 

= —g(ni(l,k)—ne(1,k))-

Similarly, the discretized electron continuity equation and electron energy balance

equation at r = O can be derived as

 

2 (3 l 2 ( ( ‘) ( 1)) '
—-J -,k + - Je 1,k+— -J¢ l,k—-
r3/2 er 2 (AZk+AZk-l) Z 2 Z 2 (A.13)

= ne(i, k) [nn(1,k)kion(l,k)—ar(l,k)ni(l,k)]

and

._2_. (a) 2 H 1)-( A»
qer 22k +(A2k+AZk-1) qez12k+2 qez l,k 2

’3/2

 

(A.14)
(1, k)

abs

=ne(l,k)2Hi(l,k)+P
,-

where Je, Jez, qe, and qeZ are deﬁned in (3.39), (3.40), (3.42), and (3.43).

LIST OF REFERENCES

[1]

[2]

[3]

[4]

[5]

[6]

[71

[8]

[9]

LIST OF REFERENCES

J. Asmussen, “Electron cyclotron resonance microwave discharges for etching and
thin-ﬁlm deposition,” J. Vac. Sci. Technol. A, vol. 7, no. 3, pp. 883-893, May/June
1989.

S. Whitechair, J. Assmussen, and S. Nakanishi, “Microwave electrothermal

thruster performance in helium gas,” J. Propul. Power, vol. 3, pp. 136-144, March]
April 1987.

J. Zhang, B. Huang, D. K. Reinhard, and J. Asmussen, “An investigation of

electromagnetic ﬁeld patterns during microwave plasma diamond thin ﬁlm
deposition,” J. Vac. Sci. Technol. A, vol. 8, no. 3, pp. 2124-2128, May/Junel990.

J. Hopwood, M. Dahimene, D. K. Reinhard, and J. Asmussen, “Plasma etching
with a microwave cavity plasma disk source,” J. Vac. Sci. Technol. B, vol. 6, no. 1,
pp268-27l, Jan/Feb 1988.

S. J. Pearton, U. K. Charkabarti, A. Katz, A. P. Perley, and W. S. Hobson,
“Comparison of CI-I4/1-12/Ar reactive ion etching and electron cyclotron resonance
plasma etching of In-based III-V alloys,” J. Vac. Sci. Technol. B, vol. 9, no. 3, pp.
1421-1432, May/June 1991.

J. R. Rogers, “Properties of steady-state, high-pressure, argon microwave
discharges,” Ph.D. Dissertation, Michigan State University, East Lansing, 1982.

M. L. Passow, M. L. Brake, P. Lopez, W. B. McColl, and T. E. Repetti,
“Microwave resonant-cavity-produced air discharges,” IEEE Trans. Plasma Sci,
vol. 19, no. 2, pp. 219-228, April 1991.

J. B. Manring, “Electromagnetic ﬁeld solutions for the natural modes of a
cylindrical cavity loaded with lossy materials,” Ph.D. Dissertation, Michigan State
University, East Lansing, 1992.

S. Offermanns, “Resonance characteristics of a cavity-operated electrodeless high-

184

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

185

pressure microwave discharge system,” IEEE Trans. Microwave Theory Tech., vol.
38, no. 7, pp. 904-911, July 1990.

S. Offermanns, “Electrodeless high-pressure microwave discharges,” J. App.
Phys, vol. 67, no. 1, pp. 115-123, Jan 1990.

K. S. Yee, “Numerical solution of initial boundary value problems involving
Maxwell’s equations in Isotropic media,” IEEE Trans. Antennas Propagat., vol.
14, no. 3, pp. 302-307, May 1966.

R. Holand, “THREDE: a free-ﬁeld EMP coupling and scattering code,” IEEE
Trans. Nuclear Sci., vol. 24, no. 6, pp. 2416-2421, Dec 1977.

A. Taﬂove, and M. E. Brodwin, “Numerical solution of steady-state
electromagnetic scattering problems using the time-dependent Maxwell’s
equations”, IEEE Trans. Microwave Theory Tech., vol. 23, no. 8, pp. 623-630, Aug
1975.

A. Taﬂove, “Application of the ﬁnite-difference time-domain method to sinusoidal
steady-state electromagnetic-penetration problems,” IEEE Trans. Electromagn.
Compat., vol. 22, no. 3, pp. 191—202, Aug 1980.

G. Mur, “Absorbing boundary conditions for the ﬁnite-difference approximation
of the time-domain electromagnetic-ﬁeld equations,” IEEE Trans. Electromagn.
Compat., vol. 23, no. 4, pp. 377-382, Nov 1981.

R. Luebbers, F. P. Hunsberger, K. S. Kunz, R. B. Standler, and M. Schneider, “A
frequency-dependent ﬁnite-difference time—domain formulation for dispersive
materials,” IEEE Trans. Electromagn. Compat, vol. 32, no. 3, pp. 222-227, Aug
1990.

C. Wang and O. P. Gandhi, “Numerical simulation of annular phased arrays for
anatomically based models using the FDTD method,” IEEE Trans. Microwave
Theory Tech., vol. 37, no. 1, pp. 118-126, Jan 1989.

D. Sullivan, “Three-dimensional computer simulation in deep regional
hyperthermia using the ﬁnite-difference ' time-domain method,” IEEE Trans.
Microwave Theory Tech., vol. 38, no. 2, pp. 204-211, Feb 1990.

J. Chen, and O. P. Gandhi, “Electromagnetic deposition in an anatomically based
model of man for leakage ﬁelds of a parallel-plate dielectric heater,” IEEE Trans.
Microwave Theory Tech., vol. 37, no. 1, pp. 174-180, Jan 1989.

X. Zhang and K. K. Mei, “Time-domain ﬁnite difference approach to the
calculation of the frequency-dependent characteristics of microstrip discontinuity,”
IEEE Trans. Microwave Theory Tech., vol. 36, no. 12, pp. 1775-1787, Dec 1988.

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

186

X. Zhang, J. Fang, K. K. Mei, and Y. Liu, “Calculations of the dispersive
characteristics of microstrips by the time-domain ﬁnite difference method,” IEEE
Trans. Microwave Theory Tech., vol. 36, no. 2, pp. 263-267, Feb 1988.

G. Liang, Y Liu, and K. K. Mei, “Full-wave analysis of coplanar waveguide and
slotline using the time-domain ﬁnite-difference method,” IEEE Trans. Microwave
Theory Tech., vol. 37, no. 12, pp. 1949-1957, Dec 1989.

D. M. Sheen, A. M. Ali, M. D. Abouzahra, and J. A. Kong, “Application of the
three-dimensional ﬁnite-difference time domain method to the analysis of planar

microstrip circuits,” IEEE Trans. Microwave Theory Tech., vol. 38, no. 7, pp. 849-
856, July 1990

D. H. Choi, and W. J. R. Hoefer, “The ﬁnite-difference-time-domain method and
its application to eigenvalue problems,” IEEE Trans. Microwave Theory Tech.,
vol. 34, no. 12, pp. 1464-1469, Dec 1986.

A. Navarro, M. J. Nunez, and E. Martin, “Study of TEO and TMO modes in
dielectric resonators by a ﬁnite difference time-domain method coupled with the

discrete fourier transform,” IEEE Trans. Microwave Theory Tech., vol. 39, no. 1,
pp. 14-17,Jan 1991.

M. F. Iskander, O. Andrade, H. Kimrey, R. Smith, and S. Lamoreaux,
“Computational techniques in modeling and quantifying microwave interactions
with materials,” Ceramic Transactions, vol. 21, Ed. Clark, Gae, and Sutton, pp.
141-147, 1991.

D. B. Graves and K. F. Jensen, “A continuum model of DC and RF discharges,”
IEEE Trans. Plasma Sci., vol. 14, no. 2, pp. 78-91, April 1986.

D. B. Graves, “Fluid model simulations of a 13.56-MHz rf discharge: time and

space dependence of rates of electron impact excitation,” J. Appl. Phys., vol. 62,
no. 1, pp. 88-94, July 1987.

J. P. Boeuf, “Numerical model of rf glow discharges,” Physical Review A, vol. 36,
no. 6, pp. 2782-2792, Sep 1987.

J. H. Tsai and C. Wu, “Two-dimensional simulations of rf glow discharges in N2
and SF6,” Physical Review A, vol. 41, no. 10, pp. 5626-5644, May 1990.

M. J. Kushner, “Monte-Carlo simulation of electron properties in rf parallel plate
capacitively coupled discharges,” J. Appl. Phys., vol. 54, no. 9, pp. 4958-4956, Sep
1983.

B. E. Thompson and H. H. Sawin, “Monte Carlo simulation of ion transport
through rf glow-discharge sheaths,” J. Appl. Phys., vol. 63, no. 7, pp. 2241-2251,

[33]

[34]

[35]

[36]

[37]

[38]

[391

[40]

[41]

[42]

[43]

[44]

187
April 1988.

J. I. Ulacia and J. P. McVittie, “A two-dimensional computer simulation for dry
etching using Monte Carlo techniques,” J. Appl. Phys., vol. 65, no. 4, pp. 1484-
1491, Feb 1989.

G. R. Govinda Raju and M. S. Dincer, “Monte Carlo simulation of electron
swarms in nitrogen in uniform E x B ﬁeld,” IEEE Trans. Plasma Sci., vol. 18, no.
5, Pp. 819-825, Oct 1990.

T. E. Sheridan, M. J. Goeckner, and J. Goree, “Model of energetic electron
transport in magnetron discharges,” Vac. Sci. Technol., vol. A8, no. 1, pp. 30-37,

Jan/Feb 1990.

C. K. Birdsall, “Particle-in-cell charged-particle simulations, plus Monte Carlo
collisions with neutral atoms, PIC-MCC,” IEEE Trans. Plasma Sci., vol. 19, no. 2,
pp. 65-85, April 1991.

D. Vender and R. W. Boswell, “Numerical modeling of low-pressure RF plasmas,”
IEEE Trans. Plasma Sci., vol. 18, no. 4, pp. 725-732, April 1990.

M. Surendra and D. B. Graves, “Particle simulations of radio-frequency glow
discharges,” IEEE Trans. Plasma Sci., vol. 19, no. 2, pp.144-157, May 1991.

H. W. Trombley, F. L. Terry, and M. E. Elta, “A self-consistent particle model for
the simulation of RF glow discharges,” IEEE Trans. Plasma Sci., vol. 19, no. 2, pp.
158-162, April 1991.

R. W. Boswell and D. Vender, “Simulation of Pulsed electropositive and
electronegtive plasmas,” IEEE Trans. Plasma Sci., vol. 19, no. 2, pp. 141-143,
April 1991.

T. A. Grotjohn, “Numerical modeling of a compact ECR ion source,” Rev. Sci.
Instruments, vol. 63, pp. 2535-2539, 1992.

M. S. Barnes, T. J. Colter, and M. E. Elta, “Large-signal time-domain modeling of
low-pressure rf glow discharges,” J. Appl. Phys., vol. 61, no. 1, pp. 81-89, Jan
1987.

M. S. Barnes, T. J. Cotler, and M. E. Elta, “A staggered-mesh ﬁnite-difference

numerical method for solving the transport equations in low pressure RF glow
discharges,” J. Comp. Phys, vol. 77, pp. 53-72, 1988.

N. Sato and H. Tagashia, “A hybrid Monte Carlo/ﬂuid model of RF plasmas in a
SiH4/I-12 mixture,” IEEE Trans. Plasma Sci., vol. 19, no. 2, pp. 102-112, April
1991.

[45]

[46]

[47]

[48]

[49]

[50]

[51]

[52]

[53]

[54]

[55]

[56]

[57]

[58]

188

R. K. Porteous and D. B. Graves, “Modeling and simulation of magnetically
conﬁned low-pressure plasmas in two dimensions,” IEEE Trans. Plasma Sci., vol.
19, no. 2, pp. 204-213, April 1991.

Y. Weng and M. J. Kushner, “electron energy distribution in electron cyclotron
resonance discharges for materials processing,” J. Appl. Phys., vol. 72 no. 1, pp.
33-42, July 1992.

R. W. Boswell and I. J. Morey, “Self-consistent simulation of a parallel-plate rf
discharge,” Appl. Phys. Lett, vol. 52, no. 1, pp. 21-23, Jan 1988.

J. Zhang, “Experimental development of microwave cavity plasma reactors for
large area and high rate diamond ﬁlm deposition,” Ph. D. Dissertation, Michgan
State University, 1993.

R. F. Harrington, Time-Harmonic Electromagnetic Fields, McGraw-Hill, 1961.

B. E. Cherrington, Gaseous Electronic and Gas Lasers, Pergamon, New York,
1979.

M. Surendra et a1, “Self-consistent dc glow-discharge simulations applied to
diamond ﬁlm deposition reactors,” J.Appl. Phys., vol. 71, no. 10, pp. 5189-5198,
May 1992.

J. D. P. Passchier and W. J. Goeheer, “A two-dimensional ﬂuid model for an argon
rf discharge,” J.Appl. Phys., vol. 74, no. 6, pp. 3744-3751, Sep 1993.

S. Ramo, J. R. Whinnery, and T. Van Duzer, Fields and Waves in Communication
Electronics, Wiley, New York, 1965.

G. L. King and T. A. Grotjohn, 40th National Symposium of the American
Vacuum Society, Orlando, 1993.

W. H. Press et a1, Numerical Recipes, Cambridge University Press, 1986.

T. A. Grotjohn et a1, “Modeling the electromagnetic excitation of a compact ECR
ion/free radical source,” Rev. Sci. Instruments, vol. 65, no. 5, pp. 1761-1765, May
1994.

R. K. Janev et a1, Elementary Processes in Hydrogen-Helium Plasmas, Springer,
Berlin, 1987.

G. King, “Temperature and concentration of ionic and neutral species in resonant
microwave cavity plasma discharges,” Ph.D. Dissertation, Michigan State
University, 1994.

189

[59] Fongray F. Young and Chwan-Hwa Wu, “Two-dimensional, self-consistent, three-
moment simulation of RF glow discharges,” IEEE Trans. Plasma Sci., vol. 21, no.
3, pp. 312-321, June 1993.

[60] Yao-Tsung Tsai, “AC simulation of ﬁeld effect transistors with a hydrodynamic
transport model,” Ph.D. Dissertation, Michigan State University, 1990.

[61] Chi-Jung Huang, “An investigation of hot electron induced degradation for silicon
bipolar transistors,” Ph.D. Dissertation, Michigan State University, 1992.

[62] Dennis M. Manos and Daniel L. Flamm, Plasma Etching, Academic press, San
Diago, 1989.

[63] R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles,
McGraw-Hill, New York, 1981.

[64] C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simulation,
McGraw-Hill, New York, 1985.

[65] M. E. Tajima-Tochiki, Computational Plasma Physics, Addison-Wesley, 1988.