A Quantum Chemical Exploration of SiC Chemical Vapor Deposition

(1)

Linköping Studies in Science and Technology Dissertation No. 1828

A Quantum Chemical Exploration of

SiC Chemical Vapor Deposition

Pitsiri Sukkaew

Department of Physics, Chemistry, and Biology (IFM) Linköping University, SE-581 83 Linköping, Sweden

(2)

Printed in Sweden by LiuTryck 2017

ISSN 0345-7524

(3)

Abstract

SiC is a wide bandgap semiconductor with many attractive properties. It has attracted particular attentions in the areas of power and sensor devices as well as biomedical and biosensor applications. This is owing to its properties such as large bandgap, high breakdown electric field, high thermal conductivities and chemically robustness. Typically, SiC homoepitaxial layers are grown using the chemical vapor deposition (CVD) technique. Experimental studies of SiC CVD have been limited to post-process measuring of the layer rather than in situ measurements. In most cases, the observations are presented in terms of input conditions rather than in terms of the unknown growth condi-tion near the surface. This makes it difficult to really understand the underly-ing mechanism of what causes the features observed experimentally. With help of computational methods such as computational fluid dynamic (CFD) we can now explore various variables that are usually not possible to measure. CFD modeling of SiC CVD, however, requires inputs such as thermochemical properties and chemical reactions, which in many cases are not known. In this thesis, we use quantum chemical calculations to provide the missing details complementary to CFD modeling.

We first determine the thermochemical properties of the halides and halohy-drides of Si and C species, SiHnXm and CHnXm, for X being F, Cl and Br

which were shown to be reliable compared to the available experimental and/or theoretical data. In the study of gas-phase kinetics, we combine ab ini-tio methods and DFTs with convenini-tional transiini-tion state theory to derive ki-netic parameters for gas phase reactions related to Si-H-X species. Lastly, we study surface adsorptions related to SiCVD such as adsorptions of small C-H and Si-C-H-X species, and in the case of C-C-H adsorption, the study was ex-tended to include subsequent surface reactions where stable surface products may be formed.

(4)

Sammanfattning

SiC är ett halvledarmaterial med stort bandgap som har många attraktiva egen-skaper. Det har fått särskilt mycket uppmärksamhet inom kraftkomponenter och sensorer, men också för tillämpningar inom biomedicin och biosensorer. Detta beror på materialets stora bandgap, förmågan att klara höga elektriska fält, goda värmeledningsförmåga och kemiska stabilitet. Epitaxiska skikt av SiC för kommersiella tillämpningar tillverkas med hjälp av kemisk ångdepo-sition (Chemical Vapor Depoångdepo-sition, CVD). Experimentella studier av SiC CVD begränsas till mätningar som görs i efterhand, snarare än under proces-sens gång. För det mesta görs observationer med utgångspunkt från procesproces-sens ingångsvärden istället för från de okända tillväxtförhållandena vid ytan. Detta gör det svårt att verkligen förstå de underliggande mekanismerna för de obser-verade experimentella resultaten. Med hjälp av beräkningsmetoder, såsom computational fluid dynamics (CFD) kan vi utforska olika variabler som van-ligtvis inte går att mäta. Dock kräver CFD-modellering av SiC CVD ingångs-data i form av termokemiska egenskaper och kemiska reaktioner, och dessa data är ofta okända. I den här avhandlingen har vi använt kvantkemiska beräk-ningar för att ta reda på de detaljer som saknas för fortsatt CFD-modellering. Till att börja med har vi bestämt de termokemiska egenskaperna för halider och hydrohalider av SiHnXm och CHnXm, där X är F, Cl och Br, och som visats

vara tillförlitliga data jämfört med tillgängliga experimentella och teoretiska data. För att studera gasfas-kinetik har vi kombinerat ab initio och DFT-beräkningar med konventionell ”transition state theory” för att härleda kine-tiska parametrar för gasfas-reaktioner för Si-H-X species. Slutligen har vi studerat ytadsorption relaterad till SiC CVD, såsom adsorption av små C-H och Si-H-X species, och som i fallet C-H adsorption utökats till att även in-kludera efterföljande ytreaktioner där stabila ytprodukter kan bildas.

(5)

Acknowledgement

I would like to acknowledge my supervisors, Dr. Örjan Danielsson, Prof. Erik Janzén, Dr. Olof Kordina and Prof. Lars Ojamäe, for the assistance and sup-port during the study, for introducing me into the field of CFD modeling and quantum chemistry, and for giving me opportunities to learn and broaden my skills in various areas. This work would not have been possible without guid-ance and encouragement from all of you.

I would also like to thank all my friends and colleagues especially Pontus, Ildiko, Henke, Daniel, Ian, Martin, Chams, Zhafira, Aek, Thang, Duc, Vo-lodymyr, Christina, Elham, Ted, Abeni, Mama Lee, Yuttapoom, Prae, Mike, Sit, Björn, Louise and Chao for interesting discussions, enjoyable dinners and lunches and for support and inspirations.

I would like to thank my Dad and Mom, my sister and my cats for being there for me, for all the laughter and for always cheering me up.

(6)

List of included publications

Paper

1

Shortcomings of CVD modeling of SiC today

Ö. Danielsson, P. Sukkaew, L. Ojamäe, O. Kordina and E. Janzén Theoretical Chemistry Accounts (2013) 132, 1398.

My contributions: I calculated all thermochemical properties, performed parts of CFD kinetic model and wrote part of the paper

Paper

2

Thermochemical properties of halides and halohydrides of

sili-con and carbon

P. Sukkaew, L. Ojamäe, O. Kordina, E. Janzén and Ö. Danielsson ECS Journal of Solid State Science and Technology (2016) 5, P27.

My contributions: I calculated all thermochemical properties and wrote the entire paper

Paper

3

Silicon chemistry in fluorinated chemical vapor deposition of

silicon carbide

P. Stenberg, P. Sukkaew, I. Farkas, O. Kordina, E. Janzén, L. Ojamäe, Ö. Da-nielsson and H. Pedersen

J. Phys. Chem. C, Just Accepted Manuscript DOI: 10.1021/acs.jpcc.6b10849

My contributions: I calculated all surface chemistry and wrote part of the pa-per

Paper

4

Adsorption of silicon halides and silicon halohydrides on

4

H–

SiC during SiC-CVD process.

P. Sukkaew, L. Ojamäe, O. Kordina, E. Janzén and Ö. Danielsson In manuscript.

My contributions: I calculated all surface chemistry and wrote the entire paper

Paper

5

Ab initio study of growth mechanism of

4

H−SiC: adsorption

and surface reaction of C

2

H

2

, C

2

H

4

, CH

4

, and CH

3

P. Sukkaew, L. Ojamäe, O. Kordina, E. Janzén and Ö. Danielsson J. Phys. Chem. C, Just Accepted Manuscript

DOI: 10.1021/acs.jpcc.6b11085.

(7)

Paper

6

Opening the black box of silicon carbide chemical vapor

depo-sition

P. Sukkaew, P. Stenberg, O. Kordina, E. Janzén, H. Pedersen, L. Ojamäe, and Ö. Danielsson

