UFDC Home  myUFDC Home  Help 



Full Text  
QUANTUM CHEMISTRY CALCULATIONS OF THE REACTIONS OF GASEOUS OXYGEN ATOMS WITH CLEAN AND ADSORBATETERMINATED Si(100)(2x1) By PAULO EMILIO HERRERAMORALES A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2005 Copyright 2005 by Paulo Emilio HerreraMorales To my late grandmother, Catalina Caal Prado. ACKNOWLEDGMENTS I wish to express my sincere appreciation to Dr. Jason F. Weaver (my supervisory committee chair) for all his encouragement and advice throughout my Ph.D. program. I would like to extend my deepest gratitude to the members of my advisory committee (Dr. Gar B. Hoflund, Dr. Vaneica Y. Young, and Dr. Fan Ren) for their assistance, technical support, and guidance. I also thank Dr. David Micha, Dr. Erick Deumens, and Dr. Adrian Roitberg (from the Quantum Theory Project) and Dr. Helena HagelinWeaver (from the Chemical Engineering Department) for their valuable support with quantum chemistry matters. I would like to thank my colleagues at Universidad Rafael Landivar in Guatemala, for their support and encouragement. I would like to express my sincere appreciation and gratitude to my parents, Victor Manuel Herrera and Emma Morales de Herrera; and to my sister, Servanda Virginia Herrera, for supporting me and my parents for the duration of this work. Special thanks go to my fiancee Sonia Maria Garcia, who supported me more than once. She possesses admirable patience. I would like to thank the Quantum Theory Project of the University of Florida for the use of its computational facilities; and Brad Shumbera, for his excellent work in setting up the computer network for our research laboratory. I also want to recognize the other Weaver Research Group members, Alex A. Gerrard, Jau Jiun Chen, Sunil Devarajan and Heywood Kan, for their help with different issues through the years. TABLE OF CONTENTS page A C K N O W L ED G M EN T S ................................................ ........................................... iv LIST OF TABLES .............................................. viii LIST O F FIG U RE S .............. ......................... ........................... ....................... .. .. .... .x A B S T R A C T .................................................................................................... ........ .. x iii CHAPTER 1 IN T R O D U C T IO N .................................................................. .. ... .... ............... 1 1.1 Research Objective ..................... .. ........ ..... .......................... 1.2 T he Si(100)(2x l) Surface .................... ................................. ................. ............... 1 1.3 Chemistry of Atomic Oxygen on Si(100)(2xl).......................... ...............3...... 1.4 Important Siliconbased Materials for Microelectronics Industry.......................4... 2 THEORETICAL BACKGROUND OF QUANTUM CHEMISTRY C A L C U L A T IO N S ......................................................................... .. . .. ...............8 2 .1 Introdu action ...................... .... ... ..................... ..................... ............ .. .8 2.2 Energy Functionals and Early Density Functional Theory (DFT) .................... 10 2.3 HohenbergKohnSham Density Functional Theory (KSDFT) .......................14 2.3.1 HohenbergKohn Lemmas and ExchangeCorrelation Functional D e fin itio n .............. .... ........ ..................................................................... 1 5 2.3.2 KohnSham Selfconsistent Field Method .............................................17 2.3.3 Exchange and Correlation Functional Approximations ..............................18 2.3.3.1 Local density approximation (LDA)........................................ 19 2.3.3.2 D ensity gradient corrections...................................... ................ 20 2 .3.3.3 H ybrid m ethods ......................................................... ................ 22 2.3.3.4 E m pirical D FT m ethods ................ .......... .......... ...................... 24 2.3.4 KohnSham Density Functional Theory Computational Chemistry ...........25 2.4 B asis Sets ................. ........................................ ............... 27 2.5 Geometry Optimization and Transition State Search Calculations...................31 2.5.1 Geometry Optimization of Energy Minima ..........................................33 2.5.2 T transition State Searches ...................................................... ................ 38 3 NITROGEN ATOM ABSTRACTION FROM Si(100)(2xl)...............................41 3 .1 In tro d u ctio n ........................................................................................................... 4 1 3.2 C om putational A pproach.................................... ....................... ................ 43 3 .3 R e su lts...................... ... .. .............................. ........................................... .. 4 5 3.3.1 Bonding Configurations of a Nitrogen Atom on Si(100)(2xl)..............45 3.3.2 Nitrogen Abstraction by a GasPhase Oxygen Atom...............................46 3.3.3 Abstraction of N Adsorbed at the Dangling Bond [R(ad)] .....................48 3.3.4 Abstraction of the Nitrogen Bonded Across the Dimer [R(db)] .............51 3.3.5 Abstraction of the Nitrogen Bonded at a Backbond [R(bb)]...................54 3.3.6 Abstraction from the N=Si3 Structure [R(sat)]..................................... 59 3 .4 D iscu ssio n ............................................................................................................. 6 3 4 CHEMISTRY ON Si(100)(2xl) DURING EARLY STAGES OF OXIDATION W ITH O (3P) .................................................................................. . .................66 4 .1 In tro d u ctio n ...........................................................................................................6 6 4 .2 T heoretical A approach .......................................... ......................... ................ 69 4 .3 R esu lts an d D iscu ssion ....................................................................... ............... 70 4.3.1 Structures with One Adsorbed Oxygen Atom (OiSi9H12) .....................77 4.3.2 Structures with Two Adsorbed Oxygen Atoms (02Si9H12)...................80 4.3.3 Structures with Three Adsorbed Oxygen Atoms (03Si9H12) .................86 4 .3 .4 O x y g en In sertio n ........................................................................................9 1 4.3.4.1 Therm odynam ic considerations ................................ ................ 91 4.3.4.2 K inetic considerations ............................................... ................ 94 5 ATOMIC OXYGEN INSERTION INTO ETHYLENE AND ACETYLENETERMINATED Si(100)(2xl).............................................100 5 .1 In tro du ctio n ......................................................................................................... 10 0 5.2 Theoretical A approach ................. ............................................................ 102 5.3 R results and D iscu ssion ........................................................... ...................... 104 5.3.1 Relative Energies of Oxidized EthyleneCovered Si(100) Clusters.........110 5.3.2 Relative Energies of Oxidized Acetyleneterminated Surfaces .............115 5.3.3 Oxygen Insertion M echanism s....... .......... ....................................... 119 5.3.3.1 Therm dynamic considerations ........................ ................... 119 5.3.3.2 K inetic considerations ....... ........... ....................................... 120 6 SUMMARY AND CONCLUSIONS.......... ..........................124 6.1 Nitrogen Atom Abstraction from Si(100)(2xl) by Gaseous Atomic Oxygen ..124 6.2 Initial Step of Si(100)(2xl) Oxidation by Gaseous Atomic Oxygen..............125 6.3 Oxidation of C2H2 and C2H4covered Si(100)(2xl) by Gaseous Atomic O x y g e n ............................................................................................................ . 12 7 APPENDIX A QUANTUM CHEMISTRY SOFTWARE ............................................................ 128 A I1 F ile F orm at .xyz .... ... ......................................... ....................... . .......... 133 A .2 F ile F orm at .zm t ............................................. .. ..................... ... ......... ............ 33 A.3 HyperChem Files (Format .hin)................ ............................................ 136 A.4 Gaussian03 Files (Formats .inp, .log and .chk) ...................... ................... 138 B EXAMPLE OF GAUSSIAN03 OUTPUT FILE ........................... ..................... 141 C CALCULATION PROTOCOL FOR INITIAL GEOMETRY FOR GEOMETRY OPTIMIZATIONS ............................. .......... ........ ............... 157 D CALCULATION PROTOCOL FOR INITIAL GEOMETRY FOR TRANSITION STATE SEARCHES ................................... .............................. 159 LIST O F R EFEREN CE S ... ................................................................... ............... 166 BIOGRAPH ICAL SKETCH .................. .............................................................. 173 LIST OF TABLES Table page 11 Selected properties of silicon dioxide ................................................... ...............6... 21 Approximations made to correct ThomasFermiDFT before the development of KohnSham DFT ...................................... ........ ................... 13 22 Variables used to express the approximate exchangecorrelation functionals E x c [p (r )] ............................................................................................................... ... 1 8 41 Suboxide penalty energies for various silicon oxidation states...............................75 42 Penalty energies of O iSi9H 12 isom ers ................................................ ................ 78 43 Penalty energies of 0 2Si9H 12 isom ers ................................................ ................ 81 44 Penalty energies of 0 3Si9H 12 isom ers .....................................................................87 45 Energy barriers for migration of oxygen atoms in OSi9H12Od structures from a dangling bond site to a backbond or dimer bond site in the spintriplet state..........96 51 Ei(SiSi) structure and energy of adsorption compared to those of reported dimerized C2H4terminated Si(100) structures.......................... ...................107 52 Ai(SiSi) structure and energy of adsorption compared to those of reported dimerized C2H2terminated Si(100) structures.......................... .................. 107 53 Calculated energy of formation of the different oxidized bridges that form when an oxygen atom inserts into the C2H2 and C2H4terminated Si9H12 clusters ...........109 54 Penalty energies of spinsinglet OC2H4Si9H12 isomers................................. 110 55 Penalty energies of spintriplet OC2H4Si9H12 isomers ..................................... 111 56 Penalty energies of spinsinglet OC2H2Si9H12 isomers..................................115 57 Penalty energies of spintriplet OC2H2Si9H12 isomers ..................................116 Ai Quantun chemistry programs used in this work......................... ...................132 A2 File formats used for proper handling of the coordinates systems and results in our geometry optimization calculations ....... ....... ...................... 132 A3 Gaussian03 input files sections ....... ........ .......... ...................... 138 C1 Calculation protocol used to determine the best initial geometry for the geometry optimization, transition state search and frequency calculations .........................157 LIST OF FIGURES Figure page 11 Si(100) 2x1 surface show ing the dim er row s. ............................... ..................... 2 12 Metaloxidesemiconductor fieldeffect transistor (MOSFET). ............................. 6 21 Flowchart of the KohnSham SCF procedure ............................... ..................... 26 22 Tw odim ensional potential energy surface ......................................... ................ 31 23 Flowcharts for a quasiNewton algorithm for geometry optimization.................. 36 31 The Si9H12 cluster used in our UB3LYP calculations.........................................47 32 Structural inform ation of N Si9H 12 clusters .................................. ..................... 48 33 Reaction pathway for the nitrogen abstraction from R(ad)................................ 49 34 Critical point structures of the nitrogen abstraction from R(ad)........................... 50 35 Reaction pathway of nitrogen abstraction from R(db) .............................................52 36 Critical point structures of the nitrogen abstraction from R(db)...........................53 37 Reaction pathways for the abstraction of the nitrogen atom inserted in a SiSi b a c k b o n d ................................................................................................................. 5 6 38 Molecular precursors formed after Ochemisorption onto R(bb) .........................57 39 Structures formed during nitrogen abstraction from MP(bbt) .............................. 58 310 Structures involved in the nitrogen abstraction from MP(bbs)............................. 59 311 Reaction pathways for Nabstraction from (and Oatom adsorption on) R(sat)......61 312 Structures formed during the direct nitrogen abstraction from R(sat) ..................62 313 Structures involved in the Oatom chemisorption on R(sat)............................... 62 41 Structural information and highest occupied molecular orbital plot of clean Si9H12 clu sters .................................................................................................... ........ .. 7 4 42 R elative energies of O 1Si9H 12 isom ers............................................... ................ 77 43 Structural characteristics of Si9H12 clusters with one oxygen atom......................79 44 R elative energies of 0 2Si9H 12 isom ers............................................... ................ 83 45 Structural information of dangling bond isomers with two oxygen atoms adsorbed o n S i(10 0 ) .............................................................................. .. ............... 8 4 46 Structural inform ation of 02Si9H 12 isom ers....................................... ................ 85 47 Relative energies of optimized structures for three oxygen atom incorporation into S i( 1 0 0 ) .................................................................................................................. ... 8 8 48 Structural information of the dangling bond 03Si9H12 isomer DB3(OSiOSiOd) ...89 49 Structural characteristics 03Si9H 12 isom ers ....................................... ................ 90 410 Qualitative representation of a minimum energy path for oxygen atom incorporation on Si(100)(2xl) based on the relative energy of the structures........ 93 411 Structural information of transition state structures for oxygen migration in 0 1Si9H 12 structures.. .................................................................................... 95 412 Structural information of transition state structures for oxygen insertion in 02Si9H 12 structures.. .................................................................................... 98 413 Structural information of transition state structure TS3(OSiOSiOdOSiOSiO) ......98 414 Proposed preferred path for initial steps of oxidation of Si(100)(2xl) by oxygen ato m s ...................................................................................................... ....... .. 9 9 51 Structural information of optimized clusters for olefincovered Si(100)(2xl).....106 52 Relative energies of spinsinglet OC2H4Si9H12 isomers....................................111 53 Relative energies of spintriplet OC2H4Si9H12 isomers..................................112 54 Structural information of OC2H4 Si9H12 isomers ..................... ...................114 55 Structural information of OC2H2Si9H12 isomers ...................... ...................117 56 Relative energies of spinsinglet OC2H2Si9H12 structures.................................118 57 Relative energies of spintriplet OC2H2Si9H12 isomers..................................119 58 Structural information of clusters involved in oxygen insertion pathways on C2H2 and C2H4terminated Si(100) surfaces ........................................121 59 Potential energy surfaces for insertion of an oxygen atom on an C2H2 covered Si(100)(2x l) surface ............. ............... ................................................ 123 510 Potential energy surfaces for insertion of an oxygen atom on an C2H4 covered Si(100)(2x l) surface ............. ............... ................................................ 123 A i H yperChem 7 interface .................. ............................................................ 129 A2 Molden visualization interface and main command screens................................130 A3 gOpenM ol visualization interface ....... ...... ...... ..................... 131 A4 Example of a .xyz file describing a Si9H1203 cluster.................. ...................134 A5 ZM atrix file describing a Si9H12+xCy cluster ....... ... ................. ................... 135 A6 HyperChem 7 files describing a Si9H1203 cluster....................... ...................137 A7 Gaussian03 file types ............................................................. 140 D1 Input file for a transition state search that uses the OPT(QST2) option ............. 160 D2 Diagram of the structures of reactant and product used to search for a transition state structure using OPT(QST2) or OPT(QST3) in Gaussian03 ........................161 D3 Zmatrix (.zmt) file generated by M olden ........................................ ............... 162 D4 Gaussian03 input file for a transition state search where the product structure had to be reordered using M olden.................................... ...................... ............... 163 D5 Process used to generate an initial guess for the transition state search using OPT(Q ST3) in G aussian03 .......................................................... 164 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy QUANTUM CHEMISTRY CALCULATIONS OF THE REACTIONS OF GASEOUS OXYGEN ATOMS WITH CLEAN AND ADSORBATETERMINATED Si(100)(2xl) By Paulo Emilio HerreraMorales May 2005 Chair: Jason F. Weaver Major Department: Chemical Engineering Reactions of gasphase radicals at solid surfaces are fundamental to the plasmaassisted processing of semiconductor materials. In addition to adsorbing efficiently, radicals incident from the gasphase can also stimulate several types of elementary processes before thermally accommodating to the surface. Reactions that occur under such conditions may be classified as nonthermal; and examples include atom insertion, directatom abstraction and collisioninduced reaction and desorption. Indeed, nonthermal surface reactions play a critical role in determining the enhanced surface reactivity afforded by plasma processing. Advancing the fundamental understanding of radicalsurface reactions is therefore of considerable importance to improving control in plasmaassisted materials processing. Our study used quantum chemical calculations to investigate the interactions of gasphase oxygen atoms [0(3P)] with clean and adsorbatemodified Si(100)(2xl) surfaces. We carefully studied reaction pathways for the oxygen insertion, adsorbate abstraction and adsorbate migration of these species on reaction with 0(3P) using density functional theory (DFT) and silicon clusters. Knowledge of these reaction pathways helps determine the viability of certain reactions on Si(100)(2xl) by gasphase oxygen atoms in plasmaassisted material processing in general. Another goal of our study was to better understand the chemistry of certain ultrathin siliconbased films. Silicon dioxide has been at the center of the microelectronics industry given its multiple electrical properties. However, as the size of the devices keeps shrinking according to Moore's law predictions of number of transistors, the demands imposed on this material necessitate a detailed understanding of its chemistry at the molecular level urging researchers to look for several alternatives. Materials such as silicon oxynitride (SiOxNy) and silicon oxycarbide (SiOxCy) have been suggested as possible silicon dioxide substitutes as gate dielectric and interface dielectric, respectively. Silicon oxynitride can be obtained by either inserting nitrogen atoms into silicon dioxide or oxygen atoms into a silicon nitride film, resulting in ultrathin layers that enhance properties of pure silicon dioxide as a gate dielectric. Silicon oxycarbide is one of the products that evolved as an alternative for silicon dioxide as a result of the recent trend of organic functionalization of the silicon surfaces (a process by which hydrocarbon molecules are chemisorbed on top of the silicon surfaces to combine the semiconductor properties of silicon with the organic functionality of carbon). Thus we used quantum chemistry calculations to study the initial steps of oxidation of the clean Si(100)(2xl) surface, of the abstraction of nitrogen atoms, and for insertion of gaseous oxygen atoms on acetylene and ethyleneterminated Si(100). CHAPTER 1 INTRODUCTION 1.1 Research Objective Our main objective was to advance fundamental understanding of the reactions of gaseous oxygen atoms with clean and adsorbatemodified Si(100)(2xl) surfaces, motivated by an interest in gaining insights into the possible reaction pathways during the initial steps of oxidation. Our work was also motivated by the potential benefits of fabricating ultrathin siliconbased films and by the need to precisely control the properties of such films. Density functional theory (DFT) quantum chemistry calculations were run to investigate the possible reactions. This chapter provides a brief introduction to the Si(100)(2xl) surface, an overview of the current state of knowledge of the oxidation process of this surface, and a brief description of the industrially important silicon dioxide, silicon oxynitride, and silicon oxycarbide, semiconductor materials whose production can be deeply impacted by knowledge of oxygenatom chemistry at the molecular level. Density functional theory and the background of the quantum chemistry calculations are explained in detail in Chapter 2. 1.2 The Si(100)(2xl) Surface Because of its numerous industrial applications and model character for semiconductor surface science, the industrially important Si(100) surface has been widely studied, dating back to the early work of Schlier and Farnsworth [1]. It is known that when first formed, the surface reconstructs by forming silicon dimmers, and that the surface adopts a (2xl) orientation (where 2x1 designates the new periodicity of the surface atoms). The surface forms these dimers so it can reduce the number of dangling bonds per surface atom, from two in the bulkterminated (lxl) surface to one in the reconstructed (2xl), and also because their formation lowers the surface energy by about 1.0 eV [2] (Figure 11). Figure 11. Si(100) (2xl) surface showing the dimer rows. It is generally accepted [38] that the dimers are asymmetric (buckled) and consist of an sp2like bonded down atom, that moves closer to the plane of its three nearest neighbors; and an up atom that moves away from the plane of its neighbors, and possesses an slike dangling bond. This rehybridization process from the bulk sp3 to the sp2like configuration of the silicon dimer atoms is accompanied by a charge transfer from the down to the up atom. To minimize this electrostatic effect (and to relieve local stress) the direction of buckling alternates within the dimer rows. The silicon dimer bonds in Si(100)(2xl) are similar to a C=C double bond [9]. The organic double bond consists of two types of interactions: a a bond with symmetry around the axis connecting the two atoms, and a sr bond with a nodal plane along the axis. Similar interactions exist in the SiSi dimer bonds, with the difference that the ar interaction is sufficiently weak that the dimer is not held in a symmetric configuration, thus adopting the buckled asymmetry that constitutes its more recognizable structural feature [10] (Figure 11). 1.3 Chemistry of Atomic Oxygen on Si(100)(2xl) Among the most reliable sources of gasphase radicals such as groundstate atomic oxygen [0(3P)] are plasmas, which typically contain a variety of highly reactive, unstable species (radicals and charged particles) that, in addition to adsorbing efficiently, can also stimulate several types of elementary processes before thermally accommodating to the surface. These elementary processes include atom insertion, directatom abstraction, and collisioninduced reaction and desorption. It is generally accepted that atom abstraction can occur in several ways, ranging from the EleyRideal mechanism (ER) of direct abstraction (in which an atom is abstracted directly from the surface in a single collision with an incident species, without reaching thermal equilibrium with the surface) [ 11]; and the LangmuirHinshelwood mechanism (LH) of collisioninduced reaction and desorption (in which the incident species reaches thermal equilibrium with the surface before reacting with the adsorbate and eventually evolving into the gasphase) [12]. An intermediate reaction mechanism for atom abstraction is the the hot atom mechanism (HA) in which the incident species experiences multiple collisions with the surface, but does not fully thermalize before the reactive encounter [13]. To study the aforementioned reaction mechanisms, it is necessary to have a convenient source for all types of radicals. Unfortunately, such a source exists only for hydrogen radicals, which can be produced efficiently by thermal dissociation of molecular H2 over a hot filament [14]. Thus, abstraction reactions on silicon have been widely investigated only for hydrogen atoms [1516]; and several kinetic models for the H/Si(100) system have been proposed, most of them favoring the ER mechanism [1718], the hot atom precursor scenario [1921], or both [2225]. Other systems, however (such as those of reactions of oxygen atoms on metal and semiconductor surfaces) have not been studied as often, despite their industrial and technological importance. Nevertheless, a few studies on oxidation of silicon surfaces by gaseous oxygen atoms have been reported. Engel et al. [2631] performed a detailed kinetic study of oxidation of clean Si(100) by gaseous oxygenatoms and reported that Oatom oxidation is facile compared with that caused by molecular oxygen, and that several monolayers of oxide can be formed efficiently via the direct insertion of O atoms into nearsurface bonds. Orellana et al. [32] studying the effect of the spinstate of the surface on the oxidation of Si(100) by 02, found that the oxidation proceeds with a high probability along the triplet potential energy surface (PES), and concluded that the spin state of the surface may be an important factor in determining the quality of the Si/Si02 interface. Despite much study, the reaction mechanism for abstracting hydrogen atoms from Si(100)(2xl) is still unclear. So the oxygenatom chemistry on Si(100) is far from understood. Thus, we used quantum chemical calculations to investigate interactions of the Oatoms with clean Si(100), spindoublet Ncovered Si(100), and spinsinglet and spintriplet C2H2 and C2H4terminated Si(100) surfaces. We focused on determining the reaction pathways leading to oxygen insertion or adsorbate direct abstraction from these surfaces, and on the effects that the spinstate has on the processes. This information may be important for the microscopic understanding of the initial steps of the silicon oxidation process by gaseous oxygen atoms. It may also be relevant for developing improved techniques for the fabrication of ultrathin layers of siliconbased industrial materials. 1.4 Important Siliconbased Materials for Microelectronics Industry Semiconductors are extremely useful for electronic purposes, because they can carry an electric current by electron propagation or hole propagation, and because this current is generally unidirectional, and the amount of current may be influenced by an external agent [33]. The most widely used insulator is silicon dioxide (SiO2), which is also one of the most commonly encountered substances in daily life. Modern integrated circuit industry was made possible by the unique properties of silicon dioxide. It is the only native oxide of silicon that is stable in water and at elevated temperatures, is an excellent electrical insulator, is a mask to common diffusing species, and is capable of forming a nearly perfect electrical interface with its substrate (Table 11). SiO2 has strong, directional covalent bonds, and has a welldefined local structure in which four oxygen atoms are arrayed at the corners of a tetrahedron around a central silicon atom. The SiOSi bond angles are essentially the tetrahedral angle (109.40), and the SiO distance is 1.61 A with very little variation. It is these siloxanebridge bonds between silicon atoms that give SiO2 many of its unique properties [33]. The Si/SiO2 interface is one of the key components of the commonly used metaloxidesemiconductor fieldeffect transistor (MOSFET), the building block of integrated circuits (Figure 12). This type of transistor continues to be the predominant device in ultralarge scale integrated circuits (USIC) because it is simple to scale down [34]. MOSFETs consists of a source and a drain (regions of doped silicon), a gate dielectric (an extremely thin silicon dioxide layer), and a gate electrode (a layer of polycrystalline silicon), rendered conductive by heavy doping, to bias the gate dielectric. Other MOSFET elements include thinmetal interconnect layers (made of Cu) to connect the transistor electrically to other parts of the circuit, and SiO2 dielectric layers to provide electrical isolation between the Cu interconnects and other devices. Basically, MOSFETs are switches that allow current to flow from the source to the drain only when the gate electrode supplies the appropriate bias voltage through the gate dielectric [34]. Table 11. Selected properties of silicon dioxide. Oxide native to silicon Low interfacial defect density (1010eVcm2) Melting point = 1713C Energy gap = 9 eV Resistivity 1015 cm Dielectric strength ~ 107 V/cm Dielectric constant relative to air (K) =3.9 Field Oxide Gate / C, / Silicon Source Drain Gate dielectric Figure 12. Metaloxidesemiconductor fieldeffect transistor (MOSFET). Silicon dioxide has a relatively lowdielectric constant (K = 3.9) [33]. Thus, since high gate dielectric capacitance is necessary to produce the required drive currents from submicron devices, and because capacitance is inversely proportional to gate dielectric thickness, the SiO2 layers have been scaled to extremely thin dimensions where the Si/Si02 interface becomes a more critical (and limiting) part of the gate dielectric. As a result, 1.0 nm Si02 layers are commonly found in modern materials. Such a thin silicon dioxide layer is mostly interface, and it contains about five layers of Si atoms (at least two of which reside at the interface). This has given rise to several electrical and performance problems (including impurity penetration through the dioxide, and a lack of convenient insulating properties) that make it necessary to look for alternate gate dielectric materials. Over the past several years, ultrathin silicon oxynitride (SiOxNy) films have been incorporated into metaloxidesemiconductor (MOS) devices as an alternative for pure SiO2 [33]. Incorporating nitrogen into SiO2 is a relatively simple method for fabricating silicon oxynitride films, since the nitrogen atoms tend to aggregate at the Si/SiO2 interface, and results in dielectric layers that enhance resistance to gate current leakage and inhibit boron penetration into the dielectric [34]. In recent years, surface functionalization also has been intensively investigated [3536]. This is the process of depositing layers of organic molecules on semiconductor surfaces to impart some property of the organic materials to the semiconductor device [3536]. The organic molecules might be designed to serve in place of gate oxides in metaloxide semiconductor fieldeffect transistor (MOSFET) devices, where the higher wire resistance of smaller metal lines and the crosstalk of closely spaced metal increase the interconnect resistance and capacitance product delay. This requires a lowdielectric constant (low k) material as the interlayer dielectric, and lowresistance conductors (such as aluminum). Until recently, silicon dioxide was the material of choice, but the decrease of the device dimensions and the resulting change of aluminum for copper have necessitated the search for alternative materials. Silicon oxycarbide (SiCxOy), a low k hydrid between organic and inorganic materials obtained by organic functionalization, is one of the most favorable candidates for interlayer dielectric in MOSFETs [3739]. CHAPTER 2 THEORETICAL BACKGROUND OF QUANTUM CHEMISTRY CALCULATIONS 2.1 Introduction From a chemistry and physics viewpoint three types of systems exist [40]: Systems that possess a very small number of particles where the particles do not interact and that can be studied classically and exactly, systems that have many particles and that require statistical methods to give a nearly exact solution of the properties under study, and systems where only a few particles exist, but where the particles interact considerably with one another. The latter cannot be solved exactly (neither classically nor by statistical modeling) so a different strategy is needed. Electron interactions with other electrons and nuclei within a molecule are an example of these complicated systems for which classical and statistical approaches are insufficient. If one considers that the molecule has N electrons with electronic coordinates ri, the total electronic energy of the system is obtained by solving the Schrodinger equation where H is the Hamiltonian operator, expressed as in Equation 21. 1 n n n 1 H= V + V(r) + (21) 2 I I j r, rj The three terms on the righthand side of Equation 21 correspond to the kinetic energy, the classical Coulomb electrostatic potential between nuclei and electrons, and the electronelectron interactions, respectively. From these, the electronelectron interaction term highly complicates the calculation, given that it describes the highly coupled motions of the electrons. If one attempts the classical approach, one can ignore the electronelectron interactions completely. This decouples the electronic motions and factorizes the problem into N completely independent problems, one for each electron. However, this approach loses essential elements of the physics of the system. If the statistical approach is attempted there is no need for intuitive simplifications, but a large supercomputer is needed to solve the Schrodinger equation directly. This requires due care in each calculation step, and has the drawback that the result is difficult to interpret and gives little physical insight. Thus, a molecular system is one of those systems for which neither classical nor statistical methods are sufficient. As indicated, these systems require a different approach, and a hybrid between the classical and statistical methods seems the most reasonable. The electronelectron interaction term in Equation 21 is replaced by a term that can be solved statistically or computationally, while the first two terms are solved classically. The available methods give a loss of few of the finer details of the system, but they factorize the original problem into N independent oneelectron problems that can be solved iteratively (given that usually the oneelectron equations depend on the energy, while the energy depends on the oneelectron equations) [40]. The two morewidely used iterative computational chemistry methods are HartreeFock (HF) and KohnSham density functional theory (KohnSham DFT) [4144]. Both methods are based on the selfconsistentfield (SCF) concept and share many conceptual and computational features. However, they are essentially different because while HF is based on singleelectron orbitals {t (r)} and many electron wavefunctions (') constructed from them, KohnSham DFT is based on the electron density p(r) and the fictitious singleparticle orbitals {(p, (r)} [4245]. This chapter presents the mostcommon energy functionals (Section 2.2), followed by a detailed description of the KohnSham SCF scheme and its implementation in Gaussian03 (Section 2.3), an introduction to the basis sets (Section 2.4) and an explanation of the geometry optimization schemes used in our study (Section 2.5). Since our main focus is on DFT calculations, we present only the main characteristics of the HartreeFock method, and do not discuss the more complex and exact postHartreeFock methods (such as CI and MP2). 2.2 Energy Functionals and Early Density Functional Theory (DFT) It is useful to express the total electronic energy of a molecular system partitioned in energy functionals, as follows E = E +E+E + Ex + Ec (22) given that the computer programs available find the approximate solutions to the Schrodinger equation by calculating each one of the functionals of Equation 22 separately. A functional is a mathematical device that maps objects onto numbers, and can be viewed as a function whose argument is also a function. Terms on the righthand side correspond to the kinetic energy of the electrons, the Coulomb energy of electrons caused by their attraction to the nuclei, the Coulomb energy that the electrons would have in their own field if they moved independently and if each electron repelled itself, a correction for electron exchange, and a correction for electron correlation, respectively. The last two terms can be combined as Exc = Ex + Ec rendering a correction for the false assumptions involved in Ej(i.e., that electrons do not perturb one another at close approach and that their motions are independent). The exchange energy Exis the electron stabilization resulting from the Fermi correlation (i.e., the dependence of motion arising from the Pauli Principle); while Ec (the correlation energy) arises from the correlation between the motions of electrons with different spins. These five terms have entirely different magnitudes, with ET, EV and Ej constituting most of E. The trend is as follows: ET > Ev > Ej >> Ex >> Ec. Additionally (as seen later), the expressions for Ev and Ej are common for all SCF methods used in quantum chemistry; while the formulas for ET, Ex and Ec are the ones that differentiate those methods [40, 45]. The HartreeFock approximation (the first SCF method developed) is based on orbital functionals, which are welldefined procedures that take the orbitals of a system and return an energy. The early orbital functionals that were suggested are the Hartree kinetic functional, the selfinteractioncorrection functional and the Fock exchange functional. In 1928, Hartree [40, 45] presented a model in which the ith electron in an atom moves completely independent of the other in an orbital y, (r). Since all electrons are represented this way, the system is then a collection of uncorrelated electrons, and the total kinetic energy of the atom ET is the sum of the kinetic energies of the electrons. H " E7H Y f V, (r)V 2, (rr (23) Molecular orbital (MO) theory originated when Hartree's model was extended to molecules by delocalizing the { V., (r) } over several atoms. However, the model neither gives the exact ET (given that the electrons do not move independently) nor excludes the selfinteractions of electrons (i.e., the exchange energy), and an additional selfinteractioncorrection functional is required. A correction functional was presented in 1930 by Fock [40, 45], who indicated that the Hartree wavefunction violates the Pauli Exclusion Principle because it is not properly antisymmetric. He then solved the problem by rendering the wavefunction antisymmetric as follows EF30 = l Idrdr2 (24). 2 2 ri r2 At this point in the history of quantum chemistry, it was realized that it is cumbersome to work with wavefunctions, given that they have essentially no physical interpretation and have units of probability density to the onehalf power. An alternative method was suggested, based on an observable that should allow the construction of the electronic Hamiltonian. Since the Hamiltonian depends on the position and atomic number of nuclei and the total number of electrons N, the best observable for the alternative method is the electron density p, given that the total number of electrons Nis obtained from N = p(r)dr (25) and, since the nuclei correspond to point charges, the nuclei positions are local maxima of p. Only the relationship between the atomic number ZA and the electron density was not clear, but it was found that ZA is related to the spherically averaged density. With this solved, Eq. 22 for the total electronic energy can be expressed in terms of density functional E[p(r)] = E, [p(r)] + E, [p(r)] + Ej [p(r)] + Ex [p(r)] (26) The simplest approximation for the total electron energy was the approach of ThomasFermi [4647], who assumed that the system had classical behavior, and then reduced Equation 26 to Equation 27. E[p(r)] = E, [p(r)] + E, [p(r)] + E, [p(r)] (27) In Equation 27, the Coulomb interaction terms in the righthand side are given by Equation 28 and 29. nuclei Z E,[p(r)] J= p(rr (28) and E [r)] p(r pr dr2 (29) 2 r, r2 E, [p(r)] and Ej [p(r)] are extremely useful for modern SCF methods, whether based on HF or DFT theories. Equation 28 is exact and is used in all the available SCF methods. Equation 29 is an approximation that considers that the electrons move independently and that each electron experiences the field caused by all the electrons (including itself). Thomas and Fermi [4647] also derived a kinetic functional based on the uniform electron gas, orjellium, which is an idealized system in the limit of N and V oo of a system of N electrons in a cubical box of volume V, that has a uniformly distributed positive charge that renders the system neutral. Thomas and Fermi chose jellium because although it is a manyelectron system, it is completely defined by its electron density p [4647]. The functional is as shown in Equation 210. 3 6 2 2/3 5/3(r r E,4p(r)] = (60f a (210) The ThomasFermi kinetic energy functional is the only density functional that has an elegant mathematical derivation, but unfortunately, it is not accurate enough to be chemically useful. Also, it was the first DFT functional as such, because it showed that nonelectrostatic energy terms can be expressed in terms of the electron density. Table 21. Approximations made to correct ThomasFermiDFT before the development of KohnSham DFT Functional Author Year Comments Hole function: Slater 1930 The hole function h(ri; r2) [ 1 p( f(r,)p ()dr2 Wigner 1933 is centered on electron E, [p(r)] = 2 r, r21 and 1 and a function of + p(rl)h(rl; r2)r1dr2 Seitz electron 2 coordinates. 2 r, _r2 Hole exchange energy ofjellium: Bloch and 1930 When combined with the 3 1/ 3 3 Dirac work of Thomas and EXD3[p(r)]= p43 (r)dr Fermi, is known as the 2 4;rj ThomasFermiDirac approximation. Slater exchange hole function: Slater 1951 Completely ignores 9(31/3 correlation effects. The Efs[p(r)]= J43 (r)dr exchange hole function 0 8 at any position is calculated as a sphere of constant potential with a radius that depends on the magnitude of the electron density at that position. Xa: Gaspar 1954 a= 4 works better than 9a 3 )1/3 P ( r Slater 1971 both 1 and 2/3. Ex[P(r)]= (r)dr and Wood with a = 1 for a = 1 for Exs51 and 2/3 for ExD3o However, these ThomasFermi functionals (known as ThomasFermi DFT) are no longer widely used, because they predict that molecules are unstable with respect to the dissociation of their constituent atoms. Table 21 shows the approximations made to correct for problems with ThomasFermiDFT before KohnSham density functional theory was developed; of these approximations, Xa is the only one that still is used in certain scientific endeavours. 2.3 HohenbergKohnSham Density Functional Theory (KSDFT) The early DFT model presented in Section 2.2 did not have widespread use because it resulted in fairly large errors in molecular calculations, and because the theories were not rigorously founded. This changed in 1964, when Hohenberg and Kohn (HK) [41] presented two key theorems to the scientific community. Before we further describe the HK theorems, we must clarify that we have limited this description to the simplest systems, that is, N nonrelativistic, interacting electrons in a nonmagnetic state with Hamiltonian E[p(r)]= ET [p(r)] + E [op(r)] + Ej [p(r)] (211) where ET[p(r)] Ey[p(r)] Zv(rj) and Ej[p(r)] 1 2 2i.jri rj For mathematical reasons, we considered a broad class of Hamiltonians with electrons moving in an arbitrary external potential v(r), besides the physically relevant Coulomb potentials due to point nuclei. 2.3.1 HohenbergKohn Lemmas and ExchangeCorrelation Functional Definition The starting point of KSDFT is the rigorous, simpleexistence lemma of Hohenberg and Kohn (HK) [41] who considered that the specification of the groundstate density, p(r), determines the external potential v(r) uniquely (to within an additive constant C). Since p(r) also determines N by integration, N[p(r)] = pc (r)dr = N, it determines the full Hamiltonian H (and, implicitly, all the properties determined by H). Hohenberg and Kohn [41, 45] proved that the groundstate density determines the external potential via reduction adabsurdum. They assumed that there were two potentials, va(r) and vb(r), with groundstate wave functions Yo,a and TO,b, respectively. They also assumed that both wave functions give rise to the same density p (r). Unless [vb(r) va(r) = constant] To,a and TO,b cannot be equal since they satisfy different Schrodinger equations. If one denotes the Hamiltonian and groundstate energies associated with Yo,a and YO,b by Ha and Hb and Eo,a and Eo,b, then one has by the minimal property of the groundstate, Eo,b = so that EO,a By interchanging a and b, HK found in exactly the same way EO,b Finally, adding Equations 212 and 213 leads to the inconsistency Eo,a + Eo,b < Eo,a + Eo,b. (214) With the help of the existence lemma, HK also showed (in a second lemma) that the density obeys a variational (or minimal) principle. For a given v(r), they defined the following energy functional of p(r) Ev(r) [P(r)] = f v(r)p(r)dr + F[p(r)], (215) where F[p(r)] ('P[p(r)IET + Ej 'P[p(r)]). (216) One should notice that F is a functional of p(r), since the wave function itself is a functional of p(r). The variational principle is then given by the expression Ev(r) [p(r)]> Ev(r) [o (r)] E (217) where po(r) and E are the density and energy of the groundstate. The equality in Eq. 217 holds only if po(r) = p(r); in other words, the variational principle states that any calculated energy will be higher than that of the groundstate. Then, Hohenberg and Kohn [41, 45] extracted from F[p(r)] its largest and elementary contributions by writing Equation 218, F[p(r)]= ET [(r)]+ I p(r ) drdr' + EXC [(r)] (218) 2 r r' where E,[p(r)] is the kinetic energy of a noninteracting system of electrons with density p(r), and the next term is the classical expression for the interaction energy. Equation 218 is the KSDFT definition of Exc[p(r)]. 2.3.2 KohnSham Selfconsistent Field Method If Exc is ignored in Eq. 218, the physical content of the theory becomes identical to that of the Hartree approximation [44]. Kohn and Sham (KS) [42] noticed that and transformed the EulerLagrange equation associated with the stationarity of Ev[p(r)] into a new set of selfconsistent equations, the KohnSham (KS) equations. 1 V + v(r)+ f (dr'+Vxc(r) j (r)= 0, (219a) N N[p(r)] = (r2 (219b) j=1 and VXC = (r) (219c). These equations are analogous to HartreeFock equations, although they also include correlation effects. They must be solved selfconsistently, like the HartreeFock equations, calculating Vxc in each cycle from Eq. 219c, with an appropriate approximation for Exc[p(r)] (Section 2.3.3). Despite the appearance of simple, singleparticle orbitals, the KS equations are in principle exact provided that the exact Exc is used in Eq. 219c; that is, the only error in the theory is due to approximations of Exc[p(r)]. Once Equations 219 are selfconsistently solved, the groundstate energy is obtained by E 2 rr r E = rcj 2r r drdr' f VXC (r)p(r)dr + Exc [o(r)] (220) 1 where e, and p are the selfconsistent quantities. The individual eigenfunctions and eigenvalues, p, and e,, of the Equations 219 have no strict physical significance. 2.3.3 Exchange and Correlation Functional Approximations In principle, the KSDFT exchangecorrelation functional accounts for the classical and quantum mechanical electronelectron repulsion and corrects for the difference in kinetic energy between the fictitious noninteracting system and the real electron system. Table 22. Variables used to express the approximate exchangecorrelation functionals Exc[p(r)] Variable Equation Comments Energy density (cxc) Ex [p(r)]= p(r),c [p(r)]dr Function of the electron density, given by the sum of individual exchange and correlation contributions in the system. Effective radius [rs(r)] r 1/3 Used when the electron r (r)= density is expressed for 4;Tp(r) exactly one electron. This electron would be contained within a perfect sphere that has the same density throughout its center. Normalized spin (p(r) p (r) Used to express spin polarization () (r) r= p(r) densities at any position Sa p n in openshell systems. By where: a corresponds to convention, the spin spinup and fs to spindown densities are given by electrons. 1 P (r)= P(r)(+1) 2 and p (r)= p(r) pa(r). However, there is no known systematic way to achieve an arbitrarily high level of accuracy for this functional. Therefore, an approximate form for the correct exchangecorrelation functional must be selected so that the KohnSham differential equations for the orbitals can be obtained by minimization of the total energy in Eq. 220. These Ex [,o(r)] approximations can be broadly classified in four major groups: the local density approximation functionals, the gradientcorrected functionals, the hybrid functionals and semiempirical DFT. Moreover, all these approximate functionals can be expressed in terms of the three variables summarized in Table 22. 2.3.3.1 Local density approximation (LDA) The designation local density approximation originally referred to any DFT functional where the energy density at some position r was computed exclusively from the value of the electron density at that particular position (i.e., the local value of p). However, since the only functionals that follow this definition are those derived based on the uniform electron gas electron density, generally the name LDA is applied only to exchange and correlations functionals derived from jellium. In the special case when one needs to calculate the exchange energy density for a spin polarized system, LDA can be modified with Equation 221. ex [p(r), = [01=(r)]+ e [p(r)] [p(r)] (1+ )42( 1 4/3 2 (221) 9Q (31/3 where c4 p(r)= 9a 3 3 [r], 810 and s\x [p(r)] is analogous to so [p(r)]. It describes a uniform electron gas composed of electrons with the same spin. This expression is the local spin density approximation (LSDA) [40, 45]. Since there has been little success in finding analytical expressions for c [p(r)], most of the work available on correlation energy density calculations is approximated by numerical techniques. One of the most widely used functionals is that of Vosko, Wilk and Nusair (VWN) [48]. After Ceperly and Alder used Monte Carlo techniques in 1980 to compute the total energy for jellium of various densities and found their correlation energy by subtracting the exchange, VWN developed a spin polarized correlation local density functional analogous to the exchange functional presented in Eq. 221 by fitting these computational results. The VWN proposed expression is Equation 222. r 2b f(4c b2)12 In + 2b tan 1 (4c b 2 A rs +b(r)/2 c 4cb2/2 2r/2 +b (222) fc(r < >:2bln  ) (222) C 2 bxo r/2 2 2(b+2xo) tan 1(4c b2 1/2 x2 +bxo +c r+br +C (4cb22 /2 2 /2 +b where i = 0 or 1 (analogous to E [p0(r)] and Ex [p(r)]), and with different empirical constant sets for gc [p(r)] and cc [p(r)]. This correlation functional, known asEcjv[p(r)], is exact for jellium but not for molecular systems and it is a good example of how many empirically optimized constants are required in DFT to approximate the unknown exactExc [p(r)]. Furthermore, the semiempirical flavor of functionals such as Ecv [p(r)] is the main inspiration for the development of the semiempirical DFT methods (Section 2.3.3.4). 2.3.3.2 Density gradient corrections Because the LDA exchange functionals were derived from a HF density matrix constructed with the plane wave orbitals of the uniform electron gas, they typically underestimate the exchange energy by roughly 10 to 15 %, and, consequently, need a correction for the nonuniformity of atomic and molecular electron densities. An improvement would be the generation of exchange and correlation functionals that depend on the local value of the electron density, and on the extent to which the density is locally changing (i.e., the density gradient). These types of functionals are called gradientcorrected functionals, and the methodology by which they are formed is known as the generalized gradient approximation (GGA). Most of these gradientcorrected functionals are constructed by adding a correction to an LDA functional energy density expression, that is, Equation 223. 'c [p(r), Vp(r)] = WLD [p(r)] + A /cx(r) (223) where x(r) is the reduced density gradient given by Eq. 224, x(r)= Vp( (224), and X/C indicates that this could be either an exchange or a correlation gradientcorrected density functional). The reduced density gradient has small values in bonding regions, larger values in core regions and very large values in Rydberg regions of molecules. The most widely used GGA exchange functional is the one developed by Becke in 1988 [49], which combines LSDA exchange with a gradient correction, and yields accurate exchange energies for atoms and correct exchange energy density in Rydberg regions. In its simplest form, the functional is given by Equation 225, 2 E88 = ED30 p )4/3 X dr (225) x 'x VJ ) 1+ 6/fxt sinh x, where the semiempirical parameter fl has a value of 0.0042, obtained by fitting to the ExF30 (HartreeFock) exchange energies of six noble gas atoms (He to Rn). Alternative GGA exchange functionals are do not include empirically optimized parameters. With respect to correlation, E'j N[p(r)] predicts correlation energies that are too large, making it a poor starting point for a gradientcorrected functional. Instead, the most popular GGA correlation functional developed by Lee, Yang and Parr (LYP) [50] computes full correlation energies. LYP abandoned jellium in favor of the He atom, the simplest system with a nonvanishing correlation energy. Their approach was based on earlier work of Colle and Salvetti [51], and it was improved for spincompensated and spinpolarized versions by Miehlich et al. [52]. The expression reported by Miehlich et al. for this gradientcorrected correlation functional, denoted as a secondorder gradient expansion is Equation 226. Ej [p,Vp]= 1a 4 ) 1 1 Vv1 i+dp 226 ab {Cp [211/3CF((a)8/3 +(P)8/3)+ 47 7 2 (226) 2 18 9 [p P 2P2Vp2 2 (pa)2)Vp2+(2p2 (P)2)V] 2} where = exp( 1)p 11/3, 1 3 + dp 1 CF (3 2)23 and the 1+ dp13 1d+dp1'3 10 parameters a, b, c and d were obtained from the work of Colle and Salvetti [51]. The values for the parameters are a = 0.04918, b = 0.132, c = 0.2533 and d= 0.349. ECL [p(r)] yields better results than EwN'[p(r)] and is very important in DFT despite some theoretical deficiencies. 2.3.3.3 Hybrid methods Given that the basis of the KohnSham approximation is the noninteractive reference system, it was suggested that one can control the conversion from a non interacting to a real interacting system. With the HollmannFeynman theorem, the DFT exchangecorrelation energy can be computed as shown in Equation 227. Exc = (Aj (A) Y (A))dA (227) where 2 is the extent of interelectronic interaction, going from 0 (noninteraction system) to 1 (real system). Development of hybrid DFT techniques resulted in methods that are significantly more efficient than ab initio methods of comparable accuracy. The relatively good accuracy and computational efficiency of hybrid methods indicate that further improvements might lead to methods that are both highly accurate and relatively efficient. All these methods include HartreeFock exchange because a small exact exchange component is a necessary constituent of any exchangecorrelation approximation aiming for accurate molecular energetic. The methods use the conventional mixing method introduced by Becke [53], which expands the functional with respect to the electron density and its gradient according to Eq. 223 and adds adjustable coeffiencts, c, yB, (228) an approach that is appealing because its application to a new molecule requires only a single KohnSham calculation. Becke developed the B3LYP exchangecorrelation functional [53], which is defined by Equation 229. EB3LY E28 +Ev +Ej +(Ic1)E 0 +c1E3 +cAE88 +(1c3)EWN +c3EL (229) where the coefficients cl, c2 and c3 are semiempirical coefficients determined by a fit to experimental data. This method cannot be derived rigorously and, at least in its current status, is basically empirical. Nonetheless, it is widely used because its ability to predict the atomization energies of normal systems is very close to their exact value. Kang and Musgrave [54] proposed a hybrid DFT method in which the exchange functional is a mix of Slater exchange and exact exchange. The KMLYP energy functional is given by Eq. 230. ELYP EH28 + EV + E +(1 a)Ef30 +aE3 + (I b)E vwN + bEL (230) where a= 0.557 and b = 0.448 and the subscripts are the same as in Eq. 229. 2.3.3.4 Empirical DFT methods These methods increased the emphasis on empirical parameterization, arguing that by modestly increasing parameterization it may be possible to obtain good functionals. Additionally, empirical DFT methods (EDF) looked for an exchangecorrelation functional that was optimized for a relatively small basis set and questioned the need for including HartreeFock exchange to obtain good agreement with experiment. Hybrid methods suppose that a small exactexchange component is a necessary constituent of any exchangecorrelation approximation aiming for accurate molecular energetic, but this introduces nonlocal effects and consequent computational complications that EDF methods tried to eliminate. Empirical density functional theory 1 (EDF1) [55] developed a DFT functional by linearly combining several existing functionals. It was optimized to yield accurate thermochemistry when used with the 63 1+G* Pople basis set (Section 2.4). Among other variational schemes, it considered linear combinations of different functionals or linear combinations of the same functionals with different parameters. The component functionals used are those of Becke [49] and Lee, Yang and Parr [50], discussed above. EDF2 [56] was developed following the method used in EDF1 [55], in which a DFT functional is obtained by linearly combining several existing orbital and density functionals, with a special focus on giving accurate harmonic frequencies when used with the ccpVTZ [57] basis sets. The expression EDF2 is as shown in Equation 231. EXF2= blEF3 []+ b2ED3 []+b3E"88 [,]+b4EDFo[p]+ bE"N []+b6E [po]++bEDF l[p] (231) where b, = 0.1695, b2 = 0.2811, b3 = 0.6227, b4 = 0.0551, b5 = 0.3029, b6 = 0.5998, and b7 = 0.0053. 2.3.4 KohnSham Density Functional Theory Computational Chemistry Figure 22 summarizes the steps involved in a KSDFT geometry optimization. For the most part, they are the same as those used in HartreeFock theory. First, a basis set is chosen to construct the KS orbitals (Eq. 219) and then an initial estimate of the molecular geometry. After that, the overlap integrals and the kinetic energy and nuclear attraction integrals are computed. The latter two kinds of integrals are called oneelectron integrals in HF theory to distinguish them from the twoelectron Coulomb and exchange integrals. (However, in KS theory all integrals can in some sense be regarded as oneelectron integrals since every one reflects the interactions of each one electron with external potentials). To evaluate the remaining integrals, one must guess an initial density, and this density can be constructed as a matrix entirely equivalent to the P HF density matrix, which describes the degree to which individual basis functions contribute to the manyelectron wave function, and thus, how important the Coulomb and exchange integrals should be. The elements of P are then computed as indicated in Equation 232. occupied P,'0 = 2 aia7 (232) where the coefficients aA, and a, specify the contribution of the basis functions to the molecular orbital i and the factor of two appears because this is a restricted calculation, meaning that it is considering only single wave functions in which all orbitals are doubly occupied [45]. Pople and Nesbet [58] presented the equation for P for radicals and excited states, which is known as unrestricted HartreeFock and which has one wave function for each electron. /Choose basis set Figure 21. Flowchart of the KohnSham SCF procedure Once the guess density is determined one can construct Vxc (Eq. 219c) and evaluate the remaining integrals in each KS matrix element. New orbitals are determined from solution of the secular equation, the density is determined from those orbitals, and it is compared to the density from the preceding iteration. Once convergence of the SCF is achieved, the energy is computed by plugging the final density into Eq. 220. At this point, the calculation proceeds according to the geometry optimization algorithm (i.e., if the geometry does not correspond to the optimal point, a new structure will be found and the KS SCF process will be run again, until the optimum is reached [45]). The geometry optimization algorithm is described in Section 2.5, but first, we present a brief introduction to the basis sets used in this work. 2.4 Basis Sets As discussed, Equations 219 are formally similar to the HF SCF ones, except that the exchange nonlocal oneelectron operator is replaced by a local exchangecorrelation operator depending on the total electron density. Thus, the KS equations can be efficiently calculated and the single electron wave functions can be represented by several basis functions such as Slatertype orbitals. A basis set is then, a set of mathematical basis functions from which the wave function is constructed in a quantum chemistry calculation. Slatertype orbitals have a number of attractive features primarily associated with the degree to which they closely resemble hydrogenic atomic orbitals. They suffer, however, from a fairly significant limitation. There is no analytical solution available for the general fourindex integral given en Equation 233 when the basis functions are STOs. The requirement that such integrals be solved by numerical methods severely limits their utility in molecular systems of any significant size. 1 2 Alternatives to the use of STOs have been proposed. One of them indicates that for there to be an analytical solution of the general fourindex integral (Eq. 233) formed from such functions is required that the radial decay of the STOs be changed from e' to er2. That is, the AOlike functions are chosen to have the form of a Gaussian function. The general functional form of a normalized Gaussiantype orbital (GTO) in atomcentered Cartesian coordinates is as shown in Equation 234. (x,y,z;a,i,j,k) 2a3/4 (8a)i+j+ki!j!k! /2 yj zk e (x2 + +z2) (234) KL) (2i) (2j) (2k)! Smalls GTO sets can be used for a wider range of chemical problems but involve some loss of flexibility in the resulting molecular orbitals. The simplest level of basis is minimal and corresponds to one basis function per atomic orbital. The next level is splitvalence or valencemultipleCin which two basis functions are used for each valence atomic orbital. This second level is known to give a better description of the relative energies and of some geometrical features of molecules. Further improvement of a basis set requires addition of functions of higher angular quantum numbers (polarization functions) [59]. The basis functions are generally contracted (i.e., each basis function is a linear combination of a number of primitive Gaussian functions). The contracted GTOs used in DFT to represent the electronic orbital wave functions, often were originally developed for HartreeFock calculations. A considerable increase in computational efficiency can be achieved if the exponents of the Gaussian primitives are shared between different basis functions [60]. At the splitvalence level, this has been exploited by sharing primitive exponents among s and p functions for the valence functions. In particular, a series of basis was defined and designated KLMG by Pople and coworkers [6168] where K, L, and M are integers. Such a basis for a firstrow element (Li to Ne) consists of an stype innershell function with K Gaussians, an inner set of valence s and ptype functions with L Gaussians, and another outer sp set with M Gaussians. Both valence sets have shared exponents. For hydrogen, only two stype valence functions (with L and M Gaussians) are used. Among this basis sets are 321G [6768], 431G [64, 69], and 631G [63, 6566]. The original 431G and 631G splitvalence basis sets were obtained by optimizing all Gaussian exponents and contraction coefficients to give the lowest spinunrestricted HartreeFock energy for the atomic groundstate. The Pople basis sets have seen widespread use among the scientific community [45]. A problem with the calculations based on KLMG basis sets is that s and p functions centered on the atoms do not provide sufficient mathematical flexibility to adequately describe the wave function for some geometries. This occurs because the molecular orbitals require more mathematical flexibility than do the atoms. Because of the utility of AOlike GTOs, this flexibility is almost always added in the form of basis functions corresponding to one quantum number higher angular momentum than the valence orbitals. Thus, for a firstrow atom the most useful polarization functions are d GTOs, and for hydrogen, p GTOs. A variety of molecular properties prove to be sensitive to the presence of polarization functions; for example, d functions on secondrow atoms are absolutely required to make reasonable predictions for the geometries of molecules including such atoms in formally hypervalent bonding situations (e.g., siliconates). Because the total number of functions begins to grow rather quickly with the addition of polarization functions, early calculations typically made use of only a single set. Pople and coworkers introduced first a simple nomenclature scheme to indicate the presence of the polarization functions, the "*" (star). Thus, 631G* implies a set of d functions added to polarize the p functions in 631G. A second star implies p functions on H and He (e.g., 6311G**) [45]. However, realizing the tendency to use more than one set of polarization functions in modern calculations, their basis set nomenclature now typically includes an explicit enumeration of those functions instead of the star nomenclature. Thus, 631G(d) is preferred over 631G* because the former generalizes to allow names such as, for example, 631G(3d2fg,2pd), which implies heavy atoms polarized by three sets of d functions, two sets of f functions, and a set of g functions, and hydrogen atoms by two sets of p functions and one of d. Finally, we must indicate that the highest energy MOs of anions and highly excited electronic states tend to be much more spatially diffuse than the molecular orbitals described so far. When a basis set does not have the flexibility necessary to allow a weakly bound electron to localize far from the remaining density, significant errors in energies and other molecular properties can occur. To address this limitation, standard basis sets are often augmented with diffuse basis functions. In the Pople family of basis sets, the presence of diffuse functions is indicated by a "+" in the basis set name. Thus, 63 l+G(d) indicates that heavy atoms have been augmented with an additional one s and one set of p functions having small exponents. A second plus indicates the presence of diffuse s functions on H [e.g., 631 l1++G(2d,p)]. For the Pople basis sets, the exponents for diffuse functions were variationally optimized on the anionic oneheavyatom hydrides and are the same for 321G, 631G and 6311G [45]. In this work we use the 321G*, the doublep plus polarization 631G(d), and the diffuse triplep plus polarization 631 l1++G(2d,p) Pople basis sets. 2.5 Geometry Optimization and Transition State Search Calculations Before we discuss how to perform geometry optimizations of molecular systems, we must understand the concept of a potential energy surface (PES) since a potential energy surface describes the energy of a particular molecule as a function of its geometry, and a geometry optimization is a process to find minima on a potential energy surface. Finding the total PES is a very complicated task, mainly because molecules have many atoms and many coordinates that describe their geometry. A simplified threedimensional PES can be represented as a topographic surface with valleys and saddle points (Figure 22). saddle point Global maximum (transition state) local minum = ;:J mnumu (reactant) 'Tr Te :tAble product) Local minimum (alternate product) Figure 22. Twodimensional potential energy surface Generally, calculation of most potential energy surfaces is based on the BornOppenheimer approximation which specifies that since the electrons are much lighter than the nuclei, the electronic part of the wave function readjusts almost immediately to any nuclear motion. Thus, potential energy surfaces can be obtained by calculating the total electronic energy for all the possible fixed nuclear positions that the molecule may adopt, an approach that is satisfactory for groundstate systems. On a PES, a reaction path is the movement from a valley of reactants to a valley of products, and a reaction mechanism is the path across the PES. Then, the important points in the calculated PES are the minima and the saddle points, with the equilibrium geometries corresponding to minima and the saddle points to the highest points on a reaction path that requires the least energy to get from the reactants to the products. Both PES minima and saddle points are called stationary points (or critical points), because for them the gradient (first derivative) is zero, and in classical mechanics, negatives of the first derivatives of the potential energy surface are the forces on the atoms in the molecule, as expressed by Eq. 234. F = VEV (234) An energy minimum must satisfy two conditions [70]: its first derivative (i.e., gradient of the forces) must be zero (critical point). its second derivative matrix must be positive definite (i.e., all eigenvalues of the Hessian must be positive or all of the vibrational frequencies must be real). Because of the difficulty of finding a complete PES, most quantum chemistry geometry optimization methods are focused on finding equilibrium structures and transition states directly. In this work, we performed geometry optimizations of Si9H12+wOxNyCz (w,x,y,z = 0,1,2, etc.) clusters using the unrestricted hybrid density functional method B3LYP combined with the diffuse tripleL plus polarization 631 l1++G(2d,p) basis set to expand the molecular orbitals of chemically active atoms (the surface Si atoms, four Si atoms of the first subsurface layer, the carbon atoms, the nitrogen atoms, and the oxygen atoms), and the doublet plus polarization 631G(d) basis set to describe the remaining subsurface silicon atoms and terminating hydrogen atoms (Section 2.4). Usually, finding equilibrium geometries and transition structures requires application of unconstrained optimization on the PES. However, we imposed constraints on the hydrogen atoms that terminate broken SiSi bonds in the clusters, to mimic the constraints that the rest of the solid should impose on the surface dimer under study. 2.5.1 Geometry Optimization of Energy Minima There are three types of geometry optimization algorithms: Methods that use only the energy, methods that use the first derivatives of the potential energy surface with respect to geometric parameters (i.e., gradient methods), and methods that require second derivatives (i.e., Newton or NewtonRaphson methods). The methods that only use the energy are the most widely applicable, but the slowest to converge, while the second derivative methods converge very fast but, given that analytic second derivatives are not readily available and more costly than the gradient methods. Thus, the gradient algorithms are the preferred methods, and (if analytical gradients are not available) it is usually efficient to calculate them numerically. In this work we performed quantum chemistry calculations with Gaussian03 [71], which uses the gradient method known as the Berny algorithm [7273] (Figure 23) as the standard method for calculations of both geometry optimizations to a local minimum and transition state searches. Gradient methods, such as the Berny algorithm, approximate the potential energy surface as a quadratic function and calculate the energy with the expressions E(r)= Eo + gAr + 1Ar'HAr (235a) 2 g = go + HAr (235b) where Ar = r ro, g is the gradient vector, and H is the Hessian matrix. Since for an energy minimum g = 0, the minimum structure can be obtained by solving linear equations Eq. 236a and Eq. 236b. g = go + HAr = 0 (236a) Ar = H g0 (236b) known as the Newton step. These first derivatives can be found analytically or numerically. The rate of convergence of the Berny algorithm depends predominantly on six factors: * Initial geometry of the structure to be optimized * Coordinate system chosen. * Initial guess for the Hessian matrix. * Hessian matrix updating method. * Step size control. The initial geometry of the structure is a critical issue. It cannot be just any structure that resembles the molecule under study as it would imply that the optimization using a large basis set would require numerous unnecessary steps. Generally, molecular mechanics and semiempirical methods are used to refine a raw structure and obtain a much better initial guess. In this work, we also ran a DFT/UB3LYP geometry optimization with the 321G* Pople basis set after the semiempirical calculations to get as close as possible to the final minimum energy structure. Early geometry optimization used nonredundant internal coordinates (e.g., the Z matrix internal coordinates), but later it was shown that Cartesian coordinates and a combination of Cartesian and internal coordinates had some advantages for particular systems [74]. Cartesian coordinates are the simplest and give an unambiguous representation of the structure, but they are also strongly coupled (i.e., to change a bond length one must change the x, y, and z positions of several atoms). Alternatively, internal coordinates (based on bond lengths, valence angles and dihedral angles) avoid problems with rigid body rotation s and translations. In internal coordinates, the coupling is much smaller than with Cartesian coordinates, and as a result the Hessian is more diagonal. (Appendix C reviews the protocol used to determine the initial coordinates for our calculations). The rate of convergence of a geometry optimization depends also on the initial guess of the nuclear Hessian since a closer guess to the actual Hessian will translate into a faster convergence. Computationally inexpensive empirical methods to generate initial nuclear Hessian have been proposed by Schlegel [75] and Fischer and Alml6f [76], who also compared the effect of different coordinate systems on the geometry optimization. An initial empirical estimate of a diagonal Hessian can be quite satisfactory for redundant internal coordinates and can be readily transformed to other coordinate systems. However, Baker [7778] showed that minimizations with Cartesian and redundant internal coordinates were similar if a good initial estimate of the Hessian is used (i.e., a molecular mechanics or semiempirical Hessian instead of a unit Hessian). The Hessian from an optimization or a frequency calculation previously performed at a lower level of theory is often a very good initial guess for an optimization at a higher level. Some difficult optimizations may require an accurate initial Hessian computed analytically or numerically at the same level of theory as that used for the optimization. The default in Gaussian03 estimates bond stretching force constants from empirical rules, and obtains average anglebending and torsion constants from vibrational spectra or theoretical calculations, using smaller and less expensive basis sets to calculate the force constants. Finally, since we are interested in constrained structures (given the rigidity that the silicon bulk would impose on the surface dimer chemistry) it is important to know that any components of the gradient vector corresponding to frozen variables are set to zero or projected out, thereby eliminating their direct contribution to the next optimization step. choose coordinate syste " //Choose coordinate system/ Figure 23. Flowcharts for a quasiNewton algorithm for geometry optimization Proper update of the Hessian is essential for efficient optimizations, as the quality of the Hessian at each point is critical for the success of the optimization. However, the calculation of the exact Hessian at each point would make the process lengthy and costly. Instead, the Hessian is adjusted in the quadratic approximation represented by Eq. 236 so that it fits the gradient gi at the current point ri and the gradient gi1 at the previous point. This leads to Equation 237. H,Ar = Ag (237) In Eq. 237 Ar = r, r, and Ag = g, g, There are numerous methods to update the Hessian, for example that of Broyden, Fletcher, Goldfab and Shanno (BFGS) [70], which can be written as expressed in Eq. 238. H, = H1 + AgAg'/Ar'Ag H,ArAr'H1 /Ar'H,1Ar (238) Eq. 238 is symmetric, positive definite, and minimizes the norm of the change in the Hessian. Schlegel et al. [74] indicate that most modem updating algorithms give similar results. Once there is a new dependable Hessian, a Newton step is taken on the model quadratic surface. The Hessian must be positive definite (i.e., all of its eigenvalues must be positive, for the step to be in the downhill direction). If the structure is far from the minimum (e.g., large gradients) or the potential energy surface is very flat (one or more small eigenvalues of the Hessian), then a simple Newton step may be too large, taking the molecule beyond the region where the model quadratic surface is valid. In this case, a shorter step must be taken. One can limit the step to be no larger than a trust radius which can be adjusted as the calculation proceeds, depending on whether the change in the predicted energy compares well with the actual calculated energy difference. In the Bemy algorithm of Gaussian03 the trust radius is updated using the method of Fletcher [79]. 2.5.2 Transition State Searches A transition structure is the highest point on the reaction path that requires the least energy to get from the reactant A to the product B. In other words, it is a stationary point that is an energy maximum in one direction and a minimum in all others. For a point to be considered to be a transition state structure the first derivatives must be zero and the energy must be a maximum along the reaction path connecting the valley of reactants with the valley of products on the potential energy surface. The transition state structure must be a critical point of index one (i.e., one of the eigenvalues of the Hessian matrix must be negative and all others must be positive). Energy, structure and vibrational information for transition state structures are obtained from the transition state theory. As the case for minima, transition structures can be found by geometry optimization, although some modifications must be incorporated into the procedure. Initial estimates for transition structures are more difficult to obtain than for equilibrium structures since transition state geometries vary much more than equilibrium geometries. Moreover, molecular mechanics generally cannot handle transition states involving the making and breaking of bonds. As a result, when searching for transition states, many quasiNewton methods need to start fairly close to the quadratic region of the transition state. Several techniques have been proposed to search for the initial structure of the transition state, including synchronous transit [80] and eigenvector following [81]. Gaussian03 combines synchronous transit and quasiNewton methods to find transition states [74, 80, 82]. This method requires two or three structures to start the optimization (one in the reactant valley, the second in the product valley, and an optional third structure as an initial guess for the TS geometry). If there is no third structure, the initial guess is obtained by interpolating between the reactants and products in redundant internal coordinates. The first few steps search for a maximum along an arc of circle connecting the reactantlike with the productlike structure, and a minimum in all other directions. The initial guess of the points along the reaction path is obtained by interpolating between two input structures; the energy and gradient are calculated at each point and an empirical estimate of the Hessian is obtained also at each point. The highest energy point on the path is chosen to be optimized to the closest TS, dividing the path into two downhill segments. In the remaining optimization steps, a quasi Newton based eigenvector following optimization is guided by the tangent to the arc of circle passing through the presumed TS and minima according to the implementation developed by Baker in 1986 [81]. As in the geometry optimizations, the use of a certain coordinate system can improve the transition state search considerably. Ayala and Schlegel [73] showed that the combined used of redundant internal coordinates and the tangent to the quadratic synchronous transit (QST) path improves considerably the transition state optimization, mainly because a better Hessian is obtained in the first few steps and improves the search direction. These redundant internal coordinates are based on the work of Pulay et al.[8384], who defined a natural internal coordinate system that minimizes the number of redundancies by using local pseudosymmetry coordinates about each atom and special coordinates for ring structures, and Peng et al. [74], who reduced the number of special cases of these coordinates by using a simpler set of internal coordinates composed of all bond lengths, bond angles and dihedral angles and applied it successfully for transitions state searches. Peng et al. [74] defined their coordinates as follows: First, the interatomic distances are examined to determine which atoms are bonded. Then, a bond angle bend coordinate is assigned for any two atoms bonded to the same third atom. Special attention must be given to linear bond angles; if the bond angle is greater than 1750, then two orthogonal linear angle bend coordinates are generated. Finally, a dihedral angle coordinate is assigned for each pair of atoms bonded to opposite ends of a bond. If one or both of the bond angles involved in a dihedral angle is linear, then the dihedral is omitted. In addition to the redundant internal coordinates generated automatically, extra stretch, bend and dihedral angle coordinates can be specified in the input. For regular transition state optimizations starting from one structure, the bonds being made or broken need to be specified in the input. Some difficult transition state searches may require an accurate initial Hessian computed analytically or numerically at the same level of theory as that used for the optimization instead of an updated one. Furthermore, Hessian update methods are suitable for finding minima for small and mediumsized molecules, but for difficult cases, such as some of the transition state structures investigated in this work, the Hessian has to be recalculated every few steps (or even every step) instead of being updated. This is equivalent to a Newton or NewtonRaphson algorithm. This is a very expensive method, and by default, Gaussian03 avoids it and updates the Hessian; however, its use has the advantage of obtaining an excellent description of the Hessian at each point of the calculation and also, since a vibrational frequency analysis is automatically done at the converged structure there is no need for an additional frequency job to determine the number of positive eigenvalues or zero point vibrational energies. The protocol that we used for transition state searches is presented in Appendix D. CHAPTER 3 NITROGEN ATOM ABSTRACTION FROM Si(100)(2xl) 3.1 Introduction The reactions of gasphase radicals at solid surfaces are fundamental to the plasmaassisted processing of semiconductor materials. In addition to adsorbing efficiently, radicals incident from the gasphase can also stimulate several types of elementary processes before thermally accommodating to the surface, including directatom abstraction and collisioninduced reaction and desorption. Directatom abstraction can occur by an EleyRideal mechanism in which an atom is abstracted from the surface in a single collision with an incident species [11] or by a hot atom mechanism in which the incident species experiences multiple collisions with the surface but does not fully thermalize before the reactive encounter [13]. Indeed, nonthermal surface reactions such as these play a critical role in determining the enhanced surface reactivity afforded by plasma processing. Advancing the fundamental understanding of radicalsurface reactions is therefore of considerable importance to improving control in plasmaassisted materials processing in addition to being of scientific interest. In the present study, we have used quantum chemical calculations to investigate the interactions of an oxygen atom with nitrogen atoms incorporated into the Si(100)(2xl) surface, focusing on pathways that lead to direct nitrogen abstraction from the surface by the formation of gaseous NO. This investigation is motivated firstly by an interest in determining the viability of direct nitrogen abstraction from Si(100) by gasphase oxygen atoms, and to gain insights into the possible pathways for this surface reaction. It is further motivated by the potential benefits that may be realized by incorporating ultrathin silicon oxynitride films into metaloxidesemiconductor (MOS) devices and the need to precisely control the properties of such films. Over the past several years, ultrathin silicon oxynitride (SiOxNy) films have been incorporated into MOS devices as an alternative for SiO2, which is no longer suitable as an insulator for MOS devices because of the decreasing dimensions of microelectronic devices imposed by Moore's law [33, 85]. Incorporating nitrogen into SiO2 is a relatively simple method for fabricating silicon oxynitride films, and results in dielectric layers that enhance resistance to gate current leakage and inhibit boron penetration into the dielectric [86]. However, as the sizes of semiconductor devices continue to decrease, it is becoming critical to develop processing methods that afford control of film properties and composition at the monolayer level. Lowtemperature remote plasma processing offers distinct advantages over other ultrathin silicon oxynitride film preparation methods [33]. For example, researchers [8788] have reported that the concentration profiles of nitrogen and oxygen atoms within SiOxNy films grown by remote plasma processing on Si(100) can vary significantly with both the feed gas composition and processing protocol that is used. Of particular interest is the observation that the exposure of Si(100) to an N20 or N2/02 plasma produces a film with a highly heterogeneous composition profile in which the nitrogen atoms accumulate at the filmSi interface and the oxygen atoms remain close to the filmvacuum interface. It was suggested that this composition profile results, at least in part, by chemical reactions in which gaseous oxygen atoms efficiently scavenge nitrogen atoms located closest to the filmvacuum interface. In a separate study, Watanabe and Tatsumi [89] also observed a decrease in surface nitrogen concentration after exposing a nitrided Si(100) surface to an oxygen plasma, and asserted that nitrogen atoms are removed from the surface through reactions with gaseous oxygen atoms. While thermally activated reactions between accommodated oxygen and nitrogen atoms could have taken place at the surface temperature of 750 C used in the work of Watanabe and Tatsumi, the other investigations cited were conducted with the substrate held at 300 C, which is well below the temperature for appreciable thermal decomposition of silicon oxide [30] and nitride surfaces [9091]. It is therefore likely that the oxygeninduced removal of nitrogen from the surfaces of these films occurs by nonthermal processes such as EleyRideal abstraction. In the present study, we find that direct nitrogen abstraction by gaseous oxygen atoms is energetically favorable when nitrogen is bound at the Si(100) surface in coordinatively unsaturated configurations, but less so when the nitrogen is triply coordinated with surrounding silicon atoms. 3.2 Computational Approach Quantum chemical calculations of oxygen induced abstraction of a nitrogen atom from the Si(100)(2x1) surface were performed using density functional theory (DFT) and cluster models of the surface. Cluster energies were computed using the unrestricted, hybrid threeparameter Becke method (UB3LYP), which combines the gradientcorrected exchange functional of Becke [49, 53] with the LeeYangParr (LYP) correlation functional [50]. An unrestricted approach was used for all calculations because the reactants, products and transition states are open shell structures. Geometry optimization methods based on the Berny algorithm were used for calculating local minima [7374], while transition state structures were obtained using a multiple step procedure. In this approach, an initial transition structure is obtained using the quadratic synchronous transitguided quasiNewton method (QSTN) [74, 82]. A frequency analysis is then performed to determine the normal mode of vibration that has the largest imaginary frequency. In the final step, this normal mode is selected to search for a transition structure using an eigenvector following method. This extensive approach was found to be more robust than the standard search routines available in the Gaussian03 package for identifying transition structures for the open shell systems we investigated. After convergence, a final frequency analysis was performed to confirm that the local minima and transition structures were zero or first order saddle points, respectively [92]. Figure 31 shows the Si9H12 onedimer cluster that was used to model the Si(100)(2xl) surface in these calculations. This cluster model has been widely used to study chemical reactions on Si(100)(2xl) using DFT [10, 93101], mainly because it is the smallest structure that adequately represents the main structural characteristics of the Si(100)(2xl) surface such as the tilted silicon dimer bond and welloriented sp3 covalent bonds. DFT calculations using the onedimer cluster have been found to accurately predict bond energies and barriers for several reactions on Si(100)(2xl) [10, 99] for which nonlocal electronic effects [98] are of minor importance. Hydrogen atoms are used to terminate the bonds of the Si cluster (i.e., hydrogen atoms are used as substituents for the bulk silicon atoms that are removed by truncation of SiSi bonds at the exterior of the cluster). These H atoms preserve the tetrahedral sp3 bonding environment of subsurface Si atoms, mimic strain that the bulk silicon atoms would impose on the boundary of the cluster, and generally have a negligible effect on the quantum chemistry calculation itself [99, 102]. Geometric constraints are imposed on the hydrogen atoms to improve the simulation of bulk strain effects as follows: All H atoms in place of bulk Si atoms are held fixed in their ideal tetrahedral configurations, while those replacing silicon atoms in neighboring dimers are fixed in positions that mimic the buckled dimer structure [96, 99]. The third and fourth layer Si atoms are not directly constrained, although their displacements are hindered by the constraints imposed on their neighbors and terminating hydrogen atoms; the fourth layer Si atom is practically held fixed in space due to the limited displacement of the third layer silicon atoms. Finally, all chemically active atoms, which include the silicon dimer atoms, the second layer Si atoms, and the nitrogen and oxygen atoms, are allowed to fully relax during the geometry optimizations. To reduce the computational expense of the calculations, a mixed basis set was used to expand the electronic wave function. A diffuse triplep plus polarization 631 l1++G(2d,p) basis set was used to describe the chemically active atoms while the remaining subsurface silicon atoms and terminating hydrogen atoms are described with a doublep plus polarization 631G(d) basis set [92]. All the structures investigated have a spin doublet multiplicity and the calculations were done using the Gaussian03 program [71]. 3.3 Results 3.3.1 Bonding Configurations of a Nitrogen Atom on Si(100)(2xl) We investigated the pathways for nitrogen abstraction from four different configurations of a nitrogen atom bonded on the Si(100)(2x1) surface to explore the influence of the local surface bonding environment on these reactions. The NSi(100) structures that we examined were recently predicted by Widjaja et al. [98] to be possible equilibrium structures resulting from incorporating a single nitrogen atom into the silicon surface. To benchmark our calculations, we optimized the NSi(100) structures using the same computational procedure as used by those authors. The structures and the corresponding energies of formation that we predicted are shown in Figure 32, where the zero of energy is taken to be the energy of the nitrogen atom and silicon cluster at infinite separation. Each structure has a spin multiplicity of two. The lowest energy structure, which has an energy of formation of 2.69 eV, is obtained when the nitrogen atom bonds to a single dimer atom [R(ad)]. The remaining structures are the N atom bonded with both dimer atoms in an epoxidelike structure [R(db)], the N atom inserted into a SiSi backbond [R(bb)] and finally the N atom bonded with one dimer atom and two second layer Si atoms [R(sat)]. The geometrical properties of these structures are in excellent agreement with those reported by Widjaja et al. [98], and the energies of structures R(ad), R(db) and R(sat) differ by less than 0.11 eV from their results. However, the energy of formation of structure R(bb) is found to be higher by 0.38 eV from the prediction of Widjaja et al. The differences in the energies predicted in these studies most likely arise from slight differences in the way geometric constraints are imposed in the calculations. It is noted that the energy of structure R(bb) is particularly sensitive to the positions of the terminating hydrogen atoms. Despite this difference, the geometries and trends in the relative energies of the NSi(100) structures are very close to those reported in the study of Widjaja et al. [98]. Finally, each of the structures (Figure 32) is predicted to be a local minimum on the doublet NSi(100) potential energy surface, and the lowest barrier for interconversion between the structures is 0.48 eV [98]. These characteristics suggest that at moderate surface temperature each structure could exist in appreciable concentrations during the initial stages of nitridation. 3.3.2 Nitrogen Abstraction by a GasPhase Oxygen Atom We investigated reactions between species only in their respective electronic groundstates. All of the reactions occur on a doublet potential energy surface and Firstlayer dimerr) Secondlayer Thirdlayer [100] Fourthlayer [110] [110] Figure 31. The Si9H12 cluster used in our UB3LYP calculations. Si atoms are represented by dark spheres and H atoms are shown as lightcolored spheres. may be represented by the general equation, 0(3P) + NSi(100) (Ms = 2)  NO(21) + Si(100) (Ms= 1), where Ms is the initial spin multiplicity of the surface cluster. For each of the NSi(100) structures investigated, nitrogen abstraction by a gasphase Oatom was found to be highly exothermic. This conclusion is reached by considering that the bond energy of NO in its doublet groundstate is about 6.5 eV, whereas the energies of formation of the NSi(100) structures range from about 2.7 to 4.7 eV (Figure 32). We initially explored direct pathways for abstraction in which the NO bond forms and the SiN bonds break in a single elementary step, as depicted by the reaction equation shown above. Extensive searches revealed a transition structure for only one singlestep abstraction process, namely, direct abstraction from the R(sat) structure. For structures R(ad), R(db) and R(bb), the abstraction pathways exhibit energy wells between the reactants and products due to the formation and interconversion of adsorbed NO species. These molecular precursor structures are produced when the 0(3P) atom attaches directly to a nitrogen dangling bond, yielding an NO species that is bound to the surface. A C B D Figure 32. Structural information of NSi9H12 clusters. A) R(ad). B) R(db). C) R(bb). D) R(sat). All distances are in A. Terminating hydrogen atoms have been excluded for clarity. Si atoms are shown as lightcolored spheres and N atoms as darker spheres. 3.3.3 Abstraction of N Adsorbed at the Dangling Bond [R(ad)] Figure 33 shows the predicted pathway by which an oxygen atom abstracts a nitrogen atom that is adsorbed on a single Si atom of the surface dimer [R(ad)]. The reaction involves the initial formation of molecular precursor MP(ad) (Figure 34B), which then decomposes to produce a gaseous NO molecule and the bare Si cluster, structure P (Figure 34C). The production of MP(ad) is highly exothermic (4.84 eV), since this reaction involves the formation of a strong NO double bond without cleavage of other bonds in the cluster. The subsequent decomposition of MP(ad) via the formation of a gaseous NO molecule is endothermic by 0.91 eV. A transition structure could not be found along the path from structure MP(ad) to P so only the thermochemical barrier must be overcome in this step. Although the molecular precursor lies at a lower energy than the abstraction product, a minimum of 3.93 eV of energy would need to be dissipated to the solid for the molecular precursor to stabilize. Moreover, since the NO species couples to the solid mainly through the SiN bond, it is likely that the SiN bond 0 ReAfrence; S0 + + : 7.53 MP(ad) Figure 33. Reaction pathway for the nitrogen abstraction from R(ad). Energiesarein units of eV and the zero of energy is taken to be the groundstate nitrogen atom, the 0(3P) atom and the singlet Si cluster at infinite separation. atom, the 0(3P) atom and the singlet Si cluster at infinite separation. will become highly excited and break during NO bond formation since only 0.91 eV are needed to detach the NO molecule from the structure. This result therefore suggests that nitrogen abstraction from structure R(ad) will be highly facile, occurring effectively in a single atomsurface collision. 0 110.81 Figure 34. Critical point structures of the nitrogen abstraction from R(ad). A) R(ad) B) MP(ad). C) P (final product). All distances are in A and the terminating hydrogen atoms have been excluded for clarity. The geometric and electronic changes that occur during the reaction shown in Figure 33 provide additional insights for understanding this reaction. In reactant R(ad), two unpaired alpha electrons are localized on the N atom, and one beta electron resides mainly on the triply coordinated Si atom of the surface dimer. Since the 0(3P) atom also has two unpaired electrons, its interaction with the adsorbed N atom results in the highly exothermic formation of an NO double bond. The NO bond in structure MP(ad) is close in length to that of the NO bond of an isolated NO molecule (1.20 vs. 1.15 A, respectively, Figure 34). The predominant change in the relative positions of the Si and N atoms that occurs during the first reaction step is elongation of the SiN bond from 1.75 to 1.91 A, which indicates a weakening of the SiN bond when the NO double bond of MP(ad) is formed. In the final reaction step, the molecular precursor decomposes by cleavage of the SiN bond, with one electron transferring to a ;r* orbital of the NO molecule and the other remaining at the surface. This remaining electron experiences a weak ;rinteraction with the initially unpaired electron on the opposing dimer atom since the bare Si surface in these calculations is taken to be the singlet, buckled dimer structure that is thought by many to be the groundstate of the Si(100)(2xl) surface [103104]. 3.3.4 Abstraction of the Nitrogen Bonded Across the Dimer [R(db)] Shown in Figure 35 is the pathway predicted for the abstraction of a nitrogen atom bonded across the surface dimer [R(db), Figure 36(A)]. The first step in this pathway is NO bond formation, resulting in structure MP(db). The exothermicity of this reaction is 3 eV, which is quite significant but is still lower by 1.84 eV than NO bond formation in the analogous reaction R(ab)  MP(ab) (Figure 33). Formation of structure MP(db) is followed by NO migration from the bridging site to the dangling bond site from which NO desorbs. The migration path is predicted to have an energy barrier of 0.48 eV, and the final desorption step presents a thermochemical barrier of 0.91 eV as discussed. Each of these barriers is significantly less than the energy released in the initial formation of the NO bond. Thus, unless energy is efficiently dissipated away from the initial collision zone, NO desorption can occur rapidly by this pathway as well. 52 0 0 + + .joev O Oxygen .'. Nitrogen 4.66 O R(db) TS(dbad) 7.52 7.66 MP(ad) MP(db) Figure 35. Reaction pathway of nitrogen abstraction from R(db). Energies expressed are in eV. The geometric changes that the cluster undergoes during these reactions steps are illustrated in Figure 36. NO bond formation to produce MP(db) causes each SiN bond to stretch from 1.68 to 1.80 A, and the SiSi dimer bond to contract from 2.55 to 2.37 A, which is indicative of weaker SiN bonds and a stronger Si dimer bond. Also, the NO bond in MP(db) is longer than that in MP(ad) (1.27 vs. 1.20 A), suggesting a weaker NO bond in the former. These changes may be understood by considering the corresponding changes in the alpha spin densities on the surface atoms. In the initial reactant R(db), an unpaired electron is distributed mainly between the Si atoms of the epoxidelike ring. The alpha spin density on these Si atoms decreases to zero upon formation of molecular precursor MP(db), and an unpaired electron is distributed almost evenly between the N 126.0e A .C .,, ,.4 B , D 1'. Figure 36. Critical point structures of the nitrogen abstraction from R(db). Structural information and relative energy. A) R(db). B) MP(db). C) TS(dbad). D) MP(ad). All distances are in A. The terminating hydrogen atoms were excluded for clarity. The final product is shown in Figure 34(C). and 0 atoms in the MP(db) structure. These changes indicate that NO bond formation leading to MP(db) involves a pairing between the lone electron on the reactant R(db) and an unpaired electron on the incident 0(3P) atom. The second lone electron on the 0(3P) atom is transferred to a 7r*like orbital between the N and 0 atoms. The resulting NO bond effectively has a bond order of 1.5 and is therefore weaker and longer than the NO double bond of MP(ad). The migration of the NO species from the bridging to the dangling bond position is activated by 0.48 eV, since it involves the simultaneous cleavage of a SiN bond and the formation of an NO double bond. In this reaction, an electron from the initial SiN bond is transferred to the NO species to form the double bond, while the second electron remains localized on the Si atom. The transition structure TS(dbad) for this reaction is shown in Figure 36(C). To reach TS(dbad) from MP(db), the NO species tilts away from one of the Si dimer atoms, causing the SiN bond to stretch to 2.44 A. At TS(dbad), the NO bond contracted to 1.21 A, and the remaining SiN bond is elongated from 1.68 to 1.85 A. On the path from TS(dbad) to MP(ad), the structure adopts a more stable configuration as the NO species tilts farther away from the opposing Si dimer atom, and the NO bond contracts to it final value of 1.20 A. 3.3.5 Abstraction of the Nitrogen Bonded at a Backbond [R(bb)] We found two pathways for nitrogen abstraction from the backbonded position [structure R(bb)] (Figure 37). For each path, the oxygen atom is predicted to first bond directly with the nitrogen atom of structure R(bb) to generate molecular precursor MP(bbs) (Figure 38B) in which the NO bond is nearly parallel to the bisector of the SiNSi bond angle. In one pathway, the 0 atom tilts toward the Si atom to form an SiO bond resulting in structure MP(bbt) which can then decompose by sequential NSi bond cleavage, forming MP(ad) and then a gaseous NO molecule. The alternative pathway involves direct NSi bond cleavage, causing the symmetric structure MP(bbs) to convert directly to MP(ad) from which the NO molecule detaches (Figure 37, inset). The reaction to produce MP(bbs) from the nitrogen backbonded structure is exothermic by 2.68 eV, and only a small energy barrier of 0.03 eV must be overcome for MP(bbs) to transform into the more stable molecular precursor MP(bbt). An energy barrier of 0.38 eV must be overcome for MP(bbt) to transform to MP(ad), while a barrier of 0.55 eV is predicted for the direct conversion of MP(bbs) to the MP(ad) structure. The overall reaction is exothermic by 2.25 eV, and the energy released during the initial NO bond formation is much greater than the barriers that must be surmounted for the NO molecule to desorb from the cluster. The general bonding characteristics of the MP(bbs) structure (Figure 38B) are similar to those of the MP(db) structure (Figure 36B) discussed above. In particular, formation of the NO bond in MP(bbs) involves a pairing between lone electrons on the N and 0 atoms, with the second lone electron of the 0 atom transferring to a z*like orbital localized along the NO bond. Despite these similarities, however, structure MP(bbs) is less energetically favorable than MP(db) by 0.61 eV. A comparison of the structures (Figures 36B and 38B) suggests that MP(bbs) is more strained since its formation is accompanied by an increase of the SiNSi bond angle and elongation of the SiN bonds. The NO bond in MP(bbs) is also longer than in MP(db), which suggests that additional strain would be imposed if the NO bond were to achieve its optimum length. One of the pathways for N abstraction involves the conversion of structure MP(bbs) to MP(bbt) (Figure 38D), which is exothermic by 0.25 eV and presents a barrier of only 0.03 eV (Figure 37). Although the energy barrier is quite small for this reaction, vibrational analysis indicates that structures MP(bbs) and TS(bbsbbt) do correspond to zero and first order saddle points, respectively. The exothermicity for this reaction step is also relatively small mainly because unfavorable structural changes offset the energy gained in forming the strong SiO bond. The predominant structural changes include stretching of the NO bond from 1.32 to 1.51 A and an increase of the SiNSi bond angle from 132.60 to 164.50. Elongation of the NO bond indicates a weakening of this bond, and the substantial increase in the Alternative path 31....r o 0 TS(bbsad)'t at D e0 7.05 " dQ. .0 MP(bbs) 7.53 4.37 ______,______ R(bb) C O 0Q T r6 62 7.05 7.02 r6 92> P MP() TS(bsbbt TS(b d P(bbt) 53 MP(a d) 0 + + :1 A.0 eoov 0 Oxygen Nitrogen Figure 37. Reaction pathways for the abstraction of the nitrogen atom inserted in a SiSi backbond. Main panel: 0 + R(bb) MP(bbs) TS(bbsbbt) MP(bbt)  TS(bbtad)  MP(ad)  P. Top inset: 0 + R(bb)  MP(bbs) TS(bbsad)  MP(ad)  P. Energies are in units of eV. SiNSi bond angle imparts strain on the SiSi bonds in underlying layers. Energy is also required for this reaction because an electron must be removed from a ;rlike bond along the Si surface dimer for the SiO bond to form. The remaining steps in this abstraction pathway include migration of the NO species of structure MP(bbt) to the dangling bond position, and then NO desorption from MP(ad). The migration reaction MP(bbt) > MP(ad) is exothermic by 0.23 eV and has an energy barrier of 0.38 eV (Figure 37). The geometry of the transition structure TS(bbtad) (Figure 39B) reveals that the NO species tilts away from the Si dimer atom, causing the SiO bond to break and the NO bond to strengthen and contract early along the reaction path. Concurrently, the SiNSi bond angle decreases to allow the SiSi backbond to begin forming and the lower SiN bond to break. Formation of the SiSi backbond and the NO bond is completed on the MP(ad) side of the barrier. In the second abstraction pathway for N(bb), the symmetric MP(bbs) structure converts directly to MP(ad) by the path shown in the inset of Figure 37. This reaction step is exothermic by 0.48 eV and has an energy barrier of 0.55 eV. The corresponding transition structure TS(bbsad) for this path is shown in Figure 310(B). The bond 0 1058, A 103.20 B D Figure 38. Molecular precursors formed after Ochemisorption onto R(bb). Structural information of the structures involved in the reaction pathway. A) R(bb). B) MP(bbs). C) TS(bbsbbt). D) MP(bbt). All distances are in A. The terminating hydrogen atoms were excluded for clarity. lengths in TS(bbsad) are nearly identical to those in TS(bbtad) (Figure 39B), which is not surprising since the reactions that proceed through these transition states have a common product. Nevertheless, the energies of these transition structures, referenced to the same initial state, differ by 0.42 eV. The most striking structural difference is that in TS(bbtad) the NO bond is nearly parallel with the plane defined by the SiNSi ring, whereas the NO bond is significantly tilted out of this plane in TS(bbsad). This suggests that overlap of the ,rlike orbitals in the NO species with those in the SiNSi ring is enhanced in the planar configuration of TS(bbtad), thereby lowering the energy of the structure relative to that of TS(bbsad). Figure 39. Structures formed during nitrogen abstraction from MP(bbt). A) MP(bbt). B) TS(bbtad). C) MP(ad). D) P. All distances are in A. A B) TS(bbsad). C) MP(ad). D) P. All distances are in A 3.3.6 Abstraction from the N=Si3 Structure [R(sat)] Figure 311 shows two possible pathways for N abstraction from the N Si3 structure [R(sat), Figure 312(A)] by an incoming 0(3P) atom, and one pathway for the adsorption of the 0 atom at a Si dangling bond site. The first abstraction pathway is a direct process wherein the incident 0 atom removes the N atom from the surface in a single step. This reaction is the only singlestep abstraction process that was identified in the study. Direct N abstraction from the N Si3 structure is predicted to be exothermic by 2.57 eV and presents a relatively small energy barrier of 0.20 eV. The N0 separation at the transition structure is 1.23 A (TS(satP) (Figure 312B), which is only 0.12 A longer than the NO bond of the isolated NO molecule. In addition, each secondlayer SiN bond has broken and two SiSi backbonds have formed on the path to the transition structure. This indicates that direct abstraction from NSi3 is a concerted process in which three SiN bonds break while two SiSi backbonds and a NO bond simultaneously form during a single gassurface collision. This reaction would require the upper Si dimer atom to move toward the second layer Si atoms by a considerable distance, and the N atom to simultaneously move out of the plane defined by these three Si atoms as the 0 atom approaches the cluster and the NO bond begins to form. The other possible pathway for N abstraction from the NSi3 structure is also a concerted process in which the interaction of the Oatom with the cluster produces molecular precursor MP(bbs) (Figure 38B), which then decomposes by the pathways shown in Figure 37. However, neither an intermediate structure nor a transition structure could be found for the reaction 0(3P) + R(sat) " MP(bbs). Thus, for this reaction to occur, an NO bond and a SiSi backbond must form, and a SiN bond must break during a single collision of the Oatom with the cluster. Although neither of these abstraction pathways is energetically prohibitive, they each involve the formation and scission of multiple bonds in a single gassurface collision. Moreover, these reactions require that the Oatom come in close proximity to the dangling bond of the Si dimer atom where it could adsorb. Indeed, adsorption of the Oatom at the dangling bond site (Figure 311) is predicted to be barrierless and highly exothermic (6.35 eV), and is therefore concluded that to be the predominant outcome of the interaction between an 0(3P) atom and the NSi3 structure. As may be seen in Figure 312, adsorption of an 0(3P) atom at the Si dangling bond causes the local NSi3 configuration to change negligibly, but it does result in contraction of the SiSi dimer 61 . 3.85 " 4 05 >TS(satP) 4.05', R(sa t) I N S /'TS(bbsad)', 6.62 7.05 \ c p MIP (bbs) 7.53 1MP(ad) Q + + J>. 6J, UOeV 0 *',Q Oxygen C Nitrogen 10.41 O(ad)R(sat) Figure 311. Reaction pathways for Nabstraction from (and Oatom adsorption on) R(sat). Energies are in units of eV and the zero of energy is taken to be the groundstate nitrogen atom, the 0(3P) atom and the singlet Si cluster at infinite separation. bond. This contraction indicates that OSi formation alleviates repulsive interactions between the two unpaired electrons that are present on the silicon dimer atoms of the N=Si3 structure. In fact, we find that the energy of Oatom adsorption at the Si dangling bond site of R(sat) is 1.9 eV higher than that for 0adsorption on the clean Si(100)(2xl) surface. This indicates that significant electronic repulsion exists along the dimer of the R(sat) structure, and that Oatom adsorption is highly favorable on this surface. IC.S* Figure 312. Structures formed during the direct nitrogen abstraction from R(sat). A) 0(3P) attacks the nitrogen atom of R(sat). B) Transition structure [TS(satP)]. C) NO molecule desorbed into the gas phase and bare Si surface (P). All distances are in A. The terminating hydrogen atoms were excluded for clarity. O 04 8, Figure 313. Structures involved in the Oatom chemisorption on R(sat). Structural information and relative energy. A) 0(3P) atom approaching structure R(sat) and B) Oatom adsorbed at the Si dangling bond site [O(ad) R(sat)] All distances are in A. The terminating hydrogen atoms were excluded for clarity. 3.4 Discussion Nitrogen abstraction by a gasphase 0(3P) atom is highly exothermic for each of the NSi(100) structures investigated in this study. However, abstraction is predicted to occur in a single step only for the reaction of 0(3P) with the coordinatively saturated N atom of the N=Si3 structure, and (in this case) singlestep abstraction appears to be much less probable than 0 adsorption due to the energetic differences in these reactions and because multiple bonds must break and form for singlestep abstraction to occur. For each of the coordinatively unsaturated NSi(100) structures, N abstraction is predicted to occur by a precursormediated pathway that is initiated by the formation of an NO bond and the release of between 2.7 and 4.8 eV into the surface. Since the subsequent elementary steps leading to NO desorption have barriers that are less than 1.0 eV, NO bond formation provides the system with excess energy that could readily promote local bond rearrangements and ultimately NO desorption. Indeed, nitrogen abstraction by such a pathway would effectively be an EleyRideal process since the NO product would evolve into the gasphase within no more than a few vibrational periods after the initial gassurface collision. Alternatively, the energy released during NO bond formation could be dissipated away from the surface bonds, allowing the adsorbed NO species to equilibrate in one of the configurations predicted. However, considering the significant amount of energy (2 to 4 eV) that would need to be channeled away from the initial collision zone for a precursor to stabilize, it is reasonable to expect NO desorption to be the more likely outcome of initial NO bond formation in these systems. Nevertheless, calculations of the dynamics for abstraction, and particularly the efficiency of energy dissipation to the solid, are needed to explore the propensity for 0(3P) atoms to abstract nitrogen from Si(100) by the pathways predicted here. Prior experimental studies have provided indirect evidence for nitrogen abstraction from Si surfaces by gaseous atomic oxygen. For example, exposure of Si(100) to plasmas containing both 0 and N atoms results in a surface that is depleted of nitrogen [87, 88, 105]. Similarly, treating nitrided Si(100) with an oxygen plasma has been reported to lower the surface nitrogen concentration [89, 106]. We recently conducted reactive scattering experiments in UHV to directly examine the abstraction of nitrogen from Si(100) by an atomic oxygen beam, but we did not observe the evolution of gaseous NO using mass spectrometry or nitrogen depletion from the surface [107]. However, the absence of measurable nitrogen abstraction in our experiments can be attributed to the bonding state of nitrogen that was investigated. In that work, nitrogen was incorporated into Si(100) by thermally decomposing NH3 on the surface at a substrate temperature of 900 K; adsorption at this high surface temperature is necessary to ensure the complete desorption of hydrogen. Above about 700 K, however, nitrogen has been found to diffuse into the subsurface region of Si(100) [108] and is apparently inaccessible for direct abstraction by gaseous Oatoms. Our experimental observations therefore suggest that the NSi(100) structures considered in the present computational study are not stable at temperatures greater than about 700 K. In fact, recent quantum chemical investigations predict that a single nitrogen atom does have an energetic preference to bond in the subsurface of Si(100) [109110]. Nevertheless, the NSi(100) structures investigated in the present study, which were reported originally by Widjaja et al. [98], may be stable at surface temperatures lower than 700 K since, to our knowledge, barriers for the migration of nitrogen into the Si(100) subsurface have not been reported, and may be high enough to enable nitrogen atoms to stabilize at the surface at low to moderate temperatures. Hence, nitrogen abstraction by gaseous Oatoms may indeed occur during plasmaenhanced oxynitridation of Si(100) if the surface temperature is maintained sufficiently low. Another possibility is that the barriers for nitrogen diffusion to the subsurface are small, which would imply that the NSi(100) structures that we investigated in this computational study (Figure 32) do not exist in appreciable concentrations at low to moderate surface temperature. In this case, a mechanism other than direct abstraction by gaseous oxygen atoms would be needed to explain observations of nitrogen depletion at the Si(100) surface by reaction with oxygen plasmas. Experiments to directly investigate nitrogen abstraction by atomic oxygen will require a method for adsorbing N atoms onto Si(100) at low surface temperature, which should be possible using an active nitrogen source such as gaseous Natoms or through nonthermal activation of nitrogencontaining adsorbates. CHAPTER 4 CHEMISTRY ON Si(100)(2x1) DURING EARLY STAGES OF OXIDATION WITH O(3P) 4.1 Introduction Oxygen adsorption on Si(100)(2xl) is an important process because it corresponds to the initial stages of the formation of the Si/Si02 interface, a system of high technological relevance in microelectronic devices, especially in the metaloxide semiconductor fieldeffect transistor (MOSFET) gate dielectrics. Because of the demands imposed by Moore's law, the thickness of the silicon dioxide gate dielectrics grown on Si(100) is crucial for the future of the microelectronic industry. As these layers approach the sub nanometer region (it is now possible to fabricate devices with Si02 gate thickness of 1.3 nm [34], which correspond to about ten silicon atoms across the Si/Si02 interface) it is important to understand the detailed mechanism of initial oxidation of silicon surfaces. This small thickness is a limit for the use of Si02 gate dielectrics given that the concentration of undesirable silicon suboxides (Si Si2+ and Si3+) in the Si/Si02 interface becomes increasingly greater. Because of this, aggressive research for alternative gate dielectric materials is taking place, but even for these new gate highK dielectric materials, knowledge of the initial steps of Si02 formation is critical given that many of them still involve one or two monolayers of Si02 deposited over the gate channel region. Several methods to fabricate silicon dioxide have been used successfully over the year [3334]. However, because of the continuous scaling of the devices, these traditional methods are becoming increasingly less efficient. Furthermore, as the size of the silicon dioxide layers approach molecular and atomic dimensions, it would be ideal to have a fabrication method based on molecular chemistry and low temperature, so the product obtained by the manufacturing process remains in a given configuration and is not affected by thermal changes of the system. One of these methods is plasma processing of silicon substrates. The interactions of radical species such as oxygen atoms, which can be present in large fractions in plasma environment, with silicon surfaces are relatively poorly understood. Engel and coworkers [2631] reported detailed kinetics studies of the 0 atom oxidation of clean silicon using a plasmabased atomic beam source. Also, Yasuda et al. [111] studied the 0 atom oxidation of hydrogenterminated silicon using a hotfilament source. These studies revealed that 0 atom oxidation is favorable compared with that due to 02, and that several monolayers of oxide can be formed efficiently via the direct insertion of 0 atoms into nearsurface bonds. We explore the structures formed during the initial steps of 0(3P) incorporation on clean Si(100)(2xl) using small Si9H12 cluster models and gradientcorrected density functional theory (DFT). From the time when early studies of asymmetric dimer descriptions [112] and initial adsorption of hydrogen on Si(100) [94] were made, the Si9H12 cluster has been proven a very useful tool for performing theoretical studies of local chemistry on Si(100)(2xl). Recent comprehensive oxidation studies which combine Si9H12 clusters with gradientcorrected density functional theory (DFT) techniques yielded accurate energetic and vibrational frequencies for several oxidation products, although none of them has treated oxidation by 0(3P) specifically. Most of previous work has focused on the waterinduced oxidation of the clean Si(100)(2x 1) surface [34, 95, 96, 104, 113, 114], and some on oxidation by molecular oxygen [100]. What has been found so far using the DFT/cluster approach can be summarized as follows: Initially, water adsorbs dissociatively on top of single silicon dimer and forms H and OH fragments which are most stable when they are bonded with the hydroxyl oriented away from the surface dimer bond [95, 104]. Thermodynamic [96] and mechanistic studies [113114] revealed a strong tendency for 0 to agglomerate on the dimers of Si(100) and predicted that the three and fiveoxygen agglomerated structures were the most stable. Formation of epoxidelike rings upon dehydrogenation of the surface in the event of oxygen agglomeration was also reported. This tendency for 0 to agglomerate was found to be thermodynamically driven since it is energetically more favorable to have one dimer with n oxygens, for n=25, and (n1) oxygenfree dimers, than it is to have n dimers each with one oxygen. In this work, we performed DFT/cluster calculations to determine minimum energy geometries to elucidate the relative thermodynamic energies of the possible oxidation products that can form by insertion of up to three 0(3P) atoms, and we use this result to predict a minimum energy reaction path. These oxidized Si9H12 clusters were classified in three different sets of isomers (OiSi9H12, 02Si9H12 or 03Si9H12) depending on the number of oxygen atoms adsorbed on the surface. The energy of these isomers was then compared to one another to determine how it was affected by the formation of the two SiO bonds in place of one SiSi bond, by the oxidation state of the surface silicon atoms, by the spinstate of the surface (spinsinglet vs. spintriplet) and by strain effects (all quantified by using a bond energy model that assumes that they are independent and additive). Finally, we investigated transition states for the insertion process, to develop a model for the preferred mechanism of the initial steps of Si(100) oxidation by 0(3P). 4.2 Theoretical Approach Our theoretical approach was based on KohnSham density functional theory (DFT) calculations of clusters of silicon that represent the local bonding arrangement of the surface. The silicon cluster used was the Si9H12 (Figure 31) which is the smallest structure that appropriately represents the main structural characteristics of the Si(100)(2xl) surface (i.e., covalent tetrahedral sp3 arrangement of SiSi bonds and a tilted dimer). This cluster has twelve H atoms to terminate all dangling bonds resulting from truncation of SiSi bonds at the exterior of the cluster; these terminating hydrogen atoms preserve the tetrahedral bonding and have a negligible effect on the predicted energies of the different clusters [102, 115]. The clusters were constrained by imposing boundary conditions that mimic the strain that the bulk silicon would impose on the surface dimer under study. In particular, the hydrogen atoms were fixed in their positions along the directions of truncated SiSi bonds that they terminate. Two types of hydrogen atoms can be found in this cluster. The bulk hydrogen atoms, fixed along tetrahedral directions, and the neighboring dimer hydrogen atoms, constrained in positions that mimic nearest silicon dimers. The third and fourth layer silicon atoms were not directly constrained, but the displacement of these atoms is minimal because they are surrounded by fixed H atoms. All the chemically active atoms including the first and second layer silicon atoms were allowed to relax completely unconstrained. KohnSham density functional theory [4144] is used for the electronic structures calculations. Specifically, we use the B3LYP hybridgradientcorrected method [53, 116] which calculates the exchange correlation term by means of a linear combination of local, gradientcorrected and exact HartreeFock exchange terms with the Becke gradientcorrected term (B88) [49], and the local and gradientcorrected correlations terms of VoskoWilkNusair (VWN) [48] and LeeYangParr (LYP) [50], respectively. The electronic wave function was expanded using mixed Gaussian basis sets. Diffuse triple(plus polarization 631 l1++G(2d,p) was used for the oxygen atoms and first and second layer silicon atoms, while subsurface silicon and terminating H atoms were expanded with a double(plus polarization 63 1G(d) basis set. This approach, used successfully in similar studies [97100], focuses basis functions on the chemically active portion of the cluster and accurately describes orbitals involved in the reaction while minimizing the computational expense. All calculations were run using an unrestricted approach to calculate the openshell structures of the spintriplet systems. Comparison between restricted and unrestricted calculation results for closed shell systems (i.e., spinsinglet state clusters) showed negligible differences. The commercial program Gaussian03 [71] was used to run all the calculations. 4.3 Results and Discussion A primary goal of this study was to assess the factors which influence the thermodynamic stability of local structures formed in the early stages of Si(100) oxidation. Toward this end, we performed energy minimization for each possible isomer in which one, two or three oxygen atoms are inserted into the first and second layer SiSi bonds of the Si9H12 cluster. Since these calculations generated many structures, a shorthand notation of capital letters followed by a subscript and a code in parenthesis is used for labeling the clusters to facilitate discussion. The first capital letters which are DB, ME or TS stand for dangling bond, minimum energy, and transition state structure, respectively. Subscripts 1 or 3 appear after the capital letter to indicate whether the structure is a spin singlet or triplet. Finally, the code in parenthesis is based on the notation used by Stefanov and Raghavachari [96] in which SiSi corresponds to a non oxidized cluster, SiOSi is a cluster with the 0 atom inserted into the dimer bond, SiSiO2 is a structure with two oxygen atoms inserted into the two backbonds of the same dimer silicon atom, OSiSiO has two oxygen atoms inserted in a backbond pair at the same side of the dimer SiSi bond, OxSiSiO has the two oxygen atoms in oppositeside backbonds, and SiSiOd is a cluster with an oxygen atom bonded to a single silicon atom of the dimer. Several factors determine the heat of formation of an oxidized cluster, including the Si oxidation states, bond strain and the spin state of the cluster. We investigated clusters in both singlet and triplet spin states since the lowest energy structures of the clean Si(100) correspond to these spin states. The spin has two primary effects on the energy of the cluster. Firstly, a singlet cluster has one extra SiSi bond compared with a triplet cluster, which tends to lower the energy of the structure. However, the formation of this extra bond also alters the geometric structure of a cluster and hence the strain energy. The strain is typically higher in a singlet cluster than in the analogous triplet structure for the onedimer model. To directly compare the energies of singlet and triplet clusters, the energy of an oxidized cluster is defined with respect to a common reference state, which was taken to be the energy of the appropriate number of isolated oxygen atoms plus that of the clean (unoxidized) cluster in its singlet groundstate, [structure MEi(SiSi)]. The singlet bare cluster was chosen as the reference state because DFT predicts that this structure minimizes the energy of the clean Si(100) surface. The dangling bonds on the Si dimer atoms are unpaired in the spintriplet state of the bare cluster, [ME3(SiSi), Figure 41(a)], but they pair to form a weak ar bond along the dimer in the spinsinglet state, [MEi(SiSi), Figure 41(B)]. This extra bond in MEi(SiSi) shortens and tilts the silicon dimer bond along the [110] direction, resulting in an asymmetric structure that is 0.37 eV more stable than its symmetric ME3(SiSi) counterpart. The highest occupied molecular orbitals (HOMO) of both MEi(SiSi) and ME3(SiSi) were also calculated (Figure 41). The 0.37 eV energy difference, labeled as AEspin, reflects the energy gained by the 7c interaction, which has been reported to be 1.24 eV [117]. It also reflects the difference in strain energies between the singlet and triplet bare clusters. Since the initial surface is taken to be the singlet bare cluster, the production of an oxidized cluster in the triplet spin state may be considered to occur in two steps, namely, conversion of the bare cluster from the singlet to the triplet groundstate, and then insertion of oxygen atoms into the triplet bare cluster. Thus, when comparing the energies of the oxidized clusters, it is important to recognize that the singlet to triplet conversion step contributes an energy penalty of 0.37 eV to the heat of formation of a triplet cluster. This energy is taken into account in our analysis of other factors that determine the cluster energies as discussed below. We quantified factors which contribute to the energy of an oxidized cluster by invoking a bond energy model, which has previously been shown to provide a reasonable representation of the heats of silicon suboxide formation [118]. In the bond energy model, bond strain and distinct SiO bond strengths for Si+1, Si+2 and Si+3 oxidation states are treated as separable effects that contribute additively to the total heat of formation. In cases for which each oxygen atom in the cluster is inserted into a SiSi bond, the bond energy model yields the following expression for the energy of an oxidized cluster, AEUB3LP = N(2so s) + AEsuboxide + 5(AEsp + AEtran (41) where AEUB3LYP is the energy of the oxidized cluster calculated by DFT, Nis the total number of oxygen atoms in the cluster, s,o is the SiO bond energy in stoichiometric ft quartz crystal and s,, is the SiSi bond energy in bulk silicon. In addition, AEuboxde is the total suboxide penalty energy, AEpn is the 0.37 eV energy difference between the triplet and singlet bare clusters, where 8 is zero for the spinsinglet and one for the spintriplet clusters, and AEexess is the excess energy that is mainly related to changes in strain that result from oxygen insertion into the bare cluster. The suboxide energy penalty is defined by the Equation 42. AEsuboxide ZNs (Aj) (42) x where Nsgi is the number of silicon atoms bonded to x oxygen atoms (with x = 1, 2 or 3), and Ax is the energy penalty of a Si atom in the +x oxidation state. The concept of suboxide penalty energies was introduced in the bond energy model of Hamman [119] to take into account the effective increase in SiO bond strength as the Si oxidation state increases. This effect originates mainly from the greater amount of ionic character in SiO bonds involving Si atoms in higher oxidation states, but quantum resonance has also been suggested to enhance the bond strengths [119]. To quantify the differences in suboxide SiO bond strengths, penalty energies were computed by Hamann [119] and Bongiorno and Pasquarello [118] for each Si suboxide state, under the assumption that the penalty energies can simply be added to the (2esioesisi) term to compute the heats of formation of suboxide structures. Since the quantity esio is defined as the SiO bond energy in crystalline SiO2, the penalty energies have positive values so that their contribution reduces the energy change in forming SiO suboxide bonds. Table 41 shows the values of the suboxide penalty energies determined in the separate studies of Hamann [119] and Bongiorno and Pasquarello [118]. Relative stability [110] L [iT0] Structural information and highest occupied molecular orbital plot of clean Si9H12 clusters. A) Symmetric spintriplet surface ME3(SiSi). B) asymmetric spinsinglet surface MEi(SiSi). MEi(SiSi) is calculated to be more stable than ME3(SiSi) by 0.37 eV. Silicon and hydrogen are represented by grey and white balls, respectively. Bond lengths are expressed in A. Figure 41. Table 41. Suboxide penalty energies for various silicon oxidation states. The superscripts denote the number of oxygen atoms bonded to a particular silicon atom. Energies are given in eV. Oxidation state Energy penalty Ref 118 Ref 119 Si+1 A1 0.50 0.47 Si+2 A2 0.51 0.51 Si+3 A3 0.22 0.24 The penalty energies determined in these studies by Hamann [119] and Bongiorno and Pasquarello [118] are in good agreement with one another and are significant in value (Table 41), which suggests that maximizing the Si oxidation states is an important driving force that determines the types of local structures that form on Si(100) during initial oxidation. We used the bond energy model given in Eq. 42 to quantify changes in strain energy resulting from the insertion of oxygen atoms into SiSi bonds of the bare cluster. We chose to use an average of the suboxide penalty energies reported by Hamann [119] and Bongiorno and Pasquarello [118] in our analysis. This approach appears justified considering the close agreement in the penalty energies determined in those investigations. In addition, the values of Eso and Es,s were taken to be 4.35 and  2.02 eV, respectively, as reported for the SiO bond energy in crystalline Si02 (/f quartz, the lowestenergy form), and for the SiSi bond energy in bulk silicon [96, 117]. Additionally, for spinsinglet structures where the oxygen atom inserts into the silicon dimer, the appropriate value of s,s, in Eq. 42 corresponds to the energy for SiSi bondd cleavage (1.24 eV) [117]. These choices introduce some uncertainty in the absolute strain energies that are calculated because Eso and Es,s are experimentally determined values and are, therefore, not subject to the systematic errors that may affect the energies predicted by DFT. We also investigated clusters in which one of the oxygen atoms adsorbs on a single Si atom of the surface dimer and forms only one SiO bond, which we label as SiOd. These structures can serve as precursors to the formation of SiOSi linkages which form when the adsorbed 0 atom inserts into SiSi bonds near the surface. In fact, a recent experimental study by Gerrard et al. [107] shows that gaseous 0(3P) atoms initially adsorb at surface dangling bond sites on Si(100), forming SiOd bonds, before incorporating into the solid. Hence, investigating these dangling bond structures is essential for developing a mechanistic understanding of Si(100) oxidation by gaseous oxygen atoms. For the dangling bond structures, the following variation of Eq. 42 was used to quantify strain energies, AUB3LYP = (N X2sSiO SiSi) + SiOd + Esuboxide + 8(AEspin)+ Aexcess (43) where all quantities have the same definitions as given above except for eSiOd which is the energy of the SiOd bonds. We find that SiOd bonds have a more delocalized character than the SiO bonds in SiOSi linkages due to the lone electron from the oxygen atom. One consequence of this delocalization is that the bond energy model with the parameters defined above, yields excessive values for the relative strain energies for the dangling bond structures (Section 4.3.2). In virtue of this, the value of eSiOd has been defined by the difference in energy between the DB3(SiSiOd) and MEi(SiSi) structures. Finally, it is important to point out that charge transfer involved in the incorporation of adsorbed species on Si(100)(2xl) can exhibit nonlocal effects (i.e., charge transfer to a neighboring dimer on the surface) which may limit the overall usefulness of the OxSi9H12 clusters in estimating heats of formation. However, Widjaja and Musgrave [100] reported that the electron density in the similar system of 02 adsorbed on Si(100)(2xl) remains localized within the onedimer environment. This observation supports the use of the small Si9H12 cluster for our investigation. 4.3.1 Structures with One Adsorbed Oxygen Atom (OiSi9H12) Since we are interested in structures formed during the early stages of oxidation of the Si(100)(2xl) surface, we considered only oxygen atom insertion into the surface dimer bond, one or more backbonds or a dangling bond. Optimized total energies of the different oxidized clusters are used to evaluate their stability relative to that of the clean clusters. For now, we postpone discussing the mechanistic aspects (i.e., reaction barriers) of these reactions and focus exclusively on the relative thermodynamic stabilities of the different isomers. 3  3.5  5.5559 5.81  6.5  Structure Figure 42. Relative energies of OiSi9H12 isomers. A) DB3(SiSiOd). B) ME3(SiSiO). C) ME3(SiOSi). D) MEi(SiOSi). E) MEi(SiSiO). Energy reference: MEi(SiSi) + 0(3P) = 0.00 eV. Oxygen and silicon atoms are represented by black and gray balls, respectively. For the minimum energy analysis, the energies of the OiSi9H12 isomers are referenced to the reaction 0(3P) + Si9H12 (Ms = 1)  OiSi9H12 (Ms = 1 or 3) (44) where Ms is the spin multiplicity. Table 42 shows the different energies that contribute to the calculated UB3LYP energies of the O1Si9H12 clusters, in order of ascendant stability. The DB3(SiSiOd) structure (Figure 43A) is the least favorable cluster because it only has one new SiO bond. According to the bond energy model for dangling bond structures (Eq. 43) and the definition of the ESiOd energy, there is not any stress associated to DB3(SiSiOd). This value is not surprising, since the only difference in the geometric characteristics of the DB3(SiSiOd) and ME3(SiSi) structures (Figure 41A) is an 0.02 A elongation of each of the backbonds that are located on the oxygen side of the dimer. Table 42. Penalty energies of O1Si9H12 isomers. Contributions of the different factors to the calculated DFT/UB3LYP energies, according to the bond energy model. Energies are given in eV. For all structures: N= 1. Isomers are listed in order of increasing thermodynamic stability. Energy reference: MEi(SiSi) + 20(3P) = 0.00 eV. Structure N(2esioesisi)a ESiOd AEsuboxide 6(AEspin) MAexcess AEUB3LYP (Eq. 43) Dangling Bond Isomers DB3(SiSiOd) 0.00 4.74 0.00 0.37 0.00 4.37 Minimum Energy Isomers ME3(SiSiO) 6.68 N/A 1.00 0.37 0.28 5.59 ME3(SiOSi) 6.68 N/A 1.00 0.37 0.50 5.80 MEI(SiOSi)b 7.46 N/A 1.00 0.00 0.56 5.90 MEi(SiSiO) 6.68 N/A 1.00 0.00 0.49 6.17 a: corresponds to (N1) (2esioesisi) for dangling bond clusters. b: uses a value of 1.24 eV for the esisi, which corresponds to the SiSi xZ bond In the remaining O1Si9H12 structures, two SiO bonds are formed at the expense of one SiSi bond, so each structure has the same suboxide energy penalty. Thus, differences in the relative energies arise from the spin penalty energy and the different amounts of strain in the structures. Firstly, the spinsinglet isomers MEi(SiOSi) (Figure 43D) and MEi(SiSiO) (Figure 43E) are predicted to be the more energetically favorable structures, indicating that strain relief in the triplet structures is insufficient to Spintriplet 2.34. ,. =97.50 . =42.8 S=39.70 [100] [110] [11iTb] Spinsinglet Structural characteristics of Si9H12 clusters with one oxygen atom. A) DB3(SiSiOd). B) ME3(SiOSi). C) ME3(SiSiO). D) MEi(SiOSi). E) MEi(SiSiO). Oxygen and silicon atoms are represented by black and gray balls, respectively. Bond lengths are expressed in A. The crystalographic coordinate origin gives the orientation of the clusters. The surface is along the [100] plane. Figure 43. compensate the spin penalty energy. Structure MEi(SiOSi) is more strained, and therefore less favorable than MEi(SiSiO) by 0.27 eV, since formation of the epoxidelike ring stretches and weakens the SiSi dimer bond. This result agrees with observations made by Weldon et al. [113114], who were the first researchers to report the epoxidelike ring structures and suggested that they were the thermodynamically favored products after oxygen agglomeration (i.e., after three or more oxygen atoms absorb on the Si(100) surface). For the triplet structures, ME3(SiOSi) is more favorable than ME3(SiSiO) by 0.22 eV since oxygen insertion breaks the dimer bond and thereby relieves the strain that the dimer bond imposes on the bare cluster. Notice in Figure 43(B) that the silicon dimer atoms are separated by 3.00 A in the ME3(SiOSi) structure. Finally, the bond energy model predicts negative strain energy (i.e. strain relief, for almost all of the O1Si9H12 structures, suggesting that each structure possesses less strain than the corresponding bare cluster). The only exception is the MEi(SiOSi) structure, in which the predominant structural change is a significant elongation of the silicon dimer bond to a length of 2.57 A, that increases the strain in the cluster by 0.56 eV. Nevertheless, the silicon dimer bond is not cleaved in this structure given that the formation of the SiO bond only disrupts the weak SiSi interaction along the dimer. This formation does not break the stronger SiSi o bond, although it elongates the bond considerably. 4.3.2 Structures with Two Adsorbed Oxygen Atoms (02Si9H12) Figure 44 shows the relative thermodynamic stability of all the isomers with two oxygen atoms (02Si9H12), and the various contributions to these relative stabilities are presented in Table 43, in order of ascendant stability. The energies of the isomers are referenced to the reaction expressed in Eq. 45. 20(3P) + Si(100) (Ms = 1)  02Si(100) (Ms = 1 or 3) (45) where, as before, Ms is the spin multiplicity. As in the case of the OiSi9H12 clusters, the dangling bond structures are the least favorable structures, because they only have three SiO bonds instead of the four that form on all the other isomers given that one of the oxygen atoms bonds to a silicon dangling bond. An interesting structural effect is predicted by the DFT/UB3LYP calculations for the dangling bond structures which significantly influences the heats of SiOd bond formation. We find that SiOd bond SiSi bond lengths and angles relative to the reactant ME3(SiSi). In contrast, formation of an SiOd bond on the ME3(SiOSi) and ME3(SiSiO) structures causes one SiSi bond length to increase, thereby weakening the structure. In particular, with respect to the reactant structures, the dimer SiSi bond is stretched by about 0.07 A in DB3(OSiSiOd) and DB3(SiSiO2d), and one of the SiSi backbonds on the Od side of the dimer is Table 43. Penalty energies of 02Si9H12 isomers. Contributions of the different factors to the calculated DFT/UB3LYP energies, according to the bond energy model. Energies are given in eV. For all structures: N= 2. Isomers are listed in order of increasing thermodynamic stability. Energy reference: MEi(SiSi) + 20(3P) = 0.00 eV AEUB3LYP Structure N(2esioesisi)a ESiOd AEsuboxide 6( AEspin) AEexcess (Eq. 43) Dangling Bond Isomers DB3(OSiSiOd) 6.68 4.74 1.00 0.37 0.23 9.82 DB3(SiSiO2d) 6.68 4.74 1.00 0.37 0.27 10.32 DB3(SiOSiOd) 6.68 4.74 1.00 0.37 0.38 10.43 Minimum Energy Isomers MEi(OSiSiO) 13.37 N/A 2.00 0.00 0.33 11.04 ME3(OxSiSiO) 13.37 N/A 2.00 0.37 0.26 11.25 MEi(OxSiSiO) 13.37 N/A 2.00 0.00 0.06 11.31 ME3(SiSiO2) 13.37 N/A 1.51 0.37 0.05 11.44 ME3(OSiSiO) 13.37 N/A 2.00 0.37 0.70 11.69 ME3(SiOSiO) 13.37 N/A 1.51 0.37 0.57 12.05 MEi(SiOSiO)b 14.15 N/A 1.51 0.00 0.48 12.15 MEi(SiSiO2) 13.37 N/A 1.51 0.00 0.30 12.16 a: corresponds to (N1) (2esioesisi) for dangling bond clusters. b: uses a value of 1.24 eV for the esisi, which corresponds to the SiSi 7r bond. stretched by 0.10 A in the DB3(SiOSiOd) structure. As expected, this SiSi bond elongation lowers the heat of SiOd bond formation. For example, the heat of the reaction ME3(SiSi) + 0(3P) = VME3(SiSiOd) is 0.50 eV greater than the heat of the reaction ME3(SiSiO) + 0(3P) = ME3(OSiSiOd). This energy difference may be attributed almost entirely to the observed SiSi bond weakening since the SiOd bonds in both structures involve a Si+1 species, and the only appreciable structural difference between reactants and products is the longer SiSi dimer bond in ME3(OSiSiOd). It is tempting to conclude that SiSi bond weakening occurs in the O2Si9H12 structures because the Si atoms in these bonds have higher partial positive charges than in ME3(SiSiOd). For example, the dimer bond in ME3(SiSiOd) involves an Sio and Si+1 species, whereas in VIME3(OSiSiOd) the dimer Si atoms are both nominally in the +1 oxidation state. However, similar SiSi bond stretching is not observed in 02Si9H12 structures in which each oxygen atom is present in an SiOSi linkage (Figure 46). Hence, delocalization of the unpaired electron of the Od atom must play an important role in weakening SiSi bonds in the O2Si9H12 dangling bond structures. The remaining 02Si9H12 isomers result from the formation of SiO bonds in siloxane bridges or epoxidelike rings in place of the dimer bond or backbonds (Figure 46). The relative energies of these structures are determined by a combination of effects, though the suboxide energy penalty generally has the largest influence. The overall trend is that the structures possessing four Si+1 species are less favorable than those with a combination of two Si+ and one Si+2 species, and that the spinsinglet structures are more favorable than their spintriplet counterparts. Indeed, three of the four isomers possessing a Si+2 atom are the most energetically favorable among the two 0atom isomers and, in order of increasing stability, the most favorable structures are ME3(SiOSiO) < MEi(SiOSiO) < MEi(SiSiO2). The ME3(SiSiO2) structure falls out of this trend as it is 0.25 eV less favorable than the ME3(OSiSiO) structure which has four Si+1 species and, consequently, should be less stable. After accounting for the suboxide penalties, the bond energy model suggests a substantial 0.75 eV difference in the excess energies of these two structures, which probably arises from a combination of strain effects associated to the bonding site of the oxygen atom. The results suggest that forming two siloxanebridges with two different silicon dimer atoms is the favored product. 9.1 A 9.6 C 9.82 . 10.0.33 10.6 10.43 E F G 11.25 1131 J K 11.6 1144  12.1  12.05 12.15 12.16 12.6  Structure Figure 44. Relative energies of 02Si9H12 isomers. A) DB3(OSiSiOd). B) DB3(SiSiO2d). C) DB3(SiOSiOd). D) MEi(OSiSiO). E) ME3(OxSiSiO). F) MEi(OxSiSiO). G) ME3(SiSiO2). H) ME3(OSiSiO). I) ME3(SiOSiO). J) MEi(SiOSiO). K) MEi(SiSiO2). Energy reference: MEi(SiSi) + 20(3p) = 0.00 eV. Oxygen and silicon atoms are represented by black and gray spheres, respectively. Another substantial effect of strain is the observation that structure ME3(OSiSiO) (Figure 46A) is significantly lower in energy than structure MEi(OSiSiO) (Figure 46E), despite the spin penalty energy that applies to the triplet structure. After taking into account the spin penalty energy, the bond energy model suggests that MEi(OSiSiO) has 1.02 eV more strain than ME3(OSiSiO). In both of these structures, the formation of two p =108.5o Side view Top view Figure 45. Structural information of dangling bond isomers with two oxygen atoms adsorbed on Si(100). A) DB3(OSiSiOd). B) DB3(SiSiO2d). C) DB3(SiOSiOd). Bond lengths are expressed in A. Oxygen atoms are represented with black balls. Silicon and hydrogen atoms are represented in gray and white balls, respectively. Side view Top view Side view Top view Spintriplet Spinsinglet Figure 46. Structural information of 02Si9H12 isomers. Structures are listed by increasing relative stability of the spinsinglet isomers. A) ME3(OSiSiO). B) ME3(OxSiSiO). C) ME3(SiOSiO). D) ME3(SiSiO2). E) MEi(OSiSiO). F) MEi(OxSiSiO). G) MEi(SiOSiO). H) MEi(SiSiO2). Bond lengths are expressed in A. Oxygen atoms are represented with black balls. Silicon and hydrogen atoms are represented in gray and white balls, respectively. siloxane bridges on the same side of the silicon dimer bond shifts the dimer out of its original symmetry plane, without changing its orientation along the [110] direction. However, these structures experience markedly different strain effects due to the different bonding interactions along their respective SiSi dimers. In ME3(OSiSiO), the formation of the two siloxane bridges alleviates strain in the remaining SiSi backbonds, which each shorten by 0.04 A. Aside from being shifted out of plane, the dimer bond remains largely unaltered by the formation of ME3(OSiSiO). These structural changes are predicted to relieve about 0.70 eV of strain from the cluster, and result in considerable stabilization of the structure, making it the most energetically favorable of all the structures with four Si+1 species. In contrast, the formation of the siloxane bridges significantly weakens the w interaction across the silicon dimer in MEi(OSiSiO), causing the dimer bond to stretch by 0.16 A, resulting in an 0.33 eV increase in strain relative to the singlet bare cluster, according to the bond energy analysis. 4.3.3 Structures with Three Adsorbed Oxygen Atoms (O3Si9H12) Figure 47 shows relative energies of isomers with three oxygen atoms and table 44 shows the different contributions to the relative energies of the structures, as determined from the bond energy model. As in the previous two cases discussed above, the energies of the isomers are defined with respect to the reaction described by Equation 46. Si9H12(Ms=l) + 30(3P)  03Si9H12 (Ms=1 or 3) (46) We calculated the relative energies and structural characteristics of the 03Si9H12 isomers (Figures 48 and Figure 49). We found that the relative energy of the structures is determined by a combination of effects, including strain and suboxide formation. 