Copyright by Sangheon Lee 2010 The Dissertation Committee for Sangheon Lee certifies that this is the approved version of the following dissertation: First Principles-based Atomistic Modeling of the Structural Properties of Silicon-Oxide Nanomaterials Committee: Gyeong S. Hwang, Supervisor James R. Chelikowsky Sanjay K. Banerjee John G. Ekerdt Charles B. Mullins First Principles-based Atomistic Modeling of the Structural Properties of Silicon-Oxide Nanomaterials by Sangheon Lee, B.S. Dissertation Presented to the Faculty of the Graduate School of The University of Texas at Austin in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy The University of Texas at Austin August 2010 Dedicated to my parents, Jaewoong Lee and Oghee Kim, and my brother Seungyeol Lee. v Acknowledgements First of all, I would like to thank my advisor, Professor Gyeong S. Hwang. While supporting me for the past six years, he gave me many important research opportunities, carefully advised me, and patiently waited for my progress. Next, I would like to thank former members in the Hwang research group, who include Taras A. Kirichenko, Devina Pillay, Decai Yu, Scott A. Harrison, Jason Kenney, Soohwan Lee, Dr. Yun Wang, and Dr. Chin-Lung Kuo. Their pioneering works in the developing research group provided me with a solid basis which enabled me to perform more advanced researches. I would like to thank current students in the Hwang research group, who include Hyun Woo Kim, Eunsu Paek, Kyoung Eun Kweon, John A. Stephens, Hyung Chul Ham, Chia-Yun Chou, Yongjin Lee, and Robert J. Bondi. This dissertation represents a team effort made by these students. In particular, Robert and I worked in close collaboration in a very complementary fashion, and we were able to successfully perform the assigned projects and publish the important findings. It is a great honor for me to have Professor Chelikowsky, Professor Banerjee, Professor Ekerdt, and Professor Mullins as my dissertation committee members. When I was accepted to the UT-Austin, Professor Ekerdt was the department chair and Professor Mullins was the graduate advisor. I would like to appreciate them for giving me a chance to pursue a PhD in the world best chemical engineering department. It was a turning point in my life. I would like to thank Tim Guinn for safely managing the UNIX/Linux system. I would also like to thank the Texas Advanced Computing Center for use of their computing resources. Lastly, I am grateful for support from the Donald D. Harrington Graduate Fellows Program. vi First Principles-based Atomistic Modeling of the Structural Properties of Silicon-Oxide Nanomaterials Publication No._____________ Sangheon Lee, Ph.D. The University of Texas at Austin, 2010 Supervisor: Gyeong S. Hwang We have developed continuous random network (CRN) model based Metropolis Monte Carlo simulation tools which are capable of predicting the structural properties of amorphous semiconductor and oxide materials as well as their interface. To bolster the reliability of the CRN model, we have developed force fields based on gradient corrected density functional theory (DFT) calculations. Our in-house CRN-MMC tools have been massively parallelized, which allows us to create fairly large model structures within a reasonable computational time. Using the integrated CRN-MMC tools, we have elucidated the complex growth and structure of self-interstitial and vacancy clusters in silicon and the effect of strain on the structure and stability of the defect clusters. Our work for vacancy clusters suggests that small vacancy defects (Vn) exclusively favor fourfold-coordination thermodynamically with no significant kinetic limitation rather than void-like structure formation, which has widely vii been adapted to explain the behavior and properties of vacancy defects. Our results also highlight the identification of stable high-symmetry fourfold-coordinated V12 and V32 clusters that could be expected to exist to a large extent in a vacancy rich region although its direct characterization appears impractical at present. Our work for self-interstitial clusters (In) provides the first theoretical support for earlier experiments which suggest a shape transition from compact to elongated structures around n = 10. When the cluster size is smaller than 10, the stable I4 and I8 compact clusters are found to inhibit the formation of elongated defects, whereas the newly discovered fourfold-coordinated I12 state is found to serve as an effective nucleation center for large extended defects. Our CRN-MMC approach also enabled us to elucidate the underlying mechanisms of synthesis and manipulation of Si rich insulators as well as the fundamental understanding of the relationship between the atomic structure and properties. We developed a valence force field based on a modified Keating model for the structure and energetics of amorphous Si rich oxide (a-SiOx, 0  x  2) materials. In particular, our work emphasizes the importance of correctly describing the wide Si-O-Si angle distribution. Our work also suggests that the relative rigidity between Si and SiO2 matrices is critical in determination of the Si/SiO2 interface structure. The present potential model coupled with the CRN-MMC method can be used to create structural models (free of coordination defects) for complex a-SiOx-based materials, which will further allow thorough studies of the properties of these materials. viii Table of Contents List of Tables xi List of Figures xii Chapter 1: Introduction 1 Chapter 2: Growth of Small Self-interstitial Clusters 7 2.1 Introduction ........................................................................................................7 2.2 Computational Approach ...................................................................................9 2.3 Atomic Structure of Compact Self-Interstitial Clusters (In, n ≤ 10) ................19 2.4 Stability of Compact Self-Interstitial Clusters .................................................20 2.5 Lowest-energy Small Compact Clusters (In, 5 ≤ n ≤ 12) .................................24 2.6 Compact-to-elongated Transition ....................................................................29 2.7 Growth to {311} Extended Defects .................................................................35 2.8 Summary ..........................................................................................................47 Chapter 3: Formation and Structure of Vacancy Defects in Silicon 50 3.1 Introduction ......................................................................................................50 3.2 Computational Approach .................................................................................52 3.3 Determination of Fourfold-coordinated Vacancy Clusters ..............................53 3.4 Relative Stability between Fourfold-coordinated and PHR Vacancy Defects (Vn, n ≤ 48) ..........................................................................................61 3.5 Extremely Stable FC V12 and V32 .....................................................................68 3.6 Kinetic Effects on the Formation of Fourfold-coordinated Vacancy Defects .............................................................................................................72 3.7 Summary ..........................................................................................................80 Chapter 4: Strain Effect on Structure and Stability of Native Defects in Silicon 83 4.1 Introduction ......................................................................................................83 4.2 Computational Approach .................................................................................84 4.3 Biaxial Strain Effects on the Structure and Stability of Self-Interstitial Clusters in Silicon ............................................................................................86 ix 4.4 Strain Effects on the Stability and Structure of Vacancy Clusters in Si ..........95 Chapter 5: Ab Initio Parameterized Valence Force Field for the Structure and Stability of Amorphous SiOx (0 ≤ x ≤ 2) Materials 99 5.1 Introduction ......................................................................................................99 5.2 Computational Approach ...............................................................................103 5.2.1 Valence Force Field Model for a-SiOx (0  x  2) ........................................103 5.2.2 Determination of Force Field Parameters ......................................................106 5.2.3 Metropolis Monte Carlo Simulations.............................................................113 5.2.4 Density-functional-theory Calculations .........................................................115 5.3 Energetics of Amorphous SiOx: Force Field Models ....................................115 5.4 Phase Separation: a-Si Cluster Embedded in a-SiO2 Matrix ........................121 5.5 Mechanical Properties ....................................................................................128 5.6 Summary ........................................................................................................132 Chapter 5: Summary 134 Appendix A: Summary of Publications 137 A.1 Structure and stability of small compact self-interstitial clusters in crystalline silicon ...........................................................................................137 A.2 Growth and shape transition of small silicon self-interstitial clusters ...........137 A.3 Integrated atomistic modeling of interstitial defect growth in silicon ...........138 A.4 Biaxial strain effects on the structure and stability of self-interstitial clusters in silicon............................................................................................138 A.5 Theoretical characterization of silicon self-interstitial clusters in uniform strain fields .....................................................................................................139 A.6 Prediction of the formulation of stable periodic self-interstitial cluster chains [(I4)m, m = 1-4] in Si under biaxial strain ............................................139 A.7 Theoretical determination of stable fourfold-coordinated vacancy clusters in silicon .........................................................................................................140 A.8 Formation and structure of vacancy defects in silicon: Combined Metropolis Monte Carlo, tight-binding molecular dynamics, and density functional theory calculations ........................................................................140 A.9 Strain effects on the stability and structure of vacancy clusters in Si: A first-principles study ......................................................................................141 x A.10 On the origin of Si nanocrystal formation in a Si suboxide matrix ...............142 A.11 Strain-induced formation of surface defects in amorphous silica: A theoretical prediction .....................................................................................143 A.12 Structure and dynamics of silicon-oxygen pairs and their role in silicon self-diffusion in amorphous silica ..................................................................144 A.13 First-principles study of the mechanical and optical properties of amorphous hydrogenated silicon and silicon-rich silicon oxide ....................145 Appendix B: CRN-MMC Simulation 147 B.1 Metropolis Sampling ......................................................................................147 B.2 Bond-switching Method.................................................................................147 Appendix C: Additional Applications of CRN-MMC Simulation 148 C.1 c-Si/a-SiO2 Interface ......................................................................................148 C.2 Hydrogenated Amorphous Silicon (a-Si:H) ..................................................148 C.3 Oxidized Silicon Nanowire ............................................................................149 C.4 Oxidized Silicon Quantum Dot ......................................................................150 C.5 Oxidized Embedded Silicon Nanoparticle .....................................................151 C.6 Fullerene ........................................................................................................152 C.6.1 Buckyball .......................................................................................................152 C.6.2 Topological Defects in Graphene and Carbon Nanotube ..............................153 Bibliography 154 Vita 166   xi List of Tables  2.1: Predicted values for the formation energies per interstitial [Êf(n)] of small self-interstitial clusters (In, n = 1 – 10) from DFT-GGA/LDA and Keating- like potential (KT) calculations as indicated. The supercells α, β, and γ consists of 192 + n, 400 + n, and 480 + n atoms, respectively, where n is the number of interstitials. The corresponding defect symmetries are also summarized. .......................................................................................................... 22 5.1: Force constants for KT(LH) potential. Force constants for other KT potentials are also summarized for reference. ..................................................... 106 5.2: Calculated equilibrium distances of Si-Si and Si-O bonds from bulk structures. The values from the cluster calculations are shown in (). ................ 107 5. 3: Si suboxide statistics averaged over four independent structures of a-SiO0.5 (64 Si and 32 O atoms), a-SiO1.0 (64 Si and 64 O atoms), and a-SiO1.5 (64 Si and 96 O atoms) used in Figure 5.4 and 5.5(a). These structures were constructed from CRN-MMC simulations based on the KT(LH) potential excluding suboxide penalty energies. ................................................................. 116 5.4: Si suboxide statistics averaged over four independent structures of a-SiO0.5 (64 Si and 32 O atoms), a-SiO1.0 (64 Si and 64 O atoms), and a-SiO1.5 (64 Si and 96 O atoms) used in Fig. 5.5(b). These structures were constructed from CRN-MMC simulations based on the KT(LH) potential including suboxide penalty energies. .................................................................................. 120 5.5: Si suboxide statistics averaged over five independent KT(LH) and KT(TT) models of np-Si/a-SiO2 (480 Si and 720 O atoms). ............................................ 123 5.6: Computed average mechanical properties based on KT(LH) and KT(TT) potentials [1] for ten independent a-Si (216 Si atoms) and a-SiO2 (216 Si and 432 O atoms) structures. Strain was applied during mechanical property calculations using the same KT(LH) and KT(TT) potentials used during the initial CRN-MMC simulations. Relevant experimental data is also summarized for comparison. ....................................................................... 129 xii List of Figures 2.1: Predicted lowest energy configurations for (a) I1, (b) I2, (c) I3, and (d) I4. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated [5]... ............. 12 2.2: Predicted minimum-energy penta-interstitial (I5) structure with C1 symmetry from three different perspectives, as indicated. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white), and the fifth interstitial is also indicated [5]... .......................................... 14 2.3: Predicted two nearly degenerate configurations for the hexa-interstitial cluster (I6) from three different perspectives, as indicated (see upper insets in Figure 2.2). Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). Also, the fifth and sixth interstitials captured by the I4 cluster are also indicated in (b). The corresponding defect symmetry is indicated [5]... ........................................................................ 15 2.4: Predicted three local minimum configurations for the hepta-interstitial cluster (I7) from three different perspectives, as indicated (see upper insets in Figure 2.2). The I7a structure is about 0.39 eV and 0.75 eV more favorable than the I7b and I7c structures, respectively. All structures are fourfold-coordinated. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The corresponding defect symmetries are also indicated [5]... ...................................................................... 17 2.5: Predicted minimum energy configuration with Ci symmetry for the octa- interstitial cluster from three different perspectives, as indicated (see upper insets in Figure 2.2). Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white) [5]... .................................................. 18 2.6: Predicted two local minimum configurations for the ennea-interstitial cluster (I9) from three different perspectives, as indicated (see upper insets in Figure 2.2). The I9a structure is about 0.09 eV more favorable than the I9b structure. In (a) the ninth interstitial captured by the I8 cluster is indicated. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The corresponding defect symmetries are also indicated [5]... ....................................................................................................... 18 2.7: Predicted two local minimum configurations for the deca-interstitial cluster (I10) from three different perspectives, as indicated (see upper insets in Figure 2.2). The I10a structure is about 0.3 eV more favorable than the I10b structure. Both structures are fourfold-coordinated. Grey (gold) balls xiii indicated more distorted atoms than the rest of the lattice atoms (in white). The corresponding defect symmetries are also indicated [5]... ............................. 20 2.8: (a) Predicted formation energies per interstitial (Êf) of small compact interstitial clusters from our DFT-GGA calculations and (b) their approximate capture radii for mobile interstitials. The capture radius is approximated as rc = (3V/4π)1/3, where V is a spherical-like volume that corresponds to the number of Si lattice atoms (excluding interstitials) whose strain energies are greater than a given value (as indicated) [5]... ............. 23 2.9: Predicted minimum-energy configurations for (a) penta-, (b) hexa-, (c) hepta-, (d) octa-, (e) ennea-, (f) deca-, (g) hendeca-, (h) dodeca-interstitial defects in Si. For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated. Grey (gold) balls represent more distorted atoms than the rest of the lattice atoms (in white) [6]... ............................................................................................ 25 2.10: Compact, minimum energy configurations with corresponding defect symmetries for (a) hendeca- and (b) dodeca-interstitial defects in Si. For each defect, the top and bottom panels show two different views, as indicated. Light grey (gold) wireframe represents the bulk Si lattice. Dark grey spheres represent interstitial atoms and their highly strained neighbors [7].... ..................................................................................................... 28 2.11: Calculated formation energies per interstitial (Êf) for interstitial clusters shown in Figure 2.1, 2.9, and F#_013 as well as chain-like elongated structures (indicated as ‘Chain-like’) as a function of cluster size (n). For the chain-like case, the atomic structures of Kim et al. [51] were recalculated within DFT-GGA. To minimize possible interactions between a defect and its periodic images, we evaluated the formation energies by scaling the supercell size as follows (‘This work’): 192+n (I1– I2), 400+n (I3-I6), 480+n (I7-I10), 560+n (I11-I13), and 672+n (I14-I16) atom supercells, where n is the number of interstitials. Likewise, the following supercells sizes were used (‘Chain-like’): 400+n (I3-I4), 480+n (I5-I6), 560+n (I7-I8), and 640+n (I9-I10) atom supercells. For n > 10, the chain- like cluster formation energy is forecasted using Êf(n) = {6.83 + (n − 2)×1.28}/n, see also Figure 2.17 [7]... .................................................................. 30 2.12: Configuration for the chain-like, elongated hexa-interstitial defect with C1 symmetry in Si. The top and bottom panels show two different views, as indicated. In both views the strain energy distribution around the cluster is binned into three different levels [high (≥ 0.25eV) = dark grey (red) balls, medium (0.15 – 0.25 eV) = grey balls, and low (< 0.15 eV) = light grey xiv (gold) wireframe, based on strain energy values from force field calculations] [7]... ................................................................................................. 31 2.13: Minimum energy structures for I12 through I16 with the I12-like core, which exhibit preferred growth along the [110] direction when n > 10. Configurations for I12 through I16 are depicted with dark grey balls showing highly distorted atoms and light grey (gold) wireframe representing the bulk lattice. Additionally, the upper right panel provides an atomic-level strain distribution profile for the I12 cluster and nearby bulk Si based on average bond length. Each atom is assigned an average bond length based on the four bonds formed with its nearest neighbors and is gradient-shaded accordingly. The color spectrum shifts from red to white to blue (RWB) as the average bond lengths shift from compressive to strain-free to tensile. The RWB spectrum range displayed here covers the Si-Si DFT equilibrium bond length of 2.36  0.07Å [7]... ............................. 34 2.14: Three different defect core structures which are infinitely long in the [110] direction: (a) I12-like (as shown for the FC I12 cluster identified in this work); (b) chain-like; and (c) {311}. For each core, the left and right panels show two different views, as indicated. Grey balls indicate more distorted atoms than the rest of the light grey (gold) wireframe lattice [7]... ....... 37 2.15: Minimum energy configurations for {311} core structures of (a) I12 and (b) I16 clusters in Si. For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated. Dark grey (red) balls indicate highly distorted atoms in the (110) interfaces, grey balls represent moderately distorted atoms, and light grey (gold) wireframe represents the bulk Si lattice [7]... .................................... 38 2.16: Minimum energy configurations for the (110) edges of (a) the elongated I12-like, (b) chain-like, and (c) {311} core configurations. Dark grey (red) balls indicate highly distorted atoms contained in the (110) interfaces, grey balls represent moderately distorted atoms, and light grey (gold) wireframe represents the bulk Si lattice [7]... ......................................................................... 39 2.17: Predicted formation energies per interstitial of I12-like, chain-like, and {311} core configuration clusters as a function of size (n) using Ef(n) = {ܧ௙௘௡ௗ + (n – n′)Ê௙௖}/n, where ܧ௙௘௡ௗ is the (110) interface energy, Ê௙௖ is the core energy per interstitial for the semi-infinite structure, and n′ is the number of interface atoms. Here, the interface energies were obtained by numerically fitting DFT-GGA values for the formation energies of small clusters: I12 and I16 were used for the I12-like and {311} cases, while I8, I9 and I10 represented the chain-like, elongated structures. The predicted xv values are based on Ê௙௖ = 1.29, 1.28, and 1.11 eV and ܧ௙௘௡ௗ = 8.74, 6.83, and 11.62 eV for the I12-like, chain-like, and {311} cases, respectively. For the number of interface atoms, n′ = 4 was used for the I12-like and {311} structures, while n′ = 2 was used for the chain-like, elongated structure [7]... ....................................................................................................................... 42 2.18: Block diagram for a typical {311} extended defect structure. Characteristic dimensions observed experimentally and alignment of the defect within the bulk Si lattice is shown. From the computational results obtained, incremental growth along <233> (small block extension) followed by more favorable growth along <110> (block arrows) is anticipated [7]... .................................................................................................... 43 2.19: (a) An atomic-level strain distribution profile throughout the {311} core structure and nearby bulk material based on average bond length. Each atom is assigned an average bond length based on the four bonds formed with its nearest neighbors and is gradient-shaded accordingly. The color spectrum shifts from red to white to blue (RWB) as the average bond lengths shift from compressive to strain-free to tensile. The RWB spectrum range displayed here covers the Si-Si DFT equilibrium bond length of 2.36  0.07Å. Additionally, the constituent rings of the model {311} core have been shaded grey to help visualize the defect cross- section in the bulk. The “I” (interstitial) and “E” (edge) ring annotations reflect the nomenclature of Kohyama and Takeda [82] to generalize the periodicity seen in {311} planar structures. From left to right, this structure would be labeled “EIIE”. The “I” repeating unit contains two net interstitials per period along <110>. From the strain profile, the structure core is compressively strained, the {311} interface planes are under slight compression, and the {233} interface planes are under tensile strain. Minimum energy configurations for (b) “{311} core + I3” and (c) “{311} core + I4” structures. In both cases, the interstitials appended themselves to the tensile {233} interface to help minimize excess strain energy. The appended tri-interstitial and tetra-interstitial configurations (bright red) show strong resemblance to the ground state I3 and I4 compact configurations of Figure 2.1 [7]... ......................................................................... 44 3.1: Predicted fourfold-coordinated configurations for (a) V3, (b) V4, (c) V5, and (d) V6 in Si. Grey (gold) balls indicate more distorted atoms than the rest of the lattice atoms (in white). For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated... ..................................................................................................... 56 xvi 3.2: Predicted fourfold-coordinated configurations for (a) V7, (b) V8, (c) V9, and (d) V10 in Si. Grey (gold) balls indicate more distorted atoms than the rest of the lattice atoms (in white). For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated... ..................................................................................................... 57 3.3: Predicted fourfold-coordinated configurations for (a) V11, (b) V12, (c) V13, and (d) V14 in Si. Grey (gold) balls indicate more distorted atoms than the rest of the lattice atoms (in white). For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated... ....................................................................................... 58 3.4: Predicted fourfold-coordinated configurations for (a) V15, (b) V16, (c) V17, and (d) V18 in Si. Grey (gold) balls indicate more distorted atoms than the rest of the lattice atoms (in white). For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated... ....................................................................................... 59 3.5: Distributions of bond length (left panels) and bond angle (right panels) deviations associated with the distorted atoms as shown in Figure 3.1, 3.2, 3.3, and 3.4. The dotted lines indicate the equilibrium values of r0 = 2.365 Å and 0 = 109.5 ° [8]... .............................................................................. 60 3.6: (a) Calculated formation energies per vacancy (Êf) and (b) binding energies of vacancy clusters as a function of cluster size (n) for both fourfold-coordinated (indicated as “Fourfold”) and PHR (“PHR”) configurations. For the PHR case, the atomic structures from previous studies [90,91] were recalculated within DFT-GGA. To minimize possible interactions between a defect and its periodic images, we carefully evaluated the formation energies by changing the supercell size; 480–n and 576–n atom supercells, where n is the number of vacancies, were used for V1-V14 and V15-V18, respectively. The inset shows a variation in the total energy difference (∆E in eV) between the fourfold- coordinated and PHR cases [8]... .......................................................................... 62 3.7: (a) DFT-predicted formation energies per vacancy (Êf) of vacancy clusters as a function of cluster size (n) for both fourfold-coordinated (indicated as FC) and “part of hexagonal ring” (PHR) configurations, together with adjusted values (FC*) for the FC configurations. For the PHR case of n > 18, the DFT formation energies per vacancy are well fitted with the power function of (−0.82 + 3.68n2/3)/n, as indicated with a dotted line. (b) Variation in the total energy difference between the FC and PHR structures from DFT calculations (indicated as ΔE), together with adjusted values (ΔE*). Refer to the text for the adjustment in the FC defect xvii formation energies using FF calculations, which is associated with incomplete relaxation of the defect system due to limited supercell size used for DFT calculations [9]... ............................................................................ 65 3.8: Comparison of the FC cluster formation energies (Ef) from between DFT (indicated as DFT) and FF (as detailed in Section 3.2) calculations. The inset shows the Ef variations of FC V12 and V32 in terms of supercell size; here N indicates the number of Si lattice atoms prior to vacancy defect creation [9]... ......................................................................................................... 67 3.9: (a) Ball-and-stick illustration of the atomic structure of V12, together with Wannier function centers as indicated by small black balls. Grey (gold) balls indicate more distorted atoms than the rest of the lattice atoms (in white), while five dark grey (red) balls indicate a well relaxed tetragon at the defect center. (b) Calculated total density of states (TDOS) for the c-Si system with (upper panel) and without (lower panel) the V12 cluster. Here, EF indicates the Fermi level [8]... ......................................................................... 70 3.10: (a) Core and strained nearest neighbors of FC V32 shown isolated from bulk Si lattice. In (a), we annotate the antiparallel edge type, where member atoms of each chain are distinguished by either “x” or “†” labels (blue or red spheres, respectively). The ten atoms comprising the adamantine cage (AC) in the center of the structure are identified as dark grey spheres. (b) Illustration of a simplified geometric diagram of FC V32 that highlights the S24 group symmetry, AC core, and the two antiparallel (double line) and four parallel (single line) cluster edges. Axis labels denote cluster orientation with respect to the Si lattice. (c) Reorientation of FC V32 to emphasize the parallel edge, where atoms comprising these chains are annotated with “x” labels (blue spheres). (d) Total density of states comparison of bulk Si to V32 near the band gap using 1000 atom basis supercells with the conduction band edges set as references [9]... .............. 71 3.11: Mono-vacancy induced structural rearrangement, together with variation in the formation energy (Ef = nÊf) as a function of annealing time for V3 + V → V4 during TBMD simulations at 1400K. All atomic configurations from the TBMD simulation were optimized within DFT-GGA calculations using one Gamma k-point sampling for 512−n atom supercells and illustrated in the top and middle panels. Grey wire frame represent bulk Si lattice and dark grey spheres represent highly strained FC atoms while brown spheres represent a Si mono-vacancy. Si dangling defects are also depicted by red (black) spheres [9]... .................................................................... 73 3.12: Mono-vacancy induced structural rearrangement, together with variation in the formation energy (Ef = nÊf) as a function of annealing time for V11 + xviii V → V12 during TBMD simulations at 1400K. All atomic configurations from the TBMD simulation were optimized within DFT-GGA calculations using one Gamma k-point sampling for 512−n atom supercells and illustrated in the top and middle panels. Grey wire frame represent bulk Si lattice and dark grey spheres represent highly strained FC atoms. Si dangling defects are also depicted by red (black) spheres [9]... ........................... 74 3.13: Mono-vacancy induced structural rearrangement, together with variation in the formation energy (Ef = nÊf) as a function of annealing time for V4 + V → V5 during TBMD simulations at 1400K. All atomic configurations from the TBMD simulation were optimized within DFT-GGA calculations using one Gamma k-point sampling for 512−n atom supercells and illustrated in the top and middle panels. Grey wire frame represent bulk Si lattice and dark grey spheres represent highly strained FC atoms. Si dangling defects are also depicted by red (black) spheres [9]... ........................... 75 3.14: Structural evolution of amorphous pockets with deficit of 12 atoms during TBMD simulations at 1400K together with variation in the formation energy (Ef = nÊf) as a function of annealing time. All atomic configurations from the TBMD simulation were optimized within DFT- GGA calculations using one Gamma k-point sampling for 512−n atom supercells. Filled circles illustrate the formation of the lowest-energy FC V12 with C2 symmetry, together with their atomic configurations as indicated. Open circles illustrate the formation of a PHR V12 structure, together with their atomic configurations as indicated; however, the PHR defect formation would be rather unlikely in a typical annealing process, as discussed in the text. Grey wire frame represent bulk Si lattice and dark grey spheres represent highly strained FC atoms. Si dangling and floating bonds are depicted by red (black) and blue (black, pointed also by the arrow) spheres, respectively [9]... ......................................................................... 76 4.1: Tensile biaxial stress/strain interaction in our model Si supercell. In the figure, applied tensile stress, σ║, in the plane of the substrate acts equally in all directions as shown by block arrows and produces a tensile strain. In response, the lattice contracts in the out-of-plane direction as shown by the solid black arrows. Under compressive strain conditions, the directions of all arrows are inverted [10]... ................................................................................ 85 4.2: Predicted lowest energy configurations for I4 under 4% compressive (top), strain-free (middle, I4g), and 4% tensile (bottom, I44%t) biaxial strain. Biaxial strain is applied to the supercell as shown in Figure 4.1. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated... ................... 87 xix 4.3: Predicted lowest energy configurations for I3 under 4% compressive (top, I34%c), strain-free (middle, I3g), and 4% tensile (bottom, I34%t) biaxial strain. Biaxial strain is applied to the supercell as shown in Figure 4.1. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated... ................... 89 4.4: Predicted lowest energy configurations for I7 under 4% compressive (top, I74%c), strain-free (middle, I7g), and 4% tensile (bottom, I74%t) biaxial strain. Biaxial strain is applied to the supercell as shown in Figure 4.1. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated... ................... 90 4.5: Predicted lowest energy configurations for I8 under 4% compressive (top), strain-free (middle, I8g), and 4% tensile (bottom, I84%t) biaxial strain. Biaxial strain is applied to the supercell as shown in Figure 4.1. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated... ................... 91 4.6: Formation energy response per interstitial compared for the two relevant orientations of I4 using 256+n supercells. Perspective views along the [001] are shown for reference along with corresponding orientations of the S4 rotation-reflection axes. Light gray (gold) wireframe represents the bulk Si lattice. Dark gray spheres denote interstitials and their highly strained neighbors [11]... ...................................................................................... 93 4.7: Predicted lowest energy configurations for I12 under 4% compressive (top, I124%c), strain-free (middle, I12g), and 4% tensile (bottom, I124%t) biaxial strain. Biaxial strain is applied to the supercell as shown in Figure 4.1. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated... ................... 94 4.8: Formation energy per vacancy as a function of biaxial strain for the strain- free, ground-state FC V3 structure and a FC trivacancy configuration (V3t) identified by formation under tensile biaxial strain conditions. Inset cluster configurations are shown along [11ത0]. Light gray (gold) wireframe represents bulk crystalline Si. Dark gray spheres represent highly strained atoms neighboring the vacancy clusters [13]... ..................................................... 97 4.9: Formation energy per vacancy as a function of biaxial strain for the strain- free, ground-state FC V5 structure and a FC pentavacancy configuration (V5t) identified by formation under tensile biaxial strain conditions. These structures are approximately degenerate under strain-free conditions. Inset cluster configurations are shown along [11ത0]. Light gray (gold) wireframe xx represents bulk crystalline Si. Dark gray spheres represent highly strained atoms neighboring the vacancy clusters [13]... ..................................................... 98 5.1: Representative cluster models used for calculating Si-Si and Si-O bond lengths. Cluster models for other oxidation states (see Table 5.2) are obtained by adjusting the number of O atoms from these models. For calculating Si-O bond lengths, the two Si atoms neighboring the central O atom retain the same oxidation state... ................................................................ 108 5.2: Variations (ΔE) in total energies (from DFT, KT(LH), and KT(TT) calculations) per bond of (a) c-Si (with 8 atoms) and (b) c-SiO2 (β- cristobalite with 8 SiO2 units) as a function of the ratio (L/L0) of the strained lattice constant (L) to the equilibrium lattice constant (L0)... ................ 109 5.3: Variations (ΔE) in total energies (from DFT, KT(LH), and KT(TT) calculations) of the cluster model (inset) as a function of Si-O-Si bond angle (θ)... ........................................................................................................... 111 5.4: Relative total energies per Si atom (ΔÊtotal) from KT(LH), KT(TT), WT1, WT2, and DFT calculations for a-SiOx (x = 0, 0.5, 1.0, 1.5, and 2.0) (64 Si and 64x O atoms) structures constructed from CRN-MMC simulations based on the KT(LH) potential without suboxide penalty energies. For each x, four independent structures are represented. For x = 0.5, 1.0, and 1.5, the distributions of Si oxidation states are summarized in Table 5.3... ....... 117 5.5: Average relative strain (ΔÊstrain) and suboxide (ΔÊsubox, inset) energies per Si based on KT(LH), KT(TT), and DFT calculations for a-SiOx (x = 0, 0.5, 1.0, 1.5, and 2.0) models constructed from CRN-MMC simulations based on the KT(LH) potential (a) excluding and (b) including suboxide penalty energies. For each x, four independent structures are represented. For x = 0.5, 1.0, and 1.5, the distributions of Si oxidation states for (a) and (b) are summarized in Tables 5.3 and 5.4, respectively... .............................................. 119 5.6: Atomic configurations for (a) KT(LH) and (b) KT(TT) models for the a-Si cluster embedded in a-SiO2 matrix (np-Si/a-SiO2). Grey wireframe represents O atoms and Si4+ states that comprise the a-SiO2 phase. Yellow, blue, red, and grey balls represent Si0, Si1+, Si2+, and Si3+ states, respectively... ...................................................................................................... 122 5.7: Profiles of (a) strain (ΔÊstrain) and (b) suboxide (ΔÊsubox) energies per Si along radial directions from cluster centers of KT(LH) and KT(TT) models for np-Si/a-SiO2 (480 Si and 720 O atoms). The cluster center is defined as the center of mass of Si0 atoms. The nominal interface radius, r0, is defined in the text. Each data point represents the average value xxi within a given concentric spherical shell (2 Å thick) sampled over four independent structures. The two solid, horizontal lines depict the calculated strain energies for bulk a-Si and a-SiO2 with 216 Si and SiO2 units, respectively. All energies are calculated with the KT(LH) potential............................................................................................................... 124 5.8: Profiles of (a) Si1+, (b) Si2+, and (c) Si3+ oxidation state distributions along radial directions from cluster centers of KT(LH) and KT(TT) models of np-Si/a-SiO2 (480 Si and 720 O atoms). The cluster center is defined as the center of mass of Si0 atoms. The nominal interface radius, r0, is defined in the text. Each data point represents the average value within a given concentric spherical shell (2 Å thick) sampled over four independent structures... .......................................................................................................... 126 5.9: Ring-size distributions for the (a) entire, (b) a-Si, and (c) a-SiO2 regions of KT(LH) and KT(TT) models of np-Si/a-SiO2 (480 Si and 720 O atoms)... ....... 127 A.1: Series of snapshots from KMC simulations of phase separation in Si suboxide with the initial Si supersaturation of (a) 10%, (b) 20%, and (c) 30% at 1100 ºC. Here, only Si0 atoms are displayed. The box size is 8.1×8.1×8.1 nm3 [189].. ...................................................................................... 140 A.2: Variation in the relative total energy (in eV) of a thin a-SiO2 slab during ab-initio MD relaxation at 1200 K as a function of annealing time (in ps). The strain induced formation of surface edge-sharing tetrahedron and silanone defects is also indicated, together with their configurations (in the insets) [17]... ....................................................................................................... 141 A.3: Atomic structure of a SiO pair in a-SiO2. The enclosed shows SiO pair insertion into the a-SiO2 matrix, resulting in the formation of a divalent Si atom. The small dark (red) and big gray (pink) balls represent O and Si atoms, respectively [18]... ................................................................................... 142 A.4: Imaginary (a) and real (b) parts of the complex dielectric function spectra for a-SiOx computed from DFT-GGA calculations. The spectra are annotated according to O stoichiometry (x) relative to Si. For each composition, only the lowest energy structure is represented in the spectra [19]. ..................................................................................................................... 144 B.1: Schematic illustration of bond switching event in which bonds AB and CD are broken and bonds AC and BD are newly formed... ...................................... 145 C.1: Planar Si/SiO2 interface configurations for (a) Si(100)/a-SiO2, (b) Si(110)/a-SiO2, and (c) Si(111)/a-SiO2. The small gray balls represent O xxii atoms, and big gold (Si0), green (Si1+), brown (Si2+), red (Si3+), and white (Si4+) balls represent Si atoms with different oxidation states... ......................... 146 C.2: Amorphous configurations of a-Si:H (64 Si and 4 H atoms). Silicon and hydrogen atoms are represented by gold (large) and blue (small) spheres, respectively... ...................................................................................................... 146 C.3: Oxidized (a) Si<100>, (b) Si<110>, and (c) Si<111> nanowire. The small gray balls represent O atoms, and big gold (Si0), green (Si1+), brown (Si2+), red (Si3+), and white (Si4+) balls represent Si atoms with different oxidation states.................................................................................................... 147 C.4: Oxidized Si quantum dot and cross sections. The small gray balls represent O atoms, and big gold (Si0), green (Si1+), brown (Si2+), red (Si3+), and white (Si4+) balls represent Si atoms with different oxidation states... ........ 148 C.5: Oxide embedded (a) amorphous Si cluster and (b) Si nanocrystal. The small gray balls represent O atoms, and big gold (Si0), green (Si1+), brown (Si2+), red (Si3+), and white (Si4+) balls represent Si atoms with different oxidation states.................................................................................................... 149 C.6: Comparison of the strain energies (ΔEstrain) for Cn buckyball structures between our force field (indicated as This work) and DFT (indicated as DFT). All buckyball structures are from CRN-MMC simulations... ................. 150 C.7: Topological defects: (a) V12 in graphene and (b) V4 carbon nanotube from CRN-MMC simulations... ................................................................................... 151 1 Chapter 1 Introduction The ability to control the synthesis and manipulation of nanostructured materials with atomic scale precision, together with an accurate assessment of the structure- property relationships, offers enormous opportunities for the development of a variety of novel electronic, chemical, and biological devices. Significant theoretical and experimental efforts have recently been undertaken to increase the fundamental understanding of the synthesis, manipulation and characterization, and in turn experimental control of the structural and electronic properties of nanomaterial-based systems. For most nanomaterials, however the fundamental mechanics of synthesis and manipulation and their structures are poorly known, hampering the development of nanotechnology-based novel devices. Our work intends to develop the theoretical foundations for the synthesis, manipulation and characterization of low dimensional Si and Ge structures encapsulated or embedded in an insulator including: i) two-dimensional (2D) Si-SiO2 multi quantum wells, ii) one-dimensional (1-D) Si nanowires encapsulated within SiO2, and iii) zero- dimensional (0-D) Si and Ge nanocrystals embedded in SiO2. The Si-Ge-Insulator nanosystems are receiving considerable attention because of not only their interesting fundamental properties but also promising applications in electronics and photonics. The electronic/optical properties of the Si-Ge-Insulator nanosystems are determined by a complex combination of i) the size and dimensionality of the Si/Ge 2 nanostructures and ii) the atomic structure, bonding, composition, and doping at Si/Ge- Insulator, Si/Ge bulk, interfaces, and insulator matrices. Experiments may provide many clues to the atomic properties and behaviors, but their interpretations often remain controversial due largely to difficulties in direct measurement. While current experimental techniques are still limited to providing complementary atomic-level, real space information, comprehensive multiscale modeling based on first principles quantum mechanics can contribute greatly to i) elucidating the underlying mechanisms of synthesis, manipulation, and characterization, ii) understanding the structure and structure-property relationships, and iii) improving existing (or developing new) process technologies. This work is grouped into three thrusts including: (1) Defect engineering in Si-Ge-Insulator nanosystems, (2) Oxidation of Si nanostructures and Si-oxide interfacial properties, and (3) Synthesis of Si and Ge nanoparticles in oxide matrices. To address these problems theoretically, we developed continuous random network (CRN) model based Metropolis Monte Carlo simulation tools. The CRN based MMC (CRN-MMC) approach has been successfully used to determine the fully relaxed structure of fully coordinated amorphous materials and their interfaces [1-4], while molecular dynamics simulations alone may not always guarantee the construction of thermally equilibrated structures due to their intrinsic time scale limitations. To bolster the reliability of the CRN model, we carefully optimized the potential parameters based on gradient corrected density functional theory (DFT) calculations. Furthermore, our in- 3 house CRN-MMC tools have been massively parallelized, which allows us to create fairly large model structures within a reasonable computational time. First our CRN-MMC approach was applied to address long time unsolved problems for the fundamental behavior of native defects in crystalline silicon (c-Si) [5- 14]. Vacancies (V) and self-interstitials (I) are two fundamental native defects in c-Si. It is well known that these native defects usually grow into larger open volume defects or extended defects via small vacancy or interstitial clusters. These native defects not only exhibit interesting physics of their own but also greatly alter material properties. They also play an important role in dopant transient enhanced diffusion and electrical activation/deactivation. However, the structure and stability of larger defect clusters were still uncertain due largely to their small size and the possible complexity in their geometries, which hamper both the current experimental and conventional theoretical studies. Our CRN-MMC simulations enabled us to theoretically determine the structure and stability of small vacancy and self-interstitial clusters. Our results highlight the identification of stable high-symmetry fourfold-coordinated i) V12 and V32 clusters which could be expected to exist to a large extent in a vacancy rich region although its direct characterization appears impractical at present and ii) I12 cluster which could be expected to serve as an effective nucleation for larger extended defects. Ultimately these structural models enabled us to reveal the possible mechanisms of how these small native defect clusters grow into larger open volume or extended defects. Recently Yu, Lee, and Hwang [189] have presented the underlying mechanism of Si nanocrystal formation in an amorphous Si suboxide matrix based on extensive first 4 principles based atomistic modeling. The results predict that the formation of oxide embedded Si clusters is primarily attributed to chemical phase separation to Si and SiO2, which is mainly driven by suboxide penalty, with a minor contribution of strain. The phase separation turns out to be primarily controlled by diffusion of O atoms rather than Si atoms. From the kinetic Monte Carlo simulations based on these fundamental findings, the authors have identified two growth mechanisms: “coalescencelike” and “pseudo- Ostwald ripening.” The former is mainly responsible for fast Si cluster growth at the early stage of annealing where the clusters are close each other, while the latter becomes important when the density of clusters is low such that they are separated by large distances. The simulation results also show that the ripening process takes place several orders of magnitude slower than the coalescencelike growth, and the prevailing coalescencelike growth behavior results in a big variation in the Si cluster size with the Si:O ratio. The results agree well with experimental observations that show strong dependence of the cluster size on the initial Si supersaturation as well as rapid formation of Si clusters at the early stage of annealing with very slow ripening. This is the first ever and only theoretical work for the origin of Si nanocrystal formation in a Si suboxide matrix. While the simple but physically sound growth model explains well the intriguing growth behavior of Si nanoclusters in amorphous Si suboxide, further investigations into Si cluster crystallization are required to develop an improved analytical or kinetic model for precise description of the structural evolution of embedded Si nanoclusters. Our CRN-MMC approach also enabled us to elucidate various scientific and technological questions associated with Si-Insulator nanosystems. While recent first 5 principles based molecular dynamics studies provide atomic-level information for a few selected cases of these nanosystems [15-18], their applications on real space systems are limited largely due to the computational time and system size limitations. The CRN- MMC simulations, based on first principles quantum mechanics, can be an alternative to the conventional first principles based molecular dynamics due to the computational efficiency in searching minimum energy structures although it is not as accurate as the first principles based calculations. We have performed an intensive research to develop multiscale modeling strategies and reliable force field models based on first principles quantum mechanics. We observed that the previously introduced Keating-like potential [1] induces a significant deviations from our first principles calculations as the O:Si ratio varies. To fix this problem, we modified the current Keating-like potential and also optimize the bond stretching and bond angle bending force constants based on the first principles calculations. The improved force field calculations significantly increase the reliability of the CRN model with respect to the first principles calculations. Our CRN- MMC simulations enabled us to theoretically determine the structure and stability of technologically-important systems containing Si and/or O, such as amorphous silicon (a- Si), hydrogenated amorphous silicon (a-Si:H), amorphous silica (a-SiO2), planar crystal silicon/oxide interface (c-Si/a-SiO2), and oxide embedded silicon nanoparticle (np-Si/a- SiO2). [19-22] Ultimately these structural models enabled us to elucidate the underlying mechanisms of synthesis and manipulation of Si rich insulators as well as the fundamental understanding of the relationship between the atomic structure and 6 properties. Our results highlight that the relative rigidity between Si and SiO2 phases is critical in determination of the Si/SiO2 interface structure. 7 Chapter 2 Growth of Small Self-interstitial Clusters 2.1 Introduction There have been significant efforts to understand the fundamental behavior of Si interstitial defects created by bombardment with energetic dopant ions, due to their crucial role in defining ultrashallow pn junctions for ever smaller semiconductor device fabrication. Single Si interstitials are highly mobile even at room temperature [23,24]. Hence, in bulk Si, self-interstitials are likely to remain in the form of clusters or interstitial-impurity complexes. The formation and structure of rod-like {311} defects have been well characterized by high-resolution transmission electron microscopy (HRTEM) [25-27]. In addition, a series of recent spectroscopy measurements [28-33] have evidenced existence of small compact self-interstitial clusters before they evolve into larger extended defects. In ultrashallow junction formation with low-energy implanted dopants, such small interstitial clusters are thought to be a main source for free interstitials responsible for dopant transient enhanced diffusion and agglomeration during post-implantation thermal treatment [34-45]. Hence, significant experimental and theoretical efforts [5,46-61] have been made to determine the structure and stability of small self-interstitial clusters as well as their growth to larger extended defects, yet still unclear. 8 Earlier inverse model studies [62,63] based on experimental observations of the spatial and temporal evolution of extended defect distribution in ion-implanted Si suggested the occurrence of structural transition between small compact clusters and larger extended defects when the cluster size is around 10 atoms. This prediction has been supported by a series of low temperature photoluminescence (PL) studies [29,31,32] which show the signatures of small compact clusters of various sizes in ion-implanted Si upon short annealing. In addition, the PL measurements demonstrate changes in the spectra from multiple sharp peaks to broad but distinct signatures during prolonged annealing, ascribed to the shape evolution of self-interstitial clusters from compact to extended forms. While the empirical studies are rather limited to explicitly show the formation and atomic structure of small compact interstitial clusters, a few previous theoretical studies [51,54,61] also predicted the possibility that a chain-like, elongated tri- interstitial cluster may play a role for the compact-to-extended transition. Only a few small compact interstitial clusters (In, n = 2-4) have been unequivocally identified by first principles calculations [48,54,64]. At present, still the atomic structure and stability of larger compact clusters (n  5) are uncertain, due largely to the possible complexity in their geometries. The lack of information hampers revealing how small compact interstitial clusters grow to larger extended defects. In this work, we present an effective computational approach to determine the structure and energetics of interstitial clusters, which combines continuous random network model, tight binding molecular dynamics, and first principles quantum mechanical calculations. Using the new approach, we have identified the minimum- 9 energy configurations and formation energies of small, compact self-interstitial clusters (n  12). We also compare our work with existing theoretical studies and experimental measurements, which shows excellent agreement. Based on these structural models, we present the growth and structural evolution of Si self-interstitial clusters in Si and discuss the growth of extended {311} defects from the small clusters. Such accurate determination of the structure and stability of small interstitial clusters will enable more rigorous examinations of temporal and spatial evolution of interstitial concentration profiles, and in turn interstitial mediated dopant diffusion and clustering during the formation of ultrashallow junctions in Si-based electronic devices. 2.2 Computational Approach Within the continuous random network (CRN) model [65], a disordered structure is generated via a large number of bond transpositions using Metropolis Monte Carlo (MMC) sampling. The CRN based MMC (CRN-MMC) approach has been successfully used to determine the fully relaxed structure of amorphous materials and their interfaces, while molecular dynamics (MD) simulations alone may not always guarantee the construction of thermally equilibrated structures due to their intrinsic time scale limitations. Likewise, we expect that the structure of self-interstitial clusters in Si can also be predicted using CRN-MMC simulations, if all atoms in the clusters are fourfold- coordinated. Indeed, recent theoretical studies [48,54] have demonstrated that such fourfold-coordinated (FC) structures are energetically favored. Our in-house CRN-MMC 10 code has been massively parallelized, which allows us to create fairly large interstitial clusters in crystalline Si within a reasonable computation time. We employed Keating-like potentials [66] which have been proven to be reliable for studying the relaxed structure of disordered Si materials. Within the Keating-like valence force model, the strain energy (Estrain) is defined as: ܧ௦௧௥௔௜௡ ൌ ଵଶ∑ ݇௕ሺܾ௜ െ ܾ଴ሻଶ ൅ ଵ ଶ∑ ݇ఏሺܿ݋ݏߠ௜௝ െ ܿ݋ݏߠ଴ሻଶ௜௜ , (2.1) where the first and second terms on the right hand side show changes in the strain energy arising from deviations in bond lengths and bond angles, respectively, from their equilibrium values, bi and b0 represent the i th bond length and the equilibrium value respectively, θij is the bond angle between bonds i and j to a common atom with the equilibrium value of θ0, and kb and kθ are force constants respectively for the two-body and three-body interactions. The potential parameters were carefully optimized based on first principles density functional theory (DFT) calculations. First, the two-body force constant (kb) was adjusted to fit DFT values for the total energy variation of crystalline Si with varying amounts of strain (from 10 % compressive to 10 % tensile). Then, the three-body force constant (kθ) was adjusted to fit DFT values for the strain energies of five different amorphous Si model structures (of each is within a 64-atom simple cubic cell). Finally, the values for kb and kθ were further refined simultaneously based on several interstitial clusters, such as I4, I7 and I8. Through such careful optimization, we have obtained a set 11 of parameters, kb = 11.976 eV/Å2 and kθ = 2.097 eV, for DFT values of b0 = 2.364 Å and θ0 =109.5°. Using the CRN-MMC approach, we first determined possible fourfold- coordinated structures for interstitial clusters of different sizes (I3 – I10), starting with various initial configurations for each case. The thermal stability of the fourfold- coordinated model structures was carefully checked using tight binding molecular dynamics (TBMD) based on highly optimized semi-empirical potentials developed by Lenosky et al. [67]. The structure and energetics of the stable clusters determined by combined CRN-MMC and TBMD simulations were further refined using first principles calculations based on DFT within the generalized gradient approximation of Perdew and Wang (GGA-PW91) [68]. While the CRN-MMC, TBMD, or DFT method alone is likely limited to sample all possible cluster configurations, the combined approach has demonstrated to be an effective mean to determine the complex minimum energy states of small self-interstitial clusters in Si. All atomic structures and energetics reported herein were calculated using the well established planewave program VASP (Vienna Ab initio Simulation Package) [69]. For the sake of reference, the formation energies of interstitial clusters were also evaluated within the local density approximation (LDA) by relaxing the GGA structures. For the GGA and LDA calculations, the supercell lattice constants were fixed at 5.460 and 5.382 Å, respectively, as obtained from careful volume optimization. A planewave cutoff-energy of 160 eV was used. Vanderbilt ultrasoft pseudopotentials [70] for core- electron interactions were employed, and Brillouin zone sampling was performed using 12 Monkhorst–Pack type k-point meshes. The mesh size was set to (222) for the 216- atom simple cubic supercell, and was properly adjusted with supercell size. For each defect system, all atoms were fully relaxed using the conjugate gradient method until residual forces on constituent atoms become smaller than 510-2 eV/Å. 2.3 Atomic Structure of Compact Self-Interstitial Clusters (In, n ≤ 10) We first calculated the structure and formation energies of the single- and di- interstitials for the sake of reference. As shown in Figure 2.1, the <110>-split interstitial (I) and the C1h di-interstitial (I2) turn out to be most stable in the neutral sate, with formation energies (per interstitial) of 3.80 eV and 2.79 eV, respectively. The results are in good agreement with those from previous DFT calculations [49, 54]. Figure 2.1: Predicted lowest energy configurations for (a) I1, (b) I2, (c) I3, and (d) I4. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated [5]. 13 It is now well accepted that the tri-interstitial (I3) and tetra-interstitial (I4) clusters preferentially form the C2 and D2d structures, respectively, where all atoms are fourfold- coordinated [48,54], as illustrated in Figure 2.1. Our CRN-MMC simulations also predict the ground-state C2 and D2d configurations. It is also worth emphasizing that the CRN- MMC approach consistently yields the lowest-energy structures, irrespective of the initial position of interstitials as long as they are placed sufficiently close to each other. For the penta-interstitial (I5) cluster, the most favorable fourfold-coordinated structure, yet highly distorted while having the trace of I4, was also determined via the CRN-MMC simulations. During TBMD simulation at 1000 K, indeed the highly strained structure was quickly reconfigured to the more stable I4 + I structure where the fifth interstitial is located near the seven-member ring of the I4 cluster (see Figure 2.2). This implies that such fourfold-coordinated geometry is not thermodynamically favorable for the I5 cluster. We have also carefully examined other possible local-minimum structures as suggested by earlier modified embedded atom method calculations [71], such as the structure of I4 + split interstitial. However, they turn out to be at least 0.6 eV less favorable than the I4 + I structure according to our DFT-GGA calculations. Figure 2.3 shows a fourfold-coordinated structure of the hexa-interstitial (I6) in which all six interstitials [in grey (gold)] are placed on a (111) plane. In addition, by introducing two interstitials to the I4 compact structure we obtain another I6 (I4 + 2I) structure with C2 symmetry in which the newly-introduced atoms are placed around the seven-member ring of the I4 cluster. Hereafter the former [(a)] and the latter [(b)] are referred to as I6a and I6b, respectively. Although yielding a highly strained atom 14 (indicated as A), our DFT-GGA total energy calculations show that the FC I6a structure is as favorable as the I6b structure. Figure 2.2: Predicted minimum-energy penta-interstitial (I5) structure with C1 symmetry from three different perspectives, as indicated. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white), and the fifth interstitial is also indicated [5]. 15 Figure 2.3: Predicted two nearly degenerate configurations for the hexa-interstitial cluster (I6) from three different perspectives, as indicated (see upper insets in Figure 2.2). Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). Also, the fifth and sixth interstitials captured by the I4 cluster are also indicated in (b). The corresponding defect symmetry is indicated [5]. For the hepta-interstitial (I7) cluster, we have identified three stable fourfold- coordinated structures. As shown in Figure 2.4, the first one (referred to as I7a) is constituted with a combination of the I3 and I4 ground state structures, lowering the total formation energy by 0.37 eV compared with the I3 and I4 ground state structures apart. 16 The second one (referred to as I7b) consists of two I3 clusters which are connected through an atom with C2 symmetry. It turns out that the center atom leads two four- member rings (associated with the I3 clusters) to two less strained five-member rings, thereby lowering defect-induced strains. As a result, the I7b formation energy per interstitial of 1.94 eV is noticeably smaller than 2.06 eV for the I3 cluster. For the other local-minimum structure (referred to as I7c), all seven interstitials are placed on a (111) plane with CS symmetry, like the I6a structure. According to our DFT-GGA calculations, the I7a is predicted to be 0.39 eV and 0.75 eV more favorable than the I7b and I7c, respectively. These results indicate that the I7 clusters are preferentially formed with the I3 and I4 clusters as building blocks. Our CRN-MMC calculations show that the octa-interstitial (I8) cluster preferentially consists of two stable I4 clusters as building blocks. Depending on the alignment of the I4 clusters, several almost degenerate FC I8 structures have been determined. Among them, as shown in Figure 2.5 the lowest-energy structure exhibits C2 symmetry, where two I4 clusters are placed next to each other while having the C2 rotation axes perpendicular to the (001) plane. The I8 cluster is likely to lower to some extent the induced strain, relative to two isolated I4 clusters. The resulting energy gain is estimated to be about 0.3 eV. 17 Figure 2.4: Predicted three local minimum configurations for the hepta-interstitial cluster (I7) from three different perspectives, as indicated (see upper insets in Figure 2.2). The I7a structure is about 0.39 eV and 0.75 eV more favorable than the I7b and I7c structures, respectively. All structures are fourfold-coordinated. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The corresponding defect symmetries are also indicated [5]. 18 Figure 2.5: Predicted minimum energy configuration with Ci symmetry for the octa- interstitial cluster from three different perspectives, as indicated (see upper insets in Figure 2.2). Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white) [5]. Figure 2.6: Predicted two local minimum configurations for the ennea-interstitial cluster (I9) from three different perspectives, as indicated (see upper insets in Figure 2.2). The I9a structure is about 0.09 eV more favorable than the I9b structure. In (a) the ninth interstitial captured by the I8 cluster is indicated. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The corresponding defect symmetries are also indicated [5]. 19 For the ennea-interstitial (I9) cluster, one can first expect a structure [referred to as I9a, Figure 2.6(a)] where the ninth interstitial is located around an I4 seven-member ring of the I8 cluster, as seen for the I5 cluster earlier. Indeed, a single interstitial can significantly be stabilized near the I8 cluster. Unlike the I5, however we also identify a stable fourfold-coordinated structure which is nearly degenerate with the I9a structure. The highly distorted structure [referred to as I9b, Figure 2.6(b)] looks like a combination of two I4 clusters through an additional interstitial. Finally, we find that the deca-interstitial (I10) cluster preferentially forms a fourfold-coordinated structure. Figure 2.7 shows two stable FC I10 configurations identified, referred to as I10a and I10b, respectively, hereafter. The I10a structure is highly disordered while the I10a shows Ci symmetry, but the former is predicted to be about 0.3 eV favorable than the latter. Note that in the I10b structure all ten interstitials are placed on a (111) plane, like the I6b and I7b structures. One might also expect the I8 + 2I structure where two single interstitials are trapped around the I8 cluster, like the I4 + 2I cluster. However, the I8 + 2I structure turns out to be about 1.7 eV less favorable than the I10a structure. It is worth noting that the energy gain by bond formation exceeds the energy loss by strain arising from the fourfold-coordination as the cluster size increases from I6 to I10. This is apparently due to more flexibility in the bond rearrangement in the larger cluster. 20 Figure 2.7: Predicted two local minimum configurations for the deca-interstitial cluster (I10) from three different perspectives, as indicated (see upper insets in Figure 2.2). The I10a structure is about 0.3 eV more favorable than the I10b structure. Both structures are fourfold-coordinated. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The corresponding defect symmetries are also indicated [5]. 2.4 Stability of Compact Self-Interstitial Clusters (In, n ≤ 10) In Table 2.1, we summarize the formation energies of the small interstitial clusters from our supercell calculations at fixed volumes (as specified) within both the PW91- GGA [see also Figure 2.8(a)] and LDA (local density approximation), together with the values attained using the Keating-like (KT) potential model in which parameters were optimized based on DFT-GGA results. As such, the KT and GGA values are in agreement. In addition, particularly for the fourfold-coordinated defect structures, we can see that the GGA and LDA values are very close. In fact, this is not surprising 21 considering that density functional theory generally works well for such simple covalent interactions. Here the defect formation energy per interstitial [ܧ෠௙ሺ݊ሻ] is given as: ܧ෠௙ሺ݊ሻ ൌ ቄܧሺܰ ൅ ݊ሻ െ ேା௡ே ܧሺܰሻቅ /݊, (2.2) where E(N+n) and E(N) are the total energies of N-atom supercells with a n-interstitial cluster and with no defect. Then the total formation energy [ܧ௙ሺ݊ሻ] is given by: ܧ௙ሺ݊ሻ ൌ ݊ܧ෠௙ሺ݊ሻ. (2.3) The result clearly demonstrates that by and large the formation energy per interstitial decreases as the number of interstitials increases, as also predicted by earlier studies [51,62,63]. From the result, one can also see that the formation energy variation exhibits a non-monotonic trend, with strong minima at I4 and I8. The results consistent with earlier inverse model studies based on experimental observations [62,63], which demonstrated that clusters containing four or eight atoms are particularly stable. The oscillating behavior in the stability of small interstitial clusters has also been predicted by other theoretical studies [72-74]. 22 Table 2.1: Predicted values for the formation energies per interstitial [Êf(n)] of small self-interstitial clusters (In, n = 1 – 10) from DFT-GGA/LDA and Keating-like potential (KT) calculations as indicated. The supercells α, β, and γ consists of 192 + n, 400 + n, and 480 + n atoms, respectively, where n is the number of interstitials. The corresponding defect symmetries are also summarized. 23 Figure 2.8: (a) Predicted formation energies per interstitial (Êf) of small compact interstitial clusters from our DFT-GGA calculations and (b) their approximate capture radii for mobile interstitials. The capture radius is approximated as rc = (3V/4π)1/3, where V is a spherical-like volume that corresponds to the number of Si lattice atoms (excluding interstitials) whose strain energies are greater than a given value (as indicated) [5]. 24 For the identified interstitial clusters, we also estimated their capture radii for mobile interstitials based on the calculation of defect-induced strain fields. For each interstitial cluster, we first counted the number of Si lattice atoms (excluding interstitials) whose strain energies are greater than a given value and then calculated the corresponding volume (V) in bulk Si. This gives a value of capture radius rc = (3V/4π)1/3, assuming a spherical volume. We admit that the approach could be oversimplified, but it should be physically sound and sufficient in approximating changes in the capture radius with cluster size, given that the formation energy of interstitials is a function of local strain [4,210,211]. As summarized in Figure 2.8(b), by and large, the predicted capture radius increases with cluster size, consistent with earlier inverse model studies [62,63]. 2.5 Lowest-energy small compact clusters (In, 5 ≤ n ≤ 12) Figure 2.9 shows minimum energy structures which we have identified for I5 through I12 using a combination of CRN-MMC, TBMD, and DFT-GGA calculations. For each cluster size, we first constructed possible fourfold-coordinated structures using CRN-MMC simulations, followed by TBMD simulations at high temperatures (> 1000 K) to check their thermal stability. Then, we used DFT-GGA calculations to refine the geometries of the identified stable clusters and compared their formation energies to determine the lowest energy configuration among them. The combined approach has proven successful in determination of minimum energy configurations for Si self- interstitial clusters (In, n  3), particularly when they prefer fourfold-coordination [5]. 25 Here, we only present the lowest energy state for each cluster size among several local minima identified from our extensive search. Other stable structures can be found elsewhere [5]. Figure 2.9: Predicted minimum-energy configurations for (a) penta-, (b) hexa-, (c) hepta-, (d) octa-, (e) ennea-, (f) deca-, (g) hendeca-, (h) dodeca-interstitial defects in Si. For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated. Grey (gold) balls represent more distorted atoms than the rest of the lattice atoms (in white) [6]. 26 Our results show that small interstitial clusters tend to favor compact structures, but the compact geometry is no longer favorable when the cluster size exceeds 10 interstitials. All atoms in the hendeca-interstitial (I11) [Figure 2.9(g)] and dodeca- interstitial (I12) [Figure 2.9(h)] clusters are fourfold-coordinated and they are elongated along the [110] direction. While the I11 cluster exhibits a somewhat asymmetric shape, the I12 structure is perfectly symmetric and less strained. The I12 structure has mirror symmetry with respect to the (110) plane (perpendicular to the C2 rotation axis along the [110] direction). Having four interstitials in each unit, the core part [see also Figure 2.14(a)] is characterized by four adjacent six-membered rings surrounded by five-, six-, and seven-membered rings. The (110) surface atoms are also fourfold-coordinated and their bond angles and lengths range from 98.8o to 128.2o and from 2.243 Å to 2.348 Å, respectively. Both parameters deviate insignificantly from their respective equilibrium values of 109.5o and 2.365 Å in crystalline Si. The I12 formation energy is predicted to be 1.60 eV per interstitial, which is substantially lower than 3.80 eV for the split-<110> mono-interstitial. We also find that the octa-interstitial (I8) cluster (comprising two stable compact I4 clusters) is very stable. For the minimum energy I8 configuration, two complete I4 core structures are adjacent along the [110] direction, as shown in Figure 2.9(d). The I4-I4 alignment with C2v symmetry yields two eight-membered rings in the interface region, which turns out to lower the induced strain relative to two isolated I4 core clusters. While there is an energy variation as the I4-I4 alignment within the Si lattice changes [5], the 27 predicted ground state structure of I8 is about 0.71 eV more favorable than the case where two I4 clusters are fully separated. As a consequence of the high thermal stability of the I4 and I8 structures, the penta-interstitial (I5) [Figure 2.9(a)] and ennea-interstitial (I9) clusters [Figure 2.9(e)] favor the I4 + I and I8 + I configurations, respectively, where the additional single interstitial is located near their host I4 or I8 cluster. For the hexa-interstitial (I6) [Figure 2.9(b)] and hepta-interstitial (I7) clusters [Figure 2.9(c)], the combined I4 + I2 and I4 + I3 configurations appear energetically favored with energy gains of 1.16 eV and 0.37 eV, respectively, over their fully-separated counterparts (isolated I2, I3, or I4 clusters). Similarly, a stable I8 + I2 configuration is identified for the deca-interstitial (I10) cluster, but turns out to be 0.36 eV less favorable than the stable fourfold-coordinated structure as shown in Figure 2.9(f). Our calculations also predict the I8 + I3 [Figure 2.10(a)] and I8 + I4 [Figure 2.10(b)] structures to be 0.93 eV and 1.59 eV less stable than the elongated, FC I11 and I12 structures, respectively. 28 Figure 2.10: Compact, minimum energy configurations with corresponding defect symmetries for (a) hendeca- and (b) dodeca-interstitial defects in Si. For each defect, the top and bottom panels show two different views, as indicated. Light grey (gold) wireframe represents the bulk Si lattice. Dark grey spheres represent interstitial atoms and their highly strained neighbors [7]. 29 2.6 Compact-to-elongated Transition Figure 2.11 summarizes both the calculated formation energies of small interstitial clusters and several configurations of elongated structures. The predicted Ef values of 3.80 eV, 2.79 eV, 2.06 eV, and 1.85 eV, respectively, for the I, I2, I3 and I4 clusters are in good agreement with previous DFT-GGA calculations [49, 54]. The predicted formation energies exhibit an oscillating trend as the cluster size varies, with strong minima at n = 4 and n = 8, consistent with inverse modeling experiments [62, 63]. Based on these calculations, we discuss the nucleation and growth of large, extended defects by evaluating three proposed models. Prior to our study, the presence of stable, chain-like, elongated, interstitial clusters [Figure 2.12] for n ≥ 3 was reported by Kim et al. [51]. These chain-like configurations are extended by concatenating dumbbell-like split-<110> interstitials to the existing core structure. While the fourfold-coordinated compact cluster was identified to be more stable than the chain-like, elongated cluster at n = 4, this study did not report the existence of larger compact clusters based on the I4 core that have lower energy configurations than the chain-like, elongated clusters. At the time, the existing calculations led to the conclusion that small compact clusters may evolve into chain-like, elongated clusters for n ≈ 5-8. This model of the compact-to-extended structural growth transition mechanism proposed the chain-like, elongated tri-interstitial cluster as a key component and seed for extended structure growth [51, 54]. 30 Figure 2.11: Calculated formation energies per interstitial (Êf) for interstitial clusters shown in Figure 2.1, 2.9, and F#_013 as well as chain-like elongated structures (indicated as ‘Chain-like’) as a function of cluster size (n). For the chain-like case, the atomic structures of Kim et al. [51] were recalculated within DFT-GGA. To minimize possible interactions between a defect and its periodic images, we evaluated the formation energies by scaling the supercell size as follows (‘This work’): 192+n (I1–I2), 400+n (I3- I6), 480+n (I7-I10), 560+n (I11-I13), and 672+n (I14-I16) atom supercells, where n is the number of interstitials. Likewise, the following supercells sizes were used (‘Chain-like’): 400+n (I3-I4), 480+n (I5-I6), 560+n (I7-I8), and 640+n (I9-I10) atom supercells. For n > 10, the chain-like cluster formation energy is forecasted using Êf(n) = {6.83 + (n − 2)×1.28}/n, see also Figure 2.17 [7]. 31 Figure 2.12: Configuration for the chain-like, elongated hexa-interstitial defect with C1 symmetry in Si. The top and bottom panels show two different views, as indicated. In both views the strain energy distribution around the cluster is binned into three different levels [high (≥ 0.25eV) = dark grey (red) balls, medium (0.15 – 0.25 eV) = grey balls, and low (< 0.15 eV) = light grey (gold) wireframe, based on strain energy values from force field calculations] [7]. 32 Our recent studies using the combined computational approach have shown that the compact configurations are comparable or more stable than the chain-like, elongated clusters for n  8, as shown in Figure 2.12. Our results also suggest that the transition from compact to well-ordered chain-like configurations would be difficult because the compact I4 core structure is relatively very stable. Our TBMD simulation at 1100K shows that the I4 + I structure remains nearly unchanged for 20 ps, whereas the compact I3 ground state cluster easily collapses when it captures an additional interstitial and would further reconfigure itself into the stable I4 core structure. This behavior suggests that the stable I4 and I8 compact clusters would kinetically and/or thermodynamically deter the formation of small chain-like clusters (n ≤ 10). For n > 10, the elongated configuration becomes more stable than the compact shape. The elongated configuration generally has strained {110} terminating surfaces or interfaces. This excess strain energy deters chain-like, elongated configurations for small n, but the terminating surface effects are diminished as n increases. This observation suggests that fourfold-coordinated compact clusters may evolve into chain-like, elongated structures for n  10 and further growth may occur by the capture of interstitials at the strained {110} terminating surfaces. Another plausible transition growth mechanism may involve the elongated, FC I11 and I12 clusters identified in our work. In the size regime where n > 10, we do not anticipate interconversion between the elongated, fourfold-coordinated clusters we have presented and the chain-like, elongated configurations proposed by Kim et al. because of 33 the well-ordered structure and high thermal stability exhibited by the elongated fourfold- coordinated clusters. Figure 2.13 (upper right panel) shows an atomic-level strain distribution profile for the FC I12 cluster based on average bond lengths. The color of each atom is gradient- shaded according to the average bond length shared with its four neighbors. Compressive strain is indicated by red, strain-free is white, and tensile strain appears as blue. Larger strain is seen near the (110) interfaces, indicating that they are more active than the {311} and {233} interfaces. Our TBMD simulations indeed demonstrate that an additional interstitial placed around the I12 cluster is preferentially captured at the (110) interfaces, as defined, to form the I13 cluster. Additionally, Figure 2.13 reveals the minimum energy fourfold-coordinated configurations of the I14, I15, and I16 clusters during the proposed growth sequence from the I12 configuration. The I13, I14, I15 and I16 formation energies per interstitial are predicted to be 1.63 eV, 1.63 eV, 1.60 eV, and 1.50 eV, respectively. Note that the I16 structure with C2h symmetry is an extension of I12 with an additional core unit added in the [110] direction. By capturing additional interstitials, the interstitial defect structure grows preferentially along the [110] direction. This behavior suggests that I12 can serve as an effective nucleation seed for the growth of larger extended defects. 34 Figure 2.13: Minimum energy structures for I12 through I16 with the I12-like core, which exhibit preferred growth along the [110] direction when n > 10. Configurations for I12 through I16 are depicted with dark grey balls showing highly distorted atoms and light grey (gold) wireframe representing the bulk lattice. Additionally, the upper right panel provides an atomic-level strain distribution profile for the I12 cluster and nearby bulk Si based on average bond length. Each atom is assigned an average bond length based on the four bonds formed with its nearest neighbors and is gradient-shaded accordingly. The color spectrum shifts from red to white to blue (RWB) as the average bond lengths shift from compressive to strain-free to tensile. The RWB spectrum range displayed here covers the Si-Si DFT equilibrium bond length of 2.36  0.07Å [7]. 35 2.7 Growth to {311} Extended Defects We also examined the possible evolution of small, fourfold-coordinated and chain-like, elongated clusters into extended defects with a {311} habit plane. In Figure 2.14, we compare the I12-like core structure [Figure 2.14(a)], the chain-like, elongated core structure [Figure 2.14(b)] of Kim et al., and a plausible {311} core structure [Figure 2.14(c)]. In Figure 2.14, all three defect cores shown are cross-sections that effectively extend infinitely in the [110] direction using periodic boundary conditions of the supercell method. The {311} core consists of four six-membered rings in the middle layer and six five- and six seven-membered rings in the outer layers, as proposed by earlier studies [75-77]. Other {311} core configurations have also been predicted [51,75- 78], but the most favorable structural sequence is still undetermined. We feel our conclusions are general enough that our observed results of the interaction between small interstitial clusters and the exemplary {311} structure in Figure 2.14(c) can be applied to other plausible {311} defect structures. For each repeating unit of the structures shown in Figure 2.14, the I12-like and {311} cores contain four interstitials, while the chain-like, elongated structure has only two. The I12-like core has a higher interstitial density than either the chain-like or {311} cores. Compared to the I12-like core, the chain-like core must be twice as long along [110] and the {311} core must be twice as long along [33ത2] to match the number of net embedded interstitials. The I12-like and {311} cores both exhibit a representative fraction of the {311}-habit plane, but the chain-like, elongated structure does not. In addition, 36 transformation from the I12-like to {311} core is plausible via structural relaxation along the <233> direction. However, it appears unlikely that the chain-like, elongated core of Figure 2.14(b) could reconfigure into the {311} core of Figure 2.14(c) because the chain- like configuration contains two less net interstitials per unit length along the <110> direction. The required structural transformation to convert the chain-like core of Figure 2.14(b) into the {311} core of Figure 2.14(c) would need the chain-like configuration to compress by 50% along <110> while simultaneously stretching by 100% along <233>. This type of physical transformation on an extended scale is improbable. The less dense {311} core has a lower Ef than the I12-like core, which suggests the {311} core is more relaxed than the I12-like core. The predicted formation energies of the {311}, I12-like, and chain-like, elongated cores are 1.11 eV, 1.29 eV, and 1.28 eV, respectively, assuming the cores are infinitely long in the [110] direction, as depicted in Figure 2.14. Our results reinforce the possibility of a compact-to-elongated transition mechanism initiated from fourfold-coordinated, elongated I11 and I12 clusters. This conclusion is consistent with earlier experimental results which suggested the incipience of a compact-to-elongated transition for n  10 [62,63]. 37 Figure 2.14: Three different defect core structures which are infinitely long in the [110] direction: (a) I12-like (as shown for the FC I12 cluster identified in this work); (b) chain- like; and (c) {311}. For each core, the left and right panels show two different views, as indicated. Grey balls indicate more distorted atoms than the rest of the light grey (gold) wireframe lattice [7]. 38 Figure 2.15: Minimum energy configurations for {311} core structures of (a) I12 and (b) I16 clusters in Si. For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated. Dark grey (red) balls indicate highly distorted atoms in the (110) interfaces, grey balls represent moderately distorted atoms, and light grey (gold) wireframe represents the bulk Si lattice [7]. 39 Figure 2.16: Minimum energy configurations for the (110) edges of (a) the elongated I12- like, (b) chain-like, and (c) {311} core configurations. Dark grey (red) balls indicate highly distorted atoms contained in the (110) interfaces, grey balls represent moderately distorted atoms, and light grey (gold) wireframe represents the bulk Si lattice [7]. 40 We also compared the relative stability of both the I12 [Figure 2.15(a)] and I16 [Figure 2.15(b)] clusters with the {311} core configuration, since the {311} core is energetically more favorable than the I12-like core. As shown in Figure 2.16(c), the (110) interfaces can be fourfold-coordinated for the {311} core structure, even though each surface initially contains two dangling bonds. The fourfold-coordination lowers the Ef of I12 by 1.1 eV over the initial structure with dangling bonds. Using CRN-MMC simulations, we also confirmed that the fourfold-coordinated interfaces can favorably evolve for larger clusters with a {311} core configuration; however, we find that the (110) interface configuration is more strained for the {311} core than for the I12-like core [Figure 2.16(a)]. Our calculations predict the {311} core structures (with Ci symmetry) to be 1.4 eV and 0.8 eV less favorable than the I12-like core containing I12 and I16 clusters, respectively. Figure 2.17 also compares the relative stability of large clusters with the I12-like, chain-like, and {311} cores which we approximate using: Ê௙ሺ݊ሻ ൌ ா೑ ೐೙೏ାሺ௡ି௡ᇲሻÊ೑೎ ௡ , (2.4) where ܧ௙௘௡ௗ is the {110} interface energy, Ê௙௖ is the core energy per interstitial for the semi-infinite elongated structure, and n′ represents the number of (110) interface terminating atoms. For the I12-like and {311} cores, we use the core formation energies of Ê௙௖ = 1.29 eV and 1.11 eV to estimate the corresponding (110) interface energies as 41 ܧ௙௘௡ௗ = 8.74 eV and 11.62 eV, respectively. Note that both cases have four (110) interface atoms (n′ = 4). Similarly, for the chain-like, elongated core, we obtain Ê௙௖ = 1.28 eV and ܧ௙௘௡ௗ = 6.83 eV, where n′ = 2 (the chain-like structure has two (110) interface atoms). From this result, the {311} core structure becomes most favorable for n  20 atoms, while the I12-like core should prevail for n < 20 atoms. This supports a possible structural relaxation from the I12-like to {311} core structure near the n ≈ 20 size regime as the interstitial cluster grows along the [110] direction. Experimentally observed {311} defects are typically as long as 1 m along the <110> direction and range from 1 to 100 nm in width along the <233> direction. Note that the <233> direction is perpendicular to the <110> and <311> directions as shown in Figure 2.18. This implies that a {311} defect gradually expands along the <233> direction while it quickly elongates along the <110> direction as it incorporates more interstitials during a thermal annealing process. This is shown schematically in Figure 2.18. We also constructed an extended defect with a {311} habit plane by adding repeating units to the {311} core of Figure 2.14(c) along the <233> direction. The predicted Ef per interstitial is significantly lower than the {311} core structure in Figure 2.14(c) with Ef () = 0.89 eV, which is in good agreement with expectations of 0.7 to 1.4 eV as reported for extended {311} defects in previous studies [209]. 42 Figure 2.17: Predicted formation energies per interstitial of I12-like, chain-like, and {311} core configuration clusters as a function of size (n) using Ef(n) = {ܧ௙௘௡ௗ + (n – n′)Ê௙௖}/n, where ܧ௙௘௡ௗ is the (110) interface energy, Ê௙௖ is the core energy per interstitial for the semi-infinite structure, and n′ is the number of interface atoms. Here, the interface energies were obtained by numerically fitting DFT-GGA values for the formation energies of small clusters: I12 and I16 were used for the I12-like and {311} cases, while I8, I9 and I10 represented the chain-like, elongated structures. The predicted values are based on Ê௙௖ = 1.29, 1.28, and 1.11 eV and ܧ௙௘௡ௗ = 8.74, 6.83, and 11.62 eV for the I12-like, chain-like, and {311} cases, respectively. For the number of interface atoms, n′ = 4 was used for the I12-like and {311} structures, while n′ = 2 was used for the chain-like, elongated structure [7]. 43 Figure 2.18: Block diagram for a typical {311} extended defect structure. Characteristic dimensions observed experimentally and alignment of the defect within the bulk Si lattice is shown. From the computational results obtained, incremental growth along <233> (small block extension) followed by more favorable growth along <110> (block arrows) is anticipated [7]. Finally, we briefly discuss the predicted growth behavior along the <233> direction based on our preliminary results. Figure 2.19(a) exhibits an atomic-level strain distribution profile throughout the model {311} core structure using the same parameters as Figure 2.13 (upper right panel). According to our analysis, the (33ത2) interface is under tensile strain (blue), while the (11ത3ത) interface is under mild compressive strain (pink), which is consistent with a previous theoretical study [77]. Strain distributions from all extended defects with a {311} habit plane we studied show a similar trend. This indicates that additional interstitials, which generally induce a localized compressive strain field, are more likely to aggregate at the (33ത2) tensile interface to quench strain rather than at the (11ത3ത) compressive interface. 44 Figure 2.19: (a) An atomic-level strain distribution profile throughout the {311} core structure and nearby bulk material based on average bond length. Each atom is assigned an average bond length based on the four bonds formed with its nearest neighbors and is gradient-shaded accordingly. The color spectrum shifts from red to white to blue (RWB) as the average bond lengths shift from compressive to strain-free to tensile. The RWB spectrum range displayed here covers the Si-Si DFT equilibrium bond length of 2.36  0.07Å. Additionally, the constituent rings of the model {311} core have been shaded grey to help visualize the defect cross-section in the bulk. The “I” (interstitial) and “E” (edge) ring annotations reflect the nomenclature of Kohyama and Takeda [82] to generalize the periodicity seen in {311} planar structures. From left to right, this structure would be labeled “EIIE”. The “I” repeating unit contains two net interstitials per period along <110>. From the strain profile, the structure core is compressively strained, the {311} interface planes are under slight compression, and the {233} interface planes are under tensile strain. Minimum energy configurations for (b) “{311} core + I3” and (c) “{311} core + I4” structures. In both cases, the interstitials appended themselves to the tensile {233} interface to help minimize excess strain energy. The appended tri-interstitial and tetra-interstitial configurations (bright red) show strong resemblance to the ground state I3 and I4 compact configurations of Figure 2.1 [7]. 45 To verify this potential interstitial growth mechanism influenced by the local strain field, we added additional interstitials near an extended defect with a {311} core. Our CRN-MMC simulations clearly show that the interstitials preferentially append to the extended defect at the (33ത2) interfaces, as shown in Figure 2.19. In Figure 2.19(b), three additional interstitials form a stable tri-interstitial cluster whose configuration is similar to the ground state I3 cluster. The appended tri-interstitial cluster is compressively strained, as indicated by the bright red color. At the (33ത2) interface, two five-membered rings of the ground state I3 cluster combine with two seven-membered rings of the {311} core, resulting in two six-membered rings. Likewise, in Figure 2.19(c), four additional interstitials form a stable tetra-interstitial cluster which resembles the ground state I4 cluster. Again, the appended tetra-interstitial cluster is compressively strain and appears bright red. At the (33ത2) interface, two seven-membered rings of the ground state I4 cluster combine with two seven-membered rings of the {311} core, resulting in two seven-membered rings. As a result, the hybrid “{311} core + I3” and “{311} core + I4” structures have binding energies of 1.1 eV and 1.2 eV with respect to the separated “{311} core + I3” and “{311} core + I4” component structures, respectively, according to our force field calculations. For comparison, we also identified the structure and stability of possible I3 and I4 clusters at the (11ത3ത) interface plane. Unlike the clusters appended to the (33ത2) interfaces, the I3 and I4 clusters do not yield similar structures to the ground state I3 and I4 clusters, respectively, when formed in the mildly compressive region along the (11ത3ത) interface. As a consequence, the binding energies of the I3 and I4 46 clusters at the (11ത3ത) interface were significantly lower than those at the (33ത2) interface by 0.6 and 0.5 eV, respectively. We also evaluated the addition of the small interstitial clusters near the (1ത1ത0) interface of the {311} defect structure. We estimated the binding energy of the tetra- interstitial cluster (when it forms another repeating unit along the [1ത1ത0] direction) at the (1ത1ത0) interface as 2.5 eV, which is twice as large as the binding energy found at the ( 233 ) interface. Given the binding energy results of our work and the typical ratio of dimensions observed experimentally for extended {311} defect structures, it is reasonable to propose that the defect growth rate along <110> is largest, followed by the growth rate along <233>, while the growth rate along <311> should be slowest. As depicted by the block diagram in Figure 2.18, a possible multi-step mechanism is shown that integrates our computational results for extended defect growth. First, newly formed small compact cluster units append to the {233} interfaces. Second, these new additional interstitial building blocks seed rapid chain growth along the <110> direction. The eventual result is the extension of the {311} defect in the <233> dimension. We also propose the possibility of relatively highly-strained small compact clusters playing a critical role as initiators for further structural relaxation of the already grown {311} defect, ultimately leading to the formation of well-relaxed, extended defects with {311} habit planes. 47 2.8 Summary In this work, we present an effective computational approach to determine the structure and energetics of compact self-interstitial clusters in Si, which combines continuous random network (CRN) model, tight binding molecular dynamics (TBMD), and density functional theory (DFT) calculations. For a given defect cluster size, we first constructed possible fourfold-coordinated structures using CRN based Metropolis Monte Carlo (CRN-MMC) simulations, followed by TBMD simulations at high temperatures to check the stability of the fourfold-coordinated structures. Then, DFT calculations were performed to refine the structures of stable defect clusters identified, and compare their formation energies to identify local minimum-energy structures. This combined approach has been demonstrated to be an efficient technique for determining the lowest- energy configurations of defect clusters (In, n  3), particularly when they prefer to form fourfold-coordinated structures. Using the combined simulations, we have identified new stable compact structures for self-interstitial clusters with n = 6 through 10. The calculated relative formation energies per interstitial of the small compact defect clusters (n  10) are also compared with those from earlier inverse model studies based on experiments, which shows good agreement. Both our calculation and the inverse modeling show that the formation energy variation exhibits a non-monotonic trend, with strong minima at I4 and I8. The improved understanding greatly assist in explaining and predicting the temporal and spatial changes of interstitial concentrations, which is essential for describing interstitial mediated dopant diffusion and clustering in 48 ultrashallow junction formation required for future generations of Si-based electronic devices. We have determined stable, compact geometries for small interstitial clusters (n < 10) and demonstrated that stable I4 and I8 compact clusters would kinetically or/and thermodynamically inhibit the formation of chain-like, elongated clusters. When the cluster size exceeds 10 atoms, our calculations show that elongated structures become more favorable energetically, although they commonly consist of highly-strained {110} interfaces which make them less favorable configurations for smaller clusters. In particular, the newly discovered FC I12 configuration is found to serve as an effective nucleation center for larger extended defects. For the first time, this theoretical work provides explicit support for earlier experiments which suggested the occurrence of the compact-to-elongated transition at n  10. In addition, the predicted formation energies per interstitial generally decrease with cluster size, but also exhibit an oscillating trend with strong minima at n = 4 and n = 8, which is consistent with earlier inverse model studies based on experiment. While the FC I12-like core is energetically less favorable than a typical {311} defect core, our theoretical results suggest the possible occurrence of further structural relaxation from the I12-like to {311} core configuration as the cluster grows larger than n ≈ 20 along the <110> direction. We also propose the possible formation of compact interstitial clusters at the {233} interfaces of extended {311} defect structures driven by minimization of local strain energy. These compact clusters formed at the {233} interfaces can potentially seed further growth along the <110> direction, ultimately leading to incremental growth along <233> and more rapid growth along 49 <110> for extended {311} defect structures. Our integrated modeling technique has many applications ranging from identification of small interstitial cluster configurations to prediction of the growth behavior of extended defect structures. 50 Chapter 3 Formation and Structure of Vacancy Defects in Silicon 3.1 Introduction Vacancies in crystalline Si have sustained interest because these native defects can greatly influence the mechanical, electrical, and optical properties of the bulk semiconductor. Vacancies are often introduced during manufacturing processes such as ion implantation and crystal growth and expected to predominantly reside in clusters or complexes with other impurities because mono-vacancies are highly mobile, even at ambient temperatures [79]. Vacancy defects can getter impurities and also recombine with self-interstitials to reduce defect concentrations in active electronic device regions [80]. While the formation of stable voids on the scale of a few nanometers was evidenced by transmission electron microscopy (TEM) [81], direct experimental characterization of small vacancy clusters is challenging because their dimensions often exceed the capabilities of TEM. Small vacancy defects have been of particular interest because they are a main source or a getter for mobile vacancies and interstitials which are largely responsible for dopant transient enhanced diffusion and electrical deactivation in ultrashallow junction formation for Si-based electronic devices [82-87]. Previous theoretical studies [88-91] proposed “part of hexagonal ring” (PHR) and spherically shaped clusters (SPC) models for the structure of vacancy clusters. According to the models, the V6, V10 and V14 PHR 51 configurations are predicted to be particularly stable, because of their relatively reduced number of Si-Si broken bonds. The PHR and SPC models have also been used to explain large open volume defects [92-97]. However, fourfold-coordinated (FC) V3, V4 and V5 clusters were recently identified to be more energetically favorable than their PHR counterparts [98]. These fourfold-coordinated structures were obtained by placing additional Si atoms to terminate dangling bonds in ring hexa-vacancies. For larger vacancy defects (Vn, n  7), no explicit theoretical account is currently available on their stable fourfold-coordinated configurations which could be too complex to be determined by simple trial and error static calculations. In particular, few studies have examined the likelihood of high-symmetry configurations for stable FC vacancy clusters, as seen in self-interstitial clusters which exhibit well-ordered extended configurations [5-7, 10- 12,51,75,76]. In addition, limited information exists concerning the relative stability between competing FC and PHR/SPC-type structures with cluster size. Finally, little is known regarding the kinetic role in formation of both FC and open-volume (PHR/SPC- type) defects. In this work, we establish the structure and formation energies of fourfold- coordinated vacancy defects in the size range of 3-48 vacancies using a combination of continuous random network model based Metropolis Monte Carlo (CRN-MMC), tight- binding molecular dynamics (TBMD) and density functional theory (DFT) calculations. The combined approach has been demonstrated to be an effective means to identify stable fourfold-coordinated defect clusters in Si [5,6]. We also examine the relative stability between small fourfold-coordinated and PHR vacancy defects (Vn, 3  n  48), showing 52 that for each cluster size the newly identified fourfold-coordinated structure is energetically more favorable than the conventional PHR structure. In addition, our results highlight the identification of an extremely stable FC V32 configuration with high symmetry, which likely influences the growth of larger vacancy structures. We also discuss kinetic feasibility of the formation of thermodynamically-favored FC vacancy clusters using TBMD simulation results. The fundamental findings will greatly contribute to a better understanding of the properties of native defects in Si, and their impact on relevant materials properties and processing. 3.2 Computational Approach All atomic structures and energies reported herein were calculated using a plane- wave basis set pseudopotential method within the generalized gradient approximation of Perdew and Wang (GGA-PW91) [68] to density functional theory (DFT), as implemented in the Vienna Ab initio Simulation Package (VASP) [69]. Vanderbilt-type ultrasoft pseudopotentials [70] were used for core-electron interactions. Outer electron wave-functions were expanded using a planewave basis set with a kinetic energy cutoff of 160 eV. For Brillouin zone sampling, we used a (2×2×2) Monkhorst-Pack k-point mesh for 480–n and 576–n atom supercells and only one Gamma k-point sampling for 1000−n atom supercells. For each defect system, all atoms were fully relaxed using the conjugate gradient method until residual forces on constituent atoms became smaller than 5×10-2 eV/Å. The supercell DFT calculations employed a fixed Si lattice constant of 5.46 53 Å as obtained from volume optimization of defect-free Si. To take into account possible incomplete relaxation of the defect system due to limited supercell size, adjustment of the FC defect formation energies was also made using force field calculations with increasing supercell size. This adjustment appears to be necessary for improved estimation of the formation energies of large FC clusters as discussed later. For TBMD simulations, semi- empirical potentials developed by Lenosky et al. [67] were used. A Keating (KT)-like valence bond model was employed for CRN-MMC calculations. Within the Keating-like valence force model, the strain energy (Estrain) is given as ܧ௦௧௥௔௜௡ ൌ ଵଶ∑ ݇௕ሺܾ௜ െ ܾ଴ሻଶ ൅ ଵ ଶ∑ ݇ఏሺܿ݋ݏߠ௜௝ െ ܿ݋ݏߠ଴ሻଶ௜௜ , (3.1) where bi is the i th bond length and θij is the bond angle between bonds i and j, and the equilibrium and force constants are b0 = 2.365 Å, θ0 = 109.5°, kb = 6.951 eV/Å2 and kθ = 1.868 eV. A detailed description of KT parameter optimization can be found in Sec. 2.2. 3.3 Determination of Fourfold-coordinated Vacancy Clusters We first determine the structure and energetics of fourfold-coordinated vacancy clusters in the size range of 3-18 vacancies. For each cluster size, we first generate possible fourfold-coordinated configurations using CRN-MMC simulations, followed by TBMD simulations at high temperatures (> 1000 K) to check their thermal stability. Then, we employ DFT-GGA calculations to refine the geometries of the stable clusters, and compared their formation energies to determine the lowest-energy structure among them. The combined approach has been proven to be an effective way to identify stable 54 fourfold-coordinated native defects in crystalline Si [5,6]. Figure 3.1, 3.2, 3.3, and 3.4 show the predicted minimum energy configurations of fourfold-coordinated vacancy defects (Vn, 3  n  18). Compared with previously reported small fourfold-coordinated vacancy clusters (Vn, 3 ≤ n ≤ 6) [98], our calculations yield the same configurations for V3, V4 and V6. The previously proposed V5 structure was obtained by relaxing its conventional PHR state [98], exhibiting a combination of three five- and two six-membered rings. On the other hand, as shown in Figure 3.1(c), the newly identified V5 structure consists of one five-, two six-, and two seven-membered rings, as typically seen for well relaxed {311} extended defects [75,76]. The new structure is predicted to be 0.33 eV more favorable than the previous one. It is evident that larger-membered rings are more flexible than smaller ones, which in turn more effectively relieves the defect-induced strain. In addition, we find that the fourfold-coordinated vacancy clusters commonly contain four- or/and five-membered rings which are preferentially placed in the center region, while surrounded by more flexible larger membered (seven- or/and eight-membered) rings. For instance, as shown in Figure 3.2(b) the predicted FC V8 structure with C2h symmetry is composed of two symmetric voids which face each other with four-membered ring as mirror plane, and each void is surrounded by one four-, six six-, and two seven- membered rings. The fourfold-coordinated vacancy clusters are tensilely strained with bond length deviations of -0.05 – 0.4 Å from the equilibrium value of ≈ 2.36 Å, see Figure 3.5. Analyzing calculated maximally localized Wannier functions for the defect 55 configurations, we see that Wannier centers insignificantly, mostly less than 0.1 Å, depart from the mid positions of two bonded atoms (see Figure 3.9(a) for the V12 case). The spread of the Wannier functions is also close to that in bulk crystalline Si, confirming fourfold-coordination of the vacancy defects. In addition, our density of state analysis shows no energy levels within the Si band gap. This indicates that the fourfold- coordinated vacancy defects are optically inactive, which may in turn impede their direct characterization using optical and electrical measurements. 56 Figure 3.1: Predicted fourfold-coordinated configurations for (a) V3, (b) V4, (c) V5, and (d) V6 in Si. Grey (gold) balls indicate more distorted atoms than the rest of the lattice atoms (in white). For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated. 57 Figure 3.2: Predicted fourfold-coordinated configurations for (a) V7, (b) V8, (c) V9, and (d) V10 in Si. Grey (gold) balls indicate more distorted atoms than the rest of the lattice atoms (in white). For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated. 58 Figure 3.3: Predicted fourfold-coordinated configurations for (a) V11, (b) V12, (c) V13, and (d) V14 in Si. Grey (gold) balls indicate more distorted atoms than the rest of the lattice atoms (in white). For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated. 59 Figure 3.4: Predicted fourfold-coordinated configurations for (a) V15, (b) V16, (c) V17, and (d) V18 in Si. Grey (gold) balls indicate more distorted atoms than the rest of the lattice atoms (in white). For each defect, the left and right panels show two different views, as indicated. The symmetry of each defect is also indicated. 60 Figure 3.5: Distributions of bond length (left panels) and bond angle (right panels) deviations associated with the distorted atoms as shown in Figure 3.1, 3.2, 3.3, and 3.4. The dotted lines indicate the equilibrium values of r0 = 2.365 Å and 0 = 109.5 ° [8]. 61 3.4 Relative Stability between Fourfold-coordinated and PHR Vacancy Defects (Vn, n ≤ 48) Figure 3.6(a) shows the calculated formation energies of small vacancy defects (Vn, n = 3 – 18) in both fourfold-coordinated and PHR configurations. Here, the formation energy per vacancy [Ê௙ሺ݊ሻ] is given by: Ê௙ሺ݊ሻ ൌ ቄܧሺܰ െ ݊ሻ െ ேି௡ே ܧሺܰሻቅ /݊, (3.2) where E(N–n) and E(N) are the total energies of N-atom supercells with a n-vacancy cluster and with no defect, respectively. Then the total formation energy [ܧ௙ሺ݊ሻ] is given by: ܧ௙ሺ݊ሻ ൌ ݊Ê௙ሺ݊ሻ. (3.3) In the inset, the total energy differences between the fourfold-coordinated and the PHR states are also presented. For the PHR vacancy clusters, care is taken to ensure possible pairings between dangling bonds created by removal of Si lattice atoms, via structural relaxation using combined CRN-MMC and DFT calculations. The rigorous ionic relaxation leads to reconfiguration of V5 to a fourfold-coordinated state as also reported by Makhov and Lewis [98]. Given the ease of reconfiguration, here we consider the fourfold-coordinated structure as the PHR V5, rather than the conventional PHR structure. Recall that the fourfold-coordinated PHR state is 0.33 eV less favorable than our newly identified fourfold-coordinated structure. For other PHR clusters, there is no significant 62 structural rearrangement. Hence, the ionic relaxation leads to insignificant variations in their formation energies due to the well known flatness of the total energy surface of vacancy defects. Figure 3.6: (a) Calculated formation energies per vacancy (Êf) and (b) binding energies of vacancy clusters as a function of cluster size (n) for both fourfold-coordinated (indicated as “Fourfold”) and PHR (“PHR”) configurations. For the PHR case, the atomic structures from previous studies [90,91] were recalculated within DFT-GGA. To minimize possible interactions between a defect and its periodic images, we carefully evaluated the formation energies by changing the supercell size; 480–n and 576–n atom supercells, where n is the number of vacancies, were used for V1-V14 and V15-V18, respectively. The inset shows a variation in the total energy difference (∆E in eV) between the fourfold-coordinated and PHR cases [8]. 63 For the small vacancy defects, our results demonstrate that the fourfold- coordinated configurations are energetically more favorable than the PHR configurations. The energy difference can be as large as 4 eV when n = 12. While the information related to the fourfold-coordinated states was lacking, the PHR states have been mostly considered for understanding the behavior and properties of vacancies. In addition, spherically shaped clusters (SPC) have also been identified to be stable. Particularly for V17, the SPC was expected to be more favorable than the PHR cluster [91]. Our DFT- GGA calculation, however, predicts the FC V17 cluster to be 2.6 eV more favorable than the SPC. According to our calculations, small vacancy defects thermodynamically favor fourfold-coordination in crystalline Si. This is apparently due to the fact that the energy gain by bond formation exceeds the strain energy associated with bond stretching and bond angle distortion. It is also well established that energetically small self-interstitial clusters (In, n  3) prefer to be fourfold-coordinated [5,6]. Considering less stiffness in the tensily strained structure compared to the compressively strained one, one can expect that the fourfold-coordination of vacancy defects (which are commonly under tensile strain) will be energetically more facile than interstitial defects (which are typically under compression). Indeed, for the Keating-like valence force model, the two-body force parameters of kb = 6.951 eV/Å2 and 11.976 eV/Å2 provide the best fit to our DFT energetics for fourfold-coordinated vacancy and interstitial defects, respectively. Note that a smaller value of force constant indicates less stiffness in the defect structure. 64 Annealing of vacancy defects may be preceded by dissociation into smaller ones in a vacancy-rich region. Thus, to examine their thermal stability we also calculate the binding energies of the fourfold-coordinated and the PHR vacancy clusters. The results are summarized in Figure 3.6(b). Here, the binding energy is given by: ܧ௕ሺ݊ሻ ൌ ܧ௙ሺ݊ െ 1ሻ ൅ ܧ௙ሺ1ሻ െ ܧ௙ሺ݊ሻ, (3.4) which represents an energy cost for single vacancy liberation from a given cluster, i.e., Vn → Vn-1 + V. For the PHR case, the binding energies show an oscillating behavior with local maxima at n = 5, 10, 14 and 18, in good agreement with previous studies [91], except for V5 and V6. The discrepancy is attributed to the fact that for V5 we consider the stable fourfold-coordinated PHR state converted from the conventional PHR (vide supra), leading to increase the binding energy of V5 while lowering the V6 binding strength accordingly. For the fourfold-coordinated vacancy clusters, their binding energies are similar in magnitude to those of the PHR clusters, although they exhibit local maxima and minima at different sizes. Figure 3.7(a) summarizes the predicted formation energies of PHR and FC vacancy clusters up to n = 48. The total energy differences between PHR and FC defects are presented in Figure 3.7(b), which clearly demonstrates that the FC configurations are energetically more favorable than the PHR configurations for 3 ≤ n ≤ 48. The energetically-favored formation of FC defects indicates that the bond energy gain by fourfold-coordination exceeds the strain energy increase via consequent lattice distortions. 65 Figure 3.7: (a) DFT-predicted formation energies per vacancy (Êf) of vacancy clusters as a function of cluster size (n) for both fourfold-coordinated (indicated as FC) and “part of hexagonal ring” (PHR) configurations, together with adjusted values (FC*) for the FC configurations. For the PHR case of n > 18, the DFT formation energies per vacancy are well fitted with the power function of (−0.82 + 3.68n2/3)/n, as indicated with a dotted line. (b) Variation in the total energy difference between the FC and PHR structures from DFT calculations (indicated as ΔE), together with adjusted values (ΔE*). Refer to the text for the adjustment in the FC defect formation energies using FF calculations, which is associated with incomplete relaxation of the defect system due to limited supercell size used for DFT calculations [9]. 66 For the DFT calculations, 480–n, 576–n, and 1000−n atom supercells were used for V1-V14, V15-V18, and V19-V48, respectively. For the PHR case, the atomic structures of n = 3-18 from previous studies [91] were recalculated within DFT-GGA while the other structures of n > 18 were constructed by optimizing initial structures which are prepared to have minimum dangling bonds. For the formation energy calculation, highly strained FC structures commonly require larger supercells than PHR structures. To avoid large and computationally expensive supercells, we adjusted the formation energies [Ef(n) = nÊf(n)] for the DFT calculation with Ef(n) from the force field (FF) calculations (as described in Ref. 8): ܧ௙ሺ݊ሻ ൌ ܧ௙ிிሺ݊, ௅ܰሻ ൅ ሼܧ௙஽ி்ሺ݊, ௌܰሻ െ ܧ௙ிிሺ݊, ௌܰሻሽ, (3.5) where superscripts FF and DFT indicate values calculated from the respective methods, respectively, while NL and NS represent large (27000 atoms which is 15×15×15 times larger than a 8-atom cubic cell) and small (≤ 1000 atoms employed for DFT calculations) supercells, respectively, both without defects. The discrepancies between the DFT (indicated as FC) and adjusted (FC*) values in Figure 3.7(a) indicate the importance of the adjustment. We also carefully checked the validity of the adjustment approach. As illustrated in Figure 3.8, the DFT [ܧ௙஽ி்ሺ݊, ௌܰሻ] and FF [ܧ௙ிிሺ݊, ௌܰሻ] values for the FC Vn defects in these small (NS ≤ 1000 atoms) supercells are almost identical, implying that the FF method is adequate to estimate the formation energies of FC defects. This is not surprising considering the well-relaxed FC structures with no significant departure from 67 crystalline Si, and thus their energetics can be described reasonably well by Keating-like potentials [5-8,10-12]. The inset in Figure 3.8 exhibits the variations of ܧ௙ிிሺ݊, ௅ܰሻ for FC V12 and FC V32 with increasing NL. The result demonstrates that the ܧ௙ிிሺ݊, ௅ܰሻ values nearly saturate when NL  27000 while the smaller case (V12) converges much faster, confirming the soundness of our choice of NL = 27000 atoms for the adjustment of the formation energies of FC defects. Figure 3.8: Comparison of the FC cluster formation energies (Ef) from between DFT (indicated as DFT) and FF (as detailed in Section 3.2) calculations. The inset shows the Ef variations of FC V12 and V32 in terms of supercell size; here N indicates the number of Si lattice atoms prior to vacancy defect creation [9]. 68 One important feature from our calculations is that the FC defect formation energy appears to saturate for n > 30 (although some variation with n still occurs), while the PHR defect formation energy monotonically decreases. As a result, the total energy difference between the FC and PHR states peaks at n = 32, then gradually decays; for n >> 50, the PHR state is expected to be energetically more favorable than the FC state. This further supports our previous conclusion [8] that the Ef(n) of void-like defects (such as PHR defects) is largely governed by the void-surface energy which is proportional to void-surface area (~ n2/3), whereas the Ef(n) of FC defects is determined by the number of strained Si atoms which is proportional to cluster size (~ n). Therefore, we expect the formation energy per vacancy in the PHR case to monotonically decrease (~ n2/3/n = n−1/3), while the FC case instead saturates and becomes size-independent (~ n/n = 1). 3.5 Extremely Stable FC V12 and V32 Among all the defect configurations considered, the FC V12 turns out to be the most favorable thermodynamically, with a binding energy of 5.36 eV. As shown in Figure 3.9(a), the V12 structure has C2 symmetry with a combination of four identical structural units. The structure of each unit consists of two five-, two six-, and two eight- membered rings, resembling a small void. The V12 structure can be viewed as the V17 SPC with a well-relaxed tetragon inside. (The V17 SPC is obtained by removing up to the second nearest neighbors around a center atom). It is worth pointing out that the six bond angles around the center atom (as indicated) are 109.5o, identical to the equilibrium 69 tetrahedral bond angle of crystalline Si. The well relaxed fourfold-coordinated structure can also be evidenced by the Wannier centers which are almost-perfectly located at bond mid-points [Figure 3.9(a)] as well as the total density of states which appears virtually identical to that for defect-free Si [Figure 3.9(b)]. Having significant stability, the FC V12 cluster could be expected to exist to a large extent in a vacancy rich region, although its direct characterization appears impractical at present. Another, perhaps most important, observation is the existence of an extremely stable FC V32 configuration with a Êf(n) of only 0.89 eV and 0.83 eV before and after adjustment, respectively. To better understand its remarkable stability we analyze the atomic configuration of FC V32 in Figure 3.10. The FC V32 structure contains an “adamantine cage” (AC) core surrounded by a distorted shell and exhibits tetrahedral geometry with four triangular (111) facets and four corners. The covalent network inside the shell region is distorted and under tension, but also shows a characteristic pattern of n-membered rings consistent with the overall cluster symmetry. The six edges along the tetrahedron boundary are each defined by two essentially parallel and approximately linear chains of five Si atoms. Figure 3.10 illustrates the subtle differences among the six cluster edges, while Figure 3.10(b) conveys the overall symmetry of FC V32 and indicates that the six edges of FC V32 are not all identical. We define four edges as parallel and two edges as antiparallel. The atomic-level detail of Figure 3.10(a) and (c) contrasts how the two edges are different. In the antiparallel edges, the bonded atoms among the five constituents of each chain are staggered. In the parallel edges, the bonded atoms among the five members of each chain are arranged side-by-side. Note that the alignment of 70 unbonded atoms between the two parallel chains is responsible for an 8-atom ring on each of these edges. If the linear arrangement of one of the five-membered chains of the antiparallel edge is reversed, a parallel edge configuration results. Despite the variations in edge structure, the four (111) faces and four corners are identical to each other, respectively. Therefore, the approximately regular tetrahedron, which has one C2 axis and two S4 rotation-reflection axes, is characterized by S24 symmetry. Figure 3.9: (a) Ball-and-stick illustration of the atomic structure of V12, together with Wannier function centers as indicated by small black balls. Grey (gold) balls indicate more distorted atoms than the rest of the lattice atoms (in white), while five dark grey (red) balls indicate a well relaxed tetragon at the defect center. (b) Calculated total density of states (TDOS) for the c-Si system with (upper panel) and without (lower panel) the V12 cluster. Here, EF indicates the Fermi level [8]. 71 Figure 3.10: (a) Core and strained nearest neighbors of FC V32 shown isolated from bulk Si lattice. In (a), we annotate the antiparallel edge type, where member atoms of each chain are distinguished by either “x” or “†” labels (blue or red spheres, respectively). The ten atoms comprising the adamantine cage (AC) in the center of the structure are identified as dark grey spheres. (b) Illustration of a simplified geometric diagram of FC V32 that highlights the S24 group symmetry, AC core, and the two antiparallel (double line) and four parallel (single line) cluster edges. Axis labels denote cluster orientation with respect to the Si lattice. (c) Reorientation of FC V32 to emphasize the parallel edge, where atoms comprising these chains are annotated with “x” labels (blue spheres). (d) Total density of states comparison of bulk Si to V32 near the band gap using 1000 atom basis supercells with the conduction band edges set as references [9]. 72 It is worth emphasizing that each of the four corner regions of the FC V32 tetrahedron resembles the C2 symmetry FC V12 structure that is comprised of four identical void-like structural units surrounding a tetragonal core [8]. Figure 3.10(d) shows that the band gap in the presence of V32 is free of additional defect states and is essentially the same as the c-Si gap. The enhanced stability of FC V32 is likely attributed to its unique crystalline core and symmetric, yet distorted, covalent bond topology in the shell region. We also find that when n > 25 the stable FC clusters identified here commonly show the core/shell-like arrangement (with an AC core), while smaller FC clusters (10 < n < 25) exhibit the trace of ordered V12 [8]. The atomic configurations of other FC clusters will be presented elsewhere. 3.6 Kinetic Effects on the Formation of Fourfold-coordinated Vacancy Defects Next, we examine if the formation of FC vacancy defects can also be kinetically favored over their PHR counterparts. Vacancy defects may grow by capturing mobile single vacancies (and possibly also di-vacancies [85,87]) from the crystalline matrix. In addition, vacancy clusters may also be formed by recrystallization of low-density amorphous regions that are deficient in Si atoms compared to the perfect crystal. To investigate vacancy cluster growth, we performed TBMD simulations at 1400K for 73 several systems of small clusters (Vn, n = 2-4, 11), each of which had an additional mono- vacancy in close proximity. Figure 3.11: Mono-vacancy induced structural rearrangement, together with variation in the formation energy (Ef = nÊf) as a function of annealing time for V3 + V → V4 during TBMD simulations at 1400K. All atomic configurations from the TBMD simulation were optimized within DFT-GGA calculations using one Gamma k-point sampling for 512−n atom supercells and illustrated in the top and middle panels. Grey wire frame represent bulk Si lattice and dark grey spheres represent highly strained FC atoms while brown spheres represent a Si mono-vacancy. Si dangling defects are also depicted by red (black) spheres [9]. 74 Figure 3.12: Mono-vacancy induced structural rearrangement, together with variation in the formation energy (Ef = nÊf) as a function of annealing time for V11 + V → V12 during TBMD simulations at 1400K. All atomic configurations from the TBMD simulation were optimized within DFT-GGA calculations using one Gamma k-point sampling for 512−n atom supercells and illustrated in the top and middle panels. Grey wire frame represent bulk Si lattice and dark grey spheres represent highly strained FC atoms. Si dangling defects are also depicted by red (black) spheres [9]. 75 Figure 3.13: Mono-vacancy induced structural rearrangement, together with variation in the formation energy (Ef = nÊf) as a function of annealing time for V4 + V → V5 during TBMD simulations at 1400K. All atomic configurations from the TBMD simulation were optimized within DFT-GGA calculations using one Gamma k-point sampling for 512−n atom supercells and illustrated in the top and middle panels. Grey wire frame represent bulk Si lattice and dark grey spheres represent highly strained FC atoms. Si dangling defects are also depicted by red (black) spheres [9]. 76 Figure 3.14: Structural evolution of amorphous pockets with deficit of 12 atoms during TBMD simulations at 1400K together with variation in the formation energy (Ef = nÊf) as a function of annealing time. All atomic configurations from the TBMD simulation were optimized within DFT-GGA calculations using one Gamma k-point sampling for 512−n atom supercells. Filled circles illustrate the formation of the lowest-energy FC V12 with C2 symmetry, together with their atomic configurations as indicated. Open circles illustrate the formation of a PHR V12 structure, together with their atomic configurations as indicated; however, the PHR defect formation would be rather unlikely in a typical annealing process, as discussed in the text. Grey wire frame represent bulk Si lattice and dark grey spheres represent highly strained FC atoms. Si dangling and floating bonds are depicted by red (black) and blue (black, pointed also by the arrow) spheres, respectively [9]. 77 Figure 3.11 illustrates V3 + V → V4, where the initial FC structure of V3 is quickly perturbed upon capturing a mono-vacancy and yields multiple coordination defects (CDs), i.e., dangling and floating bonds. The intermediate structure undergoes successive rearrangements that render a stable PHR structure around 12.5 ps, which subsequently relaxes to the lowest-energy FC structure with C2h symmetry within 15 ps. Similarly, as shown in Figure 3.12, the V11 + V → V12 case leads to the lowest-energy FC V12 structure with C2 symmetry within 7 ps. During the agglomeration, multiple CDs arise and diffuse around the cluster, which in turn enhances structural rearrangement while the majority of Si atoms in the cluster region remain fourfold-coordinated. For the V2 + V → V3 and V4 + V → V5 cases, the overall agglomeration mechanisms are similar to the V3 + V → V4 and V11 + V → V12 cases. However, after 20 ps of the TBMD simulations, the V3 and V5 aggregates did not complete fourfold- coordination. We continued the TBMD simulations for additional 20 ps (= 40 ps total). As shown in Figure 3.13, the V5 ended up with the lowest-energy FC configuration around 35 ps, but the V3 was not fully fourfold-coordinated after 40ps of the TBMD simulation. The relatively slow reconfigurations of the V3 and V5 aggregates to their lowest-energy FC configurations can be related to the small energy differences of 0.26 eV and 0.33 eV, respectively, between the FC and PHR states; in contrast, the fourfold- coordination of V4 and V12 yields significant energy gains of 1.36 eV and 3.99 eV (before adjustment) over their PHR counterparts. In particular, the V3 FC structure is highly strained and thus can be easily ruptured with generating CDs during the high temperature 78 annealing at 1400K, retarding the fourfold-coordination that leads to its lowest energy state. Using TBMD simulations at 1400 K, we proceeded to evaluate the structural evolution of amorphous pockets with atomic Si deficiencies. We employed CRN-based MMC simulations to prepare initial amorphous pocket structures embedded in c-Si supercells. For each model amorphous system in our statistical sample, twelve atoms were deficient in the amorphous region relative to c-Si, but all constituent Si atoms remained fourfold-coordinated. Consequently, complete recrystallization will leave twelve vacancies. Among four independent model systems evaluated, one case yielded the lowest-energy FC V12 with C2 symmetry within 20 ps as shown in Figure 3.14 (dark circle). Upon annealing, the amorphous pocket is relaxed via bond rearrangements assisted by coordination defects and subsequently leads to its recrystallization from the interface between the amorphous and crystalline regions. After 17 ps, the FC V12 configuration is formed and remains virtually unchanged throughout the remainder of the 60 ps simulation. For another two cases, the amorphous pockets gradually transformed into the lowest energy FC V12 state, but not completed within the 60 ps simulation duration. As shown in the inset of Figure 3.14, our TBMD simulations also show the possibility that an amorphous pocket would evolve into a PHR defect. We noticed that the PHR defect formation can take place if the rate of CD creation by breaking highly strained Si-Si bonds is higher than that of CD annealing through bond network rearrangements (i.e., when several Si-Si bonds are ruptured for a short period of time 79 within which annealing of the created CDs seldom occurs). The nearly simultaneous rupture of multiple Si-Si bonds could be possible only when the amorphous network is significantly strained and also under immediate high-temperature treatment, like the condition employed in the TBMD simulation above. However, a typical annealing process involves a ramp-up period during which we expect a highly-strained network would be gradually relaxed without the nearly simultaneous break of multiple Si-Si bonds. This leads us to speculate that the formation of open-volume PHR-type would be rather unlikely during recrystallization of an amorphous region, especially when the amorphous pocket is large enough that it can be rearranged without inducing significant local strains. While the creation of PHR defects cannot be excluded during processing, for the sake of completeness we also performed TBMD at 1400 K starting with the PHR V12 configuration, but there was no indication for PHR  FC conversion although dangling bonds diffuse along the cluster surface. This suggests that there exists a sizable barrier for the PHR  FC transformation for sufficiently large vacancy defects, which is not surprising considering the significant lattice distortions required for linking the under- coordinated Si atoms to be fourfold-coordinated. The results may also indicate that PHR defects could coexist with FC defects while the latter are likely to be prevailing. Although further studies might be necessary to verify the kinetic role in the structural evolution of various Vn clusters, our work suggests the possible formation of thermodynamically-favored FC vacancy clusters for n < 48 with no significant kinetic limitations. This suggests that a large fraction of small vacancy defects can exist in the FC state in c-Si. However, it appears to be impractical to explicitly characterize FC 80 vacancy clusters with currently available spectroscopic techniques because of their optically inactive nature (given no energy states within the Si band gap [see Figure 3.10(d)]). Nonetheless, the presence of FC defects could be expected to alter the properties and behavior of Si materials such as thermal conductivity, elastic modulus, and impurity gettering. Further studies are therefore warranted to better understand the formation and nature of FC vacancy defects together with open volume PHR-type defects, particularly their effects on the thermal/electrical transport, mechanical, and gettering properties of Si-based nanomaterials for future thermoelectric and electronic applications. 3.7 Summary Based on combined continuous random network (CRN) model, tight binding molecular dynamics (TBMD), and density functional theory (DFT) calculations, we present stable fourfold-coordination for small vacancy clusters in the size range of 3-48 vacancies. For each cluster size, we first generated possible fourfold-coordinated vacancy clusters using CRN based Metropolis Monte Carlo (CRN-MMC) simulations, followed by TBMD simulations at high temperatures to check the stability of the fourfold-coordinated structures. Then, DFT calculations were performed to refine the geometries of the stable clusters, and compared their formation energies to determine the lowest-energy structure among them. In this size regime our calculations demonstrate that FC configurations are energetically more favorable than PHR structures which have been until now considered to be prevailing; however, with increasing cluster size (up to a 81 few nanometers in diameter) less strained open volume PHR-type defects will be ultimately favored as the energy penalty from strained atoms in FC clusters overwhelms energy gains from minimized dangling bond density. The preference for fourfold- coordinated structuring is apparently attributed to the fact that the energy gain by bond formation exceeds the strain energy associated. Our results highlight very stable FC V12 and V32 structures. The FC V12 consists of four identical structural units while each unit has two five-, two six-, and two eight- membered rings. The FC V12 structure is predicted to be about 4 eV more favorable than the conventional PHR structure. The FC V32 configuration exhibits a complex but ordered tetrahedral shape containing an “adamantine cage” core surrounded by a distorted shell. When n > 25 stable FC clusters commonly show the core/shell feature, while smaller FC clusters (10 < n < 25) exhibit the trace of ordered V12 that is comprised of four identical void-like structural units surrounding a tetragonal core. Our TBMD simulations also suggest that the formation of thermodynamically favored FC vacancy clusters could be kinetically facile. While small vacancy defects thermodynamically favor fourfold-coordination in crystalline Si, our theoretical study also demonstrates that the PHR configuration will become energetically more favorable than the fourfold- coordinated configuration when the defect size is greater than a couple of nanometers in diameter. Given the significant stability, we expect that the FC V12 and V32 defect would exist to a large extent in a vacancy rich region. Our findings emphasize the importance of better understanding the nature of stable FC vacancy defects to effectively estimate and manipulate the properties and behavior of Si-based materials. The improved 82 understanding regarding the structure and stability of vacancy defects will greatly assist in better understanding the properties of native defects in Si, and explaining and predicting their impact on relevant materials properties and processing. 83 Chapter 4 Strain Effect on the Structure and Stability of Native Defects in Silicon 4.1 Introduction Strain engineering has received intense attention in the semiconductor industry in the past few years as a vehicle to extend silicon complementary metal-oxide- semiconductor (CMOS) transistor performance in modern high-performance electronics [10]. Relative to other CMOS enhancement techniques, incorporation of process strain only adds a few percent to wafer cost since only a few extra steps are required and often existing steps can simply be modified by tuning the strain of deposited thin films [99]. As a result, the semiconductor industry began widely incorporating process-induced strain into manufacturing flows with the introduction of the 90 nm node as a cost effective technique to help extend transistor performance improvement consistent with Moor’s Law [10]. There are two ways of technological relevance to apply strain to the channel of a metal-oxide semiconductor field-effect transistor (MOSFET): biaxial strain, which is sometimes called global or bulk strain because it is implemented at the substrate level, and uniaxial strain, which is sometimes referenced as local or process-induced strain [100]. Tensile biaxial strain can be implemented in a transistor channel by building the entire device in strained Si epitaxy over a thick Si1-xGex layer. Uniaxial strain is often 84 implemented at the device level in two ways: (1) selective epitaxial growth of a binary Si alloy (Si1-xGex for channel compression and Si1-xCx for channel tension) in the source/drain recessed regions to impose strain along the transistor channel direction or (2) deposition of a high-stress silicon nitride cap layer which mechanically couples the local strain of the film into the underlying transistor channel. Hydrostatic strain can be achieved in bulk Si simply by subjecting the material to a uniform pressure field, but this form is technologically less relevant than biaxial or uniaxial strain [10]. In this work, we examined how the presence of biaxial strain influences the stability of small interstitial and vacancy clusters. 4.2 Computational Approach In epitaxial growth, biaxial strain occurs in an epilayer when lattice mismatch is present between the epilayer and the substrate. If the lattice constant of the epilayer, a, is smaller than the lattice constant of the substrate, a0, then tensile biaxial strain results. Conversely, if a is larger than a0, compressive biaxial strain results. It is useful to define the lattice constant in the plane of strain as a║ and the lattice constant perpendicular to the plane of strain as a┴. In the regime of linear elastic behavior we model for Si, when a║ expands under tensile strain, a┴ simultaneously contracts. Under compressive strain, a║ contracts while a┴ expands. This phenomenon is known as the Poisson effect, as visualized in Figure 4.1 [101]. 85 Figure 4.1: Tensile biaxial stress/strain interaction in our model Si supercell. In the figure, applied tensile stress, σ║, in the plane of the substrate acts equally in all directions as shown by block arrows and produces a tensile strain. In response, the lattice contracts in the out-of-plane direction as shown by the solid black arrows. Under compressive strain conditions, the directions of all arrows are inverted [10]. We can quantify the Poisson effect for biaxially strained systems by defining a quantity, ν*, that relates the ratio of in-plane strain, ε║, and out-of-plane strain, ε┴. In our system, the values of a║ in Si under tensile strain conditions are equal to representative values of aSiGe, which is the lattice constant of a binary SiGe system. We calculate the in- plane strain as ε║ = (aSiGe - aSi)/aSi and the out-of-plane strain as ε┴ = (a┴ - aSi)/aSi. The experimental value of aSi is 5.4309 Å and aGe is 5.6461 Å, [102] so 4% tensile strain is the limiting case of Si grown over pure Ge. From linear elastic theory, [103, 101, 102, 104] the relationship between out-of-plane and in-plane strain for a cubic crystal can be expressed in terms of two elastic stiffness constants, 86 ߥכ ൌ െ ఌ఼ఌצ ൌ 2ሺ ஼భమ ஼భభሻ. (4.1) Quantity ν* is valid for deformation of a few percent strain. Using the values of 16.6×1011 and 6.4×1011 dyn/cm2 for C11 and C12, [105] respectively, the value of ν* is 0.771. Using ν* and the expression for ε║ and ε┴, we calculated the values of a┴ for each independent value of a║ studied. All results we present here for biaxially strained Si are based on ν* = 0.771. By iterating through a reasonable range of a┴ for each value of a║ studied, we numerically verified that minimum-energy supercell dimensions occur as ν* converges to 0.771. The formation energy in terms of cluster size (n) and strain condition (ε), Ef(n, ε), is given by ܧ௙ሺ݊, ߝሻ ൌ ܧ௧௢௧ሺ݊, ߝሻ െ ௡ାேே ܧ௕௨௟௞ሺߝሻ, (4.2) where ܧ௧௢௧ሺ݊, ߝሻ is the total energy of the In cluster in the n + N atom supercell, n is the size of the interstitial cluster, N is the basis number of atoms in the bulk Si supercell, and ܧ௕௨௟௞ሺߝሻ is the total energy of the N atoms supercell of crystalline Si at a given biaxial strain condition. [10] 4.3 Biaxial Strain Effects on the Structure and Stability of Self- Interstitial Clusters in Silicon To examine how the presence of biaxial strain influences the stability of small interstitial clusters, we applied the integrated atomistic modeling procedure under 4% 87 compressive and 4% tensile biaxial strain to their ground-state configurations (Ing, n < 12) under the strain-free condition. We first focused on I3 and I4 clusters as these clusters are seen to function as core building blocks for larger clusters under the strain-free condition. Figure 4.2: Predicted lowest energy configurations for I4 under 4% compressive (top), strain-free (middle, I4g), and 4% tensile (bottom, I44%t) biaxial strain. Biaxial strain is applied to the supercell as shown in Figure 4.1. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated. 88 Figure 4.2 shows the lowest-energy configurations for I4 formed under strained conditions. Under both 4% compressive and 4% tensile biaxial strain conditions, I4 yields the same configuration (I44%c and I44%t) as the I4g, but I44%c and I44%t are aligned along the different orientations. Figure 4.3 shows lowest-energy configurations for I3 formed under strained conditions. Under 4% compressive strain, I3 yields a previously unreported configuration (I34%c). The I34%c configuration has D2d symmetry, like the I4g core structure, and also exhibits similar orientation-dependent strain behavior as previously detailed for I4 [10]. Most interstitial clusters are stabilized under tensile conditions because the lattice becomes more accommodating to interstitial atoms as it stretches. The ܫଷସ%௖configuration is the only Si interstitial cluster investigated that is stabilized as strain condition becomes more tensile. Under 4% tensile strain, I3 yields the same geometry (I34%t) as the I3g but is aligned along the different orientation. Figure 4.4 shows lowest-energy configurations for I7 formed under strained conditions. Under 4% compressive strain, I7 yields a previously unreported configuration (I74%c). Interestingly, the I74%c is a combination of I44%c and I34%c. Likewise, under 4% tensile strain, I7 yields a combination of I44%t and I34%t, i.e., the same configuration as I74%g. Under biaxial strain conditions, I7 exemplifies the significance of the I3 and I4 cores in the formation of larger clusters. In addition, the orientations (I4 and I3) and configurations (I3) of the constituent cores of I7 under biaxial strain follow the energetic trends seen for the isolated core components. 89 Figure 4.3: Predicted lowest energy configurations for I3 under 4% compressive (top, I34%c), strain-free (middle, I3g), and 4% tensile (bottom, I34%t) biaxial strain. Biaxial strain is applied to the supercell as shown in Figure 4.1. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated. 90 Figure 4.4: Predicted lowest energy configurations for I7 under 4% compressive (top, I74%c), strain-free (middle, I7g), and 4% tensile (bottom, I74%t) biaxial strain. Biaxial strain is applied to the supercell as shown in Figure 4.1. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated. 91 Figure 4.5 shows lowest-energy configurations for I8 formed under strained conditions. Under 4% compressive strain, I8 yields a combination of two I44%c clusters and exhibit the same configuration as the I8g. Likewise, under 4% tensile strain, I8 yield a combination of two I44%t clusters and exhibit the same configuration as I84%g, but aligned along a different orientation. Figure 4.5: Predicted lowest energy configurations for I8 under 4% compressive (top), strain-free (middle, I8g), and 4% tensile (bottom, I84%t) biaxial strain. Biaxial strain is applied to the supercell as shown in Figure 4.1. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated. 92 Our DFT calculations reveal that both I4 and I8 become more stable under certain strain conditions [10]. The Ef(4, ε) dependence of I4 on biaxial strain is determined by the orientation of I4 with respect to the biaxial strain field as shown in Figure 4.6 [the Ef(8, ε) response of I8 is nearly identical]. Motivated by this interesting behavior observed for I4 core derivative configurations, we proceeded to search for a potential family of extended configurations composed of I4 core units using the integrated atomic modeling procedure [5-7] under strained conditions. Starting with the I8 ground state configuration, four interstitials were added in the vicinity and a conditionally stable I12 structure was discovered as depicted in Figure 4.7. Under 4% compressive strain, I12 yields a combination of three I44%c clusters. Likewise, under 4% tensile strain, I12 yield a combination of three I44%t clusters. In order to emphasize the intraconfigurational composition of I4 cores, we will refer to the I4 chain configurations as (I4)m, where m represents the number of I4 cores in the chain. Our results suggest that the (I4)m configurations may participate in the compact-to-extended transition regime (10 ≤ n ≤ 20) of self-interstitial cluster growth under sufficient strain conditions [11]. 93 Figure 4.6: Formation energy response per interstitial compared for the two relevant orientations of I4 using 256+n supercells. Perspective views along the [001] are shown for reference along with corresponding orientations of the S4 rotation-reflection axes. Light gray (gold) wireframe represents the bulk Si lattice. Dark gray spheres denote interstitials and their highly strained neighbors [11]. 94 Figure 4.7: Predicted lowest energy configurations for I12 under 4% compressive (top, I124%c), strain-free (middle, I12g), and 4% tensile (bottom, I124%t) biaxial strain. Biaxial strain is applied to the supercell as shown in Figure 4.1. Grey (gold) balls indicated more distorted atoms than the rest of the lattice atoms (in white). The symmetry for each cluster is also indicated. 95 4.4 Strain Effects on the Stability and Structure of Vacancy Clusters in Si The effect of strain on the stability and structure of small, neutral Si vacancy clusters (Vn, n ≤ 12) was investigated using periodic DFT calculations. Our results indicate that compressive strain generally stabilizes vacancy clusters, which is complementary to the conclusions of our previous work [10-12] that showed that tensile strain stabilizes Si interstitial clusters. The magnitude of stabilization provided by biaxial compression is generally greater than the magnitude of stabilization observed under uniaxial compression. Similar to our previous work of interstitials, we also observe orientation effects for vacancy clusters in uniform strain fields. The C2h-symmetry configurations of V4 [see Figure 3.1(b)] and V8 [see Figure 3.2(b)] exhibit orientation-dependent behavior under biaxial strain in correlation with the orientation-dependent behavior shown by interstitial clusters with C2h symmetry. This observation is significant because it emphasizes the critical role that cluster symmetry performs in orientation-dependent stabilization in strained Si, regardless of cluster composition. Our calculations show that strain conditions can modulate the relative stability of competing fourfold-coordinated (FC) and PHR-type configurations. While FC configurations are more favorable for small vacancy clusters in strain-free environments, PHR-type configurations become preferred under certain tensile conditions in both biaxial and uniaxial cases [13]. The relative stabilization of PHR-type configurations 96 under tension is largely attributable to the flat nature of the formation-energy dependence on strain that is generally exhibited by clusters with incomplete fourfold-coordination. In addition, our results illuminate a general orientation effect for PHR-type configurations (2 ≤ n ≤ 6) under uniaxial strain related to alignment of the parent ring structure with respect to the direction of applied strain. Our simulations also demonstrate that strain of sufficient magnitude can initiate configuration changes. We provide two examples, V3t (see Figure 4.8) and V5t (see Figure 4.9), of vacancy configurations we identified by cluster formation under 3% tensile biaxial strain. These particular structures were highlighted because they exhibited different formation energy responses as a function of strain when compared to their counterpart configurations formed under stain-free conditions. Our study furthers the understanding of Si vacancy clusters by providing insight into the influence that strain has on the stability and structure of vacancy defects, which is an important step toward property prediction and ultimately defect engineering in Si-based materials [13]. 97 Figure 4.8: Formation energy per vacancy as a function of biaxial strain for the strain- free, ground-state FC V3 structure and a FC tri-vacancy configuration (V3t) identified by formation under tensile biaxial strain conditions. Inset cluster configurations are shown along [11ത0]. Light gray (gold) wireframe represents bulk crystalline Si. Dark gray spheres represent highly strained atoms neighboring the vacancy clusters [13]. 98 Figure 4.9: Formation energy per vacancy as a function of biaxial strain for the strain- free, ground-state FC V5 structure and a FC penta-vacancy configuration (V5t) identified by formation under tensile biaxial strain conditions. These structures are approximately degenerate under strain-free conditions. Inset cluster configurations are shown along [11ത0]. Light gray (gold) wireframe represents bulk crystalline Si. Dark gray spheres represent highly strained atoms neighboring the vacancy clusters [13]. 99 Chapter 5 Ab Initio Parameterized Valence Force Field for the Structure and Energetics of Amorphous SiOx (0 ≤ x ≤ 2) Materials 5.1 Introduction Amorphous silicon-rich oxides (a-SiOx, 0 ≤ x ≤ 2) have garnered great attention not only from their unique properties, but also potential technological importance in various electronic, optoelectronic, and energy applications. Amorphous phases lack long range order or well-defined atomic structure and are thermodynamically less favorable than their corresponding crystalline phases [106,107]. As a consequence of lattice distortion and variations in composition, a-SiOx materials typically have different properties from their crystalline counterparts and also Si/SiO2 multiphase systems. The properties of a-SiOx can be controlled by varying Si:O composition ratio and incorporating various impurities. For instance, hydrogenated a-SiOx materials exhibit visible room temperature photoluminescence (PL) properties and a higher photosensitivity than other materials with comparable optical gaps, making them viable candidates for Si-based optoelectronic applications [108-111]. In addition, doping with boron or phosphorous atoms might enable realization of light-emitting diodes [64,112- 114]. Si suboxides may undergo phase separation to yield oxide-embedded Si nanoparticles during high temperature annealing; the Si nanoparticles (np-Si) can be 100 amorphous or crystalline depending on the annealing temperature [115-119]. The np- Si/a-SiO2 system emits visible PL with high efficiency at room temperature, while luminescence from bulk crystalline Si (c-Si) is negligible as a result of its indirect band structure. The discovery of Si-based luminescence has generated considerable interest in its potential application to integrated optoelectronic devices [120-124]. In addition, oxide-embedded Si nanoparticles have been envisioned as possible discrete storage elements in non-volatile flash memories [125-131]. Earlier studies have suggested that the performance of np-Si based devices would be determined by a complex combination of the following attributes: Si particle size, shape, and crystallinity; Si/SiO2 interface structure and strain; and near-interface defects (bonding, chemical, and structural) [132]. It is therefore important to develop a detailed understanding of the structure, strain, and composition of a-SiOx materials and their interfaces. An atomic level understanding of a-SiOx materials derived from experimental methods has thus far remained elusive in part from the limited resolution of common instrumentation. A complementary effort has been made in the development of nanoscale atomistic models of Si and its amorphous oxides (a-SiOx, 0  x  2), where first-principles methods have achieved widespread usage in characterization of the structure and stability of these systems. Despite substantial growth in computational power, most technologically-important systems containing Si and/or O, such as a-Si, a- SiO2, c-Si/a-SiO2, and np-Si/a-SiO2, are cost-prohibitive to address exclusively with first- principles calculations because realistic atomic structural models typically contain complex topologies requiring large numbers of atoms. As an alternative to first- 101 principles calculations, computationally less expensive empirical interatomic potentials have been widely used; various types of interatomic potentials have been developed for Si [66,133-156], SiO2 [157-163], and Si/SiO2 composites [1, 164-166]. The primary liability in application of empirical potentials is their limited transferability. Even for Si, no single empirical potential provides an adequate description of Si-Si interaction across all relevant forms (bulk crystal, amorphous and liquid phases, surfaces, and clusters) of this prototypical semiconductor; consequently, it is reasonable to conclude that generation of a single potential that comprehends all pertinent phenomena is likely an insurmountable task. Therefore, we believe it is crucial to develop application-specific potentials that can be robustly-applied to describe atomic structure in designated systems. Experimentally, a-Si [107,167,106] and a-SiO2 [168,169] are known to form random tetrahedral networks characterized by both long- range disorder and short-range order similar to that of their parent crystals. At c-Si/a- SiO2 interfaces, electrical and electron-spin-resonance measurements show extremely low densities (typically between 1010 to 1012 cm-2) [170-172] of interface defects consistent with an almost-perfect bonding network across the interface. Therefore, typical a-SiOx systems can be well-represented by fully-coordinated random tetrahedral networks in which Si and O atoms have four and two bonds, respectively, and O-O bonding is ideally absent. To generate a structural model of a-SiOx, energy minima identification at the atomic level is a fundamental first step. To this end, molecular dynamics (MD) or Monte Carlo (MC) methods coupled to various interatomic potentials have been widely used as 102 detailed in the following work: (1) ab-initio MD for a-Si [173], a-SiO2 [174-176], and planar c-Si/a-SiO2 [177]; (2) classical MD for a-Si [156], a-SiO2 [162,178,179], planar c- Si/a-SiO2 [164,166,180-182], and np-Si/a-SiO2 [206]; and (3) classical MC for a-Si [65,66], a-SiO2 [157], a-SiOx [183], planar c-Si/a-SiO2 [1,184,185], and np-Si/a-SiO2 [3,186-187]. Ab-initio MD permits accurate description of atomic motion, but its utility is restricted to small systems and short time scales because of steep computational requirements. Classical MD permits simulation of relatively large systems, but the same time scale limitations can compromise complete structural relaxation and the availability of reliable force fields for bond formation/scission is also a concern. For our work, classical MC using a CRN model is a proven approach for the construction of fully- relaxed a-SiOx structures [168]. Within the CRN model, an amorphous system is relaxed via a large number of bond transpositions using Metropolis Monte Carlo (MMC) sampling [65], where the validity of the final structure also strongly depends on application of a reliable force field. Since this approach does not require description of bond formation/scission, simple and computationally less expensive valence force field (VFF) models like the three-body, harmonic Keating-like (KT) potentials have been widely used, permitting simulation of larger systems. Overall, the VFF models are effective in prediction of minimum-energy configurations of fully-coordinated Si-based materials when the bond lengths and angles do not significantly deviate from their equilibrium values [5-8,14]; however, relatively little optimization of a-SiOx VFF models is currently available in the literature. 103 In this work, we focus on construction of a semiempirical VFF model which accurately depicts the energetics of fully-coordinated a-SiOx systems (for a wide range of local distortion from ground state tetrahedral coordination), but excludes bond formation and scission. The formulation of our VFF model is based on the modified Keating model proposed by Tu and Tersoff [1], but adjusted to fit DFT total energy calculations for various atomic structures. We investigate how the atomistic description of a-SiOx is affected by force field selection by performing CRN-MMC simulations utilizing the original Tu and Tersoff potential and our modified potential to construct various amorphous systems including bulk, homogeneous a-SiOx (x = 0, 0.5, 1.0, 1.5, and 2.0) and multiphase a-SiOx systems (np-Si/a-SiO2). 5.2 Computational Approach 5.2.1 Valence Force Field Model for a-SiOx (0  x  2) Within the VFF model, the relative energies of a-SiOx materials are evaluated in terms of the increase of total energy (ΔEtotal) with respect to the Si-Si and Si-O bond energies obtained respectively from c-Si and c-SiO2 (β-cristobalite in this work). The ΔEtotal can be given by the sum of the changes in strain energy (ΔEstrain) and suboxide energy (ΔEsubox): ΔEtotal = ΔEstrain + ΔEsubox. (5.1) 104 The suboxide (penalty) energy (ΔEsubox) represents an increase in the Si-Si and Si-O bond energies arising from the various oxidation states of Si [188]. For a given Si-rich suboxide system, ΔEsubox can be obtained by adding the suboxide penalties of individual Si atoms with intermediate oxidation states (+1, +2, +3). Using periodic c-SiOx (x = 0.5, 1.0, and 1.5) models (see Figure 2 in Ref. 189), our DFT calculations predict the suboxide energies of 0.54, 0.57, and 0.29 eV for Si1+, Si2+, and Si3+, respectively, in good agreement with previous DFT results [188,189,190]. Strain energy (Estrain) arises from lattice distortions involving bond stretching, bond angle distortion, torsion resistance, and non-bonding interactions. The structure, stability, and phonon properties of bulk disordered Si and SiO2 materials have been successfully studied using a Keating-like VFF model:    i ji ijibstrain kbbkE , 2 0 2 0 )cos(cos2 1)( 2 1  , (5.2) where kb (in eV/Å2) and kθ (in eV) refer to the bond-stretching and angle-distortion force constants, respectively, bi and b0 (in Å) are the lengths of the ith bond and the equilibrium (reference) bond, respectively, and θij and θ0 (in degrees) are the angles subtended by the ith and jth bonds and the equilibrium bond angle, respectively. The three-body harmonic potential offers a satisfactory description of the strain of Si and SiO2 materials particularly when the departure of the bond lengths and bond angles from their respective equilibrium values is insignificant [66, 83,157]. For strain and suboxide energy variations, Keating-like (KT) potentials have been applied to examine the network topology and properties of a-Si and a-SiO2 [1,3,184]. In 105 particular, the KT potential parameterized by Tu and Tersoff [referred to as KT(TT), hereafter] has been widely employed to determine the atomic structure and energetics of amorphous Si/SiO2 multiphase systems, including planar c-Si/a-SiO2 interfaces [1,2] and np-Si/a-SiO2 [3,186,191]. Recently, Lee and Hwang optimized KT potential parameters based on the geometries and energies from density-functional theory (DFT) calculations [referred to as KT(LH) to distinguish it from KT(TT)]. This optimization procedure is detailed in the following section. Earlier studies emphasized the importance of kθ for the Si-O-Si bond angle to achieve a realistic description of bulk a-SiO2. The relatively small kθ(Si-O-Si) value in the KT(TT) potential appears to result in a large discrepancy between experimental and simulation results for the Si-O-Si bond angle distribution in the highly-distorted a-SiO2 network, which was corrected by making the kθ(Si-O-Si) term stronger [157]. Compared to bulk a-SiO2, we expect the SiO2 structure near the Si/SiO2 interface to be more distorted due to strain arising from the lattice mismatch between Si and SiO2. Therefore, an improved description of Si-O-Si bond angle distortion is likely warranted to obtain more realistic structural models and energetics for Si/SiO2 interface systems. Likewise, other angle distortion force constants (such as kθ(O-Si-O), kθ(Si-Si-O), and kθ(Si-Si-Si)) might also require reexamination since bond topologies and strain energies in a-SiOx materials are mainly governed by three-body contributions. In addition, the equilibrium Si-Si and Si-O bond lengths are known to be functions of Si charge state [172,188,190,192]. It is reasonable to infer that optimization of b0(Si-Si) and b0(Si-O) values in terms of Si oxidation state might be influential; however, we will later show 106 that equilibrium bond length variations are insignificant as well as their influence on the suboxide structure. 5.2.2 Determination of Force Field Parameters The potential parameters are determined by fitting VFF total energy data to DFT values in the following sequence that corresponds to an increase in the degrees of freedom of each training set: (1) b0(Si-Si) and b0(Si-O); (2) kb(Si-Si) and kb(Si-O); (3) θ0(Si-Si-Si); θ0(O-Si-O), θ0(Si-O-Si), and θ0(Si-Si-O); (4) kθ(Si-Si-Si); (5) kθ(Si-O-Si) [and also nθ(Si-O-Si), power of the three-body term]; (6) kθ(O-Si-O); and (7) kθ(Si-Si-O). Table 5.1 summarizes the force constant values for both the KT(LH) and KT(TT) potentials along with other KT potential parameters for bulk a-SiO2. Table 5.2 likewise summarizes calculated and tabulated b0 values. Table 5.1: Force constants for KT(LH) potential. Force constants for other KT potentials are also summarized for reference. 107 Table 5.2: Calculated equilibrium distances of Si-Si and Si-O bonds from bulk structures. The values from the cluster calculations are shown in (). For determination of b0(Sim-Sin) and b0(Sim-O), where m and n indicate the oxidation states of respective Si atoms, we used periodic crystalline Si0, Si1+, Si2+, Si3+, and Si4+ lattice models (see Figure 2 in Ref. 189) as well as cluster models for all oxidation states (see Figure 5.1). When m = n, the b0 values can be obtained from periodic calculations. Assuming that the variation of b0 with Si oxidation state is 108 identical for the periodic and cluster calculations, we tabulated the b0 values when m  n based on the cluster calculation results:     ),(),( ),(),(),(),(),(),(),( nnBmmB nnBnmBmmbnmbmmBnnbnmb   , (5.3) where b and B refer to periodic and cluster calculation results, respectively. The b0(Sim- O) value decreases with increasing m, possibly attributed to a reduction in covalent bond character, consistent with previous ab initio calculations [172,188,190,192]. We observe the contribution of oxidation state to insignificantly affect the resultant a-SiOx bond topology as corroborated by minor perturbations (substantially less than 0.1 Å) in average Sim-Sin and Sim-O bond lengths. Figure 5.1: Representative cluster models used for calculating Si-Si and Si-O bond lengths. Cluster models for other oxidation states (see Table 5.2) are obtained by adjusting the number of O atoms from these models. For calculating Si-O bond lengths, the two Si atoms neighboring the central O atom retain the same oxidation state. 109 Figure 5.2: Variations (ΔE) in total energies (from DFT, KT(LH), and KT(TT) calculations) per bond of (a) c-Si (with 8 atoms) and (b) c-SiO2 (β-cristobalite with 8 SiO2 units) as a function of the ratio (L/L0) of the strained lattice constant (L) to the equilibrium lattice constant (L0). 110 For kb(Si0-Si0) and kb(Si4+-O) parameters, we calculated variations in the total energies of c-Si (with 8 atoms) and c-SiO2 (β-cristobalite with 8 SiO2 units) by varying their respective lattice constants from -5% to 5%. For DFT calculations, a Monkhorst- Pack (8×8×8) k-point mesh was used for Brillouin-zone sampling. Optimized values of kb(Si0-Si0) = 9.08 eV/Å2 and kb(Si4+-O) = 31.90 eV/Å2 are close to corresponding KT(TT) values as shown in Table 5.1. For c-Si and c-SiO2, the variation of ΔE computed by KT(LH), KT(TT), and DFT calculations is shown in Figs. 5.2(a) and 5.2(b), respectively, as a function of the magnitude of bond strain. For both model systems, the VFF values are in good agreement with DFT values near equilibrium. The potential dependence of kb on Si oxidation state was also examined using cluster model calculations. From our results, kb appears inversely proportional to b0; however, the magnitude of kb variation is sufficiently small, so we can safely disregard the oxidation effect. With the equilibrium bond angle of θ0(Si-Si-Si) = 109.5°, we optimized the force constant kθ(Si-Si-Si) using four independent 64-atom a-Si supercells. The optimal value was obtained through minimization of the cross-validation error () which is given by:    N n n FF n DFT EEN 1 2)()(2 )(1 , (5.4) where )(nDFTE and )(n FFE refer to the DFT and FF energies, respectively, of the n th of N total a-Si models in the training set; in this case, the energies were evaluated based on fully- relaxed structures (with the same network) from each calculation. The same procedure was applied in optimization of other kθ values, unless stated otherwise. Our optimized kθ(Si-Si-Si) value of 1.795 eV is only half of the corresponding KT(TT) value of 3.58 eV. 111 For the remaining three-body force constants, equilibrium bond angles were set at θ0(O-Si-O) = 109.5°, θ0(Si-Si-O) = 109.5°, and θ0(Si-O-Si) = 180°, which are well- established for the Si-O system. Given that both Si-O-Si and O-Si-O bond angle distortions contribute to the energetics of a-SiO2, we first determined kθ(Si-O-Si) using a cluster model structure [see Figure 5.3 inset], then computed kθ(O-Si-O) using four independent, periodic a-SiO2 model structures (each containing 64-SiO2 units). Figure 5.3: Variations (ΔE) in total energies (from DFT, KT(LH), and KT(TT) calculations) of the cluster model (inset) as a function of Si-O-Si bond angle (θ). 112 From our DFT cluster calculations [see Figure 5.3], the total energy only slightly changes as θ(Si-O-Si) is reduced from 180° to 150°, but rapidly increases for θ(Si-O-Si) < 120°. The KT(TT) values of kθ(Si-O-Si) = 0.75 eV and nθ(Si-O-Si) = 2 (power of the corresponding three-body term) show reasonable agreement with DFT results for 150o < θ(Si-O-Si)  180 o, yet exhibit significant underestimation for θ(Si-O-Si) < 150o. Many previous studies [183,193-195] have demonstrated that the amorphous silica structure has a wide Si-O-Si angle distribution that may vary from 120° to 180°. Note that θ(Si-O-Si) can be around 130° and 160°, respectively, in three- and four-membered rings in a-SiO2. To more rigorously describe the strain energy variation associated with the wide distribution of the O-subtended bond angle, we adjusted not only kθ(Si-O-Si), but also nθ(Si-O-Si), in the three-body term,  nk )cos(cos 0 . By fitting VFF total energies to DFT values from the cluster and subsequent periodic calculations, we obtained kθ(Si-O- Si) = 2.62 eV and nθ(Si-O-Si) = 2.2. Note that these parameters are comparable to kθ(Si- O-Si) = 2.0 eV with nθ(Si-O-Si) = 2 from a careful reoptimization for a-SiO2 bulk phases by Alfthan et al. [157]. Concomitantly, we obtained kθ(O-Si-O) = 10.25 eV, which is also much larger than 4.32 eV from KT(TT), using four independent, periodic a-SiO2 model structures. We also adjusted kθ(Si-O-Si) and kθ(O-Si-O) simultaneously using four independent, periodic a-SiO2 model structures, but found that the optimized values were essentially unchanged. Finally, we determined kθ(Si-Si-O) = 4.165 eV using four independent, periodic a-SiO models comprised of 64 Si-O units, which is close to the KT(TT) value of 3.93 eV. 113 5.2.3 Metropolis Monte Carlo Simulations All a-SiOx structures we present were generated by CRN-MMC simulations combined with either KT(LH) or KT(TT) potentials. The atomic structure of each model system evolves toward thermodynamic equilibrium though MC bond-switching moves [65], which we implemented using the extended WWW (Wooten-Winer-Weaire) bond transposition scheme [66]. A bond switching move involves two bonds, A-B and C-D, across four unique atoms (A, B, C, and D) and forms two new bonds B-D and A-C by severing bonds A-B and C-D. A sampling process selects one of five different combinations of four distinct atoms (A, B, C, and D): Si(A)-Si(B)-Si(C)-Si(D); O(A)- Si(B)-O-Si(C)-O(D); Si(A)-Si(B)-Si(C)-O(D); Si(A)-Si(B)-O(C)-Si(D); and O(A)-Si(B)- Si(C)-O(D), where atoms A and C, as well as atoms B and D, must not be directly connected by a bond prior to the switching maneuver. For O(A)-Si(B)-O-Si(C)-O(D), the O atom between atom B and atom C is first selected randomly, then the remaining atoms are randomly identified. For the remaining combinations, either atom B or C is first selected randomly, then the remaining atoms are randomly identified. The acceptance or rejection of each bond-switching move is determined using probability P = min[1, exp(- ΔE/kBT)], where ΔE is the change in ΔEtotal resulting from the bond-switching move. Before and after each bond-switching move, the system is relaxed by Polak and Ribiere’s conjugate-gradient method [196]. During the MMC simulation, we included an additional repulsive term (Er) in ΔEtotal to effectively prevent nonbonded atoms from interacting [1,66]. Inclusion of Er is 114 particularly important in a-SiOx topological determination, likely because the flexible Si- O-Si linkages permit much more structural degrees of freedom than fourfold-coordinated a-Si. The repulsive contribution is given by:   mn mnr rdE 3 2 )( , (5.5) where m and n denote atoms which are neither 1st nor 2nd neighbors in the network, rmn is the distance between two atoms (evaluated only for rmn < d2), and d2 is a cutoff distance. We used the following parameters: d2(Si-Si) = 3.84 Å, d2(Si-O) = 3.2 Å, d2(O-O) = 2.61 Å, and γ = 0.5 eV/Å3, referring to Ref. 66 and 157. The Er term becomes negligible for the well-relaxed a-SiOx models presented in this paper. The following procedure was used for construction of each a-SiOx (0 ≤ x ≤ 2) structure model. First, we began with a randomized Si configuration in a periodic supercell with volume (V) given by V = VSi×NSi, where NSi denotes the number of Si atoms and VSi is the unit volume of a-Si. The randomized Si configuration was sequentially relaxed at temperatures of 5000, 4000, 3000, 2000, and 1000 K with approximately 1000NSi trials for each temperature. Next, NO (= xNSi) O atoms were randomly incorporated into Si-Si bonds in the a-Si model, resulting in an intermediate a- SiOx model with volume (V) given by V = VSi×(NSi – NO/2) + VSiO2×NO/2, where VSiO2 denotes the unit volume of a-SiO2. VSi and VSiO2 were extracted from corresponding experimental densities of 2.28 g/cm3 [167] and 2.2 g/cm3 [174], respectively. This intermediate configuration was further relaxed in a thermal sequence of 5000, 4000, 3000, 115 2000 and 1000 K with approximately 200(NSi + NO) trials for each temperature. Each time the simulation temperature was decremented, the lowest-energy configuration from the completed temperature step was selected as the initial configuration for the ensuing simulation step. 5.2.4 Density-functional-theory Calculations All DFT calculations herein were performed using the well-established planewave program, VASP [69], within the generalized gradient approximation of Perdew and Wang (GGA-PW91) [68]. Vanderbilt-type ultrasoft pseudopotentials [70] were adopted to describe the interaction between ion cores and valence electrons. Valence electron wave-functions were expanded using a planewave basis set with a kinetic-energy cut-off of 400 eV. For Brillouin zone sampling, we used a (222) Monkhorst-Pack k-point mesh for all periodic a-SiOx supercell models with 64 Si and 64x O atoms (sufficient for disordered systems) and -point sampling for cluster models, unless noted otherwise. All structures were fully-relaxed using the conjugate gradient method until residual forces on constituent atoms became smaller than 510-2 eV/Å. 5.3 Energetics of Amorphous SiOx: Force Field Models We evaluated the reliability of the force fields considered in this work for the energetics of a-SiOx materials by comparison with DFT results. Besides KT(LH) and KT(TT), we also looked at extended Stillinger-Weber (SW) potentials without and with 116 considering the suboxide penalty as proposed by Watanabe et al. [164,165]; for convenience, the former and latter are referred to as WT1 and WT2, hereafter. First, we prepared model structures for a-SiOx (x = 0, 0.5, 1, 1.5, and 2) using MC simulations based on the KT(LH) potential without including ΔEsubox to avoid suboxide phase separation into Si and SiO2. For x = 0.5, 1, and 1.5, the prevailing Si oxidation states are +1, +2, and +3, respectively, as listed in Table 5.3; in the model structures, O atoms are almost evenly distributed. For each x, we considered four independent structures each of which contains 64 Si atoms with 64x O atoms. Table 5.3: Si suboxide statistics averaged over four independent structures of a-SiO0.5 (64 Si and 32 O atoms), a-SiO1.0 (64 Si and 64 O atoms), and a-SiO1.5 (64 Si and 96 O atoms) used in Figure 5.4 and 5.5(a). These structures were constructed from CRN- MMC simulations based on the KT(LH) potential excluding suboxide penalty energies. 117 Figure 5.4: Relative total energies per Si atom (ΔÊtotal) from KT(LH), KT(TT), WT1, WT2, and DFT calculations for a-SiOx (x = 0, 0.5, 1.0, 1.5, and 2.0) (64 Si and 64x O atoms) structures constructed from CRN-MMC simulations based on the KT(LH) potential without suboxide penalty energies. For each x, four independent structures are represented. For x = 0.5, 1.0, and 1.5, the distributions of Si oxidation states are summarized in Table 5.3. Figure 5.4 shows the variations of ΔÊtotal (= ΔEtotal per Si atom) with x from DFT, KT(LH), KT(TT), WT1, and WT2 calculations. The DFT result (distribution) in Figure 5.4 resembles a parabola with maximum at x  1, driven mainly by suboxide penalty contribution (later demonstrated in Figure 5.5). Among the four classical potentials, KT(LH) exhibits the best agreement with DFT for all a-SiOx models. KT(TT) tends to 118 overestimate and underestimate the total energies of a-Si and a-SiO2, respectively. As expected, WT1 (with no suboxide penalty contribution) yields no significant variation in ΔÊtotal with x, while WT2 (whose pair-like interaction term was modified in order to describe the suboxide penalty [165,197]) significantly overestimates the suboxide contribution. In addition, compared to DFT, both WT1 and WT2 are likely to underestimate ΔÊtotal in a-Si, while showing relatively good agreement in a-SiO2. Figure 5.5(a) presents the average strain energies per Si (ΔÊstrain) from KT(LH), KT(TT), and DFT calculations, which were obtained by subtracting the average suboxide penalty energies (in the inset) from the average total energies (in Figure 5.4). For each x, average values represent four independent structures considered. Overall, KT(LH) and DFT are in good agreement. Compared to KT(LH) and DFT, KT(TT) yields a noticeably larger ΔÊstrain value in a-Si, where ΔÊstrain monotonically decreases with increasing O content and becomes smallest in a-SiO2. The ΔEstrain overestimation of KT(TT) for a-Si is mainly attributed to the larger kθ(Si-Si-Si) value of 3.58 eV relative to 1.795 eV in KT(LH), while the underestimated ΔEstrain in a-SiO2 is due to the smaller kθ(Si-O-Si) and kθ(O-Si-O) values of 0.75 eV and 4.32 eV relative to the respective values of 2.62 eV (with nθ = 2.2) and 10.25 eV in KT(LH). Note that the two-body force constants, kb(Si-Si) and kb(Si-O), are comparable for KT(TT) and KT(LH). 119 Figure 5.5: Average relative strain (ΔÊstrain) and suboxide (ΔÊsubox, inset) energies per Si based on KT(LH), KT(TT), and DFT calculations for a-SiOx (x = 0, 0.5, 1.0, 1.5, and 2.0) models constructed from CRN-MMC simulations based on the KT(LH) potential (a) excluding and (b) including suboxide penalty energies. For each x, four independent structures are represented. For x = 0.5, 1.0, and 1.5, the distributions of Si oxidation states for (a) and (b) are summarized in Tables 5.3 and 5.4, respectively. 120 We repeated the above procedure using model structures with partial phase separation which were obtained from MC simulations including ΔEsubox in KT(LH). For each x, four independent structures were considered. As summarized in Table 5.4, the Si oxidation state statistics clearly indicate the formation of Si and SiO2 phases in the suboxide systems. Figure 5.5(b) shows the variations of ΔÊstrain (and ΔÊsubox in the inset) as a function of Si:O ratio (x) from KT(LH), KT(TT), and DFT calculations. For the suboxide systems, KT(LH) shows excellent agreement with DFT, but the KT(TT) values significantly deviate from the DFT and KT(LH) values. Table 5.4: Si suboxide statistics averaged over four independent structures of a-SiO0.5 (64 Si and 32 O atoms), a-SiO1.0 (64 Si and 64 O atoms), and a-SiO1.5 (64 Si and 96 O atoms) used in Fig. 5.5(b). These structures were constructed from CRN-MMC simulations based on the KT(LH) potential including suboxide penalty energies. It is worth noting that ΔEstrain increases while ΔEsubox drops in the phase separation of suboxides, as seen from the separation-induced changes of ΔÊstrain and ΔÊsubox [Figure 5.5(a) vs. Figure 5.5(b)]. For instance, the phase separation results in an increase in 121 ΔÊstrain from 0.127 to 0.225 eV/Si when x = 1 (i.e., a-SiO), while ΔÊsubox reduces by 0.227 eV/Si. These results suggest that the role of strain might be important in determining the atomic configurations, particularly in the Si/SiO2 interface region, although the phase separation is mainly driven by the reduction of suboxide penalty energy [189]. 5.4 Phase separation: a-Si Cluster Embedded in a-SiO2 Matrix In this section, we examine how the atomic-level description of phase separation in a-SiO2 is affected by the choice of force fields. In particular, based on the KT(LH) and KT(TT) potentials, we attempt to assess the role of strain in determination of the atomic configuration near the Si/SiO2 interface. For both KT(LH) and KT(TT) potentials, we constructed five independent phase-separated model structures using CRN-MMC simulations. The structure generation procedure adopted the following steps for each model: (i) construction of a 480-atom a-Si supercell; (ii) insertion of 720 O atoms into Si- Si bonds from the supercell perimeter inward with concurrent volume compensation (following V = VSi×(NSi – NO/2) + VSiO2×NO/2 from Sec. 5.2.3 C); (iii) execution of O atom hopping moves at 100 K over 200(NSi + NO) trials to induce further phase separation (only ΔEsubox was considered, not ΔEstrain, to expedite phase separation); (iv) implementation of bond-switching moves within the oxide phase through a thermal sequence of 5000, 4000, 3000, 2000, and 1000K over approximately 200(NSi + NO) trials for each temperature; and (v) completion of bond-switching maneuvers throughout the supercell (both phases) in consecutive thermal stages of 3000, 2000, and 1000K over 122 approximately 200(NSi + NO) trials for each temperature. Each time the simulation temperature changed, the lowest-energy configuration from the prior simulation was adopted as the initial configuration for the subsequent simulation stage. This extensive approach provides a thorough description of phase separation in the a-SiO1.5 suboxide that leads to the formation of an a-Si cluster embedded in a a-SiO2 matrix. Example configurations from our simulations are presented in Figure 5.6. Figure 5.6: Atomic configurations for (a) KT(LH) and (b) KT(TT) models for the a-Si cluster embedded in a-SiO2 matrix (np-Si/a-SiO2). Grey wireframe represents O atoms and Si4+ states that comprise the a-SiO2 phase. Yellow, blue, red, and grey balls represent Si0, Si1+, Si2+, and Si3+ states, respectively. Comparing the phase-separated structures from KT(LH) and KT(TT) potential- based simulations (referred to as KT(LH) and KT(TT) models, hereafter), we find important discrepancies in the degree of phase separation (identifiable by the distribution of intermediate Si oxidation states) as well as the distribution of ring sizes. In Table 5.5, we summarize the relative concentrations of Si oxidation states for the KT(LH) and 123 KT(TT) models. While Si3+ is the dominant suboxide state in both models because of its low suboxide energy (0.29 eV) relative to those of Si1+ (0.54 eV) and Si2+ (0.57 eV), the overall concentration of suboxide states (Si1+, Si2+, Si3+) is higher in the KT(LH) model than the KT(TT) model. In addition, Table 5.5 shows that each suboxide state is more abundant in the KT(LH) model relative to the KT(TT) model. To provide some quantification of the abruptness of the phase transition interface regions, we calculated ratios of Si/SiO2 states (Si0, Si4+) to suboxide states (Si1+, Si2+, Si3+) as 0.78 and 0.81 for the KT(LH) and KT(TT) models, respectively. These results suggest that the KT(LH) model should yield more graded Si/SiO2 interface profiles with smaller a-Si cluster phases than KT(TT) models. Table 5.5: Si suboxide statistics averaged over five independent KT(LH) and KT(TT) models of np-Si/a-SiO2 (480 Si and 720 O atoms). 124 Figure 5.7: Profiles of (a) strain (ΔÊstrain) and (b) suboxide (ΔÊsubox) energies per Si along radial directions from cluster centers of KT(LH) and KT(TT) models for np-Si/a-SiO2 (480 Si and 720 O atoms). The cluster center is defined as the center of mass of Si0 atoms. The nominal interface radius, r0, is defined in the text. Each data point represents the average value within a given concentric spherical shell (2 Å thick) sampled over four independent structures. The two solid, horizontal lines depict the calculated strain energies for bulk a-Si and a-SiO2 with 216 Si and SiO2 units, respectively. All energies are calculated with the KT(LH) potential. 125 To further characterize the suboxide transition interface, we provide energy and Si suboxide distribution profiles in Figure 5.7 and 5.8, respectively, along radial directions from the a-Si cluster centers for both models. As shown in Figure 5.7, the a-SiO2 region in the KT(TT) model exhibits much higher ΔÊstrain values, but less Si suboxide penalty contributions, than the KT(LH) model; however, the KT(TT) model a-Si region exhibits lower ΔÊstrain values than observed in the KT(LH) model. We also observe that both a-Si and a-SiO2 regions in the proximity of the Si/SiO2 interfaces yield higher strain energies than bulk a-Si and a-SiO2. In Figure 5.8, the radial profiles of Si suboxide distribution clarify the inferences about phase transition abruptness extracted from the suboxide distribution results compiled in Table 5.5. For each model, we define a nominal interface radius, r0, that effectively defines a reference for the Si/SiO2 interface, )()( )()( 21 21 0      SinSin SirSir r , (5.6) where )( mSir is the distance of a Si atom with oxidation state m from the cluster center, )( mSin is the number of Sim, and the summations are conducted over all four independent samples studied. The Si1+ and Si2+ states can be interpreted as perimeter Si atoms of the a-Si phase with one and two O neighbors, respectively. For the Si1+ and Si3+ oxidation states, the increased spread in radial distribution of the KT(LH) model over the KT(TT) model is readily apparent in Figure 5.8. For both models, we observe a prominent peak in the Si3+ oxidation state just outside of r0. The KT(LH) Si3+ distribution also exhibits a more graded phase transition interface since the KT(LH) model has both a lower peak 126 and a more significant distribution tail on the a-SiO2 side (r – r0 > 5Å) when contrasted to the KT(TT) distribution. Figure 5.8: Profiles of (a) Si1+, (b) Si2+, and (c) Si3+ oxidation state distributions along radial directions from cluster centers of KT(LH) and KT(TT) models of np-Si/a-SiO2 (480 Si and 720 O atoms). The cluster center is defined as the center of mass of Si0 atoms. The nominal interface radius, r0, is defined in the text. Each data point represents the average value within a given concentric spherical shell (2 Å thick) sampled over four independent structures. 127 Figure 5.9: Ring-size distributions for the (a) entire, (b) a-Si, and (c) a-SiO2 regions of KT(LH) and KT(TT) models of np-Si/a-SiO2 (480 Si and 720 O atoms). 128 In Figure 5.9, we present ring-size distributions for the (a) total, (b) a-Si, and (c) a-SiO2 components of the KT(LH) and KT(TT) model structures. For the a-Si and a- SiO2 cases, the paths comprised solely of Si0 and Si4+ atoms are counted as rings, respectively; for comparison, the ring-size distributions of bulk a-Si and a-SiO2 are also provided. For a-Si, the embedded phase contains more five-membered rings in both models, rather than the energetically-favored six-membered rings that are most frequently observed in bulk a-Si. Likewise, the a-SiO2 phase in the two-phase system yields broader ring size distributions than in bulk a-SiO2 for both models. This indicates that the phase-separated Si and SiO2 structures are more strained than their bulk counterparts. We also notice in the a-Si phase that the KT(LH) model structures tend to contain more five-membered rings than the KT(TT) model structures; on the other hand, the latter generally exhibit broader ring size distributions than the former in the a-SiO2 phase. This is not surprising considering that the KT(LH) potential over- and underestimates lattice strain in a-SiO2 and a-Si, respectively, compared to the KT(TT) potential. 5.5 Mechanical Properties Our calculations suggest that the relative rigidity between Si and SiO2 matrices is critical in determination of the Si/SiO2 interface structure. Elastic (or Young’s) modulus (Y) and bulk modulus (B) are important metrics of the rigidity of an elastic response. For various a-SiOx compositions, these two moduli were successfully evaluated by first- principles calculations using a statistical approach in our previous work [19] and the 129 endpoint cases (x = 0 and 2) have been well-characterized through experimental measurements [198-205]. Additional mechanical properties, such as the Poisson ratio (ν) and shear modulus (G), can be calculated once Y and B are known because only two of these four quantities are independent in isotropic materials. [19] We apply our previously reported moduli calculation method to VFF total energy data to evaluate Y and B based on both KT(LH) and KT(TT) potentials for a-Si and a-SiO2 in order to quantify the degree of rigidity in respective a-Si and a-SiO2 matrices. Table 5.6: Computed average mechanical properties based on KT(LH) and KT(TT) potentials [1] for ten independent a-Si (216 Si atoms) and a-SiO2 (216 Si and 432 O atoms) structures. Strain was applied during mechanical property calculations using the same KT(LH) and KT(TT) potentials used during the initial CRN-MMC simulations. Relevant experimental data is also summarized for comparison. 130 Table 5.6 provides a summary of our mechanical property calculations along with relevant experimental data for comparison. The KT(LH) and KT(TT) potential-based calculations exhibit significant differences in Y and B values for both a-Si and a-SiO2, where the former provides better agreement with experimental data than the latter. The similar B values for a-Si from KT(LH) and KT(TT) calculations can be attributed to nearly identical kb values for the two potentials, suggesting that the bulk modulus is nearly unaffected by bond angle (Si-Si-Si) distortions. The remaining disparities for Y and B values in both a-Si and a-SiO2 between KT(LH) and KT(TT) potentials can be clearly explained by the aforementioned differences in kb and k. Considering the disparate nature of a-Si and a-SiO2 in the local proximity of a Si/SiO2 interface, the relative rigidity of SiO2 to Si should be an important factor in structural determination of these interfaces. Since both bond (Si-Si and Si-O) stretching and angle (Si-Si-Si, O-Si-O, and Si-O-Si) distortion contribute to Young’s modulus for a- Si and a-SiO2, we attempted to quantify the relative rigidity between a-Si and a-SiO2 using Y, rather than B. We evaluate the following dimensionless number, , as a measure of relative rigidity: Sia SiOa Y Y   2 , (5.7) 131 where Y(a-Si) and Y(a-SiO2) are the Young’s moduli for bulk a-Si and bulk a-SiO2, respectively. Our calculations show that a-Si (Y = 158.4 GPa) from KT(TT) is slightly more rigid than a-Si (Y = 124.3 GPa) from KT(LH), while a-SiO2 (Y = 46.3 GPa) from KT(TT) is far less rigid than a-SiO2 (Y = 98.9 GPa) from KT(LH). From these Y values, we obtain γKT(LH) = 0.8 and γKT(TT) = 0.3 for the KT(LH) and KT(TT) potentials, respectively. This indicates that the relative rigidity of SiO2 to Si is significantly underestimated by KT(TT). The smaller γKT(TT) implies that application of the KT(TT) potential will likely lead to structural rearrangement in the a-SiO2 phase driven by minimization of strain exerted on the a-Si phase, ultimately resulting in excess distortion in the a-SiO2 structure. In contrast, the larger γKT(LH) implies that a similar driving force for a-SiO2 structural distortion is significantly reduced for the KT(LH) potential. This provides a plausible explanation for the contrasting strain energy profiles of the KT(LH) and KT(TT) potentials as depicted in Figure 5.7(a). The occurrence of relatively more graded (abrupt) Si/SiO2 interfaces for the KT(LH) (KT(TT)) model structures can be explained by the difference in rigidity between a-Si and a-SiO2 phases. Phase separation of a-SiOx into Si and SiO2 phases is driven by minimization of the suboxide energy, but it concurrently creates additional distortion from lattice mismatch between Si and SiO2; as a result, the increase of strain energy from lattice mismatch tends to temper the formation of abrupt boundaries. In application of the KT(TT) potential, the excessively pliable a-SiO2 phase permits disproportionate lattice distortion on the a-SiO2 side of the interface which leads the system to form relatively abrupt Si/SiO2 interfaces. In contrast, in application of the 132 KT(LH) potential, the relatively more rigid a-SiO2 side of the interface is more resistive to accommodation of lattice distortion, so formation of relatively graded Si/SiO2 interfaces is favored (see Figure 5.8). 5.6 Summary We present a valence force field based on a modified Keating model for the structure and energetics of amorphous Si rich oxide (a-SiOx, 0  x  2) materials. The potential parameters for the strain energy contribution were optimized to fit DFT results for various cluster and periodic model structures. Suboxide energies were determined using DFT calculations of periodic c-SiOx (x = 0.5, 1.0, and 1.5) models, which are 0.54, 0.57, and 0.29 eV for Si1+, Si2+, and Si3+, respectively. We particularly focused on precise optimization of bond angle force constants such as kθ(Si-O-Si), kθ(O-Si-O), kθ(Si- Si-O), kθ(Si-Si-Si) since bond topologies and strain energies in a-SiOx are mainly governed by the three-body contributions. In this work, to more rigorously describe the strain energy variation associated with a wide Si-O-Si angle distribution (particularly in a highly strained Si/SiO2 composite system), we adjusted not only kθ(Si-O-Si), but also nθ(Si-O-Si), in the three-body term,  nk )cos(cos 0 . We also considered variations in the equilibrium bond lengths such as b0(Si-Si) and b0(Si-O) in terms of Si oxidation state, but the contribution of oxidation state turns out to insignificantly affect the resultant a- SiOx bond topology. For the energetics of various a-SiOx (0  x  2) systems, the present potential model agrees well with DFT for all O/Si composition ratios, while earlier 133 Keating-like and modified Stillinger-Weber potential models exhibit significant deviations from the present model and DFT. These results emphasize the importance of correctly describing the wide Si-O-Si angle distribution by making the corresponding bond-bending term stronger as well as softening of the Si lattice in the amorphous phase by making Si-Si-Si bond-bending term weaker. We also find that phase separation in a- SiOx results in an increase in the strain energy while the suboxide penalty energy reduces. Although the phase separation is mainly driven by the reduction of suboxide penalty energy, our calculations demonstrate that the role of strain is important in determining the atomic configurations particularly in the highly strained Si/SiO2 interface region. Our study also suggests that the relative rigidity between Si and SiO2 matrices is critical in determination of the Si/SiO2 interface structure. As such, as a measure of relative rigidity we introduced and evaluated a dimensionless number SiaSiOa YY  /2 , where Y(a-Si) and Y(a-SiO2) are the Young’s moduli for bulk a-Si and bulk a-SiO2, respectively. From the present potential model, the value of  is estimated to be 0.8 in the a-Si/a-SiO2 system, and decreases in the c-Si/a-SiO2 case. A smaller γ implies larger structural rearrangement in the SiO2 part driven by minimization of strain exerted on the Si part, ultimately resulting in more distortion in the a-SiO2 structure with a broader ring size distribution as well as less graded Si/SiO2 interface layer with a lower concentration of suboxide states (Si1+, Si2+, Si3+). The present potential model coupled with the CRN-MMC method can be used to create structural models (free of coordination defects) for complex a-SiOx- based materials, which will further allow thorough studies of the properties of these materials. 134 Chapter 6 Summary We present an effective computational approach to determine the structure and energetics of native defect clusters in Si, which combines continuous random network (CRN) model, tight binding molecular dynamics (TBMD), and density functional theory (DFT) calculations. This combined approach has been demonstrated to be an efficient technique for determining the lowest-energy configurations of defect clusters, particularly when they prefer to form fourfold-coordinated structures. Using the integrated approach, for the first time we have elucidated the complex growth and structure of self-interstitial and vacancy clusters in silicon and the effect of strain on the structure and stability of the defect clusters. The improved understanding regarding the structure and stability of native defects will greatly assist in better understanding the properties of native defects in Si, and explaining and predicting their impact on relevant materials properties and processing. We also present a valence force field based on a modified Keating model for the structure and energetics of amorphous Si rich oxide (a-SiOx, 0  x  2) materials. The potential parameters for the strain energy contribution were optimized to fit DFT results for various cluster and periodic model structures. Although the phase separation is mainly driven by the reduction of suboxide penalty energy, our calculations demonstrate that the role of strain is important in determining the atomic configurations particularly in the highly strained Si/SiO2 interface region. Our study also suggests that the relative rigidity between Si and SiO2 matrices is critical in determination of the Si/SiO2 interface 135 structure. The present potential model coupled with the CRN-MMC method can be used to create structural models (free of coordination defects) for complex a-SiOx-based materials, which will further allow thorough studies of the properties of these materials. 136 Appendices 137 Appendix A Summary of Publications A.1 Structure and stability of small compact self-interstitial clusters in crystalline silicon Sangheon Lee and Gyeong S. Hwang, Phys. Rev. B 77, 085210 (2008) We have determined the atomic structures and energies of small, compact self- interstitial clusters (In, n  10) in Si using a combination of continuous random network model, tight binding molecular dynamics, and density functional theory calculations. We present new local-minimum configurations for compact self-interstitial clusters with n = 6 through 10, together with well defined smaller clusters (n  5) for comparison. The cluster formation energies per interstitial exhibit strong minima at n= 4 and 8, in good agreement with earlier inverse model studies based on experimental observations. A.2 Growth and shape transition of small silicon self-interstitial clusters Sangheon Lee and Gyeong S. Hwang, Phys. Rev. B 78, 045204 (2008) The growth behavior of small self-interstitial clusters in crystalline Si is presented based on extensive combined Metropolis Monte Carlo, tight-binding molecular dynamics and density functional theory calculations. New stable structures for small interstitial clusters (In, 5  n  16) are determined, showing that the compact geometry appears favored when the cluster size is smaller than 10 atoms (n < 10). The fourfold-coordinated dodeca-interstitial (I12) structure with C2h symmetry is identified to serve as an effective nucleation center for larger extended defects. This work provides the first theoretical support for earlier experiments which suggest a shape transition from compact to elongated structures around n =10. 138 A.3 Integrated atomistic modeling of interstitial defect growth in silicon Sangheon Lee, Robert J. Bondi, and Gyeong S. Hwang, Mol. Simulat. 35, 867 (2009). A new theoretical approach which combines Metropolis Monte Carlo, tight- binding molecular dynamics, and density functional theory calculations is introduced as an efficient technique to determine the structure and stability of native defects in crystalline silicon. Based on this combined approach, the growth behavior of self- interstitial defects in crystalline Si is presented. New stable structures for small interstitial clusters (In, 5  n  16) are determined and show that the compact geometry appears favored when the cluster size is smaller than 10 atoms (n < 10). The fourfold- coordinated dodeca-interstitial (I12) structure with C2h symmetry is identified as an effective nucleation center for larger extended defects. This work provides the first theoretical support for earlier experiments which suggest a shape transition from compact to elongated structures around n = 10. We also provide some theoretical evidence that suggests that {311} extended defects grow slowly along <233> and relatively faster along <110>, which is consistent with typical defect aspect ratios observed through transmission electron microscopy. A.4 Biaxial strain effects on the structure and stability of self- interstitial clusters in silicon Robert J. Bondi, Sangheon Lee, and Gyeong S. Hwang, Phys. Rev. B 79, 104106 (2009). Using first-principles density-functional theory calculations, we examine variations in the structure and stability of small self-interstitial clusters (In, n ≤ 10) in crystalline silicon across a range of biaxial strain conditions (-3% ≤ ε ≤ 3%) on Si(100). Under the strain conditions considered, there is no significant deviation in the ground- state configuration of any cluster from the strain-free case. However, the relative stability of I4 and I8 is significantly increased under both compressive and tensile conditions, while other cluster sizes generally show less sensitivity to changes in strain. This suggests that I4 and I8 likely play an even larger role in the clustering/dissolution of interstitial defects in strained Si relative to strain-free Si. We find that the noteworthy strain dependence of I4 and I8 is attributed to the unique shape and symmetry of the I4- like core which allows reorientation within the lattice that is dependent on the compressive/tensile nature of biaxial strain. 139 A.5 Theoretical characterization of silicon self-interstitial clusters in uniform strain fields Robert J. Bondi, Sangheon Lee, and Gyeong S. Hwang, Phys. Rev. B 80, 125202 (2009) We use first-principles density-functional theory calculations to evaluate the orientation-dependent stability of small, neutral self-interstitial clusters (In, n ≤ 4) in crystalline Si across a range of uniform strain conditions (-4% ≤ ε ≤ 4%) in both uniaxial and biaxial strain fields on Si(100). Comprehending the behavior of these small clusters under strain is important in extending our understanding of the evolutionary cycle of interstitial defects during the ion implantation and annealing processes that occur during semiconductor manufacturing. Our calculation results suggest that strain of sufficient magnitude can contribute to significant ground-state structural distortion and even generation of different cluster configurations. Our study also indicates that the relative stability change per unit change in applied strain is greater in the biaxial case than the uniaxial case for interstitial clusters. We provided localized strain-distribution profiles and modification of bulk Si density of states to characterize the extent to which interstitial clusters modulate crystalline Si structure. A.6 Prediction of the formulation of stable periodic self-interstitial cluster chains [(I4)m, m = 1-4] in Si under biaxial strain Robert J. Bondi, Sangheon Lee, and Gyeong S. Hwang, Appl. Phys. Lett. 94, 1 (2009) Using density functional theory calculations, we examined the structure and stability of extendable self-interstitial cluster configurations (In, n = 12, 16) with four- atom periodicity in crystalline silicon under biaxial strain (-4% ≤ ε ≤ 4%) on Si(100). In the absence of strain, the ground state configurations of I12 and I16 share a common structure (I12-like) with C2h symmetry and a four-atom repeating unit; however, we identified an extended configuration based on I4 (D2d symmetry) cluster aggregates [(I4)m (m = 3,4)] along <110> that is more favorable under certain magnitude of strain. While both the I12-like and (I4)m configurations exhibit relative stabilities that are a function f both strain and orientation, the larger (I4)m orientation effect is the primary reason that these structures are preferred in both highly tensile and highly compressive environments. This suggests that I4 derivatives may participate in the growth transition of Si self- interstitial clusters in the compact-to-extended size regime (10 ≤ n ≤ 20) under stain. 140 A.7 Theoretical determination of stable fourfold-coordinated vacancy clusters in silicon Sangheon Lee and Gyeong S. Hwang, Phys. Rev. B 78, 125310 (2008) We have identified stable fourfold-coordinated vacancy clusters (Vn, 3  n  18) in Si using a combination of Metropolis Monte Carlo, tight binding molecular dynamics, and density functional theory calculations. Our calculations show that the small vacancy defects exclusively favor fourfold-coordination thermodynamically, rather than hexagonal ring-like structure formation which has widely been adapted to explain the behavior and properties of vacancy defects. Among those examined, the fourfold- coordinated V12 cluster with C2 symmetry is identified to be the most stable, yielding a formation energy of 1.16 eV per vacancy. The fourfold-coordinated V12 structure is about 4 eV more favorable than the conventional hexagonal ring structure. We also discuss how the relative stability between the fourfold-coordinated and hexagonal ring configurations will change as the cluster size increases to greater than a few tens of vacancies. A.8 Formation and structure of vacancy defects in silicon: Combined Metropolis Monte Carlo, tight-binding molecular dynamics, and density functional theory calculations Sangheon Lee, Robert J. Bondi, and Gyeong S. Hwang, Phys. Rev. B 80, 245209 (2009) We present the formation and structure of vacancy clusters (Vn, n  48) in crystalline Si based on combined Metropolis Monte Carlo, tight binding molecular dynamics, and density functional theory calculations. In this size regime, vacancy clusters are predicted to favor fourfold-coordination by nullifying dangling bonds created by Si lattice atom removal. Our results also highlight the identification of a stable high- symmetry V32 configuration that exhibits a complex, but ordered tetrahedral core/shell shape. When n > 25, fourfold-coordinated (FC) clusters commonly show the core/shell figure while smaller FC clusters (10 < n < 25) exhibit the trace of the high-symmetry V12 structure that exhibits four identical void-like structural units surrounding a tetragonal core. In addition, our study reveals that the formation of thermodynamically favored FC clusters can be kinetically facile. 141 A.9 Strain effects on the stability and structure of vacancy clusters in Si: A first-principles study Robert J. Bondi, Sangheon Lee, and Gyeong S. Hwang, Phys. Rev. B 81, 245206 (2010) Using first-principles density-functional theory calculations, we investigate the influence of both biaxial and uniaxial strain (-4% ≤ ε ≤ 4%) on the stability and structure of small, neutral vacancy clusters (Vn, n ≤12) on Si(100). A thorough understanding of vacancy clusters under stain is an important step toward elucidation of the evolutionary life cycle of native defects, especially during semiconductor manufacturing. Fourfold- coordinated (FC) structures are more favorable than “partial hexagonal ring” (PHR) structures in the size regime of our study under stain-free conditions; however, FC structures are also more rigid and consequently more sensitive to strain. Our calculation results indicate that PHR structures can be thermodynamically more favorable than FC structures in the presence of specific strain conditions. In addition, we identify orientation effects in which the cluster symmetry and its alignment within the strain field dictate cluster stability; in consequence, both configurations and orientations are essential factors in the identification of minimum-energy vacancy structures in strained Si. Furthermore, highlights of our simulation results suggest that minimum-energy cluster configurations formed under stain are often different than minimum-energy cluster configurations formed in the absence of strain. 142 A.10 On the origin of Si nanocrystal formation in a Si suboxide matrix Decai Yu, Sangheon Lee, and Gyeong S. Hwang, J. Appl. Phys. 102, 084309 (2007) We examined mechanisms underlying Si nanocrystal formation in Si-rich SiO2 using a combination of quantum mechanical and Monte Carlo (MC) simulations. We find that this process is mainly driven by suboxide penalty arising from incomplete O coordination, with a minor contribution of strain, and it is primarily controlled by O diffusion rather than excess Si diffusion and agglomeration. The overall behavior of Si cluster growth from our MC simulations based on these fundamental findings agrees well with experiments. Figure A.1: Series of snapshots from KMC simulations of phase separation in Si suboxide with the initial Si supersaturation of (a) 10%, (b) 20%, and (c) 30% at 1100 ºC. Here, only Si0 atoms are displayed. The box size is 8.1×8.1×8.1 nm3 [189]. 143 A.11 Strain-induced formation of surface defects in amorphous silica: A theoretical prediction Chin-Lung Kuo, Sangheon Lee, and Gyeong S. Hwang, Phys. Rev. Lett. 100, 076104 (2008) We present a prediction for the formation of surface defects in a thin amorphous silica layer during relaxation of externally imposed stresses and strains, based on extensive ab initio molecular dynamics calculations. Our calculations show that the application of a biaxial compressive stress leads to the creation of edge-sharing tetrahedron and/or silanone defects at the silica surface, which turns out to facilitate strain relief with irreversible structural changes in the silica layer. We also discuss a possible correlation between the predicted formation of surface defects and the observed enhanced surface reactivity of amorphous silica under compressive strain conditions. Figure A.2: Variation in the relative total energy (in eV) of a thin a-SiO2 slab during ab- initio MD relaxation at 1200 K as a function of annealing time (in ps). The strain induced formation of surface edge-sharing tetrahedron and silanone defects is also indicated, together with their configurations (in the insets) [17]. 144 A.12 Structure and dynamics of silicon-oxygen pairs and their role in silicon self-diffusion in amorphous silica Chin-Lung Kuo, Sangheon Lee, and Gyeong S. Hwang, J. Appl. Phys. 104, 045906 (2008) Based on gradient corrected periodic density functional theory calculations, we present the formation, structure, and diffusion of SiO pairs in a-SiO2. We find that a SiO pair preferentially undergoes transformation into an O vacancy through a twofold- coordinated Si atom. We determine the pathways for SiO pair → divalent Si → O vacancy transformation and divalent Si diffusion, along with O vacancy diffusion. Based on these results, we also discuss how the presence of SiO pairs can enhance Si self- diffusion in a-SiO2. Figure A.3: Atomic structure of a SiO pair in a-SiO2. The enclosed shows SiO pair insertion into the a-SiO2 matrix, resulting in the formation of a divalent Si atom. The small dark (red) and big gray (pink) balls represent O and Si atoms, respectively [18]. 145 A.13 First-principles study of the mechanical and optical properties of amorphous hydrogenated silicon and silicon-rich silicon oxide Robert J. Bondi, Sangheon Lee, and Gyeong S. Hwang, Phys. Rev. B 81, 195207 (2010) We use first-principles density-functional theory calculations to predict mechanical and optical property variation with composition for a-Si:H (at.% H = 0, 5.9, 11.1, and 15.8) and a-SiOx (x = 0, 0.5, 1.0, 1.5, and 2.0). A better understanding of the properties of a-Si:H and amorphous silicon oxide (a-SiOx) is technologically important, particularly for photovoltaic and optoelectronic device applications, respectively. However, relatively little reliable property information is available for these amorphous materials except for the well-studied endpoint cases of a-Si and a-SiO2. Our DFT calculations within the generalized gradient approximation predicts that addition of H to a-Si monotonically reduces the elastic modulus (Y) by 18 % and bulk modulus (B) by 16 % as H incorporation increases to 15.8 at.% in a-Si:H. Similarly, addition of O to a-Si monotonically reduces Y by 35 % and B by 38 % as x increases to 2.0 in a-SiOx. Our optical spectra for the complex dielectric function, (), exhibit intensity reduction in the E2 transition peak of Im[()] and reduction in the low-frequency dielectric constant {ߝ଴ ൌ limఠ՜଴ ܴ݁ሾߝሺ߱ሻሿ} as either H or O are added to a-Si while the a-SiOx spectra additionally resolve a vivid blueshift of both the fundamental absorption edge and E2 transition energy as O content increases. Considering the large variation in reported experimental measurements and the limited availability of previous computational results, our property predictions provide valuable insight into the mechanical and optical behavior of a-Si:H and a-SiOx materials. 146 Figure A.4: Imaginary (a) and real (b) parts of the complex dielectric function spectra for a-SiOx computed from DFT-GGA calculations. The spectra are annotated according to O stoichiometry (x) relative to Si. For each composition, only the lowest energy structure is represented in the spectra [19]. 147 Appendix B CRN-MMC Simulation B.1 Metropolis Sampling The most significant breakthrough in Monte Carlo modeling took place when Metropolis et al. [207] described an approach where ‘instead of choosing configurations randomly, then weighting them with exp(-E/kBT), we choose configurations with a probability exp(-E/kBT) and weight them evenly’ [208]. The Metropolis prescription dictates that we choose points with a Boltzmann- weighted probability. The typical approach is to begin with some ‘reasonable’ configuration q1. Then q1 is randomly perturbed to give a new configuration q2. The probability p of ‘accepting’ point q2 is ݌ ൌ ݉݅݊ ቂ1, ୣ୶୮ ሺିாమ/௞ಳ்ሻୣ୶୮ ሺିாభ/௞ಳ்ሻቃ, (B.1) Thus, if the energy (E2) of point q2 is not higher than that (E1) of q1, the point is always accepted. If the energy of the second point is higher than the first, p is compared to a random number z between 0 and 1, and the move is accepted if p ≥ z [208]. B.2 Bond-switching Method In our CRN-MMC simulations, the aforementioned random perturbation occurs in the form of a bond-switching event, as depicted in Figure B.1. Wooten et al. generated a CRN model structure of amorphous Si and Ge with periodic boundary conditions in 1985 [65] by employing a basic bond-switching method. Subsequently, Tu et al. extended the bond-switching method into amorphous SiO2 [66]. Figure B.1: Schematic illustration of bond switching event in which bonds AB and CD are broken and bonds AC and BD are newly formed. 148 Appendix C Additional Applications of CRN-MMC Simulation C.1 c-Si/a-SiO2 Interface Figure C.1: Planar Si/SiO2 interface configurations for (a) Si(100)/a-SiO2, (b) Si(110)/a- SiO2, and (c) Si(111)/a-SiO2. The small gray balls represent O atoms, and big gold (Si0), green (Si1+), brown (Si2+), red (Si3+), and white (Si4+) balls represent Si atoms with different oxidation states. C.2 Hydrogenated Amorphous Silicon (a-Si:H) Figure C.2: Amorphous configurations of a-Si:H (64 Si and 4 H atoms). Silicon and hydrogen atoms are represented by gold (large) and blue (small) spheres, respectively. 149 C.3 Oxidized Silicon Nanowire Figure C.3: Oxidized (a) Si<100>, (b) Si<110>, and (c) Si<111> nanowire. The small gray balls represent O atoms, and big gold (Si0), green (Si1+), brown (Si2+), red (Si3+), and white (Si4+) balls represent Si atoms with different oxidation states. 150 C.4 Oxidized Silicon Quantum Dot Figure C.4: Oxidized Si quantum dot and cross sections. The small gray balls represent O atoms, and big gold (Si0), green (Si1+), brown (Si2+), red (Si3+), and white (Si4+) balls represent Si atoms with different oxidation states. 151 C.5 Oxidized Embedded Silicon Nanoparticle Figure C.5: Oxide embedded (a) amorphous Si cluster and (b) Si nanocrystal. The small gray balls represent O atoms, and big gold (Si0), green (Si1+), brown (Si2+), red (Si3+), and white (Si4+) balls represent Si atoms with different oxidation states. 152 C.6 Fullerene C.6.1 Buckyball Figure C.6: Comparison of the strain energies (ΔEstrain) for Cn buckyball structures between our force field (indicated as This work) and DFT (indicated as DFT). All buckyball structures are from CRN-MMC simulations. 153 C.6.2 Topological Defects in Graphene and Carbon Nanotube Figure C.7: Topological defects: (a) V12 in graphene and (b) V4 carbon nanotube from CRN-MMC simulations. 154 Bibliography 1. Y. Tu and J. Tersoff, Phys. Rev. Lett. 84, 4393 (2000). 2. Y. Tu and J. Tersoff, Phys. Rev. Lett. 89, 086102 (2002). 3. G. Hadjisavvas and P.C. Kelires, Phys. Rev. Lett. 93, 226104 (2004). 4. T.A. Kirichenko, D. Yu, S.K. Banerjee, and G.S. Hwang, Phys. Rev. B 72, 035345 (2005). 5. S. Lee and G.S. Hwang, Phys. Rev. B 77, 085210 (2008). 6. S. Lee and G.S. Hwang, Phys. Rev. B 78, 045204 (2008). 7. S. Lee, R.J. Bondi, and G.S. Hwang, Mol. Simulat. 35, 867 (2009). 8. S. Lee and G.S. Hwang, Phys, Rev. B 78, 125310 (2008). 9. S. Lee, R.J. Bondi, and G.S. Hwang, Phys. Rev. B 80, 245209 (2010). 10. R.J. Bondi, S. Lee, and G.S. Hwang, Phys. Rev. B 79, 104106 (2009). 11. R.J. Bondi, S. Lee, and G.S. Hwang, Appl. Phys. Lett. 94, 264101 (2009). 12. R.J. Bondi, S. Lee, and G.S. Hwang, Phys. Rev. B 80, 125202 (2009). 13. R.J. Bondi, S. Lee, and G.S. Hwang, Phys. Rev. B 81, 245206 (2010). 14. S. Lee, R.J. Bondi, and G.S. Hwang, Phys. Rev. B 80, 245209 (2009). 15. D. Yu, G.S. Hwang, T.A. Kirichenko, and S.K. Banerjee, Phys. Rev. B 72, 205204 (2005). 16. C.-L. Kuo and G.S. Hwang, Phys. Rev. Lett. 97, 066101 (2006). 17. C.-L. Kuo, S. Lee, and G.S. Hwang, Phys. Rev. Lett. 100, 076104 (2008). 18. C.-L. Kuo, S. Lee, and G.S. Hwang, J. Appl. Phys. 104, 054906 (2008). 19. R.J. Bondi, S. Lee, and G.S. Hwang, 81, 195207 (2010). 20. S. Lee, R.J. Bondi, and G.S. Hwang, to be published. 21. S. Lee, R.J. Bondi, and G.S. Hwang, to be published. 22. R.J. Bondi, S. Lee, and G.S. Hwang, to be published. 23. G.D. Watkins, Phys. Rev. B 12, 5824 (1975). 155 24. E.J.H. Collart, K. Weemers, N.E. B. Cowern, J. Politiek, P.H.L. Bancken, J.G.M. van Berkum, and D.J. Gravesteijn, Nucl. Instrum. Methods Phys. Res. B 139, 98 (1998). 25. S. Takeda, Jpn. J. Appl. Phys., Part 2 30, L639 (1991). 26. A. Agarwal, T.E. Haynes, D.J. Eaglesham, H.-J. Gossmann, D.C. Jacobson, J.M. Poate, and Y.E. Erokhin, Appl. Phys. Lett. 70, 3332 (1997). 27. A. Claverie, B. Colombeau, B. de Mauduit, C. Bonafos, X. Hebras, C. Ben Assayag, and F. Cristiano, Appl. Phys. A: Mater. Sci. Process. 76, 1025 (2003). 28. J.L. Benton, S. Libertino, P. Kringhøj, D.J. Eaglesham, J.M. Poate, and S. Coffa, J. Appl. Phys. 82, 120 (1997). 29. D.C. Schmidt, B.G. Svensson, M. Seibt, C. Jagadish, and G. Davies, J. Appl. Phys. 88, 2309 (2000). 30. S. Libertino, S. Coffa, and J.L. Benton, Phys. Rev. B 63, 195206 (2001). 31. M. Nakamura and S. Murakami, J. Appl. Phys. 94, 3075 (2003). 32. P.K. Giri, Semicond. Sci. Tech. 20, 638 (2005). 33. D. Pierreux and A. Stesmans, Phys. Rev. B 71, 115204 (2005). 34. D.J. Eaglesham, P.A. Stolk, H.-J. Gossmann, and J.M. Poate, Appl. Phys. Lett. 65, 2305 (1994). 35. L.H. Zhang, K.S. Jones, P.H. Chi, and D.S. Simons, Appl. Phys. Lett. 67, 2025 (1995). 36. G.S. Hwang and W.A. Goddard, Phy. Rev. Lett. 89, 055901 (2002). 37. G.S. Hwang and W.A. Goddard, Appl. Phys. Lett. 83, 1047 (2003). 38. G.S. Hwang and W.A. Goddard, Appl. Phys. Lett. 83, 3501 (2003). 39. C.-L. Kuo, W. Luo, and P. Clancy, Mol. Simulat. 29, 577 (2003). 40. S. Solmi, M. Ferri, M. Bersani, D. Giubertoni, and V. Soncini, J. Appl. Phys. 94, 4950 (2003). 41. M.Y.L. Jung, C.T.M. Kwok, R.D. Braatz, and E.G. Seebauer, J. Appl. Phys. 97, 063520 (2005). 42. S.A. Harrison, T.F. Edgar and G.S. Hwang, Appl. Phys. Lett. 87, 231905 (2005). 156 43. S.A. Harrison, T.F. Edgar and G.S. Hwang, Phys. Rev. B 72, 195414 (2005). 44. S.A. Harrison, T.F. Edgar and G.S. Hwang, Phys. Rev. B 74, 195202 (2006). 45. S.A. Harrison, T.F. Edgar and G.S. Hwang, Electrochem. Solid St. 9, G354 (2006). 46. J. Zhu, T. D. dela Rubia, L.H. Yang, G. Mailhiot, and G.H. Gilmer, Phys. Rev. B 54, 4741 (1996). 47. P.B. Rasband, P. Clancy, and M.O. Thompson, J. Appl. Phys. 79, 8998 (1996). 48. N. Arai, S. Takeda, and M. Kohyama, Phys. Rev. Lett. 78, 4265 (1997). 49. J. Kim, F. Kirchhoff, W.G. Aulbur, J.W. Wilkins, F.S. Khan, and G. Kresse, Phys. Rev. Lett. 83, 1990 (1999). 50. B.J. Coomer, J.P. Goss, R. Jones, S. Öberg, and P.R. Briddon, Physica B 274, 505 (1999). 51. J. Kim, F. Kirchhoff, J.W. Wilkins, and F.S. Khan, Phys. Rev. Lett., 84, 503 (2000). 52. S.K. Estreicher, M. Gharaibeh, P.A. Fedders, and P. Ordejόn, Phys. Rev. Lett. 86, 1247 (2001). 53. M.P. Chichkine and M.M. de Souza, Phys. Rev. B 66, 045205 (2002). 54. D.A. Richie, J. Kim, S.A. Barr, K.R.A. Hazzard, R. Hennig, and J.W. Wilkins, Phys. Rev. Lett. 92, 045501 (2004). 55. G.M. Lopez and V. Fiorentini, Phys. Rev. B 69, 155206 (2004). 56. L.A. Marqués, L. Pelaz, P. Castrillo, and J. Barbolla, Phys. Rev. B 71, 085204 (2005). 57. M. Cogoni, B.P. Uberuaga, A.F. Voter, and L. Colombo, Phys. Rev. B 71, 121203(R) (2005). 58. M. Posselt, F. Gao, and D. Zwicker, Phys. Rev. B 71, 245202 (2005). 59. A. Carvalho, R. Jones, J. Coutinho, and P.R. Briddon, Phys. Rev. B 72, 155208 (2005). 60. S.A. Centoni, B. Sadigh, G.H. Gilmer, T.J. Lenosky, T.D. de la Rubia, and C.B. Musgrave, Phys. Rev. B 72, 195206 (2005). 157 61. Y.A. Du, R.G. Hennig, T.J. Lenosky, and J.W. Wilkins, Eur. Phys. J. B 57, 229 (2007). 62. N.E.B Cowern, G. Mannino, P.A. Stolk, F. Roozeboom, H.G.A. Huizing, J.G.M. van Berkum, F. Cristiano, A. Claverie, and M. Jaraíz, Phys. Rev. Lett. 82, 4460 (1999). 63. C.J. Ortiz, P. Pichler, T. Fühner, F. Cristiano, B. Colombeau, N.E.B. Cowern, A. Claverie, J. Appl. Phys. 96, 4866 (2004). 64. M. Fujii, K. Toshikiyo, Y. Takase, Y. Yamaguchi, and S. Hayashi, J. Appl. Phys. 94, 1990 (2003). 65. F. Wooten, K. Winer, and D. Weaire, Phys. Rev. Lett. 54, 1392 (1985). 66. Y. Tu, J. Tersoff, G. Grinstein, and D. Vanderbilt, Phys. Rev. Lett. 81, 4899 (1998). 67. T.J. Lenosky, J.D. Kress, I. Kwon, A.F. Voter, B. Edwards, D.F. Richards, S. Yang, and J.B. Adams, Phys. Rev. B 55, 1528 (1997). 68. J.P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 69. G. Kresse and J. Furthmuller, VASP the guide (Vienna University of Technology, Vienna, 2001). 70. D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). 71. S. Birner, J. Kim, D.A. Richie, J.W. Wilkins, A.F. Voter, and T. Lenosky, Solid State Commun. 120, 279 (2001). 72. L. Colombo, Physica B 273-274, 458 (1999). 73. M. Gharaibeh, S.K. Estreicher, and P.A. Fedders, Physica B 273-274, 532 (1999). 74. M.P. Chichkine, M.M. De Souza, and E.M. Sankara Narayanan, Phys. Rev. Lett. 88, 085501 (2002). 75. M. Kohyama and S. Takeda, Phys. Rev. B 46, 12305 (1992). 76. J. Kim, J.W. Wilkins, F.S. Khan, and A. Canning, Phys. Rev. B 55, 16186 (1997). 77. P. Alippi and L. Colombo, Phys. Rev. B 62, 1815 (2000). 158 78. J.P. Goss, T.A.G. Eberlein, R. Jones, N. Pinho, A.T. Blumenau, T. Frauenheim, P.R. Briddon, and S. Öberg, J. Phys.: Condens. Matter 14, 12843 (2002). 79. G.D. Watkins, J.W. Corbett, and R.M. Walker, J. Appl. Phys. 30, 1198 (1959). 80. V.C. Venezia, D.J. Eaglesham, T.E. Haynes, A. Agarwal, D.C. Jacobson, H.-J. Gossmann, and F.H. Baumann, Appl. Phys. Lett. 73, 2980 (1998). 81. J.S. Williams, M.J. Conway, B.C. Williams, and J. Wong-Leung, Appl. Phys. Lett. 78, 2867 (2001). 82. R.A. Brown, D. Maroudas, and T. Sinno, J. Cryst. Growth 137, 12 (1994). 83. V.C. Venezia, T.E. Haynes, A. Agarwal, L. Pelaz, H.-J. Gossmann, D.C. Jacobson, and D.J. Eaglesham, Appl. Phys. Lett. 74, 1299 (1999). 84. M.T. Zawadzki, W. Luo, and P. Clancy, Phys. Rev. B 63, 205205 (2001). 85. G.S. Hwang and W.A. Goddard III, Phys. Rev. B 65, 233205 (2002). 86. T.A. Frewen, S.S. Kapur, W. Haeckl, W. von Ammon, and T. Sinno, J. Cryst. Growth 279, 258 (2005). 87. D.A. Abdulmalik and P.G. Coleman, Phys. Rev. Lett. 100, 095503 (2008). 88. D.J. Chadi and K.J. Chang, Phys. Rev. B 38, 1523 (1988). 89. J.L. Hastings, S.K. Estreicher, and P.A. Fedders, Phys. Rev. B 56, 10215 (1997). 90. A. Bongiorno, L. Colombo, and T. Diaz de la Rubia, Europhys. Lett. 43, 695 (1998). 91. T.E.M. Staab, A. Sieck, M. Haugk, M.J. Puska, Th. Frauenheim, and H.S. Leipner, Phys. Rev. B 65, 115210 (2002). 92. M. Saito and A. Oshiyama, Phys. Rev. B 53, 7810 (1996). 93. A. La Magna , S. Coffa, and L. Colombo, Phys. Rev. Lett. 82, 1720 (1999). 94. G. Amarendra, R. Rajaraman, G.V. Rao, K.G.M. Nair, B. Viswanathan, R. Suzuki, T. Ohdaira, and T. Mikado, Phys. Rev. B 63, 224112 (2001). 95. S. Chakravarthi and S.T. Dunham, J. Appl. Phys. 89, 4758 (2001). 96. B.P. Haley, K.M. Beardmore, and N. Grønbech-Jensen, Phys. Rev. B 74, 045217 (2006). 159 97. S. Dannefaer, V. Avalos, and O. Andersen, Eur. Phys. J. Appl. Phys. 37, 213 (2007). 98. D.V. Makhov and L.J. Lewis, Phys. Rev. Lett. 92, 255504 (2004). 99. S.E. Thompson, G. Sun, Y.S. Choi, and T. Nishida, IEEE Trans. Electron Devices 53, 1010 (2006). 100. K. Derbyshire, Solid State Technol. 50, 38 (2007). 101. P. Bhattacharya, Semiconductor Optoelectronic Devices (Prentice Hall, Upper Saddle River, NJ, 1997). 102. J. Singh, Physics of Semiconductors and Their Heterostructures (McGraw-Hill, New York, 1993). 103. L. Lin, T. Kirichenko, B.R. Sahu, G.S. Hwang, and S.K. Banerjee, Phys. Rev. B 72, 205206 (2005). 104. C.G. Van de Walle and R.M. Martin, Phys. Rev. B 34, 5621 (1986). 105. M.E. Levinshtein, S.L. Rumyantsev, and M.S. Shur, Handbook Series on Semiconductor Parameters (World Scientific, London, 1996), Vol. 1, p. 29. 106. J. Fortner and J.S. Lannin, Phys. Rev. B 39, 5527 (1989). 107. S. Kugler, G. Molnar, G. Petö, E. Zsoldos, L. Rosta, A. Menelle, and R. Bellissent, Phys. Rev. B 40, 8030 (1989). 108. S. Cheylan and R.G. Elliman, Appl. Phys. Lett. 78, 1225 (2001); S. Cheylan and R.G. Elliman, Nucl. Instrum. Methods Phys. Res. B 148, 986 (1999). 109. K.S. Min, K.V. Shcheglov, C.M. Yang, H.A. Atwater, M.L. Brongersma, and A. Polman, Appl. Phys. Lett. 69, 2033 (1996). 110. S.P. Withrow, C.W. White, A. Meldrum, J.D. Budai, D.M. Hembree, Jr., and J.C. Barbour, J. Appl. Phys. 86, 1 (1999). 111. E. Neufeld, S. Wang, R. Actz, Ch. Buchal, R. Carius, C.W. White, and D.K. Thomas, Thin Solid Films 294, 238 (1997). 112. A.I. Belov, V.A. Belyakov, V.A. Burdov, A.N. Mikhailov, and D.I. Tetelbaum, Journal of Surface Investigation-X-ray Synchrotron and Neutron Techniques 3, 627 (2009). 160 113. M. Fujii, Y. Yamguchi, Y. Takase, K. Ninomiya, and S. Hayashi, Appl. Phys. Lett. 85, 1158 (2004). 114. D.I. Tetelbaum, O.N. Gorshkov, V.A. Burdov, S.A. Truchin, A.N. Mikhaylov, D.M. Gaponova, S.V. Morozov, and A.I. Kovalev, Phys. Solid State 46, 17 (2004). 115. F. Iacona, C. Bongiorno, C. Spinella, S. Boninelli, and F. Priolo, J. Appl. Phys. 95, 3723 (2004). 116. D. Barba, F. Martin, and G.G. Ross, Nanotechnology 19, 115707 (2008). 117. Y.Q. Wang, R. Smirani, and G.G. Ross, J. Cryst. Growth 294, 486 (2006). 118. G.A. Kachurin, S.G. Cherkova, D.V. Marin, A. Misiuk, Z.S. Yanovitskaya, J. Jedrzejewsky, and I. Balberg, Phys. Status Solidi A 206, 78 (2009). 119. X.J. Hao, A.P. Podhorodecki, Y.S. Shen, G. Zatryb, J. Misiewicz, and M.A. Green, Nanotechnology 20, 485703 (2009). 120. W.L. Wilson, P.F. Szajowski, and L.E. Brus, Science 262, 1242 (1993). 121. J. Linnros, N. Lalic, A. Galeckas, and V. Grivickas, J. Appl. Phys. 86, 6128 (1999). 122. Y. Kanemitsu, T. Ogawa, K. Shiraishi, and K. Takeda, Phys. Rev. B 48, 4883 (1993). 123. L. Pavesi, L. Dal Negro, C. Mazzoleni, G. Franzo, and F. Priolo, Nature (London) 408, 440 (2000). 124. P. Pellogrino. B. Garrido, C. Garcia, J. Arbiol, J.R. Morante, Mmelchiorri, N. Daldosso, L. Pavesi, E. Scheid, and G. Sarrabayrouse, J. Appl. Phys. 97, 074312 (2005). 125. M. Carrada, A. Wellner, V. Paillard, C. Bonafos, H. Coffin, and A. Claverie, Appl. Phys. Lett. 87, 251911 (2005). 126. S. Tiwari, F. Rana, H.I. Hanafi, A. Hartstein, E.F. Crabbe, and K. Chan, Appl. Phys. Lett. 68, 1377 (1996). 127. H.I. Hanafi, S. Tiwari, and I. Khan, IEEE Trans. Electron Devices 43, 1553 (1996). 161 128. J. von Borany, T. Gebel, K.-H. Stegemann, H.-J. Thees, and M. Wittmaack, Solid-State Electronics 46, 1729 (2002). 129. S. Tiwari, F. Rana, K. Chan, L. Shi, and H. Hanafi, Appl. Phys. Lett. 69, 1232 (1996). 130. R.J. Walters, P.G. Kik, J.D. Casperson, H.A. Atwater, R. Lindstedt, M. Giorgi, and G. Bourianoff, Appl. Phys. Lett. 85, 2622 (2004). 131. T.Z. Lu, J. Shen, B. Mereu, M. Alexe, R. Scholz, V. Talalaev, and M. Zacharias, Appl. Phys. A: Mater. Sci. Process. 80, 1631 (2005). 132. L.C. Kimerling, L. Dal Negro, S. Saini, Y. Yi, D. Ahn, S. Akiyama, D. Cannon, J. Liu, J.G. Sandland, D. Sparacin, J. Michel, K. Wada, and M.R. Watts, Silicon Photonics (Springer-Verlag, Berlin, 2004). 133. P.N. Keating, Phys. Rev. 145, 637 (1966). 134. A.M. Stoneham, V.T.B. Torres, P.M. Masri, and H.R. Schober, Philos. Mag. A 58, 93 (1988). 135. J.M. Murrell and J.A. Rodrigues-Ruiz, Mol. Phys. 71, 823 (1990). 136. A.R. Al-derzi, R.L. Johnston, J.N. Murrell, and J.A. Rodrigues-Riuz, Mol. Phys. 73, 265 (1991). 137. F.H. Stillinger and T.A. Weber, Phys. Rev. B 31, 5262 (1985). 138. R.L.C. Vink, G.T. Barkema, N. Mousseau, and W.F. van der Weg, J. Non-Cryst. Solids 282, 248 (2001). 139. J.A. Hauch, D. Holland, M.P. Marder, and H.L. Swinney, Phys. Rev. Lett. 82, 3823 (1999). 140. J. Tersoff, Phys. Rev. B 37, 6991 (1988). 141. J. Tersoff, Phys. Rev. Lett. 56, 632 (1986). 142. J. Tersoff, Phys. Rev. B 38, 9902 (1988). 143. B.W. Dodson, Phys. Rev. B 35, 2795 (1987). 144. K.E. Khor and S. Das Sarma, Phys. Rev. B 38, 3318 (1988). 145. K.E. Khor and S. Das Sarma, Phys. Rev. B 39, 1188 (1989). 146. K.E. Khor and S. Das Sarma, Phys. Rev. B 40, 1319 (1989). 162 147. D.W. Brenner, Phys. Rev. Lett. 63, 1022 (1989). Brenner shows how the Tersoff format is equivalent to the embedded atom method of 148 and 149. 148. M.I. Baskes, Phys. Rev. Lett. 59, 2666 (1987). 149. M.I. Baskes, J.S. Nelson, and A.F. Wright, Phys. Rev. B 40, 6085 (1989). 150. J. Wang and A. Rockett, Phys. Rev. B 43, 12571 (1991). 151. B. Bolding and H.C. Andersen, Phys. Rev. B 41, 10568 (1990). 152. R. Biswas and D.R. Hamann, Phys. Rev B 36 (1987) 6434. 153. R. Biswas and D.R. Hamann, Phys. Rev. Lett. 55, 2001 (1985). 154. E. Kaxiras and K.C. Pandey, Phys. Rev. B 38, R12736 (1988) 155. M.Z. Bazant, E. Kaxiras, and J.F. Justo, Phys. Rev. B 56, 8542 (1997). 156. J.F. Justo, M.Z. Bazant, E. Kaxiras, V.V. Bulatov, and S. Yip, Phys. Rev. B 58, 2539 (1998). 157. S. von Alfthan, A. Kuronen, and K. Kaski, Phys. Rev. B 68, 073203 (2003). 158. B.W.H. van Beest, G.J. Kramer, and R.A. van Santen, Phys. Rev. Lett. 64, 1955 (1990). 159. S. Tsuneyuki, M. Tsukada, and H. Aoki, Y. Matsui, Phys. Rev. Lett. 61, 869 (1988). 160. P. Tangney and S. Scandolo, J. Chem. Phys. 117, 8898 (2002). 161. E. Demiralp, T. Çağin, and W.A. Goddard III, Phys. Rev. Lett. 82, 1708 (1999). 162. P. Vashishta, R.K. Kaila, J.P. Rino, and I. Ebbsjø, Phys. Rev. B 41, 12197 (1990). 163. B.P. Feuston and S.H. Garofalini, J. Chem. Phys. 89, 5818 (1988). 164. T. Watanabe, H. Fujiwara, H. Noguchi, T. Hoshino, and I. Ohdomari, Jpn. J. Appl. Phys., Part 2 38, L366 (1999). 165. T. Watanabe, D. Yamasaki, K. Tatsumura, and I. Ohdomari, Appl. Surf. Sci. 234, 207 (2004). 166. A.C.T. van Duin, A. Strachan, S. Stewman, Q. Zhang, X. Xu, and W.A. Goddard, J. Phys. Chem. A 107, 3803 (2003). 163 167. K. Laaziri, S. Kycia, S. Roorda, M. Chicoine, J.L. Robertson, J. Wang, and S.C. Moss, Phys. Rev. B 60, 13520 (1999). 168. W.H. Zachariasen, J. Am. Chem. Soc. 54, 3841 (1932). 169. B.E. Warren, H. Krutter, and O. Morningstar, J. Am. Ceram. Soc. 19, 202 (1936). 170. E.H. Poindexter, P. Caplan, B. Deal, and R. Razouk, J. Appl. Phys. 52, 879 (1981). 171. R. Helms and E.H. Poindexter, Rep. Prog. Phys. 83, 2449 (1998). 172. H. Fukuda, M. Yasuda, T. Iwabuchi, S. Kaneko, T. Ueno, I. Ohdomari, J. Appl. Phys. 72, 1906 (1992). 173. I. Štich, R. Car, and M. Parrinello, Phys. Rev. B 44, 11092 (1991). 174. J. Sarnthein, A. Pasquarello, and R. Car, Phys. Rev. Lett. 74, 4682 (1995). 175. J. Sarnthein, A. Pasquarello, and R. Car, Phys. Rev. Lett. 52, 12690 (1995). 176. A. Pasquarello and R. Car, Phys. Rev. Lett. 80, 5145 (1998). 177. A. Pasquarello, M.S. Hybertsen, and R. Car, Nature (London) 396, 58 (1998). 178. S. Susman, K.J. Volin, D.L. Price, M. Grimsditch, J.P. Rino, R.K. Kalia, and V. Vashishta, Phys. Rev. B 43, 1194 (1991). 179. L.V. Woodcock, C.A. Angell, P. Cheeseman, J. Chem. Phys. 65, 1565 (1976). 180. K. Tatsumura, T. Watanabe, D. Yamasaki, T. Shimura, M. Umeno, and I. Ohdomari, Phys. Rev. B 69, 085212 (2004). 181. K. Tatsumura, T. Watanabe, D. Yamasaki, T. Shimura, M. Umeno, and I. Ohdomari, Jpn. J. Appl. Phys., Part 1 42, 7250 (2003). 182. K. Tatsumura, T. Watanabe, D. Yamasaki, T. Shimura, M. Umeno, and I. Ohdomari, Jpn. J. Appl. Phys., Part 1 43, 492 (2004). 183. V.M. Burlakov, G.A.D. Briggs, A.P. Sutton, and Y. Tsukahara, Phys. Rev. Lett. 86, 3052 (2001). 184. K.O. Ng and D. Vanderbilt, Phys. Rev. B 59, 10132 (1999). 185. L. Kong and L.J. Lewis, Phys. Rev. B 77, 085204 (2008). 164 186. G. Hadjisavvas, I.N. Remediakis, and P.C. Kelires, Phys. Rev. B 74, 165419 (2006). 187. F. Djurabekova and K. Nordlund, Phys. Rev. B 77, 1153254 (2008). 188. D.R. Hamann, Phys. Rev. B 61, 9899 (2000). 189. D. Yu, S. Lee, and G.S. Hwang, J. Appl. Phys. 102, 084309 (2007). 190. A. Bongiorno and A. Pasquarello, Phys. Rev. B 62, R16326 (2000). 191. G. Hadjisavvas and P.C. Kelires, Physica E 38, 99 (2007). 192. T. Hottori, Crit. Rev. Solid State Mater. Sci. 20, 399 (1995). 193. J.R.G. DaSilva, D.G. Pinatti, C.E. Anderson, and M.L. Rudee, Phil. Mag. 31, 713 (1975). 194. F. Mauri, A. Pasquarello, B.G. Pfrommer, Y.-G. Yoon, and S.G. Louie, Phys. Rev. B 62, R4786 (2000). 195. R. Mozzi and B. Warren, J. Appl. Crystallogr. 2, 164 (1969). 196. E. Polak, Computational Methods in Optimization, Academic Press, 1971. 197. J. Samela, K. Nordlund, V.N. Popok, E.E.B. Campbell, Phys. Rev. B 77, 075309 (2008). 198. M. Szabadi and P. Hess, Phys. Rev. B 58, 8941 (1998). 199. D.M. Follstaedt, J.A. Knapp, and S.M. Myers, J. Mater. Res. 19, 338 (2004). 200. T. Rouxel, J. Am. Ceram. Soc. 90, 3019 (2007). 201. H. Ni, X. Li, and H. Gao, Appl. Phys. Lett. 88, 043108 (2006). 202. C. Tsou, Y.S. Huang, and H.C. Chang, NSTI-Nanotech. 3, 339 (2005). 203. B. Bhushan, Handbook of Nanotechnology (Springer, Berlin, 2004), p. 773. 204. X. Li, B. Bhushan, K. Takashima, C.-W. Baek, and Y.-K. Kim, Ultramicroscopy 97, 481 (2003). 205. © Corning Incorporated, HPFS® Fused Silica KrF Grade (2003). 206. J. Dalla Torre, J.-L. Bocquet, Y. Limoge, J.-P. Crocombette, E. Adam, G. Martin, T. Baron, P. Rivallin, and P. Mur, J. Appl. Phys. 92, 1084 (2002). 207. N. Metropolis, A.E. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller, J. Chem. Phys. 21. 1087 (1953). 165 208. C.J. Cramer, Essentials of Computational Chemistry 2nd ed, Wiley. 209. P.A. Stolk, H.-J. Gossmann, D.J. Eaglesham, D.C. Jacobson, C.S. Rafferty, G.H. Gilmer, M. Jaraíz, J.M. Poate, H.S. Luftman, and T.E. Haynes, J. Appl. Phys. 81, 6031 (1997). 210. T.A. Kirichenko, S.K. Banerjee, and G.S. Hwang, Phys. Rev. B 70, 045321 (2004). 211. T.A. Kirichenko, S.K. Banerjee, and G.S. Hwang, Phys. Status Solidi. B 241, 2303 (2004). 166 Vita Sangheon Lee was born in Dalseong-gun, Gyeongsangbuk-do, Korea on 30 March, 1978, the second son of Jaewoong Lee and Oghee Kim. He attended Kyungwon high school in 1993-1995. He received a Bachelor of Science in Chemical Engineering from Seoul National University in August 2003. While in college, he completed his compulsory military service in 1999-2001. In 2004, he worked as a process engineer for Honam Petrochemical Corporation in Korea. In August 2004, he entered the graduate school at the University of Texas to pursue a PhD and began research under the guidance of Professor Gyeong S. Hwang in the Department of Chemical Engineering. Permanent Address: Bosung-town 105-103, Cheonnae-ri 117, Hwawon-up, Dalseong-gun, Daegu-city 711-830, Republic of Korea This dissertation was typed by the author.