Modeling of moderate pressure H2/CH4 microwave
discharge used for diamond deposition
K. Hassouni, F. Silva, G. Haagelar, X. Duten, G. Lombardi, A. Gicquel(1)
T. A. Grotjohn(2)
M. Capitelli(3)
J. Röpcke(4)
(1)LIMHP,
CNRS-UPR 1311, Univeristé Paris 13, 99 Avenue J. B. Clément
93430 Villetaneuse, France
(2)Department
of Electrial Engineering, Michigan State University, USA East Lansing, USA
(3) Dipartimento
(4)
di Chimica – IMIP CNR Univeristy of Bari Italy
INP Greifswald
Investigated Device : microwave cavity
coupling system + belljar quartz vessel
Usual Experimental Conditions :
Feed Gas : H2/CH4 - %CH4 < 5
Pressure : 20-200 Torr
Input Power : 0.4-4 kW
Power density : 6-100 W/cm3
Flow Rates : 100 - 300 sccm
Substrate Temperature : 1000 K - 1300 K
Orders of magnitude :
Plasma height 3.5 cm (for 2.45 GHz)
Tg > 2000 K - 1011 cm-3 < ne <1013 cm-3
Transport and collisional phenomena in the plasma
Wave-plasma interaction
Electron heating
(E,H)
e
e + AB
AB(v), AB(r), AB
AB+, A, B
Electron/heavy species
collisions
Energy transfer, ionization, etc...
Heavy species-heavy
species collisions
Intermode energy transfer, chemistry
Convection, Diffusion
Energy and mass
transfer
EEDF, <ee>
Energy and mass transport
Plasma-surface interaction
Mass and energy transfer
VDF, Tv
Tg
ni, i=1-n
Substrate (Ts, cs)
What information do we need ?
Strongly linked
I. For a given coupling configuration, What’s the optimal reactor configuration
that insures :
1- enhanced density of active species at the growing diamond film,
2- a good thermal stability of the reactor (during up to several weeks )
process engineering :
mass and energy transport, gas phase and surface chemistry
II. What is the optimal electromagnetic coupling conditions in term of :
1- power
2- pressure
2- frequency
Electrical engineering :
resonant electromagnetic modes, wave-plasma interaction
Modeling approach
Detailed collisional model in moderate pressure H2 Plasmas
* Determination of the main chemical processes and physical phenomena in the reactor.
* Set up a satisfactory and useful Physical plasma model
2D Self-consistent model (Only for H2 discharges)
 Determine self consistent plasma, electomagnetic field and absorbed power distributions.
for axisymmetrical configuration
 May be used for power deposition optimisation and for scaleup