Submitted.

My contributions: I calculated gas-phase kinetic reactions of the F, Si-H-Cl and Si-H-Br systems

(8)

Chapter 1 Introduction

1.1 Silicon Carbide

SiC is a wide bandgap semiconductor with many attractive properties. As the Si technology reaches its maturity, the investment has moved on to developing new materials with physical properties superior to Si. SiC has attracted partic-ular attentions in the areas of power and sensor devices as well as biomedical and biosensor applications. This is owing to its properties such as large bandgap, high breakdown electric field, high thermal conductivities and chem-ically robustness 1–8_{. Table 1 shows a summary of SiC properties at room}

tem-perature in comparison to Si and other wide bandgap semiconductors. The usage of wide bandgap semiconductors in power device applications enhance the efficiencies in energy transformations, while making it possible to produce devices in smaller sizes and be more robust in operations at high temperatures

1,9,10_{. It is also expected that SiC unipolar power devices will replace Si bipolar}

devices over the blocking range of 300 – 6500 V, and there is a promising fu-ture for these devices in the ultrahigh (> 10 kV) voltage range 1_{. The chemical}

robustness of Si-C bond and large Young’s modulus also make SiC attractive in sensor applications featuring high temperature and hostile environments

4,5,11_{as well as biomedical and biosensor applications}12,13_.

Table 1 Properties of SiC in comparison to Si and other wide bandgap

semiconduc-tors.

Properties SiC Si GaN Diamond

