THEORETICAL STUDY OF THERMAL PROPERTIES AND THERMAL CONDUCTIVITIES OF CRYSTALS 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 classified information. Xiaoli Tang Certificate of Approval: Minseo Park Jianjun Dong, Chair Associate Professor Associate Professor Physics Physics Yu Lin John Williams Professor Professor Physics Physics George T. Flowers Interim Dean Graduate School THEORETICAL STUDY OF THERMAL PROPERTIES AND THERMAL CONDUCTIVITIES OF CRYSTALS Xiaoli Tang A Dissertation Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy Auburn, Alabama August 9, 2008 iii THEORETICAL STUDY OF THERMAL PROPERTIES AND THERMAL CONDUCTIVITIES OF CRYSTALS Xiaoli Tang Permission is granted to Auburn University to make copies of this dissertation at its direction, upon the request of individuals or institutions and at their expense. The author reserves all the publication rights Signature of Author Date of Graduation iv VITA Xiaoli Tang, daughter of Jiakuan Tang and Lifang Gu, was born on January 9 th , 1980, in the city of Dongtai, Jiangsu province, People?s Republic of China. She entered University of Shanghai for Science and Technology in September, 1998, and graduated with the Bachelor of Engineering degree in Fluid Mechanics and Engineering in July, 2002. She entered the graduate program in Physics Department of Auburn University in August, 2003 and joined Condensed Matter Theory Group in January 2005. Since then she is dedicated to her Ph.D research program. She married Yi Zhou in Opelika, Alabama on October 1 st , 2003, and her son, Kevin Zhou, was born on March 11 th , 2007 in East Alabama Medical Clinic, Opelika, Alabama. v DISSERTATION ABSTRACT THEORETICAL STUDY OF THERMAL PROPERTIES AND THERMAL CONDUCTIVITIES OF CRYSTALS Xiaoli Tang Doctor of Philosophy, August 9, 2008 (B.E. University of Shanghai for Science and Technology, 2002) 172 Typed Pages Directed by Jianjun Dong Standard first-principles total energy and force calculations were extended to the calculations of harmonic force constant matrices and third order lattice anharmonicity tensors with an efficient super-cell finite-difference algorithm. Phonon spectra calculated within this algorithm are in excellent agreement with other theoretical results calculated with density perturbation functional theory, as well as the available experimental measurements. The newly proposed algorithm for lattice anharmonicity was implemented with both empirical Tersoff potentials and first-principle density functional theory methods. A self testing scheme for the validity of 3 rd order lattice anharmonicity vi was also proposed and sumrule enforcement has been investigated to ensure the numerical accuracy. Statistical ensemble of phonons was then adopted to calculate and simulate the equilibrium thermal properties of solid materials. With the forces calculated from first- principle theory, fundamental thermal properties such as heat capacity, thermal expansion were calculated within the quasi-harmonic approximation. Kinetic theory was implemented to predict the non-equilibrium thermal transport properties such as phonon life time and thermal conductivity. With the newly developed computational method, we have studied the thermal and thermal transport properties of two material systems Si 136 and MgO. Our calculation predicted that a negative thermal expansion exist in Si 136 at temperature lower than 124K, and was then confirmed by experimental measurement. Green-Kubo calculation yielded 90% reduction of thermal conductivity in Si 136 compared with diamond structured Si. Cause of this reduction was then investigated using kinetic transport theory. For MgO, the pressure dependence of lattice anharmonicity was studied. Both intrinsic anharmonicity and extrinsic isotope induced phonon scattering have been considered. The isotope effect on the lattice thermal conductivity was discussed. Preliminary results of lattice thermal conductivity at a wide range of temperature were then presented. At room temperature, our theory calculated lattice thermal conductivity is 51 W/K/m, in a good agreement with experimental measurement 54 W/K/m. In addition, two models were proposed to estimate the pressure dependence of the lattice thermal conductivity. vii ACKNOWLEDGEMENTS This dissertation could not have been finished without the help and support of many people who are gratefully acknowledged here. First of all, I would like to express my deepest gratitude to my dedicated research supervisor, Prof. Jianjun Dong, who has not only offered me great suggestions and ideas regarding my research work, but also encouraged and challenged me through my Ph.D program. I would also like to thank him for his patience and kindness, especially his trust on me working remotely during my last year of Ph.D study. His great effort on helping me through the dissertation writing is greatly appreciated. I would like to thank Dr. An-Ban Chen, who initiated my research interest in theoretical/computational condensed matter physics and gave me a lot great advice on my research study. And I would also like to thank Prof. John R. Williams and Prof. Yu Lin for their valuable suggestions on the dissertation writing. Great thanks to my graduate friends, Bin Xu, for valuable research discussions and encouragement and Zengjun Chen for helping me with the paper work toward graduation while I am not in Auburn. I would also like to express my sincere thanks to my parents, my sister, my parents-in-law, my husband and my little baby Kevin for their continuous love and support. Last but not least, I would like to thank NSF (EPS-0447675 and HRD-0317741) and DOE (DE-FG02-03ER46060) for their financial support to my Ph.D research study. viii Style manual or journal used: Bibliography conforms to those of the Physics Review B Computer software used: Microsoft Word 2003 for Windows ix TABLE OF CONTENTS LIST OF FIGURES .......................................................................................................... xii LIST OF TABLES........................................................................................................... xvi CHAPTER 1 INTRODUCTION........................................................................................ 1 CHAPTER 2 FIRST-PRINCIPLES ATOM-SCALE SIMULATION AND MODELING OF CRYSTALLINE SOLIDS ............................................................................................ 8 2.1 First-Principles Density Functional Theory (DFT)............................................... 8 2.2 Atom-Scale Simulation and Modeling................................................................ 12 2.2.1 Structural Models and Crystal Symmetries ............................................. 12 2.2.2 Energy Minimization and Equation of State:........................................... 19 2.2.3 Lattice Vibration .......................................................................................21 CHAPTER 3 STATISTICAL THEORIES OF PHONONS............................................. 37 3.1 Quantum Theory of Lattice Vibration ................................................................ 37 3.1.1 Phonons: Harmonic Frequencies and Group Velocities ......................... 37 3.1.2 Quasi-Harmonic Approximation and Mode Gr?neisen Parameters ........ 39 3.1.3 Anharmonicity Induced Phonon-Phonon Scatterings.............................. 40 3.2 Equilibrium Thermal Properties ......................................................................... 43 3.2.1 Statistical Ensemseble Theory of Independent Phonons ......................... 43 3.2.2 Macroscopic Thermal Properties within the QHA .................................. 44 x 3.3 Non-Equilibirum Thermal Transport Theories................................................... 46 3.3.1 Kinetic Transport Theory......................................................................... 46 3.3.2 Green-Kubo Formula............................................................................... 51 CHAPTER 4 THERMAL PROPERTIES OF TYPE II CLATHRATE Si 136 .................. 54 4.1 Introduction......................................................................................................... 54 4.2 Crystal Lattices and Static Equation of State...................................................... 57 4.3 Lattice Phonon Spectra ....................................................................................... 59 4.4 Thermodynamic Potentials and T-P Phase Relations ......................................... 69 4.5 Thermal Properties.............................................................................................. 71 4.6. Conclusions........................................................................................................ 74 CHAPTER 5 THERMAL CONDUCTIVITY OF TETRHEDRALLY BONDED SILICON CRYSTALS ? THE GREEN-KUBO APPROACH ........................................ 76 5.1 Introduction......................................................................................................... 76 5.2 Empirical Si-Si Potential..................................................................................... 79 5.3 Results and Discussions...................................................................................... 85 5.4 Conclusion .......................................................................................................... 95 CHAPTER 6 LATTICE ANHARMONICITY OF TETRAHEDRALLY BONDED SILICON CRYSTALS ..................................................................................................... 97 6.1 Introduction......................................................................................................... 97 6.2 Single Effective Phonon Life Time Approximation........................................... 98 6.3 Lattice Anharmonicity ...................................................................................... 100 6.4 Conclusion ........................................................................................................ 112 CHAPTER 7 LATTICE THEMRAL CONDUCTIVITY OF MgO .............................. 113 xi 7.1 Introduction....................................................................................................... 113 7.2 Crystal Structure and Equation of State............................................................ 116 7.3 Phonon and Gr?neisen Parameters ................................................................... 117 7.4 Lattice Anharmonicity ...................................................................................... 121 7.5 Isotope Effect.................................................................................................... 123 7.6 Thermal Conductivity ....................................................................................... 126 7.7 Conclusion ........................................................................................................ 129 CHAPTER 8 SUMMARY AND FUTURE WORK ...................................................... 131 REFERENCES .............................................................................................................. 140 APPENDIX A NOTATION FOR SPACE GROUP AND POINT GROUP ................. 143 APPENDIX B NORMAL MODE SYMMETRY REPRESENTATION ANALYSIS.. 149 xii LIST OF FIGURES Figure 2.1 Comparison of Gr?neisen parameters calculated from FDA (blue circle) and AA, which include both sumrule un-enforced (red diamond) and sumrule enforced (black square) results.. ................................................................................................................. 41 Figure 4.1 LDA calculation of (a) phonon dispersion relations, (b) vibrational density of states, and (c) mode Gr?neisen parameters of d-Si at the equilibrium volume. ............... 59 Figure 4.2 LDA calculation of (a) phonon dispersion relations, (b) vibrational density of states, and (c) mode Gr?neisen parameters of Si 136 clathrate at the equilibrium volume. 59 Figure 4.3 (a) 4-fold ?folded? phonon dispersion plots of d-Si, calculated with unit-cell size four times that of the primitive cell ( 4 cell primitive aa= ?=2.16 nm), along with (b) the re-plotted phonon dispersion of Si 136 ( clathrate a = 1.46 nm)................................................ 60 Figure 4.4 The theoretical mode Gr?neisen parameters ? i along three high-symmetry directions. The ? i of LA (cross) and optic (triangles) phonon modes are all positive; the LA branches have values that lie within a narrow range (~+0.75) and those associated with optic modes range between +0.8 to +1.6. ................................................................. 63 Figure 4.5 Raman spectrum obtained for Si 136 at ambient P and T. There is a weak peak at ~520 cm-1 that corresponds mainly to a trace of d-Si impurity in the sample. However, there is also a calculated mode at 516 cm-1 for Si 136 clathrate at this position.. .............. 45 Figure 4.6 Selected Raman spectra of Si 136 collected at compression. Asterisk marks diamond-phase silicon.. .................................................................................................... 46 Figure 4.7 Pressure dependence of Raman shifts of Si 136 and d-Si. ................................ 47 Figure 4.8 The LDA predicted vibrational entropies in d-Si (solid line) and Si 136 (dashed line). .................................................................................................................................. 69 Figure 4.9 The theoretically predicted equilibrium T-P phase boundary between the ground state diamond phase and the Si 136 clathrate phase. Note that the transition to the ?expanded? polymorph of the element occurs within the tensile regime, at negative pressure. ............................................................................................................................ 70 xiii Figure 4.10 The calculated and measured specific isobaric heat capacities in d-Si and Si 136 . .................................................................................................................................. 71 Figure 4.11 The theoretically predicted linear coefficients of lattice thermal expansion at P=0 in d-Si and Si 136 . ........................................................................................................ 73 Figure 4.12 High-resolution X-ray powder diffraction data for Si 136 illustrated at three temperatures (5K ? bottom, 145K - middle, 275K - top). The three high angle reflections used in the data fitting are marked with an asterisk and represent (in order) the 733, 660 and 751 reflections............................................................................................................ 73 Figure 4.13 Low temperature variation (T = 5 - 275 K) variation in unit cell lattice constant (?) for clathrate-structured Si 136 obtained by fitting to high angle powder X-ray diffraction data. The line drawn through the data points is a guide to the eye. ................ 74 Figure 5.1 Phonon dispersion of d-Si using Tersoff empirical potential.......................... 84 Figure 5.2 Phonon dispersion of Si 136 using Tersoff empirical potential. ........................ 84 Figure 5.3 Macroscopic heat current fluctuations with time of d-Si system .................... 86 Figure 5.4 The normalized time correlation function ()g t at room temperature for d-Si. (a) an overall look; (b) close-up look at the beginning of the MD run; (c) close-up look at decaying of ()g t after a long time MD run....................................................................... 87 Figure 5.5 Macroscopic heat current fluctuations with time of Si 136 system.................... 89 Figure 5.6 The normalized time-dependent correlation function ()g t at room temperature for Si 136 . (a) an overall look; (b) close-up look at the beginning of the MD run; (c) close- up look at decaying of ()g t after a long time MD run...................................................... 89 Figure 5.7 Power Spectrum of heat current autocorrelation function for d-Si. ................ 91 Figure 5.8 Power Spectrum of heat current autocorrelation function for Si 136 ................. 92 Figure 5.9 Normalized the heat current autocorrelation of d-Si before (black line) and after (red line) the low-pass filter with filtering frequency 0.5Hz.................................... 94 Figure 5.10 Normalized the heat current autocorrelation of Si 136 before (black line) and after (red line) the low-pass filter with filtering frequency 0.5Hz.................................... 94 Figure 5.11 log plot of the normalized low-pass (0.5Hz) filtered heat current autocorrelation of d-Si (black line) and Si 136 (red line). ................................................... 95 xiv Figure 6.1 Fitted effective phonon life time eff ? of d-Si to its experimental heat conductivity. Temperature dependence of Si 136 was estimated assuming the equal eff ? ... 99 Figure 6.2 Atomic model of d-Si .................................................................................... 101 Figure 6.3 Atomic model of Si 136 . Three different types of atoms have been indicated by different color: 8a (red), 32e (blue) and 96g (gold)........................................................ 102 Figure 6.4 A schematic show of coordinates transformation to yield a bond-stretching A ijk .................................................................................................................................. 104 Figure 6.5 Any bond between two atoms that are tetrahedrally bonded can be rotated to z? direction like what we showed here................................................................................ 105 Figure 6.6 Pressure dependence of bond stretching Aijk terms in d-Si.......................... 106 Figure 6.7 Pressure dependence of 8a-32e bond stretching A ijk terms in Si 136 . ............. 106 Figure 6.8 Pressure dependence of 32e-96g bond stretching Aijk terms in Si 136 ........... 107 Figure 6.9 Pressure dependence of 96g-96g type (a) bond stretching A ijk terms in Si 136 . ......................................................................................................................................... 107 Figure 6.10 Pressure dependence of 96g-96g type (b) bond stretching A ijk terms in Si 136 . ......................................................................................................................................... 108 Figure 6.11 Mode Gr?neisen parameters of d-Si calculated from FDA and AA. .......... 111 Figure 6.12 Mode Gr?neisen parameters of Si 136 calculated from FDA and AA........... 111 Figure 7.1 LDA calculated harmonic phonon spectra, phonon density of state and group velocity of MgO at 0 GPa. Experiment data of phonon dispersion measured by Peckham (1967) and Sangster et al. (1970) are also shown for comparison.................................. 118 Figure 7.2 LDA calculated phonon spectra, phonon density of state and group velocity of MgO at 68 GPa. .............................................................................................................. 119 Figure 7.3 LDA calculated Gr?neisen parameter as a function of temperature at 0 GPa. ......................................................................................................................................... 119 Figure 7.4 The pressure dependence of the average and difference (inset) of the two largest anharmonic tensor elements................................................................................ 122 xv Figure 7.5 Comparison of Phonon mode Gr?neisen parameters calculated with both finite difference method and the 3 rd order lattice anharmonicity method. ............................... 123 Figure 7.6 Temperature dependence of lattice thermal conductivity ? at 0 GPa; ? was calculated using First-Principles method with a 444? ? q-point sampling over the Brillouin zone. Experimental measurements are also cited for comparison................... 127 Figure 7.7 Estimated pressure dependence of lattice thermal conductivity in MgO at 500K for model I and model II. ...................................................................................... 127 Figure 7.8 Estimated pressure dependence of lattice thermal conductivity in MgO at 2000K for model I and model II. Recent experimental results are shown in symbol plus (+).................................................................................................................................... 128 xvi LIST OF TABLES Table 2.1 One-to-one correspondence of the index for the displacement to the index of atom being displaced......................................................................................................... 26 Table 2.2 One-to-one correspondence of index for the pair displacement to that for the single atom displacement.................................................................................................. 32 Table 4.1 LDA calculated static (T=0K) Birch-Murnaghan equation of state of the ground state diamond phase Si (d-Si) and the meta-stable type-II clathrate phase of Si (Si 136 ), and available experimental parameters. The experimental data are for T = 298 K and they are taken from Ref [64] for d-Si and Ref. [65] for Si 136 . ........................................................ 58 Table 4.2 Comparison of acoustic velocities (m/s) for d-Si (experimental vs. LDA calculated) and Si 136 . The experimental data labeled with (a) and (b) are taken from Ref. and respectively. ..............................................................................................................61 Table 4.3 LDA calculated and experimentally measured Raman frequencies (?) and mode Gr?neisen parameters (? i ) for clathrate-structured Si 136 .......................................... 68 Table 5.1 Fitting parameters for silicon to be used in Equ. (5.1). .................................... 80 Table 5.2 Comparison between VASP results and Tersoff MD results: Equilibrium Birch- Murnaghan equation of state fitting parameters for d-Si and Si 136 together with experimental data.............................................................................................................. 80 Table 5.3 Si 136 ? point phonon frequencies from both Tersoff MD and VASP calculation. ........................................................................................................................................... 82 Table 6.1 List of coordinates of all 4 first nearest neighbors of two representative atoms at Wyckoff site 8a in our chosen super-cell of d-Si........................................................ 101 Table 6.2 List of coordinates of all 4 first nearest neighbors of one representative atom for each Wyckoff site in our chosen super-cell of Si 136 .................................................. 102 Table 6.3 Averaged site anharmonicity and overall bond anharmonicity of Si 136 .......... 109 xvii Table 6.4 Comparison of relative anharmoncity between d-Si and Si 136 ....................... 110 Table 7.1 Our LDA calculated static equilibrium properties of MgO fitted by 3 rd order Birch-Murnaghan equation of state, compared with previous theoretical results and experimental measurement. Experimental data was taken at room temperature............ 117 Table 7.2 Mass deviation of atom oxygen (g O ) and magnesium (g Mg )........................... 125 Table 7.3 Effects of isotope on the lattice thermal conductivity at 500K and 2000K.... 125 Table 7.4 List of ? in ()/( 0) 1PP P?? ?= =+ estimated from two models at 500K and 2000K respectively. ........................................................................................................ 129 1 CHAPTER 1 INTRODUCTION Atom-scale numerical simulation techniques, especially those based on the first- principles electronic-structure theories, have become powerful tools to understand structures, formations, dynamics, and many other physical and chemical properties of real and complex materials systems. First-principles simulations and calculations are capable of predicting various materials properties when no or limited experimental data are available, because they directly solve the electronic structures at the level of quantum mechanics and do not depend on any empirical parameters that are either suggested based on a priori theoretical assumptions or fitted with a set of experimental data. Therefore, first-principles techniques are ideal to (1) construct realistic atomistic structural models for complex materials systems, provide unbiased interpretation of complicated experimental data, and gain insights on chemical trends of materials properties, (2) predict materials properties at conditions where experimental measurements are not yet feasible, for example extreme high pressure, and (3) design novel artificial materials with optimized physical/chemical properties. 2 Despite rapid progresses on numerical algorithms and significant increases in computer speed and memory during the last two decades, the intensive computational loads limit the first-principles simulation studies to models less than a couple of hundred atoms and over periods no more than a few of thousand simulation steps on a typical single CPU workstation. Nevertheless, the recent development in parallel computers, especially the relatively low-cost Beowulf type computer clusters, provides new opportunities. A major part of the research work in this dissertation is to implement and further develop efficient parallel algorithms to calculate harmonic phonon spectra and lattice anharmonicity using first-principles density functional theory (DFT). We have successfully developed a real-space super-cell based algorithm that can accurately calculate both harmonic force constant matrices and 3 rd order lattice anharmonicity tensors for a wide range of material systems. To the best of our knowledge, our new numerical algorithm for the 3 rd order lattice anharmonicity calculation is the first of this kind. Group theory has been adopted to ensure the rotational symmetry relation among the tensor elements, and sumrule enforcement to the tensor are made to ensure the translational symmetry of the crystal system. In addition, evaluation of Gr?neisen parameters has been proposed to justify the accuracy of this complex tensor. Our first- principles-method calculated total energies, phonon spectra, and lattice anharmonicity are then adopted to predict (1) thermodynamic potentials of crystals at different (T,P) conditions within the statistical quasi-harmonic approximation (QHA), and (2) lattice thermal conductivity based on the kinetic transport theory. An alternative simulation 3 approach based on the Green-Kubo theory is also implemented to evaluate the lattice thermal conductivity as comparison. Results of first-principles calculation of two material systems are presented in this dissertation. The first is nano-open crystalline type-II silicon clathrate (Si 136 ) materials. This new Si allotrope has an open-cage structure, isostructural with low- density inclusion compounds of H 2 O-ice. It has a cubic framework in which each cubic unit-cell contains sixteen 20-atom ?cages? (dodecahedra) and eight 28-atom ?cages? (hexakaidecahedra). In addition to the elemental ?guest-free? form of Si 136 , various ?guest? atoms, including alkali or alkaline earth metals, or halogens, can be incorporated inside the atomic cages to form binary or ternary compounds. Both pristine and guest- encapsulated clathrate materials have significant technological potential because they exhibit a very wide spectrum of materials properties. In 1995, Slack predicted that open framework structures containing encapsulated rattling guest atoms may exhibit lowered ?glass-like? thermal conductivity due to scattering of acoustic heat-carrying phonons by the guest atoms, while leaving the electrical conductivity via the framework channels largely unaffected. 1 And then G.S. Nolas confirmed this prediction by measuring the conductivity for Sr 8 Ga 8 Ge 30 , which has thermal conductivity 2 orders of magnitude smaller than diamond structured Ge (d-Ge). 2 In 2000, J. Dong et al. 3 performed a theoretical calculation for both guest-free and guest-encapsulated type I Ge 46 , revealing that one order of reduction comes from the rattlers inside the cages, and the other order of reduction comes from open framework itself. After the guest-free type II clathrate was successfully synthesized by Cryko 4 , Nolas 5 measured the thermal conductivity for guest- free Si 136 in 2003, and the measurement showed that guest-free Si 136 itself exhibits a 4 rather low thermal conductivity, which suggests the open framework is the key reason for the reduction of the thermal conductivity. Meanwhile, little theoretical work has been devoted to the study of these unusual expanded-framework semiconducting crystals, including studies of both fundamental thermal properties and thermal conductivity. As a first step in this area, we use first-principles theoretical methods to predict the measurable thermal properties (such as heat capacity and thermal expansion) of the guest-free clathrate Si 136 . Our results are analyzed and then compared with previous data on the well-known ground state diamond-structured phase of this element (d-Si) 6 . Despite noticeable differences in materials density, compressibility, and electronic structures, we find that the two phases have very comparable heat capacities and thermal expansibilities. One important prediction of our calculations is that the clathrate-structured polymorph Si 136 should exhibit a region of negative thermal expansion below 124K, like the diamond-structured phase. This prediction has been confirmed by an experimental measurement. Then we calculate the time-correlation functions of heat current in both d- Si and Si 136 using classical molecular dynamics (MD) with an empirical Tersoff potential at equilibrium micro-canonical (N,E,V) conditions, and derive the non-equilibrium thermal transport properties based on statistical fluctuation-dissipation theory 7 (the Green-Kubo formula). The calculation indicates that thermal conductivity of Si 136 is only about 10.6% of that of d-Si. We further adopt the kinetic transport theory to quantitatively analyze the contributions from two major effects to this observed large reduction. Our calculation shows that the major contribution to the reduction of thermal conductivity in Si 136 comes from the flattering in phonon dispersion, which reduces the thermal conductivity by a factor of 0.8. The increase of lattice anharmonicity in Si 136 5 contributes to a further reduction by a factor of 0.3. The overall estimated reduction of lattice thermal conductivity reduction based on the kinetic theory is 14% ((10.8)(10.3)=? ?? , which is in consistent with the results obtained with the Green- Kubo Formula. The second material system studied is magnesium oxide (MgO), which is considered to be an end-member component of the lower mantle minerals. Thermal conductivity (? ) data of Earth?s constituent minerals are important for understanding any geophysical process that involves heat. 8,9,10 Probing the lattice anharmonicity in MgO is a precursor to studying more complex mineral structures and compositions relevant to the Earth. Although several rapid developments in experimental techniques were reported in recent years 11,12,13,14,15 , some pressure (P) and temperature (T) conditions of the Earth?s interiors (for example, T > 2300K or P > 100GPa) remain inaccessible for accurate measurement of ? at the current stage. Furthermore, the issue of contact associated errors for the thermal transport measurements has been raised and discussed. 14 The systematic errors of this type are especially important for accurately determining the pressure dependence in thermal transport properties. At the same time, little theoretical effort has been devoted to the first-principles calculation of this important thermal transport property of minerals, including ideal crystalline minerals (i.e. containing no isotope/composition disorder, no isolated or extended defects, or no finite-size grain boundaries). Current understanding on lattice anharmonicity and its pressure dependence is limited. Recently, Oganov and Dorogokupets reported a study on the anharmonicity effects on the thermodynamic potentials of MgO using a first-principles method 16 . In addition to the conventional quasi-harmonic approximation (QHA) results (including 6 both harmonic and anharmonic contributions), an additional correction term, whose magnitude scales as a function of T 2 , was estimated using MD simulations. Individual interatomic anharmonicity terms were not explicitly evaluated. The authors reported that at ambient pressure, the lattice anharmonicity evaluated with the QHA approach led to a noticeable overestimation in the lattice thermal expansion, an equilibrium thermal property that is believed to be closely related to lattice anharmonicity. In this work, we provide a first-principles calculation of harmonic phonon spectra, 3 rd order lattice anharmonicity in MgO. Explicit calculation of phonon relaxation time in both intrinsic anharmonicity and extrinsic isotope induced phonon scattering processes have been calculated. Lattice thermal conductivity at a wide temperature range and ambient pressure has been calculated with 444? ? q-point grids for the Brillouin zone integration. Within the single relaxation time approximation (SRTA), two models are proposed to estimate the pressure dependence of lattice thermal conductivity. The rest of this dissertation is organized as following: Chapters 2 and 3 review the fundamental theories adopted in our theoretical studies of materials and present some of the implementation details of computational methodologies. Chapter 4 presents a detailed first-principles prediction of equilibrium thermal properties of Si 136 clathrate. Chapter 5 reports an empirical (Tersoff) potential based molecular dynamics (MD) simulation study and an evaluation of lattice thermal conductivity based on the Green- Kubo theory of two tetrahedrally bonded Si crystals: d-Si and Si 136 . Chapter 6 studies the 3 rd order lattice anharmonicity and 1 st order approximation of phonon-phonon scattering rates using the same empirical (Tersoff) potentials and re-exams the ratio of lattice thermal conductivity between Si 136 and d-Si based on the kinetic transport theory. 7 Chapter 7 presents the first-principles prediction of pressure dependence of lattice anharmonicity in MgO and the temperature and pressure dependence of the lattice thermal conductivity. Isotope effect on the lattice thermal conductivity of MgO is also discussed. Finally, Chapter 8 concludes the key results of our studies and suggests related future research topics. 8 CHAPTER 2 FIRST-PRINCIPLES ATOM-SCALE SIMULATION AND MODELING OF CRYSTALLINE SOLIDS 2.1 First-Principles Density Functional Theory (DFT) An accurate total energy theory, which predicts the energy of an N-atom material system at a given structural configuration, is the foundation for all atom-scale simulation and modeling. Because of the large mass ratios between nuclei and electrons inside an atom and the fact that electronic structures of core-electrons are not environment-sensitive, an atom inside a solid are often considered as a positively-charged ion (including nuclei and core electrons) surrounded by a group of negatively charged valence electrons, and ions and electrons move at different time-scales. Born and Oppenheimer 17 (BO) proposed an adiabatic approximation, which assumes electrons respond instantaneously to the motion of ions. Within the BO approximation, the electrons are moving in an external potential from a configuration of static ions, and the total energy can be written as: ({ }) ({ }) ({ , }) ion ion I e e i ion e I i E ERErERr ??? =++ v v v v (2.1) where I R v is the position of th I ion with I from 1 to ion N , and i r r is the position of th i 9 valence electron with i from 1 to N e . While ion ion E ? is simply Coulomb energy. The exact solution of the last two terms, often referred as electronic energy, for a real bulk solid is a formidable task because it requires solving an interacting many-body Hamiltonian equation. Meanwhile, a wide range of total energy theories, from simple empirical two- body force field types to computationally intensive quantum Monte Carlo methods, have been developed to approximate the electronic energies of various solid-state materials systems. First-principles electronic-structure based total energy theories refer to the class of theories that (1) explicitly describe the motion of electrons with quantum mechanics and (2) adopt no empirical fitting parameters. Most of the calculations reported in this dissertation are calculated within the Density Functional Theory (DFT), which have been successfully adopted to predict structural, elastic and dynamical properties for a wide range of solids in the past 20 to 30 years. The DFT was proposed by Hohenberg and Kohn (HK) 18 , who showed that the total energy of N e electrons is a unique functional of the total electron density, and the minimum value of the total energy functional is the ground state energy, the density that yields the minimum energy is the ground state density. Kohn and Sham (KS) 19 later demonstrated that the problem of strongly interacting electrons can be mapped to a rather simple problem of single electron moving in an effective potential, if the exchange- correlation potential is known. The KS energy functional is: [ ] 3 () (()) ()() (()) (()) electronic ion hatree xc E nr Tnr V rnrdr E nr E nr=+ + + ? vv vv v v (2.2) where ()nr v is electron charge density, T is the kinetic energy of the system, () ion Vr v is 10 the ionic potential, 2 33 ()(') (()) ' 2' hartree enrnr E nr drdr rr = ? ? v v v vv vv is the long rang coulomb interaction between the electrons. (()) xc E nr v is the exchange-correlation energy. For a given nuclei configuration { } I R v , the KS Schr?dinger equation is: 2 2 (()) () () 2 eff i i i Vnr r r m ??? ?? ??+ = ?? ?? h v vv , (2.3a) where ( ( )) ( ( )) ( ( )) ( ( )) eff hatree xc ion V nr V nr V nr V nr= ++ vvvv , (2.3b) and 2 () () i i nr r?= ? v v . (2.3c) The KS theory provides a single-electron approach to exactly calculate the total energy of a many-electron system if the exchange correlation energy is known. To study real solids, further approximations, such as local density approximation 19 (LDA) or generalized gradient approximation 20 (GGA), are often used to approximate the exchange-correlation interactions of many electrons. Unless specified otherwise, results reported in this dissertation are obtained with the LDA. To self-consistently solve the KS equation, the electronic wave-functions are expanded with a chosen basis sets. The particular implementation of DFT theory is the Vienna ab initio Simulation Package (VASP) 21 , which adopts the planewave basis. And the j-th eigen-function of the k-point is written as: () ,, () iG k r j jk G k G rce? + ? + = ? vv v rvv v r , (2.4) 11 where G v is the reciprocal lattice vectors, k v is wave-vector. The KS equation in reciprocal space can be derived by plugging Equ. (2.4) into Equ. (2.3a): 2 2 ,' ',' ' (') 2 eff i iGG kG ikG G kG VGGc c m ? ? ++ ?? ++? = ?? ?? ? vv v vvv v v v h . (2.5) Solution can be found by diagonalizing the Hamiltonian matrix whose matrix element is 2 2 ',' (') 2 effGGkGkG H kG VGG m ? ++ =+ +? vv vvv v vvv h . Since terms with large kinetic energy 2 2 2 kG m + v v h are relatively small, all the terms whose kinetic energy larger than a chosen cutoff energy, ie., 2 2 2 cut kG E m +? v v h , (2.6) are set as zeros in order to reduce computation loads. To further reduce computation loads, only wave functions of (bond-forming) valence electrons are solved explicitly, and the effects of (environment-insensitive) core electrons are approximated with the so-called pseudo-potentials (PP). Two widely adopted pseudo-potential types are norm conserving pseudo potentials 22 and ultra-soft pseudo-potentials 23 (USPP). Unless stated otherwise, the studies reported here are performed with the USPPs. To evaluate the forces on the th i ion, we need to simply calculate the first-order derivatives the total energy with respect of the corresponding position I R v : II I I I I FE H H H H=?? =?? ? ? =??? ? ??? ? ?? ?? v . (2.7). The first term on the right is called Feynman-Hellmann (F-H) force, and the remaining 12 two terms are called Pulay forces. The Feynman-Hellmann theorem 24 states that Pulay forces vanish when the calculation has reached self-consistency and the basis set orthornomality persists and are independent of the atomic positions. Therefore, in the plane-wave implemented pseudo potential calculation, the F-H forces are exact forces if the plane wave basis is complete and the electron configuration has been relaxed sufficiently so that the calculated wave functions is as close as possible to that of the real eigenstates. 2.2 Atom-Scale Simulation and Modeling 2.2.1 Structural Models and Crystal Symmetries Solid is a macroscopic state of matter and a bulk solid usually consists ~10 23 atoms ?linked? with their neighbors with relatively strong chemical bonds. Because of the rigidity of the inter-atomic interactions, as the first order approximation, atoms inside a solid are considered to be ?fixed? in space. A crystal is a solid whose constituent atoms/molecules arranged in an ordered fashion. An ideal crystal can be described as a three-dimensional (3D) periodic array of points in real space (lattice) with a set of atoms arranged at each lattice point (atomic basis). Each lattice points are equivalent, and a vector that links two lattice points are called a lattice vector. Because of lattice periodicity, i.e. a crystal is indistinguishable when it is shifted by a lattice vector, we can define a unit-cell, which is a 3D box that can fill in all the lattice space by translation without leaving any gaps. A unit cell box is specified by three unit-cell vectors 123 ,,aaa vvv , with 123 (, , )aaa representing the lengths of lattice vectors and (, ,)? ?? representing the 13 angles between two lattice vectors. Any lattice vector can be indexed with 123 (, , )nnn : 123 11 22 33 (, , ) lattice R nnn na na na=++ r r rr . There are seven distinct crystal systems, i.e. cubic (c), hexagonal (hex), tetragonal (tet), rhombohedral (rhm), orthorhombic (orth), monoclinic (mono), and triclinic (tric). For example, a lattice with unit-cell vectors being 1 (,0,0)aa= r , 2 (0, ,0)aa= r , 3 (0,0, )aa= r is cubic lattice because 123 aaaa=== and 90???=== o . Similarly, a lattice with 1 (,0,0)aa= r , 2 (/2,3/2,0)aa a=? r , and 3 (0,0, )ac= r is hexagonal because 123 aaa= ? , 90??== o , and 120? = o . Additionally, three centering arrangements, i.e. base-centered (bac), face-centered (fc), or body centered (bc), exist. Using the combination of seven crystal systems and three lattice centering arrangements, we can categorize all 3D lattices into 14 Bravais lattice types. For example, a cubic lattice with the face-center arrangement belongs to the face- centered-cubic (fcc) Bravais lattice type, and its primitive unit-cell lattice vectors are: 1 (0,0.5 ,0.5 )aaa= r , 2 (0.5 ,0,0.5 )aaa= r , 3 (0.5 ,0.5 ,0)aaa= r , with 123 2 2 aaa a=== , and 60???=== o . The volume of the fcc primitive unit cell is ? of that of the conventional cubic cell. Three material systems studied in this dissertation, d-Si, Si 136 , and B1-MgO, happen to belong to the same fcc Bravais lattice type. The equilibrium phases of silicon and magnesium oxide at ambient conditions adopted the A4 (diamond type) lattice and the B1 (NaCl type) lattice respectively, and each contains two atoms per fcc unit-cell. Si 136 crystals adopt the type-II clathrate lattice with 34 atoms per fcc unit- cell (or 136 atoms per cubic cell), 14 The arrangement of atoms in the unit-cell can be uniquely specified with the locations of atoms inside a unit-cell using a set of internal coordinates (,, ) ii i uvw , where 1, , a iN= L , and a N is the number of atoms per unit cell. For example, the two Si atoms in a unit-cell of d-Si locate at (0, 0, 0) and (14,14,14) sites, while the one Mg atom and one O atom in the B1-structured MgO locate at (0,0,0) and (12,12,12) sites respectively. The position of each atom in the crystal can then be specified with its lattice indices 123 (, , )nnn and basis index i . For example, the position of i-th atom in the 123 (, , )nnn cell is 123 123 1 1 2 2 3 3 (, , ,) (, , ) ( ) ( ) ( ) lattice i i i i R nnni R nnn r n ua n va n wa=+=+++ rr r rr r . (2.8) As atoms of same type are indistinguishable, a group of sites occupied by the same type of atoms might be symmetrically equivalent, i.e. the crystal remains unchanged after a point-group (PG) (rotation) operation and/or a space-group (SG) (rotation and gliding) operation. The symmetry of a crystal is determined by the numbers and types of its PG and SG symmetry operators. For 3D crystals, there are 32 point groups and 230 space groups 25 . For example, both d-Si and B1-MgO crystals belong to the h O point group, while their space group symmetries are 3Fd m (#227 in the 230 SG list) and 3Fm m (#225 in the 230 SG list) respectively. Symmetry, being a dominant feature of crystals, should be and has been greatly taken advantage of in our calculations. We have implemented a FORTRAN code (FindSG.f90) to examine all the space group symmetry operators by combining all the possible point-group symmetry operations with possible glide plane and screw axis. We 15 also constructed an atom-atom mapping under each discovered symmetry operation, which is implemented within a code called Find1on1.f90. For example, after a given operation, atom #1 at position #1 moves to the position originally occupied by atom 5. We refer this as atom #1 is mapped to atom #5 with the specified operation. The output file of code Find1on1.f90 is called 1on1map.dat, which contains every possible symmetry operator and the atom-atom mapping under its operation. This information will be utilized in our derivation of irreducible single displacement moves in phonon calculations and paired displacement moves in anharmonicity calculation. Here we attached an example of 1on1map.dat file: 116 index of symmetry operator -1.0000000000 0.0000000000 0.0000000000 0.0000000000 -1.0000000000 0.0000000000 rotational operator 0.0000000000 0.000000000 -1.0000000000 0.8125000000 0.062500000 0.5625000000 glide operator 1 115 2 114 3 113 4 116 5 127 6 126 7 125 atom-atom map under current operator 118 14 119 13 120 16 121 11 122 10 123 9 124 12 125 7 126 6 127 5 128 8 16 The periodicity and symmetry of a crystal lattice can be detected by experimental techniques, such as X-ray diffraction. According to the Laue?s diffraction condition 26 or the equivalent Braggs? Law 27 , the atoms, being the scatters in the crystal, will scatter the incoming wave with wave-vector k v to all the directions, only those outgoing waves 'k v satisfying the Laue?s diffraction condition, 'kkG? = v v v , will interfere constructively. Here G v is a lattice vector in the reciprocal space, or momentum space (k-space), of a crystal. The reciprocal lattice vectors of a crystal are defined based on the corresponding real space unit-cell vectors: 23 1 123 31 2 231 12 3 312 2 () 2 () 2 () aa b aaa aa b aaa aa b aaa ? ? ? ? = ?? ? = ?? ? = ?? v v v (2.9) Obviously, 2 ij ij ab ???= v v for i and j being 1, 2, or 3. In the case of the fcc unit cell, the reciprocal lattice vectors are 1 4 111 (,,) 222 b a ? =? v , 2 4111 (, ,) 222 b a ? =? v and 3 4111 (,, ) 22 2 b a ? =? v , which form a body-centered cubic (bcc) structure with lattice constant 4 a ? . A reciprocal lattice vector is labeled with (,,)hkl index: 123hkl Ghbkblb=++ vvvv with ,,hkl?integers. Each lattice vector hkl G v in reciprocal space is associated with a set of parallel planes whose miller indexes are ()hkl in direct space. hkl G v is normal to these planes and the inter-plane spacing hkl d is equal to 2/ hkl mG? v , 17 where m is the common divider of (,,)hkl . In the case of elastic scattering, 'kk= vv , and the Braggs? Law can be easily derived from the Laue?s law as: 2sin hkl d ? ?= (2.10) where ? is the wavelength of the incoming wave, and ? is the incident angels with respect to the family of planes ()hkl . Studying reciprocal lattice vectors is of close relevance to the Fourier decomposition of lattice-periodic functions. For example, the potential ()Ur v for electrons in a crystal is a periodic over the lattice: () ( )Ur Ur R= + v v v , where R v is the direct lattice vectors. Decomposing it with a Fourier series based this translational symmetry, we have ik r ik R ik r kk kk eU e eU ??? = ?? vvvv ? () iG r G G Ur e U ? = ? v v v v v (2.11) From above equation, it is easy to derive that 1 ik R e ? = v v , or equivalently 2kR n??= v v where n is an integer. If k v is the reciprocal lattice vector, 123 khbkblb G= ++? v vvvv , then the above requirement is automatically fulfilled. Like direct lattices, reciprocal lattices also have the translational symmetry; any vector k v is equivalent to vector kG+ v v , therefore one unit cell in the reciprocal space will be enough to describe the lattice. The often used reciprocal unit cell is Wigner-Seitz primitive cell, which is also called first Brillouin zone (BZ), and often only k v within the BZ are calculated explicitly. When a infinitely large crystal is approximated with a large supercell containing 123cell NNNN=?? unit-cells ( 123 ,,NNN are the numbers of unit cells along lattice vector direction 123 ,,aaa v vv respectively) satisfying the Born-von Karman 18 periodic boundary condition, the number of k-points within BZ is cell N being the number of unit cells in direct space. When cell N ?? , for example in a real solids is N cell in order of 10 23 , the grids would be so dense that we can consider the distribution of k vectors inside is quasi-continuous. The evaluation of many physical properties involves integrating the periodic functions of k-vector in the BZ, therefore, the decent dense k-point sampling is crucial in order to yield accurate results, usually convergence test is desired in real calculation. Among all the uniform k-points sampled in BZ, many of them are dependent on the others due to the lattice point symmetry. All the independent k-points as well as their equivalent ones can be found by applying them the symmetry operations of the point group that the direct lattice belongs to, and only evaluation of the periodic function at these independent ones are necessary, which is central idea of Monkhorst-Pack 28 k-point sampling method. We have been able to calculate a lot of physical quantities, such as mode Gr?neisen parameters, phonon frequencies, mode group velocities, etc., over these reduced k-point set and then assign the calculated physical quantities to all the other equivalent k-points according to symmetry map. If the quantity is a scalar, we can simply assign the same value at an independent k-point to all the other equivalent k-points. However, if the physical quantity is a vector (for example, group velocity g V v ) or higher- order tensors, they might be different at two equivalent k points. This is because scalars stay invariant under the point group symmetry operation, while vectors and tensor usually don?t. 19 Finding the number of equivalent k-points and the map between the reduced k- point set and full coarse grids is done by code named as distinctK.f90. It reads in the unit cell information and 123 ,,NNN which determines the meshes, and will generate the coarse grid data, the reduced k-point set and also the map between them. Physical quantities over the coarse grids can be fully recovered by the code named Rebuilt.f90. In summary, unit-cell models are adopted to describe lattices and atomic basis of crystals. Based on its unit-cell lattice parameters and atomic arrangements inside the unit cell, a crystal is structurally categorized in term of space groups, point groups, and /or Bravais lattice types etc. International Tables for Crystallography is a useful reference for detailed crystallography information of lattices, such as the space-group symmetry operators and all the symmetry-allowed Wyckoff sites for atoms in lattices belonging to each of 230 space groups. 2.2.2 Energy Minimization and Equation of State In total energy calculation, cell shape and internal coordinates are allowed to relax if necessary and the total energy of the system will be calculated at each molecular dynamic step. The whole process is an iteration process. Iteration stops when the minimum energy is found, and corresponding structure is the equilibrium or optimized structure. The completion of a total energy calculation requires massive coding, and now it has become a standard technique and has been encoded into a package. There are many packages available, such as ABINIT 29 , VASP 30 , SIESTA 31 , QUANTUM ESPRESSO 32 etc. Our calculation were all carried out by VASP, which performs ab initio quantum- 20 mechanical molecular dynamics using density functional theory with pseudo-potential and plane wave basis set. For more details of this package, one can refer to VASP website 30 . The most important files in a VASP calculation are: INCAR, POSCAR, POTCAR and KPOINTS. And they can be customized according to the needs. INCAR contains information such as relaxation scheme, molecular dynamics iteration steps, energy cutoff for plane wave basis set, etc.; POSCAR contains the structural information, such as lattice vectors and internal atomic coordinates; POTCAR contains the pseudo-potential of nuclei; And KPOINTS contains the k-point sampling over the BZ for energy integration. Even though VASP allows the relaxation of the structure, including cell volume, cell shape, and internal coordinates, volume relaxation is not recommended in VASP when one tries to find the equilibrium structure and the minimum energy 30 . It is due to the fact that the arising of Pulay stress induced from the volume relaxation intends to result in an underestimation of the equilibrium volume, and then provide the minimum energy with rather large error bar. 30 However, this can be avoided by doing volume conserving calculations. The equilibrium structure can be found by doing constant- volume relaxation for multiple volume points with the same energy cutoff, and then fitting the energy-volume to an equation of state to get reliable lattice parameters and bulk modulus even with rather small energy cutoffs. There are several versions of equation of state; the one we adopted in our calculations is 3 rd order Birch-Murnaghan equation of state 33 . ()() 2 23 23 000 0 9 () ( ) 1 1(4 ')(1( ) )2 8 EV E KV V V K V V=+ ? +? ? , (2.12) 21 where ,,E VK and 'K are energy, volume, bulk modulus and its pressure derivative respectively. Those with subscript 0 indicate their equilibrium value. 2.2.3 Lattice Vibration Chemical bonds that hold atoms together to form a solid are NOT 100% rigid. In reality, atoms always oscillate around their minimal energy positions. As the amplitudes of atomic motion are relatively small compared to the inter-atomic spacing, we can expand the total energy tot E about its equilibrium value using Taylor expansion: 00 0 233 00 13 // 1,1 13 33 / ,, 1 1 ( ,..., ) ( ) ( ) 2 ( ,..., ) 1 () . 6 nn nn nn NN static N i i j xx xx iijj tot N N ijk xx ijk ijk EE Exx x xx x Ex x E xx x ho xx x == == = = ???? + ?+ ??+ ?? ? = ?? ? ??? + ??? ?? ?? ? (2.13) Here, static E is the minimal energy when all atoms are at their equilibrium positions. 0 n x and n x are the n th equilibrium and instantaneous atomic coordinates respectively, N is the total number of atoms, 3 represent (, ,)xyz three coordinate directions; and 0 nnn x xx? =? is the position deviation from equilibrium. The first order derivatives of the total energy 0 / () nn xx i E x = ? ? are associated with the forces on the atoms, which are zero at an equilibrium configuration. If we consider the atoms are connected with ?springs?, the 2 nd order derivatives 0 2 / () nn xx ij E xx = ? ?? can be interpreted as the harmonic 22 spring force constants. All higher-order derivatives represent the anharmonic spring forces among atoms. Within the harmonic approximation which neglects all the 3 rd and higher order terms in Equ. (2.13), the classical Hamiltonian of this strongly coupled anharmonic N- atom system is now simplified to a harmonic system: 233 00 1313 13 1,1 11 ( ,..., ; ,..., ) ( ,..., ) 22 NN i N N static N ij i j iij i p H xxpp Exx xx m ? == ?? = + + ?? ?? . (2.14) For a crystal, instead of dealing with force constant matrix of an infinitely large periodic system in real space, we can take advantage of lattice periodicity and define the dynamic matrices D of a unit-cell in the BZ of the reciprocal space: 0, ,0 ( ) exp[ ( )] ijl ij jl i l ij D qiqR mm ?? ?? ? =???? ? r r rr , (2.15) where q v is a q-point (k-point is conventionally used for electronic energy calculation, and q-point for vibrational energy calculation) in the BZ, ,1,,i j n= L are index of atoms in the unit cell, ,1,2,3? ? = are index of three coordinates. , ij mm are masses of atom i and j respectively, jl R r stands for the spatial location of atom j in the th l unit cell relative to the origin of the reference cell: 0 th . The summation is taken over all the unit cells under consideration, depending on the super cell size chosen for the system. Super cell should be big enough such that the interaction between the atom inside the super cell with all its images outside the super cell are negligible, but it is computationally unnecessary to use a cell larger than the one that already converges the 23 energy. Element of force constant matrix 2 0, 0 tot ijl ijl E x x ?? ?? ? ?= ?? is the second order total energy derivative with respect to the displacement of atom i in the reference cell along ? direction (0)ui ? and of atom j in th l unit cell along ? direction ()ujl ? away from their equilibrium positions. The eigen-frequency (?) and eigen-vector e of the ? th vibration eigen-mode at a reciprocal space q-point with wave vector q v can be calculated by solving the following eigen problem: 2 ()(,) (,)(,)qq q q? ?? ?=De e v vvv , (2.16) Thus the harmonic vibration in a perfect crystal can be considered as a superposition of normal modes (,)q ? v with various frequencies. Each normal mode represents a special type of vibration whose characteristic frequency is (,)qv? v , where q v is the wave vector whose magnitude is the reciprocal of the wavelength of the lattice wave, and v is the index that specifies the polarization of the wave. At each q v point, there are in total 3n number of normal modes, three of which are acoustic vibration modes and the rest 3n-3 are optical modes (n is the number of atoms in the unit cell). The vibration frequency of an acoustic mode approaches zero at the long wavelength limit, ie., 0q ? v , and these low-frequency acoustic vibration correspond to the sound waves in the lattice. We have implemented and further developed a real-space super-cell technique to calculate the force-constant matrix ? with first-principles methods. Additional LOTO splitting is then added for correction in ionic systems. 34 We first start with a fully relaxed structure, i.e. each atom is at its zero-force, minimal energy position, and then calculate 24 atomic forces on atoms when a single atom is displaced by a small, yet finite displacement. From Equ. (2.13), the ( , ,0)i? component of the forces can be derived as: 0 0 0 '' 23 0/ /' '00 0' 4 /'' ' 0'' 1 () ( ) 2 1 () . 6 nn nn nn ml tot tot tot ijljlk xx xx jl kl jliijl ijlk tot jl kl ml xx kl jl ijlklml EE E Fxx xxx xxx E xxx ho xxxx ? ???? ?????? ??? ?? ? ?? ???? == = ?? ? =? =? ? ? ? ? ? ? ? ? ??? + ???? ??? ??? ,(2.17) where l and 'l are index of unit cell, 0 indicates the reference unit cell, , , 1, ,ijk n= L are index of atoms in the unit cell with n being the number of atoms in the unit cell, ,, 1,2, 3or? ??= are three Cartesian component. The 4 th order anharmonicity is also included, but our algorithm will make it cancelled off, thus the energy actually is cut off at 5 th and beyond. According to Equ. (2.17), if we displace only one atom, say j, in the th l cell by ? , ie., jl x ? ? =+?, the ? component of the force felt by atom i in the reference cell is going to be: 23 0 0, 0,, 0,,, 11 26 I i i jl i jl jl i jl jl jl FAB ? ?? ??? ???? =?? ?? ? ? ? , (2.18a) Similarly, if we displace atom j in the th l cell by ?? , ie., jl x ? ? =??, the ? component of the force felt by atom i in the reference cell will be: 23 0 0, 0,, 0,,, 11 26 II i i jl i jl jl i jl jl jl FAB ? ?? ??? ???? =+? ?? ? + ? , (2.18b) In order to cancel off the fourth order anharomonicity contribution, we also displace atom j in the th l cell by 2+? and 2??, the forces are: 23 00, 0, 0,, 48 2 26 III i i jl i jl jl i jl jl jl FAB ? ?? ??? ???? =?? ?? ? ? ? , (2.18c) and 25 23 0 0, 0,, 0,,, 48 2 26 IV i i jl i jl jl i jl jl jl FAB ????? ??? =+? ?? ? + ? (2.18d) respectively. Now combine Equ. (2.18a) and Equ. (2.18b), we get 3 0 0 0, 0, , , 1 2 3 III i i i jl i jl jl jl FF B ? ? ?? ???? ?=???? ?. (2.19a) and combine Equ. (2.18c) and (2.18d), we get 3 0 0 0, 0, , , 8 4 3 III IV i i i jl i jl jl jl FF B ? ? ?? ???? ?=???? ?. (2.19b) From Equ. (2.19a) and (2.19b), we have the 2 nd order force constant as: 0000 0, 88 12 IIIIIV iiii ijl FFFF ???? ?? ?++? ?= ? . (2.20) As one can see from Equ. (2.20), each matrix element 0,ijl? ? ? can be achieved by extracting the HF force element 0i F ? in the total energy calculation for a super cell with atom j in the th l unit cell displaced by ? along ? direction, while others remain in their equilibrium positions. For a 33NN? force constant matrix, in principle, one should complete 26N? total number of total energy calculation in order to get the full matrix. Even for small systems, it is quite time consuming to complete all these calculations. However, the symmetry contained in the structure indicates that a lot of calculations are redundant; one should only carry out the calculation for those independence ones, and then recover all the others by proper symmetry operations. It will be shown that taking symmetry into consideration can greatly reduce the computation load and thus save a great amount of time. All the irreducible single atom displacement can be found by Code named Moves_Analysis.f90. 1on1map.dat, which contains the symmetry information for given 26 system, is the only input required for this code. This code generates a list of irreducible moves as well as the mapping information for reconstruction. The detail algorithm is explained in the next 5 paragraphs. In Cartesian setting, each atom can be displaced in 6 directions x+ , y+ , z+ , x? , y? , z? , and they are defined as the bases. In total we have 6 bases ( 6Nbasis = ), and they have the vector forms in the following and are labeled from 1 to 6 in the order as they appear: 100 10 0 0,1,0,0 , 1,0 0010 0 1 ? ??????? ?? ?? ? ??????? ?? ?? ? ? ??????? ?? ?? ? ??????? ?? ?? ? ? ??????? ?? ?? ? For a system containing N atom in the cell, there are in total 6N number of such single moves. During the initialization process, they will be labeled from 1 to 6N, and each of them has independence flag 1 indicating they are temporarily independent of one another. The index of each move directly determines the atom index and direction in the cell and vise versa. _(_1) _Index move Index atom Nbasis Index basis=??+, (2.21) _ _int 1 Index move Index atom Nbasis ?? = + ?? ?? , (2.22) __(_1)Index basis Index move Index atom Nbasis=???. (2.23) Table 2.1 One-to-one correspondence of the index for the displacement to the index of atom being displaced. 1 2 3 4 5 6 L j L 6N-5 6N-4 6N-3 6N-2 6N-1 6N atom #1 atom #N 27 The one-to-one correspondence between the index of moves for single-atom displacement and the atom itself can be seen above. We are aiming to find independent single ones among these 6N moves in order to reduce computation load. We start a loop over all the moves first, for th j move with index _imove(Note, we use _i and _f to indicate initial and final respectively), we found its atom index _iatom and basis index _i basis by using Equ. (2.22) and Equ. (2.23), and then run a loop over all the operators. Under each operator in 1on1map.dat, we can find the atom index _f atom which is equivalent to _iatomaccording to the atom-atom mapping in the file, and also, the basis index _f basis can be found by operating the operator to the initial basis: __f basis Operator i basis= ? . (2.24) With the knowledge of _f atom and _f basis , move index _f move can be determined by Equ. (2.21). If index _f move is larger than _imove and its independence flag is still 1, then we updated its independence flag to 2, and record the index _imove as well as the operator that correlates them. Of course, if index _f move is smaller than _imove or its independence flag is 2, we don?t need to take action to it, since it has been dealt with or it is no longer independent. By doing a loop over all the single moves and then a loop over all the operators, we have finished the update of independence flag for each move, and those remaining as 1 are the independent ones. With the independent move index, we know what atom and in which direction we should displace to carry out the F-H force calculation. Of course, for those dependent ones, their dependence and the corresponding operator have been recorded so that we can reconstruct the full force matrix later. 28 How far should one displace the atom from its equilibrium position in the force constant calculation? We intend to choose as small displacement as possible, since the irreducible moves are arrived by assuming the displacement does not change the crystal symmetry, which, in fact, does. Most often a displacement of 0.2% work rather well and the convergence test is desired for the exceptions. Now we have got the force constant for all the irreducible moves. However, the dynamics matrix is built upon the full force constant matrix (see Equ. (2.15)). Thus a construction of full force constant matrix based on the irreducible move results shall be performed. This is done by our reconstruction code named Rebuild_FullFM.f90. Let?s explain the major algorithm for force reconstruction for single moves. We start with a loop over moves, for a given move index _f move , if its independence flag is 1, we don?t need to do anything since the forces are calculated directly; but if the flag is 2, it means that this move is dependent on some other move with index _imove through an operator. Index _imove must be smaller than index _f move since this is the way we saved them in Move_Analysis.f90 code, and because of this, the forces on all the atoms due to displacement _imove must have been available first. Knowing the forces on all the atoms due to the displacement _imove, we can easily find the forces due to _f move by symmetry operation. Suppose the forces exerting on the atom _iatom due to the displacement _imove is _iforce , since _f move is dependent upon _imove through an Operator , the force _i force felt by _iatomin _imove calculation will become the force _fforce felt by _f atom in _f move . The Operator is used to figure out what 29 _fforce and _f atom are. _f atom is the atom that the _iatom corresponds to in the atom-atom mapping under the current Operator , and _f force can be calculated by: __fforce Operator i force= ? . (2.25) Since the atom-atom mapping under each operator is one on one, a loop over all the atoms in the cell will generate the full force matrix due to the _f move . By far, we have been able to use symmetry to figure out irreducible moves for force constant calculation and also reconstructed force constant matrix. It is obvious that symmetry has greatly reduced our calculation load. In addition, it also reduces the numerical error in the sense that recovered forces from symmetry are the same as its dependent if they are indeed the same, but if we carry out F-H force calculation for both, we might end up with slightly different force vectors between these calculations purely due to numerical error. This is especially important if there is a cancel-out effect. We have further developed a paired-displacement algorithm to evaluate third- order anharmonicity tensor A with element 3 0, , ' 0' ijlkl ijlkl E A x xx ??? ??? ? = ??? . We propose to displace a pair of coordinates in the following fashion. Instead of displacing one atom, we now displace a pair of atoms at the same time. If we displace both atom j in th l unit cell and atom k in ' th l unit cell by ? , ie., ' , jl kl xx ?? ? =+?? =+?, the ? component of the force felt by atom i in reference cell is going to be: 2 0 0, 0,' 0,, 0,',' 0,,' 3 0,,, 0,',',' 0,,',' 0,',, 1 ()( 2) 2 1 (2 6 V i i jl i kl i jl jl i kl kl i jl kl i jljljl i klklkl i jlklkl i kl jljl FAAA BB B B ????? ???? ??? ???? ???? ???? ???? =?? +? ?? + + ? ? ++ + ? , (2.26a) 30 Similarly, if ' , jl kl xx ?? ?=+??=??, then 2 0 0, 0,' 0,, 0,',' 0,,' 3 0, , , 0, ', ', ' 0, , ', ' 0, ', , 1 ()( 2) 2 1 (2 6 VI i i jl i kl i jl jl i kl kl i jl kl i jl jl jl i kl kl kl i jl kl kl i kl jl jl FAAA BB B B ????? ???? ??? ???? ???? ???? ???? =?? ?? ?? + ? ? ? ?+ ? ? , (2.26b) if ' , jl kl xx ?? ?=???=+?, then 2 0 0, 0,' 0,, 0,',' 0,,' 3 0,,, 0,',',' 0,,',' 0,',, 1 ()( 2) 2 1 (2 6 VII i i jl i kl i jl jl i kl kl i jl kl i jl jl jl i kl kl kl i jl kl kl i kl jl jl FAAA BB B B ????? ???? ??? ???? ???? ???? ???? =?? +? ?? + ? ? ? ?+? + ? , (2.26c) and if ' , jl kl xx ?? ? =?? ? =??, then 2 0 0, 0,' 0,, 0,',' 0,,' 3 0, , , 0, ', ', ' 0, , ', ' 0, ', , 1 ()( 2) 2 1 (2 6 VIII i i jl i kl i jl jl i kl kl i jl kl i jl jl jl i kl kl kl i jl kl kl i kl jl jl FAAA BB B B ????? ???? ??? ???? ???? ???? ???? =+? +? ?? + + ? + ++ + ? . (2.26d) Note that we consider j k= being a generalized pair of atoms. Combine Equ. (2.26a) and Equ. (2.26d), we get 22 2 0 0 0,, 0,',' 0,,' 2 V VIII i i i jl jl i kl kl i jl kl FF A A A ? ? ??? ??? ??? +=? ?? ?? ?, (2.27a) Combine Equ. (2.26b) and Equ. (2.26c), we get 22 2 0 0 0,, 0,',' 0,,' 2 VI VII i i i jl jl i kl kl i jl kl FF A A A ? ? ??? ??? ??? +=? ?? ?+ ?, (2.27b) From Equ. (2.27a) and Equ. (2.27b), we will have: 000 0 0, , ' 2 4 V VI VII VIII iii i ijlkl FFFF A ???? ??? ?++? = ? . (2.28) 31 Special case like ,,'j kl l? ?=== is considered separately. 0, ,ijljl A ? ?? can be determined by forces from single atom displacement. Combine Equ. (2.18a) and Equ. (2.18b), we get 2 00 0,, III ii ijljl FF A ?? ??? + =? ? . (2.29a) and combine Equ. (2.18c) and Equ. (2.18d), we get 2 00 0,, 4 III IV ii ijljl FF A ?? ??? + =? ? . (2.29b) From both Equ. (2.29a) and Equ. (2.29b), we have 0000 0, , 2 3 I II III IV iiii ijljl FFFF A ???? ??? +?? = ? . (2.30) Equ. (2.28) and Equ. (2.30) have clearly related the 0, , 'ijlkl A ? ?? terms to the F-H forces due to displacement of single and pair atoms which can be calculated through ab initio method. Each element of 0, , 'ijlkl A ? ?? can be achieved by extracting the F-H force element 0i F ? in the total energy calculation for a super cell with atom j in the th l unit cell displaced by ? , ??, 2? or 2? ? along ? direction and atom k in the 'th l unit cell displaced by ? , ?? , 2? or 2? ? along ? direction, while others remain in their equilibrium positions. Compared with force constant matrix, the full tensor evaluation requires 2 6N C more calculation due to pair moves. Again we took advantage of crystal symmetry. Moves_Analysis.f90 generates not only independent single moves for force constant calculation, but also independent pair moves for 3 rd order anharmonicity calculation. 32 Single independent move analysis and the full force matrix reconstruction have been explained in detail previously in this section; here we add the explanation for the derivation of irreducible pair moves as well as the reconstruction. First of all, we assign each pair moves ( , )i j a unique number which is defined as : _(,)(1)Index pair i j iNmovesj= ?? + (2.31) where i, and j are index of moves, as defined in Equ. (2.21), and Nmoves Nbasis Natom=?, Nbasis is 6 as defined earlier, and Natom is the number of atoms in the super cell. Since the forces induced by displacement of pair ( , )i j and pair (,)j i are the same, we only need to consider those terms with i j? , ie., the upper triangle of Table 2.2. Table 2.2 One-to-one correspondence of index for the pair displacement to that for the single atom displacement. (1 1) 6 1N?? + 2 L L 61N ? 6N (2 1) 6 1N?? + 62N + L L 26 1N? ? 26N? M M M M M M M M M (1)6iNj? ?+ M M M M M M M M (6 1) 6 1NN?? + L L L 661NN? ? 66NN? As one can see, the pair index number, acts like an identity, uniquely defines the pair moves ( , )i j . Like we do to single moves, we also assign each pair move an i-move j-move 33 independent flag and they are initially assuming to be 1, meaning they are independent of one another. First, we start with a loop over i-move from 1 to 6N, and another loop over j-move from i-move to 6N, note that only those equal to or larger than i-move are considered. Then, we add in another loop over all the operators in 1on1map.dat. For each operator, we need to find pair move (, )i j the symmetry equivalence (,)mn under its operation. The process of finding m-move for i-move and n-move for j-move is the same as that we have addressed for single atom displacement when we construct the force constant matrix. Basically, we figure out correspondent displaced atom (i-atom and j- atom) and direction (i-basis and j-basis) for pair move (, )i j ; then find correspondent pair atom (m-atom and n-atom) for (i-atom and j-atom), and displacement direction (m-basis and n-basis) for (i-basis and j-basis), and thus get the index move m-move and n-move (,)mn . Recall that only pairs i j? are considered, we need to swap the order of m and n if mn> . If _(,) _(,)Index pair m n Index pair i j> , we will update its independence flag _(,)Independence pair m n to 2, and record the symmetry operator _ ( , )OP pair m n for reconstruction use. There is one special case we need to consider separately. It is when pair move happens on the same atom while in opposite directions, it means that the atom is not displaced at all, therefore forces felt by all the atoms should be zero. We set independence flag as 0 for such pair moves. After we complete all the loops, all the independent pair moves will have independence unchanged, i.e., the independence flag remains as 1, and we will only carry out H-F force calculation for those independent pair moves. Forces due to the dependent 34 ones can be derived from symmetry and it is also done by the reconstruction code named Rebuild_FullFM.f90. This code serves two purposes, it rebuilds the forces due to all the single moves based on the knowledge of those from the irreducible singe moves, and also the forces due to all the pair moves based on those from the irreducible pair moves. The force construction of pair moves is similar to that of single moves except that we now record everything for _Index pair instead of _Index move . If the independence flag is 0, then all the forces will be set as zero. If the independence flag is 1, we will keep the F-H forces unchanged; if the independence flag is 2, we rebuild all the forces with the knowledge of mapping information and operator index. As defined earlier, tensor A is a 333NNN? ? matrix. With N being the number of atoms in the supercell, usually in the order of 100, we are dealing with a tensor containing around 6 2.7 10? number of elements. Since individual ijk A term represents the response of the system to the displacement of atoms in the supercell ( , , )i j k , and the crystal symmetry indicates that such response can be due to other possible combinations (', ', ')i j k . Therefore, the full tensor A can be represented by only rather small number of independent ijk A terms after incorporating with the symmetry. How to find those independent ijk A terms and then the symmetry map from those independent terms to full tensor has been implemented in the code named Get_Aind_map.f90. In d-Si system, 128- atom supercell was used in calculation, our calculation with Tersoff empirical potential shows that there are only 28 nonzero ijk A terms with the approximation that only terms involving 2 nd nearest neighbors are considered. 35 Another important issue is the sumrule enforcement. Numerical calculation might not preserve the intrinsic physics as desired, especially translational symmetry of the crystal. Those divergences could be avoided by the sumrule enforcement. For 3 rd order anharmonicity, its sumrule is enforced by setting the following: 0, '', '''' '''' '' ' ' 1,2,3 0 kklklkl kl k l Ar ?? ? ? ?= = ??? , (2.32) where kl r ? is the th ? component of the vector locating the th k lattice atom in the th l unit cell. If ? ?= , Equ. (2.32) corresponds to the sumrule for self terms, and if they are different, it then corresponds to the sumrule for cross terms. Both should be enforced in order to avoid the divergence due to the numerical error. Figure 2.1 shows that the sumrule enforcement has indeed improved our calculation. We plotted the Gr?neisen dispersion using two difference methods: Finite Difference Approach (FDA) and Anharmonicity Approach (AA) (will see in section (3.1.2)). Before we confined the ijk A data to satisfy the sumrule, the Gr?neisen data from AA match well with FDA results at all the q-points except those transverse modes near ? point. They basically diverge upon approaching gamma point. On the other hand, the sumrule enforced ijk A data has successfully eliminated the divergence near gamma point; the whole set of data sit right on top of the FDA data. 36 Figure 2.1 Comparison of Gr?neisen parameters calculated from FDA (blue circle) and AA, which include both sumrule un-enforced (red diamond) and sumrule enforced (black square) results. 37 CHAPTER 3 STATISTICAL THEORIES OF PHONONS 3.1 Quantum Theory of Lattice Vibration 3.1.1 Phonons: Harmonic Frequencies and Group Velocities Adopting proper normal modes of lattice vibration () v Qq v and () v Pq v as a new basis set: 3 * 1 () (,) () n v Qq meqvxq ?? ? ?= = ? v vv v (3.1) 3 * 1 () (,) () n v Pq me qvp q ?? ? ?= = ? v vv v (3.2) where ( , )eqv vv is the eigen-vector of normal mode ( , )qv v , ()x q ? v is the Fourier transformation of atomic coordinates x ? , m ? is the mass of the atom. The Hamiltonian of harmonic lattice vibration can be written as a summation of 3N simple harmonic oscillators: 3 22 2 1 11 () ( () (,) ()) 22 n vv v H qPq qvQq? = =+ ? vv v (3.3) 38 Each normal mode (,)qv v is characterized with its eigen-vector representing the vibration pattern and the corresponding vibration frequency () v q? v . The square of the amplitude of a vibration mode is proportional to the vibration energy of the given vibration mode. Within the classical mechanics, such vibration energy is a continuous variable, while quantum mechanics only allows a set of quantized energy level for the vibration mode with frequency ?: 1 () 2 n n? ?=+h , where 0,1, 2 .n = L (3.4) Total harmonic vibration energy is the summation of energy of all the normal modes, ie.,: 11 () () 2 lattice v v qv q E qnq N ? ?? =+ ?? ?? ?? v v v h . (3.5) where h is the Planck constant, 1 () 2 v q? v h is the zero energy of the normal mode, () v nq v is the vibration energy level index of the normal mode. The particle representation of this quantization of lattice vibration energy is often called phonons, and () v nq v can be correspondingly interpreted as the number of phonons occupied at the normal mode (,)qv v . At thermal equilibrium, a solid now can be treated as a system containing a collection of independent (or weakly coupled) phonons. The group velocity of a phonon is defined as: (,) (,) gq vqi qi?=?v v vv (3.6) Since 2 11 (,) (,) (,) () (,) 22 qq q qi qi eqi q eqi?? ?? ?=? =? Dvv v v , the mode group velocity can be derived using the Feynman-Hellmann theorem : 39 1 (,) (,) () (,) 2 gq v qi eqi q eqi ? =?Dv vv vv v v (3.7) 3.1.2 Quasi-Harmonic Approximation and Mode Gr?neisen Parameters When lattice anharmonicity contribution to the vibration energy is relatively small, the effects of anharmonicity can be approximated with the quasi-harmonic approximation (QHA), where the harmonic oscillator model is valid, yet the oscillation frequencies become volume-dependent. The mode Gr?neisen parameter, defined as , , ln ( ) ln qv qv dV dV ? ? ?? v v . (3.8) is introduced to quantify the effect of anharmonicity on each phonon mode. Apply the Feynman-Hellmann theorem to Equ. (3.8), we can evaluate each mode Gr?neisen parameters with the calculated volume derivatives of dynamical matrix ()/ q ddVDv and the eigen-vectors of phonon modes: , 2 , () () 2 q qv v v q d V qq dV ? ? ? =? D ee v v v v v . (3.9) where V is volume of the unit-cell, ,q ? ?v and () v qe r v are eigen-frequency and eigenvector of ?-th eigen-mode at the q v point in the BZ respectively, and ()/ q ddVDv is the volume derivative of the dynamical matrix, and the first order derivative is approximated by finite difference method. On the other hand, volume derivatives of dynamical matrix can also be evaluated using the 3 rd order lattice anharmonicity based on a perturbation approach: 35 40 *' (') , 0, ' ', '' '' '' '' ' ' '' '' ' 1()( 6() kk iq R l vv qv k kl k l k l kklklv kk eqeq A er q MM ?? ?? ? ? ??? ? ? ? =? ???? v v v v v v , (3.10) where kl r ? is the ? component of the vector locating the th k lattice atom in the th l unit cell. Equ. (3.10) further provides an indirect way to verify the calculated 3 rd order lattice anharmonicity tensors by comparing the ,qv ? v results obtained using Equ. (3.9) and Equ. (3.10). Unlike 2 nd force constant matrix, which has phonon calculation to examine whether or not those calculated elements are reasonable values, 3 rd order tensor A is not known to have direct relation to any measurable physical quantities. And the precise representation of both phonon modes and anharmonic forces are required in order to calculate accurately the phonon life time and the intrinsic lattice thermal conductivity. The consistency between two methods will definitely provide us confidence of applying our anharmonicity data to compute phonon-phonon scattering rate in phonon life time and thermal conductivity calculation. 3.1.3 Anharmonicity Induced Phonon-Phonon Scatterings Any lattice imperfection (such as finite-size grain boundaries, defects, and isotopes etc) or anharmonicity might cause interactions among phonons. In this work, we focus on the effects of lattice anharmonicity. As the first step, we ignore 4 th and higher order anharmonicity, and investigate only the three-phonon scattering mechanism. Four- phonon scattering may be important at high temperature, which will be explored in future work. 41 In general, three-phonon processes can be categorized into two types. One is the annihilation process, in which one phonon is annihilated; the other one is creation process, in which one additional phonon is created. Two types of scattering scheme in 3-phonon processes are shown in Figure 3.1. Figure 3.1 Schematic figures for two types of 3-phonon processes. The concept of scattering carries with it the implication that the rate of scattering is relatively small and for the interaction that produces scattering to be regarded as a small perturbation 36 . Scattering rate, which determines the phonon life time, can be obtained using standard quantum-mechanical time-dependent perturbation theory as embodied in Fermi Golden Rule 37 . 2 2 () f ifi PfHi ? ? ??=? ? h (3.11) where i is initial state with energy i ? and f is the final state with energy f ? under perturbation H? . In the case of lattice anharmonicity induced 3-phonon scattering Type I q q' q'' Type II q q' q'' 42 mechanism, { }(,) i inqv= , { } (,) f f nqv= , and the perturbation is the 3 rd order lattice anharmonic energy, ie., 0, , ' ' 1 (0) ( ) ( ') 6 ijlkl ijlkl H Axixjlxkl ??? ? ? ? ??? ?= ? ? ? ??? (3.12) With normal coordinates as the new basis, H? can be written as: ' '', , , ''' 11 ( , ') ( ) ( ') ( '') 6 qq qg i j k i j k qqq H A qqX qX qX q N ??? ? ? ? ??? ? ++ ?= ???? vv vv vvv v vvvv (3.13) where ('') ,, 0, , ' .' (, ') iql ql ijk i jlkl ll Aqq A e ??? ? ? ? ? ?+ ? = ? v v v v v v (3.14) and (,) () (,) v v eqv Xq Qqv m ? ? = ? v v (3.15) Combine Equ. (3.5), (3.11), (3.13), and (3.14) and (3.15), we have transition rates '' ,' q qq P v v v for type I scattering event and ', ''qq q P v v v for type II scattering event as: '' '' ,' ,' ( ) ( ')( ( '') 1) qoo o q qq i j k qq Pnqnqnq S=+ v v v vv v vv (3.16) where 2 2 '' ,' 2 12 ( , ') ( ) ( ') ( '') 24 (,)(',)('',) ( ) ij k q qq A qqe qe qe q h S Nqiqjqk mmm ??? ? ? ? ??? ??? ?? ? ? = ? v vv v vvv vv v ' '', ( (,) (',) ('',) qq qg qi q j q k? ?? ?? +? ?+? v vvv v vv hh h and ', '' '. '' ()((')1)((')1) qq o o o qq qij k Pnqnq nqS=++ vv vv v vv (3.17) 43 where 2 2 '. '' 2 12 ( , ') ( ) ( ') ( '') 24 ( , ) ( ', ) ( '', ) ( ) ij k qq q A qqe qe qe q h S Nqiqjqk mmm ??? ? ? ? ??? ??? ?? ? ? = ? vv v v vvv vv v ''', ( ( , ) ( ', ) ( '', )) qq qg qi q j q k? ?? ?? ?? ??? v vvv v vv hh h 3.2 Equilibrium Thermal Properties 3.2.1 Statistical Ensemseble Theory of Independent Phonons As discussed Section 3.1, the thermal excitation of a perfect crystal is now described as a system of independent phonons, which are indistinguishable quantum particles that obey Bose-Einstain statistics. As a result, ,qv nv , the equilibrium average number of phonon at the (,)qv v mode, is found to be , , / 1 (,) 1 qv B qv kT nT e ? ? = ? v v h , and the equilibrium averaged vibration energy is ,, , 11 (,) 2 vibration q v q v qv q E nT N ?? ?? =+ ?? ?? ? vv vr h . The same result can be derived based on the equilibrium ensemble theory. Within the harmonic approximation, the canonical partition function / jB EkT j Ze ? = ? can be easily derived as: , 12 3 , 12 /2 ()/ / , 00 0 1 Bkv NB Bkv N kT kT kT kv nn n e Ze e ? ?? ? ? ? ?? ? ?+++ ? == = ? ? ?? ? v v h L v h L . (3.18) Any other thermodynamic quantities, such as the total energy tot E , entropy S , and Helmholtz free energy F of a harmonic crystal, can be derived from the partition function Z. For example: 44 , , , / ,, 11 1 2 1 qv B qv tot static vibration static q v kT qv qvqq EE E E NNe ? ? ?=+ =+ + ? ??v v v h vvrr h h , (3.19) , , , , ln 1 1 (ln ( ) ) ln(1 ) (1) qv B qv B qv kT BV B qvq kT Z Sk ZT k e TNT e ? ? ? ? ?? ? ?? =+ = ? ? ? ?? ? v v h v h vr h , (3.20) ( ) , / , ,, 11 ln ln 1 2 qv B kT B static B static q v qv qvqq kT FE kTZE e NN ? ? ? =? =+ + ? ?? v h v vvrr h . (3.21) As shown in Equ. (3.21), the inputs required to calculate the thermodynamic potentials are the static total energy () static E V and harmonic phonon spectra ,qv ? v , both can be calculated with the first-principles total energy theory discussed in previous section. 3.2.2 Macroscopic Thermal Properties within the QHA To connect the first-principles calculated thermodynamic potentials with experimental measurements, we can further evaluate the thermal properties based on the macroscopic thermodynamic theories. For example, to predict the iso-thermal equations of state, we simply perform volume derivative of the Helmholtz free energy ( T F P V ? ?? =? ?? ? ?? ): , , , 11 (( , ) ) 2 qv static qv qvq d dE PnT dV N dV ? ? ?? =? ? + ?? ?? ? v v vv h . (3.22) Based on the QHA and the notation of mode Gr?neisen parameters ( ,qv ? v ) introduced in Section 3.1.2, we can re-write the Equ. (3.22) as: ,,, , () 11 1 (,) () (,) (( ,) ) 2 static static thermal q q q qvq dE V PTV P V P TV n T dV V N ? ?? ? ??=+ =? + + ? v vv vv h .(3.23) 45 With predicted P(V,T) equation of state, we can readily calculate all the measurable thermal quantities, such as isothermal compressibility: 1 T T V VP ? ? ?? =? ?? ? ?? , (3.24) or coefficient of thermal expansion ( ) () () 11 V T V P T PT V PT VT VPV ?? ?? ? ?? ==? =? ?? ??? ?? . (3.25) The specific heat capacity of the system is the derivative of the total energy: , , , ,,/ 2 /2 , 2 / , 11 1 1 1 qv B qv B qv B tot VqvqvkT qv qvVVqq kT kv B kT qvqB E Cc TN N Te q ke NkT e ? ? ? ? ? ? ? ?? ?? ? ? == = ?? ? ? ??? ? ? ?? = ?? ???? ? ?? ?? ? v v v vv h vv h v hv v h v h , (3.26) where , , 2 /2 ,, ,, 2 / (,) () 1 qv B qv B kT qv qv qv V qv B kT B nT e ck TkT e ? ? ?? ? ? ? ? ?? == ?? ? ? ??? ? ? ? v v h vv vv h h h (3.27) is the mode heat capacity. At low temperature B kT?hnull , V C is proportional to 3 T , it corrects the classical behavior at low T; and at high temperature B kT?hnull , it saturates to a constant which is close to Dulong-Petit limit of 3/ B katom. 38 Another important thermodynamic parameter that links thermal and mechanical properties is the (bulk) Gr?neisen parameter, which is defined as: ln () ( ) ln gV S VT PTV V EVC ? ? ? ? ? ?=?= ?? . (3.28) Plug Equ. (3.25), (3.26), and (3.27) into Equ. (3.28), we get 46 , , () q g Vqv qVV c VP CT C ? ? ? ? == ? ? v v v . (3.29) This result shows that the bulk Gr?neisen parameter g ? is simply the heat capacity weighted average mode Gr?neisen parameter ,qv ? v . 3.3 Non-Equilibirum Thermal Transport Theories 3.3.1 Kinetic Transport Theory At thermal equilibrium, thermal energy flows insides a system randomly from one part to another at the microscopic level. Yet, there is no net heat current at the macroscopic level, i.e. the statistically averaged heat current is exactly zero. On the other hand, the Fourier?s Law states a net heat current J r appears wherever a spatial gradient (? r ) of temperatures (T) exists, and the amount of the current is proportional to T? r : JT?=?? rr t . (3.30) Here ? t is the non-equilibrium transport coefficient called thermal conductivity. Any thermally excitatble particles inside a solid can contribute to the thermal currents. For semiconducting/insulating materials, lattice vibration is the dominant factor to the heat conduction, and phonons are considered as the the main carriers of thermal energy in this case. According to the simple kinect transport model, the microscopic heat current can be expressed as : , 1 (,) (,) (,) g qiq J nqi qiv qi N ?= ? v v r rv vv h (3.31) 47 where the (,) (,)nqi qi? vv h term represents the thermal energy of a phonon mode (,)qi v , and (,) g vqi r v is the corresponding group velocity of those phonons. Within this kinetic model, the lattice thermal conductivity can be simply expressed as: , 1 (,) (,) (,)(,) Vg g qiq cqiv qiv qi qi N ?? ?? ??= ? vr v vvv (3.32) where (,) V cqi v , (,) g vqi ? v , and (,)qi? v are the heat capacity, ? component of group velocity , and phonon life time of the normal mode (,)qi v respectively. Heat capacity contribution from each mode is given by Equ. (3.27) and group velocity of each phonon mode is defined in Equ. (3.6). The challenge in kinetic method is to calculate the phonon life time for each normal mode. In pure harmonic crystal, phonon life time is infinity; the phonon spectrum consists of several pure ? -function peaks. The anharmonic forces in the crystal causes the phonon spectrum to shift and broaden, leading to finite phonon life time, which is the reciprocal of the full-width at half-maximum (FWHM) of the peak. Those line shifts and line width can be measured by a variety of experimental techniques, such as neutron scattering. Phonon life time can be determined by phonon scattering rate. The variation of the phonon occupation number with time happens in the diffusion process, scattering process or by adding external field. () () () scattering diffusion field dn dn dn dn dt dt dt dt =++ (3.33) 48 If there is no external field in a thermal equilibrium system, 0 dn dt = and () 0 field dn dt = , thus () () scattering diffusion dn dn dt dt =? . Phonon life time ? , the time gap between two consecutive phonon collisions, is introduced in the relaxation time approximation which states that: () o diffusion dn n n n dt ? ? ? ? ?= (3.34) where n is the phonon occupation number at any condition, while o n is the equilibrium phonon occupation number, which follows the well-known Bose-Einstein distribution: 1 1 B o kT n e ? = ? h . (3.35) There are several phonon scattering mechanism that contributes to term () scattering dn dt , here we only consider the 3 rd order anharmonicity induced 3-phonon process in the phonon- phonon scattering mechanism, assuming other contributions are negligible. Under the two approximations just made, Equ. (3.33) for each normal mode characterized by q-vector q v and polarization index i reduces to: 3 () () () ii phononi nq dnq qdt? ? ? ?? =? ?? ?? v v v (3.36) Due to the fact that the energy excitation in oscillators can only happen in two consecutive energy levels, change of phonon occupation number with respect to the time in the 3-phonon process can be written as: 49 ()()()()()() ()()()()()() 3 '' ,' ', '' ''' () ( ) 1 ( ') 1 ( '') ( ) ( ') ( '') 1 1 ( ) 1 ( ') ( '') ( ) ( ') 1 ( '') 1 2 i phonon q ij kijk qq qq ijkij k q dn q dt nq nq nq nq nq nq S nq n q n q nq nq n q S ? ?? ?? ?? ?? ?? ++ ? ++ ?? ?? = ?? +?+ ?? ?? v vv vv vv v v vv vvvv vvvvv v (3.37) where the summation is taken over all the possible combination of two other normal modes (',)q j v and ( '', )qk v other the one under consideration (,)qi v in the 3-phonon scattering process. () i nq v , (') j nq v and ( '') k nq v are the phonon occupation number of mode (,)qi v , (',)q j v and ( '', )qk v respectively. The first term represents the type I scattering, negative sign in the first bracket indicates a reversed type I scattering; similarly to second term which represents the type II scattering, 1/2 comes from the over- counting of ( ', '')qq vv pair since they are equivalent in the type II scattering mechanism. Energy conservation in both processes shows that the following relationships hold: ()( )()( )( )( ) ( ) 1 ( ') 1 ( '') ( ) ( ') ( '') 1 oo oooo ij kijk nq n q n q nq n q n q++ = + vv vvvv for type I scattering, and ()()()( )( )( ) ( ) 1 ( ') ( '') ( ) ( ') 1 ( '') 1 ooooo o ijkij k nq n q n q nq n q n q+?++ vvvvv v for type II scattering. We now adopt the Single-Mode Excitation Approximation (SMEA), which states that () i nq v is the only one that deviates from the equilibrium value () o i nq v , i.e.: () () () (') (') ('') ('') o ii i o jj o kk nq n q nq nq n q nq n q =+? = = v vv vv vv (3.38) 50 After considering energy conservation and further neglecting the 2 nd and higher order of n? , the first term in the Equ. (3.37) is reduced to ( ) '' ,' ( ) ( '') ( ') ooq ik j nq n q n q S?? ? v vv v vv and the second term () ', '' 1 ( ) ( ') ( '') 1 2 oo q ij k nq n q n q S?? ? + + v v v vv v . Now we have: ()() '' ', '' ,' '''3 () 1 ( ) ( '') ( ') ( ') ( '') 1 2 ooqoo qi ikj jk qqphonon dn q nq n q n q S n q n q S dt ? ? ??? =? ? ? ? + + ? ??? ???? ?? vvv vv v vv v v vv .(3.39) Note that the common term () i nq? v has been factored out since it is independent of the summation. Combine Equ. (3.36) and (3.39), we have reached an expression for mode phonon life time: ()() '' ', '' ,' ''' 1 1 ( '') ( ') ( ') ( '') 1 () 2 ooqoo q kj jk i qq nq nqS nq nq S q? ?? =??++ ?? ?? ?? vvv vv v vv vv vv v (3.40) where o n can be calculated by knowing the phonon spectra, and '' ,' q qq S v vv and ', ''qq q S vv v are related to both phonon spectra and lattice anharmonicity through Equ. (3.16) and Equ. (3.17) respectively. In summary, phonon life times can be calculated by knowing phonon spectra and lattice anharmonicity. Combining mode phonon life time, mode heat capacity and mode group velocity, we have the mode contribution to the thermal conductivity. A decent k- grid sampling is desired to get an accurate result for thermal conductivity. Tetrahedron method provides a more accurate result than simple summation in terms of BZ integration. 51 3.3.2 Green-Kubo Formula Non-equilibirum transport properties are also associated with the fluctuation phenomena at thermal equilibirum. According to the Fluctuation-Dissipation Theorem (FDT) of Green-Kubo 39 , thermal conductivity (? ) is associated with the autocorrelation function of the fluctuated heat current at thermal equilibrium: 2 0 () (0) B JtJ dt kT ?? ? ? ? ? ? = ? (3.41) where ,,,xyz? ? = , ? is the volume of the system under consideration, B k is Boltzmann constant, T is temperature, J ? is the ? component of the heat current, L represents an ensemble average. The integration is taken over the time period from zero to infinity. The most fundamental assumption of the FDT is that the mechanism that allows a system at a stable thermodynamic equilibrium resorting from an external disturb is the same as the mechanism that create a steady-state flow at a non-equilibrium condition. The linear response theory assumes that the response of the system is linearly proportional to the external disturb, and then the transport coefficients are proportional to the magnitude of the fluctuations at equilibrium. A statistical approach to define fluctuation is the time-correlation function of two given physical quantities (For example A and B), which can be obtained by taking time dependent quantity A(t) at given time t, and another quantity B(t?) at different time t?, and averaging the product of these two time dependent quantities over some equilibrium ensemble. If A is the same as B, then it 52 is called autocorrelation function. The choice of the time scale origin is arbitrary and the ensemble average is invariant under time displacement, therefore: () (') (0) (' )A tBt A Bt t= ? (3.42) Just like Equ. (3.41) shows, thermal conductivity, the linear response property of the system to the temperature gradient, can be evaluated from the autocorrelation function of the system heat current in thermodynamic equilibrium. In practice, most auto- correlation functions lie between two extremes: strong correlation functioncos( )t? with () sin( )A tt?= and weak correlation function exp( / )t ?? with () 1At = or -1. Now the problem is how to calculate the time dependent heat current for a thermodynamically equilibrium system. Here we follow J.R. Hardy?s notation 40 , under the assumption that the Hamiltonian is in the form: 2 2 i i i i p H V m ?? =+ ?? ?? ? , where i p v is the momentum, i m is the mass, i V is the potential energy associated with the th i particle and the summation is taken over all the particles in the system. The heat current, or energy flux as used in Ref. (40), is then expressed as: 22 11 () , . 22 2 ii i iij j j ii i pp p JVqVHc mm i m ?? ?? ?? ?? =++?+ ???? ?? ? ?? ???? ?? v v vv h (3.43) where ? is the volume of the system, h is Planck constant, i q v is the position of the th i particle. H.c. stands for the Hermitian Conjugate. In the classical limit, heat current can be reduced to: 1 ()() i ii i i j j j j V JvEEqqv q ?? ?? ??? =?+?? ???? ?? ???? ?? v vv v . (3.44) 53 here i v v is the velocity , i E is instantaneous energy and i E is the averaged energy of the th i particle. The evaluation of time-dependent heat current is obtained by performing molecular dynamics (MD) simulations. Empirical Tersoff potential has been used in our study of silicon materials. The Wiener-Khintchine theorem enables the power spectrum of the correlation function to be obtained without having to construct the correlation function itself. If we have the Fourier transform of the time-dependent heat current: () ()exp( )FJtitdt??=? ? (3.45) then the power spectrum of the heat current autocorrelation could be obtained by: 2 () ()ZF? ?= (3.46) Code Current_Correlation.f90 is the implementation for self correlation of the heat current in Equ. (3.41). Multiple runs have been performed to reduce the numerical error. The ensemble averaged heat current autocorrelation function, one the one hand, can estimate the thermal conductivity through Equ. (3.41); on the other hand, its Fourier transformation gives the power spectrum of heat current, which can tell the distribution of the vibrational frequencies (phonon frequencies) that contributes to the reduction of thermal conductivity. Power spectrum analysis is done by code Powerspectrum.f90 54 CHAPTER 4 THERMAL PROPERTIES OF TYPE II CLATHRATE Si 136 4.1 Introduction Modern microelectronics technology is based largely on the semiconducting, diamond-structured silicon (d-Si). Several metastable forms, including amorphous (a- Si) 41 and nano-structured (nano-Si) silicon 42 , also have a variety of industrial applications. High-density polymorphic forms of Si, which are synthesized under high-pressure conditions, have been explored in great detail both theoretically and experimentally. 43 A common feature of the observed high-pressure phases is the increase of the coordination numbers of Si atoms from four (as in the ground state d-Si and meta-stable a-Si or nano- Si) to six or eight. This structural change is accompanied by a change in the electronic properties from semiconducting to metallic. Unfortunately, the technological applications of these octahedrally or eight coordinated Si materials are limited since none of those Si solids is recovered metastably to ambient conditions. Recently, a novel low-density crystalline form of elemental Si was synthesized by removal of Na atoms from the clathrate-structured framework compounds Na x Si 136 44,45,4 prepared by controlled thermal decomposition of NaSi. In contrast to the dense high-pressure phases, the atoms in this 55 new low density Si material are fully tetrahedrally bonded, and it is a wider-gap semiconductor (E g = 1.9 eV) 4 . The new Si allotrope has the type-II clathrate structure, isostructural with low-density inclusion compounds of H 2 O-ice, that has a cubic framework in which each cubic unit-cell contains sixteen 20-atom ?cages? (dodecahedra) and eight 28-atom ?cages? (hexakaidecahedra) (For example, Figure 1 in Ref. 44). In addition to the elemental ?guest-free? form of Si 136 , various ?guest? atoms, including alkali or alkaline earth metals, or halogens, can be incorporated inside the atomic cages to form binary or ternary compounds. Another closely related crystal structure is the type-I clathrate structure. Although various series of guest-encapsulated binary or ternary compounds involving Si framework atoms have been prepared in many laboratories, no one has yet reported a guest-free elemental solid with the type-I clathrate structure, and Si 136 remains as the only low density allotrope of the element that is metastably available at ambient conditions. Both pristine and guest-encapsulated clathrate materials have significant technological potential because they exhibit a very wide spectrum of materials properties. For example, electrical conductivity of Si clathrates ranges from wide-gap semiconducting 4 to metallic 46 and even superconducting (T c =6-8 K) in (Na,Ba)- containing Si clathrates. 47,48 In recent years, the potential applications of Si, Ge, and Sn- based clathrate-structured materials in thermoelectric devices have led to intensive research. 49,50,51,52 In 1995, Slack predicted that open framework structures containing encapsulated rattling guest atoms may exhibit lowered ?glass-like? thermal conductivity due to scattering of acoustic heat-carrying phonons by the guest atoms, while leaving the electrical conductivity via the framework channels largely unaffected. 1 Such materials 56 were described in terms of the ?phonon-glass electron-crystal? paradigm (PGEC). 53,54 This model then raised the possibility of ?tuning? the host-guest chemistry of clathrates independently to result in a desirable combination of electrical and thermal properties. Promising results were demonstrated by Nolas et al, 49,50,52,55 and they have since been extended by other authors. 56,57,58,59 However, the ?rattling? of guest atoms within the cages is not the only mechanism to reduce the thermal conductivity of clathrate materials. Using molecular dynamics simulations, Dong et al showed that the lattice thermal conductivity (?) of guest-free Ge 46 clathrates could be lowered by at least one order of magnitude compared with that of the corresponding diamond-structured Ge crystals. 3 Recently, Nolas et al. 52 demonstrated experimentally that the guest-free clathrate material Si 136 had an extremely low thermal conductivity, as low as that of amorphous SiO 2 , and the smallest value recorded among crystalline solids. This major difference between the two crystalline polymorphs of Si can be caused by either significant increase in lattice anharmonicity in the clathrate-structured material, or the flattening in the phonon dispersion relations associated with formation of the low-density structure, that is correlated with the presence of five-membered rings. Most work on clathrate materials to date has focused on their synthesis and structure characterization, 60 , 61 , 62 and transport measurements, 46, 63 however, studies of their fundamental thermal properties remain limited. It is especially important to explore any special features in the thermal expansion, heat capacity and lattice thermal conductivity of these unusual of expanded-framework semiconducting crystals. As a first step in this area, we used first-principles theoretical methods to predict the measurable thermal properties (such as heat capacity and thermal expansion) of the guest-free 57 clathrate Si 136 . Our results are analyzed and compared with previous data on the well- known ground state diamond-structured phase of this element. 64 Despite noticeable differences in materials density, compressibility, and electronic structures, we find that the two phases have very comparable heat capacities and thermal expansibilities. One important prediction of our calculations is that the clathrate-structured polymorph Si 136 should exhibit a region of negative thermal expansion below 140K, like the diamond- structured phase. This prediction is validated by our experimental measurements. 4.2 Crystal Lattices and Static Equation of State The two Si polymorphs happen to belong to the same face-centered-cubic (fcc) 3Fd m symmetry space group (#227). The Si atoms in the diamond structure occupy a single 8a site with T d symmetry, without any internal degrees of freedom. However, the type-II clathrate lattice contains three distinct sites; 8a (T d ) site, 32e (C 3v ), and 96g (C s ), that have independent sets of internal coordinates (x 32e , x 96g , z 96g ). For static equation of state (EOS) calculations of Si 136 , the internal coordinates were optimized by minimizing the total static energy while maintaining the external lattice vectors at fixed values. The calculated lattice total static energies at various volumes of d-Si and Si 136 crystals are fitted to a 3 rd -order Birch-Murnaghan equation of state (BM-EoS), and the fitting parameters are listed in Table 4.1, along with experimental data. 64,65 The LDA calculations show that the minimal static energy of Si 136 is only 79 meV/atom (i.e. 917 K in the temperature unit) higher than that of the ground state d-Si. Such small energetic differences suggest the clathrate phases are energetically accessible phases. If there exist energetic barriers to prevent back transformation to the ground state 58 diamond phase, the clathrate phases are possible meta-stable phases. This result is consistent with the observed metastability and recovery to ambient conditions of the guest-free clathrate phase Si 136 , synthesized experimentally. 4 The calculated (static) equilibrium cubic lattice constants are 0.5397 nm and 1.455 nm for the 8-atom unit-cell of d-Si and 136-atom unit-cell of Si 136 respectively. These are in excellent agreement with experimental measurements of 0.5431 nm and 1.4644 nm for d-Si 64 and Si 136 65 respectively. Because of its more complex crystal structure, the lattice periodicity in Si 136 clathrate is approximately 2 ? 3 times of that of d-Si. The mass density of the guest-free clathrate is 13% less than that of d-Si. This decrease of density is accompanied by a 17% increase of the compressibility (?=1/K o ) in Si 136 . Table 4.1: LDA calculated static (T=0K) Birch-Murnaghan equation of state of the ground state diamond phase Si (d-Si) and the meta-stable type-II clathrate phase of Si (Si 136 ), and available experimental parameters. The experimental data are for T = 298 K and they are taken from Ref [64] for d-Si and Ref. [65] for Si 136 . d-Si LDA EXPT. Si 136 LDA EXPT. E 0 (eV/atom) -5.954 - -5.875 - V 0 (? 3 /atom) 19.650 20.024 22.650 23.091 B 0 (GPa) 95.54 98 81.877 90 B? 3.936 - 3.937 5.2 59 4.3 Lattice Phonon Spectra Figure 4.1: LDA calculation of (a) phonon dispersion relations, (b) vibrational density of states, and (c) mode Gr?neisen parameters of d-Si at the equilibrium volume. Figure 4.2: LDA calculation of (a) phonon dispersion relations, (b) vibrational density of states, and (c) mode Gr?neisen parameters of Si 136 clathrate at the equilibrium volume. 60 Figure 4.3 (a) 4-fold ?folded? phonon dispersion plots of d-Si, calculated with unit-cell size four times that of the primitive cell ( 4 cell primitive aa= ?=2.16 nm), along with (b) the re-plotted phonon dispersion of Si 136 ( clathrate a = 1.46 nm). Figure 4.1 and 4.2 shows the LDA calculated phonon spectra of the two Si crystals at their respective static equilibrium volumes, at T = 0 K. Since both Si polymorphs are based on tetrahedrally-coordinated Si atoms with similar sp 3 Si-Si bonding, we expect the local force constants and vibrational properties to have analogous behavior. Three regions within the calculated ? i ( q r ) relations in d-Si (Figure 4.2(a) and 4.3(b)) can be identified as due to (a) TA branches (low frequency), (b) low/medium frequency modes due to LA branches, and (c) high frequency optic modes. At first sight, the phonon dispersion relations of Si 136 (Figure 4.2(a)) appear to be considerably more complex than those for d-Si (Figure 4.1(a)). To better illustrate the ?folding? effects due 61 to the enlarged unit cells of clathrate materials, we plotted the ?folded? phonon dispersion relation of d-Si in Figure 4.3(a) using a unit-cell of four times of the primitive unit-cell, along with the re-plotted phonon dispersion relation of Si 136 clathrate in Figure 4.3(b). Clearly, some phonon modes of the clathrate phase can be approximately mapped on to the phonon branches of d-Si, with the dispersion relations ?folded? back towards the first Brillouin zone of d-Si. This observation is important, because it means that the elastic properties and phonon propagation relations within d-Si and Si 136 at small wave-vectors are similar to each other. To further illustrate this similarity, we listed the calculated group velocities of acoustic phonon modes in d-Si and Si 136 in Table 4.2. Table 4.2 Comparison of acoustic velocities (m/s) for d-Si (experimental vs. LDA calculated) and Si 136 . The experimental data labeled with (a) and (b) are taken from Ref. 66 and 67 respectively. Direction Mode d-Si EXPT. d-Si LDA Si 136 LDA TA 5843 a 5109 4800 [001] LA 8433 a 8268 8346 TA1 4673 a 4090 4783 TA2 5830 a 5126 4860 [011] LA 9134 a 8782 8271 TA1 5099 b 4462 4859 TA2 5099 b 4462 4859 [111] LA 9245 b 8982 8283 62 The acoustic phonon velocities in Si 136 differ from their counterparts in d-Si only by a few percent. However, some ?gaps? are developed within the Si 136 ?s ? i (k) relations and they cause important modifications of the Umklapp processes, that results in a dramatic reduction of the thermal conductivity for the Si 136 polymorph. Meanwhile, a prominent ?flattening? can be seen in several of the phonon branches of Si 136 compared with d-Si: especially in the low-frequency regime between 150-200 cm -1 , where many branches are nearly dispersionless (i.e., nearly zero group velocities). The overall features of the vibrational density of states (VDOS) of Si 136 and d-Si appear quite similar, on initial inspection (Figure 4.1(b) and 4.2(b)). However, there are notable differences between the two that can be attributed to specific structural features. First, the Si 136 structure contains 5-membered rings; these cause the expansion of the unit cell motif and the decrease in density, compared with d-Si. This is manifested in the nature and appearance of van Hove singularities within the VDOS. Such singularities are intrinsic to highly symmetric crystals with small unit cells, such as d-Si, and they are readily visible in the g(?) functions. They are absent in the VDOS of Si 136 ; the g(?) plots show instead a ?saw tooth? behavior in the middle of the strong bands in the VDOS (Figure 4.2(b)). Second, due to the flattening in the phonon dispersions discussed in the above text, a sharp peak in the VDOS g(?) of Si 136 emerges at around 175 cm -1 . Third, a blue shift of about 30 cm -1 is observed among the highest frequency transverse optic (TO) bands of Si 136 , with respect to d-Si. Our theoretical results of the phonon spectra in Si clathrates are in agreement with a previous theoretical study using empirical potentials. 68 However, those calculations were based on a Tersoff potential that did not reproduce the 63 30 cm -1 shift in TO frequencies predicted by the first-principles methods. Also, the nearly dispersionless phonon branches occurring near 200 cm -1 predicted by the first-principles calculations were more ?spread out? in frequency, in the empirical study. Figure 4.4 The theoretical mode Gr?neisen parameters ? i along three high-symmetry directions. The ? i of LA (cross) and optic (triangles) phonon modes are all positive; the LA branches have values that lie within a narrow range (~+0.75) and those associated with optic modes range between +0.8 to +1.6. We also find overall similarities but certain essential differences among the Gr?neisen parameters of d-Si vs. Si 136 (i.e. ln ln i i V ? ? ? =? ? ), corresponding to the phonon modes in the two structures (Figure 4.1(c), 4.2(c)). The calculated mode shifts are compared with experimental values in Table 4.3. To compare our calculations with previous predictions of Wei et al. 69 for d-Si, we have plotted i ? along three high- 64 symmetry directions (Figure 4.4). Our calculated Gr?neisen parameters indicate that i ? ranges between -2 and +1.6 cm -1 /GPa, for all phonon branches. The i ? values for TA modes (circles) are negative, with the most negative values located at the Brillouin zone boundary (X or L points). The implications of such negative Gr?neisen parameters for TA phonons on the thermal properties of d-Si and Si 136 at low T are discussed further below. The i ? values of LA (cross) and optic (triangles) phonon modes are all positive; the LA branches have values that lie within a narrow range (~+0.75) and those associated with optic modes range between +0.8 to +1.6. The distinct character of the TA vs. LA branches, and their derived optic phonons, are illustrated in Figure 4.1(c), using q r -points sampled uniformly over a 10 10 10?? grid. The i ? vs. ? i plot for Si 136 clathrate is shown in Figure 4.2(c), using q r -points sampled uniformly over a 5 5 5? ? grid. We can identify similar regions in the corresponding plot for d-Si (Figure 4.1(c)). Although the boundaries between the three regions are less well defined for Si 136 compared with d-Si, 99% of the optic phonon branches for Si 136 can be characterized as derived from TA-like, LA-like, or optic branches issued from d-Si. The values of the Gr?neisen parameters i ? for Si 136 for these groups of modes are similar to those for d-Si; the highest values in both phases are at around +1.6, but the lowest value of i ? in Si 136 is ~ -1.6, as compared with -2 in d-Si. The calculated zone-center frequencies and their pressure derivatives (reported as mode Gr?neisen parameters) are generally in excellent agreement with experimentally measured Raman data (Figure 4.5 and Table 4.3). Selected Raman spectra of Si 136 collected at compression are shown in Figure 4.6 and the corresponding Raman shifts vs. 65 pressure are plotted in Figure 4.7. Experimental results from our collaborators agree with previous data, for the pressure shifts measured for three Raman active modes of Si 136 70 . The Gr?neisen parameters are defined as ln1 () i iT T P ? ? ? ? = ? , where ? T is the compressibility, and its value for Si 136 was taken from Ref. 65. The two lowest frequency modes (T 2g and E g ) are derived directly from the top of the acoustic branches, folded due to the lattice expansion within the Brillouin zone. These appear as distinct peaks in the Raman spectra of Si 136 , and they exhibit a marked frequency decrease with increasing pressure, as predicted by the theoretical calculations. 71 The high-frequency modes that are associated with Si-Si stretching vibrations (>400 cm -1 ) all exhibit positive Gr?neisen shifts also as predicted. Figure 4.5 Raman spectrum obtained for Si 136 at ambient P and T. There is a weak peak at ~520 cm-1 that corresponds mainly to a trace of d-Si impurity in the sample. However, there is also a calculated mode at 516 cm-1 for Si 136 clathrate at this position. 66 Figure 4.6 Selected Raman spectra of Si 136 collected at compression. Asterisk marks diamond-phase silicon. There are some discrepancies between the theoretical and experimentally observed pressure shifts, mainly occurring in the mid-frequency region between 200 cm -1 and 400 cm -1 . The T 2g , E g and A 1g modes at 325, 360 and 387 cm -1 are observed to have negative pressure shifts, although they are all predicted to have positive Gr?neisen parameters. However, that spectral region contains mixed contributions from Si-Si stretching and bending vibrations and it is expected that small deviations from the actual potential energy surface could cause large changes in the coupling between the two types 67 of motion and the resulting pressure shifts of individual modes. However, there is good overall agreement between the observed and calculated mode frequencies at the Brillouin zone center, and with the overall pattern of Gr?neisen shifts as a function of pressure. Pressure, GPa 02468 Waven u mber, cm -1 100 200 300 400 500 T 2g E g T 2g T 2g T 2g E g A 1g T 2g A 1g T 2g T 2g E g Si Figure 4.7 Pressure dependence of Raman shifts of Si 136 and d-Si. 68 Table 4.3 LDA calculated and experimentally measured Raman frequencies (?) and mode Gr?neisen parameters (? i ) for clathrate-structured Si 136 . Mode Frequency cm -1 (Theory) Frequency cm -1 (Experiment) Gr?neisen parameters (? i ) (Theory) Gr?neisen parameters (? i ) (Experiment) T 2g 121 117 -1.17 -1.32?0.04 E g 130 130 -0.71 -0.86?0.04 T 2g 176 184 -1.22 -- T 2g 267 271 0.94 1.01?0.04 A 1g 316 -- 1.04 -- T 2g 325 324 0.93 -1.17?0.04 E g 360 360 1.18 -0.66?0.04 A 1g 397 387 1.08 -0.89?0.05 T 2g 406 401 1.49 0.97?0.04 A 1g 458 454 1.22 1.33?0.04 E g 463 1.26 T 2g 466 466 1.09 -- T 2g 473 472 1.42 -- E g 483 480 1.33 1.21?0.07 T 2g 487 488 1.04 1.39?0.07 (Si) 516 520 0.94 0.95?0.005 2-phonon 920 2-phonon 965 69 4.4 Thermodynamic Potentials and T-P Phase Relations Use Equ. (3.21), we calculated the Helmholtz free energy F over a fine grid of T points, and fitted the V dependency of free energy F at each temperature using the 3 rd - order BM-EOS. Then, we constructed the Gibbs free energy G(T,P) via the Legendre transformations. We derived all other thermodynamic quantities as derivatives of the F(T,V) or G(T,P) functions. We plot the vibrational entropies, i.e. 1 0 () () () ( 1) ln(1 ) BB kT kT PVB B GF Skgdee TT kT ?? ? ?? ? ? ? ?? ?? =? =? = ? ? ??? ?? ?? ? hh h , of the d-Si and Si 136 phases, as functions of T at their respective static equilibrium volumes (Figure 4.8). The metastable Si 136 polymorph possesses slightly larger entropy (S vib ) than d-Si; more exactly by 0.132 k B /atom at 300 K, and 0.156 k B /atom at 600 K. Figure 4.8 The LDA predicted vibrational entropies in d-Si (solid line) and Si 136 (dashed line). 70 This small entropy difference (?S) between the Si 136 and d-Si polymorphs results in a positive and steep dP tr /dT slope for the equilibrium phase boundary between the two phases (Figure 4.9). The phase transition is predicted to lie at negative pressure (-P), i.e., in the tensile stressed regime. 65,72 The Si 136 phase is predicted to become more stable than d-Si at P = ?3.84 GPa, at T=0 K. The transition pressure P tr becomes ?3.11 GPa at 1200 K, and the Clapeyron slope (dP tr /dT) is predicted to be 0.0007 GPa/K, at T=600 K. The steep slope of the transition pressure indicates that the metastable phase boundary between the Si 136 and d-Si polymorphs of elemental Si intercepts the melting curve at negative pressure, so that Si 136 never becomes an equilibrium phase at positive pressure values. 72 Figure 4.9 The theoretically predicted equilibrium T-P phase boundary between the ground state diamond phase and the Si 136 clathrate phase. Note that the transition to the ?expanded? polymorph of the element occurs within the tensile regime, at negative pressure. 71 4.5 Thermal Properties We also calculated the thermal properties for the Si 136 clathrate-structured polymorph at P=0. We first predicted the specific heat capacity () pP S CT T ? = ? (Figure 4.10) under isobaric conditions based on the calculated vibrational entropies as a function of temperature at zero pressure (Figure 4.10). The discrete symbols in Figure 4.10 correspond to experimental data d-Si. 52 Our quasi-harmonic computational results are clearly valid for d-Si; the theoretical prediction lies within a few percent of experimental data. 73,74 The calculated C p of Si 136 is ~0.05 k B /atom greater than that of d-Si at room temperature. The measured C p of Si 136 has been reported by Nolas et al. 52 Although the reported values are noticeably higher than our calculation results, the same group re- measured the C p data recently and the updated experimental data is in much closer agreement with our predicted C p data. 75 Figure 4.10 The calculated and measured specific isobaric heat capacities in d-Si and Si 136 . 72 We then used our LDA results to predict the ?(T) relations at P=0. It is well known that d-Si exhibits a region of negative thermal expansion (NTE) between T =0- 150 K, that is often ascribed to anharmonic phonon-phonon interactions, but that is well reproduced by quasi-harmonic ab initio calculations 69 as well as MD simulations. 76 In our study, the phonon properties as a function of the temperature were evaluated within the LDA, according to the quasi-harmonic model. The predicted ?(T) behavior at low T for d-Si and Si 136 calculated at the same level of theory is shown in Figure 4.11. The result for d-Si indicates a marked NTE region at low temperature, as observed experimentally. Si 136 also reveals itself as a NTE material, with a minimum in the ?(T) function appearing between 0-125K. Low temperature X-ray diffraction (Figure 4.12) measurements from our collaborators are in agreement with this predicted NTE behavior for Si 136 . The variation in the cubic unit cell parameter (a o ) with T (Figure 4.13) was determined from three high angle reflections (733, 660, 751) in Figure 4.12. The precision in the experimental data (shown as error bars in Figure 4.13) is ultimately limited by the intrinsic peak widths of the Si 136 sample; this peak broadening is attributed to defects associated with the Na-loss occurring during formation of the guest-free clathrate. However, the data points collected are sufficient to show at least a flattening in the ?(T) data between 5-200 K, with a minimum apparent near 90-120 K, that may be statistically significant (Figure 4.13). Following this, the data exhibit a positive ?(T) at higher T. The LDA calculated temperature dependence of the coefficient of thermal expansion (?) is in excellent agreement with the experimental result. 73 Figure 4.11 The theoretically predicted linear coefficients of lattice thermal expansion at P=0 in d-Si and Si 136 . Figure 4.12 High-resolution X-ray powder diffraction data for Si 136 illustrated at three temperatures (5K ? bottom, 145K - middle, 275K - top). The three high angle reflections used in the data fitting are marked with an asterisk and represent (in order) the 733, 660 and 751 reflections. 74 Figure 4.13 Low temperature variation (T = 5 - 275 K) variation in unit cell lattice constant (?) for clathrate-structured Si 136 obtained by fitting to high angle powder X-ray diffraction data. The line drawn through the data points is a guide to the eye. The origin of the low-T NTE behavior in d-Si can be associated with the negative Gr?neisen parameter for a TA phonon that reaches the Brillouin zone boundary at ~ 200 cm -1 . The NTE in Si 136 can be traced to a similar cause, as a Raman active mode near this frequency unusually has a large negative Gr?neisen parameter that is both calculated theoretically and observed experimentally (Figure 4.2 and Table 4.3). 4.6. Conclusions We have studied the thermal properties of the novel guest-free clathrate polymorph of silicon (Si 136 ) based on first-principles calculations combined with experimental X-ray and Raman scattering measurements. The Si 136 clathrate is metastable 75 compared with the d-Si phase at ambient P and T. Theory indicates that it becomes thermodynamically stable within a negative pressure regime, at P = -2 to -4 GPa. The dP tr /dT Clapeyron slope is estimated as 7 ? 10 -4 GPa/K from ab initio calculations. Although it has been shown previously that some properties, such as the electronic band gap, are critically dependent upon the lattice expansion between the diamond-structured and ?expanded-volume? clathrate polymorphs, our current studies reveal that the thermal properties of the two phases involving long-wavelength phonons are similar to each other. The vibrational properties of Si 136 phonons are similar to those of d-Si, and they can be understood in terms of Brillouin zone reduction following the unit cell expansion between d-Si and Si 136 . The phonon modes in the two phases also have very similar characteristics. We find that the coefficients of thermal expansion in the two Si phases are comparable in our studies, which suggests that it is less likely that the significant reduction of lattice thermal conductivity in clathrate materials is mainly caused by any large increases of anharmonic lattice interactions in clathrate systems. 76 CHAPTER 5 THERMAL CONDUCTIVITY OF TETRHEDRALLY BONDED SILICON CRYSTALS ? THE GREEN-KUBO APPROACH 5.1 Introduction Thermoelectricity is a class of fundamental physical phenomena that involves direct conversion between thermal and electrical and thermal energies. Thermoelectric devices have been developed for applications such as power generation with waste heat in automobile engines, or solid-state cooling for microelectronic devices. The key limiting factor of current thermoelectric technologies is the lack of high-performance thermoelectric materials, whose efficiency is described with a dimension-less figure of merit 2 /ZT S T? ?= , where T, S, ? and ? are temperature, Seebeck coefficient, electrical and thermal conductivity respectively. For many practical applications, the ZT of a thermoelectric material should at least be 3 or larger. As electrical and thermal transport properties of solids are correlated through phonon scattering at microscopic level, large /? ? ratio is rare in conventional materials systems. An ideal high ZT thermoelectric material must simultaneously have a low thermal conductivity? and a high Seebeck coefficient S and electrical conductivity ? . Usually, electrical transport 77 has a strong dependence on carrier concentration. For example, the electrical conductivity in semiconductors can be altered by several orders of magnitudes by doping. In contrast, the lattice thermal conductivity of a material is difficult to control at a given (T,P) condition. In 1995, Slack predicted that crystals with open framework structures and encapsulated rattling guest atoms may exhibit ?glass-like? low thermal conductivity due to the scattering of acoustic heat-carrying phonons by the guest atoms, while leaving the high electrical conductivity via the framework channels largely unaffected. 1 Such materials were described in terms of the ?phonon-glass-electron-crystal? paradigm (PGEC). 53,54 This model then raised the possibility of ?tuning? the host-guest chemistry in Si, Ge, or Sn based clathrate materials to achieve a desirable combination of electrical and thermal properties. Later, promising results have been demonstrated by Nolas et al, 49,77 ,78,79 with type-I Ge clathrates with Sr atoms at the guest sites and Ga and Ge atoms at framework sites (Sr 8 Ga 8 Ge 38 ), and their studies have since been extended by other authors. 80,81,82,83 The measured lattice thermal conductivity of Sr 8 Ga 8 Ge 38 is two-order of magnitude lower than that of d-Ge. In addition to the importance to thermoelectric applications, such glass-like low thermal conductivity in a crystalline system is of great interests from materials theory point of view. Several years ago, Dong et al. predicted that ?rattling? of guest atoms within the cages are not the only contributors to this significant reduction in thermal conductivity of the measured clathrate materials. 3 Using statistical Green-Kubo theory and molecular dynamics simulations, they demonstrated that the lattice thermal conductivity (? ) of guest-free Ge 46 clathrates could be lowered by about one order of 78 magnitude compared with that of the corresponding diamond-structured Ge crystals. 3 They reported an unexplained oscillation feature in the time correlation function of heat current in Ge 46 clathrate, and this feature is of close relevance to the factor of ten reduction in ? . The encapsulated guest atoms provide an additional one order of magnitude reduction in lattice conductivity of Ge clathrates. This prediction of significant reduction in guest-free Ge clathrates remains un-verified since no laboratories have successfully removed encapsulated guest atoms from the Ge-based clathrates. The only pristine clathrate system that has been synthesized is the type-II silicon clathrate (Si 136 ). Si 136 , which is very similar to d-Si in the aspects of local Si-Si bonding environment, thermal properties, and mechanic properties, is an ideal material system for studying the ?structural effects? of lattice thermal conductivity. Recently, Nolas et al. 52 demonstrated experimentally that the guest-free Si 136 had an extremely low thermal conductivity, as low as that of amorphous SiO 2 . With the two order magnitude reduction of ? from the d- Si, the ? of Si 136 crystals is the smallest value recorded among crystalline solids. This experimental measurement is in consistent with the prediction of Dong et al. on the significant reduction of ? for guest-free clathrate systems. Yet, the measured reduction is one order of magnitude larger. The current study aims to (1) theoretically predict the ratio of lattice thermal conductivity of d-Si and Si 136 , and (2) quantitatively estimate the effect of each factor that contributes to the reduction of lattice thermal conductivity in Si 136 . Although the previous study has demonstrated a factor of ten in ? reduction in clathrate, the calculations were based on the Ge-based type-I clathrate. The current study directly examines the type-II Si clathrates. The simulation techniques adopted here are similar to those of the previous 79 study on the Ge 46 clathrate. The molecular dynamics (MD) simulations were performed based on an empirical potential due to the requirement of large numbers of MD steps ( 6 1.5 10? steps) and large super cell models. The accuracy of an empirical interatomic potential depends on the set of experimental data adopted in fitting. Typically, the potentials are fitted with measured structural, dynamic, and elastic properties. On the other hand, the lattice thermal conductivity is strongly influenced by lattice anharmonicity, which was not included in the fitting processes of the adopted empirical potentials. Nevertheless, the ratio between the calculated lattice conductivity of the two Si phases that contain similar the local bonding environment is expected to be insensitive to the exact value of lattice anharmonicity. Compared to the previous studies of Ge systems, we adopted a much larger super cell models (2744 atoms for d-Si and 3672 atoms for Si 136 ) to further reduce correlation artifacts due to the periodic boundary conditions. Dependence of thermal conductivity on the super cell size for bulk silicon has been studied by Volz and Chen 84 and their results are served as the guide for our choice of super cell. 5.2 Empirical Si-Si Potential The empirical potential adopted in this study is developed by Tersoff 85 . The interatomic interaction among Si atoms in a tetrahedral bonding environment is modeled as in Equ. (5.1) and fitting parameters for silicon are listed in Table 5.1. 80 ( )[ ( ) ( )]; () exp( ); () exp( ); 1, 11 () cos ( )/( ), . 22 0, C ij R ij ij A ij Rij ij ijij Aij ij ijij ij ij Cij ijij ijij ijijij ij ij Vij f r f r b f r fr A r fr B r rR f rrRSRrS rR ? ? ? =+ =? =? ? < ? ? ? ??=+ ? ? << ? ?? ? > ? ? (5.1) where 1/2 , 22 2 2 2 1/2 1/2 1/2 1/2 (1 ) ; () ( ); ()1 / / ( cos); ()/2; ()/2; (); (); (); (). ii i nn n ij ij i ij ij C ij ik ijk kij ijk i i i i i ijk ij i j ij i j ij i j ij i j ij i j ij i j b fr g gcdcdh AAA BBB RRR SSS ??? ??? ?? ??? ??? ? ? =+ = ? ?=+ ? + ? ? ? =+ =+ = = = = ? Table 5.1 Fitting parameters for silicon to be used in Equ. (5.1). A(eV) 3 1.8308 10? B(eV) 2 4.7118 10? 1 ()? ? ? o 2.4799 1 ()? ? ? o 1.7322 ? 6 1.1000 10 ? ? n 1 7.8734 10 ? ? c 5 1.0039 10? d 1 1.6217 10? h 1 5.9825 10 ? ?? ()R ? o 2.7 ()S ? o 3.0 81 Table 5.2 listed the calculated 3 rd order Birch-Murnaghan equations of state of the two Si phases, the diamond-structure Si (d-Si) and type-II clathrate (Si 136 ), using the empirical Tersoff potential described above. As a comparison, the results of our first- principles LDA calculations and experiment (more details can be found in Chapter 4) are also listed. The Tersoff potential predicted energetic difference between the Si 136 and d-Si phases is in excellent agreement with that predicted by LDA calculations. Note that the Tersoff potentials are fitted with the measured structural and elastic data of d-Si, and consequently its calculated V 0 and B 0 are nearly identical to those of room temperature experimental data of d-Si. The more importance issue that the Tersoff potential correctly predicted the differences in structural and elastic properties of the two tetrahedrally bonded Si phases. As shown in Table 5.2, the Tersoff potential predicted increase in volume and reduction in bulk modulus are about 14.72% and 13.18% respectively. The same numbers predicted by the LDA calculations are 15.27% and 14.30% respectively, while the experimental data at room temperature are15.32% and 8.16% respectively. 82 Table 5.2 Comparison between VASP results and Tersoff MD results: Equilibrium Birch- Murnaghan equation of state fitting parameters for d-Si and Si 136 together with experimental data. d-Si Si 136 Tersoff LDA Expt. Tersoff LDA Expt. E 0 (eV/atom) -4.630 -5.954 - -4.591 -5.875 - V 0 (? 3 /atom) 20.035 19.650 20.024 22.984 22.650 23.091 B 0 (GPa) 97.560 95.54 98 84.701 81.877 90 B? 4.294 3.936 - 4.280 3.937 5.2 Next, we compared the calculated phonon spectra by two theoretical methods. It has been shown that, in Si 136 , 102 eigenmodes have only 42 distinct frequencies. The irreducible representation of Si 136 given by: 12 1 212 1 2 31 45 8 3 48 5 Raman IR g ggg gu uuu u A AET TAAET T?= + + + + + + + + + where the superscript Raman and IR indicate that the mode is Raman active and IR active respectively. Phonon mode representation at Gamma point can be identified by mode symmetry analysis (Appendix I). Table 5.3 listed the results of calculated Gamma point phonon frequencies of Si 136 by the Tersoff potential and LDA methods. 83 Table 5.3 Si 136 ? point phonon frequencies from both Tersoff MD and VASP calculation. Mode representation Frequency (cm -1 ) (Tersoff) Frequency (cm -1 ) LDA IR active Raman active T2g 135 117 T1g 141 127 yes Eg 148 128 T2u 165 143 T1g 175 147 yes Eu 189 157 T1u 206 165 yes T2u 209 169 T1g 241 169 yes A2g 243 173 T1u 248 174 yes T2g 249 186 T1u 278 263 yes T2g 278 269 A2u 289 278 Eu 313 289 A1g 334 314 T2g 341 322 Eg 383 359 T2u 394 365 T1u 410 376 yes A2u 426 391 A1g 434 393 T2g 434 404 A2u 449 415 T1g 426 420 yes T1u 463 424 yes A1g 483 456 Eu 485 446 T1u 503 454 yes T2u 504 459 Eg 504 459 T2g 506 465 Eu 511 460 T2g 518 471 A1u 523 477 T1g 524 474 yes T1u 525 474 yes Eg 528 481 T2g 529 482 T2u 532 484 84 The full phonon dispersions of d-Si and Si 136 from this empirical potential calculation are shown in Figure 5.1 and 5.2 respectively. Figure 5.1 Phonon dispersion of d-Si using Tersoff empirical potential. Figure 5.2 Phonon dispersion of Si 136 using Tersoff empirical potential. 85 Over all, the results of Si 136 and d-Si calculated with the Si Tersoff potential are in good agreement with those obtained with the first-principles DFT method. 5.3 Results and Discussions In the current study, a 2744-atom super cell, i.e. 7 7 7? ? of the 8-atom cubic unit cell, and a 3672-atom super cell (a 333? ? of the 136-atom cubic cell) have been chosen d-Si and type-II Si 136 clathrate respectively. In the MD simulation, we chose the single simulation time-step as 1 fs, and the simulations have been run for 2 21 and 2 20 steps for d- Si and Si 136 respectively. The powers of 2 were chosen for the time steps because it is computationally more efficient during the Fourier transformation (FT) in the correlation function calculation. Before initiating each micro-canonical (NVE) MD simulation, each system was equilibrated for a period of 100 ps with a constant temperature Gaussian thermostat 86 , which rescales the velocities of atoms at each time step to achieve a constant system temperature; and the results have been averaged with multiple starting configurations to minimize the statistical error. Thermal conductivity of d-Si and Si 136 at room temperature was calculated by using Equ. (3.41). For a cubic system, ?? ? ?= if ? ?= and zero otherwise, where ? is: 2 00 () () (0) 3 B Gt J tJ dt kT ?? ? ? ?? ? == ? ?? . (5.2) The direct calculation of the time-averaged correlation function G(t) is time-consuming. Instead, we first do Fourier-transformation for the heat current function J(t), then 86 calculate the power-spectrum of the () (0)JtJ ?? based on the convolution theorem 87 , and finally carry out an inverse Fourier transformation of the power spectrum. The better illustrate the time-decay in the correlation functions, we define a normalized correlation function ()g t as: () () (0) () (0) (0) (0) Gt JtJ gt GJJ < > == < > (5.3) For a simple system whose heat current autocorrelation function decays simply as an exponential function of time, i.e. () t g te ? ? = , the calculated ? is (0)G ? . Compared with the simple kinetic model given by Equ. (3.32), the (0)G and ? can be interpreted as the averaged 2 1 3 Vg CV and phonon life time respectively. Figure 5.3 Macroscopic heat current fluctuations with time of d-Si system 87 As described in previous section, the total macroscopic current J(t) fluctuates around zero at equilibrium. Our MD calculated fluctuation of macroscopic heat current in d-Si is illustrated in Figure 5.3. The magnitude of this fluctuation is proportional to G(0), as defined Equ. (5.3). At T=300K, the calculated G(0) for d-Si is 34.45 W/K/M/ps. Figure 5.4 shows our calculated normalized heat current autocorrelation in d-Si as a function of time for a time period of 400ps. The ()g t starts at a value of unity, and falls rapidly during the first 0.03ps to about 0.32. Then the value starts to decay at a relatively slow rate with a small oscillation added to the descending function of time, and goes to zero beyond 400ps. The small oscillation is not due to statistics as the results are obtained with ensemble averaged over at least 8 different initial conditions. Figure 5.4 The normalized time correlation function ()g t at room temperature for d-Si. (a) an overall look; (b) close-up look at the beginning of the MD run; (c) close-up look at decaying of ()g t after a long time MD run. 88 It is also worth to point out that the g(t) correlation function does NOT decay exponentially. Therefore the averaged phonon life time can not be uniquely defined. Fitting the data with the simple exponential function for some short time periods, we get the fitted ?? ? ranging from above 200ps around t=5ps to less than 50ps for t around 100ps. Nevertheless, we define int ? as the integral of the normalized correlation function g(t) from 0t = to t ??, and our int ? at T=300K is estimated to be 14ps for d-Si. Based on this definition, we obtained the calculated lattice thermal conductivity of d-Si at 300K is around 480 W/K/M. This value is about twice larger than the experimental measured value (283 W/K/M) of isotope enriched pure d-Si crystal 88 . The factor of two overestimation of lattice thermal conductivity based on the empirical Tersoff potentials is considered reasonable as the intrinsic lattice anharmonicity is not tested in the fitting process of Tersoff potentials. The under-estimation of lattice anharmonicity by the adopted Si Tersoff potential is likely to be the same in both d-Si and Si 136 . In the rest of the study, we limit our discussions primarily on the ratio between the two calculated lattice thermal conductivities, instead of on their absolute values. 89 Figure 5.5 Macroscopic heat current fluctuations with time of Si 136 system Figure 5.6 The normalized time-dependent correlation function ()g t at room temperature for Si 136 . (a) an overall look; (b) close-up look at the beginning of the MD run; (c) close- up look at decaying of ()g t after a long time MD run. 90 A small segment of J(t) function of Si 136 is plotted in Figure 5.5. The overall feature is similar to that of d-Si. More interestingly, the G(0) of Si 136 is estimated to be 36.66 W/K/M/ps, which is very close to the value in d-Si. On the other hand, the correlation function of the heat current in Si 136 is dramatically different from that in d-Si. As shown in Figure 5.6, the g(t) oscillates between positive and negative one, while the envelope of the oscillation decays in a fashion similar to that in g(t) of d-Si. This oscillation feature was first found in the g(t) of Ge 46 by Dong et al. 3 . Due to the fact that the thermal conductivity is proportional to the integration of g(t), i.e. the area beneath the ()g t curves, it is obvious that the oscillation appearing in Si 136 system will cause the ?cancellation effect?. The effect of this oscillation can be illustrated with the following model: () cos( ) t t j j j j g tAe Be t ? ? ? ? ? =+ ? (5.4) Obviously the time-integration of all the B j related terms gives zero, and the integration of g(t) over time produces: int A? ?= . (5.5) Effectively, there is factor of ?A? reduction in the calculated lattice conductivity. For the results of Si 136 at 300K, the calculated int ? is as small as 1.4ps, about 10% of that calculated for d-Si. As the G(0) are comparable in d-Si and Si 136 , our Green-Kubo calculations predicted that the lattice thermal conductivity of Si 136 (? at 300K ~ 51 W mK ) is only about 10.6% of that d-Si (482 W mK ). Our predicted one-order of magnitude reduction in type-II Si clathrates is comparable to the previous findings in 91 type-I Ge clathrates, while the experiment measurement of Nolas et al. 5 reported a much larger two-order of magnitude reduction in Si 136 . To reveal the origins of the oscillations terms in the g(t) of Si 136 , we performed the power spectrum analysis for both d-Si and Si 136 , and the results are plotted in Figure 5.7 and 5.8 respectively. To compare with the phonon spectra, the frequency is converted in the unit of cm -1 . As expected, the high-frequency components in d-Si are negligibly small, while there are several strong high-frequency peaks in the power spectrum of Si 136 . The intensity of the peaks represents the relative magnitude of j B terms in the model described by Equ. (5.4). A closer examination further reveals that the strong peaks in the power spectrum of Si 136 coincide only with the vibration frequencies of the ?-point T 1u phonon modes, whose calculated frequencies are 206, 248, 278, 410, 463, 503 and 526 cm -1 . The four strongest high-frequency peaks in the power spectrum appear at the frequencies of 206, 248, 410 and 463 cm -1 . There are also two weak peaks around 278 and 503 cm -1 . 0 100 200 300 400 500 600 700 800 900 1000 0 50 100 150 200 250 300 I n tens i t y (ar b . unit) Frequency (cm -1 ) dSi Figure 5.7 Power Spectrum of heat current autocorrelation function for d-Si. 92 0 100 200 300 400 500 600 700 800 900 1000 0 20 40 60 80 100 I n t ensit y (arb. unit) Frequency (cm -1 ) Si 136 Figure 5.8 Power Spectrum of heat current autocorrelation function for Si 136 The observed T 1u mode associated oscillation in current correlation function g(t) is consistent with group theory. If a crystal is thermally excited in a single vibration eigen mode, the kinetic and potential energy will be transferred among atoms periodically. Equivalently, this can be viewed as having non-zero ?heat current? at a given site at an instance. Yet, because the current that is associated with the vibration is also periodic, there is no net thermal energy flux at any atom site over a vibration period. Note that this periodic site heat current is distinctly different from those heat currents caused by a statistical random process. Although a site heat current can be non-zero for many vibration modes, the only non-vanishing TOTAL heat currents are those associated with the T 1u vibration modes. This is a general result for any crystals. We now predict that the oscillation feature in the correlation functions of heat current exists in all the crystals that contain T 1u vibration modes. 93 The int 136 int () () Si dSi ? ? ? ratio is around 1/10, and according to Equ. (5.5), the reduction might related to either the reduction of A, i.e. the increase of component of oscillation terms in g(t), or the reduction of? , i.e. the increase of decay rate due to phonon-phonon scattering. To quantitatively estimate the contribution from each term, we performed the low-pass filtering on the g(t) plot to eliminate the high-frequency components, which include and are not limited to the major oscillation terms. The results are plotted in Figure 5.9 and Figure 5.10 for d-Si and Si 136 respectively. The filtered g(t) functions are plotted as the red solid lines. The effective (0) eff G are estimated 0.256 (0) 8.82 dSi W G mKps ? ?= and 136 0.045 (0) 1.65 Si W G mKps ?= for d- Si and Si 136 respectively. As a result, we estimate a factor of 0.187 reduction in the effective G(0) of Si 136 , comparing with d-Si which contains similar tetrahedral Si-Si bonds, which is associated with the oscillation feature in Si clathrates. 94 Figure 5.9 Normalized the heat current autocorrelation of d-Si before (black line) and after (red line) the low-pass filter with filtering frequency 0.5Hz. Figure 5.10 Normalized the heat current autocorrelation of Si 136 before (black line) and after (red line) the low-pass filter with filtering frequency 0.5Hz. 95 The log(g(t)) plots are shown in Figure 5.11. The overall decay rates in the two plots are comparable. Again, neither plot is a simple linear function; hence it is difficult to uniquely define an effective eff ? simply from these plots. Based on the model of int (0) (0) eff eff GG?? ?== , we estimate the ratio of eff ? in Si 136 and d-Si is around 0.57, which suggests about 43% reduction in the effective phonon life time in Si 136 . Figure 5.11 log plot of the normalized low-pass (0.5Hz) filtered heat current autocorrelation of d-Si (black line) and Si 136 (red line). 5.4 Conclusion We have carried out a series of large super-cell MD simulations of d-Si and Si 136 to evaluate the autocorrelation function of heat current when the systems are at thermal equilibrium. The lattice thermal conductivity is calculated based on Green-Kubo formula. We find that ? in Si 136 is only 10.4% of that in dSi, and the reduction is contributed by 96 an 81.3% reduction of (0) eff G and a 43% reduction of eff ? in Si 136 compared to d-Si. Our predicted 89.4% reduction of ? in type-II Si clathrate is close to the predicted reduction in type-I Ge clathrate. Meanwhile, Nolas reported a 96.7% reduction in Si 136 in 2001, which is significantly larger than the theoretical results presented here. Further experimental and theoretical studies are needed to resolve this discrepancy. 97 CHAPTER 6 LATTICE ANHARMONICITY OF TETRAHEDRALLY BONDED SILICON CRYSTALS 6.1 Introduction In chapter 5, we have calculated the thermal conductivity for d-Si and Si 136 using Green-Kubo formula with empirical Tersoff potential. The calculation revealed that the thermal conductivity of Si 136 is only about 10.6% of that of d-Si at 300K. According to kinetic theory, this reduction could be caused by noticeable decrease of heat capacity, or significant increase in lattice anharmonicity in the clathrate-structured material, or the flattening in the phonon dispersion relations associated with formation of the low-density structure, which is correlated with the presence of five-membered rings. However Green- Kubo formulism does not provide such detail information but a final thermal conductivity. It is the purpose of this chapter to investigate the individual contribution to the reduction of lattice thermal conductivity in Si 136 . To yield an apple to apple comparison, Tersoffpotential was again used in this study. Both heat capacity and group velocity can be calculated with harmonic phonon spectra, while phonon life time has to take consideration of lattice anharmonicity. The 98 algorithm to calculate third order anharmonicity has been presented in Chapter 3, where symmetry of crystal system has been greatly taken advantage of. 128-atom d-Si and 136- atom Si 136 super cells were used in the calculation for both harmonicity and anharmonicity. Harmonic phonon calculation is the same as what we do for first principles method in chapter 4 and will not be repeated here. The dimension of anharmonicity tensor A is 384 384 384? ? and 408 408 408? ? for d-Si and Si 136 respectively. Among all the elements of tensor A , most of them are near zero and thus negligible compared to major terms, since the effect of displacement of atom i to atom j is almost none if, as far as atom i is concerned, atom j is farther than its second or third nearest neighbors. Thus a cutoff distance can be chosen for the truncation of these terms. Our calculation for both d-Si and Si 136 showed the distance between an atom and its second nearest neighbor is large enough to serve as the cutoff. 6.2 Single Effective Phonon Life Time Approximation Exact evaluation of the kinetic transport model in Equ. (3.32) requires the information of life-time of each phonon mode. As the first-order approximation, where we assume all phonon modes have the same effective life time eff ? , the kinetic transport equation can be simplified as: 2 , 11 [(,)(,)] 3 Vg ef qvq CqvV qv N ??= ? vv vv (6.1) Under this single effective phonon life time approximation, we are able to calculate term 2 , 11 (,) (,) 3 Vg qvq CqvV qv N ? v v vv using the calculated phonon spectra, and estimate the value of 99 eff ? based on the experimentally measured heat conductivity for d-Si. We have found the temperature dependence of eff ? (blue-plus-dashed line in Figure 6.1) by performing a calculation for mode heat capacity and group velocity at various temperature points for d- Si. Assuming Si 136 and d-Si have the same effective phonon life time: 136Si dSi eff eff ??= , we have calculated heat conductivity for Si 136 which is shown as the red-diamond-dashed square in Figure 6.1. Mode heat capacity and group velocity for both systems were calculated with k-point grids 16 16 16? ? and 888? ? for d-Si and Si 136 respectively. The experimental data are based on natural Si crystals, which include isotope disorder. At T=300K, the isotope-enriched d-Si is about 60% higher in ? , comparing with natural d- Si. 0 200 400 600 800 1000 1200 1400 1600 1800 0 50 100 150 200 250 300 0 20 40 60 80 100 Kappa (W /K/m) Temperature (K) dSi Expt. Kappa Data Si 136 estimated Kappa Data Phonon Life Time ? Tao (ps) Figure 6.1 Fitted effective phonon life time eff ? of d-Si to its experimental heat conductivity. Temperature dependence of Si 136 was estimated assuming the equal eff ? . 100 Again, we focus on the ratio between the ? of Si 136 and d-Si. Figure 6.1 showed that thermal conductivity of Si 136 is much lower than that of d-Si assuming that they have the same effective phonon life time. At 300K, the effective eff ? is about 25.92 ps. ? of d- Si is 156 W/Km and ? of Si 136 is 31.8 W/Km, yielding a conductivity ratio 136 20% Si dSi ? ? = . This ratio is in good agreement of the estimated ratio in G(0) eff in our Green-Kubo study. 6.3 Lattice Anharmonicity To estimate the ratio of eff ? , we directly calculated the 3 rd order lattice anharmonicity tensor A using our newly developed finite difference algorithm, in which atoms are only allowed to move along x-direction, y-direction, or z-direction by ? of the same magnitude, the absolute value of each individual term is thus dependent on the setting of the Cartesian coordinates. Table 6.1 and Table 6.2 list the coordinates of some representative atoms in the super cell of d-Si and Si 136 respectively. The coordinates of these atoms are relative to the Cartesian axes showed in Figure 6.1 and Figure 6.2 for d- Si and Si 136 respectively, which provide a direct visualization of their atomic model. In d- Si, all the atoms are positioned at the same Wyckoff site 8a, each silicon atom is neighbored by 4 equivalent silicon atoms. However, in Si 136 , the neighbors are not all equivalent, each atom at 8a site is neighbored to 4 atoms which are all at 32e site; each atom at 32e site is neighbored to 1 atom at 8a site and 3 atoms at 96g site; and each atom at 96g is neighbored to 1 atom at 32e site and 3 atoms at 96g site. 101 Table 6.1 List of coordinates of all 4 first nearest neighbors of two representative atoms at Wyckoff site 8a in our chosen super-cell of d-Si. Site (example atom index) Neighbor Site (example atom index) Relative coordinates of the neighbors in our chosen cell Bond Length 8a (#65) (1.345,1.345,1.345) 2.352 8a (#68) (1.345,-1.345,-1.345) 2.352 8a (#77) (-1.345,1.345,-1.345) 2.352 8a (#1) 8a (#113) (-1.345,-1.345,1.345) 2.352 8a (#1) (-1.345,-1.345,-1.345) 2.352 8a (#2) ( 1.345, 1.345,-1.345) 2.352 8a (#5) ( 1.345,-1.345, 1.345) 2.352 8a (#65) 8a(#17) (-1.345, 1.345, 1.345) 2.352 Figure 6.2 Atomic model of d-Si. 102 Table 6.2 List of coordinates of all 4 first nearest neighbors of one representative atom for each Wyckoff site in our chosen super-cell of Si 136 . Site (example atom index) Neighbor Site (example atom index) Relative coordinates of the neighbors in our chosen cell Bond Length other 32e (#9) (1.345,1.345,1.345) 2.330 32e (#16) (1.345,-1.345,-1.345) 2.330 32e (#19) (-1.345,1.345,-1.345) 2.330 8a (#1) 32e (#22) (-1.345,-1.345,1.345) 2.330 8a (#1) (-1.345,-1.345,-1.345) 2.330 96g (#42) (-0.503, -0.503, 2.244) 2.354 96g (#60) (2.244, -0.503, -0.503) 2.354 32e (#9) 96g (#75) (-0.503, 2.244, -0.503) 2.354 32e (#9) (0.503, 0.503, -2.244) 2.354 96g (#54) (-1.684,-1.684, 0.000) 2.382 Type (a) 96g (#114) (-0.908, 1.971, 0.908) 2.352 Type (b) 96g (#42) 96g (#126) (1.971,-0.908, 0.908) 2.352 Type (b) Figure 6.3 Atomic model of Si 136 . Three different types of atoms have been indicated by different color: 8a (red), 32e (blue) and 96g (gold). 103 Under current Cartesian setting, there are in total 28 nonzero independent ijk A terms for d-Si. Si 136 structure is more complex, hence the number of independent terms is much more, our calculation yields 426 nonzero independent terms. Our data suggests that those ijk A which associate with bond stretching are the major terms, which we think is very reasonable in physics. Interatomic potential is the most unsymmetrical along that direction, on the one side, atoms is strongly ?bonded? with another atom; on the other side, there is no atom at all. Since every atom in both d-Si and Si 136 is tetrahedrally bonded with the other four silicon atoms, comparison of ijk A along bond direction is desired. However current Cartesian setting does not provide a direct evaluation of these terms, but they can be derived using a coordinate transformation. Figure 6.4 shows such a transformation schematically. Given a bond linked by 2 atoms (#1 and #2), In the original Cartesian setting (x,y,z), there are in total 216 ijk A terms associated with these two atoms, let?s call them (, , ) old A ijk , a 666?? tensor where , , 1,2, 3ijk or= represents the displacements of atom #1 along original x, y, and z direction respectively, and ,, 4,5, 6i j kor= the displacements of atom #2. What we are looking for are those ijk A in the new coordinates setting (x?,y?,z?) with ,, 3, ,6i j kor= , meaning the displacements from these two atoms along bond direction z?, which corresponds to the bond stretching. The transformation is such that: (, , ) ijk il jm kl old lmn A TT T A lmn= ? (6.2) 104 and cos cos sin cos sin 0 0 0 sin cos 0 0 0 0 cos sin sin sin cos 0 0 0 000coscossincossin inco0 0 0 0 cos sin sin sin cos ?? ?? ? ?? ?? ?? ? ? ??? ? ?? ? ??? ? ? ? = ? ? T where old ? and ? are both 6 6 6? ? tensor involving all the terms associated with these two atoms before and after the transformation. T is the transformation matrix, ? is the angle between the bond direction (z?) and z axis of original Cartesian setting, and ? is the angel between the projected bond in x-y plane and x axis. After rotation, any bond between two atoms that are tetrahedrally bonded can be rotated to z? direction like what we showed here in Figure 6.4, where the new z axis (z?) is along bond direction. Both new x (x?) and new y (y?) axes are not specified because choice of them does not affect the bond-stretching terms as long as they both are perpendicular to z?. Figure 6.4 A schematic show of coordinate transformation to yield a bond-stretching A ijk . ? ? y z x z? #1 #2 105 Figure 6.5 Any bond between two atoms that are tetrahedrally bonded can be rotated to z? direction like what we showed here. For each bond there are four possible combinations of ,,i j k that corresponds to a bond stretching and they are 333 336 366 ,,A AA and 666 A , while the order of ( , , )i j k does not change the value. The calculation shows that if a bond is connected by two atoms at the same Wyckoff site, then 333 666 A A=? and 336 366 A A=? . The coordinates and bonding of both systems are listed in Table 6.1 and Table 6.2 showing that there is only one type of bond in d-Si (8a-8a), while there are four types of bonds in Si 136 , and they are 8a-32e, 32e-96g, 96g-96g(a), 96g-96g(b). Pressure dependence of bond-stretching ijk A for all these 5 types of bond has been investigated and they are shown in Figure 6.6 for d-Si and Figure 6.7-6.10 for Si 136 respectively. 106 Figure 6.6 Pressure dependence of bond stretching A ijk terms in d-Si. Figure 6.7 Pressure dependence of 8a-32e bond stretching A ijk terms in Si 136 . 1 65 -y 107 Figure 6.8 Pressure dependence of 32e-96g bond stretching A ijk terms in Si 136 . Figure 6.9 Pressure dependence of 96g-96g type (a) bond stretching A ijk terms in Si 136 . 108 Figure 6.10 Pressure dependence of 96g-96g type (b) bond stretching A ijk terms in Si 136 . As one can see from Figure 6.7-6.10, the bond stretching anharmonicity varies for different types of bonds, and they all increase with pressure in a slightly different rate. d- Si has only one type of bond, while Si 136 has four different types of bond. To compare the anharmonicity contribution to lattice thermal conductivity of d-Si with that of Si 136 , we first calculated the averaged bond-stretching anharmonicity for Si 136 . According to the bonding nature of Si 136 , we define the anharmonicity associated with each site as: 832 (8 ) ae ijk ijk AaA ? = (6.3a) 32 8 2 32 96 2 ()3( ) (32 ) 4 ea e g ijk ijk ijk AA Ae ?? +? = (6.3b) 96 32 2 (96 96 ) 2 (96 96 ) 2 ()( ( (96 ) 4 ge gga ggb ijk ijk ijk ijk AA A Ag ?? ? ++? = (6.3c) 109 The overall averaged ijk A for Si 136 is then given by: 8 (8 ) 32 (32 ) 96 (96) 136 ijk ijk ijk ijk Aa A e A A ?+? +? = (6.4) The averaged site anharmonicity and overall bond anharmonicity at zero pressure is listed in Table (6.3) Table 6.3 Averaged site anharmonicity and overall bond anharmonicity of Si 136 . A ijk (8 )A a (32 )A e (96 )Ag 333 A A 333 42.731 40.009 37.682 38.527 A 336 42.360 40.037 38.505 39.092 The absolute bond anharmonicity represents the contribution of anharmonic energy to the total energy. In order to compare anharmonicity contribution of Si 136 and d- Si to the lattice thermal conductivity, we propose a dimensionless quantity: relative anharmonicity ? : ijk A B ? = (6.5) where B is the bulk modulus. ijk A is in the unit of eV/? 3 ,which is equivalent to 160 GPa, and B is in the unit of GPa. 2 ? is proportional to the phonon scattering rate, and thus inversely proportional to the phonon life time. Table (6.4) listed ? for both d-Si and Si 136 110 and it is shown that 136 118.3% Si dSi ? ? = if ? is calculated from A 333 , leading to ~30% reduction of effective phonon life time in Si 136 ( 136 70% Si dSi ? ? ? ). Table 6.4 Comparison of relative anharmoncity between d-Si and Si 136 System B(GPa) 333 A ( eV/? 3 ) 333 /AB 336 A ( eV/? 3 ) 336 /AB d-Si 97.560 38.527 63.19 39.522 64.82 Si 136 84.701 39.561 74.73 39.092 73.85 As mentioned in chapter 3, the reliability of 3 rd order force constants are tested by calculating Gr?neisen parameters using Equ. (3.10), and then comparing them with those calculated using finite difference method Equ. (3.9), our calculated Gr?neisen parameters from both methods for d-Si and Si 136 are presented in Figure 6.11 and Figure 6.12 respectively. Excellent match between the Gr?neisen parameters calculated via anharmonicity approach and finite difference approach shows that our calculated third anharmonicity tensor A is trustable. This is a very necessary first step in order to accurately calculate the phonon life time and lattice thermal conductivity. 111 Figure 6.11 Mode Gr?neisen parameters of d-Si calculated from FDA and AA. Figure 6.12 Mode Gr?neisen parameters of Si 136 calculated from FDA and AA 112 6.4 Conclusion In a summary, we have investigated the causes of reduction of lattice thermal conductivity in Si 136 compared to d-Si. Section 6.2 shows that the flattering of the phonon dispersion (lowed phonon group velocity) in Si 136 leads to a reduction of thermal conductivity to 20% of d-Si. Section 6.3 calculated the lattice anharmonicity and analyzed the bond-stretching anharmonicity. The comparison of the relative anharmonicity of the two systems implied that reduction of phonon life time in Si 136 is also another cause (phonon life time of Si 136 is reduced by anther factor 70%). The Kinetic theory suggests that the lattice thermal conductivity of Si 136 is about 14% ( 20% 70%=?) of that of d-Si. And note that this ratio is a ball-park estimation, it provides a upper limit for ratio 136 Si dSi ? ? . And detail phonon life calculation is desired to give more rigorous comparison. 113 CHAPTER 7 LATTICE THEMRAL CONDUCTIVITY OF MgO 7.1 Introduction Thermal conductivity (? ) data of Earth?s constituent minerals are important for understanding any geophysical process that involves heat 89,90,91 . Although several rapid developments in experimental techniques were reported in recent years 92,93,94,95,96 , some pressure (P) and temperature (T) conditions of the Earth?s interiors (for example, T > 2300K or P > 100GPa) remain inaccessible for accurate measurement of ? at the current stage. Furthermore, the issue of contact associated errors for the thermal transport measurements has been raised and discussed 14 . The systematic errors of this type are especially important for accurately determining the pressure dependence in thermal transport properties. At the same time, little theoretical effort has been devoted to the first principles calculations of this important thermal transport property of minerals, including ideal crystalline minerals (i.e. containing no isotope/composition disorder, no isolated or extended defects, or no finite-size grain boundaries). Current understanding on lattice anharmonicity and its pressure dependence is limited. Using the statistical linear response theory (Green-Kubo formula), Cohen reported the first theoretical calculation of lattice 114 thermal conductivity in MgO 97 . The preliminary results reported in that study revealed an unusual pressure dependence of ? , and Cohen suggested that the phonon-phonon scattering at compression may behave differently at low and high pressure ranges. The interatomic potential adopted in that study was a non-empirical ionic VIB model. The equilibrium properties of MgO, such as thermal equation of state, predicted by this model are in reasonable agreement with experiment 98 . While equilibrium structural and energetic properties are largely determined by the harmonic inter-atomic interactions, (non-equilibrium) thermal transport properties are intrinsically influenced by lattice anharmonicity. Yet, the VIB model has not been tested systematically for the anharmonicity. Neither lattice anharmonicity nor phonon-phonon scattering is explicitly evaluated in the Green-Kubo type of calculations. Their effects are interpreted through the calculated correlation function of heat currents. To ensure convergence of time integrals of the correlation functions, the calculations of this class often require molecular dynamics (MD) simulation time of order of 1000 ps, a formidable computation task for accurate first principles techniques, such as fully self-consistent planewave methods. Recently, Oganov and Dorogokupets reported a study about the anharmonicity effects on the thermodynamic potentials of MgO using a first-principles method 99 . In addition to the conventional quasi-harmonic approximation (QHA) results (including both harmonic and anharmonic contributions), an additional correction term, whose magnitude scales as a function of T 2 , was estimated using MD simulations. Individual interatomic anharmonicity terms were not explicitly evaluated. Comparing to the calculation of heat current correlation functions, significantly less MD simulation time is needed to reasonably approximate the ensemble average. The authors reported that at ambient 115 pressure, the lattice anharmonicity evaluated with the QHA approach led to a noticeable overestimation in the lattice thermal expansion, an equilibrium thermal property that is believed to be closely related to lattice anharmonicity. This work provides a first-principles calculation of harmonic phonon spectra and 3 rd order lattice anharmonicity in MgO. Explicit calculation of both anharmonicity and mass-disorder induced phonon life time have been carried out. Lattice thermal conductivity at a wide range of temperature has then been calculated using kinetic transport theory. Finally two models have been proposed to estimate the pressure dependence of lattice thermal conductivity of MgO. The simple oxide MgO is considered to be an end-member component of the lower mantle. Probing the lattice anharmonicity and thermal conductivity of MgO is a precursor to studying more complex mineral structures and compositions relevant to the Earth. The direct derivation of harmonic force constant matrices and the 3 rd order anharmonicity tensors were carried out numerically using an efficient supercell finite-displacement (SFD) technique. The calculated lattice anharmonicity were carefully tested, and the pressure dependence is predicted. Results of the 3 rd order lattice anharmonicity are the necessary inputs for obtaining phonon-phonon scattering rates, and such rates are needed for evaluating phonon life time and lattice thermal conductivity based on non-equilibrium transport theories, such as kinetic transport theory or more sophisticated Boltzmann transport equation (BTE). We have adopted kinetic transport theory in our current study of lattice thermal conductivity. The preliminary results of lattice thermal conductivity over a wide range of temperature at ambient pressure are based on a 444? ? q-point sampling over the BZ. In addition to the anharmonicity induced phonon-phonon scattering, we also considered the mass disorder 116 (isotope effect) induced phonon-phonon scattering, and isotope effects on the lattice thermal conductivity at difference temperature will also be discussed. Explicit calculation of lattice thermal conductivity at high pressure is not presented but will be reported in the future. In this dissertation, we simply discussed the implication for the pressure dependence of lattice thermal conductivity based on the single relaxation time approximation (SRTA). A characteristic feature of these SRTA models is that the thermal conductivity of a mineral is simplified into two contributing terms: the harmonic dynamics related heat capacity and group velocity term and the anharmonic dynamics related relaxation time term. The pressure dependency of the harmonic term is evaluated with our first-principles phonon data, whereas the pressure effects on relaxation time are predicted based on some simple assumptions. These models will serve as the baseline models for our future data obtained based on the explicit kinetic transport calculations. 7.2 Crystal Structure and Equation of State Two-atom face-centered-cubic (fcc) unit-cell models were used in the EoS calculation. The Brillouin zone integration for electronic energy was calculated using the tetrahedron method over a 24 24 24? ? Monkhorst-Pack grid. Wave functions of valence electrons in Mg (2p 6 3s 2 ) atoms and O (2s 2 2p 4 ) atoms were expanded in a set of plane waves, while the core electrons (i.e. 1s 2 2s 2 in Mg and 1s 2 in O) were replaced with the Projector Augmented Wave (PAW) potentials. To achieve better numerical accuracy, the semi-core 2p electrons in Mg atoms were treated as valence electrons in our calculations. The adopted kinetic energy cutoff value for plane waves was 400eV. 117 The fitting parameters to 3 rd Birch-Murnaghan equation of states are listed in Table 7.1, experimental measurements and other theory results are also listed for comparison. Comparing to the experimental data at room temperature 100 , our static calculation underestimates the V 0 at 300K by 3%, while it overestimates the bulk modulus at 300K by 8%. The errors reported here are typical of calculations of this class, and our results are consistent with previous theoretical results 16,101,102,103 (Table 7.1). We further find that including free energy contribution due to lattice vibration further improves the agreement between LDA calculated and experimental V 0 and B 0 . 104 Table 7.1 Our LDA calculated static equilibrium properties of MgO fitted by 3 rd order Birch-Murnaghan equation of state, compared with previous theoretical results and experimental measurement. Experimental data was taken at room temperature. EOS Parameters LDA(this work) LDA a LDA b LDA c GGA d Expt. e V 0 (? 3 /MgO) 18.1 19.05 18.8 19.2 18.1 18.7 B 0 (GPa) 172.7 172.6 159 159.7 172 160.2 B? 4.2 4.0 4.30 4.26 4.09 3.99 7.3 Phonon and Gr?neisen Parameters A 128-atom super cell model with only ?-point sampling in the Brillouin zone was used in the calculation for H-F forces. Phonon calculation was based on the finite a Reference 99 b Reference 101 c Reference 102 d Reference 103 e Reference 100 118 difference method, which, however, can not handle well the optical modes for ionic systems. Thus additional LOTO correction has been made for optic phonons. The results of our LDA calculated phonon dispersion relation at 0 GPa and 68 GPa are plotted in Figure 7.1a and 7.2a respectively, and they are in excellent agreement with the available experimental data 105,106 , as well as two previously published theoretical result 101,107 using density perturbation functional theory (DPFT). Other phonon related properties, such as phonon density state (phDOS, Figure 7.1b and 7.2b), or phonon group velocity (Figure 7.1c and 7.2c), can be readily derived from the calculated phonon dispersion curves and/or dynamical matrices. Figure 7.1 LDA calculated harmonic phonon spectra, phonon density of state and group velocity of MgO at 0 GPa. Experiment data of phonon dispersion measured by Peckham (1967) and Sangster et al. (1970) are also shown for comparison. 119 Figure 7.2 LDA calculated phonon spectra, phonon density of state and group velocity of MgO at 68 GPa. Figure 7.3 LDA calculated Gr?neisen parameter as a function of temperature at 0 GPa. 120 Equilibrium thermodynamic properties of MgO at high-pressure have also been calculated within the statistical QHA. In this dissertation, I will focus only on Gr?neisen parameter ( th ? ), a dimensionless parameter that is commonly adopted in empirical models of lattice thermal conductivity because of its intrinsic connection with lattice anharmonicity 8,108 . The ? th of a mineral can be obtained based on the measured P-T-V EOS: () () th V V V PVP V UCT ? ? ? ?= ?? , (7.1) where () () VV V US CT TT ?? == is the heat capacity. Within the QHA, the bulk th ? can also be expressed in term of averaged phonon mode Gr?neisen parameters log ( ) () log i i dV gV dV ? ?? : (, ) () (, ) (, ) ii i th V cTVg V TV CTV ? = ? , (7.3) where i ? is the vibration frequency of the i th phonon mode , and / / 2 2 1 [( ( , ) ) ] 2 (, ) ( ) (1) iB kT iB kT ii i i B B dn T e cTV k dT k T e ? ??? ? + == ? h hh h represents the mode heat capacity of the corresponding phonon, and the total heat capacity Vi i Cc= ? sums contributions from each phonon. Our calculated th ? at the equilibrium volume is shown in Figure 7.3. The high temperature limit of th ? is estimated to be 1.39, about 10% smaller than the reported experimental value of 1.54 109 . The good agreement between our LDA result and 121 experimental measurement indicates that the LDA theory is capable to quantitatively describe lattice anharmonicity. Our theoretical results are not sensitive to the numerical approximations that we adopted to calculate mode i g parameters. 7.4 Lattice Anharmonicity Omitting all the 4 th and higher order of lattice anharmonicity, we further derived the full 3 rd order anharmonicity tensors at various chosen volumes/pressures. Because of the lattice symmetry of MgO, all the single-atom anharmonicity terms, i.e. the (, , )i j k indices of the tensor element ijk A are all from one single atom, are zeros. Two largest ijk A terms, around 24 3 eV ? at equilibrium volume, are found to be the bonding-stretching anharmonicity between two neighboring Mg-O atoms along the bonding direction. We refer the absolute values of these two largest anharmonicity terms as 2 OMg A and 2 OMg A respectively. Although it is well known that the charge distribution around the O 2- ions are much more sensitive to the surrounding ion configuration than that around the Mg 2+ ions, our results show that 22 1 () 2 OMg OMg AA A?= ? is only about 0.2 3 eV ? , i.e., less than 1% of the averaged value. Therefore, we define 22 1 () 2 stretching OMg OMg AAA?+ to represent the dominant term in the 3 rd order anharmonicity tensor. As shown in Figure 7.4, stretching A in MgO increases almost linearly with pressure (P), and at zero pressure log( ) stretch dA dP is estimated from linear fitting as 0.0153 GPa -1 . We find that A? decreases 122 with increasing pressure and turns negative above 30 GPa. The ||/ stretching AA? ratio remains relatively small for all the pressures studied in this paper. Figure 7.4 The pressure dependence of the average and difference (inset) of the two largest anharmonic tensor elements. Other anharmonic ijk A elements of MgO are at least one order of magnitude smaller than stretching A , and all the terms involving the atoms and their 3 rd and beyond nearest neighbors are negligible compared to the major terms. As a simplification, we only keep all the terms coming from the atoms and their 1 st and 2 nd neighbors in the calculated A tensors. To test the accuracy of the anharmonicity tensors calculated with our new algorithm, we calculated the phonon mode Gr?neisen parameters from both finite difference approach (Equ. (3.9)) and anharmonicity approach (Equ. (3.10)), and their comparison is shown in Figure 7.5. Even though the match is not as good as d-Si 123 and Si 136 , but this is based on the F-H force calculated from first-principles, thus we consider it quite satisfactory. Figure 7.5 Comparison of Phonon mode Gr?neisen parameters calculated with both finite difference method and the 3 rd order lattice anharmonicity method. 7.5 Isotope Effect Phonon-scattering processes in general can be divided into intrinsic processes and extrinsic processes. Lattice anharmonicity induced phonon scattering is the dominant intrinsic process in insulators and most semiconducting materials. The extrinsic processes include the phonon scattering at all sorts of crystal defects and crystal surfaces. In this dissertation, lattice anharmonicity for MgO has been presented and discussed in section 7.4, and as far as the extrinsic processes are concerned, only isotope effects will be considered. Many studies have shown that there is a significant contribution to the scattering from the variation in isotope mass, which can modify the thermal conductivity 124 by an appreciable degree. The first experiment studying the isotope effect was carried by Geballe and Hull in 1958, who observed an increase in thermal conductivity by 15% at room temperature when the concentration of the 74 Ge isotope was increased from 36.5% (natural) to 95.8%. 110 Isotope effect on the thermal conductivity in silicon was then studied by Capinski et al. in 1997. 111 They suggested an increase of thermal conductivity in isotope pure silicon by 60%, as compared to silicon of natural isotopic abundance. In this dissertation, we will also discuss the isotope effect of MgO on its thermal conductivity at both low temperature and high temperature. The scattering rate or the reciprocal of the relaxation time iso ? of the phonons due to single scattering by the isotopes is given by 112 : 2 1 22 ', 2(,) (,) ( (,) (',)) | (,)| | (',)| 2 iso qj qi qi qi q jgeqi eqj ?? ? ? ?? ??? ? ?? =??? ?? ?? ?? v v hvvvvv hh h ,(7.4) and 2 ()[1 ()/ ()] kk k gf mm ? ? ??=? ? , (7.5) where h is Planck constant, (,)qi? v and ( , )eqi ? v are the eigen-frequency and eigen- vector of phonon mode (,)qi v , () k f ? is the fraction of th k isotope of atom ? that has mass () k m ? , and ()m ? is the average mass of atom ? . g ? for oxygen and magnesium have been calculated and listed in Table 7.2. 125 Table 7.2 Mass deviation of atom oxygen (g O ) and magnesium (g Mg ). g O g Mg 0.336?10-4 7.398?10-4 The overall phonon relaxation time is then given by: 11 1 (,) (,) (,) anh iso qi qi qi?? ? ?+ v vv . (7.6) In order to estimate the isotope effects on the thermal conductivity of MgO, we approximate the phonon life time from anharmonicity induced phonon-phonon scattering to be eff ? such that 11 1 (,) (,) eff iso qi qi??? ?+ v v . (7.7) eff ? was calculated via 2 , 1 (,) (,) 3 eff Vg qiq cqiV qi N ? ? = ? vv v v , where ? was using experiment measured thermal conductivity for natural MgO 8 . Both V c and g V were calculated at a 20 20 20??q-grids. eff ? at 500K and 2000K are reported in Table 7.3. Table 7.3 Effects of isotope on the lattice thermal conductivity at 500K and 2000K. T (K) ? natural (W/K/m) ? effective (ps) ? isotope-pure (W/K/m) % increase 500 35 3.015 39.4 12.6 2000 7 0.495 7.15 2.1 126 With ( , ) iso qi? v , which can be evaluated from Equ. (7.4), and eff ? , we then calculated the overall phonon relaxation time, and further predicted the thermal conductivity of isotope-pure MgO ( isotope pure ? ? ). At 500K, isotope-pure MgO has 12.6% higher thermal conductivity than natural MgO, while at 2000K only 2.1% increase in the thermal conductivity (see Table 7.3). Our results indicate that the isotope effect on the thermal conductivity is more important at low temperature, while it is almost negligible at high temperature. 7.6 Thermal Conductivity Explicit calculations for both anh ? and iso ? at 444? ? were carried out using first- principles methods. By considering both intrinsic phonon-phonon scattering and extrinsic isotope induced phonon scattering, Thermal conductivity over a wide rang of temperature at ambient pressure has been calculated using kinetic transport theory. At room temperature, our calculated thermal conductivity is about 51 W/K/m, in good agreement with experimental measurement 54 W/K/m by Slack in 1962. 113 Figure 7.6 shows the temperature dependence of our calculated thermal conductivity of MgO at ambient pressure as well as some available experimental work 114,113,115 . Overall, our theory results are comparable with experimental measurements. 127 Figure 7.6 Temperature dependence of lattice thermal conductivity ? at 0 GPa; ? was calculated using First-Principles method with a 444? ? q-point sampling over the Brillouin zone. Experimental measurements are also cited for comparison. Figure 7.7 Estimated pressure dependence of lattice thermal conductivity in MgO at 500K for model I and model II. 128 Figure 7.8 Estimated pressure dependence of lattice thermal conductivity in MgO at 2000K for model I and model II. Recent experimental results are shown in symbol plus (+). We further provide a preliminary estimation on the pressure dependence of lattice thermal conductivity in MgO based on a simple kinetic transport model within single relaxation time approximation (SRTA) (See Equ. (6.1), which characterizes the effect of anharmonicity-induced phonon scattering with just one single parameter: the effective phonon life-time (or relaxation time) eff ? . Two models are proposed here for () (0) P P ? ? = . Model I assumes that eff ? is pressure independent, i.e., () 1 (0) eff eff P P ? ? = = the pressure dependence of the thermal conductivity comes from the harmonic terms and the isotope 129 effects. Model II suggests eff ? is inversely proportional to 2 (/) stretching A B , i.e., 2 2 ( ) ( ( 0) / ( 0)) (0) ( ()/() eff stretching eff stretching PAPBP? ? == = = , where B is the Bulk Modulus. Pressure dependence of thermal conductivity in MgO proposed from these two models at 500K and 2000K are plotted in Figure 7.7 and 7.8. eff ? at zero pressure were taken from Table 7.3. Results of recent experiment at 2000K 15 are better fitted to Model II. Both models provide a rather linear relationship for () (0) P P ? ? = , ie., () 1 (0) P P P ? ? ? =+ = . Slope ? of these two models at 500K and 2000K are listed in Table 7.4. Table 7.4 List of ? in ()/( 0) 1PP P?? ?= =+ estimated from two models at 500K and 2000K respectively. T(K) Model ? T(K) Model ? I 0.0059 I 0.0071 500K II 0.0244 2000 II 0.0296 7.7 Conclusion Using an efficient and accurate super-cell finite-displacement algorithm, we studied the harmonic and anharmonic lattice dynamics in MgO at high-pressure. The explicitly calculated 3 rd anharmonic interaction terms were reported, and the calculated ijk A can well reproduce the mode Gr?neisen parameters. The estimated pressure coefficient of the dominant stretch A term is about 0.0153 GPa -1 . In addition to the intrinsic 130 anharmonicity induced phonon-phonon scattering, variation in the mass of isotopes also scatters phonons, which is also a factor causing the reduction of thermal conductivity. Such isotope effect on the lattice thermal conductivity has been studied, isotope-pure MgO will have 12.6% higher thermal conductivity than natural MgO at 500K. Our results also indicate that isotope effect is important at low temperature and becomes negligible at high temperature. With both anharmonicity and isotope effects considered, we presented the preliminary results of temperature dependence of lattice thermal conductivity in MgO at ambient pressure, which shows good agreement with previous experimental measurements. Two models were then proposed for the estimation of the pressure dependence of lattice thermal conductivity. Explicit calculation of lattice thermal conductivity at high pressure shall be carried out, and a finer q-grid sampling over the BZ is desired to more accurately integrate the BZ. 131 CHAPTER 8 SUMMARY AND FUTURE WORK In this dissertation I have presented the details of our recent methodology development on first-principles materials simulation and modeling and the simulation results on two interesting materials systems. The methodology development consists two categories. First, we have extended the standard first-principle total energy and force calculations to the calculations of harmonic (2 nd order) force constant matrices and 3 rd order lattice anharmonicity tensors using an efficient super-cell finite-displacement algorithm. We have carried a series of detailed comparison between our calculated (harmonic) phonon spectra and the results reported by other groups using linear response density perturbation functional theory (DPFT) method, as well as available experimental measurements. We find that our results are in excellent agreement with experiments and results calculated with the DPFT method. Our SC-FD algorithm is currently implemented with the VASP code, but it can also be easily implemented with any other first-principle codes with any forms of pseudo-potentials and electron wave basis sets. This algorithm is ideal for parallel computation using low-cost Beowulf computer clusters. We have further proposed a paired displacement approach to calculate 3 rd order lattice 132 anharmonicity. The new algorithm has been implemented with both simple empirical Tersoff potentials and the first-principles DFT method. We also proposed a robust testing scheme and investigated the sumrule to ensure the numerical accuracy of the calculated 3 rd order anharmonicity tensor. Next, we have adopted statistical ensemble theory of phonons to calculate and simulate materials properties at finite temperature. The results derived from the first- principles calculated force constant matrices are utilized within the quasi-harmonic approximation (QHA) to predict fundamental thermodynamic potentials as functions of temperature (T) and pressure (P). The consequently derived (T,P) phase diagrams and equilibrium thermal properties, such as lattice thermal expansion, and/or heat capacity, were compared directly with experiments. Meanwhile, we have also applied the quantum scattering theory to explicitly calculate the phonon life time, and further implemented a simple kinetic transport theory to predict non-equilibrium thermal transport properties. Investigation of contribution of harmonic and anharmonic components in kinetic transport equation to the lattice thermal conductivity has been decoupled under the single phonon relaxation time approximation. Using the newly developed computational methods, we have studied the thermal properties of the novel guest-free clathrate polymorph of silicon (Si 136 ) based on first- principles calculations combined with experimental X-ray and Raman scattering measurements. The Si 136 clathrate is metastable compared with the d-Si phase at ambient P and T. Theory indicates that it becomes thermodynamically stable within a negative pressure regime, at P = -2 to -4 GPa. The dP tr /dT Clapeyron slope is estimated as 7?10 -4 GPa/K from ab initio calculations. Although it has been shown previously that some 133 properties, such as the electronic band gap, are critically dependent upon the lattice expansion between the diamond-structured and ?expanded-volume? clathrate polymorphs, our current studies reveal that the thermal properties of the two phases involving long- wavelength phonons are similar to each other. The vibrational properties of Si 136 phonons are similar to those of d-Si, and they can be understood in terms of Brillouin zone reduction following the unit cell expansion between d-Si and Si 136 . The phonon modes in the two phases also have very similar characteristics. We find that the coefficients of thermal expansion in the two Si phases are comparable in our studies, which suggests that it is less likely that the significant reduction of lattice thermal conductivity in clathrate materials is mainly caused by any large increases of anharmonic lattice interactions in clathrate systems. We also adopt the Green-Kubo formulism to calculate the thermal conductivity of d-Si and Si 136 using classical molecular dynamics with empirical Tersoff potential. Our calculation reveals that the thermal conductivity of Si 136 is about 10% of that of d-Si. According to the kinetic theory, this reduction can be attributed to lower heat capacity, or lower group velocity (flattering phonon spectra), or lower phonon life time (large lattice anharmonicity). To investigate the contribution of the low conductivity in Si 136 , third order lattice anharmonicity of both d-Si and Si 136 was further studied using the same empirical potential. Our results shows that phonon life time of Si 136 is about 70% of that of d-Si, and the reduced thermal conductivity in Si 136 mainly comes from the phonon flattering or reduced group velocity (a reduction of 80%) as suggested from first- principles thermal calculations. 134 Using the efficient and accurate super-cell finite-displacement algorithm, we then studied the harmonic and anharmonic lattice dynamics in MgO at high-pressure. The explicitly calculated 3 rd anharmonic interaction terms are reported, and the calculated ijk A can well reproduce the mode Gr?neisen parameters. The estimated pressure coefficient of the dominant stretch A term is about 0.0153 GPa -1 . Both intrinsic lattice anharmonicity and extrinsic isotope induced phonon scattering have been taken into consideration during the calculation of lattice thermal conductivity. Explicit calculation of phonon life time from two scattering mechanisms has been carried out. Kinetic transport equation was then adopted to evaluate the lattice thermal conductivity over a wide range of temperature at ambient pressure. Preliminary results of lattice thermal conductivity were based on 444?? q-grids BZ integration. At room temperature, our calculated lattice thermal conductivity is 51 W/K/m, in a good agreement with experimental measurement 54 W/K/m. The overall temperature trend is consistent with available experiments. Explicit calculation of thermal conductivity at high pressure is not yet done in this dissertation. Instead, two models have been proposed to estimate the pressure dependence. At 2000K, model II suggests the pressure coefficient to be 0.0296, while model I suggests a rate 0.0071, and model II is better fitted with the experiment. The immediate extension of the research work discussed in this dissertation is to calculate the thermal conductivity at high pressure explicitly. And the convergence of predicted lattice thermal conductivity on the size of q-point sampling deserves a much more careful investigation. 135 REFERENCES 1 G. A. Slack, Thermoelectric Materials-New Directions and Approaches, edited by T. M. Tritt, G. Mahan, H. B. lyon, Jr., and M. G. Kanatzidis(Materials Research Society, Pittsburgh, PA, 1997), Vol. 407, p. 47; G. A. Slack, in CRC Handbook of Thermoelectric, edited by D. M. Rowe( Chemical Rubber Corp., Boca Raton, FL, 1995), p. 407. 2 G. S. Nolas, J. L. Cohn, G. A. Slack, S. b. Schujman, Appl. Phys. Lett. 73, 178 (1998). 3 J. Dong, O. F. Sankey, and C. W. Myles, Phys. Rev. Lett. 86, 2361 (2001). 4 J. Gryko, P. F. McMillan, R. F. Marzke, G. K. Ramachandran, D. Patton, S. K. Deb, O. F. Sankey, Phys. Rev. B 62, R7707 (2000). 5 G. S. Nolas, M. Beekman, J. Gryko, G. A. Jr. Lamberton, T. M. Tritt, P. F. McMillan, Appl. Phys. Lett. 82, 910 (2003). 6 O. Madelung, Semiconductors - Basic Data (New York, 1996). 7 H. B. Callen and T. A. Welton, Phys. Rev. 83, 34 (1951). 8 A.M. Hofmeister, Science 283, 1699 (1999). 9 A.M. Hofmeister, Thermal conductivity of the Earth. In: Schubert, G. (Ed). Treattise on Geophysics, Vol 2: Mineral Physics (2007). 10 A.M. Hofmeister and D.A. Yuen, J. Geodyn. 44, 186 (2007). 11 S. Andersson and G. Backstrom, Rev. Sci. Instrum. 57, 1633 (1986). 12 T. Katsura, Geophys. J. Int. 122 (1), 63 (1995). 136 13 A.M. Hofmeister, Phys. Chem. Minerals 33, 45 (2006). 14 A.M. Hofmeister, PNAS 104, 9192 (2007). 15 P. Beck, A.F. Goncharov, Viktor. Struzhkin, B. Militzer, H.-K. Mao and R.J. Hemley, Appl. Phys. Lett. 91, 181914 (2007). 16 A.R. Oganov and P.I. Dorogokupets, Phys. Rev. B 67, 224110 (2003). 17 M. Born and R. Oppenheimer, Annalen der Physik, 84, 457 (1927). 18 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 19 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 20 D. C. Langreth and J. P. Perdew, Phys. Rev. B 21, 5469 (1980); D. C. Langreth and M. J. Mehl, Phys. Rev. B 28, 1809 (1983); J. P. Perdew and Y. Wang, Phys. Rev. B 33, 8800 (1986); J. P. Perdew, Phys. Rev. B 33, 8822 (1986). 21 G. Kresse and J. Furth?mller, Comput. Mater. Sci, 6, 15 (1996); G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993); G. Kresse and J. J. Furth?mller ibid. 54, 11169 (1996). 22 R. Hamann, M. Schluter, and C. Chiang, Phys. Rev. Lett. 43, 1494 (1979). 23 D. Vanderbilt, Phys. Rev. B 41, 7892 (1990); K. Laasonen, R. Car, C. Lee, and D. Vanderbilt, ibid. 43, 6796 (1991); K. Laasonen, A. Pasquarello, R. Car, C. Lee, and D. Vanderbilt, ibid. 47, 10142 (1993); G. Kesse and J. Hafner, Phys. Rev. B 48, 13115 (1993); J. Phys.: Condens. Matter 6, 8245 (1994). 24 For example, J. Ihm, A. Zunger and M.L. Cohen, J. Phys. C 12, 4409 (1979). 25 International Tables for Crystallography, Volume A, edited by Th. Hahn. Reidel, Publishing Company, Dordrecht, Boston (1996). 137 26 M. Laue, Ann. Phys. 41, 989 (1913). 27 W.L. Bragg, Proceedings of the Cambridge Philosophical Society, 17, 43 (1914). 28 H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976). 29 http://www.abinit.org/. 30 http://cms.mpi.univie.ac.at/vasp/vasp/vasp.html. 31 http://www.uam.es/departamentos/ciencias/fismateriac/siesta/. 32 http://www.quantum-espresso.org/. 33 F. Birch, Phys. Rev. 71, 809 (1947); F.D. Murnaghan, Proceedings of the National Academy of Sciences, 30, 244 (1944). 34 K. Kunc and R.M. Martin, Phys. Rev. Lett. 48, 406 (1982). 35 D. A. Broido, A. Ward and N. Mingo, Phys. Rev. B 72 014308 (2005). 36 B.K. Ridley, Theory of Transport Properties of Semiconductor Nanostructures, edited by Eckehard Schc6ll, Springer (1997), pp 357. 37 E. Fermi, Nuclear Physics. University of Chicago Press (1950). 38 For example, M. T. Dove, Introduction to Lattice Dynamics. Cambridge University Press (1993). 39 R. Kubo. Rep. Prog. Phys. 49, 255 (1986). 40 R.J. Hardy, Phys. Rev. 132, 168 (1963). 41 R.A Street, Hydrogenated Amorphous Silicon. Cambridge University Press (1991). 42 S. H. Tolber, A. B. Herhold, L. E. Brus, and A. P. Alivisatos, Phys. Rev. Lett. 76, 4384 (1996). 138 43 J. Z. Hu, L. D. Merkle, C. S. Menoni, I. L. Spain, Phys. Rev. B 34, 4679 (1986); S. J. Duclos, Y. K. Vohra, and A. L. Ruoff, Phys. Rev. B 41, 120021 (1990); M. I. McMahon, and R. J. Nelmes, Phys. Rev. B 47, 8337 (1993); M. Imai, M. Mitamura, K, Yaoita, and K. Tsuji, High Pres. Res. 15, 167 (1996). 44 J. S. Kasper, P.Hagenmuller, M. Pouchard, and C. Cros. Science 150, 1713 (1965). 45 G. K. Ramachandran, J. Dong, J. Diefenbacher, J. Gryko, R. F. Marzke, O. F. Sankey and P. McMillan, J. Solid State Chem. 145, 716 (1999). 46 S. B. Roy, K. E. Sim, and A. D. Caplin, Philos, Mag. B 65, 1445 (1992). 47 H. Kawaji, H. Horie, S. Yamanaka, and M. I shikawa, Phys. Rev. Lett. 74, 1427 (1995). 48 D. Conn?table, V. Timoshevskii, B. Masenelli, J. Beille, J. Marcus, B, Barbara, A. M. Saitta, G. ?M. Rignanese, P. M?linon, S. Yamanaka, and X. Blase, Phys. Rev. Lett. 91, 247001 (2003). 49 G. S. Nolas, J. L. Cohn, G. A. Slack, S. b. Schujman, Appl. Phys. Lett. 73, 178 (1998). 50 J. L. Cohn, G. S. Nolas, V. Fessatidis, T. H. Metcalf, G. A. Slack, Phys. Rev. Lett. 82, 779 (1999). 51 G. S. Nolas, Proc. Chemistry, Physcis, and Materials Science of Thermoelectric Materials: Beyond Bismuth Telluride ( Traverse, MI, United States, 2003). 52 G. S. Nolas, M. Beekman, J. Gryko, G. A. Jr. Lamberton, T. M. Tritt, P. F. McMillan, Appl. Phys. Lett. 82, 910 (2003). 53 T. M. Tritt, Science 283, 804 (1999). 54 D. G. Cahill and R. O. Pohl, Annu. Rev. Chem. 39, 93 (1988). 139 55 G. S. Nolas, J. ?M. Ward, J. Gryko, L. Qiu, and M. A. White, Phys. Rev. B 64, 153201 (2001). 56 A. Bentien, M. Christensen, J. D. Bryan, A. Sanchez, S. Paschen, F. Steglich, G. D. Stucky and B. B. Iversen, Phys. Rev. B 69, 045107 (2004). 57 V. Keppens, M. A. McGuire, A. Teklu, C. Laermans, B. C. Sales, D. Mandrus, B. C. Chakoumakos, Physica B 316, 95 (2002). 58 J. S. Tse, K. Uehara, R. Rousseau, A. Ker ,C. I. Ratcliffe, M. A. White and G. Mackay, Phys. Rev. Lett. 85, 114 (2000); ibid, Phys. Rev. Lett. 86, 4980 (2001). 59 L. Qiu, I. P. Swainson, G. S. Nolas and M. A. White, Phys. Rev. B 70, 035208 (2004). 60 J. S. Kasper, P.Hagenmuller, M. Pouchard, and C. Cros. Science 150, 1713 (1965). 61 G. K. Ramachandran, J. Dong, J. Diefenbacher, J. Gryko, R. F. Marzke, O. F. Sankey and P. McMillan, J. Solid State Chem. 145, 716 (1999). 62 S. Bobev, S. C. Sevov, J. Am. Chem. Soc., 123, 3389 (2001). 63 G. K. Eamachandran, P. F. McMillan, J. Dong, J. Gryko, and O. F. Sankey, Mat. Res. Soc. Symp. Proc. 626, Z13.2 (2000). 64 O. Madelung, Semiconductors - Basic Data (New York, 1996). 65 G. K. Ramachandran, P. F. McMillan, S. K. Deb, M. Somayazulu, J. Gryko, J. Dong and O. F. Sankey, J. Phys.: Condens. Matter 12 (2000). 66 J.J. Hall, Phys. Rev. 161, 756 (1967). 67 S. Wang, Y. Hsu, J. Pu, J. C. Sung, and L. G. Hwa, Mater. Chem. Phys. 85, 432 (2004). 68 K. Moriguchi, S. Munetoh, A. Shintani and T. Motooka, Phys. Rev. B 64, 195409 (2001). 140 69 S. Wei, C. Li and M. Y. Chou, Phys. Rev. B 50, 14587 (1994). 70 Y. Guyot, L. Grosvalet, and B. Champagnon, Phys. Rev. B 60, 14507 (1999). 71 J. Dong, O. F. Sankey, and G. Kern, Phys. Rev. B 60, 950 (1999). 72 M. Wilson, P. F. McMillan, Phys. Rev. Lett. 90, 135703 (2003). 73 H. R. Shanks, P. D. Maycock, P. H. Sidles, and G. C. Danielson, Phys. Rev. 130, 1743 (1963). 74 P. Flubacher, A. J. Leadbetter and J. A. Morrison, Philos. Mag. 4, 273 (1959). 75 G. Nolas and T. Tritt, private communication. 76 C. Z. Wang, C. T. Chan and K. M. Ho, Phys. Rev. B 39, 8586 (1989). 77 J. L. Cohn, G. S. Nolas, V. Fessatidis, T. H. Metcalf, G. A. Slack, Phys. Rev. Lett. 82, 779 (1999). 78 G. S. Nolas, M. Beekman, J. Gryko, G. A. Jr. Lamberton, T. M. Tritt, P. F. McMillan, Appl. Phys. Lett. 82, 910 (2003). 79 G. S. Nolas, J. ?M. Ward, J. Gryko, L. Qiu, and M. A. White, Phys. Rev. B 64, 153201 (2001). 80 A. Bentien, M. Christensen, J. D. Bryan, A. Sanchez, S. Paschen, F. Steglich, G. D. Stucky and B. B. Iversen, Phys. Rev. B 69, 045107 (2004). 81 V. Keppens, M. A. McGuire, A. Teklu, C. Laermans, B. C. Sales, D. Mandrus, B. C. Chakoumakos, Physica B 316, 95 (2002). 82 J. S. Tse, K. Uehara, R. Rousseau, A. Ker ,C. I. Ratcliffe, M. A. White and G. Mackay, Phys. Rev. Lett. 85, 114 (2000); ibid, Phys. Rev. Lett. 86, 4980 (2001). 83 L. Qiu, I. P. Swainson, G. S. Nolas and M. A. White, Phys. Rev. B 70, 035208 (2004). 141 84 S. G. Volz and G. Chen, Phys. Rev. B 61, 2651 (2000). 85 J. Tersoff, Phys. Rev. B 39, 5566 (1989). 86 . Clarke D. Brown. Mol. Phys., 51:1243, 1983. 87 R. Bracewell, The Fourier Transform and Its Applications, 3rd ed. New York: McGraw-Hill, pp. 108-112 (1999). 88 W.S. Capinski, H.J. Maris, E. Bauser, I. Silier, M. Asen-Palmer, T. Ruf, M. Cardona and E. Gmelin, App. Phys. Lett. 71, 2109 (1997). 89 A.M. Hofmeister, Science 283, 1699 (1999). 90 A.M. Hofmeister, Thermal conductivity of the Earth. In: Schubert, G. (Ed). Treattise on Geophysics, Vol 2: Mineral Physics (2007). 91 A.M. Hofmeister and D.A. Yuen, J. Geodyn. 44, 186 (2007). 92 S. Andersson and G. Backstrom, Rev. Sci. Instrum. 57, 1633 (1986). 93 T. Katsura, Geophys. J. Int. 122 (1), 63 (1995). 94 A.M. Hofmeister, Phys. Chem. Minerals 33, 45 (2006). 95 A.M. Hofmeister, PNAS 104, 9192 (2007). 96 P. Beck, A.F. Goncharov, Viktor. Struzhkin, B. Militzer, H.-K. Mao and R.J. Hemley, Appl. Phys. Lett. 91, 181914 (2007). 97 R.E. Cohen, R. E. Reviews of High Pressure Science and Technology 7, 160 (1998). 98 I. Inbar and R.E. Cohen, Geophys. Res. Lett. 22, 1533 (1995). 99 A.R. Oganov and P.I. Dorogokupets, Phys. Rev. B 67, 224110 (2003). 100 S. Speziale, C.-S. Zha, T.S. Duffy, R.J. Hemley and H.-K Mao, J. Geophys. Res. [Solid Earth] 106, 515 (2001). 142 101 B.B. Karki, R.M. Wentzcovitch, S. de Gironcoli and S. Baroni, Phys. Rev. B 61, 8793 (2000). 102 M. J. Mehl, R.E. Cohen and H. Krakauer, J. Geophys. Res. 93, 8009 (1988). 103 P.I. Dorogokupets and A. R. Oganov, Phys. Rev. B 75, 024115 (2007). 104 X. Tang, J. Dong, P. Hutchins, O. Shebanova, J. Gryko, P. Barnes, J.K. Cockcroft, M. Vickers and P.F. McMillan, Phys. Rev. B 74: 014109 (2006). 105 G. Peckham Proc. Phys. Soc. London 90, 657 (1967). 106 M.J.L. Sangster, G. Peckham and D.H. Saunderson, J. Phys. C; Solid State Phys. 3, 1026 (1970). 107 A.R. Oganov, M.J. Gillan and G.D. Price, J. Chem. Phys. 118, 10174 (2003). 108 J. Fabian and P.B. Allen, Phys. Rev. Lett. 79, 1885 (1997). 109 D. Isaak D. Anderson O.L., T. Goto, Phys. Chem.. Miner. 16, 704 (1989). 110 T. H. Geballe and G. W. Hull, Phys. Rev. 110, 773 (1958). 111 W. S. Capinski, H. J. Maris, E. Bauser, I. Siljer, M. Asen-Palmer, T. Ruf, M. Cardona, and E. Gmelin, Appl. Phys. Lett. 71, 2109 (1997). 112 Shi-ichiro Tamura, Phys. Rev. B 27, 858 (1983). 113 Glen. A. Slack. Phys. Rev. 126 427 (1962). 114 F. R. Charvat and W. D. Kingery, J. Am. Ceram. Soc. 40, 306 (1957). 115 J. F. Schatz and G. Simmons, J. Geophys. Res. 77, 6966 (1972). 143 APPENDIX A 3 NOTATION FOR SPACE GROUP AND POINT GROUP Possible point operations are rotation, n C in Sch?nflies notation (S notation) or n in Hermann-Mauguin notation (H-M notation), representing a rotation of 360/n degree about an axis; n can only be 1,2,3,4 and 6 in crystal structure in order for the rotation to be compatible with the translational symmetry; reflection, denoted by , , hvd ? ?? in S notation or m in H-M notation, representing a mirror operation about a plane. The subscript h means the mirror plan which is perpendicular to the principle rotation axis, v the mirror plan includes principal rotation axis, and d the mirror plane including the principal rotation axis and also bisects the angle between the other two 2-fold rotation axes; improper rotation, denoted by n S in S notation or n in H-M notation, representing a combination of rotation n C by 360/n degree followed by reflection h ? in a plan normal to the rotation axis. The other two operations are Identity and Inversion, which are special cases (n=1) of rotation and improper rotation respectively. The point, axis, and plane mentioned above are called symmetry elements. Each point group possesses different numbers of symmetry operation, it can has as many as 48 symmetry operations ( h O ), and 144 as low as 1 symmetry operation ( 1 C ). All the possible symmetry operations of each point group can be found in International Table of Crystallography. Two additional symmetry operations are the combination of rotation or reflection with a translation less than the unit cell size. A glide plane is a reflection in a plane followed by a translation parallel to the plane, denoted by ,,,abcn or d , where ,,abc represents the gliding direction along axis a, b, c by half of the cell length in that direction, n represents the gliding along the half of the diagonal of a face, and d represents the gliding along a quarter of the face or space diagonal of the unit cell. A screw axis is a rotation about an axis, followed by a translation along the direction of the axis. It is denoted by a number n followed by a subscript m , n represents the degree of right-hand rotation 360 / n about the rotation axis, and m represents a following translation of /mn of cell length along that axis direction. Since the translation for screw axis is within the unit cell, m must be smaller than n . Usually we say space group #227 is 3Fd m , which belongs to point group 7 h O , we are talking about the notation of the space group. Space groups commonly use HM notation while point groups often use S notation. HM notation for space groups starts with a Bravais lattice descriptor indicating the centering of the lattice followed by up to 3 set of notation for possible point operations, glide planes or screw axes. The centering of the lattice is expressed by P for primitive cell (R for rohmbohedra cell only) , C for base- centered cell, F for face-centered cell and I for body-centered cell. The order of point symmetry operation in the HM notation is not random. The first symbol denotes the symmetry along the major axis, and the second along axis of secondary importance, and the third along the tertiary important direction (see Table A.I). For cubic system, the 145 primary direction is [100]/[010]/[001] since they are equivalent, the secondary direction is the diagonal direction [111] and the tertinary direction is [110]. 3Fd m is the short HM notation of #227, the full notation is 1 42 3F dm , the symmetry operation along each direction can been easily seen from the full notation. 7 h O is the S notation for its point group. One can also derive the HM notation of the point group from its space group by replacing the glide plane by mirror and screw axis by rotation. Here we have d glide in space group notation, after replace it by mirror; we have the HM notation 3mmfor point group h O . Table A.I Primary, secondary, and ternary directions commonly used for 7 crystal systems. Symmetry direction Crystal system Primary Secondary Ternary Cubic [100]/[010]/[001] [111] [110] Hexagonal/Trigonal [001] [100]/[010] [120]/ [1 10] Tetragonal [001] [100]/[010] [110] Orthorhombic [100] [010] [001] Monoclinic [010] Triclinic None 146 Knowing the meaning of the space group symbolism, we have better idea about the International Tables for Crystallography, which contains full information about the 230 space groups. The more important information listed in this table is the Wyckoff section for each space group. As an example, a complete list of Wyckoff sites for space group 227 has been cited here for reference (Table A.II), where the coordinates are partially listed. Note that conventional cubic unit cell is used. Table A.II Wyckoff sites for space group 227 (only partial coordinates are listed) Coordinates Multiplicity Wyckoff letter Site symmetry (0,0,0)+ (0,1/2,1/2)+ (1/2,0,1/2)+ (1/2,1/2,0)+ 192 i 1 (x,y,z) ?. 96 h ..2 (1/2,y,-y+1/4), ?. 96 g ..m (x,x,z), ?? 48 f 2.mm (x,0,0), ?? 32 e .3m (x,x,x), ?? 16 d .-3m (5/8,5/8,5/8), ?? 16 c .-3m (1/8,1/8,1/8), ?? 8 b -43m (1/2,1/2,1/2) (1/4,1/4,1/4) 8 a -43m (0,0,0) (3/4,1/4,3/4) As one can see from Table A.II, the Wyckoff site contains the information of multiplicity, Wyckoff letter, site symmetry and also coordinates. Multiplicity represents 147 the number of equivalent atoms per unit cell, Wyckoff letter is simply a letter, starting with a at the bottom position and continuing upwards alphabetically; the site symmetry gives all symmetry operations that map a point onto itself; it is the subgroup of the point group to which the space group under consideration belongs. Site symmetry symbol displays the same sequence of symmetry directions as the space group symbol. Coordinates of all the points that go back to themselves under this site symmetry operation are listed in the next column. For centered space groups, the centering transformations are listed above the coordinate triplets. There are 9 Wyckoff sites in space group 227. Any structure belonging to this space group will have one or more of these Wyckoff sites. For example, d-Si and its isotope phase silicon clathrate type II (Si 136 ) both belongs to space 227, but they have different combination of Wyckoff sites. In d-Si system, there is one type of site (8a) occupied by Si, while Si in Si 136, occupies 8a, 32e and 96g three different types of site, which lead to 136 as the total number of atoms in the unit cell. Let?s take 8a site as an example to illustrate the meaning of Wyckoff site symmetry. As showed in table, 8a site has site symmetry 43m? , representing an improper 4-fold rotational symmetry along primary direction (001), a 3 fold rotational symmetry along secondary direction [111] and mirror symmetry in the ternary direction which is automatically satisfied. The complete Wyckoff position for this site can be obtained by adding four centering translations (0,0,0), (0,1/2,1/2), (1/2,0,1/2) and (1/2,1/2,0) listed on the top of the Table A.II to the two coordinate triplets (0,0,0) and (3/4,1/4,3/4) in the 8a row of the Table A.II. Thus we have positions for 8 atoms in the unit cell, any of which will go back to itself after the site symmetry operation 43m? . 148 Similarly, we can derive the positions for 32e and 96g sites. Note that these two sites have unknown parameter, which can be determined by X-ray diffraction experiments. 149 APPENDIX B NORMAL MODE SYMMETRY REPRESENTATION ANALYSIS Phonon calculation predicts the dispersion curve, as well as the phonon modes at Gamma point. However, not all the gamma phonon modes can be detected in experiments, such as Infrared spectroscopy and Raman Scattering, which both characterize the vibrations of chemical bonds. The mechanism of these experiments determines the vibrational selection rules for crystals. What is of our interest is how to characterize all the vibrational modes calculated by our first-principles method, and in the meantime, determine which modes are IR active and which modes are Raman active. For a vibration to be Raman active modes, the polarizability of the crystal will change with the vibration motion, for a vibration to be IR active, the dipole moment of the crystal will change with the vibration motion. By doing this analysis, we are not only able to compare our phonon frequency with existing experimental values, but also predict those without experiment results. 150 In order to get the irreducible representation of the phonon modes, and further the optically active vibration modes, Ref.I is a good source of help. Here is a simple illustration about how to derive them for crystals provided Site Symmetry Table (Appendix I in Ref. I), Character Table (Appendix II in Ref. I) and Correlation Table (Appendix III in Ref. I). In some situation, Table 14 in Ref. I could also be used. System Si 136 will be used as an example in the demonstration, hopefully it will serve as a no- brain procedure in determination of irreducible representation and optically active modes without going through too much group theory. The first step is to gather the structure information for the system under consideration. Here we know that Si 136 belongs to space group 227, and silicon atoms are located at three different Wyckoff sites, 8a, 32e, and 96g. Note this is the notation for conventional cubic cell. In primitive cell, they can be regarded as 2a, 8e, and 24g sites respectively. We shall consider each site separately. For system with multiple types of atoms, each site of every type of atoms should be considered separately. However the idea is the same for every one of them. Then we use the Site Symmetry Table to find the point symmetry operations of each Wyckoff site. The symmetry table lists the space group, point group, and also the point symmetry operations in S notation for all the Wyckoff positions the space group has and they are listed in the alphabetical order starting from a. Site symmetry table gives the following information for space group 227: I W. G. Fateley, F. R. Dollish, N. T. McDevitt and F. F. Bentley, Infrared and Raman Selection Rules for Molecular and Lattice Vibrations: The correlation Method (John Wiley & Sons, Inc. 1972). 151 Table B.I Site symmetries of space group 227 Space group Site symmetries 227 Fd3m O h 7 2T d (2); 2D 3d (4); C 3v (8); C 2v (12); C s (24); C 2 (24); C 1 (48) T d , D 3d , C 3v , C 2v , C s , C 2 and C 1 are subgroups of point group O h . Since the sites symmetries are ordered alphabetically (see Table B.II), it is rather easy to find the site symmetry for any Wyckoff site we desire. The site symmetry of 8a, 32e and 96g are T d , C 3v and C s respectively. Let?s consider 8a (T d ) site first. Table B.II Wyckoff letter representation of each site symmetry of space group 227 Wyckoff Position 8a 8b 16c 16d 32e 48f 96g 96h 96i Site symmetry T d T d D 3d D 3d C 3v C 2v C s C 2 C 1 We then go to Character Table of site symmetry T d and identify all the symmetry representations of the translations T x ,T y and T z . From Ref. I appendix II, we can find the character table for site symmetry T d . Table B.III list all the symmetry operations (E, C 3 , C 2 , S 4 , ? d ) of the point group T d and also its subsets (A 1 , A 2 , E, T 1 , T 2 ), which are often called symmetry representations and they are orthogonal to one another. Note that representations like A or B has single degeneracy, E has double degeneracy and T has triple degeneracy. The numbers in Table B.III are called character, indicating the effect of an operation in a given representation. The righteous two columns are function operators. 152 We can see that T 2 contains translational operators, similarly, A 1 +E for C 3v and 2A?+A?? for C s . Table B.III Character table for T d point group. T d E 8C 3 3C 2 6S 4 6? d Linear, Roations Quadratic A 1 1 1 1 1 1 X 2 +y 2 +z 2 A 1 1 1 1 -1 -1 E 2 -1 2 0 0 (2z 2 -x 2 -y 2 ,x 2 -y 2 ) T 1 3 0 -1 1 -1 (R x ,R y ,R z ) T 2 3 0 -1 -1 1 (T x ,T y , T z ) (xy,xz,yz) Table B.IV Correlation table for O h point group (only T d , C 3v and C s are listed) O h T d C 3v C s (? h ) C s (? d ) A 1g A 1 A 1 A? A? A 2g A 2 A 2 A? A?? E g E E 2A? A?+A?? T 1g T 1 A 2 +E A?+A?? A?+2A?? T 2g T 2 A 1 +E A?+2A?? 2A?+A?? A 1u A 2 A 2 A?? A?? A 2u A 1 A 1 A?? A? E u E E 2A?? A?+A?? T 1u T 2 A 1 +E 2A?+A?? 2A?+A?? T 2u T 1 A 2 +E 2A?+A?? A?+2A?? 153 Character table only gives the irreducible representations for point group T d . However we want to have the irreducible representations for point group O h which has higher symmetry, thus we need to transform them into those in O h representations. Such transformation is listed in Correlation Table. Ref. I appendix III lists the correlation tables for all the point groups. Here we cited the correlation between O h and its subgroups we needed for Si 136 in Table B.IV: Notice that in the Table B.IV, we found two columns for C s , one tagged with ? h and the other with ? d . Which one shall we use? Table 14 in Ref. I tells the answer. In Tables 14, we find the following for space group 227: Space group number C 2 ?, ? h C 2 , ? d ? d 227 O h 7 h f g Table 14 says that if the atom under consideration is in Wyckoff site h, the column with (C 2 ?, ? h ) should be chosen; if it is f, then the column with (C 2 , ? d ) should be chosen; if it is g, the column with (? d ) should be chosen. In current case, the Wyckoff site is g, thus we shall use the column tagged with ? d . The correspondence between the representations of O h and those of its subgroup T d , C 3v and C s is summarized in Table B.V. Therefore the irreducible representation for gamma phonon modes is: 12 1212 12 31 458 3 485 g ggg gu uuuu A AETTAAETT?= + + + + + + + + + (B.I) In crystal, there are 3N degrees of vibrational freedom, where N is the number of atoms in the Bravais cell. For Si 136 , we used the 34-atom primitive cell. We are expecting 154 3N=102 degrees of vibrational freedom. This can be checked at this point by adding all the contributions in the representation. 3111 4 2 53 8 311 31 4 2 8 3 53 102?=?+?+?+?+?+?+?+?+?+?= Table B.V Summary of symmetry representations for Si 136 Point group Symmetry representation Symmetry representation of point group Oh T d T 2 T 2g +T 1u A 1 A 1g +T 2g +A 2u +T 1u C 3v E E g +T 1g +T 2g +E u +T 1u +T 2u 2A? 2(A 1g +E g +T 1g +2T 2g +A 2u +E u +2T 1u +T 2u ) C s A?? A 2g +E g +2T 1g +T 2g +A 1u + Eu +T 1u +2T 2u Among 3N degrees of vibrational freedom, three are acoustical vibrations. And when we consider only vibrations at the Brillouin zone center, their frequencies are zero, and they are of no physical interest, these three acoustical vibrations should be subtracted from the irreducible representations. The symmetry representations of acoustical vibrations are those correspond to translations T x , T y and T z in point group character table, and they are also IR active modes. Raman active modes are those corresponding to rotation R x , R y , and R z . O h character table tells that acoustic modes are T 1u mode, and it is also the IR active mode; Raman active modes are T 1g mode. Subtract one T 1u mode from Equ. (B.I), we then get 155 the optical irreducible representation for Si 136 in Equ. (B.II). IR and Raman active modes are indicated by superscript. 12 1 212 1 2 31 45 8 3 47 5 optic Raman IR g ggg gu uuu u A AET TAAET T?= + + + + ++ + + + (B.II)