1D transport model of the plasma on the stagnation line
for both H2 and H2/CH4 discharges
* Investigate the coupling between Chemistry, Energy Transfer and Transport Phenomena
* Determine the behavior of the plasma temperatures and species.
Detailed Collisional model – H2 (1/2)
I. Vibrational Kinetics
Larichuta, Celiberto et al.
* e- + H2(X,v) ==> e- + H2(X,w) |w-v| < 4
* e- + H2(X,v) ==> e- + H2(B1Su,v’) ===> e- + H2(X,w)
* vv relaxation : H2(v) + H2(w) ===> H2(v1) + H2(w  1)
* v-t exchange : H2(v) + H2 (or H) ===> H2(v1) + H2 (or H)
* Dissociation from upper v by heavy species collisions:
H2(v) + H2 (or H) ===> 2H + H2 (or H)
II. Ground States Kinetics
C. D. Scott et al. J. Thermophysics, Vol. 10 1996, p 426
* e- + H2(v=4-10) ==> H + H* H2+ + H2 ==> H3+ + H
* Recombination and mutual neutralization of ions
For more detailed description : K. Hassouni, A. Gicquel, M. Capitelli and J. Loureiro
Plasma source Science and Technology, 8(3), 494 (1999)
.../..
Detailed Collisional model
III. Electronically Excited States Kinetics : H2 and H
Larichuta, Celiberto et al.
H2 : K. Sawada et al. J. Appl. Phys., 78 (5), (1995), p. 2913
H : J. A. Kunc et al., Phys. Fluids, Vol. 30(7), (1987), p. 2255
* e- + H2(v) ==> e- + H2* and
e- + H(n) ==> e- + H(m)
* e- + H2 (v) = [H2**]=> H2+ + 2e- and
e- + H(n) ==> 2e- + H+
* e- + H2 (v) = [H2**]=> 2H + e* H2* + M ==> H2* or 2H
* H(n) + H ==> H(m) + H
* M* ==> M*’ + hn (M=H or H2)
Optically thick plasma for Lyman radiations
M. Glass-Maujean Phys. Rev. Lett., Vol. 62 (2), (1989), p. 144
H(n=2) + H2 ==> H3+ + e- (s = 15-30 A2)
H(n>2) + H2 ==> H3+ + e- - Assumption in this work For more detailed description : K. Hassouni, A. Gicquel, M. Capitelli ad J. Loureiro
Plasma source Science and Technology, 8(3), 494 (1999)
Simulation procedure
Homogeneous
Plasma
Use the effective field assumption
=> stationary situation
EEDF
VDF
H2*(n=1-36)
H*(n=1-40)
H+, H2+, H3+, H-, eTg
D ata
M W P i n p , P , dt , ds, w/ N
I n i ti a l G u e s s
E / N , P l a s m a c o m p o si ti o n , T g
Correct E/N according to (MWPc - MWPinp)
Thin
Boundary Layer
Substrate
Use of the detailed kinetics in the frame
of Quasi-homogeneous plasma model
E l e c t r o n B o l t z m a n n E q u a ti o n
T w o ter m e x p a nsio n
E E D F , R a t e C o n s t a n t s , n e -s
C a l c u l a te t h e t r a n s p o r t C o e f f i c i e n t s
S p e ci e s a n d e n e r g y L o s s e s a t t h e w a ll
S p e ci e s k i n e ti c s E q u a ti o n s
N e w P l a s m a c o m p o s i ti o n
T o t a l E n e r g y E q u a ti o n
N e w Tg
C alc u la te th e n e w a bs o rb e d M W p o w e r
No
| M W P c - M W Pin p| < e
Y es
T g , E / N , E E D F , V D F a n d S p e c i e s d e n s i ti e s
Quasi-homogenous plasma model
Most significant results : vibrational distribution/
MWPDav = 4.5 W/cm
-1
10
MWPDav = 9.0 W/cm
MWPDav = 15.0 W/cm
-2
10
MWPDav = 22.0 W/cm
-3
10
VDF - (Nv/N0)
3
MWPDav = 30.0 W/cm
-4
10
-5
10
-6
10
-7
10
-8
10
Vibrational quantum
13
-9
10
0
1
2
3
4
5
6
7
8 9 10 11 12 14
-10
10
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5
Vibrational energy (eV)
Densité de puissance dissipée (W/cm3)
0
10
102
Exc. Vib.
101
Exc. trans-rot.
100
Dissociation
Exc. Elec. H2
10-1
Exc. Elec. H
10-2
10-30
5
10 15 20 25
MWPDav (W/cm3)
30
Quasi-homogenous plasma model : most significant results
EEDF behavior and electron impact rate constants
rate constants of Electron-impact process
H2 dissociation
Bimodal distribution
1E+00
1E-02
Te-l = 19200 K
FDEE (eV-3/2)
Kd (cm^3/s)
1E-04
1E-06
Te-h = 8500 K
1E-08
1E-10
1E-9
1E-10
1E-11
1E-12
1E-13
1E-14
1E-15
1E-16
1E-17
1E-18
1E-19
1E-20
1E-21
1E-22
[H]=0
[H]=0.1
[H]=0.2
[H]=0.3
0
1E-12
0
0,5
1
1,5
2
2,5
3
3,5
Electron Average Energy (eV)
5
10
15
20
25
ee(eV)
FDEE Bimodale
Te-h = f(Te-l) (univoque)
Energy balance  Te-l  Te-h
Rate constants only depend on Te-l
 We can use a scalar model to describe the electron kinetics