4H 6H 3C Energy gap 1,4,14,15 _3.2 _3.0 _2.3 _1.12 _3.4 _5.5 Thermal conductivity (W cm-1 _K-1₎14,16 4.9 4.9 3.6 1.3 1.3 20 Electron mobility (cm2_V-1_s-1₎4,14,17 1000 400 800 1500 900 2200 Breakdown field (MV cm-1₎4,7,14 ~3 2 3.2 0.3‒0.6 ~3 10 Electron saturation velocity (cm s-1₎14,17 2×107 _2×107 _2×107 _1×107 _2×107 _3×107 Melting Point (K) 4,16 ₃₁₀₃ ₃₁₀₃ ₃₁₀₃ ₁₆₉₀ ₂₆₇₃ ₁₆₇₃ Young’s modulus (GPa) 4 300 to 500 130‒180 200‒300 1000 Relative dielectric constant 6,11,14 10 9.7 9.6 11.9 8.9 5.7

(12)

SiC crystallizes in both cubic zincblende and hexagonal wurtzite structures and exhibits more than 200 polytypes 18–20_{. These polytypes, distinguished by}

the lattice types and stacking orders, are shown to have rather similar proper-ties regarding their physical, chemical and mechanical aspects 15_{. The}

differ-ences are seen in the band gap and band structures and are shown to vary as functions of hexagonalities 15_{, which are defined as the percentages of atoms}

in hexagonal positions in the unit cells. Among SiC many polytypes, the ma-jority of research are focused on the 6H, 4H and 3C polytypes, shown in Fig-ure 1.1. The first two polytypes are hexagonal and the last is cubic. The 6H and 4H polytypes have been used as substrates for the group III-N heteroepi-taxial growth 21_{, graphene}22_{as well as homoepitaxial growth of SiC. For GaN}

epitaxy, 6H and 4H SiC substrates have advantages in having small lattice constant mismatch to GaN and good thermal conductivity 21_{. At present both}

4H- and 6H-SiC substrates are commercially available in large size and good quality. For electronic applications, the 4H polytype is the most common ow-ing to its technological status and physical properties. More reviews on the SiC high power applications may be found in Ref. 1,2,6,23–25_{. The}₃_{C polytype,}

4H-SiC A B C B A 3C-SiC A B C A B A B C A C B A 6H-SiC

Figure 1.1 Stacking orders of 3C, 4H, and 6H polytypes of SiC. The letter A, B and C refer to the three stacking positions of a closed packed structure. The 3C polytype has 3 repeated units, ABC. The 4H polytype has 4 repeated units, ABCB, and the 6H poly-type has 6 repeated units, ABCACB.

(13)

on the other hand, is preferable for micro electro mechanical systems (MEMS) sensors for its ability to grow on Si wafers and hence compatible with conven-tional Si based MEMS processes 4,5_.

1.2 Chemical Vapor Deposition (CVD) of SiC

The progress in SiC technology has sped up drastically since the emergence of high quality and large SiC wafers grown by seed sublimation method 26,27_.

Active layer of SiC is commonly carried out homoepitaxially using chemical vapor deposition (CVD) technique 28,29_{. Chemical vapor deposition (CVD) as}

its name suggests is a process where deposition occurs via chemical reactions from vapor phase precursors 30_{. The process has been utilized to grow various}

kinds of materials as reviewed by Ref. 31_{. The quality of homoepitaxial layers}

depends directly on the ability to control the uniformity, layer thickness and doping density while at the same time suppressing the formation of defects and, particularly in SiC homoepitaxy, avoiding polytype mixing. The polytype mixing can be reduced or completely eliminated by a method called “step-controlled epitaxy” 32_{. This technique uses surface steps to replicate the}

sub-strate polytype to the layers. By cutting a subsub-strate at a slight off angle from the usual c-plane, more steps are produced on the surface with narrower ter-races lying in between. As the widths of terter-races become narrower, the species adsorbed somewhere on the terraces can now migrate on the surface to reach the step edges and the growth occurs steadily from the step edges as a result. This way of growth ensures that the stacking order of the substrate is passed on to the growing layer resulting in a single crystal. In this way it is possible to achieve stable growth at lower temperature of polytypes that otherwise would require very high temperature such as 4H and 6H at a lower tempera-ture 32_.

CVD of 4H polytypes is usually carried out at the temperature between 1500 to 1600 ℃ and the pressure of 100-300 mbar. Typically, the process proceeds using precursors such as silane (SiH4) and light hydrocarbons such as ethylene

(C2H4) or propane (C3H8). Alternatively, usage of precursors such as

methyl-trichlorosilane (CH3SiCl3) have also been reported 28,29. During the process,

the precursors are pre-diluted with a carrier gas (usually hydrogen) and fed into the heated growth chamber where the substrate is placed and well heated. When heated sufficiently, the precursor molecules start to react via pyrolysis reactions producing reactive gaseous species which when reaching the surface will get adsorbed into a new layer. The standard chemistry of SiC CVD is carried out with only the precursors and the carrier gas in the gas flow.

(14)

How-ever, problems occur when attempting to increase the growth rate by feeding in more precursors, i.e. by increasing the precursor concentrations. It is ob-served that when high concentrations of precursors have been used so-called silicon droplets start to form in the gas and damage the growing surface when they get into contact 29_{. This can be suppressed by adding halogen into the}

system to dissolve the droplets. With halogen addition, the growth rate can be increased significantly, as high as a few hundred micrometers per hour 29_{. In}

most cases, Cl has been used for this purpose. Recent research also found that Br and F chemistries are applicable to SiC CVD 33–35_.

Experimental studies of SiC CVD have been limited to post-process measur-ing of the layer rather than in situ measurements. In most cases, the observa-tions are presented in terms of input condiobserva-tions at the inlet and reactor setup such as precursor flow rates, C/Si ratios, doping agent flow rates, process temperature etc. rather than in terms of the actual condition on the surface. It should be noted that the inlet conditions are not necessary the same as the condition just above the surface. Most importantly, the conditions at the inlet and on the surface might not be related via linear relationships 36_{. Without}

knowing the actual surface condition, it is difficult to elucidate the underlying mechanism. This can be remedied by using various modeling methods.

1.3 Present Status of SiC-CVD Process Modeling

Computational modeling can be used to study what is happening inside the CVD growth chamber. Modeling of CVD-related phenomena has been carried out by many different methods, depending on the time and length scales of the phenomena. For calculations in macroscale (reactor scale), computational fluid dynamic (CFD) is a powerful method which has been applied to various types of CVD processes and materials allowing predictions of heat and flow distributions inside the growth chamber, types and concentration profiles of gaseous molecules and radicals, profiles of growth rates and doping 37–40_.

Us-age of CFD ranges from designing of reactor geometry to process control. CFD studies of SiC CVD were reported in Ref 38,41–43_.

For CVD calculations, CFD requires detailed input of all aspects ranging from material properties in heat and mass transports to elementary chemical reac-tions occurring in the process. The methods for solving the mass and heat transport equations are well developed and the calculations are reliable 38,42_.

(15)

part and a surface part. For SiC deposition, the gas phase part involves C-Si-H-halide containing species. The reactions involving C-H and Si-H species, both experimental and theoretical data, are available in the literature 38,41,44–47

and recently Cl-related reactions 48–51_{have also become available based on ab}

initio and density functional theory (DFT) calculations. The data have been used extensively in CFD modeling 36,51–53_{. As the calculation cost is high, it is}

preferable to preselect only the important reactions rather than using the entire extensive set. In this way the computational cost can be reduced without los-ing important details. However, for SiC CVD such reduced reaction sets are shown to work well only within limited conditions and need to be generalized to cover wider process conditions in order to be useful 53_{. Alternatively, it is}

possible to predict the gas phase compositions directly from their thermody-namic properties 39,54,55_{. The derived compositions are the condition as the}

system reaches its thermodynamic equilibrium and hence when its total Gibbs free energy is the lowest. It is a useful way to explore the effects of changes in process conditions on the species formations both in gas phase and condensed phase. For this to be reliable it is important to include all important species including the short-lived radicals where thermodynamic properties are usually not available in the commonly used databases. Quantum chemical calculations such as ab initio, DFT and their composites are great tools to derive such properties, especially for short-lived species which are hard to measure exper-imentally 55–57_.

The details of surface reactions are much less studied than the gas phase. Due to necessity surface models for CFD in the literature usually combine experi-mental results carried out under different conditions than the CVD 58–63_.

Com-putational studies on adsorptions on SiC surfaces are also rare and unsystem-atic 64–70_{. A useful model for thin film growth needs to include kinetics}

occur-ring at different positions (step edge, terrace, surface configuration, etc.). Sur-face adsorption studies, as well as studies of kinetics and thermochemistry, have been performed extensively using DFT. DFT has its strength in its accu-racy and costs compared to other ab initio methods 71,72_{. Time and length}

scales for DFT are in the order of 1 ps to 1 ns and 1 Å to a few nm respectively

73_{. With DFT kinetic barrier heights can be obtained with high accuracy and}

together with transition state theory the actual rate can be calculated. It should also be noted that each DFT calculation provides information for a certain reaction on a certain surface configuration. To build up a surface model for CFD requires multiple DFT calculations covering a wide range of reactions and surface configurations.

(16)

Alternative methods to simulate micro to mesoscale such as molecular dynam-ics (MD) and kinetic Monte Carlo (kMC) can also be applied to study surface phenomena 65,72,74–76_{. These types of methods allow much larger coverage in}

time and length scales than DFT. MD length scale ranges from a few nm to 100 nm and time scale of a few µs to 1 ms 73_{, while kMC time scale is in the}

order of a few ms to seconds and at a length scale of 100 nm to 1 µm 73_{. Unless}

combined with atomic level method such as DFT as in the case of ab initio MD, empirical data will be used to determine the time evolution in such stud-ies. Evolution of atoms and molecules in MD is directed by the use of effec-tive force fields. However, at the present time, MD empirical force fields are not yet suitable for semiconductor growth applications due to their lack of parameters for the details to predict correct chemical reactions involved in the process. A compromising way to solve this problem is to add a quantum me-chanical method such as DFT to calculate chemical reactions at a certain, fixed interval of time. This ab initio MD however has its drawback in its lim-ited time and length scales due to the usage of DFT. As the time scale is short-ened, the chance of actually observing chemical reactions happening during the run is also small, especially for reactions with large energy barriers. kMC can also be empirical based but different from MD in the methods describing time evolution. Using kMC the system is allowed to evolve based on probabil-ity and transition rate of events 73_{. The calculations can be performed at a time}

scale more appropriate to growth and at the same time cover a large number of atoms and molecules. In this way, surface morphology can be taken into ac-count and surface evolution at each positions can be observed. But a drawback with atomistic simulation techniques such as MD and kMC are that lengthy separate simulations need to be performed for e.g. each temperature or pres-sure. More analytic kinetic simulation techniques like CFD are exceedingly more efficient when it comes to e.g. predicting concentrations and their time-dependence of different species, which are the quantities most of interest in order to understand the CVD process.

(17)

Chapter 2 Quantum Chemical Calculations

2.1 Introduction

Quantum chemical calculations apply quantum theory to study chemical relat-ed processes. The goal is to determine electronic structures of atoms and mol-ecules, from which the energies and equilibrium geometries can be derived. At atomic level, the particles are described in term of wave functions which are governed by the Schrödinger equation. The non-relativistic, time-independent Schrödinger equation in atomic units is written as

(2.1) where

The Hamiltonian is the energy operator containing (listed from left to right) the electron kinetic term ( ), the nuclei kinetic term ( ), the electron-nuclei interaction ( ), the electron-electron interaction ( ) and lastly the nuclei-nuclei interaction ( ). is the wave function describing the system.

The Schrödinger equation is a many-body equation and thus it is impractical to solve directly unless the system under study is extremely small. To remedy this, approximations are usually made. The first approximation considered here is the Born-Oppenheimer (BO) approximation, which states that the elec-trons adjust themselves to the motion of the nuclei or in other words the nu-cleus movement is separable from that of the electrons. The BO assumption is based on the fact that the nuclei are much heavier than the electrons. It has been shown to work well unless two (or more) solutions come close in ener-gies. With BO approximation, the Schrödinger equation can be separated into two parts, the first describing the electronic wave function and the latter

(18)

de-scribing the nucleus wave function. The electronic Schrödinger equation is written as

(2.2a)

where ,

and the nucleus Schrödinger equation is

(2.2b)

where and .

The Eq. 2.2a and 2.2b satisfy the total Schrödinger equation,

(2.2c)

where ,

It should be noted that, in the electronic Schrödinger equation, Eq. 2.2a, the nuclei positions are treated as fixed parameters, while in the nuclear wave function equation, Eq. 2.2b, the presence of electrons are treated as a potential acting on the nuclei. It is clear by comparing Eq. 2.1 and Eq. 2.2a - 2.2c that the BO approximation removes the coupling terms between the kinetics of nuclei and electrons.

2.2 Hartree-Fock Approximation

The Hartree-Fock (HF) method 77–80_{reduces the N-particle Schrödinger}

equa-tion, Eq. 2.2a, into a one-particle equation by treating as an average field acting on each electron. The HF method uses a Slater determinant as a trial wave function,

(19)

. (2.3)

Here is a single-electron wave function, called spinorbital, which is a prod-uct of a spatial orbital and a spin function. is the number of electrons in the system. The determinant formulation ensures that the wave function is anti-symmetric with respect to an interchange between two electron coordinates. With the wave function written this way, it is inherently assumed that the wave functions can be separated into a product of single-electron spinorbitals. Moreover, by using a single determinant as a solution, HF entirely neglects the effects of electron correlation.

With the Slater determinant ansatz, Eq. 2.3, the electron-electron repulsion in Eq. 2.2a is separated into two terms,

,

.

(2.4)

The first term ( ) is the classical Coulomb interaction, and the second term ( ) is the exchange interaction between two electrons. By choosing to be orthonormal to each other, the exchange term is nonvanished only for the spinorbitals with the same spins. Using variation method, the one-electron HF equation is written as

(2.5)

where ,

,

(20)

and ,

where and refer to the spatial coordinate and the spatial plus the spin

co-ordinate respectively. The operator is called the Fock operator. The operator is a one-electron Hamiltonian containing the electron kinetics and the po-tential from the nuclei.

With the HF orbitals given by Eq. 2.5, the HF wave function can be formulat-ed using Eq. 2.3 with the energy derived from the electronic Schrödinger equa-tion Eq. 2.2a,

.

(2.6a) Alternatively can be written in term of HF orbital energy, , as

. (2.6b)

It should be noted that the total energy, , is not equal to the sum of the orbital energies, . The sum of the orbital energies, , has in fact count-ed the repulsion effects twice and thus necount-eds to be correctcount-ed with the term

.

HF requires an initial guess of orbitals to start with, and from there the calcu-lation continues iteratively until the wave function and energy remain un-changed, or in other words until the self-consistency is reached. As it is based on the variational method, HF energy is guaranteed to be either equal to or larger than the exact value. Despite HF’s limitation in neglecting electron correlations, HF orbitals can be used to build up a trial wave function for other methods such as MP2 and CC, as well as to provide an exchange energy for hybrid DFTs.

From the HF determinant, a new trial wave function can be built by adding more determinants to the HF wave function. This can be achieved by, for ex-ample, using a linear combination of determinants as in configuration

(21)

interac-tion (CI) method. The addiinterac-tional determinants can be built by replacing one or more orbitals in the HF determinant with virtual ones, and they are thus called excited Slater determinants. In a “frozen core” approximation, only the va-lence electrons are allowed to be excited. The excited determinants are classi-fied from the number of excited electrons. When one, two or three electrons are excited, the determinants are called respectively as singly, doubly, triply excited determinants. The virtual orbitals are generally obtained from an HF calculation whenever a basis set larger than a minimum set has been used. To get a reasonable description of the wave function, however, the number of virtual orbitals should not be too small. From the new trial wave function, electron correlations can now be taken into account by solving for a wave function that are a solution to the Schrödinger equation.

2.3 Møller-Plesset Many-Body Perturbation Theory

Correlation effects can be recovered using many-body perturbation theory 77

where the Hamiltonian, , is written in terms of a reference Hamiltonian, , and a perturbed Hamiltonian, , with being a parameter determin-ing the strength of perturbation,

. (2.7) The perturbed Schrödinger equation is

. (2.8) with the energy ) and wave function ( ) written using Taylor expansion as follows,

, (2.9a)

and . (2.9b)

(22)

. (2.10) By mapping the terms with the same power of , we obtain the nth_-order

per-turbation equations. These equations provide formula for the nth_{-order energy}

corrections and the nth_{-order wave function corrections.}

In the nth_{-order Møller-Plesset (MP) many-body perturbation theory, MPn, the}

reference Hamiltonian is chosen to be the sum of the Fock operators 81_,

. (2.11)

With the parameter set to 1, the perturbed Hamiltonian is

. (2.12)

Here is the electron Hamiltonian Eq (2.2a). The 0th_{-order wave function,}

, is the HF determinant and the ith_{-order wave function,} _{, in Eq.}_2.9_{b is}

written in terms of the excited Slater determinants. For MP0 and MP1 levels, the energies are simply and respectively. The correlation correction starts from MP2 level with the energy correction,

. (2.13)

Here, is an excited determinant in which and occupied orbitals are replaced by and virtual orbitals. The summations and run over all occupied and virtual orbitals respectively. The reader is referred to a review paper by D. Cremer 82_{for explicit forms of MPn. It should be noted}

that the energies derived from MPn do not always converge monotonically toward the limiting value, as shown in Fig. 2.1a. Instead, it could oscillate around the value as n increases, as shown in Fig. 2.1b 82,83_{. Nevertheless, for a}

(23)

well behaved system, acceptable results can be expected between MP3 and MP477_.

Figure 2.1 Correlation (single point) energies of the MPn methods in the unit of Hartree of a) CH4 and b) HF calculated using cc-pVnZ basis sets, for n = D,

T, Q and 5. The correlation energies are calculated based on equilibrium struc-tures optimized using CCSD(T) level and cc-pVDZ basis set.

2.4 Coupled Cluster Theory

The coupled-cluster (CC) theory 77,84_{introduces an ansatz of the form,}

, (2.14)

where .

The operator is called the cluster operator and written as . The operator produces the ith_{-excited Slater}

determinant for the HF wave function, , such as

b)

