Theoretical Study of Pressure-Induced Phase Transitions and Thermal Properties for Main-Group Oxides and Nitrides Except where reference is made to the work of others, the work described in this dissertation is my own or was done in collaboration with my advisory committee. This dissertation does not include proprietary or classi?ed information. Bin Xu Certi?cate of Approval: Minseo Park Associate Professor Physics Jianjun Dong, Chair Associate Professor Physics Yu Lin Professor Physics Satoshi Hinata Professor Physics George T. Flowers Dean Graduate School Theoretical Study of Pressure-Induced Phase Transitions and Thermal Properties for Main-Group Oxides and Nitrides Bin Xu A Dissertation Submitted to the Graduate Faculty of Auburn University in Partial Ful?llment of the Requirements for the Degree of Doctor of Philosophy Auburn, Alabama August 10, 2009 Theoretical Study of Pressure-Induced Phase Transitions and Thermal Properties for Main-Group Oxides and Nitrides Bin Xu Permission is granted to Auburn University to make copies of this dissertation at its discretion, upon the request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iii Vita Bin Xu, son of Qingshan Xu and Xiufen Cai, was born on April 21, 1980 in Shijiazhuang, Hebei, P. R. China. He attended University of Science and Technology of China (USTC) in 1998, and graduated with a Bachelor of Science in Physics in 2003. He entered Graduate School at Auburn University to pursue a Doctoral degree in Physics in September 2003. He was married to Shi Chen, daughter of Taixuan Chen and Junhua Wang, on April 25, 2007 in China. iv Dissertation Abstract Theoretical Study of Pressure-Induced Phase Transitions and Thermal Properties for Main-Group Oxides and Nitrides Bin Xu Doctor of Philosophy, August 10, 2009 (B.S., University of Science and Technology of China, 2003) 211 Typed Pages Directed by Jianjun Dong Main group nitrides and oxides are important solid compounds with applications in ?elds ranging from structural ceramics to catalysts and electronic materials. We have theoretically investigated the pressure-induced phase transitions and thermal properties for a series of oxides, nitrides and oxynitrides of IIIB, IVB and IIIB group, including Al2O3, AlN, Si3N4, Ga2O3, and Ga3O3N. In this dissertation, thermodynamic potentials at ?nite temperatures were calculated within the quasi-harmonic approximation (QHA). Structural optimization and total energy calculation of unit cell models were carried out based on ?rst-principles density functional theory (DFT) within the local density approximation (LDA). Vibrational spectra were calculated using the real-space supercell force-constant (SC-FC) method. For pressure- induced phase transitions, we have studied the equilibrium transition conditions, as well as the kinetic process, including the investigation of transition pathways, kinetic barriers and the softening-phonon induced displacive transitions. In our studies of equilibrium transition conditions, we calculated the T-P phase dia- grams for Al2O3, AlN, Si3N4, and Ga2O3. Our predicted transition pressures are 84 and 134 GPa at 300 K for corundum?Rh2O3(II)?postperovskite transitions in Al2O3, 9.9 GPa at v 0 K for wurtzite?rocksalt transition in AlN, 7.5/7.0 GPa at 300 K for ?/??? transitions in Si3N4 and 0.5 and 39 GPa at 300 K for ????Rh2O3(II) transitions in Ga2O3. In our studies of the pressure-induced reconstructive phase transition pathways, we ex- amined the corundum-to-Rh2O3(II) transition in Al2O3 and the wurtzite-to-rocksalt tran- sition in AlN. We showed that the rhombohedral corundum phase and the orthorhombic Rh2O3(II) phase are related by intermediate structures with monoclinic symmetry (P2/c). Using the proposed transition pathway, we calculated the kinetic barriers for the forward (C-to-R) and backward (R-to-C) transitions and further predict the meta-stabilities of the two phases. For AlN, di?erent transition pathways were previously proposed. We reinter- preted the bond-preserving paths with long-range patterns of the ?transition units?, and our calculated kinetic barriers indicate that the long-range pattern is less important. We also showed that the bond-breaking paths are not energetically favored. In addition, based on the pressure dependencies of the barrier heights we explained the discrepancy of tran- sition pressure between the room temperature observation and the calculated equilibrium result. In our studies of displacive transitions in Si3N4, although ? phase is dynamically stable at low pressure, two competing phonon-softening mechanisms are found under high pressure. If the ??? transition is bypassed due to kinetic reasons at lower temperatures, the ? phase is predicted to undergo a ?rst-order ??P3 transition at about 38.5 GPa. This predicted metastable high-pressure P3 phase is structurally related to ?-Si3N4. We have also calculated the thermodynamic and elastic properties of these systems, and selected results are presented. Our predictions are in good agreement with available experimental and other theoretical data. Furthermore, we studied the shifted Raman scattering and its correlation with the growth direction in Ga2O3 nanowires. And collaborated with an experimental study that synthesized spinel-structured gallium oxynitride from Ga2O3+GaN mixtures at high pres- sure and high temperature, we showed that the optimal synthesis pressure is predicted to be close to the ?-to-? transition pressure of Ga2O3. vi Acknowledgments Although only my name appears as the author of this dissertation, its completion is truly contributed from many great people. I would like to thank everyone and everything came along into my life, even without mentioning here. I express my sincere gratitude to my research advisor, Dr. Jianjun Dong, for his continuous encouragement, support and advice throughout my graduate years. Dr. Dong taught me not only the knowledge and techniques in this ?eld, but also how to question thoughts and express ideas. I enjoy the great bene?t from his instruction, especially the emphasis on details. His tolerance and being understanding is greatly appreciated. I would like to thank Dr. An-Ban Chen, who introduced me to Auburn University. He is always there to help with my research and career. And I appreciate the helpful suggestions from the committee members, Dr. Yu Lin, Dr. Satoshi Hinata, Dr. Minseo Park, and the outside reader, Dr. Orlando Acevedo. When I was working as a TA, Dr. Satoshi Hinata gave me great advice on how to teach as a ?rst-time lecturer. And I thank Dr. Robert Boivin for our discussion about recitations and labs. I am also thankful to our IT management professional Philip Forrest. He ensured the computational ability of the computer cluster which is important for my research. And I would like to acknowledge my graduate friend Xiaoli Tang, for the research discussion and the in?uence of her diligent attitude, Zengjun Chen for his knowledge on Latex. Beyond all the above gratitude, the support and care from my family is invaluable. I especially would like to thank my parents and my wife, Shi Chen, for all their understanding and care. Finally, I appreciate the ?nancial support from DOE (DE-FG02-03ER46060) and NSF (HRD-0317741) that funded part of the research in this dissertation. vii Style manual or journal used Bibliography follows American Physical Society (APS) style Computer software used The document preparation package TEX (speci?cally LATEX) together with the departmental style-?le auphd.sty. viii Table of Contents List of Figures xii List of Tables xviii 1 INTRODUCTION 1 1.1 Pressure-Induced Phase Transitions and Thermodynamic Properties of Main- Group Oxides and Nitrides . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Al2O3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2 AlN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.3 Si3N4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.4 Ga2O3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.1.5 Ga3O3N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2 First-Principles Theoretical Studies . . . . . . . . . . . . . . . . . . . . . . . 10 1.3 Outline of Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2 THEORY OF PHASE TRANSITIONS IN SOLIDS 15 2.1 Phases and Crystal Symmetries . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Equilibrium Thermodynamic Theory of Phase Transitions . . . . . . . . . . 19 2.2.1 Thermodynamic Stability . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.2 Phase Diagram and Classi?cation of Phase Transitions . . . . . . . . 21 2.3 Landau Theory for Phase Transitions of Second Kind . . . . . . . . . . . . 25 2.3.1 The Order Parameter and Landau Free Energy . . . . . . . . . . . . 25 2.3.2 Dynamic Lattice Stability and Soft Phonon Modes . . . . . . . . . . 27 2.4 Transition Paths of Reconstructive Phase Transitions . . . . . . . . . . . . . 28 3 FIRST-PRINCIPLES CALCULATIONS OF THERMODYNAMIC PROP- ERTIES OF CRYSTALS 33 3.1 First-Principles Total Energy Theory . . . . . . . . . . . . . . . . . . . . . . 33 3.1.1 Density Functional Theory (DFT) . . . . . . . . . . . . . . . . . . . 33 3.2 Statistical Theory of Bulk Crystals . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.1 Phonon Theory of Lattice Dynamics . . . . . . . . . . . . . . . . . . 36 3.2.2 Statistical Harmonic Approximation . . . . . . . . . . . . . . . . . . 41 3.2.3 First-Principles Phonon Calculations . . . . . . . . . . . . . . . . . . 42 3.3 Thermodynamic Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.1 Quasi-Harmonic Approximation . . . . . . . . . . . . . . . . . . . . 46 3.3.2 Equation of State (EOS) Models . . . . . . . . . . . . . . . . . . . . 50 3.3.3 Case Study: MgO . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 ix 4 PRESSURE-INDUCED PHASE TRANSITIONS IN ALUMINIUM OXIDE AND ALUMINIUM NITRIDE 59 4.1 Aluminium Oxide: Al2O3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.1.2 Crystal Structures, Total Energies and Vibrational Properties . . . . 62 4.1.3 T-P Phase Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.1.4 Transition Pathways in the ?-to-Rh2O3(II) Transition . . . . . . . . 72 4.1.5 Thermal Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.1.6 High Pressure Elasticity . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.1.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.2 Aluminum Nitride: AlN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.2.2 Total Energy and Equilibrium Phase Transition . . . . . . . . . . . . 94 4.2.3 Microscopic Mechanism for the Wurtzite-to-Rocksalt Transition . . . 95 4.2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5 STUDY OF HIGH-PRESSURE PHASE TRANSITIONS IN SILICON NI- TRIDE 102 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.2 Crystal Structures, Static Binding Energies, and Vibrational Spectra . . . . 106 5.3 Equilibrium Thermodynamic Stability and Phase Transitions . . . . . . . . 114 5.4 Phonon-Softening Induced Structural Instability in ?-Si3N4 at High Pressures 116 5.5 Room Temperature Metastable P3 Phase . . . . . . . . . . . . . . . . . . . 121 5.6 Thermodynamic Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 5.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 6 FIRST-PRINCIPLES STUDY OF GALLIUM OXIDE AND GALLIUM OXYNITRIDE 132 6.1 Gallium Oxide: Ga2O3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 6.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 6.1.2 Total Energy and Vibrational Properties . . . . . . . . . . . . . . . . 136 6.1.3 T-P Phase Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 6.1.4 Blue-Shifted Raman Scattering in Gallium Oxide Nanowires . . . . . 141 6.1.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 6.2 Gallium Oxynitride: Ga3O3N . . . . . . . . . . . . . . . . . . . . . . . . . . 145 6.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 6.2.2 Total Energy Calculations of GaN and Ga3O3N . . . . . . . . . . . . 147 6.2.3 Theoretical Study of the Synthesis, Structure, and Stability of Ga3O3N Spinel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 6.2.4 Electronic Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 6.2.5 Phonon Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 6.2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 7 CONCLUSIONS AND FUTURE WORK 159 Appendices 180 x A Total Energy Calculation with VASP 181 B Elasticity Calculation 182 C Beyond Harmonic Approximation: Perturbation Theory of Lattice An- harmonicity 188 xi List of Figures 1.1 Photo of recently constructed departmental computer cluster. . . . . . . . . 12 2.1 Phase diagram with two solid phases, one liquid phase and one gas phase. Two triple points and one critical point are present. . . . . . . . . . . . . . 23 2.2 Variation of the Landau free energy with positive and negative A coe?cients. 26 3.1 LDA (with US-PP) calculated (a) phonon dispersion relation, (b) vibrational density of states of ?-Al2O3 at zero pressure. High symmetry points are A:parenleftbig 0, 12,0parenrightbig, ?: (0,0,0) Z: parenleftbig12, 12, 12parenrightbig, and D: parenleftbig12,0, 12parenrightbig. Lines denote theoretical spectrum and discrete squares denote experimental data. . . . . . . . . . . . 47 3.2 Potential energy vs. interatomic distance curve. Dashed line denotes har- monic potential curve. Minimum potential energy is provided at r = a . . . 48 3.3 LDA calculated thermal free energy at several volume points. (a) At T = 0 K. (b) At T = 2000 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.4 LDA calculated Fvib?V data at T = 2000 K are ?tted to several EOS models. 55 3.5 Thermal expansion of MgO at zero pressure. Lines represent this work from LDA calculation. Discrete symbols are reported measurements. . . . . . . . 56 3.6 Thermal expansion of MgO at zero pressure. Solid line denotes our calcula- tion with Debye model. Discrete symbols are reported measurements. . . . 58 4.1 Crystal structures of four polymorphs of Al2O3: (a) Corundum, (b) Rh2O3- (II), (c) Pbnm perovskite and (d) post-perovskite. . . . . . . . . . . . . . . 63 4.2 LDA calculated phonon dispersion curves and vibrational density of state of ?-Al2O3 at zero pressure using (a) US-PP and (b) PAW. Discrete squares denote experimental data172. . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.3 Calculated Raman-active frequencies of ?-Al2O3 as a function of pressure. (a) US-PP, (b) PAW. Solid and dashed lines represent Eg and A1g modes from linear ?tting with data sets below 20 GPa. . . . . . . . . . . . . . . . . 67 xii 4.4 Calculated IR-active frequencies of ?-Al2O3 as a function of pressure. (a) TO modes with US-PP, (b) TO modes with PAW, (c) LO modes with US-PP, (d) LO modes with PAW. Solid and dashed curves denote Eu and A2u modes from 2nd order polynomial ?tting. . . . . . . . . . . . . . . . . . . . . . . . . 69 4.5 Calculated Raman-active frequencies of Rh2O3(II)-Al2O3 as a function of pressure. (a), (c), (e), (g) US-PP, (b), (d), (f), (h) PAW. Solid curves are from 2nd order polynomial ?tting. . . . . . . . . . . . . . . . . . . . . . . . . 71 4.6 Calculated Raman-active frequencies of pPV-Al2O3 as a function of pressure. (a) US-PP, (b) PAW. Solid, dashed, dotted and dash-dotted curves represent Ag, B1g, B2g and B3g modes from 2nd order polynomial ?tting. . . . . . . . 72 4.7 Static enthalpies of Al2O3 polymorphs as a function of pressure relative to that of corundum phase using (a) US-PP and (b) PAW. The static ? ?Rh2O3(II) and Rh2O3(II)?pPV transition pressure are present. The PV phase is metastable at all pressures considered among the four polymorphs. 73 4.8 Calculated T-P phase diagram of Al2O3. (a) US-PP, (b) PAW. Dashed lines denote possible phase boundaries. . . . . . . . . . . . . . . . . . . . . . . . . 74 4.9 Comparison of structures of (a) corundum and (b) Rh2O3(II)-Al2O3. The corundum phase is viewed from an angle in which it appears ?orthorhombic?- like. The di?erences between two polymorphs are highlighted. . . . . . . . . 75 4.10 Al-O bond lengths as a function of the transition parameter at 84.0 GPa. As the transition parameter changes from 0 to 1, corundum structure transforms smoothly to the Rh2O3(II) phase. There are 12 distinct bonds and only one bond breaks and reforms during the transformation. . . . . . . . . . . . . . 77 4.11 Volume per atom for the corundum?Rh2O3(II) transformation at several pressures. Same volume ranges are used to illustrate the pressure e?ect. . . 77 4.12 Enthalpies for the corundum?Rh2O3(II) transformation at several pressures. 78 4.13 Enthalpy barrier height for the corundum?Rh2O3(II) and Rh2O3(II)?corundum transformations as a function of pressure. . . . . . . 79 4.14 Calculated phonon dispersion curves of corundum-Al2O3 at 204 GPa. . . . 79 4.15 Calculated phonon dispersion curves of Rh2O3(II)-Al2O3 at 0 GPa. . . . . . 79 4.16 Comparison of the present theoretical calculation with measured thermal ex- pansion coe?cients of ?-Al2O3 as a function of temperature at zero pressure. Solid and dashed lines represent our LDA calculations using US-PP and PAW method, respectively. Discrete symbols are reported experimental data.16?21 81 xiii 4.17 Theoretical thermal expansion coe?cients of ?-Al2O3 as a function of tem- perature at several pressures from using (a) US-PP, (b) PAW. . . . . . . . . 82 4.18 Comparison of the present LDA predicted adiabatic bulk modulus as a func- tion of temperature of ?-Al2O3 with Voigt-Ruess-Hill data by Te?et23, by Goto et al.25 and the polycrystalline data by sound velocity measurements by Chung and Simmons24. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.19 Comparison of calculated isobaric heat capacity and entropy as a function of temperature of ?-Al2O3 with experimental data at zero pressure. Discrete open circles denote measured data from Furukawa22. . . . . . . . . . . . . . 83 4.20 LDA calculated thermal expansion coe?cients of Rh2O3(II)-Al2O3 at several pressures as a function of temperature up to 3000 K. (a) US-PP, (b) PAW . 84 4.21 LDA calculated thermal expansion coe?cients of pPV-Al2O3 at several pres- sures as a function of temperature up to 3000 K. (a) US-PP, (b) PAW . . . 84 4.22 Pressure dependence of the elastic constants of corundum. (a) US-PP, (b) PAW. Discrete symbols represent directly calculated data. Curves are from the 2nd-order polynomial ?tting. . . . . . . . . . . . . . . . . . . . . . . . . 85 4.23 Pressure dependence of the elastic constants of Rh2O3(II)-Al2O3. (a) US- PP, (b) PAW. Discrete symbols represent directly calculated data. Curves are from the 2nd-order polynomial ?tting. . . . . . . . . . . . . . . . . . . . 87 4.24 Pressure dependence of the elastic constants of pPV-Al2O3. (a) US-PP, (b) PAW. Discrete symbols (except crosses) represent directly calculated data from this work. Crosses denote GGA calculation from Stackhouse et al.36 at 136 GPa. Curves are from the 2nd-order polynomial ?tting. . . . . . . . . . 88 4.25 Phase transformation along TP1. (a) Unit cell of B4 structure, (b) Unit cell of B1 structure, (c), (d) Transformation pattern of ?transition units? in B4 and B1 phases. Dashed lines denote primitive unit cell. (Origin shifted onto the position of an atom for better illustration of the unit cell.) . . . . . . . 92 4.26 Phase transformation along TP2, TP3, TP4 and TP5. (a) Unit cell of B4 structure, (b) Unit cell of B1 structure, (c), (d) Transformation pattern of ?transition units? in B4 and B1 phases. (Origin shifted onto the position of an atom for better illustration of unit cell.) . . . . . . . . . . . . . . . . . . 93 4.27 Four-atom ?transition unit? from B4 to B1 structure. . . . . . . . . . . . . 94 4.28 Cohesive energy per atom as a function of volume for B4 and B1-AlN. En- thalpy variation of both phases as a function of pressure is also shown. . . . 96 xiv 4.29 Enthalpy as a function of transition parameter relative to B4 phase for ?ve TPs at six pressures. (a) at 0 GPa, (b) at 5 GPa, (c) at the predicted transition pressure 9.85 GPa, (d) at 15 GPa, (e) at 20 GPa and (f) at 30 GPa.100 4.30 Forward (B4-to-B1) and backward (B1-to-B4) activation barrier height (en- thalpy) as a function of pressure for ?ve TPs. . . . . . . . . . . . . . . . . . 101 5.1 Polymorphs of Si3N4 and synthesis conditions . . . . . . . . . . . . . . . . . 104 5.2 Crystal structures of (a),(b) ?-, (c),(d) ?-, and (e),(f),(g) ?-Si3N4. In the panel of ?- and ?-Si3N4, the ?rst graph illustrates the unit-cell model and the second graph is the 2?2?1 supercell model viewed in the direction of c axis. In the panel of ?-Si3N4, the ?rst graph shows the conventional cubic cell of the spinel structure and the following two graphs show the fourfold and sixfold coordinated Si units(SiN4 and SiN6) with tetrahedra and octahedra, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.3 Energy-volume curves for ?-, ?- and ?-Si3N4 in the scale of per atom. ? phase has an equilibrium energy of 3 meV lower than that of ? phase. E0 of ? phase is 93 meV higher than ? phase. . . . . . . . . . . . . . . . . . . . . 110 5.4 Phonon dispersion curves and vibrational density of states (VDOS) of (a) ?-, (b) ?-, and (c) ?-Si3N4 at zero pressure. . . . . . . . . . . . . . . . . . . . . 112 5.5 Calculated dispersion curves (scattered circles) of mode Gr?uneisen parameter of (a) ?-, (b) ?-, and (c) ?-Si3N4 at zero pressure. Red horizontal line is present to separate the positive and negative values. . . . . . . . . . . . . . 113 5.6 Gibbs free energy of ?-Si3N4 relative to that of ? phase as a function of temperature. Solid, dashed and dotted lines represent the pressure of 0, 5 and 10 GPa, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.7 T-P phase diagram of Si3N4. Solid curve denotes the phase boundary be- tween ?- and ?-Si3N4. Dashed curve denotes the phase boundary between ?- and ?-Si3N4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.8 (a) Raman, (b) IR and (c) silent mode frequencies as a function of pressure up to 60 GPa for ?-Si3N4. Experimental pressure dependence of Raman modes up to 30 GPa is also presented in discrete symbols as a comparison200. Solid squares denote measurements upon pressure increase and open squares denote measurements upon pressure decrease. Several low-frequency modes are found to decrease with increasing pressure. One Bu branch of silent modes is found dropping to zero at about 60 GPa. . . . . . . . . . . . . . . 117 xv 5.9 Phonon dispersion of ?-Si3N4 at a pressure of 48 GPa. Two competing soft phonon modes are found: one TA branch at M point and one optic branch at ? point. No LOTO splitting correction is added for the interests of low- frequency modes only. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.10 The square of vibration frequency (?2) as a function of pressure for two competing soft phonon branches: one TA branch at M point and one Bu branch at ? point. Solid squares and circles represent data from calculation. Solid and dashed lines are from a linear ?tting. . . . . . . . . . . . . . . . . 119 5.11 Ball-stick models of (a) P63/m, (b) P?6, (c) P21/m and (d) P3 structures viewed along the c axis. Balls in dark color represent N atoms and Si atoms are in light color. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.12 The total energy of P63/m (?), P?6, P21/m and P3 structures as a function of volume. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.13 Enthalpy landscape and its contour plot as a function of fxy and fz at the transition pressure of 38.5 GPa. . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.14 Enthalpy barrier (relative to ? phase) as a function of linearly interpreted transition parameter at 30 GPa (thinner) and the transition pressure of 38.5 GPa (thicker). Solid curves denote the ??P3??P3 path and the dashed curves denote direct ??P3 path. Horizontal axis is de?ned as qualitative structural similarity. The left end represents ? structure and the right end represent P3 structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.15 Temperature dependence of volume thermal expansion coe?cient of bulk ?- Si3N4 at zero pressure. Solid and dashed lines show present work with static energies ?tted to the 2nd-order BM-EOS and 3rd-order BM-EOS, respectively, and both thermal free energies ?tted to the 2nd-order BM-EOS. Discrete symbols denote experimental data83?86. . . . . . . . . . . . . . . . . . . . . 127 5.16 Temperature dependence of bulk Gr?uneisen parameter of ?-Si3N4 at zero pressure. Solid and dashed lines show present work with static energies ?tted to the 2nd-order BM-EOS and 3rd-order BM-EOS, respectively, and both thermal free energies are ?tted to the 2nd-order BM-EOS. Discrete symbols denote experimental data85. . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.17 Temperature dependence of volume thermal expansion coe?cient of bulk ?- Si3N4 at pressures of 0, 10, 20 and 30 GPa. . . . . . . . . . . . . . . . . . . 129 5.18 Temperature dependence of isobaric heat capacity of ?-Si3N4 at zero pressure. Solid and dashed lines show present work with static energies ?tted to the 2nd-order BM-EOS and 3rd-order BM-EOS, respectively, and both thermal free energies are ?tted to the 2nd-order BM-EOS. Discrete symbols denote experimental data85?88. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 xvi 5.19 Temperature dependence of volume thermal expansion coe?cient of ?-, ?- and ?-Si3N4 at zero pressure. Discrete symbols represent experimental data89,90. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 6.1 Crystal structure of monoclinic ?-Ga2O3 phase. Note that Ga3+ cations (light color) occupy both tetrahedral and octahedral interstices within the ccp lattice of O2? ions (dark color). . . . . . . . . . . . . . . . . . . . . . . 133 6.2 Enthalpy di?erences relative to ? phase for ?, ? and Rh2O3(II)-Ga2O3 as a function of pressure up to 60 GPa. The ? and ? curves cross at 0.5 GPa, the ? and Rh2O3(II) curves cross at 40.9 GPa. . . . . . . . . . . . . . . . . . . 137 6.3 T-P phase diagram for ?, ? and Rh2O3(II)-Ga2O3 polymorphs. . . . . . . . 140 6.4 LDA-calculated formation enthalpy ?H (i.e., Gibbs free energy of formation ?Gformation at zero temperature) as a function of pressure using (a) LDA methods and (b) the GGA approach. The calculated positive ?H suggests an endothermic formation. The solid plots are calculated assuming the oxyni- trides are synthesized from ?-Ga2O3 and wurtzite GaN, and the dashed plots are calculated assuming the oxynitrides are synthesized from ?-Ga2O3 and wurtzite GaN. The solid plot and the dashed plot cross at the pressure of the ?-to-? phase transition in Ga2O3 (predicted to be 0.5 and 6.6 GPa by LDA and GGA methods, respectively). . . . . . . . . . . . . . . . . . . . . . . . . 151 6.5 (a) Electronic band structure and (b) Electronic density of states for Ga3O3N calculated using ?rst-principles (DFT) methods within the LDA. . . . . . . 155 6.6 Photoluminescence spectrum of the Ga2.8N0.64O3.24 sample obtained via high-P,T synthesis in the multianvil experiments, using UV laser excitation (325 nm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 6.7 (a) Raman spectrum collected for the Ga2.8N0.64O3.24 sample at an excita- tion wavelength of 514.5 nm. The bold solid lines on the frequency scale below indicate the positions of the Raman bands for the analogous spinel form of ?-Ge3N4. (b) Phonon density of states (VDOS) calculated for the R?3m pseudocubic Ga3O3N phase, predicted as ?model I? in the enthalpy calculations (Table 6.7 ). The dashed lines indicate the Raman-active modes for that phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 xvii List of Tables 2.1 Point and space groups of Bravais lattices and crystal structures . . . . . . 17 4.1 Space groups, formula units per primitive unit cell Z and Wycko? sites of four polymorphs of Al2O3: Corundum, Rh2O3-(II), Pbnm perovskite and post-perovskite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.2 Third-order BM-EOS parameters for Al2O3 polymorphs . . . . . . . . . . . 65 4.3 The theoretical data of Raman frequencies at zero pressure and their pressure dependencies, and comparison with other reported experimental results on corundum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.4 IR-active frequencies of corundum at zero pressure . . . . . . . . . . . . . . 69 4.5 Elastic constants of corundum at zero pressure and their pressure dependencies 86 4.6 Comparison of elastic moduli between corundum and Rh2O3(II) phase. The corundum phase is treated as a monoclinic crystal with 20 atoms per unit cell. 88 4.7 Third-order BM-EOS parameters, zero pressure structural parameters for B4 and B1-AlN and B4-to-B1 transition pressure Pt. Parameters listed for B1 phase is for the 8-atom conventional cubic unit cell. . . . . . . . . . . . . . 95 4.8 Bond-preserving TPs for the B4-to-B1 phase transition. G denotes the space group of two end structures and G? denotes the common subgroup. Z rep- resents the number of formula units per unit cell along the TP. The trans- formation matrices are given with respect to the primitive unit cell lattice vectors and the fractional atomic coordinates are in terms of the setting of G?. u ? 3/8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.1 Space groups, formula units per primitive unit cell and Wycko? sites of hexag- onal ?-, ?-Si3N4 and cubic ?-Si3N4. . . . . . . . . . . . . . . . . . . . . . . 102 5.2 Summary of calculated and measured crystal parameters of ?-, ?- and ?- Si3N4. V0 is the equilibrium volume per atom, B is the bulk modulus and B? is the ?rst-order pressure derivative. Measurements were made at room temperature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 xviii 5.3 Summary of phase transition pressure and temperature for synthesizing ?-Si3N4115 6.1 Third-order BM-EOS parameters for Ga2O3 polymorphs . . . . . . . . . . . 136 6.2 Calculated and experimental zone-center Raman peak positions and Gr?uneisen parameter for the ?-Ga2O3 phase . . . . . . . . . . . . . . . . . . 138 6.3 Calculated and experimental zone-center Raman peak positions for the ?- Ga2O3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 6.4 Estimated internal strains . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 6.5 Raman mode frequencies and frequency shifts in ?-Ga2O3 nanowires with the [40?1] and [110] growth directions. Overall, excellent agreement between the observed and calculated shifts is seen for all mode frequencies except the one marked with an *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 6.6 Third-order Birch-Murnaghan EOS parameters, zero pressure structural pa- rameters for wurtzite GaN . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 6.7 Parameters of the third-order Birch-Murnaghan equation of states of three unit-cell-based atomic models of the spinel-structured Ga3O3N that have the lowest energy calculated within the LDA. . . . . . . . . . . . . . . . . . . . 150 6.8 Analysis of local coordination ordering schemes within the three unit-cell models of the spinel-structured Ga3O3N in terms of the ratio of various types of AX4 tetrahedral and AX6 octahedral units. . . . . . . . . . . . . . . . . . 153 A.1 Value and meaning of ISIF tag . . . . . . . . . . . . . . . . . . . . . . . . . 181 B.1 Independent components of Cij for all the crystal classes . . . . . . . . . . . 184 B.2 Correction term to Cij due to pressure . . . . . . . . . . . . . . . . . . . . . 187 xix Chapter 1 INTRODUCTION 1.1 Pressure-Induced Phase Transitions and Thermodynamic Properties of Main-Group Oxides and Nitrides Main group nitrides and oxides are important solid compounds with applications in ?elds ranging from structural ceramics to catalysts and electronic materials. In this disser- tation, I report our studies of ?ve material systems of these groups: Al2O3, AlN, Si3N4, Ga2O3, and Ga3O3N. 1.1.1 Al2O3 Aluminium oxide (Al2O3) is often referred to as alumina, sapphire or aloxite in the ce- ramic, mining and materials science communities. It is commonly used as an abrasive due to its hardness and as a refractory material due to its high melting point. The naturally- occurring crystalline form of Al2O3 at ambient condition is primarily corundum. Rubies (Cr+3 doped) and sapphires are gem-quality forms of corundum. In high-pressure experi- ments, ruby usually serves as a standard pressure gauge (ruby scale) in diamond anvil cell (DAC)1 and sapphire is used as window material in shock wave experiments2. Al2O3 is also one of the major constituents of the Earth?s lower mantle. High-pressure behaviors and thermal properties of alumina are important for both experimental research and better understanding of the interior of the earth. At ambient conditions, The crystalline form of Al2O3 is corundum (?-Al2O3, space group R?3c). This structure is known to exist over a wide range of pressure and tem- perature conditions. In the past twenty years theoretical3?7 and experimental8?11 studies showed that ?-Al2O3 transforms into the Rh2O3(II) structure (space group Pbcn) around 1 90 GPa, and further transforms into the post-perovskite structure (space group Cmcm) at ?130 GPa. Several calculations predicted a transition from Rh2O3(II) structure to the orthorhombic perovskite (PV) structure (space group Pbnm) at even higher pressures5,12. However, it has never been observed in experiments. Recent calculations showed that Pbnm perovskite is not thermodynamically favored with respect to the corundum, Rh2O3(II) and post-perovskite phases6,7. Here we report our LDA calculated equilibrium T-P phase di- agram with both ultra-soft pseudopotentials (US-PP) and PAW method. The results are compared with other reported calculations3?6,6,7,13. In the corundum?Rh2O3(II) transformation, in situ heating is found necessary. At room temperature, X-ray di?raction experiments showed that corundum phase is stable up to 175 GPa14,15, which implies the existence of a large kinetic barrier for the corundum to transform into the Rh2O3(II) phase. On the other hand, Lin et al. found that the high-pressure Rh2O3(II) phase can be seen as low as 85 GPa (Pt = 96 GPa) on decompres- sion after laser heating10, indicating the Rh2O3(II)?corundum transition is also sluggish. No experiment has quenched the Rh2O3(II) phase to ambient conditions. We proposed a transformation pathway for the corundum?Rh2O3(II) transition and evaluate the kinetic barrier based on the proposed pathway. We further predict the meta-stability of two phases. Thermodynamic properties, such as thermal expansion coe?cient (TEC), heat capac- ity CP, entropy and adiabatic bulk modulus have long been studied experimentally16?25 . Almost a decade ago Hama et al. calculated the thermal properties of corundum phase by extending the formalism of Thomsen and combining the results with the Vinet model and the Debye model for lattice vibrtations26. To date, no ?rst-principles studies have been reported to predict the thermodynamic properties of Al2O3. And, both experimental and theoretical data of thermal properties of the high-pressure phases are lacking. We thus present our calculated high-pressure TEC and bulk modulus for Rh2O3 (II) and pPV phases as a function of temperature up to 3000 K. Before 2004, experimental elastic constants of ?-Al2O3 were reported from as early as 1950s25,27?29. And many e?orts of ab initio calculations were made to predict the elastic 2 constants independently30,31. However, some calculations showed con?ict in the sign of C1432. In 2004, Gladden et al.33 reported their reexamination of the elasticity of ?-Al2O3 using resonance ultrasound spectroscopy (RUS) and con?rmed that C14 is positive rather than negative. Gladden?s conclusion was later con?rmed by both measurement with a di?erent technique34 and ?rst-principles calculations35. Although measurements of Cij for the high-pressure phases of Al2O3 are still lacking, Duan et al.31 and Stackhouse36 have calculated the elastic constants of Rh2O3(II) (from 75 GPa to 300 GPa) and post-perovskite phases (at 136 GPa), respectively. Here we report our predicted Cij of ?-, Rh2O3(II)- and pPV-Al2O3 with pressure dependencies, which will be compared with other available results. Elastic properties of both corundum and Rh2O3(II) phases have been studied, but a direct comparison between the two phases is not accessible due to di?erent crystal classes. Here we compare the elastic properties of corundum and Rh2O3(II) phase by treating both phases as the common-subgroup monoclinic lattices. 1.1.2 AlN AlN is an important semiconductor primarily due to its wide band-gap and thermal properties. The structural changes and properties at high pressures are of special interests. Experiments37?40 and calculations41?46 revealed that a pressure-induced ?rst-order phase transition from wurtzite (B4) structure to rocksalt (B1) structure happens for AlN. On the experimental side, without heating, the lowest pressure at which the rocksalt structure started to show up is 14 GPa39 and the B4-to-B1 transition was observed to complete at 20-31.4 GPa39,40. Xia et al. also found that the rocksalt phase is quenchable to ambient conditions39. However, theoretical calculations consistently predicted a transition pressure that is lower than the experimental values. Most of the recent ?rst-principles predicted static Pt is less than 10 GPa. The discrepancy implies the existence of a hysteresis for the forward and backward transitions which is caused by the activation barrier at the transition. To calculate the activation barrier, we ?rst examine the previously proposed transition paths. Using the computer program COMSUBS 47, Stokes et al. proposed ?ve TPs by 3 preserving the nearest neighbors along the transition pathway and limiting the size of unit cell to no more than four di?erent Wycko? positions48. He also found a common bilayer sliding mechanism along ?ve paths, which had been pointed out by Zahn et al.49 and Sowa50. Another work adopted the systematic approach is from Capillas et al.51, which was performed using the databases and tools provided by the Bilbao Crystallographic Server52. They proposed eight possible paths with di?erent orthorhombic and monoclinic symmetries by setting the maximum k-index equal to 4, strain tolerance Stol < 0.15 and maximum atomic displacement ?tol< 2 ?A. The intermediate structures along all the eight paths have eight atoms per unit cell. The di?erence between these two reports is that the nearest neighbours are not preserved along all the eight paths proposed by Capillas et al. DFT calculations from Shimojo et al. showed that the enthalpy barrier of transforma- tion is independent of the three paths (Cmc21, Pna21 and P21) for CdSe. Cai suggested, without calculation, that the B4-to-B1 transition is characterized by the transformation of the four-atom ?transition unit?, while the long-range pattern may be less important. On the other hand, using MD simulations, Zahn et al. pointed out that the favored paths have a tendency to avoid excess strains during the transformation. It would be interesting to investigate the proposed TPs with another material via ?rst- principles method. Previous studies53?55 suggested that the energetically favored TPs are bond preserving. In this paper, we studied all ?ve bond-preserving TPs proposed by Stokes et al. and one bond-breaking TP proposed by Capillas et al. from an energetic point of view for the B4-to-B1 transition in AlN. The correlation of the enthalpy barrier with di?erent TPs and strains will be discussed. We also relate the bilayer sliding mechanism to the long-range patterns of ?transition units?, and Cai?s hypothesis will be examined. 1.1.3 Si3N4 For Si3N4, it is widely used in cutting tools and anti-friction bearings due to its excel- lent mechanical properties, low mass density, and thermal stability56. It is also used as an insulator layer or as an etch mask because of its dielectric properties and a better di?usion 4 barrier against impurities in microelectronics57. For its technological importance, the me- chanical and thermal properties of silicon nitride at ambient pressure has been investigated extensively by both experiment and theory56,57. In contrast, its properties at high-pressure is less known. ? (P31c) and ? (P63/m) phases are the only two bulk polymorphs of Si3N4 known at ambient pressure. Both phases can be synthesized by nitriding pure silicon. In 1999, ?-Si3N4 (or c-Si3N4, Fd?3m) with the cubic spinel structure was synthesized at high pressure and high temperature58. Despite intensive research e?orts in searching the ?post-spinel? phases in Group-IVB nitrides, the spinel structured ? phase remains as the only experimentally identi?ed high-pressure phase. Phase transitions in Si3N4 have drawn extensive attention for more than a decade. The relative phase stability between ? and ? phases has been a topic of investigation for many years. Direct measurements of energetics of Si3N4 were reported by Liang et al.59. However, the di?erence in formation enthalpies between ?- and ?-Si3N4 was founded to be less than the intrinsic experimental uncertainty of ?22 kJ/mol (?32.6 meV/atom). Nevertheless, the ? phase is believed to be the ground state in Si3N4 because no ??? transition is ever observed. The stability condition for ? phase has been experimentally studied at temperatures of 1300??1800?C and pressures up to 60 GPa60?67. A solution-precipitation mechanism was proposed for the ??? transformation67. The observed liquid phase on the ?-Si3N4 surfaces was believed to lower the activation energy of atomic transportation. The stability of pristine ?-Si3N4 at high temperatures is ascribed to the extremely high value of the activation energy with clean surfaces. On the theory side, several studies con?rmed that the static bonding energy of ? phase is slightly higher than that of ? phase68?71. Wendel et al.70 and Kuwabara et al.71 carried out statistical QHA calculations, and they both found that the ? phase remains metastable in the temperature range from 0 to 2000 K at ambient pressure. Yet, pressure e?ects on the relative thermodynamic stability between ? and ? phases was not addressed in previous studies. Our study is to understand the relative thermodynamic stability at high pressures. 5 The spinel structured ?-Si3N4 was ?rst synthesized by Zerr et al. with laser-heated diamond anvil cell (LH-DAC)58. Later experiments showed that ?-Si3N4 can be obtained from both ? and ?-Si3N4 upon compression and simultaneous in-situ heating58,72?75. We predict the equilibrium phase boundaries for the ??? and ??? transitions. The transition pressures for the ??? transition is about 0.5 GPa lower than that of ??? transition. The ? phase is quenchable to the ambient condition, and it remains stable at temperatures ranging up to about 1670 K at ambient pressure76,77. When ?-Si3N4 ?decomposes? at ambient pressure upon heating, the samples may consist of both ? and ?-Si3N476. The in-situ heating to high temperature is found to be necessary to form the ?-Si3N4 at high pressures. At room temperature the ??? transition is, however, by-passed. Zerr found that ?-Si3N4 exists up to 34 GPa and it then transforms into a new phase (labeled as ?-phase) under further compression78. This phase transition was identi?ed by Raman spectroscopy and energy dispersive X-ray powder di?raction (EDXD). But the structure of ? phase was not fully determined. Zerr proposed three possible unit-cells based on the EDXD pattern: two tetragonal and one orthorhombic. He further suggested that the ?-Si3N4 should be considered as a metastable intermediate stage in the ??? transition. Kroll has proposed a metastable willemite-II-Si3N4 phase which is an intermediate between ? and ?-Si3N4 in both energetics and density79. However, the wII phase is unlikely to be the experimentally observed unknown phase at high pressure and room temperature. Because 1) the wII phase, which is structurally closely related to the spinel ?-Si3N4, has been shown to have a signi?cantly lower activation barrier for the ??wII transformation, comparing to that of ??? transformation79. Although the activation barrier of the ??wII transformation is unknown, it is more likely to be high enough to exclude the room temperature transition. 2) The calculated Raman frequencies of wII-Si3N4 could not match many strong peaks appeared in the measurements, e.g., two observed peaks at about 500 cm?1 and 550 cm?1 are absent in the calculation. A recent experimental work from McMillan et al. reproduced Zerr?s ?ndings on ?-Si3N4, but excluded the wII cubic structure80. 6 Meanwhile, ?-Ge3N4 is found to transform into the metastable polymorph ?-Ge3N4 with hexagonal P3 symmetry at room temperature by Soignard et al.81 Ab initio calculation from Dong et al. showed that a ??P?6?P3 transition sequence could occur in Ge3N4 at the pressure of about 20 GPa and 28 GPa82, which are of second-order that driven by soft phonons. If ?-Ge3N4 directly transforms into the P3 structure, the transition was predicted to be ?rst-order and Pt = ?23 GPa. Dong also pointed out that the ??P?6 transition is originated from a soft silent Bu mode. Room temperature experimental study by Soignard et al. con?rmed the direct ??P3 transition associated with a 5-7% volume reduction81. The Raman data they observed excludes the intermediate P?6 structure. Based on the density consideration, Soignard et al. suggested that the new polymorph is a ?post-phenacite? phase, in stead of ?post-spinel?. Comparison of the X-ray di?raction and Raman data between Ge3N4 and Si3N4 shows similarity which may suggest a P3 structure for ?-Si3N4. It is still unclear whether there are intrinsic di?erences between the HP-RT behaviors of Si3N4 and Ge3N4, or the experimental results may be interpreted di?erently. Our study is to theoretically investigate structural instabilities and possible metastable phase transitions in ?- and ?-Si3N4 at high pressures and room temperature. We found no sign of dynamical instability in the ? phase at high-pressure. On the other hand, we predicted a phonon- softening related ?rst-order phase transition at about 38.5 GPa in ? phase. At this pressure, the density of the proposed high-pressure phase is 4.16 g/cm3 which is larger than that of ?-Si3N4 (3.71 g/cm3), yet smaller than that of ?-Si3N4 (4.53 g/cm3, calculated). We further estimated the kinetic barrier heights for our proposed ??P3 transition, which is only 67.23 meV/atom at 38.5 GPa. Despite being of ?rst-order phase transition, the small barrier height suggests that the P3 phase is unlikely to be recovered below 38.5 GPa. We also performed a series of systematical calculations of thermodynamic properties of Si3N4, such as thermal expansion coe?cient (TEC), heat capacity and bulk Gr?uneisen parameter, and compared our results with available experimental data83?90 and some pre- vious calculations70,71,91,92. The overall good agreement with experiment validates the 7 adopted statistical quasi-harmonic approximation (QHA) and the Birch-Murnagahn equa- tion of states (EOS) models. Our results support the prediction from Kuwabara et al.71 on the negative TEC of ? and ? phases at temperatures below 100 K. We attributed the origin of the negative TEC to the low-frequency phonon modes with the negative mode Gr?uneisen ratios in the two phases. 1.1.4 Ga2O3 Monoclinic gallium oxide (Ga2O3) is usually known as a wide-band-gap semiconduc- tor (Eg = 4.9 eV); however, the conductivity can be varied from insulating to conducting behavior depending upon the preparation conditions93. It is well known that Ga2O3 can exist in several forms, including ?, ?, ?, ?, and ? polymorphs that all have di?erent struc- ture types94. Of these, the most stable form at ambient conditions is determined to be ?-Ga2O3. It is of great interest to determine the pressure-induced phase transformations among Ga2O3 polymorphs in order to establish the stable and metastable phase relations between di?erent crystalline modi?cations, and to evaluate their production under di?erent synthesis conditions. The relative densities of ?- and ?-Ga2O3 are 5.94 and 6.48 g?cm?3, respectively95, indicating that a ? ? ? transformation should occur at high pressure. Nanocrystalline ?-Ga2O3 particles embedded in a glassy matrix were also studied at high pressure using energy-dispersive x-ray di?raction96. In that work, a ?-to-? phase transfor- mation was found to be initiated at 6 GPa, but the process was not completed by 15 GPa, the highest pressure achieved in the study. However, it is known that the silica glass host matrix undergoes important structural and density changes within this pressure range97,98, so that it is not yet known if the structural changes are intrinsic to the ?-Ga2O3 material presumably in?uenced by the nanocrystalline nature of the sample , or are promoted by anomalous densi?cation among the SiO2 matrix. These results prompted us to theoretically investigate the high-pressure behavior of the phase-pure bulk ?-Ga2O3. One-dimensional nanostructured forms of ?-phase of gallium oxide (?-Ga2O3) such as nanotubes, nanobelts, and nanowires, have attracted recent interest due to enhanced 8 optical properties99,100. Recently, Choi et al.101 synthesized ?-Ga2O3 nanowires (diameter range of 15?45 nm) with a [001] growth direction using an arc-discharge method. Gao et al.102 synthesized [40?1] ?-Ga2O3 nanowires with diameters ranging from ?10?100 nm in a vertical radio-frequency furnace. Interestingly, the Raman mode frequencies of the [001] ?-Ga2O3 nanowires coincide with the corresponding frequencies in bulk ?-Ga2O3101. On the other hand102, the Raman mode frequencies of the [40?1] ?-Ga2O3 nanowires are redshifted relative to corresponding frequencies in bulk ?-Ga2O3 by 4?23 cm?1. Using plasma-enhanced chemical vapor deposition, Rao et al. have synthesized ?-Ga2O3 nanowires whose growth is along the [110] direction103, and the Raman spectrum is signi?cantly blueshifted in frequency104. Here we focus on the ?rst-principles calculations of the Raman mode frequencies under internal strains. Our calculated Raman frequency shifts suggest that the observed shifts in the nanowires with the [40?1] and [110] growth directions can be explained in term of di?erent internal strains, in contrast to the previously suggested quantum con?nement e?ects and defect-induced e?ects. 1.1.5 Ga3O3N The group 13 oxynitride materials have many useful properties related to their elec- tronic structure. ?-Ga2O3 with the corundum structure is conveniently alloyed with Al2O3 to provide selective reduction catalysts for gaseous NOx105, and various other Ga2O3 phases have been proposed as gas sensors, and in nanoscale structures as electron emitters and mag- netic memory materials106. Within the Al2O3?AlN system, several important AlxOyNz ceramic alloys and compounds are known. At high AlN contents, layered forms based on hexagonal/cubic intergrowths are present. As the Al2O3 content is increased, cubic spinel-structured materials begin to appear. A large family of defect spinels (?-Al2O3, AlxOyNz) contain vacancies on both cation and anion sites107. A stoichiometric oxynitride spinel-structured compound is obtained at the Al3O3N composition, in which Al3+ ions are present on the octahedral and tetrahedral sites, and O2? and N3? occupy tetrahedral anion sites108,109. Among the related nitride compounds Si3N4 and Ge3N4, high-pressure 9 synthesis has recently resulted in formation of a new class of spinel structures, that contain Si4+ and Ge4+ cations on both tetrahedral and octahedral sites110?114. Gallium oxyni- tride (Ga3O3N) has been predicted to form a new spinel-structured compound within the Ga2O3?GaN system, with potentially useful electronic properties115,116. It is predicted to be a direct wide bandgap semiconductor, comparable with GaN115. There has previously been an experimental report of a cubic gallium oxynitride phase with composition close to Ga2.8O3.5N0.5, that formed metastably during GaN thin ?lm synthesis from chemical precur- sors117,118. Here, we report our ?rst-principles theoretical study of the formation energetics, stability, and electronic properties of the Ga3O3N spinel-structured phase, combined with experiments using a combination of high pressure-high temperature techniques to establish the formation and stability of spinel-structured Ga3O3N from Ga2O3+GaN mixtures, and to determine the chemical composition, structure and properties of the resulting materials. 1.2 First-Principles Theoretical Studies Simulation of macroscopic properties of a physical system typically involves solving an ordinary or partial di?erential equation over large numbers of degrees of freedom. Even if a precise mathematical theory is available, it is only in very few cases that analytical solutions are possible. Computational physics is the study and implementation of numerical algorithm to solve problems in physics where a quantitative theory exists. Before the prevalence of powerful computers, empirical and semi-empirical approaches are often adopted, which rely on the phenomenological model or parameters ?tted from measured data. Nowadays, ab initio (or ?rst-principles) calculations are routinely performed in the ?elds of computational physics and chemistry. Compared with the empirical approach, predictions from ?rst- principles method can provide unbiased comparison and interpretation to the experimental data. It is also advantageous in calculating properties of materials at conditions where no or limited experimental data is available, and designing novel materials with promising properties. As one of the ab initio approaches, density function theory (DFT) is extremely successful in the electronic structure calculations for solids, which results from the work 10 of Hohenberg, Kohn and Sham119,120. Beside bulk systems, density functional theory has also become popular for complex materials such as nanostructures. In this dissertation, within the DFT, the many-electron exchange-correlation is approximated with the local density approximation (LDA). For part of the study that is associated with small energy di?erences, the LDA results are compared with calculations using the generalized gradient approximation (GGA). To improve numerical e?ciency, core electrons were approximated with ultrasoft pseudopotentials (US-PP)121. In the case of Al2O3, we also conducted a parallel comparative study using Projector Augmented-Wave (PAW) method. As mentioned above, the DFT based calculation is only possible and useful with the rapid development of high performance computer technologies, especially the parallel com- puting techniques. Considering the current computational e?ciency, the DFT calculation can only deal with atomistic models of no more than a few hundred atoms on a single- core CPU. Our calculations are performed on the departmental computational resource, a distributed memory computer system (Beowulf cluster). This cluster is comprised of 96 AMD Athlon MP CPUs running Red Hat Enterprise Linux 5. A photo of our recently con- structed cluster system is shown in Figure 1.1. We have implemented parallel algorithms to distribute each job into a set of calculations. Each corresponds to one DFT calcula- tion. The most computationally expensive part in this study is to calculate the real-space force-constant matrix using supercell models of ?100-200 atoms. I have been focused on adopting and further developing computational methods based on density functional theory and statistical theory to study the behavior of materials, in particular the pressure induced solid-solid phase transitions and thermodynamic properties, over a wide range of temperature-pressure (T-P) conditions. To predict the equilibrium transition pressure at a ?nite temperature and thermodynamic properties (such as thermal expansion coe?cient (TEC), heat capacity, entropy and bulk modulus, etc.), the calculation of free energy is required. We adopt quasiharmonic approximation (QHA) to model the lattice dynamics. Our approach of calculating phonon frequencies belongs to a method of the direct approach: a ?rst-principles real-space supercell force-constant (SC-FC) method 11 Figure 1.1: Photo of recently constructed departmental computer cluster. which calculates the phonon frequencies from the forces obtained via the Hellmann-Feynman theorem. This technique is proved to be e?cient and successful in predicting the full phonon spectra of many materials122?129. For pressure-induced phase transitions, we are interested in investigating both equi- librium transition conditions and the kinetic process. Regarding reconstructive transitions with no group-subgroup relation, the equilibrium phase boundary is determined by equat- ing the Gibbs free energy of the two phases. However, in general, there exists a kinetic barrier at the transition, which can be overcome by the thermal activation energy from the environment (small barrier) or heating process (large barrier). Based on the pressure dependencies of the forward and backward barrier heights, we can predict the metastability of polymorphs. The kinetic barrier can be estimated from the knowledge of the micro- scopic mechanism of the transition. For transformations with a group-subgroup relation, a transition path (TP) can be easily de?ned by a set of continuous atomic displacements 12 and/or lattice strains that the system transforms from one phase into another. For some simple reconstructive transformations with no group-subgroup relation, although nucleation processes may occur, di?usionless collective atomic displacements can still characterize the transformation on a local basis. The concept of transition pathway is then possible to describe the transformation from the starting phase to the ending phase in a continuous manner. Among the in?nite number of ways to transform one structure into another, the theoretical studies are restricted only to those most possible paths, e.g., preservation of bonds, less strains, etc. The pressure dependence of phonon spectrum can also reveal the information of struc- tural stability under compression. In the case when soft phonon happens, a vanishing phonon frequency indicates the disappearance of the restoring forces related to the corre- sponding normal mode, and the structure consequently undergoes a continuous transition to a lower-symmetry phase, which can be found according to the vibrational pattern based on the associated eigenvectors. However, the transition happened at room or higher tem- peratures is often ?rst-order, where no phonon has become soft yet. But this ?rst-order phase transition is driven by the phonon mode with softening tendency. 1.3 Outline of Dissertation The rest of this dissertation is organized as follows. In Chapter Two I ?rst review the fundamental theory of solid-solid phase transitions, which include the reconstructive phase transitions and the soft-phonon driven continuous transitions. Then our computational ap- proaches of investigating the transformation mechanism and phase metastability/instability are introduced. In Chapter Three I discuss the ?rst-principles methodologies we adopted in this study, i.e., total energy calculation based on density functional theory and our methods to calculate the ?nite-temperature thermodynamic potentials. The in?uence of equation of state models on the prediction of thermal properties is also presented. In Chapter Four we show our ?rst-principles studies on the pressure-induced phase transitions in Al2O3 and AlN. Our focus is on the microscopic mechanism of the transition from an energetic point 13 of view. The thermal and elastic properties of Al2O3 are also studied. In Chapter Five we predict a metastable high-pressure transition for ?-Si3N4. We also predict the thermody- namic properties of three known polymorphs of Si3N4. In Chapter Six I present our study of pressure-induced phase transitions in bulk Ga2O3 and the blue-shifted Raman frequencies in Ga2O3 nanowires, as well as the theoretical optimal synthesis condition, electronic prop- erties and phonon spectrum of the spinel-structured gallium oxynitride (Ga3O3N). Finally in Chapter Seven we summarize our key results and discuss the future works. 14 Chapter 2 THEORY OF PHASE TRANSITIONS IN SOLIDS 2.1 Phases and Crystal Symmetries Before we discuss the phase transitions, it is necessary to elucidate the concept of phase. The word ?phase? may have di?erent meanings in di?erent disciplines. In thermodynamics, a phase of a macroscopic system refers to a type of internal structure of the system when the composition of the system is speci?ed. Throughout a single phase, all the physical properties of the system is uniform at the macroscopic level. These properties include chemical composition, density, heat capacity, index of refraction and so on. For a material system, the internal structure is simply its atomic con?guration if magnetic or electronic degrees of freedom can be ignored. A macroscopic solid material is composed of numerous atoms or molecules, usually in the order of ? 1023. The atoms or molecules that compose the solid are closely packed together and the chemical bonds between them are relatively strong. The short-range structural orders are largely controlled by the chemistry of the constituent elements, while statistics a?ects the long-range structural order. According to their long-range atomic order/disorder, solids are divided into two categories: crystals and amorphous solids. Real solids contain imperfections, such as surfaces, grain boundaries, and defects. In this dissertation, my study is limited to only ideal crystals, which are good approximations to real bulk crystalline materials that contain only small amount of defects and impurities. Symmetry is an essential character of crystals. Atoms in a crystal vibrate around their equilibrium positions, which form a periodical and ordered pattern in the three dimension space. The pattern is made up of a group of equilibrium positions of atoms, which is called the atomic basis. The collection of repetition of identical structural units with translational 15 symmetry is called Bravais lattice. crystal structure = Bravais lattice + atomic basis (2.1) The location of each Bravais lattice point, in another word, the origin of each repeated unit in a crystal, can be described by R = n1a1 + n2a2 + n3a3 (2.2) where the 3-D vectors ai (i = 1,2,3) are linearly independent unit-cell vectors. n1, n2 and n3 are integers. The smallest repeatable unit cell of a lattice is called primitive unit cell, and the corresponding unit vectors are primitive unit-cell vectors. Note that the choice of unit-cells in a crystal is not unique. There are seven distinct crystal systems, i.e., cubic, hexagonal, rhombohedral (also known as trigonal), tetragonal, orthorhombic, monoclinic, and triclinic. By considering ad- ditional possible lattice points at body centers, face centers, or base centers, whereas not to reduplicate, there are precisely 14 distinct Bravais lattices, e.g., the orthorhombic system includes four di?erent Bravais lattices, i.e., simple orthorhombic, body-centered orthorhom- bic, face-centered orthorhombic, and base-centered orthorhombic lattices; and for the cubic system, there are three distinct Bravais lattices, i.e., simple cubic (sc), body-centered cubic (bcc), and face-centered cubic (fcc) lattices. In the case of body-centered, face-centered, or base-centered Bravais lattice, the primitive unit cell is di?erent from its conventional unit cell. And accordingly, the primitive unit vectors are di?erent from those conventional unit vectors. As an example, for face-centered cubic lattice (fcc), the conventional unit cell is constructed by having atoms on the corners of a cube and an atom at the center of each 16 Table 2.1: Point and space groups of Bravais lattices and crystal structures Bravais lattice Crystal structure (Basis of spherical symmetry) (Basis of arbitrary symmetry) Number of point group 7 crystal systems 32 crystallographic point group Number of space group 14 Bravais lattices 230 space groups face. The conventional unit vectors and primitive unit vectors are ? ???? ? ??? ?? a1 = (a,0,0) a2 = (0,a,0) a3 = (0,0,a) and ? ???? ? ??? ?? a1 =parenleftbig0, a2, a2parenrightbig a2 =parenleftbiga2,0, a2parenrightbig a3 =parenleftbiga2, a2,0parenrightbig (2.3) respectively. The volume of the fcc primitive unit cell is 1/4 of the conventional cubic cell. In addition to translational symmetry as described by the Bravais lattice, a crystal is also invariant under a set of point group symmetry operations, such as rotation, re?ection and inversion. There are totally 32 crystallographic point groups. Each point group can be classi?ed into one of the 7 crystal systems. Point group is also called crystal class. The combination of translational symmetry operations and point group operations in a crystal is referred as the crystal?s space group symmetry. There are a grand total of 230 space groups130. Table 2.1 shows the relation between point and space groups of Bravais lattices and crystal structures. The nomenclature of space group is not unique. One commonly adopted notation is listed in the International Union of Crystallography, which assigns each space group with a number (#1 to #230). The Hermann-Mauguin notation and Sch?on?ies notation are also commonly used in crystallography community. For example, ?-Al2O3 belongs to the space group #167, which can be equivalently noted as R3C or D63d. Experimentally, the crystal structure can be determined by either atomic scale imag- ing or di?raction techniques, such as X-ray di?raction (XRD). To understand the general principles involved in solving di?raction data, it is necessary to introduce the concept of the reciprocal lattice. The reciprocal lattice plays a fundamental role in studies of functions 17 with the periodicity of a Bravais lattice. Once a Bravais lattice in real space is given, one can construct the reciprocal lattice correspondingly. The primitive unit vectors in reciprocal space are de?ned as ? ??? ?? ??? ?? b1 = 2? a2?a3a1?(a2?a3) b2 = 2? a3?a1a2?(a3?a1) b3 = 2? a1?a2a3?(a1?a2) (2.4) where ai (i = 1,2,3) are primitive unit vectors of real space Bravais lattice. Under this de?nition, ai ?bj = 2??ij (i,j = 1,2,3). Any reciprocal lattice vector can be written as G = k1b1 + k2b2 + k3b3 (2.5) where ki (i = 1,2,3) are integers. Then G?R = 2? (n1k1 + n2k2 + n3k3) (2.6) and eiG?R = 1 (2.7) where R = n1a1 + n2a2 + n3a3 is any real space lattice vector. For any family of lattice planes separated by a distance d, there are reciprocal lattice vectors normal to the planes. The shortest one has a length of 2?/d. It is convenient to describe the orientation of lattice planes based on this relation. The commonly used notation is Miller indices. The Miller indices for a family of lattice planes are the integral coe?cients of the shortest reciprocal lattice vector perpendicular to the planes in terms of the primitive unit vectors in reciprocal space, e.g., the reciprocal lattice vector hb1 + kb2 + lb3 determines a family of planes with Miller indices (h,k,l). The di?raction pattern of X-rays incident on a crystal provides information on inter- planar spacing, and ultimately the space group and structure of the crystal. Generally, there are two equivalent ways to explain the phenomenon of X-ray di?raction by a perfect 18 crystal, i.e., the Bragg law and von Laue condition. W. L. Bragg simply treated the crystal as parallel planes of atoms and the incident waves are specularly (mirror-like) re?ected from these planes. The condition for a constructive interference leads to the famous Bragg law 2dsin? = n? (2.8) where d is the distance between two adjacent parallel planes, ? is the angle measured from the plane (90? minus angle of incidence), integer n is the order of interference and ? is the wavelength of incident radiation. The total de?ected angle measured from the incident light is 2?. The von Laue approach regards the crystal as identical lattice points (a set of atoms) which can reradiate the incident wave in all directions. Sharp peaks will be observed at directions where constructive interference happens. The Laue condition can be written as ?k = G (2.9) i.e., the change in wave vector equals to a reciprocal lattice vector. Under the assumption of elastic scattering (the magnitude of the wave vector does not change), the Laue condition yields Bragg law. The practical and detailed experimental methods to determine the crystal structure are beyond the scope of this dissertation. 2.2 Equilibrium Thermodynamic Theory of Phase Transitions 2.2.1 Thermodynamic Stability The occurrence of phase transitions can be interpreted synonymous to changes of atomic structures of matters. This phenomenon has long been studied and many natural forms of transitions are well observed in our everyday life, e.g., water freezes into ice below freezing temperatures. Here we focus our attention to the solid-solid phase transitions, especially the 19 pressure-induced structural phase transitions. Within this limitation, some important as- pects of the complete theory of phase transitions, such as critical phenomena, melting, glass transition, ferromagnetism, superconductivity, etc., are not considered in this dissertation. The ?rst set of questions are: (1) why does a phase transition happen? (2) when does it happen? and (3) how does it happen? The ?rst question is related to the interactions between a solid and its surrounding media. For simplicity, we ignore magnetic and electric interactions, and focus only on mechanical and thermal interactions. According to the 1st and 2nd law of thermodynamics, the most thermodynamically stable state of an isolated state is the state that minimizes the energy (E). When the system is interacting with the surrounding media through mechanical work, the system reaches a mechanical equilibrium with the media when its pressure (P) equals the pressure of the media, and it becomes thermodynamically stable against any spontaneous ?uctuations when its enthalpy (H), de- ?ned as H = E + PV , becomes minimized. Here, V represents the volume of the system. Similarly, when the system exchanges its energy with the media through heating, the ther- mal equilibrium is reached when the temperature (T) of the system equals the temperature of the media, and it becomes thermodynamically stable when its Helmholtz free energy F, de?ned as F = E ?TS, becomes minimized. Here, S is the entropy of the system. In our studies, a solid is considered to be in a thermal equilibrium at temperature T and a mechanical equilibrium at pressure P with the surrounding media. At a given (T,P) condition, the thermodynamically stable state is reached when: ?E + P?V ?T?S ? 0 (2.10) where ?E, ?V , and ?S are virtual variations of total internal energy, volume, and entropy respectively. Given a su?ciently long period of time and allowing all possible ?uctuations, a solid at a given (T,P) condition should reach its thermodynamically stable state, which minimizes its Gibbs free energy G, de?ned as G = E + PV ? TS. Equation 2.10 is the Gibbs-Duhem stability criterion which is equivalent to ?G ? 0. 20 A thermodynamically stable phase of a solid at a given (T,P) condition corresponds to a global minimal of Gibbs free energy for all atomic con?gurations. Correspondingly, metastable phases are referred to those associated local minimums of G. Altering T or P will change the Gibbs free energy landscapes and may in turn shu?e the relative orders among the local G minimums. The emergence of a new global minimal of G will lead to a phase transition in the solid. The answer of when and how a phase transition happens requires detailed knowledge of not only the thermodynamic properties of individual phases but also the kinetic process that connect the initial and ?nal phases. At equilibrium condition, a phase transition becomes possible when a new phase has lower Gibbs free energy than the initial phase. The initial phase then can be considered as a metastable phase at the new (T,P) condition. However, the lifetime of a metastable state is controlled by its local kinetic barrier heights. In reality, metastable states can exist for a long period of time and an equilibrium phase transition can be hindered when there exist large kinetic barriers that against the ?uctuations. A well-known example is the two crystalline forms of element carbon: diamond and graphite. Graphite is known to be the ground state at ambient conditions whereas diamond is also ?stable? in a wide range of usual environment unless enough activation energy (such as heating) is provided. 2.2.2 Phase Diagram and Classification of Phase Transitions When two phases of a single-component system coexist at equilibrium, three conditions must be satis?ed between two phases: 1. The temperatures of the two phases must be equal. T1 = T2 (2.11) 21 2. The pressures in two phases must be equal. P1 = P2 (2.12) 3. The chemical potentials of the two phases must be equal. ?1 = ?2 (2.13) Since the chemical potential is the molar Gibbs free energy, we can rewrite the phase equilibrium conditions in the following form if both phases have the same number of atoms. G1 (T,P) = G2 (T,P) (2.14) The basic concepts of an equilibrium phase diagram is sketched in Figure 2.1. Each point in the two-dimensional (T,P) space denotes one equilibrium state of the system. In most the T-P region of the phase diagram, only a single phase is thermodynamically stable. At the boundary of two adjacent regions, two phases are equally stable along the T-P phase boundary lines. Phase transitions happen when the phase boundaries are crossed. At triple points, such as the (T1,P1)and (T2,P2) points in the plot, three phases coexist. For a single component system, no more than three phases can reach simultaneous thermodynamic equilibria. There is a special point C where the phase boundary curve between the liquid and gas phases come to an end. This point is named as critical point, at and beyond which the liquid and gas phases become identical. It should be mentioned that the critical point can exist only for phases which are quantitatively di?erent, e.g., a liquid and a gas. On the contrary, solid phases have certain internal symmetries. The di?erence between two solid phases of a substance, or between the liquid and solid phases, are qualitatively di?erent. During a quasi-static process of phase transition, the system may absorb or release the so-called latent heat (L). The latent heat is associated with the change of entropy of the 22 P T C Liquid Gas Solid 1 Solid 2 (T1,P1) (T2,P2) (TC,PC) Figure 2.1: Phase diagram with two solid phases, one liquid phase and one gas phase. Two triple points and one critical point are present. system. L = T?S (2.15) where ?S is the entropy di?erence between the ?nal and initial phases. ?S may be zero in continuous phase transitions. If we di?erentiate both sides of equation 2.14 with respect to temperature, the Clapeyron- Clausius equation can be derived. dT dP = V2 ?V1 S2 ?S1 = T?V L (2.16) Clapeyron-Clausius equation gives the slope of the phase boundary in the phase dia- gram which is directly related to the change of volume and latent heat. We can estimate the amount of heat being absorbed/released during the transition from equation 2.16, nev- ertheless both ?V and L are zero for continuous phase transitions (see below). 23 As discussed above, while Gibbs free energy is always continuous during equilibrium phase transitions, other thermodynamic variables might be discontinuous. I now introduce the concept of order of a phase transition based on the Gibbs free energy. 1. First-order transitions are accompanied with a discontinuous change in the ?rst-order derivatives of the thermodynamic potentials, such as volume and entropy. Thus a density jump and a latent heat are presented during a ?rst-order phase transition. Transitions of this type is referred as the ?rst kind or discontinuous phase transitions. V = parenleftbigg?G ?P parenrightbigg T S = ? parenleftbigg?G ?T parenrightbigg P (2.17) 2. In second-order transitions, the ?rst-order derivatives of Gibbs free energies are con- tinuous. Yet, the transitions are accompanied with discontinuity in second-order derivatives of thermodynamic potentials, such as isobaric speci?c heat capacity CP and isothermal compressibility ?T. CP = T parenleftbigg?S ?T parenrightbigg P = ?T parenleftbigg?2G ?2T parenrightbigg P (2.18) ?T = ? 1V parenleftbigg?V ?P parenrightbigg T = ?1/V parenleftbigg?2G ?2P parenrightbigg T (2.19) Higher order transitions can be speci?ed in a similar fashion. Within the framework of Landau theory, second and higher order transitions belong to the second-kind (continuous). For these transitions, system continuously passes from one phase to another without any abrupt changes in volume, entropy and consequently no latent heat. In the case of solid-solid phase transitions in crystals, reconstructive transitions reform the crystal lattices from one type to a distinctly di?erent one without any direct structural relations, e.g., no group-subgroup relations. This type of transitions are always ?rst-order. On the other hand, during displacive phase transitions, atoms gradually deviate from their original equilibrium positions to their corresponding new positions in the ?nal phase in 24 a collective fashion. Although the displacements might be small, crystal symmetry is al- tered. Since the displacements happen in a continuous manner, one structure must has a higher symmetry and the other structure has a lower symmetry which is a subgroup of the spacegroup of the high-symmetry phase. The displacive transitions can be any order 2.3 Landau Theory for Phase Transitions of Second Kind 2.3.1 The Order Parameter and Landau Free Energy Landau theory is a general phenomenological theory for phase transitions of second kind. Landau introduced a quantity order parameter to describe the changes in structure during the phase transition. The order parameter is zero in the high symmetry phase and takes non-zero value in the low-symmetry phase. For displacive phase transitions, the order parameter can be easily chosen as the atomic displacement from the equilibrium sites of the initial high-symmetry phase. Usually the free energy is independent of the sign of the order parameter that it only contains even-power terms. G(T,P,?) = G0 + A?2 + B?4 + ??? (2.20) where G0 is the free energy of the high symmetry phase. Apparently this expansion is valid only for small values of the order parameter, i.e., close to the transition point. For simplicity, we consider equation 2.20 up to the fourth-order term, the sign of coe?cient A will determine whether it leads to a single minimum at ? = 0 or a local maximum at ? = 0. Plot of the free energy as a function of order parameter is shown in Figure 2.2 for these two cases. If A > 0, the global minimum is at ? = 0 which corresponds to the high symmetry phase, while if A < 0, the high symmetry phase become energetically unstable and the system is stable at non-zero values of the order parameter (??0). In this dissertation, we have studied the pressure induced displacive phase transition from phenacite to post-phenacite phase in Si3N4. Because the transition is largely pressure driven, we neglect the temperature dependence of the landau free energy. Therefore at the 25 G(?) A>0 A<0 -?0 ?0 ? Figure 2.2: Variation of the Landau free energy with positive and negative A coe?cients. transition pressure the coe?cient A changes sign, so that it is positive for the high symmetry phenacite phase below the transition pressure and is negative for the post-phenacite phase above the transition pressure. The simplest form of this condition is G(P,?) = G0 + a(Pt ?P)?2 + b?4 (2.21) where a and b both have positive values. We assume A = a(Pt ?P) where Pt is the transition pressure, and B = b. When P > Pt, the equilibrium condition ?G?? = 0 leads to the non-zero value of the order parameter. ?0 = radicalbigg a(P ?Pt) 2b (2.22) Note that the ?rst order derivative of Gibbs free energy with respect to pressure is continuous at P = Pt. 26 2.3.2 Dynamic Lattice Stability and Soft Phonon Modes As a phenomenological approach, Landau theory is mathematically simple and ?exible as the theoretical background of many studies of phase transitions. However, from the ab initio simulation point of view, the new low-symmetry phase is determined from a di?erent approach. For instance, the low-symmetry structure can be derived by relaxing the high- symmetry structure to equilibrium at the conditions where it is unstable. Usually a slight distortion is needed to break the symmetry. For phonon instability (elastic instability is less common), the distortion is a set of collective atomic displacements, which is associated with the soft phonon theory. At ?nite temperature, atoms inside a solid never ?freeze? at their equilibrium posi- tions. Instead, they oscillate around their respective equilibrium positions. The atomic displacements of atoms can be analyzed in term of the vibrational normal modes. The quasi-particle representation of quantization in harmonic lattice vibration is phonon. The squares of phonon (vibrational) frequencies are normally positive. A vanishing ?2 indicates that the restoring forces related to that normal mode disappear. The corresponding vibra- tion is called soft mode or soft phonon. A crystal structure containing soft phonon modes are considered as dynamically unstable because atomic displacements along the eigenmode pattern will lead to a new structure with lower energy. Note that the new structure belongs to a lower symmetry whose space group is a subgroup of the original space group, and the eigenmodes of soft phonons are often adopted as the order parameters in Landau theory discussed above. The simplest case of soft phonon modes is that only one phonon mode becomes softened. If the soft mode drops to zero at ?-point (center of Brillouin zone), the periodicity of the low symmetry structure will be of the same size of unit cell as the high symmetry phase. If phonon softening happens at the zone boundary, the periodicity of the new phase will be di?erent from the high symmetry phase. The unit cell of the low symmetry phase, in general, is a super cell of the previous primitive unit cell. The size of the super cell is determined by the symmetry of the system and the k-point that has soft phonon, e.g., if the soft phonon 27 happens at X-point (1,0,0), the new structure will have a double-sized primitive unit cell vector along a1 compared to the previous phase. Real space displacements of atoms from high symmetry equilibrium sites can be obtained based on the eigenvector of the soft mode. By breaking the symmetry with atomic displacements (usually small), the low-symmetry structure can be found by relaxing the distorted high-symmetry con?guration energetically at the condition that it is unstable. In the case that more than one mode (at di?erent k-points) involve soft phonons, in principle, the real space atomic displacements that lead to the low-symmetry stable phase can be derived from a linear combination of the eigenvectors of the soft modes. In practice, we usually treat the soft phonon modes in serial. For each soft phonon, a new low-symmetry structure can be derived from total energy calculation by imposing its symmetry. In the next step, we calculate the phonon dispersion for the new phase. We repeat the process until all soft phonon modes disappear. The energetically favored structure among these stable or metastable phases is the one with the lowest energy, which can be determined from the E-V curves. 2.4 Transition Paths of Reconstructive Phase Transitions Many structural phase transitions in solids belong to the ?rst-order reconstructive category, i.e., the symmetries of the two structures have no group-subgroup relation and typically bonds breaking/forming process occurs during the transition. The dynamics of reconstructive phase transitions is not yet well understood for most known solid-solid phase transitions. It is known that for some transitions of this kind, an interface forms between the initial and ?nal phases, and atoms go through a long range di?usion process during the transformation. In our study, we examine the transformation mechanism from a di?erent angle, assuming the initial and ?nal structures are related by a di?usionless transition path (TP). Intermediate structures along the path can be given with a crystalline character ap- proximately. Transition pathway studies are helpful for understanding the mechanism of transitions at a microscopic level, and predicting the kinetic barrier height of the proposed 28 path(s) from the energy landscape. However, there are in principle an in?nite number of pathways for one structure to transform into another. If neglecting defects and sur- face e?ects that exist in reality, we assume a certain symmetry is maintained when atoms are displaced collectively along the TP. For di?usionless reconstructive transitions, it is known that the nucleation free path is de?ned by is a common subgroup of both end struc- tures131?134. Although the number of common subgroups is in?nite, it is possible to obtain a ?nite set of common subgroups with certain constrains, such as the size of unit cell, bonds preserving, maximum strains etc. Among these TPs, the energetically favorable path(s) is/are that/those associated with the lowest kinetic energy barrier(s). From the theoretical point of view we could investigate the most possible paths and ?nd the one(s) with the least barrier height. Stokes and Hatch have implemented a com- puter program called COMSUBS which systematically ?nds maximal common subgroups of two structures. The major constraint for COMSUBS is the bond condition, since many studies suggested that the most energetically favored TPs are those preserve the number of bonds53?55. Other constrains that a?ect the number of common subgroups are the minimum and maximum size of unit-cell, strain tolerance, maximum atomic displacement, minimum distance allowed between nearest neighbors, etc. This method has been adopted in the stud- ies of B1-to-B2 transitions in sodium chloride (NaCl) and lead sul?de (PbS)135, zinc-blende to rocksalt transitions in silicon carbide (SiC)136, and wurtzite-to-rocksalt transition48. For each possible transition path, in principle, the energy barrier height is determined by the calculated n-dimensional potential-energy surface (PES), where n is the degree of freedom of the intermediate state. For example, in the case of wurtzite-to-rocksalt transition in AlN, the degree of freedom is 6 for the Pna21 path, i.e., there are 6 free parameters that can be adjusted independently. The complete information of the wurtzite-to-rocksalt transformation path can be revealed from a 6-dimensional PES. However, ?rst-principles calculation of a 6D PES is not an easy task. To calculate the barrier height, two reasonable approaches have been investigated. The ?rst approach is named as bow-function method, as adopted in the study of zinc-blende to rocksalt transition in SiC136. The ?rst step within 29 this method is to calculate the enthalpy along a linear TP. Here ?linear? means that all the structural parameters, including both external parameters (lattice parameters) and internal coordinates, vary according to a single transition parameter. xm = (1 ??)xmi + ?xmf (2.23) where ? is the transition parameter which varies from 0 to 1, xmi and xmf are the initial and ?nal mth structural parameter, respectively. In many cases we have studied, the peak of the barrier height locates at around ? = 0.5. In this way, although we obviously overestimate the enthalpy barrier, it is e?cient for eliminating some TPs that give much larger barriers. With the number of TPs being reduced, we further develop a numerical algorithm to possibly lower the barrier height. On the basis of previous linear TP, we assume the saddle point (which gives the enthalpy barrier height) of the ?true? TP is not far from the maximum point of the linear one. By adding a quadratic term to equation 2.23 we change the linear function to a bow function. xm = (1 ??)xmi + ?xmf + wm? (1 ??) (2.24) where wm is the ?weight? of the bow function for the mth structural parameter. Then we minimize the enthalpy barrier height with respect to wm. This is done by searching for the lowest height of the peak with di?erent weights (positive or negative) for each structural parameter. If more than one peak are present, we will track the highest peak across the entire TP. If the change on wm decreases the barrier height, we keep the change, otherwise we undo the change. This process is done when no changes on wm will further lower the height of the peak. At this situation, saddle point of the true TP is found and so for the true barrier height. We should stress that our aim is to ?nd the barrier height and our proposed TP only matches the saddle point and two end points of the true TP. The quadratic bow function that links these three points is generally di?erent from the true TP. 30 However, the bow-function algorithm may fail if the landscape of the enthalpy function is complex, e.g., if the saddle point is surrounded by peaks and the true path is o? the linear path, the barrier height found by bow function can be higher than the saddle point. In all the concerns our estimated enthalpy barrier height will always be an upper limit compared to the true barrier. Besides the bow function method, sometimes a reconstructive phase transition can be characterized by only a few structural parameters. In this case a proper choice of the struc- tural parameter is essential. Other ?unimportant? parameters will be fully relaxed during the DFT calculation. In the simplest case, if only one structural parameter is primarily responsible during the transition, the transition parameter is de?ned based on the changes of this parameter. If two equally important structural parameters both cause signi?cant changes in enthalpy from one structure to the other, we can calculate the enthalpies as a function of two transition parameters (a 2D PES). Calculations involve more than two transition parameters can be cumbersome in the sense of both computational load and analysis. Let?s take the single parameter case as an example to illustrate the procedure of tran- sition pathway calculations. The particular structural parameter can be found either by comparison of two end structures or by a testing calculation as follows. If we choose the transition parameter at values which are uniform samplings between two end phases (equa- tion 2.23), the saddle point normally locates at ? = 0.5 or nearby (this testing is not applicable to the case where a metastable state exists around ? = 0.5). A series of to- tal energy calculation can be done by ?xing only one structural parameter at a time and fully relaxing the others starting with ? = 0.5. If the parameter we ?xed is unimportant, the resulting energy will be close to the energy of one end phase, whereas the important structural parameter will yield an energy apparently di?erent from the end phases. Next we sample the transition parameter in the [0,1] region and calculate the energies at each sampling point by ?xing the chosen parameter and relaxing the other parameters. To get the enthalpy under a given pressure, we preform the energy calculations at several volume 31 points and obtain the enthalpy based on the equation of states, H = E + PV . Finally a relative enthalpy as a function of transition parameter can be plotted and the barrier height is estimated from that. We emphasize that the calculated barrier height is approximate and should be understood as an upper limit of the true one. For technical reason, it is easy to ?x the internal parameter(s), however, di?cult to do so for the external ones. So the second approach currently only works for transitions that characterized by changes in the internal structural parameter(s). 32 Chapter 3 FIRST-PRINCIPLES CALCULATIONS OF THERMODYNAMIC PROPERTIES OF CRYSTALS 3.1 First-Principles Total Energy Theory 3.1.1 Density Functional Theory (DFT) Atomic-scale theory and simulation of solids rely on accurate evaluation of the total energy of an N-atom system with a speci?ed atomic con?guration. Many theories, evolving from simple empirical models in early years to quantum Monte Carlo method, have been developed to estimate the total electronic energy of a given system. Bene?ting from the increasing power of computational technology and to our purpose, in this study, we adopt density functional theory (DFT). DFT is one of the most successful ?rst-principles quan- tum mechanical theories for atoms, molecules and solids. Nowadays, DFT calculations are routinely performed in the ?elds of materials physics and chemistry. In quantum mechanics, a system is described by the wave function ? of its Hamiltonian. Within Born-Oppenheimer approximation137 (adiabatic approximation), one assumes that the electrons are in their ground state for the instantaneous ionic con?guration at any moment. Therefore, we are able to treat ionic motion with classical mechanics, and only the electronic Hamiltonian is treated with quantum mechanics. A stationary electronic state of a N-electron Hamiltonian is described by a wave function ?(r1,...,rN) satisfying the Schr?odinger equation: ?H? =bracketleftBig?T + ?V + ?UbracketrightBig? = ? ? Nsummationdisplay i ? planckover2pi1 2 2m? 2 i + Nsummationdisplay i v(ri) + summationdisplay i 800 K, Henderson?s data has a nearly linear temperature dependence which 126 is questionable. In most temperature range, our calculated thermal expansion coe?cient from EOS-I is consistent with Bruls? results. s48 s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s48s46s48 s48s46s50 s48s46s52 s48s46s54 s48s46s56 s49s46s48 s49s46s50 s49s46s52 s49s46s54 s32 s32 s84 s104 s101 s114 s109 s97 s108 s32 s101 s120 s112 s97 s110 s115 s105 s111 s110 s32 s40 s49 s48 s45 s53 s32 s75 s45 s49 s41 s84s101s109s112s101s114s97s116s117s114s101s32s40s75s41 s69s120s112s101s114s105s109s101s110s116 s32s72s101s110s100s101s114s115s111s110s32s40s49s57s55s53s41 s32s66s114s117s108s101s115s32s40s50s48s48s49s41 s32s82s101s101s98s101s114s32s40s50s48s48s53s41 s32s83s99s104s110s101s105s100s101s114s32s40s49s57s57s52s41 s84s104s105s115s32s119s111s114s107 s32s69s79s83s45s73 s32s69s79s83s45s73s73 Figure 5.15: Temperature dependence of volume thermal expansion coe?cient of bulk ?- Si3N4 at zero pressure. Solid and dashed lines show present work with static energies ?tted to the 2nd-order BM-EOS and 3rd-order BM-EOS, respectively, and both thermal free energies ?tted to the 2nd-order BM-EOS. Discrete symbols denote experimental data83?86. The negative thermal expansion coe?cient below 150 K can be related to the negative bulk Gr?uneisen parameter ? in that temperature range. ? = ?CV /(BTV ). Figure 5.16 shows our calculated Gr?uneisen parameter of ?-Si3N4 together with reported experimental data. The present temperature dependence of bulk Gr?uneisen parameter is in good agree- ment with Bruls? measured data in the temperature range between 300 K and 1300 K. The estimated percentage di?erence between experiment and calculation is within 10% for 300 K< T <500 K and the di?erence is gradually reduced to about 2% at T = 1300 K. Also, one can notice that the Gr?uneisen parameter is not sensitive to the two EOS schemes we used. The EOS-I scheme yield a larger thermal expansivity and larger volume at high temperatures, while the EOS-II scheme yield a larger isothermal bulk modulus, and both EOS schemes give very similar CV curves at all temperatures. The cancelation e?ect leads to the similarity in the bulk Gr?uneisen parameter. The bulk Gr?uneisen parameter is the weighted average of mode Gr?uneisen ratios. At low temperature, only low-frequency phonons, which mostly have negative mode Gr?uneisen 127 s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s45s49s46s48 s45s48s46s56 s45s48s46s54 s45s48s46s52 s45s48s46s50 s48s46s48 s48s46s50 s48s46s52 s48s46s54 s48s46s56 s49s46s48 s69s120s112s101s114s105s109s101s110s116 s32s66s114s117s108s115s32s40s50s48s48s49s41 s84s104s105s115s32s119s111s114s107 s32s69s79s83s45s73 s32s69s79s83s45s73s73 s71 s114 s252 s110 s101 s105 s115 s101 s110 s32 s112 s97 s114 s97 s109 s101 s116 s101 s114 s84s101s109s112s101s114s97s116s117s114s101s32s40s75s41 Figure 5.16: Temperature dependence of bulk Gr?uneisen parameter of ?-Si3N4 at zero pressure. Solid and dashed lines show present work with static energies ?tted to the 2nd- order BM-EOS and 3rd-order BM-EOS, respectively, and both thermal free energies are ?tted to the 2nd-order BM-EOS. Discrete symbols denote experimental data85. ratios, are thermally excited. This yields the negative overall bulk Gr?uneisen parameters and consequently it leads to the negative thermal expansion coe?cients. The two branches that corresponding to the most negative mode Gr?uneisen parameters are found to be the softening M-point TA and ?-point Bu modes, which are responsible for the instability of ?-Si3N4 at high pressures. Figure 5.17 shows the calculated TEC of ?-Si3N4 based on EOS-II as a function of temperature at several pressures up to 30 GPa. As pressure increases from 0 to 30 GPa, the TEC decreases from 1.11 ? 10?5 to 0.69 ? 10?5 K?1 at 2000 K temperature. The negative TEC range extends from below 150 K at 0 GPa to 220 K at 30 GPa. The most negative TEC value also decreases from ?3.11 ? 10?7 to ?5.09 ? 10?7 K?1. This pressure e?ect of TEC in ?-Si3N4 is in agreement with the calculated pressure e?ect on low frequency phonon modes and the soft-phonon associated structural instability discussed in earlier sections. Isobaric heat capacity per atom as a function of temperature at zero pressure is shown in Figure 5.18. First, both EOS schemes yield very alike curves throughout the plotted tem- perature region, which implies that CP is not sensitive to the equation of states adopted. 128 s48 s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s48s46s48 s48s46s50 s48s46s52 s48s46s54 s48s46s56 s49s46s48 s49s46s50 s49s46s52 s49s46s54 s32s80s32s61s32s32s32s48s32s71s80s97 s32s80s32s61s32s49s48s32s71s80s97 s32s80s32s61s32s50s48s32s71s80s97 s32s80s32s61s32s51s48s32s71s80s97 s32 s32 s84s101s109s112s101s114s97s116s117s114s101s32s40s75s41 s84 s104 s101 s114 s109 s97 s108 s32 s101 s120 s112 s97 s110 s115s105 s111 s110 s32 s40 s49 s48 s45 s53 s32 s75 s45 s49 s41 Figure 5.17: Temperature dependence of volume thermal expansion coe?cient of bulk ?- Si3N4 at pressures of 0, 10, 20 and 30 GPa. Below room temperature our calculations are in excellent agreement with all shown exper- imental results, and the agreement is still good above room temperature except Reeber?s data. Note that Reeber?s data are considerably di?erent from other experimental results and theoretical calculations in both thermal expansivity and isobaric heat capacity. At the same time, we ?nd a persistent agreement with the Bruls? measurements for both CP and TEC. Figure 5.19 shows the temperature dependencies of TEC for ?-, ?- and ?-Si3N4 at zero pressure. The equation of state scheme we adopted is EOS-II. We also present the experimental TEC of ?-Si3N4 in the same plot. Although the two experimental works do not agree well with each other, our prediction is consistent with Paszkowicz et al. below room temperature and the agreement is still reasonable in the temperature range from 300 K to 1000 K. More high quality measurements are required to assess the TEC of ?- Si3N4. Thermal expansion coe?cients of ?- and ?-Si3N4 are similar in magnitude over all temperatures and their di?erence at 2000 K is less than 2?10?6 K?1. In Kuwabara et al.?s calculation, TEC of ? phase was predicted to be slightly smaller than that of ? phase. This might be ascribed to the sensitive nature of TEC upon the numeric ?uctuation of thermal data. Both ? and ? phases present negative TEC at low temperature which is consistent 129 s48 s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s48s46s48 s48s46s53 s49s46s48 s49s46s53 s50s46s48 s50s46s53 s51s46s48 s32 s32 s67 s80 s32 s40 s107 s66 s47 s97 s116 s111s109 s41 s84s101s109s112s101s114s97s116s117s114s101s32s40s75s41 s69s120s112s101s114s105s109s101s110s116 s32s75s117s114s105s121s97s109s97s32s40s49s57s55s56s41 s32s75s111s115s104s99s104s101s110s107s111s32s40s49s57s56s50s41 s32s66s114s117s108s115s32s40s50s48s48s49s41 s32s82s101s101s98s101s114s32s40s50s48s48s53s41 s84s104s105s115s32s119s111s114s107 s32s69s79s83s45s73 s32s69s79s83s45s73s73 Figure 5.18: Temperature dependence of isobaric heat capacity of ?-Si3N4 at zero pressure. Solid and dashed lines show present work with static energies ?tted to the 2nd-order BM- EOS and 3rd-order BM-EOS, respectively, and both thermal free energies are ?tted to the 2nd-order BM-EOS. Discrete symbols denote experimental data85?88. with the ?rst-principles calculation from Kuwabara et al. The lowest TEC for the ? phase is about ?1.5 ? 10?7 K?1 at 90 K, and for the ? phase it is ?3.1 ? 10?7 K?1 at 100 K. Unlike ?- and ?-Si3N4, TEC of ? phase is about twice as large as the that of ?- or ?-Si3N4 in the temperature range from room temperature to 2000 K. At 2000 K, we predict a TEC of 2.3 ? 10?5 K?1 which is in good agreement with Kuwabara et al.?s ab initio calculated 2.2 ? 10?5 K?1. Paszkowicz et al. pointed out that their measured TEC tends to vanish for T < 100 K90. However, our calculation is not clearly supportive. Similarly, two other ?rst-principles calculations from Paszkowicz et al. and Kuwabara et al. do not show the negative or vanishing feature either. 5.7 Conclusions In summary, we have theoretically studied phase transitions in silicon nitride (Si3N4) at high pressure using a ?rst-principles density functional theory method. We ?nd that ?-Si3N4 remains as a metastable phase at temperatures up to 2000 K and pressures up to 10 GPa. The equilibrium ??? transition pressure is predicted as 7.5 GPa at 300K 130 s48 s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s48s46s48 s48s46s53 s49s46s48 s49s46s53 s50s46s48 s50s46s53 s32 s32 s84 s104 s101 s114 s109 s97 s108 s32 s101 s120 s112 s97 s110 s115s105 s111 s110 s32 s99 s111 s101 s102 s102 s105 s99 s105 s101 s110 s116 s32 s40 s49 s48 s45 s53 s32 s75 s41 s84s101s109s112s101s114s97s116s117s114s101s32s40s75s41 s69s120s112s101s114s105s109s101s110s116s32s40 s45s83s105 s51 s78 s52 s41 s84s104s101s111s114s121s32s40s116s104s105s115s32s119s111s114s107s41 s32s72s105s110s116s122s101s110s32s40s50s48s48s51s41 s32 s32s112s104s97s115s101 s32s80s97s115s122s107s111s119s105s99s122s32s40s50s48s48s52s41 s32 s32s112s104s97s101 s32 s32s112s104s97s115s101 Figure 5.19: Temperature dependence of volume thermal expansion coe?cient of ?-, ?- and ?-Si3N4 at zero pressure. Discrete symbols represent experimental data89,90. and it increases to 9.0 GPa at 2000K. Both ?- and ?-Si3N4 are dynamically stable at low pressure. However, two competing phonon-softening mechanisms are found in the ? phase at high pressures. At room temperature, ?-Si3N4 is predicted to undergo a ?rst-order ??P3 transition above 38.5 GPa, while ?-Si3N4 shows no signs of structural instability. The predicted metastable high-pressure P3 phase is structurally related to ?-Si3N4. The enthalpy barrier height is estimated as only 67.23 meV/atom. Our LDA predicted thermal expansion coe?cient, heat capacity and bulk Gr?uneisen parameter are in good agreement with Bruls? measured results. We ?nd relatively large discrepancies between our calculation with experimental data from Reeber. And we attribute the cause of predicted negative TEC at low temperatures in ? and ?-Si3N4 to the low-frequency phonon modes that have negative mode Gr?uneisen ratios. 131 Chapter 6 FIRST-PRINCIPLES STUDY OF GALLIUM OXIDE AND GALLIUM OXYNITRIDE 6.1 Gallium Oxide: Ga2O3 6.1.1 Introduction The oxides of group 13 elements (Al, Ga, In) are important solid-state compounds with applications in ?elds ranging from structural ceramics to catalysts and electronic materi- als201. Monoclinic gallium oxide (Ga2O3) is usually known as a wide-band-gap semiconduc- tor (Eg = 4.9 eV) ; however, the conductivity can be varied from insulating to conducting behavior depending upon the preparation conditions93. Due to its tunable optical and elec- tronic properties, ?-Ga2O3 is being developed for use in a wide variety of applications, for instance, as optical windows202, in high-temperature chemical gas sensors203, as a mag- netic memory material204, and for dielectric thin ?lms205. Recently, considerable e?ort has been devoted to the study of low-dimensional Ga2O3 materials, and ?-Ga2O3 nanowires have been obtained through physical evaporation and arc-discharge methods101. ?-Ga2O3 has also attracted recent interest as a phosphor host material for applications in thin ?lm electroluminescent displays206,207. Due to its chemical and thermal stability, ?-Ga2O3 may emerge as a useful alternative to sul?de based phosphors208. It is well known that Ga2O3 can exist in several forms, including ?, ?, ?, ?, and ? poly- morphs that all have di?erent structure types94. Of these, the most stable form at ambient conditions is determined to be ?-Ga2O3 (monoclinic C2/m, Figure 6.1)94. However, other metastable varieties can be prepared and they have been characterized at ambient pressure and temperature. This is an important observation, because the di?erent forms have dra- matically di?erent optoelectronic properties. For example, the band gap of the ?-Ga2O3 132 Figure 6.1: Crystal structure of monoclinic ?-Ga2O3 phase. Note that Ga3+ cations (light color) occupy both tetrahedral and octahedral interstices within the ccp lattice of O2? ions (dark color). polymorph that is isostructural with corundum (?-Al2O3) is 2.41 eV, much narrower than that of ?-Ga2O3209. It is of great interest to determine the pressure-induced phase transformations among Ga2O3 polymorphs in order to establish the stable and metastable phase relations between di?erent crystalline modi?cations, and to evaluate their production under di?erent synthe- sis conditions. It is particularly important to understand the role of di?erential mechanical stresses that are present in creation of nanoparticles or nanowires, in promoting the forma- tion of speci?c polymorphic forms. The high-pressure behavior of Al2O3 compounds has been studied extensively, par- ticularly the corundum-structured ?-Al2O3 phase, because of its importance as a mineral structure within the deep Earth and also due to the widespread use of ruby (Cr3+-doped ?-Al2O3) as a luminescent pressure gauge for in situ high-pressure experiments in the di- amond anvil cell1. Cr3+-doped ?-Ga2O3 has likewise been proposed as a pressure gauge material. The R1 luminescence line in this phase shows a pressure shift nearly three times that of ruby, indicating that it would make a more sensitive pressure sensor that is especially 133 useful in the lower pressure range210. In situ high-pressure and high-temperature measure- ments on ?-Al2O3 using synchrotron x-ray di?raction in a diamond anvil cell, combined with ab initio theory predictions, have now been used to characterize a transition into the Rh2O3-II structure occurring at P ? 100 GPa and T >? 1000 K5,8. The high-pressure behavior of Ga2O3 has received much less attention. The vari- ous low-density Ga2O3 structures encountered at low pressure contain the Ga3+ cations in tetrahedral coordination (i.e., GaO4 species). The thermodynamically stable ?-Ga2O3 poly- morph is isomorphous with the metastable ?-Al2O3 structure, which represents a key phase achieved during metastable transformations among various partially dehydrated ?transi- tional? aluminas as they evolve towards corundum211. ?-Al2O3 constitutes an intermediate structure between the cubic close packing of anions achieved within the low-temperature metastable aluminas, and hexagonally close-packed ?-Al2O3 corundum (isomorphous with ?-Ga2O3). In a recent study using synchrotron energy-dispersive x-ray di?raction techniques in the diamond anvil cell, it was reported that a sample of ??-Ga2O3? transformed to a structure assigned to be tetragonal at a pressure of approximately 13.3 GPa.212 However, the x-ray di?raction pattern of the starting material most strongly resembled that of ?- Ga2O3, rather than the ?-form, and a mixture of phases was present. Commercial Ga2O3 samples usually consist mainly of ?-Ga2O3, along with some ?-Ga2O3; that phase can be removed by heat treatment94. The relative densities of ?- and ?-Ga2O3 are 5.94 and 6.48 g?cm?3, respectively95, indicating that a ??? transformation should occur at high pressure. Nanocrystalline ?-Ga2O3 particles embedded in a glassy matrix were also studied at high pressure using energy-dispersive x-ray di?raction96. In that work, a ?-to-? phase transformation was found to be initiated at 6 GPa, but the process was not completed by 15 GPa, the highest pressure achieved in the study. However, it is known that the silica glass host matrix undergoes important structural and density changes within this pressure range97,98, so that it is not yet known if the structural changes are intrinsic to the ?-Ga2O3 material presumably in?uenced by the nanocrystalline nature of the sample , or 134 are promoted by anomalous densi?cation among the SiO2 matrix. These results prompted us to theoretically investigate the high-pressure behavior of the phase-pure bulk ?-Ga2O3, accompanied with an experimental study using Raman spectroscopy and high-resolution synchrotron x-ray di?raction angle dispersive techniques. Recently, two further pressure-induced phase transitions have been con?rmed from both ?rst-principles calculations and high-pressure x-ray di?raction measurements using laser-heated diamond-anvil cell (LHDAC). The further transition sequence found in Ga2O3 is the same as that in Al2O3, i.e., ??Rh2O3(II)?postperovskite (Cmcm CaIrO3-type). The experimentally determined Pt for these two transitions are ?37 GPa at 2000?100 K and 164 GPa at 1300?500 K. One-dimensional nanostructured forms of ?-phase of gallium oxide (?-Ga2O3) such as nanotubes, nanobelts, and nanowires, have attracted recent interest due to enhanced optical properties99,100. Recently, Choi et al.101 synthesized ?-Ga2O3 nanowires (diameter range of 15?45 nm) with a [001] growth direction using an arc-discharge method. Gao et al.102 synthesized [40?1] ?-Ga2O3 nanowires with diameters ranging from ?10?100 nm in a vertical radio-frequency furnace. Interestingly, the Raman mode frequencies of the [001] ?-Ga2O3 nanowires coincide with the corresponding frequencies in bulk ?-Ga2O3101. On the other hand102, the Raman mode frequencies of the [40?1] ?-Ga2O3 nanowires are redshifted relative to corresponding frequencies in bulk ?-Ga2O3 by 4?23 cm?1. Using plasma-enhanced chemical vapor deposition, Rao et al. have synthesized ?-Ga2O3 nanowires whose growth is along the [110] direction103, and the Raman spectrum is signi?cantly blueshifted in frequency104. Here we focus on the ?rst-principles calculations of the Raman mode frequencies under internal strains. Our calculated Raman frequency shifts suggest that the observed shifts in the nanowires with the [40?1] and [110] growth directions can be explained in term of di?erent internal strains, in contrast to the previously suggested quantum con?nement e?ects and defect-induced e?ects. 135 Table 6.1: Third-order BM-EOS parameters for Ga2O3 polymorphs Source Phase V0 (?A3/atom) B (GPa) B? This work ? 10.134 177 3.84 Calculation214 10.228 142 4.1 Experiment213 202(7) 2.4(6) Experiment96 191(5) 8.3(9) Experiment214 10.462 134(12) 4 This work ? 9.375 226 3.971 Calculation214 9.448 243 4.0 Experiment213 ?250 Experiment214 9.632(1) 223(2) 4 This work Rh2O3(II) 9.130 234 3.98 Calculation214 9.188 244 4.3 Experiment214 9.244(4) 271(10) 4 6.1.2 Total Energy and Vibrational Properties During the course of the experimental studies from our collaborators104,213, we began a parallel theoretical investigation of the structures, relative energetics, and properties of Ga2O3 phases. Our calculations are based on the ?rst-principles density functional theory within the local density approximation (LDA). The calculations were carried out with the VASP codes139?142, using planewave basis sets and ultrasoft pseudopotentials (US-PP)121. In this study, both valence (4s4p) and semi-core (3d) electrons in Ga atoms were treated explicitly, while the core electrons were approximated with the US-PP. The energy cuto? of the plane-wave basis was chosen as 396 eV. The Brillouin zone integration of total energy of the unit cells was carried out using 6 ? 6 ? 6 grids for both phases. Equilibrium volume V0, bulk moduli B and B? for ?, ? and Rh2O3(II) phases are obtained from the least-square ?tting to the 3rd-order Birch-Murnaghan equation of state. Fitting parameters V0, B and B?, are listed in Table 6.1. Our LDA results agree reasonably with other calculated and experimental data. Compared with the experimental values of V0, our prediction has an underestimation of less than 3%. The static enthalpies of three polymorphs relative to the ? phase are plotted in Figure 6.2 as a function of pressure. It 136 s48 s49s48 s50s48 s51s48 s52s48 s53s48 s54s48 s45s48s46s49s48 s45s48s46s48s53 s48s46s48s48 s48s46s48s53 s48s46s49s48 s48s46s49s53 s48s46s50s48 s48s46s50s53 s69 s110 s116 s104 s97 s108 s112 s121 s32 s68 s105 s102 s102 s101 s114 s101 s110 s99 s101 s32 s40 s101 s86 s47 s97 s116 s111 s109 s41 s80s114s101s115s115s117s114s101s32s40s71s80s97s41 s32 s32s112s104s97s115s101 s32 s32s112s104s97s115s101 s32s82s104 s50 s79 s51 s40s73s73s41 s71s97 s50 s79 s51 Figure 6.2: Enthalpy di?erences relative to ? phase for ?, ? and Rh2O3(II)-Ga2O3 as a function of pressure up to 60 GPa. The ? and ? curves cross at 0.5 GPa, the ? and Rh2O3(II) curves cross at 40.9 GPa. can be seen that ? phase is predicted to transform into ? phase at 0.5 GPa and ? phase further transforms into Rh2O3(II) phase at 40.9 GPa. The phonon dynamical matrices Dij (k) were constructed at the Brillouin zone center (?-point: k = 0) using a (realspace) force constant matrix ?ij (r) by calculating forces on each atom as it is slightly displaced from its equilibrium position (e.g., by 0.015 ?A). Further approximations were adopted to calculate the Dij (k) matrix at a general k-point. In the case of ?-Ga2O3, we ?rst obtained the real space ?ij (r) matrix using a 120-atom supercell model. Because of the large size of the supercell model, and the fact that the material is insulating wide-gap semiconducting , we can safely neglect interatomic interactions between atoms separated by >50% of the supercell lattice constants, and thus obtain Dij (k) by Fourier transformation of the real-space matrix elements ?ij (r). Fifteen Raman-active modes are expected for the ?-Ga2O3 structure (point symmetry C32h) from symmetry analysis: ?Raman = 10Ag + 5Bg (6.1) 137 Table 6.2: Calculated and experimental zone-center Raman peak positions and Gr?uneisen parameter for the ?-Ga2O3 phase Frequency (cm?1) Gr?uneisen ratio calculated measured Mode This work Empirical This work symmetry (LDA) calculation215 Machon et al.213 Rao et al.104 (LDA) (measured213) Ag 104 113 110.2 1.39 Bg 113 114 113.6 -0.7 Bg 149 152 144.7 144 1.53 1.97(8) Ag 165 166 169.2 169 1.00 0.35(3) Ag 205 195 200.4 200 1.30 0.98(2) Ag 317 308 318.6 317 1.13 0.95(1) Ag 346 353 346.4 344 1.83 1.52(1) Bg 356 360 1.47 Ag 418 406 415.7 416 0.58 0.78(4) Ag 467 468 472 1.26 Bg 474 474 473.5 1.14 1.27(9) Ag 600 628 629 1.70 Bg 626 644 628.7 0.8 1.54(3) Ag 637 654 652.5 654 1.39 1.39(2) Ag 732 760 763.9 767 1.23 1.11(1) The calculated frequencies and their mode Gr?uneisen parameters are reported in Ta- ble 6.2. The Raman-active modes of ?-Ga2O3 can be classi?ed into three groups: high- frequency stretching and bending of GaO4 tetrahedra (?770?500 cm?1), midfrequency de- formation of Ga2O6 octahedra (?480?310 cm?1) , and lowfrequency libration and transla- tion (below 200 cm?1) of tetrahedra-octahedra chains215. The listed experimental values of Machon et al. are from an unpolarized Raman spectrum of ?-Ga2O3 recorded from pow- dered material obtained by annealing a commercial sample213. Our calculated frequencies agree well with the observed values, to within 0.1%?6%, which is typical for LDA calcula- tions (Table 6.2). The Bg mode predicted at 356 cm?1 was likely unresolved from the Ag mode at 346 cm?1 in the measured Raman spectrum213; however, a weak peak at this fre- quency was recorded by Dohy et al.215. The band observed at 474 cm?1 also likely contains contributions from the calculated Ag and Bg modes at 467 and 474 cm?1. Surprisingly, we found no experimental evidence for the predicted Ag mode at 600 cm?1. We have no explanation for that observation. Symmetry analysis indicates that seven Raman active modes (2A1g + 5Eg) are expected for the corundum structure. Our theoretically calculated frequencies agree to within 0.7%?5.5% with the observed values (Table 6.3). 138 Table 6.3: Calculated and experimental zone-center Raman peak positions for the ?-Ga2O3 Frequency (cm?1) Mode symmetry (calculated) (measured) A1g 215 217.4 Eg 239 240.8 Eg 281 286.1 Eg 344 328.7 Eg 410 432.2 A1g 551 573 Eg 680 688.1 6.1.3 T-P Phase Diagram The monoclinic ?-Ga2O3 structure is the stable polymorph at ambient pressure and temperature. In this phase, the O2? anions form a slightly distorted fcc lattice and cations occupy tetrahedral and octahedral interstices (Figure 6.1). This structure is quite di?erent from that of the ?-Ga2O3 phase (corundum structure), which is based on a distorted hcp O2? sublattice with 2/3 of the octahedral interstices occupied by Ga3+ ions . The ??? transition is expected to result from increasing the pressure, from the observed density relationships between the two phases. The transformation involves a change in the O2? packing from cubic to hexagonal, accompanied by a shift in Ga3+ ions between tetrahedral and octahedral sites. The reconstructive nature of the transition indicates that it is ther- modynamically of the ?rst order, and it might be expected to involve a large activation energy, which gives rise to slow transformation kinetics at low temperature. In addition, the high-density structure has signi?cant possibilities for disorder among the Ga3+ positions on octahedral sites within the hcp O2? sublattice, which might not be readily detected by x-ray di?raction8. ?-? and ?-Rh2O3(II) phase boundaries of Ga2O3 are plotted in Figure 6.3. The cal- culated transition pressure from ? to ? phase is 0.3 GPa at 0 K and 1.6 GPa at 2000 K. The Clapeyron slope for this transition is positive which has the value of about +0.6 MPa/K at 1000 K. Our predicted Pt is consistent with reported LDA calculation from Yusa et al 214. On the experimental side, ?-Ga2O3 was synthesized from ?-Ga2O3 at 4.4 GPa 139 s48 s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s48 s49s48 s50s48 s51s48 s52s48 s53s48 s54s48 s32 s32 s80s114s101s115s115s117s114s101s32s40s71s80s97s41 s84 s101 s109 s112 s101 s114 s116 s117 s114 s101 s32 s40 s75 s41 s82s104 s50 s79 s51 s40s73s73s41 s71s97 s50 s79 s51 Figure 6.3: T-P phase diagram for ?, ? and Rh2O3(II)-Ga2O3 polymorphs. and 1000 ?C95,216. However, the room temperature compression experiments obtain much higher transition pressures. Tu et al. have reported a ?-to-? phase transition in Ga2O3 at 13.3 GPa212. But there is some uncertainty in the nature of the starting material used in that work. Recently, x-ray di?raction and Raman scattering results from Machon et al.213 on ?-Ga2O3 clearly indicate that a pressure-induced phase transformation occurs within the P = 20 ? 22 GPa range, and perhaps as low as P = 18.5 GPa. Analysis of the x-ray di?raction data suggest that the high-density phase corresponds to corundum-structured ?-Ga2O3. However, broadening observed both in the x-ray di?raction peaks and in the Raman spectra indicate that the material is structurally disordered. The large discrepancy on Pt may be caused by a large kinetic barrier, which is signi?cantly lowered under com- pression. Our calculated ?-to-Rh2O3(II) transition pressure is 40 GPa at 0 K and 37 GPa at 2000 K, which is in excellent agreement with the experimental value (about 37 GPa at 2000?100 K)214. The phase boundary we predict shows negative Clapeyron slope which is about -1.65 MPa/K at 1000 K. Tsuchiya et al. reported their calculated Clapeyron slope of the ?-Rh2O3(II) boundary, which is -2.2 MPa/K at 1000 K217. The discrepancy may be ascribed to the adopted computational methodologies. 140 6.1.4 Blue-Shifted Raman Scattering in Gallium Oxide Nanowires Rao et al.104 reported the Raman and Fourier transform infrared spectra of ?-Ga2O3 nanowires with [110] growth direction which is blueshifted relative to the bulk spectra by 10?40 cm?1. The blueshift in phonon frequencies of low-dimensional materials are often attributed to the size-con?nement e?ect218,219. However, the average diameter of the ?- Ga2O3 nanowires is around 25 nm. It is unlikely that the quantum size con?nement at this length scale is signi?cant enough to cause the phonon shifts as large as 50 cm?1. Furthermore, three distinctly di?erent shift patterns have been experimentally observed for the ?-Ga2O3 nanowires of di?erent growth directions. In contrast to the blueshift in the Raman and FTIR spectra reported by Rao et al., Choi et al. showed that their Fourier transform Raman spectrum of [001] ?-Ga2O3 nanowires is identical to that of bulk ?- Ga2O3101, while Gao et al. exhibited a redshift of 4?23 cm?1 in the Raman peak frequencies of their [40?1] ?-Ga2O3 nanowires relative to the corresponding Raman frequencies in bulk ?-Ga2O3102. The size con?nement e?ect is clearly insu?cient to explain the diversity of the observed shift patterns. On the other hand, the redshift in the phonon frequencies has also been attributed to the presence of impurities and defects, such as point defects, twins, and stacking faults220. These defects are also likely to be responsible for additional vibrational modes observed in the Raman spectra (and to a small extent in the FTIR spectrum) of nanowires. From a detailed high-resolution transmission electron microscopy (HRTEM) study, Gao et al. con?rmed the presence of twins and edge dislocations in their nanowires102. Dai et al.221 also proposed that the O vacancies and the stacking faults caused an abnormality in the Ga?O bond vibration and led to redshift in the Raman frequencies. Although this simple hypothesis is plausible, there is one obvious weakness, i.e., lack of close correlation between defect types and the growth directions. Presumably, similar defects might exist in the nanowires with di?erent growth directions. It is also not clear which types of defects will lead to a blueshift in vibrational frequencies. Moreover, di?erent regions in the nanowires contain di?erent defects which would imply that di?erent shifts in the Raman and/or IR spectra 141 should be observed when di?erent regions of the same nanowire are probed. However, this does not seem to be the case and instead overall distinct blueshifts or redshifts have been observed for a given nanowire. Therefore, alternative models that are capable of describing these diverse peak-shift patterns in a consistent fashion are needed. Based on a ?rst-principles calculation which we describe next, we propose that the phonon frequencies in di?erent ?-Ga2O3 nanowires are shifted as a result of internal strains in the nanowire. The basic assumption of our model is the presence of non-negligible internal strains in the nanowires due to their large surface/volume ratio. Di?erent growth direc- tions will cause di?erent surface reconstruction, and consequently lead to internal strains of di?erent magnitudes and directions. This model provides a consistent explanation for all three aforementioned Raman spectra. Direct ?rst principles calculations of phonon frequencies of 25-nm-diam nanowires is a computationally challenging task as large supercell models (of at least tens of thousands of atoms) are needed. Instead, our current computation study focuses on providing a quan- titative estimation of the internal strains which can account for the observed blue- and redshifts in the Raman frequencies for [110] and [40?1] ?-Ga2O3 nanowires, respectively. We have calculated the strain dependencies of the bulk ?-Ga2O3 using a density functional theory (DFT) method. The internal strains of the nanowires were estimated based on the least-squares ?tting of the experimentally observed Raman frequency shifts with theoret- ically predicted linear strain coe?cients d?/d?ij, where ? and ?ij are Raman frequencies and components of strain tensors, respectively. The ?-point phonons of bulk ?-Ga2O3 were calculated with a real-space ?nite displace- ment technique129. Because of its C2h space-group symmetry, the LO-TO splitting in the optic modes of ?-Ga2O3 only exist for the infrared (Au and Bu) phonon modes, not the Raman active (Ag and Bg) phonon modes. Therefore, all our Raman frequency calculations of ?-Ga2O3 were carried out with 10-atom base-centered monoclinic unit-cell model without the correction for the macroscopic interaction. As shown in Table 6.2, the theoretical data for bulk ?-Ga2O3 matches well with 13 out of the 15 Raman active (10Ag + 5Bg) modes 142 Table 6.4: Estimated internal strains [110] [40?1] Strain nanowire nanowire ?11 -0.0077 0.0029 ?22 0.0180 -0.0064 ?33 -0.0311 0.0106 ?13 0.0233 -0.0256 ?V/V -0.0208 0.0071 observed experimentally104, as well as those of the previous study of Dohy et al 215. In both cases of the unobserved Raman modes, there is another Raman active mode in the close proximity. For example, our LDA calculations predicted two Raman modes at 469 and 474 cm?1, and two Raman modes at 601 and 629 cm?1. This suggests that it is possible that the two ?missing? Raman modes are hidden by the stronger adjacent Raman modes. The strain tensor of this monoclinic crystal has six independent elements, ?11, ?22, ?33, ?23, ?13 and ?12. For simplicity, we restricted this study to linear e?ects, i.e., ? (?) ? ?0 +summationtext(d?/d?ij)??ij. This approximation is valid for small strains. We further neglected the strain of ?23 or ?12 because their d?/d?ij coe?cients are zeroes due to the monoclinic lattice symmetry. For each of four remaining types of strains (?11, ?22, ?33, and ?13), the Raman frequencies were calculated for ?ve ?nite strain values between -0.02 and +0.02. The calculated frequencies were then ?tted with a polynomial function to obtain the linear strain coe?cients. Fitting the experimental data within our strain-induced phonon shifts model, we predict the internal strains in nanowires which showed the three distinct Raman spectra. The results of the [40?1] and [110] nanowires are listed in Table 6.4, and our model predicts the strain tensor for the [001] nanowires contain non-negligible ?11, ?22, ?33, and ?13 components. As shown in Table 6.5, we obtain overall excellent ?ts for both the redshifted and blueshifted Raman spectra, with exception of the 134 cm?1 Bg mode in the [40?1] nanowire (Gao et al.). Our calculation shows that the [110] nanowire is compressed along its a and c axis, and stretched along its b axis. The strain in the [40?1] nanowire exhibits a contrasting pattern and its strain magnitude is only about 1/3 of that evaluated for the [110] nanowire. In both 143 Table 6.5: Raman mode frequencies and frequency shifts in ?-Ga2O3 nanowires with the [40?1] and [110] growth directions. Overall, excellent agreement between the observed and calculated shifts is seen for all mode frequencies except the one marked with an *. Gao et al.102 This work [40?1] growth direction [110] growth direction Frequency Calculated Frequency Calculated Bulk Nanowire shifts frequency Bulk Nanowire shifts frequency frequency frequency ?? shifts frequency frequency ?? shifts (cm?1) (cm?1) (cm?1) (cm?1) (cm?1) (cm?1) (cm?1) (cm?1) -4.8 12.0 1.7 -3.1 142 134 -8* -1.0 144 3.8 167 160 -7 -7.2 169 180 +11 10.6 198 194 -4 -6.5 200 213 +13 13.5 320 -3.9 317 9.5 344 332 -12 -6.7 344 16.8 -1.7 0.7 415 409 -6 -5.3 416 428 +12 12.3 473 -4.4 472 492 +20 21.1 -6.2 18.5 627 -5.7 629 645 +16 17.8 -0.3 3.6 651 641 -10 -14.8 654 697 +43 36.4 765 742 -23 -20.6 767 810 +43 47.0 cases, the a axis has the smallest change (Table 6.4). The strain-induced volume changes are predicted to be -2% and 0.7% for the [110] and [40?1] nanowires, respectively. Seo et al.222 studied the internal strains of GaN nanowires using x-ray measurements and they reported the strains of ?xx=2.3%, ?yy=-0.734%, and ?zz=-0.4% based on their experimental x-ray measurement. The magnitudes of our predicted strains of ?-Ga2O3 are comparable to those of GaN nanowires. 6.1.5 Conclusions For bulk Ga2O3, we have calculated the Raman frequencies and their mode Gr?uneisen parameters for ? and ? phases. Good agreement is achieved between our results with experiments and other calculations. We also predict the equilibrium T-P phase diagram of Ga2O3 consists of ?, ? and Rh2O3(II) phases. Our LDA calculated transition pressure from ? phase to ? phase is 0.3 GPa at 0 K and 1.6 GPa at 2000 K. The Clapeyron slope 144 for this transition is positive which has the value of about +0.6 MPa/K at 1000 K. This is consistent with previous calculations. However, the experimental Pt at room temperature is much larger (20?22 GPa). The discrepancy can be attributed to the existence of a large kinetic barrier, which may be signi?cantly lowered under compression. Our calculated ?-to- Rh2O3(II) transition pressure is 40 GPa at 0 K and 37 GPa at 2000 K, which is in excellent agreement with the experimental value (about 37 GPa at about 2000 K). In the case of nanowires, based on a comparison of the experimental Raman mode frequencies with our ?rst-principles calculations, we ?nd compelling evidence for growth direction induced internal strains in ?-Ga2O3 nanowires which signi?cantly in?uence the vibrational mode frequencies. Within the linear model approximation, the observed blue and redshifts of peak frequencies in the micro-Raman spectra of the ?-Ga2O3 nanowires with di?erent growth directions can be attributed to two small anisotropic internal strains: one compressive strain of 2% volume change, and the other tensile strain of 0.7% volume change. The overall high quality of the ?tted models to available experimental data suggests a strong correlation between the shifts in Raman mode frequencies and the growth direction-induced internal strains in the Ga2O3 nanowires. 6.2 Gallium Oxynitride: Ga3O3N 6.2.1 Introduction Semiconducting materials with a wide bandgap are of interest for applications in high- temperature electronics and in optoelectronic devices, particularly when the gap is direct. The oxides and nitrides of gallium are wide-gap semiconductors that give rise to materials that are useful in the blue to UV range at short wavelength. The cubic (sphalerite) and hexagonal (wurtzite) forms of GaN have bandgap energies of 3.3223 and 3.4 eV224, respec- tively. When alloyed with In and Al, hexagonal GaN-based materials with a direct gap have been developed for use in light-emitting diodes and lasers at wavelengths extending 145 from the blue to the ultraviolet range225. The thermodynamically stable ?-Ga2O3 poly- morph has a direct gap of 4.7 eV, and it has also been proposed for development as a solid-state LED material for UV applications226. The group 13 oxynitride materials have other useful properties related to their electronic structure. ?-Ga2O3 with the corundum structure is conveniently alloyed with Al2O3 to provide selective reduction catalysts for gaseous NOx105, and various other Ga2O3 phases have been proposed as gas sensors, and in nanoscale structures as electron emitters and magnetic memory materials106. Within the Al2O3?AlN system, several important AlxOyNz ceramic alloys and compounds are known. At high AlN contents, layered forms based on hexagonal/cubic intergrowths are present. As the Al2O3 content is increased, cubic spinel-structured materials begin to appear. A large family of defect spinels (?-Al2O3, AlxOyNz) contain vacancies on both cation and anion sites107. A stoichiometric oxynitride spinel-structured compound is obtained at the Al3O3N composition, in which Al3+ ions are present on the octahedral and tetrahedral sites, and O2? and N3? occupy tetrahedral anion sites108,109. Among the related nitride compounds Si3N4 and Ge3N4, high-pressure synthesis has recently resulted in formation of a new class of spinel structures, that contain Si4+ and Ge4+ cations on both tetrahedral and octahedral sites110?114. The new solid-state compounds are recoverable to ambient conditions, and they possess high hardness and low compressibility, comparable with materials such as Al2O3-corundum107. The new group 14 nitride spinels are also predicted to be wide direct band semiconductors, with band gaps calculated to lie within the range 2.2-4.0 eV82,227. The analogous spinel-structured compound Sn3N4 has also been made at ambient pressure228,229, raising the possibility of future preparation of ?-(Si,Ge)3N4 ?lms via metastable synthesis routes, such as chemical vapor deposition, to yield materials compatible with optoelectronics applications, for example. Gallium oxynitride (Ga3O3N) has been predicted to form a new spinel-structured com- pound within the Ga2O3?GaN system, with potentially useful electronic properties115,116. It is predicted to be a direct wide bandgap semiconductor, comparable with GaN115. There 146 has previously been an experimental report of a cubic gallium oxynitride phase with com- position close to Ga2.8O3.5N0.5, that formed metastably during GaN thin ?lm synthesis from chemical precursors117,118. Here, we report our ?rst-principles theoretical study of the formation energetics, stability, and electronic properties of the Ga3O3N spinel-structured phase, combined with experiments using a combination of high pressure-high tempera- ture techniques to establish the formation and stability of spinel-structured Ga3O3N from Ga2O3+GaN mixtures, and to determine the chemical composition, structure and proper- ties of the resulting materials. 6.2.2 Total Energy Calculations of GaN and Ga3O3N In order to obtain the Gibbs free energy of formation of gallium oxynitride, Gibbs free energies of Ga2O3, GaN and Ga3O3N are required. As an approximation, assuming a can- celation e?ect on the vibrational entropies between Ga2O3+GaN and Ga3O3N, the static formation enthalpy ?H (P) is adopted as an estimation of ?G(T,P). Enthalpies of Ga2O3 polymorphs as a function of pressure has been shown in section 6.1.2. For GaN, previous studies showed that the ground wurtzite structure undergoes a phase transition into the rocksalt structure at high pressure. The experimentally observed transition pressure is at 37-52 GPa230?232 and the theoretically calculated Pt is 35-52 GPa43,233?235. In this study, the synthesis pressure of Ga3O3N (several GPa) is much lower than the WZ-to-RS Pt of GaN. Therefore the only polymorph we considered for GaN is the wurtzite phase, which is one of the starting materials (99.99% Ga2O3, containing a mixture of ?+? phases, and 99.99% pure GaN) in the experimental work236. Here, we ?rst present our total energy calculation of wurtzite-GaN and spinel-Ga3O3N. The calculation was on the basis of ?rst- principles density functional theory (DFT), implemented by the VASP code139?142. The many-electron exchange-correlation interaction was approximated within the local density approximation (LDA). For parts of the study that were associated with small energy dif- ferences, we compared the LDA results with calculations using the generalized gradient approximation (GGA). To improve numerical e?ciency, core electrons were approximated 147 Table 6.6: Third-order Birch-Murnaghan EOS parameters, zero pressure structural param- eters for wurtzite GaN Source B0 (GPa) B? a (?A) c/a u This work (LDA) 198 4.01 3.144 1.631 0.376 This work (GGA) 175 4.15 3.187 1.633 0.375 Calculation43 196 4.3 3.180 1.632 0.376 Calculation184 172 5.11 1.632 0.376 Experiment230 245 0.377 Experiment231 188 3.2 3.191 1.626 Experiment232 237 4.3 3.191 1.627 0.377 with ultrasoft pseudopotentials (US-PP)121, and only the s and p valence electrons in these elements and the semicore 3d electrons in Ga were treated explicitly. The wave functions of the valence and semi-core electrons were expanded using a planewave basis, with a kinetic energy cut-o? set at 348 eV for GaN and 396 eV for Ga3O3N. The calculated static energies at various volumes are ?tted to the 3rd-order Birch- Murnaghan EOS. General agreement of ?tting parameters, B0 and B?, and zero pressure structural parameters is achieved between our calculation with other reported values (Table 6.6). GGA calculated bulk modulus is about 12% less than that of the LDA value, which is also the case from other calculations. The experimental B0 is scattered from 188 to 245 GPa, as well. Despite of the underestimation of B0, the lattice parameter a by GGA gives a better agreement with other reported values. Results for Ga-O-N phases relevant to the synthesis of the spinel-structured Ga3O3N oxynitrides and the stability and properties of that compound are described below. 6.2.3 Theoretical Study of the Synthesis, Structure, and Stability of Ga3O3N Spinel We performed ab initio calculations of the atomic structures and energetic properties of the Ga3O3N systems to understand the formation and thermodynamic stability of the spinel-structured phase. The ideal spinel crystal has an A3X4 stoichiometry (A and X represent cations and anions, respectively) with two molecular equivalents per primitive unit 148 cell corresponding to the fcc structure. Assuming a Ga:O:N ratio of 3:3:1 that corresponds to the ideal stoichiometry , a solid of n Ga3O3N molecular units contains n tetrahedrally coordinated Ga atoms (labeled as IV Ga), 2n octahedrally coordinated Ga atoms (V IGa), 3n O atoms, and n N atoms. The Gibbs free energy of formation of the oxynitride from a mixture containing the corresponding oxide and nitride is de?ned as ?Gformation (T,P) = GGa3O3N (T,P) ? [GGa2O3 (T,P) + GGaN (T,P)] (6.2) Here G represents the Gibbs free energies for each molecular unit, and a negative ?Gformation corresponds to a driving force for the formation of the oxynitride. Because the oxynitride spinel can have O or N atoms distributed among the anion sites, it is likely that there will be a large con?gurational entropy term contained within ?Gformation. A complete statistical modeling of the Gibbs free energy of Ga3O3N material that might contain such oxygen/nitrogen compositional disorder requires calculating a very large number of atomic con?gurations using supercell models. In the ?rst stage of our theoretical investigation, we performed the ab initio energetic calculations with a limited number of atomic con?gura- tions using (pseudo) face-centered-cubic unitcell models. We interpreted our data using a simpli?ed model that breaks the Gibbs energy into two terms: GGa3O3N (T,P) = Hground (P) + Galloy (T,P) (6.3) Here Hground is the temperature-independent enthalpy of the ground-state (lowest energy) con?guration, which allows us to gain insights into the energetically favored local coordina- tion states and their O/N ordering, and to study pressure e?ects on the Gibbs free energy of formation. The temperature e?ects are then described by the second term, Galloy, which models the contributions related to the O/N disorder. In this study, we ignored the entropy contributions due to the lattice vibrations. We ?rst studied the pressure e?ects on formation of Ga3O3N by setting T = 0 K; the Gibbs free energy of formation ?Gformation is thus identical to the enthalpy of formation 149 Table 6.7: Parameters of the third-order Birch-Murnaghan equation of states of three unit- cell-based atomic models of the spinel-structured Ga3O3N that have the lowest energy calculated within the LDA. E0 V0 a0 B0 B? model (eV/Ga3O3N) (?A3/Ga3O3N) (?A) (GPa) I -48.109 (0.000) 69.5424 8.2246 210 4.13 II -47.966 (0.143) 69.6928 8.2357 208 4.17 III -47.854 (0.255) 69.6181 8.2275 209 4.14 ?H (P) = ?E + P?V . We ?rst carried out a search for the lowest energy unit-cell con- ?guration of the spinel-structured Ga3O3N. Among the 8!/(6!2!) = 28 con?gurations, there are three crystallographically distinct O/N arrangements for (Ga3O3N)2. We used the LDA to calculate lattice parameters at P = 1 atm and also the V (P) relations for these three models: the calculated third-order Birch-Murnaghan equation of state parameters for each of these three unit-cell models are listed in Table 6.7. The con?guration labeled as model I has a rhombohedral symmetry (R?3m). This has the lowest equilibrium energy, and it is considered to provide the ground-state con?guration of Ga3O3N in our theoretical study. The rhombohedral distortion from the ideal cubic structure is very small; i.e., the angles between the pseudocubic lattice vectors are 89.33? within the equilibrium con?guration. The calculated equilibrium volume and the bulk modulus of this ground-state con?guration compare favorably with the experimental measurements (see below). At zero pressure, our LDA calculations predict an endothermic enthalpy of formation, with ?H = +272 meV per Ga3O3N unit (26.2 kJ/mol) with respect to Ga2O3 (monoclinic ? phase) + GaN (hexag- onal wurtzite structure). The LDA calculations predict a negative d?H/dP slope (solid line in Figure 6.4a), so the Gibbs free energy of formation becomes negative at P > 17 GPa. However, Ga2O3 does not always remains in the monoclinic ? phase at high pressure. According to our (static) LDA calculations, a ?-to-? phase transition in Ga2O3 takes place at P ? 0.5 GPa. More importantly, the slope of d?H/dP becomes positive when Ga2O3 is in the ? phase (dashed line plot, Figure 6.4a). The experimental phase transition has been reported to occur between 0.1 MPa and 4.4 GPa237,238; however, a direct transition pressure was not recorded in those studies213. Recently, Tu et al. have reported a ?-to-? 150 Figure 6.4: LDA-calculated formation enthalpy ?H (i.e., Gibbs free energy of formation ?Gformation at zero temperature) as a function of pressure using (a) LDA methods and (b) the GGA approach. The calculated positive ?H suggests an endothermic formation. The solid plots are calculated assuming the oxynitrides are synthesized from ?-Ga2O3 and wurtzite GaN, and the dashed plots are calculated assuming the oxynitrides are synthesized from ?-Ga2O3 and wurtzite GaN. The solid plot and the dashed plot cross at the pressure of the ?-to-? phase transition in Ga2O3 (predicted to be 0.5 and 6.6 GPa by LDA and GGA methods, respectively). phase transition in Ga2O3 at 13.3 GPa212. However, there is some uncertainty in the nature of the starting material used in that work, and Machon et al. repeated the study: he found that the ?-to-? transition occurs at P = 22 GPa213, without heating. It is obvious that the nature of the Ga2O3 starting material, and how it transforms under high-P, T conditions, will a?ect the energetics of the high-pressure synthesis experiment. Within the present study, we repeated our calculations within the GGA (Figure 6.4b). The GGA calculations predict a phase transition pressure for Ga2O3 of 6.6 GPa. The GGA formation enthalpy for the Ga3O3N synthesis reaction is predicted as +370 meV (35.6 kJ/mol) at P = 1 atm, and this value is reduced to its minimum of +245 meV (23.6 kJ/mol) at the Ga2O3 transition pressure. Despite the quantitative di?erences between the two sets of calculations, the GGA results agree qualitatively with the ?ndings from the LDA calculations. Both the LDA and GGA studies predict that the optimal pressure for 151 synthesizing the spinel-structured Ga3O3N from Ga2O3 and GaN mixtures is around that of the ?-to-? phase transition in Ga2O3 (P ? 6.6 GPa, according to GGA calculations). Next, we investigated the temperature e?ects on the stability of the oxynitrides systems. A simple ideal solution model has been previously adopted to estimate the e?ects of O/N disorder in spinel-structured oxynitrides116,239. Such a simple statistical model is valid only in the cases where all the atomic con?gurations have very similar energies, and it approximates the additional Gibbs free energy term with a contribution due to the alloy disorder with a pure entropic term: Galloy (T,P) = ?4kBT [xlnx + (1 ?x) ln (1 ?x)] (6.4) (here the factor 4 is due to the presence of four possible di?erent anion sites in the molecular unit that are assumed to be equally accessible). It is known that such simple models can signi?cantly underestimate the alloy formation temperature if some of the atomic con?g- urations are energetically inaccessible at the temperature of the experiment. To obtain a better estimate of the entropic contribution to the Gibbs free energy that is more relevant to the experimental results, we studied the correlation between the energetic properties of the oxynitride materials and their local atomic ordering schemes among the three unit-cell models (Table 6.8). The spinel structure is described as a packing of AX4 tetrahedra and AX6 octahedra present in a 1:2 ratio. In the case of Ga3O3N, there are ?ve possible types of AX4 tetrahedra, i.e., IV GaO4, IV GaO3N, IV GaO2N2, IV GaON3, and IV GaN4. Similarly, there exist seven types of AX6 octahedra: V IGaO6, V IGaO5N, V IGaO4N2, V IGaO3N3, V IGaO2N4, VIGaON5, and VIGaN6. The distribution of various types of AX4 and AX6 units within the three models studied is listed in Table 6.8. No IV GaON3, IV GaN4, V IGaON5, or V IGaN6 species were considered to simplify the statistical analysis; such N-rich AX4 or AX6 species are expected to have low concentration because the average O:N ratio is 3:1. The LDA calculations show a clear energetic preference for the structure containing IV GaO3N sites over those with the combination 50% IV GaO4 + 50% IV GaO2N2 152 Table 6.8: Analysis of local coordination ordering schemes within the three unit-cell models of the spinel-structured Ga3O3N in terms of the ratio of various types of AX4 tetrahedral and AX6 octahedral units. tetrahedral Ga (%) octahedral Ga (%) model IV GaO4 IV GaO3N IV GaO2N2 VIGaO6 V IGaO5N V IGaO4N2 I 100 25 75 II 100 50 50 III 50 50 50 50 and the 25% V IGaO6 +75% VIGaO4N2 combination over that with 50% V IGaO5N + 50% V IGaO4N2. We then constructed a three-energy-level model to investigate the consequences of this anion site ordering on the formation energetics of the oxynitride spinel, using the ground- state energies of the three lowest energy models found above (Table 6.7): ?Gformation (T,P) = ?H (P) ? 12kBT ln bracketleftBig 4 + 12e?2??1/kBT + 12e?2??2/kBT bracketrightBig (6.5) Here ?H (P) is the static formation enthalpy as discussed above, the values 4, 12, and 12 are the degeneracies of the levels corresponding to each model, and ??1 and ??2 are the energies of the two ?excited? states (i.e., models II and III) compared to the ground-state energy (as listed in column 2 of Table 6.7). At a given pressure, the alloy formation tem- perature corresponds to that at which ?Gformation becomes zero. On the basis of our LDA calculations, this condition occurs at Talloy = 2800 K. This result means that, according to our model, a stoichiometric Ga3O3N spinel phase would become thermodynamically stabi- lized with respect to other ordering schemes and could be synthesized from Ga2O3 + GaN mixtures above P ? 6?7 GPa and T = 2800 K. This temperature estimate is considerably higher than that used experimentally to synthesize an oxynitride spinel-structured material from the component oxides and nitrides (synthesis temperatures as low as 1200 ?C at 5 GPa; see below). The main reasons for the discrepancy are that the oxynitride materials obtained experimentally contain vacancies on the Ga3+ and perhaps also on the anion sites and the O:N ratio obtained is larger than the ideal stoichiometry that was modeled. Both 153 considerations will have a large e?ect on the relative energies of ground-state and exper- imentally accessible ?excited-state? models. Also, it is not yet clear that it is justi?ed to exclude high-N-content environments such as V GaON3, IV GaN4, V IGaON5, or V IGaN6 on purely statistical grounds from the thermodynamic treatment. Such species could have special stability due to local bonding environments, including bond valence constraints128. That will have to be tested in future investigations of local coordinations in GaxOyNz mate- rials using appropriate experimental probes of the local structural environments (e.g., X-ray absorption spectroscopy/EXAFS, NMR, etc.). 6.2.4 Electronic Properties Using ab initio LDA methods, we calculated the electronic band structure within the LDA on the basis of our atomic model of the ground-state con?guration of Ga3O3N (model I). The LDA-predicted band dispersion is plotted in Figure 6.5a from the center of the Bril- louin zone (the ? point) along three directions to the F, T, and L points at the zone bound- aries. The electronic density of states function is shown in (b). Our calculations indicate that spinel-structured gallium oxynitrides are direct wide band gap semiconductors, with optoelectronic properties that are similar to those of wurtzite- and sphalerite-structured GaN and gallium oxides. Because of well-known limitations of the LDA for such electronic structure calculations, the predicted magnitude of the band gap (2.1 eV) is likely to be underestimated in this study. The same LDA methods underestimate the band gaps of the GaN (wurtzite) and ?-Ga2O3 phases by 1.3 and 2.3 eV, respectively; we thus expect the experimental band gap of Ga3O3N (spinel) to lie around 4 eV. Our collaborator obtained room-temperature photoluminescence spectra of the Ga2.8N0.64O3.24 sample obtained by high-P,T synthesis, using 325 nm laser excitation (Figure 6.6). The onset of the photolu- minescence signal begins just below 2.5 and extends to 1.5 eV. Because the experimentally synthesized material contains a large quantity of defects on the Ga3+ sites, and also per- haps on the anion sites, along with O/N disorder, it is unlikely that the photoluminescence feature corresponds to excitations across the band gap. Instead, the observed PL band is 154 Figure 6.5: (a) Electronic band structure and (b) Electronic density of states for Ga3O3N calculated using ?rst-principles (DFT) methods within the LDA. likely to arise from defect related transitions between states mainly within the gap, so that the intrinsic band gap for a stoichiometric ordered material would lie considerably above 2.5 eV, as predicted by theory. 6.2.5 Phonon Spectrum The Raman spectrum of the new oxynitride spinel phase contains several broad bands (Figure 6.7), indicating substantial disorder among the O and N atoms on the 32e sites in the spinel structure and/or the presence of cation (Ga3+) or anion vacancies. The broad bands have maxima near 700 and 800 cm?1, and also near 300 cm?1, that correspond generally to the positions of Raman-active modes within the analogous spinel-structured compound ?-Ge3N4240. To aid in the interpretation of our experimentally obtained Raman spectra of GaxOyNz phases, we calculated the ? point phonon frequencies of the 14-atom rhombo- hedral unit cell of Ga3O3N (model I structure) using ?rst-principles LDA methods. The same technique was previously used to predict the positions of the Raman-active vibrational modes in ?-Ge3N4114,240. The pseudocubic R?3m model structure for Ga3O3N is predicted to have nine Raman-active modes (4A1g + 5Eg), with zone-center frequencies calculated 155 Figure 6.6: Photoluminescence spectrum of the Ga2.8N0.64O3.24 sample obtained via high- P,T synthesis in the multianvil experiments, using UV laser excitation (325 nm) at 213, 219, 367, 379, 499, 512, 634, 647, and 782 cm?1. These calculated frequencies are denoted by dashed lines in Figure 6.7. Within an anion-disordered structure, as expected for the real GaxOyNz spinel, he predicted zone-center modes act as poles for interpreting the broadened spectra that approach the full vibrational density of states (VDOS). The expected mode frequencies are grouped around the ?ve frequency values (i.e., 216, 373, 506, 640, and 782 cm?1) that are associated with the ideal spinel structure. By analogy with the Raman spectrum of ?-Ge3N4, we expect the lowest frequency peak (216 cm?1) and the two highest frequency peaks (640 and 782 cm?1) to have the strongest intensities, whereas the two intermediate frequency peaks (373 and 506 cm?1) are relatively weak. We cannot yet directly calculate the e?ects of the O/N disorder on the broadening patterns observed in our Raman spectra. However, it is likely that the observed Raman spectrum provides a ?rst view of the VDOS functions of Ga3O3N and also ?-Ge3N4 spinels. To provide a semiquantitative estimation of the widths of the broad peaks, we cal- culated the full phonon dispersion using a 112-atom supercell model. The Born e?ective charge induced LO-TO splitting in the ionic compounds are corrected on the basis of the interplanar force constant model proposed by Kunc and Martin156. As shown in Figure 6.7, 156 Figure 6.7: (a) Raman spectrum collected for the Ga2.8N0.64O3.24 sample at an excitation wavelength of 514.5 nm. The bold solid lines on the frequency scale below indicate the positions of the Raman bands for the analogous spinel form of ?-Ge3N4. (b) Phonon density of states (VDOS) calculated for the R?3m pseudocubic Ga3O3N phase, predicted as ?model I? in the enthalpy calculations (Table 6.7 ). The dashed lines indicate the Raman-active modes for that phase. the highest frequency strong Raman peak (i.e., the Ag mode near 782 cm?1) is expected to be the sharpest one. The other strong peak near 640 cm?1 is expected to have a broader width (40-60 cm?1). 6.2.6 Conclusions Our theoretical study showed that the most stable structure for Ga3O3N corresponds to a rhombohedral distortion of the ideal spinel structure. The formation of Ga3O3N is endothermic at ambient pressure and low temperature, and the optimal synthesis pressure is predicted to lie close to that for the ?-to-? phase transition in Ga2O3 (around 6.6 GPa according to our GGA calculations). The calculated direct band gap energy for a stoichio- metric oxynitride spinel was estimated to be around 4 eV. This value is larger than that obtained from photoluminescence data collected on the experimentally synthesized sample, which likely contains Ga3+ vacancies and other structural defects. The synthesis of this new Ga-O-N phase makes contact with the important optoelectronic materials known to exist 157 in the (Ga,Al,In)N system that provide light-emitting diodes and solid-state lasers in the blue to UV range. A cubic Ga3O3N material similar to the compound synthesized here has recently been prepared in thin ?lm form via chemical precursor techniques241. That result indicates that the new materials could be developed for use within novel optoelectronic devices. 158 Chapter 7 CONCLUSIONS AND FUTURE WORK In this dissertation, I have adopted and further developed a series of ?rst-principles computational techniques to theoretically investigate the pressure-induced phase transitions and thermodynamic properties for several representative main-group oxides and nitrides. For Al2O3, we have systematically investigated the pressure-induced phase transforma- tions. The sequence of transitions under compression, i.e., corundum?Rh2O3(II)? pPV, and the T-P phase diagram we obtained are consistent with previous theoretical and ex- perimental studies. Results using US-PP and PAW are presented. By ?nding a transition path that links the corundum and Rh2O3(II) phases with monoclinic intermediate struc- tures which has a space group P2/c, we proposed a single-bond breaking-and-reforming (SB-BAR) mechanism to describe the transformation between corundum and Rh2O3(II) phase. Total energy calculation using PAW method shows that the enthalpy barrier height of the forward corundum-to-Rh2O3(II) transition is around 130 meV/atom and remains unchanged with varied pressures. However, the barrier height of the backward Rh2O3(II)- to-corundum transition decreases signi?cantly under decompression, which indicates that the Rh2O3(II) phase may not be quenchable to ambient conditions. In addition, great simi- larity is found in the elastic constants between corundum and Rh2O3(II) phase, except C33. The larger C33 of Rh2O3(II) phase means that its c axis is less compressible than that of corundum. Zero pressure ?-point phonon frequencies and their pressure dependencies are calculated for three stable polymorphs. Comparison with measured data of corundum and theoretical Raman frequencies of pPV phase suggest that, within LDA, predictions using PAW method have better results than US-PP. This conclusion is also con?rmed from the study of elastic properties. We have calculated the elastic constants and their pressure dependencies for corundum and two high-pressure phases. PAW calculated Cij yield better 159 agreement with available experimental data of ?-Al2O3. The temperature dependencies of thermal expansion coe?cient from low to high pressures have been predicted for ?-, Rh2O3(II)-, and pPV-Al2O3. For corundum at zero pressure, our US-PP data is in good agreement with Amatuni et al.?s result measured from 300 K to 2000 K, and with Aldebert et al.?s result at temperatures above 1000 K. Our PAW data lies in the middle of the ex- perimental data and agrees with Schauer?s data below 700 K and Wachtman et al.?s data above 1200 K. Our calculated heat capacity CP, entropy and adiabatic bulk modulus of corundum phase also agree well with measured results. For AlN, we presented our ab initio calculation of activation barriers of the B4-to-B1 transition for ?ve TPs proposed by Stokes et al. We showed that the ?ve bond-preserving paths can be interpreted as transformation of di?erent long-range patterns of the ?transition units? (two di?erent orientations). The transformation of ?transition unit? is equivalent to the path along TP1 (with Cmc21 symmetry). Our calculated kinetic barriers are comparable for all ?ve paths at pressures from 0 GPa to 30 GPa, which indicate that the wurtzite-to- rocksalt transition is characterized by the transformation of the ?transition unit?, while the long-range pattern is less important. And the di?erence in strains of di?erent TPs is not a major factor for at least the transition from wurtzite to rocksalt phase in AlN. In addition, the bond-breaking path is not energetically favored compared with the bond-preserving paths. Besides the studies of di?erent TPs, our estimated forward and backward barrier heights are consistent with experimental observation and previous calculations. For Si3N4, in summary, we have theoretically studied phase transitions in silicon nitride (Si3N4) at high pressure using a ?rst-principles density functional theory method. We ?nd that ?-Si3N4 remains as a metastable phase at temperatures up to 2000 K and pressures up to 10 GPa. The equilibrium ??? transition pressure is predicted as 7.5 GPa at 300K and it increases to 9.0 GPa at 2000K. Both ?- and ?-Si3N4 are dynamically stable at low pressure. However, two competing phonon-softening mechanisms are found in the ? phase at high pressures. At room temperature, ?-Si3N4 is predicted to undergo a ?rst-order ??P3 transition above 38.5 GPa, while ?-Si3N4 shows no signs of structural instability. 160 The predicted metastable high-pressure P3 phase is structurally related to ?-Si3N4. The enthalpy barrier height is estimated as only 67.23 meV/atom. Our LDA predicted thermal expansion coe?cient, heat capacity and bulk Gr?uneisen parameter are in good agreement with Bruls? measured results. We ?nd relatively large discrepancies between our calculation with experimental data from Reeber. And we attribute the cause of predicted negative TEC at low temperatures in ? and ?-Si3N4 to the low-frequency phonon modes that have negative mode Gr?uneisen ratios. For bulk Ga2O3, we have calculated the Raman frequencies and their mode Gr?uneisen parameters for ? and ? phases. Good agreement is achieved between our results with experiments and other calculations. We also predict the equilibrium T-P phase diagram of Ga2O3 consists of ?, ? and Rh2O3(II) phases. Our LDA calculated transition pressure from ? to ? phase is 0.3 GPa at 0 K and 1.6 GPa at 2000 K. The Clapeyron slope for this transition is positive which has the value of about +0.6 MPa/K at 1000 K. This is consistent with previous calculations. However, the experimental Pt is much larger. The large discrepancy can be attributed to a large kinetic barrier, which may be signi?cantly lowered under compression. Our calculated ?-to-Rh2O3(II) transition pressure is 40 GPa at 0 K and 37 GPa at 2000 K, which is in excellent agreement with experimental value (about 37 GPa at about 2000 K). In the case of nanowires, based on a comparison of the experimental Raman mode frequencies with our ?rst-principles calculations, we ?nd compelling evidence for growth-direction-induced internal strains in ?-Ga2O3 nanowires which signi?cantly in?uence the vibrational mode frequencies. Within the linear model approximation, the observed blue and redshifts of peak frequencies in the micro-Raman spectra of the ?-Ga2O3 nanowires with di?erent growth directions can be attributed to two small anisotropic internal strains: one compressive strain of 2% volume change, and the other tensile strain of 0.7% volume change. The overall high quality of the ?tted models to available experimental data suggests a strong correlation between the shifts in Raman mode frequencies and the growth-direction-induced internal strains in the Ga2O3 nanowires. 161 For gallium oxynitride, our theoretical study showed that the most stable structure for Ga3O3N corresponds to a rhombohedral distortion of the ideal spinel structure. The formation of Ga3O3N is endothermic at ambient pressure and low temperature, and the optimal synthesis pressure is predicted to lie close to that for the ?-to-? phase transition in Ga2O3 (around 6.6 GPa according to our GGA calculations). The calculated direct band gap energy for a stoichiometric oxynitride spinel was estimated to be around 4 eV. This value is larger than that obtained from photoluminescence data collected on the experimen- tally synthesized sample, which likely contains Ga3+ vacancies and other structural defects. The synthesis of this new Ga-O-N phase makes contact with the important optoelectronic materials known to exist in the (Ga,Al,In)N system that provide light-emitting diodes and solid-state lasers in the blue to UV range. A cubic Ga3O3N material similar to the com- pound synthesized here has recently been prepared in thin ?lm form via chemical precursor techniques. That result indicates that the new materials could be developed for use within novel optoelectronic devices. The above results show that DFT calculations within the frame of quasi-harmonic approximation (QHA) is successful in calculating the ?nite-temperature thermodynamic potentials for hard materials. However, to improve the prediction of thermodynamic prop- erties for normal materials at high temperatures or strongly anharmonic crystals even at low temperatures, we are looking forward to developing and implementing an algorithm based on the perturbation theory to add anharmonic correction to the QHA vibrational free energy, as shown in Appendix C. The ?rst order correction is contributed from the 3rd and 4th order lattice anharmonicity tensors. With these anharmonicity tensors, the temperature dependence of elasticity and phonon frequencies are also within the scope. And, for the metastable P3 phase we found in Si3N4, our collaborator who did ex- periments showed that there is a mismatch between the observed X-ray di?raction (XRD) pattern after phase transition at high pressure and the calculated XRD pattern from P3 structure at the same pressure. But the unknown ? phase is believed to be structurally close 162 to the P3 phase. We will look at the low-frequency phonon modes with softening tendency for ?, P?6, P3 and P21/m phases, which may suggest other possible structures. 163 Bibliography 1. H. Mao, P. Bell, J. Shaner, and D. Steinberg, J. Appl. Phys. 49, 3276 (1978). 2. W. Nellis and C. Yoo, J. Geophys. Res. 95, 21749 (1990). 3. H. Cynn, D. Isaak, E. Cohen Ronald, M. Nicol, and O. Anderson, Am. Mineral. 75, 439 (1990). 4. F. Marton and R. Cohen, Am. Mineral. 79, 789 (1994). 5. K. Thomson, R. Wentzcovitch, and M. Bukowinski, Science 274, 1880 (1996). 6. A. Oganov and S. Ono, PNAS 102, 10828 (2005). 7. R. Caracas and R. Cohen, Geophys. Res. Lett 32, L06303 (2005). 8. N. Funamori and R. Jeanloz, Science 278, 1109 (1997). 9. T. Mashimo, K. Tsumoto, K. Nakamura, Y. Noguchi, K. Fukuoka, and Y. Syono, Geophys. Res. Lett 27, 2021 (2000). 10. J. Lin, O. Degtyareva, C. Prewitt, P. Dera, N. Sata, E. Gregoryanz, H. Mao, and R. Hemley, Nat. Mater. 3, 389 (2004). 11. S. Ono, A. Oganov, T. Koyama, and H. Shimizu, Earth Planet. Sci. Lett. 246, 326 (2006). 12. S. Jahn, P. Madden, and M. Wilson, Phys. Rev. B 69, 20106 (2004). 13. K. Umemoto and R. Wentzcovitch, PNAS 105, 6526 (2008). 14. A. Jephcoat, R. Hemley, H. Mao, and K. Goettel, Physica B 150, 115 (1988). 164 15. J. Mougin, T. Le Bihan, and G. Lucazeau, J. Phys. Chem. Solids. 62, 553 (2001). 16. J. Wachtman, T. Scuderi, and G. Cleek, J. Am. Ceram. Soc. 45, 319 (1962). 17. A. Schauer, Can. J. Phys. 43, 523 (1965). 18. A. Amatuni, T. Malyutina, V. Chekhovskoi, and V. Petukhov, High Temp-High Press 8, 565 (1976). 19. G. White and R. Roberts, High Temp-High Press 15, 321 (1983). 20. P. Aldebert and J. Traverse, High Temp-High Press 16, 127 (1984). 21. G. Fiquet, P. Richet, and G. Montagnac, Phys. Chem. Miner. 27, 103 (1999). 22. G. Furukawa, T. Douglas, R. McCoskey, and D. Ginnings, J. Res. Nat. Bur. Stand. 57, 67 (1956). 23. W. Te?t, Natl Bur Stand 70, 277 (1966). 24. D. Chung and G. Simmons, J. Appl. Phys. 39, 5316 (1968). 25. T. Goto, O. ANDERSON, I. Ohno, and S. Yamamoto, J. Geophys. Res. 94, 7588 (1989). 26. J. Hama and K. Suito, Phys. Chem. Miner. 28, 258 (2001). 27. W. Mayer and E. Hiedemann, J. Acoust. Soc. Am. 30, 756 (1958). 28. J. Wachtman, W. Te?t, D. Lam, and R. Stinch?eld, J. Am. Ceram. Soc. 43, 334 (1960). 29. J. Gieske and G. Barsch, Phys. Stat. sol. (b) 29, 121 (1968). 30. R. Cohen, Geophys. Res. Lett. 14, 37 (1987). 31. W. Duan, B. Karki, and R. Wentzcovitch, Am. Mineral. 84, 1961 (1999). 165 32. Y. Le Page, P. Saxe, and J. Rodgers, Phys. Status Solidi C 229, 1155 (2002). 33. J. Gladden, J. So, J. Maynard, P. Saxe, and Y. Le Page, Appl. Phys. Lett. 85, 392 (2004). 34. D. Hovis, A. Reddy, and A. Heuer, Appl. Phys. Lett. 88, 13 (2006). 35. S. Shang, Y. Wang, and Z. Liu, Appl. Phys. Lett 90, 101909 (2007). 36. S. Stackhouse, J. Brodholt, and G. Price, Geophys. Res. Lett. 32, 13305 (2005). 37. H. Vollst?aDt, E. Ito, M. Akaishi, S. Akimoto, and O. Fukunaga, Proc. Jpn. Acad., B 66, 7 (1990). 38. M. Ueno, A. Onodera, O. Shimomura, and K. Takemura, Phys. Rev. B 45, 10123 (1992). 39. Q. Xia, H. Xia, and A. Ruo?, J. Appl. Phys. 73, 8198 (1993). 40. S. Uehara, T. Masamoto, A. Onodera, M. Ueno, O. Shimomura, and K. Takemura, J. Phys. Chem. Solids 58, 2093 (1997). 41. P. E. Van Camp, V. E. Van Doren, and J. T. Devreese, Phys. Rev. B 44, 9056 (1991). 42. N. E. Christensen and I. Gorczyca, Phys. Rev. B 50, 4397 (1994). 43. J. Serrano, A. Rubio, E. Hern?andez, A. Mu?noz, and A. Mujica, Phys. Rev. B 62, 16612 (2000). 44. A. Saitta and F. Decremps, Phys. Rev. B 70, 35214 (2004). 45. C. Silva, H. Leite Alves, L. Scolfaro, and J. Leite, Phys. Status Solidi C 2, 2468 (2005). 46. J. Cai and N. Chen, Phys. Rev. B 75, 134109 (2007). 47. H. Stokes and D. Hatch, Isotropy Software Package (2004), URL http://stokes.byu.edu/isotropy.html. 166 48. H. Stokes, J. Gunter, D. Hatch, J. Dong, H. Wang, and J. Lewis, Phys. Rev. B 76, 012102 (2007). 49. D. Zahn, Y. Grin, and S. Leoni, Phys. Rev. B 72, 64110 (2005). 50. H. Sowa, Acta Crystallogr., Sect. A: Found. Crystallogr. 61, 325 (2005). 51. C. Capillas, J. Perez-Mato, and M. Aroyo, J. Phys.: Condens. Matter 19, 275203 (2007). 52. M. Aroyo, J. Perez-Mato, C. Capillas, E. Kroumova, S. Ivantchev, G. Madariaga, A. Kirov, and H. Wondratschek, Z. Kristallogr. 221, 15 (2006). 53. H. Sowa, Acta Crystallogr., Sect. A: Found. Crystallogr. 57, 176 (2001). 54. S. Limpijumnong and W. Lambrecht, Phys. Rev. Lett. 86, 91 (2001). 55. F. Shimojo, S. Kodiyalam, I. Ebbsj?o, R. Kalia, A. Nakano, and P. Vashishta, Phys. Rev. B 70, 184111 (2004). 56. K. Komeya and M. Matsui, Materials Science and Technology, vol. 11 (Wiley-VCH, Weinheim, 1994). 57. F. Schr?oder, Gmelin Handbook of Inorganic and Organometallic Chemistry, vol. Si Suppl. B 5c of Silicon Nitride in Electronics (Springer-Verlag, Berlin, 1991). 58. A. Zerr, G. Miehe, G. Serghiou, M. Schwarz, E. Kroke, R. Riedel, H. Fue?, P. Kroll, and R. Boehler, Nature 400, 340 (1999). 59. J. Liang, L. Topor, A. Navrotsky, and M. Mitomo, J. Mater. Res 14, 1959 (1999). 60. C. Greskovich and S. Prochazka, J. Am. Ceram. Soc. 60, 471 (1977). 61. L. Bowen, R. Weston, T. Carruthers, and R. Brook, J. Mater. Sci. 13, 341 (1978). 62. M. Shimada, N. Ogawa, M. Koizumi, F. Dachille, and R. Roy, Am. Ceram. Soc. Bull. 58, 519 (1979). 167 63. L. Gauckler, H. Hohnke, and T. Tien, J. Am. Ceram. Soc. 63, 35 (1980). 64. K. Jack, Progress in Nitrogen Ceramics, Martinus Nijho?, The Hague, Netherlands pp. 45?60 (1983). 65. J. Park, J. Kim, and C. Kim, J. Am. Ceram. Soc. 70, 240 (1987). 66. H. Hirai and K. Kondo, J. Am. Ceram. Soc. 77, 487 (1994). 67. H. Suematsu, M. Mitomo, T. Mitchell, J. Petrovic, O. Fukunaga, and N. Ohashi, J. Am. Ceram. Soc. 80, 615 (1997). 68. R. Gr?un, Structural Crystallography and Crystal Chemistry 567, 7408 (1979). 69. W. Y. Ching, L. Ouyang, and J. D. Gale, Phys. Rev. B 61, 8696 (2000). 70. J. Wendel and W. Goddard, J. Chem. Phys. 97, 5048 (1992). 71. A. Kuwabara, K. Matsunaga, and I. Tanaka, Phys. Rev. B 78, 064104 (2008). 72. M. Schwarz, G. Miehe, A. Zerr, E. Kroke, B. Poe, H. Fuess, and R. Riedel, Adv. Mater. 12, 883 (2000). 73. T. Sekine, H. He, T. Kobayashi, M. Zhang, and F. Xu, Appl. Phys. Lett. 76, 3706 (2000). 74. E. Soignard, M. Somayazulu, J. Dong, O. Sankey, and P. McMillan, J. Phys: Condens. Matter 13, 557 (2001). 75. Y. Zhang, A. Navrotsky, and T. Sekine, J. Mar. Res. 21, 41 (2006). 76. J. Jiang, F. Kragh, D. Frost, K. Stahl, and H. Lindelov, J. Phys: Condens. Matter 13, 515 (2001). 77. T. Sekine and T. Mitsuhashi, Appl. Phys. Lett. 79, 2719 (2001). 78. A. Zerr, Phys. Stat. sol. (b) 227, R4 (2001). 168 79. P. Kroll, J. Solid State Chem. 176, 530 (2003). 80. P. McMillan, O. Shebanova, D. Daisenberger, R. Cabrera, E. Bailey, A. Hector, V. Lees, D. Machon, A. Sella, and M. Wilson, Phase Transitions 10, 1003 (2007). 81. E. Soignard, P. McMillan, C. Hejny, and K. Leinenweber, J. Solid State Chem. 177, 299 (2004). 82. J. Dong, O. Sankey, S. Deb, G. Wolf, and P. McMillan, Phys. Rev. B 61, 11979 (2000). 83. C. M. B. Henderson and D. Taylor, Trans. J. Brit. Ceramic Soc. 74, 49 (1975). 84. J. Schneider, F. Frey, N. Johnson, and K. Laschke, Z. Kristallogr. 209, 328 (1994). 85. R. J. Bruls, H. T. Hintzen, G. de With, R. Metselaar, and J. C. van Miltenburg, J. Phys. Chem. Solids. 62, 783 (2001). 86. R. R. Reeber, Therm. Conduct. 27, 525 (2005). 87. M. Kuriyama, Am. Ceram. Soc. Bull. 57, 1119 (1978). 88. V. Koshchenko and G. Ya, Inorg. Mater. Inorg. Mater. 18, 903 (1982). 89. H. T. Hintzen, M. R. M. M. Hendrix, H. Wondergem, C. M. Fang, T. Sekine, and G. de With, J. Alloys Compd. 351, 40 (2003). 90. W. Paszkowicz, R. Minikayev, P. Piszora, M. Knapp, C. Bahtz, J. M. Recio, M. Mar- ques, P. Mori-Sanchez, L. Gerward, and J. Z. Jiang, Phys. Rev. B 69, 52103 (2004). 91. C. Fang, G. de Wijs, H. Hintzen, et al., J. Appl. Phys. 93, 5175 (2003). 92. L. Wang, J. Sun, W. Yang, and R. Tian, Acta Phys. Pol., A 114, 807 (2008). 93. M. Fleischer and H. Meixner, J. Mater. Sci. Lett. 11, 1728 (1992). 94. R. Roy, V. Hill, and E. Osborn, J. Am. Chem. Soc. 74, 719 (1952). 95. J. Remeika and M. Marezio, Appl. Phys. Lett. 8, 87 (1966). 169 96. K. Lipinska-Kalita, B. Chen, M. Kruger, Y. Ohki, J. Murowchick, and E. Gogol, Phys. Rev. B 68, 035209 (2003). 97. M. Grimsditch, Phys. Rev. Lett. 52, 2379 (1984). 98. C. Meade, R. Hemley, and H. Mao, Phys. Rev. Lett. 69, 1387 (1992). 99. C. Liang, G. Meng, G. Wang, Y. Wang, L. Zhang, and S. Zhang, Appl. Phys. Lett. 78, 3202 (2001). 100. G. Gundiah, A. Govindaraj, and C. Rao, Chem. Phys. Lett. 351, 189 (2002). 101. Y. Choi, W. Kim, Y. Park, S. Lee, D. Bae, Y. Lee, G. Park, W. Choi, N. Lee, and J. Kim, Adv. Mater. 12, 746 (2000). 102. Y. Gao, Y. Bando, T. Sato, Y. Zhang, and X. Gao, Appl. Phys. Lett. 81, 2267 (2002). 103. S. Sharma and M. K. Sunkara, J. Am. Chem. Soc. 124, 1228812293 (2002). 104. R. Rao, A. Rao, B. Xu, J. Dong, S. Sharma, and M. Sunkara, J. Appl. Phys. 98, 094312 (2005). 105. M. Haneda, Y. Kintaichi, H. Shimada, and H. Hamada, Chem. Lett. 27, 181 (1998). 106. E. Aubay and D. Gourier, J. Phys. Chem. 96, 5513 (1992). 107. N. Corbin, J. Eur. Ceram. Soc. 5, 143 (1989). 108. H. Willems, G. De With, R. Metselaar, R. Helmholdt, and K. Petersen, J. Mater. Sci. Lett. 12, 1470 (1993). 109. V. Dravid, J. Sutli?, A. Westwood, M. Notis, and C. Lyman, Philos. Mag. A 61, 417 (1990). 110. A. Zerr, G. Miehe, G. Serghiou, M. Schwarz, E. Kroke, and R. Riedel, Nature 400, 340 (1999). 170 111. K. Leinenweber, M. O?Kee?e, M. Somayazulu, H. Hubert, P. McMillan, and G. Wolf, Chem. Eur. J. 5, 3076 (1999). 112. G. Serghiou, G. Miehe, O. Tschauner, A. Zerr, and R. Boehler, J. Chem. Phys. 111, 4659 (1999). 113. H. He, T. Sekine, T. Kobayashi, and K. Kimoto, J. Appl. Phys. 90, 4403 (2001). 114. E. Soignard, P. McMillan, and K. Leinenweber, Chem. Mater 16, 5344 (2004). 115. J. Lowther, T. Wagner, I. Kinski, and R. Riedel, J. Alloys Compd. 376, 1 (2004). 116. P. Kroll, R. Dronskowski, and M. Martin, J. Mater. Chem. 15, 3296 (2005). 117. S. Wolter, J. DeLucca, S. Mohney, R. Kern, and C. Kuo, Thin Solid Films 371, 153 (2000). 118. M. Puchinger, D. Kisailus, F. Lange, and T. Wagner, J. Cryst. Growth 245, 219 (2002). 119. P. Hohenberg and W. Kohn, Phys. Rev 136, B864 (1964). 120. W. Kohn, L. Sham, et al., Phys. Rev 140, A1133 (1965). 121. D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). 122. W. Frank, C. Els?asser, and M. F?ahnle, Phys. Rev. Lett. 74, 1791 (1995). 123. G. Kresse, J. Furthm?uller, and J. Hafner, Europhys. Lett 32, 729 (1995). 124. A. Garc??a and D. Vanderbilt, Phys. Rev. B 54, 3817 (1996). 125. J. Dong, O. Sankey, and G. Kern, Phys. Rev. B 60, 950 (1999). 126. J. Dong and O. Sankey, J. Phys: Condens. Matter 11, 6129 (1999). 127. J. Dong, J. Tomfohr, and O. Sankey, Phys. Rev. B 61, 5967 (2000). 171 128. J. Dong, J. Deslippe, O. Sankey, E. Soignard, and P. McMillan, Phys. Rev. B 67, 94104 (2003). 129. Z. Feng, Sic Power Materials: Devices and Applications (Springer, 2004). 130. T. Hahn and C. A. J., International Tables for Crystallography, Volume A: Space-group Symmetry (Kluwer Academic Publishers, Dordrecht, 1992), 3rd ed. 131. V. Dmitriev, S. Rochal, Y. Gufan, and P. Toledano, Physical review letters 60, 1958 (1988). 132. C. Mailhiot and A. McMahan, Physical review. B, Condensed matter 44, 11578 (1991). 133. A. Christy, Acta Crystallographica Section B: Structural Science 49, 987 (1993). 134. M. O?Kee?e, B. Hyde, and W. Baur, Crystal Structures I. Patterns and Symmetry (Mineralogical Society of America Washington, DC, 1996). 135. H. Stokes, D. Hatch, J. Dong, and J. Lewis, J. Phys. Chem. Solids Phys Rev B 69, 174111 (2004). 136. D. Hatch, H. Stokes, J. Dong, J. Gunter, H. Wang, and J. Lewis, Phys. Rev. B 71, 184109 (2005). 137. M. Born and R. Oppenheimer, Ann. Phys. 389, 457 (1927). 138. L. Thomas, in Proc. Cambridge Philos. Soc (1927), vol. 23, pp. 542?548. 139. G. Kresse and J. Furthm?uller, Phys. Rev. B 54, 11169 (1996). 140. G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993). 141. G. Kresse and J. Furthm?uller, Comput. Mater. Sci 6, 15 (1996). 142. G. Kresse and J. Hafner, J. Phys.: Condens. Matter 6, 5 (1994). 143. P. Bl?ochl, Phys. Rev. B 50, 17953 (1994). 172 144. G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 145. D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). 146. J. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). 147. D. Langreth and J. Perdew, Phys. Rev. B 21, 5469 (1980). 148. D. Langreth and M. Mehl, Phys. Rev. B 28, 1809 (1983). 149. J. Perdew, Phys. Rev. B 33, 8822 (1986). 150. J. Perdew and W. Yue, Phys. Rev. B 33, 8800 (1986). 151. J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 152. R. King-Smith and R. Needs, J. Phys.: Condens. Matter 2, 3431 (1990). 153. K. Karch, P. Pavone, W. Windl, O. Sch?utt, and D. Strauch, Phys. Rev. B 50, 17054 (1994). 154. O. Sch?utt, P. Pavone, W. Windl, K. Karch, and D. Strauch, Phys. Rev. B 50, 3746 (1994). 155. M. Yin and M. Cohen, Phys. Rev. B 26, 3259 (1982). 156. K. Kunc and R. Martin, Phys. Rev. Lett. 48, 406 (1982). 157. G. Lehmann and M. Taut, Physica Status Solidi(b) 54, 469 (1972). 158. P. Bl ?ochl, O. Jepsen, and O. Andersen, Physical Review B 49, 16223 (1994). 159. F. Birch, J. Appl. Phys. 9, 279 (1938). 160. F. Birch, Phys. Rev. 71, 809 (1947). 161. P. Vinet, J. Ferrante, J. Rose, and J. Smith, J. Geophys. Res. 92, 9319 (1987). 173 162. P. Vinet, J. Rose, J. Ferrante, and J. Smith, J. Phys.: Condens. Matter 1, 1941 (1989). 163. F. Stacey, B. Brennan, and R. Irvine, Surv. Geophys. 4, 189 (1981). 164. O. Anderson, Equations of State of Solids for Geophysics and Ceramic Science (Oxford University Press, USA, 1995). 165. S. Porto and R. Krishnan, J. Chem. Phys. 47, 1009 (1967). 166. G. Watson Jr, W. Daniels, and C. Wang, J. Appl. Phys. 52, 956 (1981). 167. J. Xu, E. Huang, J. Lin, and L. Xu, Am. Mineral. 80, 1157 (1995). 168. A. Barker, Phys. Rev. 132, 1474 (1963). 169. F. Gervais and B. Piriou, J. Phys. C 7, 2374 (1974). 170. M. Schubert, T. Tiwald, and C. Herzinger, Phys. Rev. B 61, 8187 (2000). 171. S. Shin, F. Pollad, and P. Paccah, in Proceedings of the third international conference on light scattering in solids, edited by M. Balkanski, R. Leite, and S. Porto (Wiley, New York, 1976), pp. 401?405. 172. H. Schober, D. Strauch, and B. Dorner, Z. Phys. B: Condens. Matter 92, 273 (1993). 173. R. Heid, D. Strauch, and K. Bohnen, Phys. Rev. B 61, 8625 (2000). 174. Z. suppressLodziana and K. Parli?nski, Phys. Rev. B 67, 174106 (2003). 175. B. Montanari, B. Civalleri, C. Zicovich-Wilson, and R. Dovesi, Int. J. Quantum Chem. 106, 1703 (2006). 176. D. Schiferl, W. Denner, H. Schulz, and W. Holzapfel, J. Appl. Phys. 49, 4411 (1978). 177. S. Tolbert and A. Alivisatos, J. Chem. Phys. 102, 4642 (1995). 178. S. Sharma and Y. Gupta, Phys. Rev. B 58, 5964 (1998). 174 179. M. Knudson and Y. Gupta, Phys. Rev. Lett. 81, 2938 (1998). 180. M. Knudson, Y. Gupta, and A. Kunz, Phys. Rev. B 59, 11704 (1999). 181. M. Wilson and P. Madden, J. Phys: Condens. Matter 14, 4629 (2002). 182. S. Limpijumnong and S. Jungthawan, Phys. Rev. B 70, 054104 (2004). 183. Y. Cai, S. Wu, R. Xu, and J. Yu, Phys. Rev. B 73, 184104 (2006). 184. C. Stamp? and C. G. Van de Walle, Phys. Rev. B 59, 5521 (1999). 185. S. Saib and N. Bouarissa, Eur. Phys. J. B 47, 379 (2005). 186. H. Suematsu, J. Petrovic, and T. Mitchell, Mater. Res. Soc. Symp. Proc. 287, 449 (1992). 187. A. Togo and P. Kroll, J. Comput. Chem. 29, 2255 (2008). 188. P. Kroll and J. Von Appen, Phys. Stat. sol. (b) 226, 6 (2001). 189. K. Tatsumi, I. Tanaka, H. Adachi, F. Oba, and T. Sekine, J. Am. Ceram. Soc. 85, 7 (2002). 190. S. Ogata, N. Hirosaki, C. Kocer, and Y. Shibutani, Acta Mater. 52, 233 (2004). 191. M. Yashima, Y. Ando, and Y. Tabira, J. Phys. Chem. B 111, 3609 (2007). 192. H. Toraya, J. Appl. Cryst. 33, 95 (2000). 193. M. B. Kruger, J. H. Nguyen, Y. M. Li, W. A. Caldwell, M. H. Manghnani, and R. Jeanloz, Phys. Rev. B 55, 3456 (1997). 194. R. F. Zhang, S. H. Sheng, and S. Veprek, Appl. Phys. Lett. 90, 191903 (2007). 195. D. Du Boulay, N. Ishizawa, T. Atake, V. Streltsov, K. Furuya, and F. Munakata, Acta Crystallogr., Sect. B: Struct. Sci. 60, 388 (2004). 175 196. Y. Li, M. Kruger, J. Nguyen, W. Caldwell, and R. Jeanloz, Solid State Commun. 103, 107 (1997). 197. H. He, T. Sekine, T. Kobayashi, H. Hirosaki, and I. Suzuki, Phys. Rev. B 62, 11412 (2000). 198. A. Zerr, M. Kempf, M. Schwarz, E. Kroke, M. Goken, and R. Riedel, J. Am. Ceram. Soc. 85, 86 (2002). 199. J. W. Swegle, J. Appl. Phys. 68, 1563 (1990). 200. A. Zerr, G. Miehe, G. Serghiou, M. Schwarz, E. Kroke, R. Riedel, and R. Boehler, Science and Technology of High Pressure AIRAPT-17, 914 (2000). 201. G. Schmitz, P. Gassmann, and R. Franchy, J. Appl. Phys. 83, 2533 (1998). 202. M. Passlack, E. Schubert, W. Hobson, M. Hong, N. Moriya, S. Chu, K. Konstadinidis, J. Mannaerts, M. Schnoes, and G. Zydzik, J. Appl. Phys. 77, 686 (1995). 203. M. Fleischer and H. Meixner, Sens. Actuators B 4, 437 (1991). 204. E. Aubay and D. Gourier, Phys. Rev. B 47, 15023 (1993). 205. M. Passlack, N. Hunt, E. Schubert, G. Zydzik, M. Hong, J. Mannaerts, R. Opila, and R. Fischer, Appl. Phys. Lett. 64, 2715 (1994). 206. T. Xiao, A. Kitai, G. Liu, A. Nakua, and J. Barbier, Applied Physics Letters 72, 3356 (1998). 207. T. Miyata, T. Nakatani, and T. Minami, J. Lumin. 87, 1183 (2000). 208. J. Hao and M. Cocivera, J. Phys. D: Appl. Phys. 35, 433 (2002). 209. H. Kim and W. Kim, J. Appl. Phys 62, 2000 (1987). 210. T. Beales, C. Goodman, and K. Scarrott, Solid State Commun. 73, 1 (1990). 176 211. R. Zhou and R. Snyder, Acta Crystallographica Section B: Structural Science 47, 617 (1991). 212. B. Tu, Q. Cui, P. Xu, X. Wang, W. Gao, C. Wang, J. Liu, and G. Zou, J. Phys.: Condens. Matter 14, 10627 (2002). 213. D. Machon, P. McMillan, B. Xu, and J. Dong, Phys. Rev. B 73, 094125 (2006). 214. H. Yusa, T. Tsuchiya, N. Sata, and Y. Ohishi, Physical Review B 77, 064107 (2008). 215. D. Dohy, G. Lucazeau, and A. Revcolevschi, J. Solid State Chem. 45, 180 (1982). 216. M. Marezio and J. Remeika, J. Chem. Phys. 46, 1862 (1967). 217. T. Tsuchiya, H. Yusa, and J. Tsuchiya, Phys. Rev. B 76, 174108 (2007). 218. J. Rubio, H. Van Der Meulen, N. Mestres, J. Calleja, K. Wang, P. Ils, A. Forchel, N. Gippius, and S. Tikhodeev, Solid-State Electron. 40, 707 (1996). 219. A. Arora, T. Ravindran, G. Reddy, A. Sikder, and D. Misra, Diamond Relat. Mater. 10, 1477 (2001). 220. S. Zhang, B. Zhu, F. Huang, Y. Yan, E. Shang, S. Fan, and W. Han, Solid State Communications 111, 647 (1999). 221. L. Dai, X. Chen, X. Zhang, A. Jin, T. Zhou, B. Hu, and Z. Zhang, J. Appl. Phys. 92, 1062 (2002). 222. H. Seo, S. Bae, J. Park, H. Yang, K. Park, and S. Kim, J. Chem. Phys. 116, 9492 (2002). 223. H. Okumura, K. Ohta, K. Ando, W. R?uhle, T. Nagatomo, and S. Yoshida, Solid-State Electron. 41, 201 (1997). 224. H. Maruska and J. Tietjen, Appl. Phys. Lett. 15, 327 (1969). 225. S. Mohammad, A. Salvador, and H. Morkoc, Proc. IEEE 83, 1306 (1995). 177 226. H. Tippins, Phys. Rev. 140, 316 (1965). 227. S. Mo, L. Ouyang, W. Ching, I. Tanaka, Y. Koyama, and R. Riedel, Phys. Rev. Lett 83, 5046 (1999). 228. N. Scotti, W. Kockelmann, J. Senker, S. Trassel, and H. Jacobs, Z. Anorg. Allg. Chem 625, 1435 (1999). 229. M. Shemkunas, G. Wolf, K. Leinenweber, and W. Petuskey, J. Am. Ceram. Soc. 85, 101 (2002). 230. P. Perlin, C. Jauberthie-Carillon, J. Itie, A. San Miguel, I. Grzegory, and A. Polian, Phys. Rev. B 45, 83 (1992). 231. H. Xia, Q. Xia, and A. Ruo?, Phys. Rev. B 47, 12925 (1993). 232. M. Ueno, M. Yoshida, A. Onodera, O. Shimomura, and K. Takemura, Phys. Rev. B 49, 14 (1994). 233. A. Munoz and K. Kunc, Phys. Rev. B 44, 10372 (1991). 234. P. Van Camp, V. Van Doren, and J. Devreese, Solid State Commun. 81, 23 (1992). 235. R. Pandey, J. Ja?e, and N. Harrison, J. Phys. Chem. Solids 55, 1357 (1994). 236. E. Soignard, D. Machon, P. McMillan, J. Dong, B. Xu, and K. Leinenweber, Chem. Mater 17, 5465 (2005). 237. L. Foster and H. Stumpf, J. Am. Chem. Soc. 73, 1590 (1951). 238. J. Remeika and M. Marezio, Appl. Phys. Lett. 8, 87 (1966). 239. C. Fang, R. Metselaar, H. Hintzen, and G. With, J. Am. Ceram. Soc. 84, 2633 (2001). 240. S. Deb, J. Dong, H. Hubert, P. McMillan, and O. Sankey, Solid State Commun. 114, 137 (2000). 178 241. I. Kinski, G. Miehe, G. Heymann, R. Theissmann, R. Riedel, and H. Huppertz, Z. Naturforsch., B: Chem. Sci. 60, 831 (2005). 242. URL http://cms.mpi.univie.ac.at/vasp/vasp/vasp.html. 243. T. H. K. Barron and M. L. Klein, Dynamical Properties of Solids, vol. 1 (North-Holland Publishing Company American Elsevier Publishing Company, Inc., 1974). 179 Appendices 180 Appendix A Total Energy Calculation with VASP The VASP package242 that we adopt in our study is implemented with ultra-soft pseu- dopotentials (US-PP) or the projector-augmented wave (PAW). The four essential input ?les needed for a VASP calculation are INCAR, POTCAR, KPOINTS and POSCAR. The INCAR ?le contains a large number of parameters and most of these parameters should be left at their default values in our calculations. One important parameter worth mentioning here is the ISIF-tag. ISIF controls which degree of freedom (ions, cell volume, cell shape) are allowed to change. With di?erent crystal symmetries and purposes of calculations, one can choose a proper value. Table A.1 shows the meaning of each ISIF tag. Taking B1 NaCl as an example, which has cubic symmetry, if we are interested in a constant-volume total energy calculation, ISIF=0 is the choice. The POTCAR ?le contains the pseudopotential of each type of atoms. The KPOINTS ?le contains the k-point sampling scheme. And the POSCAR ?le basically includes information of cell geometry and the atomic coordinates. Table A.1: Value and meaning of ISIF tag ISIF calculate force calculate stress tensor relax ions change cell shape change cell volume 0 yes no yes no no 1 yes trace only yes no no 2 yes yes yes no no 3 yes yes yes yes yes 4 yes yes yes yes no 5 yes yes no yes no 6 yes yes no yes yes 7 yes yes no no yes Although VASP allows relaxation of cell volume, it is not recommended in an energy- volume calculation. The incompleteness of the plane wave basis set with respect to changes of volume will cause Pulay stress, which will underestimate the equilibrium volume unless a large plane wave cuto? is used. Nonetheless, if only volume conserving relaxations are carried out, the Pulay stress can often be neglected. Furthermore, errors from Pulay stress can also be reduced and reliable lattice constant, bulk modulus and other elastic properties can be yielded by doing calculations at di?erent volumes with the same energy cuto? for each calculation and then ?tting the total energies to an equation of state. But the interval of volume changes should be in the order of 5%-10% to average out the errors due to the basis set incompleteness. 181 Appendix B Elasticity Calculation A solid responds to stress by changing its shape. Elasticity is the physical property to represent the ability of a material to recover when the stress is released. ?Elastic? implies that the solid will restore to its original equilibrium shape if the stress is released. ?ij = Sijkl?kl (B.1a) ?ij = Cijkl?kl (B.1b) where ?ij and ?kl are strain and stress tensors, Sijkl and Cijkl are called the elastic com- pliance constants and elastic sti?ness constants (or simply called elastic constants, sti?ness constants). ?ij is dimensionless. ?ij has the dimension of [force]/[area] or [energy]/[volume]. Because each su?x runs from 1 to 3 (x, y, z), there are totally 81 elastic constants or com- pliance constants. Elastic constants are measurable and characteristic properties of a solid. Without considering any temperature e?ect (entropy is zero), the change of the internal energy density (energy per unit volume) purely caused by strain is a quadratic function of ?ij: E/V = 12 summationdisplay ijkl Cijkl?ij?kl (B.2) Since ?ij and ?ji always appear together, we can de?ne a symmetrical strain tensor ? as follows: ?ij = 12 (?ij + ?ji) (B.3) where ?ij = ?ji and the energy density is still in the previous form. E/V = 12 summationdisplay ijkl Cijkl?ij?kl (B.4) From the above equation it is easy to prove that Cijkl = Cijlk and Cijkl = Cjikl. These two conditions reduce the number of independent elastic constants from 81 to 36. However, it is cumbersome to deal with a fourth-rank tensor. Fortunately, the symmetry makes it possible to use a matrix form. The most commonly accepted way of translation is Voigt?s notation, by reducing the ?rst two su?xes to one and the last two su?xes in the same way. Accordingly, the stress and strain tensor (second-order) components will be represented 182 with a single su?x. ? ? ?11 ?12 ?13 ?21 ?22 ?23 ?31 ?32 ?33 ? ? ? ? ? ?1 ?6 ?5 ?6 ?2 ?4 ?5 ?4 ?3 ? ? (B.5a) ? ? ?11 ?12 ?13 ?21 ?22 ?23 ?31 ?32 ?33 ? ? ? ? ? ?1 12?6 12?5 1 2?6 ?2 1 2?41 2?5 1 2?4 ?3 ? ? (B.5b) The purpose of introducing a factor of 12 to the o?-diagonal strain components is to keep equation B.4 in the same format. Besides the su?xes abbreviation, the elastic compliance constants are changed by a factor of 1, 2, or 4 in the following way: ?? ? Sijkl = Smn if m and n are 1, 2 or 3 2Sijkl = Smn if either m or n are 4, 5 or 6 4Sijkl = Smn if both m and n are 4, 5 or 6 (B.6) While at the same time, for Cijkl there are no such factors changes. Cijukl = Cmn (B.7) where i,j,k,l = 1,2,3 and m,n = 1,??? ,6. And the energy density due to strain is written as: E/V = 12 summationdisplay ij Cij?i?j (B.8) In equation B.8, exchanging of su?x i and j does no e?ect, which leads to Cij = Cji. Thus the number of independent Cij is further reduced to 21, and this is the case for a crystal with triclinic symmetry, which is the lowest symmetry of all. Any material with higher symmetry has even less number of independent elastic constants. Results for all the crystal classes are given in table B.1. 183 Table B.1: Independent components of Cij for all the crystal classes Crystal class Independent Cij Triclinic All 21 Monoclinic unique axis b C11, C22, C33, C12, C13 , C23, C44, C55, C66, C15, C25, C35, C46 unique axis c C11, C22, C33, C12, C13, C23 , C44, C55, C66, C16, C26, C36, C45 Orthorhombic C11, C22, C33, C12, C13, C23, C44, C55, C66 Cubic C11, C12, C44 Tetragonal Classes 4, ?4, 4/m C11, C33, C12, C13, C44, C66, C16 Classes 4mm, ?42m, 422, 4/mmm C11, C33, C12, C13, C44, C66 Trigonal Classes 3, ?3 C11, C33, C12, C13, C44, C14, C15 Classes 32, ?3m, 3m C11, C33, C12, C13, C44, C14 Hexagonal C11, C33, C12, C13, C44 Isotropic C11, C12 18 4 Equation B.8 is the foundation for evaluating the elastic constants. The procedure of our calculation can be summarized as three steps. 1) Based on the symmetry of the material, we generate a series of strained POSCAR ?les. 2) We calculate the total energy for each strained POSCAR. 3) Obtain the independent elastic constants. Diagonal and o?-diagonal terms of Cij are treated di?erently. Taking cubic system for illustration, there are three independent elastic constants: C11, C12 and C44. For C11 and C44, we introduce single strain ?1 and ?4, respectively, to the unstrained system. Then total energies of the strained systems are calculated. Because C11 and C44 are both diagonal terms, they are obtained in exactly the same way. Let?s take C11 as an instance. First, assuming optimized CONTCAR ?le at a given volume is acquired in previous static total energy calculations, it serves as the unstrained POSCAR in the elasticity calculation. When single strain is applied to the system, the total energy will be changed due to strain and can be expanded as a Taylor?s series of strain, provided the strain is small (usually between 0.5% and 1.5%). E = a + b? + 12c?2 + 16d?3 + 124e?4 + ??? (B.9) where a, b, c, d and e are coe?cients of the ?rst ?ve terms in Taylor?s series. We keep up to the ?fth term as a good approximation. C11 = c (B.10) C11 is the coe?cient of the quadratic term because the elastic energy due to strain ?1 is U = 12C11?21. We have implemented a ?six-point? method to get this coe?cient. Six single ?1-strains with the amount of 3?, 2?, ?, ??, ?2?, ?3? are applied to the original POSCAR and total energy calculations are carried out for each of them. ? ??? ??? ? ??? ??? ? E1 = a + 3b? + 92c?2 + 92d?3 + 278 e?4 E2 = a + 2b? + 2c?2 + 43d?3 + 23e?4 E3 = a + b? + 12c?2 + 16d?3 + 124e?4 E4 = a?b? + 12c?2 ? 16d?3 + 124e?4 E5 = a? 2b? + 2c?2 ? 43d?3 + 23e?4 E6 = a? 3b? + 92c?2 ? 92d?3 + 278 e?4 (B.11) By canceling the terms of no interests, we have C11 = c = ?3E1 + 16E2 ? 13E3 ? 13E4 + 16E5 ? 3E624?2 (B.12) For C12, we introduce two strains ?1 and ?2, at the same time, to the unstrained system. The related elastic energy is U = C12?1?2 + C11?21 + C22?22 (B.13) C12 is the coe?cient in the front of the ?1?2 term. With pair strains, the total energy of the system can also be expressed in a Taylor?s series, being cut at the fourth-order terms 185 as an approximation. E = a + b?1 + c?2 + d?1?2 + 12e?21 + 12f?22 + 16g?31 + 12h?21?2 + 12i?1?22 +16j?32 + 124k?41 + 16l?31?2 + 14m?21?22 + 16n?1?32 + 124o?42 (B.14) where C12 = d. To eliminate the unrelated coe?cients, an ?eight-point? algorithm is implemented. The eight pair strains for (?1,?2) we apply are: (2?,2?), (2?,?2?), (?2?,?2?), (?2?,2?), (?,?), (?,??), (??,??) and (??,?). ? ??? ??? ???? ???? ??? ???? ???? ??? ??? ??? ???? ???? ??? ???? ???? ??? E1 = a + 2b?1 + 2c?2 + 4d?1?2 + 2e?21 + 2f?22 + 43g?31 + 4h?21?2 + 4i?1?22 +43j?32 + 23k?41 + 83l?31?2 + 2m?21?22 + 83n?1?32 + 23o?42 E2 = a + 2b?1 ? 2c?2 ? 4d?1?2 + 2e?21 + 2f?22 + 43g?31 ? 4h?21?2 + 4i?1?22 ?43j?32 + 23k?41 ? 83l?31?2 + 2m?21?22 ? 83n?1?32 + 23o?42 E3 = a? 2b?1 ? 2c?2 + 4d?1?2 + 2e?21 + 2f?22 ? 43g?31 ? 4h?21?2 ? 4i?1?22 ?43j?32 + 23k?41 + 83l?31?2 + 2m?21?22 + 83n?1?32 + 23o?42 E4 = a? 2b?1 + 2c?2 ? 4d?1?2 + 2e?21 + 2f?22 ? 43g?31 + 4h?21?2 ? 4i?1?22 +43j?32 + 23k?41 ? 83l?31?2 + 2m?21?22 ? 83n?1?32 + 23o?42 E5 = a + b?1 + c?2 + d?1?2 + 12e?21 + 12f?22 + 16g?31 + 12h?21?2 + 12i?1?22 +16j?32 + 124k?41 + 16l?31?2 + 14m?21?22 + 16n?1?32 + 124o?42 E6 = a + b?1 ?c?2 ?d?1?2 + 12e?21 + 12f?22 + 16g?31 ? 12h?21?2 + 12i?1?22 ?16j?32 + 124k?41 ? 16l?31?2 + 14m?21?22 ? 16n?1?32 + 124o?42 E7 = a?b?1 ?c?2 + d?1?2 + 12e?21 + 12f?22 ? 16g?31 ? 12h?21?2 ? 12i?1?22 ?16j?32 + 124k?41 + 16l?31?2 + 14m?21?22 + 16n?1?32 + 124o?42 E8 = a?b?1 + c?2 ?d?1?2 + 12e?21 + 12f?22 ? 16g?31 + 12h?21?2 ? 12i?1?22 +16j?32 + 124k?41 ? 16l?31?2 + 14m?21?22 ? 16n?1?32 + 124o?42 (B.15) To solve for C12, we do addition or subtraction with these equations and ?nd that C12 = d = ?E1 + E2 ?E3 + E4 + 16E5 ? 16E6 + 16E7 ? 16E848?2 (B.16) So far we have discussed the derivation and computation procedure of the elastic con- stants for an unstressed solid. However, elastic constants will change under stress. Here we follow Barron & Klein?s derivation and care only about the case of isotropic stress (hydro- static pressure). When pressure is applied to the solid, Barron showed that the strain-energy density relation becomes: Cijkl = 1V ?E?? ij??kl + 12P (2?ij?kl ??il?jk ??ik?jl) (B.17) where V is the volume at that pressure, P is the pressure and ?ij is the Kronecker delta (?ij = 1 when i = j; ?ij = 0 when i negationslash= j). The ?rst term in equation B.17 is obtained from the strain-energy density relation as if the solid is not under stress. The second term is a correction because of pressure. In Voigt?s notation, the correction terms for each elastic constant is shown in table B.2. 186 Table B.2: Correction term to Cij due to pressure Cij Correction term Cij Correction term Cij Correction term C11 0 C13 P C26 0 C22 0 C14 0 C34 0 C33 0 C15 0 C35 0 C44 ?12P C16 0 C36 0 C55 ?12P C23 P C45 0 C66 ?12P C24 0 C46 0 C12 P C25 0 C56 0 Again, if we use cubic crystal as an example, and let ?C represent the ?rst term in equation B.17, the three independent elastic constants are: ? ? ? C11 = ?C11 C12 = ?C12 + P C44 = ?C44 ? 12P (B.18) At the current stage, we only calculate the elastic constant and its pressure dependence at zero temperature condition. 187 Appendix C Beyond Harmonic Approximation: Perturbation Theory of Lattice Anharmonicity The quasi-harmonic approximation (QHA) has achieved a great success in predicting thermodynamic properties for many hard materials at temperatures not too close to the melting point. However, to improve the prediction of thermodynamic properties of normal materials at high temperatures or strongly anharmonic crystals even at low temperatures, it is desired to add the anharmonicity contribution to the free energy as a next level ap- proximation. This correction is usually small compare to the harmonic contribution. With the third and fourth order anharmonicity, we can evaluate the lowest order of anharmonic contribution to the harmonic free energy via perturbation theory. Within this theory, we treat the anharmonic terms in the hamiltonian as perturbations. The total hamiltonian can be written as: H = HH + HA = H0 + H2 + H3 + H4 + ??? (C.1) where HH = H0 + H2 which is the harmonic part and HA = summationtext?n=3 Hn is the anharmonic part. Hn is the nth order term in the expansion of hamiltonian. As the anharmonic e?ect is more important at high temperatures, we deal this problem within the classical limit. Then the partition function of the system is: Z = h?3N integraldisplay dx1 ??? integraldisplay dx3N integraldisplay dp1 ??? integraldisplay dp3N ? e??(HH+HA) = ZH ? integraltext dx 1 ??? integraltext dx 3N integraltext dp 1 ??? integraltext dp 3N ? e??HH ? e??HAintegraltext dx1 ???integraltext dx3N integraltext dp1 ???integraltext dp3N ? e??HH = ZH ?ZA (C.2) where ZH = h?3N integraltext dx1???integraltext dx3N integraltext dp1 ???integraltext dp3N ?e??HH is the partition function of the harmonic hamiltonian and ZA is the ensemble average of the quantity e??HA over the 188 unperturbed harmonic system. Expanding e??HA into Taylor series, we have ZA = angbracketleftBig e??HA angbracketrightBig 0 = angbracketleftBigg 1 + ?summationdisplay n=1 (?1)n n! ? n ?Hn A angbracketrightBigg = 1 + ?summationdisplay n=1 (?1)n n! ? n ??Hn A?0 (C.3) Since we have expressed the total partition function as a product of the harmonic partition function and ZA, the anharmonic contribution to the Helmholtz free energy is: FA = ?kBT lnZA = ???1 ln parenleftBigg 1 + ?summationdisplay n=1 (?1)n n! ? n ??Hn A?0 parenrightBigg (C.4) Assuming the anharmonic terms are small compared to unity, the logarithm function in equation C.4 can be expanded into series following ln (1 + x) = x? x22 + x33 ???? when ?1 < x < 1. The expansion is expressed in terms of the (??)nn! and the coe?cient of (??)nn! is called the nth cumulant ?HnA?0,C FA = ???1 ?summationdisplay n=1 (?1)n n! ? n ??Hn A?0,C (C.5) Comparing equation C.5 with equation C.4, we have ?HA?0,C = ?HA?0 (C.6a) angbracketleftbigH2 A angbracketrightbig 0,C = angbracketleftbigH2 A angbracketrightbig 0 ? (?HA?0) 2 (C.6b) angbracketleftbigH3 A angbracketrightbig 0,C = angbracketleftbigH3 A angbracketrightbig 0 ? 3 angbracketleftbigH2 A angbracketrightbig 0?HA?0 + 2 (?HA?0) 3 (C.6c) According to pairing theorem243, the lowest order anharmonic contribution from ?HA?0,C term is ?H4?0. The lowest order anharmonic contribution from angbracketleftbigH2Aangbracketrightbig0,C term is angbracketleftbigH23angbracketrightbig0 ? (?H4?0)2. Other higher order terms are ignored as an approximation. The anharmonicity correction to the free energy under this approximation is: FA = ?H4?0 ? ?2 bracketleftBigangbracketleftbig H23angbracketrightbig0 ? (?H4?0)2 bracketrightBig (C.7) 189 Using normal coordinates, the third and fourth order anharmonic terms of hamiltonian can be written as: H3 = summationdisplay q,q?,q?? summationdisplay ?,??,??? W3 parenleftbigg q,q?,q?? ?,??,??? parenrightbigg Q? (q)Q??parenleftbigq?parenrightbigQ???parenleftbigq??parenrightbig (C.8a) H4 = summationdisplay q,q?,q??,q??? summationdisplay ?,??,???,???? W4 parenleftbigg q,q?,q??,q??? ?,??,???,???? parenrightbigg Q? (q)Q??parenleftbigq?parenrightbigQ???parenleftbigq??parenrightbigQ????parenleftbigq???parenrightbig (C.8b) where W3 and W4 are de?ned as W3 parenleftbigg q,q?,q?? ?,??,??? parenrightbigg = 16 summationdisplay ??? summationdisplay ijk N?12?q+q?+q??,G?A?i,?j,?kparenleftbigq?,q??parenrightbig ?e?,i (q,?)e?,jparenleftbigq?,??parenrightbige?,kparenleftbigq??,???parenrightbig (C.9a) W4 parenleftbigg q,q?,q??,q??? ?,??,???,???? parenrightbigg = 124 summationdisplay ???? summationdisplay ijkl N?1?q+q?+q??+q???,G?B?i,?j,?k,?lparenleftbigq?,q??,q???parenrightbig ?e?,i (q,?)e?,jparenleftbigq?,??parenrightbige?,kparenleftbigq??,???parenrightbige?,lparenleftbigq???,????parenrightbig (C.9b) In equation C.9, A?i,?j,?k (q?,q??) is the third order anharmonic dynamical tensor, which can be derived from the real-space third order force-constant tensor through a Fourier transformation. B?i,?j,?k,?l (q?,q??,q???) is the fourth order anharmonic dynamical tensor. A?i,?j,?kparenleftbigq?,q??parenrightbig = 1?m imjmk summationdisplay h,h? A?i,?j,?kparenleftbig0,h,h?parenrightbig? ei(q??h+q???h?) (C.10a) B?i,?j,?k,?lparenleftbigq?,q??,q???parenrightbig = 1?m imjmkml summationdisplay h,h?,h?? B?i,?j,?k,?lparenleftbig0,h,h?,h??parenrightbig? ei(q??h+q???h?+q????h??) (C.10b) With the pairing theorem we can evaluate the averaged hamiltonian over the unper- turbed canonical ensemble. ?H4?0 = 3?3 summationdisplay q1,q2 summationdisplay ?1,?2 W4 parenleftbigg q 1,?q1,q2,?q2 ?1,?1,?2,?2 parenrightbigg ?2?1 (q1)?2?2 (q2) (C.11a) angbracketleftbigH2 3 angbracketrightbig 0 = 6 ?3 summationdisplay q1,q2,q3 summationdisplay ?1,?2,?3 W3 parenleftbigg q 1,q2,q3 ?1,?2,?3 parenrightbigg W3 parenleftbigg ?q 1,?q2,?q3 ?1,?2,?3 parenrightbigg ?2?1 (q1)?2?2 (q2)?2?3 (q3) + 9?3 summationdisplay q1,q2,q3 summationdisplay ?1,?2,?3 W3 parenleftbigg q 1,?q1,q3 ?1,?1,?3 parenrightbigg W3 parenleftbigg q 2,?q2,?q3 ?2,?2,?3 parenrightbigg ?2?1 (q1)?2?2 (q2)?2?3 (q3) (C.11b) 190 where ?? (q) is the phonon frequency of the ?th normal mode at reciprocal lattice point q. We have proposed an algorithm to evaluate the third order force-constant tensor A?i,?j,?k (?,??,???). In order to get all the tensor elements, we apply a pair of displace- ments to the system. Here we keep up to the 4th order term of the hamiltonian and the amount of each displacement is ? which is much smaller than the interatomic distance. The pair-displacement scheme we adopt is to deviate the jth atom in ? direction and the kth atom in ? direction by (?,?), (?,??), (??,?), and (??,??), respectively. The ? component of the H-F forces on the ith atom are: F++?,i (?) = ?bracketleftbig??i,?jparenleftbig?,??parenrightbig+ ??i,?kparenleftbig?,???parenrightbigbracketrightbig? ? ?12 bracketleftbigA?i,?j,?jparenleftbig?,??,??parenrightbig+ A?i,?k,?kparenleftbig?,???,???parenrightbig+ 2A?i,?j,?kparenleftbig?,??,???parenrightbigbracketrightbig? ?2 ?16 bracketleftbigg B ?i,?j,?j,?j (?,??,??,??) + B?i,?k,?k,?k (?,???,???,???) +2B?i,?j,?j,?k (?,??,???,????) + 2B?i,?j,?k,?k (?,??,???,????) bracketrightbigg ? ?3 (C.12a) F+??,i (?) = ?bracketleftbig??i,?jparenleftbig?,??parenrightbig? ??i,?kparenleftbig?,???parenrightbigbracketrightbig? ? ?12 bracketleftbigA?i,?j,?jparenleftbig?,??,??parenrightbig+ A?i,?k,?kparenleftbig?,???,???parenrightbig? 2A?i,?j,?kparenleftbig?,??,???parenrightbigbracketrightbig? ?2 ?16 bracketleftbigg B ?i,?j,?j,?j (?,??,??,??) ?B?i,?k,?k,?k (?,???,???,???) ?2B?i,?j,?j,?k (?,??,???,????) + 2B?i,?j,?k,?k (?,??,???,????) bracketrightbigg ? ?3 (C.12b) F?+?,i (?) = ?bracketleftbig???i,?jparenleftbig?,??parenrightbig+ ??i,?kparenleftbig?,???parenrightbigbracketrightbig? ? ?12 bracketleftbigA?i,?j,?jparenleftbig?,??,??parenrightbig+ A?i,?k,?kparenleftbig?,???,???parenrightbig? 2A?i,?j,?kparenleftbig?,??,???parenrightbigbracketrightbig? ?2 ?16 bracketleftbigg ?B ?i,?j,?j,?j (?,??,??,??) + B?i,?k,?k,?k (?,???,???,???) +2B?i,?j,?j,?k (?,??,???,????) ? 2B?i,?j,?k,?k (?,??,???,????) bracketrightbigg ? ?3 (C.12c) F???,i (?) = +bracketleftbig??i,?jparenleftbig?,??parenrightbig+ ??i,?kparenleftbig?,???parenrightbigbracketrightbig? ? ?12 bracketleftbigA?i,?j,?jparenleftbig?,??,??parenrightbig+ A?i,?k,?kparenleftbig?,???,???parenrightbig+ 2A?i,?j,?kparenleftbig?,??,???parenrightbigbracketrightbig? ?2 +16 bracketleftbigg B ?i,?j,?j,?j (?,??,??,??) + B?i,?k,?k,?k (?,???,???,???) +2B?i,?j,?j,?k (?,??,???,????) + 2B?i,?j,?k,?k (?,??,???,????) bracketrightbigg ? ?3 (C.12d) From equation C.12, simple derivation yields A?i,?j,?kparenleftbig?,??,???parenrightbig = ?F ++ ?,i (?) + F +? ?,i (?) ?F ?+ ?,i (?) + F ?? ?,i (?) 4?2 (C.13) For a 3N?3N?3N tensor, there are C26N +6N = 18N2+3N H-F forces required. This number can be reduced by noticing that when (?,j,??) = (?,k,???), A?i,?j,?j (?,????) can be obtained from the single-displacement calculations which have been done when acquiring 191 the second order force-constant matrix. A?i,?j,?jparenleftbig?,??,??parenrightbig = F + ?,i (?) + F ? ?,i (?) ?F 2+ ?,i (?) ?F 2? ?,i (?) 3?2 (C.14) However, there are still C26N forces needed which is not an easy task to calculate for a supercell. Again, the crystal symmetry plays an important role to greatly reduce the number of direct calculations. A code named as Moves Analysis.x has been implemented to ?gure out the irreducible number of pair-moves. Let us take the 120-atom supercell of ?-Al2O3 as an example. The total number of pair-moves is 258840 initially. After the symmetry analysis, the number of independent pair-moves is reduced to 5168. Although this number has been greatly deducted, it is still impractical to carry out ?rst-principles calculations for such a large number. At this stage, an approximation could be made that the tensor element is typically around zero if two atoms in the atomic indices of A?i,?j,?k (?,??,???) are far away from each other. If we drop all the calculations involving the distance between the jth and kth atoms that beyond the second nearest neighbor, we can reduce the number of direct pair-move calculations down to 632, which is now a feasible job. After the direct calculations, other tensor elements can be either reconstructed from the irreducible ones or equal to zero. 192