4
Quasi-homogenous plasma model
Most significant results
108
100
2
3
106
10-1
10-2
H(n=2-3) + H2 => H3+ + e-
Densité(cm-3)
ionisation rate (mol. m-3.s-1)
Ionization kinetics & H-atom excited states kinetics
e-+H2 => 2e- + H2+
H2+ + H2 => H3+ + e-
10-3
e-+H => 2e- + H+
10-4
10
15
20
25
30
MWPDav (W/cm3)
Main ionization channel : quenching of
H(n=3) and H(n=2) states
35
4
5
6
104
102
100
9.0 W/cm3
10-2
30.0 W/cm3
20
10-4
10
11
12
En (eV)
13
Excited state distribution
coupling between excited levels H*
H(n=2-3) kinetics do not depend on
H(n>3) states
14
Quasi-homogenous plasma model
Most significant results
From full model to simplified model :
9 species [H2,H,H(n=2),H(n=3),H+,H2+,H3+, e-] /2 energy modes [Tg, eedf]
12
1x10
Simplified model
0.18
Simplified model
11
9x10
0.16
11
Electron density (cm
)
0.14
-3
H-atom Relative density
0.2
0.12
0.1
0.08
0.06
0.04
8x10
11
7x10
11
6x10
11
5x10
11
4x10
11
3x10
0.02
Detailed model
11
Detailed model
2x10
0
0
10
20
30
Power density (W/cm^3)
40
11
1x10
0
5
10 15 20 25 30 35
3
MWPDav (W/cm )
2D Self-consistent model for
moderate pressure plasma flow
2 Modules
EM field simulation
domain
A plasma module
simulates the discharge in the low
pressure vessel
An electromagnetic module
simulates the EM field in the whole cavity
Hassouni, Grotjohn and Gicquel, J. Appl. Phys. (1999)
Hagelaar, Hassouni and Gicquel, J. Appl. Phys. (2004)
Plasma
simulation
domain
Plasma module
2 or 3-temperature thermo-chemically non-equilibrium flow models
ni
+   (- Di ni
t
=>
Continuity => species density
u ni )- W = 0
i
=>
+
Navie-Stockes => Flow velocity
Ee
+ div(uehe-leTe) –PEM + Qe-t + Qe-v + Qe-C =0
t
Electron energy => Te
EM module
H2-vibration mode => Tv
Total energy => Tg
E
t
Ev-m
+ div(umEv-m-lv-mTv-m) + Qv-t - Qe-v + Qv-C = 0
t
+ div[Suihi –lt-rT -Slv-mTv-m –leTe] - PEM + Qrad = 0
EM module
Coupling nodes
Electromagnetic module
A time-domain model where the plasma is considered as a high
frequency conductor
Maxwell Curl equations
 E =
H
t
  H = JHF + e0
E
t
J HF = -qe ne vHF
plasma model
Momentum equation for the HF component of the electron flow velocity
me
dvHF = - - u
qE me eff vHF
dt
plasma model
Hassouni, Grotjohn and Gicquel, J. Appl. Phys. (1999)
Hagelaar, Hassouni and Gicquel, J. Appl. Phys. (2004)
2D Self-consistent model for moderate
pressure plasma flow
Numerical method : Finite Difference Time Domain Method
Er(i+1/2,j+1,k)
 Explicit time integration
H
t
rot ( H ) = e 0
E
+J
t
Are naturally satisfied at the grid points
Ez(i+1,j,k+1/2)
rot ( E ) = - 0
Hr(i,j+1/2,k+1/2)
Hz(i+1,j+1/2,k)
Staggered grids for the different
components of E and H
Ez(i,j,k+1/2)
Leap-Frog shift between E and H, to acheive second order
accuracy
Er(i+1/2,j,k)
Only TM modes are considered in this calculation
 2 kinds of boundary conditions are used :
Perfect conductor: Et=0 and Hp=0
Nonreflecting boundary condition at the cavity inlet
 Define an excitation plane in