(24)

and .

with the summations and running over all occupied and virtual orbitals respectively.

With the CC ansatz the Schrödinger becomes . As it is not possible to include all terms in the expansion, the operators which cre-ate higher excited stcre-ates are usually abandoned in the calculation. When includes only the single ( ) and double ( ) excitations, the method is called CCSD. In CCSD, the cluster operator produces

(2.15)

Figure 2.2 Correlation (single point) energies of a) HF and b) CH4

calculat-ed using CCSD in comparison to MP4. The correlation energies (in Hartree) are shown as functions of cc-pVnZ basis sets, for n = D, T, Q and 5. The correlation energies are calculated based on equilibrium structures optimized using CCSD(T) level and cc-pVDZ basis set.

(25)

The exponential function thus includes not only the single and double exci-tations, and , but also their disconnected products , , etc. An addition of connected to CCSD results in CCSDT but this is usually too costly for most applications. Nevertheless, connected excitation is usually important 77_{. As shown in Figure}_2.2_{a and}_2.2_{b, the correlation energy from}

CCSD is inadequate in comparison to the correlation energies derived from MP4. There are various methods derived to approximate an addition of the triple excitation. Among these, the most popular one is the CCSD(T) 85_{, in}

which the triple excitation is included by perturbation. CCSD(T) is shown to give a good compromise between the cost and accuracy, and in fact found to be more accurate than CCSDT by error cancelations between the overestima-tion of triple contribuoverestima-tion and the error of truncaoverestima-tions of higher excitaoverestima-tions

