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(v1) + H2(w 1) * v-t exchange : H2(v) + H2 (or H) ===> H2(v1) + 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-leTe) –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-mTv-m) + Qv-t - Qe-v + Qv-C = 0 t + div[Suihi –lt-rT -Slv-mTv-m –leTe] - 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