the computation domain
2D Self-consistent model
Iteration scheme
Simulation procedure
Te, Tg, ne
Plasma model
12 Transport equations
9 espèces -Tg, Tv, <ee>
300-400 iterations
Grid interpolation
wp, nQ-M  ep
15-20
Iterations
Grid interpolation
Electromagnetic module
Maxwell Curl equations
HF Momentum equation for electrons
E, H, Je-HF, MWPD
15-25 microwave periods
Microwave power density
E.Je-HF
Validation of the electromagnetic module
Electric field
intensity in the
cavity
Self consistent model
without plasma
CST Microwave Studio
(Commercial code)
Some results : 600 W – 25 mbar
MWPDav=8W/cm-3
Power density (W.cm-3)
Tg(K)
2000 K
18 W/cm3
ne(K)
4x1011 cm-3
Te(K)
20000 K
Optimal power deposition at 50 mbar
||E||
1 ball regime
ne
nH
Tg
MWPinp↗ Vplasma ↗
then
MWPinp Transition 1B to 2B
500 W
600 K
Gas heating is responsible for
Discharge regime transition
Ignition for maximum E/N
Tg ↗  E/N distribution changes
800 W
1200 K
Very similar to some explenation
given for streamer to arc transition
phenomenon
1000 W
2D-self consistent model
Flow effect at high power density
With free convection
Strong free convection
Without free convection
Inlet
Radial uniformity
Simulation Domain
1D transport model for H2 plasma flow
•2 Momentum Equations => V=dv/dr|r=0 and u
dV
d  dV
r u dz + rV2 - dz  dz +L =0


et
du
u dr
+
2V
+
dz
r dz =0
•continuity equations for species :
d Ms  d  Ms  qk dTg dxs
ru dz  M xs + dz r M Ds  T dz - dz  - Ws =0


g
Substrate
• 2 energy equations : Tg, Te-l
2D SC model
2 d<ee> + d  2 d<ee> 
re cp-e (u + ue) 3k dz
dz - 3k le dz 
MWPD + Qe-v + Qe-t + Qe-C
=0
1D transport model for H2/CH4 plasma flow
31 species - 134 reactions model
Three reactions Groups
H2/CH4 (Thermal Hydrocraking : C1-2H1-6 species) :
* B. W. Yu and S. L. Girshick, J. Appl. Phys., 75(8), 1994, p. 3914
* C. T. Bowman et al., http://www.me.berkley.edu/gri_mech/
Pressure correction on three body recombination reaction
CxHy Charged species kinetics :
Ionization kinetics : e- + CxHy ==> 2e- + CxHy+
•H.Tawara et al., Research Report NIFS-DATA,
•Charge transfer kinetics : H3+ + CxHy ==> H2 + CxHy+1+
* H. Tahara et al., Jpn. J. Appl. Phys., 34.
* Dissociative Recombination : CxHy+ + e- ==> CxHy-1 + H
Lahfaoui et al. J. Chem. Phys., 106 (13)
1D transport model for H2 plasma flow
Some results
Comparison with experiments
MWPD = 9 W/cm -3 (MWPinp = 600 W and P = 25 mbar)
1,2E-02
2200
T(K)
2000
Tg-Model
Tv-model
Tg-CARS
Tv-CARS
1800
Substrate
1600
1400
Inlet (10 cm)
1,0E-02
8,0E-03
model
Actinometry
TALIF
6,0E-03
4,0E-03
Substrate
H-atom Mole Fraction
2400
2,0E-03
Inlet (10 cm)
0,0E+00
1200
0
1
2
z(cm)
3
4
0
Conclusion :
* Good Agreement on temperatures
* H-atom is overestimated in the boundary layer ( 50 %)
1
2
z(cm)
3
4
 Temperature
Fraction molaire d'hydrogène atomique
-3
30 W.cm
-3
9 W.cm
-3
12 W.cm
-3
15 W.cm
-3
23 W.cm
-3
30 W.cm
Gas temperature [K]
3000
2500
2000
-3
9 W.cm
1500
(a)
1000
1.8x10
-1
1.6x10
-1
1.4x10
-1
1.2x10
-1
1.0x10
-1
8.0x10
-2
6.0x10
-2
4.0x10
-2
2.0x10
-2
-3
30 W.cm
1
2
-3
9 W.cm
1
2
3
4
5
6
7
8
9
-3
9 W.cm
-3
12 W.cm
-3
15 W.cm
-3
23 W.cm
-3
30 W.cm
0.0
0
0
Substrate
 H-atom