86,87_.

2.5 Density Functional Theory (DFT)

DFT is similar to HF in reducing the many-body equation into a single-particle equation. It is based on the theorem by Hohenberg and Kohn 88_(HK)

which demonstrates that the ground state energy of electronic system, as well as other properties at ground state, are unique functionals of the electron den-sity, . In the same article 88_{, HK also introduced a “universal”}

functional of electron density which can be used together with any external potential in order to derive the electronic energy and demonstrated that the correct density is obtained by minimizing the energy functional.

Kohn and Sham 89_{(KS) expanded the work of HK and formulated a set of}

one-particle Schrödinger equations similar to the Hartree approach. They pro-posed that the problem of N interacting particles can be solved using an auxil-iary system of N uninteracting particles subjected to a potential V. The poten-tial V is in fact a term describing the particle interactions. The auxiliary sys-tem can then be described using,

(26)

where is the kinetic energy of non-interacting particles and is the potential resulting from the particles interactions.

The potential is composed of 1) which is an external potential from the nuclei, 2) the term which is a potential created by an electronic density and 3) which is a new potential containing all N-particle interactions. In this way, Eq. 2.16 puts all the complexities into the unknown exchange-correlation term, 90_{. Up to this point, KS formulation is exact}

and with the correct exchange-correlation term it provides the density. Similar to HF, DFT has the advantage in its low computational cost in comparison to wave function-based electron correlation methods such as MPn and CC. The quest in DFT formulation has been to obtain the unknown Exc. KS89

pre-sented the first Exc using the homogeneous electron gas assumption, which is

valid only for a slowly varying density . The thus depends only on the density and is called the local density approximation (LDA). The spin generalized case of LDA, called the local spin density approximation (LSDA), assumes that is a function of both spin densities, and

. This type of typically underestimates the exchange part of Exc by

5-10% 91_{. An improvement by including density gradients} _{leads to the}

sec-ond type of functionals called the generalized gradient approximation or GGA. GGA functionals are shown to provide improvement over LDA in the energy accuracy, but similar to LDA, GGA suffers from overbinding tendency but with smaller amounts than LDA 91_{. There are many gradient corrected}

functions available in the literature. The famous GGA functionals are B8892_,

PW8693_{, PBE}94_{and LYP}95_.

For quantum chemical calculations in which chemical bonding and thermo-chemistry are the main focus, hybrid DFTs and their descendants such as meta hybrids have been used extensively and with high success. Hybrid Exc includes

the HF exact exchange, , to LDA or GGA functional to reduce the overbinding tendency 96_{. The most popular hybrid function is B}₃_{LYP, which}

uses the B3 exchange functional from 96_{and the correlation}

func-tional from LYP 95_{. B}₃_{LYP has been used intensively in quantum chemistry}

with more than 60,000 citations 91,97_{. It was shown to be reliable for}

(27)

molecules 98,99_{. At the present, there are new functionals available which}

out-perform B3LYP in some aspects especially for kinetic barrier heights 100,101_but

at the same time more expensive 102_{. For a system with many electrons, a}

compromise then has to be made between the cost and accuracy.

2.6 Basis Sets

In solving the Schrödinger equation using ab initio methods, an orbital is usu-ally written as a linear combination of known functions, called basis sets. For these basis functions to describe an orbital correctly, they should form a com-plete set. In practice, a finite basis set is used and the accuracy thus depends on how well these functions can represent the orbitals and how many func-tions are included. In general, there are two error sources in ab initio methods: 1) the method inherited error and 2) the error due to lack of basis set conver-gence.

The basis functions used in this thesis are based on the Gaussian type orbitals (GTOs),

, (

2.17

a)

or

. (

2.17

b)

In Eq. 2.17a, the function is written in spherical coordinates. The integers , , are the principle, angular and magnetic quantum numbers respectively. is the exponent parameter determining the radial extent (how loose or tight) of the orbital. is spherical harmonic functions. In Eq. 2.17b, the function is written in Cartesian coordinate and are parameters determining the types of orbital, i.e. s, p, d, etc.

The smallest possible set of basis functions has one s-function for the first row atoms, two s-functions for Li and Be, and two s-functions and a set of p-functions (2px, 2py, 2pz) for the rest of the second row, etc. If the number of

the functions is doubled, tripled or quadrupled, with different exponents in each set, the new basis sets are called a double zeta (DZ), a triple zeta (TZ), and a quadruple zeta (DZ), respectively. As it is known that the core electrons are much less involved in chemical reactions than the valence, an addition of

(28)

Figure 2.3 Basis set convergence for correlation (single point) energies of CCSD, CCSD(T) and MP4. The figures a) to d) show the correlation ener-gies of 3_C,4_N,3_{O, and}2_{F respectively. The correlation energies are shown}

as functions of cc-pVnZ basis sets, for n = D, T, Q and 5 together with the values at the complete basis set limits, shown as DTQ and TQ5. The basis set limits are extrapolated using the relation

from Peterson et al. 199429_{. The value at DTQ uses the}

ener-gies from the cc-pVDZ, cc-pVTZ and cc-pVQZ basis sets and TQ5 uses the energies from the cc-pVTZ, cc-pVQZ and cc-pV5Z basis sets.

(29)

basis functions can be done only for the valence electrons, which is called a spilt valence basis.

The functions with larger exponents, , are appropriate in describing tighter bonds such as bonds, while those with smaller exponents allow the orbital to be more diffuse and are needed when the electrons are loosely bound. In most cases, basis functions with higher angular momentum are also needed to properly describe electron distribution along the bond, i.e. describing polariza-tion. They are especially important in describing electron correlapolariza-tion.

In adding more basis functions into the set, it is also important to consider the balance of basis sets. As a rule of thumb, the number of added functions should be maximum at the function with the lowest angular momentum, i.e. s. The number of the next lowest angular momentum function, p, should be one less and so on. For example, a balance basis is 5s4p3d2f1g. But as the func-tions f and g are much less important, they might also be left out, thus leaving the basis simply as 5s4p3d.

In most cases, the computational speed can be improved by contracting the primitive basis sets into contracted basis sets. Contraction combines the primi-tive basis functions into smaller groups of fixed linear combinations. A linear combination of these groups are then used in place of the primitive set. This introduces a new set of expansion coefficients, the number of which are much smaller than those of the primitive set. With the number of variables (expan-sion coefficients) reduced, so does the computational cost. Contraction is found to be useful especially for inner electrons which are much less involved in the chemical process and thus less effected by the environment. Contracted basis sets are for example Pople basis sets and Dunning-Huzinaga basis sets. An important consideration for choosing a basis set is the error it may produce in comparison to the value at the complete basis set (CBS) limit. Without a well-define path to the CBS value, it is difficult to predict the error and deter-mine from where it derives. This leads to constructions of correlation con-sistent basis sets and their descendants 103–105_{. With these types of basis sets,}

the correlation energies converge smoothly to a limiting value with increasing in the size of the basis sets 87_{, as shown in Figs.}_2.3_{. The smooth, converging}

behavior of the correction energies allow mathematical prediction of the CBS values, as shown in Figs. 2.387,106_.

(30)

2.7 Composite Methods

In general, despite the demand in accuracy, it is impractical to perform a high level calculation together with a basis set close to the basis set limit. Alterna-tively, a combination of different levels of theory (i.e. different methods and basis sets) can be used to improve the accuracy but in a way that reduces the cost in comparison to a direct calculation at a high level. For example, the geometry optimization is less sensitive to the level of theory than the energy

77_{. Hence, to reduce the calculation cost, a lower level of theory can be used}

for geometry optimization and from there the energy is then calculated. Vari-ous composite methods have been proposed in the literature, the choice of which depends on the level of accuracy and cost. Among these are the Gaussi-an-n theory 107–110_{, the complete basis-set method (CBS)}111–115_{, HEAT}116,117

and the Weizmann-n theories 118–120_{. In this thesis, the Gaussian-}₄_{theory and}

the Weizmann-1 theory are used to calculate the thermochemical properties of gases.

The Gaussian-4 theory (G4) 110_{provides predictions of thermochemical data}

by combining a series of ab initio calculations with empirical corrections. The method was designed to be computationally relatively light so that its applica-tion covers not only small molecules but also those of intermediate sizes. The ab initio parts of the theory contain a series of single point calculations at the CCSD(T), MP4, MP2 and HF levels and various basis sets. The HF energy is extrapolated to the value at the CBS limit. The G4 uses empirical corrections (called higher level corrections) to take into account the remaining deficien-cies. These corrections were derived from experimental data (the Gn sets) containing enthalpies of formation, ionization energies and electron and ton affinities from various types of molecules. The theory was shown to pro-duce an average absolute deviation of 0.83 kcal/mol for 454 energies in the test set 110_.

The Weizmann-1 theory (W1), unlike the Gaussian-n theories, include no em-pirical correction terms. It has been designed for the first- and second-row molecules and has been shown to produce an average absolute deviation of 0.30 kcal/mol 121_{. The method}118,121_{uses extrapolations to the CBS values for}

the HF, CCSD as well as (T) corrections. W1 also includes ab initio calcula-tions for core correlation and scalar relativistic effects 118,121,122_{. With a CBS}

extrapolation at the CCSD(T) level, a large portion of correction effects is expected to be recovered. This however causes the method to be computation-ally much heavier than the Gn.

(31)

Chapter 3 Quantum Chemistry in CVD Modeling

3.1 Optimization and Normal Mode Vibration

Applications of quantum chemistry to the study of CVD allows plenty of structural and dynamical information to be obtained. Main tasks are to deter-mine the geometrical structure and the energy of the molecule, as well as those of the transition state. The equilibrium structure of the molecule and the struc-ture of the transition state are present on a potential energy surface (PES) as stationary points. Mathematically, they are described by the vanishing of the energy gradients with respect to all geometrical parameters. For a molecule, the stationary point forms a minimum on the PES, while for a transition state it forms a saddle point which is characterized by a maximum along the reac-tion coordinate linking two minima which are the reactant and the product states.

A method to locate a minimum on a PES is called a geometry optimization. An optimization process requires a guess of an initial structure from which the first and second derivatives of energy, called the gradient and the force matrix or Hessian, are derived with respect to the internal coordinates such as bond lengths, or alternatively to the Cartesian coordinates of each constituent at-oms. During the optimization, a new structure is created from the initial ent and Hessian matrix. This new structure is then used to create a new gradi-ent and Hessian, which then leads to a prediction of a new structure. The pro-cess continues iteratively until the stationary point is located 77_.

Once the structure is found, the Hessian is calculated, diagonalized and mass weighted. From the Hessian, normal mode frequencies can be derived, which for a stationary point being a minimum on a PES should all be positive values. For a transition state, on the other hand, one frequency is required to be imag-inary due to the force constant being negative. For a correct transition state, the directional movement associated to the imaginary frequency should depict structural change along the intrinsic reaction coordinate to the product state or in the opposite direction to the reactant state. Once the normal mode frequen-cies and the geometrical structure are known, we can then proceed to calculate the partition function, which is a central function to thermodynamics.

(32)

3.2 Partition Functions

It is known from statistical thermodynamics 123_{that the probability, , of}

finding the system in the state with an energy is determined by the Boltz-mann factor, , where and are the Boltzmann constant and tem-perature respectively. The normalization constant to is called the partition function, given by It is which determines other thermody-namic properties and thus lies at the heart of thermochemistry.

For an ideal gas system, the total energy, , can be written as a sum of the energies from each molecule, . This implies that for an ideal gas system, can also be written as a function of the partition functions of individual mole-cules, For a system of N indistinguishable molemole-cules, 123_where

the factor derives from the indistinguishability of the molecules. The mo-lecular partition function, , can also be separated into four parts according to the molecular degrees of freedom, i.e. translational, , rotational, , vibra-tional, , and electronic, . We now write the partition function, , follow-ing the method used in the Gaussian software 124_as123,125

, (3.1)

where ,

,

for a linear molecule,

for a nonlinear molecule,

.

Here , , , are the mass of the molecule, the volume, the pressure and the Planck constant. and are the molecule moment of inertia and symmetry

(33)

index respectively. is the vibrational temperature of vibrational mode , and is the spin multiplicity of the molecule. In the derivation of , we have assumed that the molecule is an ideal gas, a rigid rotor and a harmon-ic oscillator. The model is defined with respect to the bottom of the well, and thus includes the zero-point vibrational mode by definition. The approxi-mation derives from the assumption that the energy difference be-tween the ground state ( defined as zero) and the first excited state, , is much larger than . Hence, all excited states are inaccessible and only the ground state energy, , contributes to . This then reduces the

summa-tion to

Once the partition function is obtained, other thermodynamic functions can be calculated using the following relations 123

(3.2) , (3.3) , (3.4) and , (3.5) , (3.6) , (3.7) . (3.8)

Here , , , , , and are the energy, entropy, enthalpy, heat capacity, Helmholtz free energy, Gibbs free energy and chemical potential respectively.

(34)

3.3 Thermochemistry

From the partition functions, , we can now derive thermochemical properties related to the SiC-CVD process. Most reactor scale modeling requires the enthalpy, entropy and heat capacity of all species involved in the process. These are usually written in the form of 7-coefficient NASA polynomials fit-ted over a specified temperature range,

, (3.9)

, (3.10)

. (3.11)

where to are the fitting coefficients, and is the absolute enthalpies

defined as

The entropy and heat capacity are readily available from ab initio calculations using Eq. 3.3 and 3.5. The enthalpy function, , however requires more cal-culations. This is due to the presence of , the standard enthalpy of formation of the molecule at 298 K which requires empirical atomic proper-ties. The standard enthalpy of formation, , is the enthalpy change at the standard condition (1 bar) when the molecule is formed from its constituent atoms in their reference states, i.e. their most stable form at 1 bar. is derived from the standard enthalpy of formation at 0 K, or , and the thermal correction at 298 K, or as