3500
simulation domain
Radiale Uniformity
Axial profiles of temperature and hydrogen mole
fraction
3
4
5
6
7
8
9
10
Position axiale [cm]
10
Axial position [cm]
1/ H-atom relative density varies between 1 et 13 % when MWP 9 à 30 W.cm-3
2/ The discharge transitions from a cold non-equilibrium plasma to a thermal
plasma
Axial distribution of hydrocarbon species
30 W.cm-3 (120 mbar / 2 kW)
C 2 Hy
Tgmax = 3200 K
10
-1
10
-1
10
-2
10
-2
10
-3
10
-3
10
-4
10
-4
10
-5
10
-5
10
-6
10
-6
10
-7
10
-7
10
-8
10
-8
10
-9
10
-9
10
CH4
CH
CH3
C
1
CH2
CH2
5% CH4
Substrat
-10
0
1
Fraction molaire
Fraction molaire
CHx
2
Plasma
3
4
5
6
7
8
9
10
Position axiale [cm]
10
C2H2
C2H
C2H4
C2
C2H6
C2H3
C2H5
5% CH4
Substrat
-10
0
1
2
Plasma
3
4
5
6
7
8
9
10
Position axiale [cm]
1/ C2H2 is the major species (Actually we have H2/C2H2 plasma)
2/ Strong density variations in the reacting boundary layer (spatial stiffness)
3/ Significant amount of C and C2 species
4/ Caution : Interpretation of line of sight measurements  outside or inside the discharge ?
Example : validation of predicted hydrocarbon
species densities
TDLAS measurements (INP : Greifswald – J. Ropcke)
60 cm
Tg~3000 K
2 cm
Tg=300 K
Tg=600 K
5 cm
Fraction molaire de C2H6
Fraction molaire de CH3
10
CH3
-3
Code Radial moyenné
Mesure IR moyennée
10
-4
0
5
10
15
20
25
30
-3
Densité de puissance [W.cm ]
-3
C2H6
10
-4
Code Radial moyenné
Mesure IR moyennée
10
-5
0
5
10
15
20
Bras optique
0
10
-1
10
-2
10
-3
10
-4
10
-5
10
-6
10
-7
10
-8
10
-9
10
-10
10
C2H2
Fraction molaire
Absorption measurements yield :
=> CH3 in the discharge
C2H6 outside the discharge
The model allows to get the spatial distribution
and the validation takes place on spatially averaged values
10
25
30
-3
Densité de puissance [W.cm ]
25 cm
2.5 cm
CH4
CH3
C2H4
0
5
CH4
CH3
C2H2
C2H4
C2H6
C2H6
10
15
20
25
Position radiale [cm]
30
1D transport model for H2/CH4 plasma flow
Comparison with IR absorption measurements
Low power density – 9 W/cm-3
C2H2+ :
e-+C2H2 => 2e- + C2H2+
e-+ C2H2+ => C2H + H
CH4 + C2H2+ => C2H3+ + CH3
C2H5+ :
C2H4 + C2H3+ => C2H5+ + C2H2
e-+ C2H5+ => C2H4 + H
C2H3+ :
CH4 + C2H2+ => C2H3+ + CH3
H3+ + C2H2 => C2H3+ + H
e-+ C2H3+ => C2H2 + H
High power density – 30 W/cm-3
C2H2+ :
e-+C2H2 => 2e- + C2H2+
e-+ C2H2+ => C2H + H
C2H3+ :
H3+ + C2H2 => C2H3+ + H
e-+ C2H3+ => C2H2 + H
C2H5+ :
C2H4 + C2H3+ => C2H5+ + C2H2
e-+ C2H5+ => C2H4 + H
Open problems in H/C moderate pressure microwave
discharges for carbon films deposition
Soot
 Treatment of more complex discharges ….
Sooting discharges for nano and ultranano-crystalline diamond
deposition (PAH and soot formation) – undergoing work
 Self consistent treatment of pulsed regime – undergoing work
Electron velocity distribution function : anisotropy effect
(nel  wexc)
 EM field simulation : Resonance region and spatial stiffness
 More detailed investigation of the effect of gas heating on the
establishment of discharge regimes