(3.12)

(35)

and . (3.14) The thermal correction at 298 K, , can be derived from ab initio calculations, Eq. 3.4. is the atomization energy at 0 K, which is defined as the difference in the energy at 0 K between the constituent atoms and the mol-ecule and can be derived from the partition functions, , of the molmol-ecules and the relevant atoms. The data of atomic and atomic thermal correction

are usually available from databases such as NIST-JANAF 126_and

NIST-CCCBDB 127_{as presented in Table}_3.1_.

Table 3.1 Atomic and - relevant to SiC CVD obtained from databases and literature.

at 0 K (kJ/mol) - (kJ/mol) Si 446 ± 8 126_{, 448.5 ± 0.8}128 _3.218126 C 711.19 ± 0.46 126 _1.051126 H 216.035 ± 0.006 126 _4.234126 F 77.28 ± 0.30 126 _4.413126 Cl 119.621 ± 0.006 126 _4.591126 Br 117.92 ± 0.06 126 _12.255126

3.4 Gibbs Free Energy Minimization

From the Gibbs free energies of the molecules, we can predict the gas phase species composition as the system reaches equilibrium by minimizing the total Gibbs free energy. This method is simple as it does not require any kinetic inputs. To obtain reliable results, all relevant species should be included in the minimization. The derived gas phase species composition is thus the composi-tion which minimizes the Gibbs free energies of the system at the temperature and pressure of consideration. By repeating the process over a range of tem-peratures (and/or pressures), a diagram of compositional change can be ob-tained as a function of temperature and pressure.

(36)

In general, the change in the total Gibbs free energy is written as 123

, (3.15) with , , , , and being the entropy, temperature, volume, pressure, chemical potential and the amount (mole) of species , respectively. The min-imization of the Gibbs free energy with respect to all species must be subject-ed to the conservation of species, , where is the number of atoms of element contained in species , and is the total number (mole) of element in the system. For the process that takes place at constant pressure and temperature, the change in the Gibbs free energy becomes simply

.

3.5 Conventional Transition State Theory

In this section, we summarize the calculation method for kinetic parameters using ab initio calculations together with conventional transition state theory (TST). TST provides a method to estimate reaction rate constants when there exists a transition state between the reactants and products. The transition state (TS) is a saddle point on the PES with a maximum along the reaction coordinate between the reactant and product as depicted in Fig. 3.1. The TS structure produces one imaginary normal mode frequencies, , corresponding to the directional movement along the intrinsic reaction coordinate to the product state or in the opposite direction to the reactant state as shown in Fig. 3.2. TST assumes that a thermal equilibrium is maintained between the reac-tant state and the transition state, which at high pressure with many collisions, can be fulfilled. At other pressures, the theory cannot account for the pressure dependent effects on chemical reactions. It is also assumed that once the sys-tem has crossed the barrier at TS it will reach the product state without re-crossing back to the reactant state. When rere-crossing happens, TST predicts an overestimation of the true rate.

3.5.1 Gas Phase Kinetics

For a reaction , the assumption of equilibrium between the reactants and TS leads to 129

(37)

. (3.16)

where is the Gibbs free energy of activation of the reac-tion and = 1 bar. For = 1 mole, and using the molar partition function,

, the equilibrium constant can be written in terms of the molar partition functions, , which is shortened from , as 129

En er gy Reaction coordinate TS Reactants

Products

Figure 3.1 Illustration of potential energy from the reactant state to the transition state (TS) and to the product state of the reaction, SiCl3 (g) + HCl (g) ↔ SiCl4 (g) + H (g). The transition state corresponds to a saddle point in energy, in contrast to the reactant and product states which are minima on the PES.

Figure 3.2 Illustration of vibrational movement corresponding to the imaginary frequency of the transition state of the reaction, SiCl3 (g) + HCl (g) ↔ SiCl4 (g) + H (g).

(38)

, (3.17)

where . The rate of converting into the

product depends on the frequency, , and the probability of its success, , . (3.18)

This gives the forward rate as where . The

factor derives from the relation . The probability is called the transmission coefficient. It introduces tunnelling effects and recrossing probability into the rate constant. It is sometimes assumed to be 1. It can be shown that for , the rate constant can be written as 129

, (3.19)

or . (3.20)

where in Eq. 3.19, the vibrational mode corresponding to is discarded from and the same is applied to in Eq. 3.20. The rate constant from Eq. 3.19 and 3.20 can also be fitted to the Arrhenius form, with the parameter being the activation energy.

3.5.2 Surface Kinetics

A solid surface can be modeled from a surface of a cluster, provided that the cluster is not too small. To model a solid surface, we follow the method of Reuter and Scheffler 130_{where the cluster translational and rotational degrees}

of freedom are removed, and so as their partition functions, and . This adjustment is applied to the transition states as well as to the cluster with ad-sorbed species. Hence, only two parts, i.e. the vibrational and electronic parts, of the partition function of the cluster are included in the calculation,

(39)

where and are derived from Eq. 3.1. Following Reuter and Scheffler 130_{, the partition function of gaseous molecules related to a surface}

reactions are written as

. (3.22) It should be noted that is similar to Eq. 3.1 except that the translational degree of freedom along the reaction coordinate has been removed, leaving only 2 out of 3 translational degrees of freedom. The 2 dimensional, transla-tional partition function is written as

= . (3.23)

Here is the area in which the gas species may reside in order to react suc-cessfully with the surface. In our study, the value of is assumed equal to the area per one surface site or equally the area per one lattice site. For a sur-face site on the (0001) surface of 4H−SiC, the area is equal to 8.178 Å2_,

calcu-lated using the 4H lattice constant of 3.073 Å 131_.

Adsorption rate. Using TST, we can derive adsorption rate of gaseous

spe-cies , corresponding to the reaction , to be 130

, (3.24)

Here is written in the unit of molecule sites-1_s-1_{. Alternatively, the}

ad-sorption rates can also be rewritten in the unit of mole area-1_time-1_{by using}

and also in the unit of length time-1_{by using}

respec-tively. Here is the Avogadro constant and is the height of 1 SiC bi-layer. is the activation energy of adsorption process at 0 K. and are defined similar to Eq. 3.21 and is defined by Eq. 3.22. Here refers to the area per one lattice site on the surface and is the fraction of surface occupied by . is the impingement rate of species from the kinetic theory,

(40)

where is the mole fraction of gaseous species and is the total pressure such that the product produces the partial pressure of species , and is the mass of gaseous species .

We can also define the sticking coefficient of species , or , using the

rela-tion, which produces

. (3.26)

The sticking coefficient of species describes the probability of the species to stick on the surface per strike.

Desorption rate. Using TST, the desorption rate, , of an adsorbed species described by a reaction, ‒‒ , is shown to be

, (3.27)

This reaction is the reverse of adsorption process. and are defined similar to Eq. 3.21, and is the activation energy of desorption process. is the fraction of surface occupied by adsorbed species .

Surface Reactions. For an on-surface reaction, , the reaction rate derived from TST is

, (3.28)

where , are defined following Eq. 3.21. Here is the par-tition function of a cluster on which an adjacent pair of and resides.

is the fraction of surface occupied by such a pair. is the activation energy of the reaction at 0 K.

(41)

Chapter 4 Results and Summary

An understanding of SiC CVD is crucial to the improvement and design of process conditions. Computational modeling has become a useful tool for such purposes due to a lack of available in situ measurements. A reactor scale mod-eling such as computational fluid dynamic (CFD) requires details of kinetics and thermochemistry in order to calculate the gas mixtures, the dominant gas-phase species as well as the growth and doping profiles. The needed infor-mation however is not always available as we have pointed out in Paper 1.

Firstly, thermochemical properties of organosilicons are usually missing from the common databases. When they are available, they sometimes contain rela-tively large errors. We showed this by comparing the data provided from the CHEMKIN database to our results derived from ab initio calculations. We showed that the heat capacity ( ) of CH3SiH2SiH(g) from the CHEMKIN

database was highly erroneous and unphysical, and that the use of this data could lead to the conclusion that this particular species was the main Si growth contributor at atmospheric pressure. Instead, when using the results from ab initio calculations we predicted that there was no significant growth contribution from this particular organosilicon species based on the kinetic model being used. However, whether it is possible to generalize and conclude that the organosilicon does not take part in the growth at all depends directly on how correct and how well the kinetic model is in explaining the growth mechanism of SiC CVD.

Most kinetic models in the literature at the time were focused on the Si-limited growth of SiC CVD. Thus, the kinetics of the hydrocarbon chemistry was less well investigated. Moreover, experimental results performed by our group 132_,

showed that there were difference in the growth rates and surface roughness when different hydrocarbons were used as precursors. This led us to study the gas-phase kinetics of hydrocarbons presented in the literature. In this study, we selected two models from the literature, one from Allendorf and Kee 41_and

the other from de Persis et al. 133_{, as well as a model compiled by one of the}

authors. The simulation conditions were kept the same for all three and no depositions were allowed. We chose to compare the results when three differ-ent hydrocarbon precursors, CH4, C2H4, C3H8, were used. The studies showed

that under the same condition the three models predicted relatively different profiles of concentrations of the most abundantly produced hydrocarbon spe-cies and atomic H. Thus, it is clear that more improvement must be made in kinetic models before we can expect reliable and realistic results.

(42)

Researchers, due to lack of kinetic data, usually chose to adopt and modify the surface kinetic model presented by Allendorf and Kee 41_{. This model however}

was shown to predict C-limited growth which conflicted to the Si-limited growth observed in experiments. Kinetic studies performed by one of the au-thors using the gas phase and surface kinetics from Allendorf and Kee 41

showed that the C-limited mode was constantly predicted for at least up to the C/Si ratio of 2.0, with no transition to the Si-limited mode near C/Si =1 as observed in experiments. Thus, it is clear that not only the gas phase model but also the surface kinetics are needed in order to build a more reliable mod-el.

In Paper 2, we derived thermochemical properties of the halides and halohy-drides of Si and C species, SiHnXm and CHnXm, for X being F, Cl and Br.

These properties are needed for CFD modeling as well as for thermodynamic equilibrium calculation of halide-assisted SiC-CVD. We chose two composite methods, namely the Gaussian-4 (G4) 110_{and the Weizmann-1 theories}118,121

as modified by Barnes et al. 122_(W₁_{RO), for the calculations. Due to limits in}

the design of the method, W1RO was not applied to the study of Br-containing species. With these two methods, we derived the atomization energies, en-thalpies of formation and entropies of the molecules listed above. We showed that the experimental data from the databases are highly inconsistent to each other as well as to the theoretical data. For example, we observed discrepan-cies in the atomization energies as large as 15‒40 kJ/mol between the experi-mental and theoretical data for SiF, CHF, SiH3F and CH3F. On the other hand,

theoretical data from the literature and our calculations agreed well. Both methods worked nicely for the molecules in the set.

Paper 3 combined the growth studies of SiC-CVD, in which SiF4 was used as

the Si precursor, together with theoretical studies of thermodynamic equilibri-um calculations and surface kinetic ab initio calculations. The results from thermochemical calculations performed under the same condition as the ex-periments suggested that the main candidates for Si growth contributors are SiF2, SiF and SiHF. This was concluded from the concentration levels of these

species and whether they are capable of contributing to the growth at the rate observed experimentally. Secondly, we studied the change in the species con-centrations vs. temperatures as well as the reactivity of the molecules. Based on the experimental growth rate profiles, these two points narrowed down the set of candidates into the three species listed above.

(43)

The part which I was mainly concerned with was the ab initio calculations of the surface kinetic of the three species. The study was performed on the (0001) surface of 4H-SiC modeled by a cluster of Si13C13, shown in Fig. 4.1a. The

cluster was completely passivated by hydrogen atoms on all sided to preserve the SiC bulk symmetry. The choice of surface passivation was justified by the typical operating conditions of SiC-CVD, i.e. a temperature of ~ 1600 ℃ in a hydrogen rich atmosphere. We chose two types of surface sites on which the Si-species were assumed to adsorb. The first site was a fully coordinated C site and the second was a C site with the presence of one dangling bond. The reason that we chose C sites as the adsorption sites, instead of Si sites, was because we were more interested in the growth with correct stoichiometry. It should be noted that the (0001) surface of 4H-SiC is a Si face. Thus, these carbon sites were obtained by pre-attaching a hydrocarbon group on a Si face before adsorption occurred. The choice of surface C groups to act as the ad-sorption sites were made by choosing the simplest and smallest hydrocarbon group. Thus, a methyl group, CH3, was chosen to represent a fully-bonded C

site and a methylene group, CH2, was chosen for a C site with the presence of

one dangling bond.

We then used DFTs to determine for the Gibbs free energies of adsorption, , of the three Si-species on these two sites over a temperature range well covering the experimental conditions. The calculated at the temperature corresponding to the experimental condition were shown to be positive for all cases, except for the SiF adsorption on a methylene site. This suggested, with the exception of SiF on a methylene site, that the products of adsorptions were less stable than the reactants. The adsorptions on the methyl site were shown to occur with the presence of transition state barriers for all cases. These tran-sition barriers were shown to be relatively high and led to sticking coefficients which were too low to explain the experimental results. Thus we concluded that it was unlikely for the three Si-species to contribute to the growth by ad-sorbing on a fully bonded C site. On the other hand, the adsorptions on the methylene showed no transition barriers between the reactant and product states. Furthermore, from the negative we concluded that SiF was likely to contribute to the growth. The positive of SiHF and SiF2 made the two

species less likely to contribute in comparison to SiF. And since was less positive for SiHF in comparison to SiF2, we concluded that the chance of

SiHF to overcome the positive is higher than for SiF2.

Paper 4 was an extension of the study from Paper 3 where we included the adsorption of SiH, SiH2, SiX, SiHX, SiX2 for X being F, Cl and Br. Also, both

A Quantum Chemical Exploration of SiC Chemical Vapor Deposition