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 firstprinciples total energy and force calculations were extended to the
calculations of harmonic force constant matrices and third order lattice anharmonicity
tensors with an efficient supercell finitedifference 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 firstprinciple 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 quasiharmonic approximation. Kinetic theory was
implemented to predict the nonequilibrium 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. GreenKubo 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. AnBan 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 parentsinlaw, my husband and my little baby Kevin for their
continuous love and support.
Last but not least, I would like to thank NSF (EPS0447675 and HRD0317741)
and DOE (DEFG0203ER46060) 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 FIRSTPRINCIPLES ATOMSCALE SIMULATION AND MODELING
OF CRYSTALLINE SOLIDS ............................................................................................ 8
2.1 FirstPrinciples Density Functional Theory (DFT)............................................... 8
2.2 AtomScale 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 QuasiHarmonic Approximation and Mode Gr?neisen Parameters ........ 39
3.1.3 Anharmonicity Induced PhononPhonon 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 NonEquilibirum Thermal Transport Theories................................................... 46
3.3.1 Kinetic Transport Theory......................................................................... 46
3.3.2 GreenKubo 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 TP Phase Relations ......................................... 69
4.5 Thermal Properties.............................................................................................. 71
4.6. Conclusions........................................................................................................ 74
CHAPTER 5 THERMAL CONDUCTIVITY OF TETRHEDRALLY BONDED
SILICON CRYSTALS ? THE GREENKUBO APPROACH ........................................ 76
5.1 Introduction......................................................................................................... 76
5.2 Empirical SiSi 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 unenforced (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 dSi 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) 4fold ?folded? phonon dispersion plots of dSi, calculated with unitcell
size four times that of the primitive cell ( 4
cell primitive
aa= ?=2.16 nm), along with (b) the
replotted phonon dispersion of Si
136
(
clathrate
a = 1.46 nm)................................................ 60
Figure 4.4 The theoretical mode Gr?neisen parameters ?
i
along three highsymmetry
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 cm1 that corresponds mainly to a trace of dSi impurity in the sample. However,
there is also a calculated mode at 516 cm1 for Si
136
clathrate at this position.. .............. 45
Figure 4.6 Selected Raman spectra of Si
136
collected at compression. Asterisk marks
diamondphase silicon.. .................................................................................................... 46
Figure 4.7 Pressure dependence of Raman shifts of Si
136
and dSi. ................................ 47
Figure 4.8 The LDA predicted vibrational entropies in dSi (solid line) and Si
136
(dashed
line). .................................................................................................................................. 69
Figure 4.9 The theoretically predicted equilibrium TP 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 dSi and
Si
136
. .................................................................................................................................. 71
Figure 4.11 The theoretically predicted linear coefficients of lattice thermal expansion at
P=0 in dSi and Si
136
. ........................................................................................................ 73
Figure 4.12 Highresolution Xray 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 clathratestructured Si
136
obtained by fitting to high angle powder Xray
diffraction data. The line drawn through the data points is a guide to the eye. ................ 74
Figure 5.1 Phonon dispersion of dSi 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 dSi system .................... 86
Figure 5.4 The normalized time correlation function ()g t at room temperature for dSi.
(a) an overall look; (b) closeup look at the beginning of the MD run; (c) closeup 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 timedependent correlation function
()g t
at room temperature
for Si
136
. (a) an overall look; (b) closeup 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 dSi. ................ 91
Figure 5.8 Power Spectrum of heat current autocorrelation function for Si
136
................. 92
Figure 5.9 Normalized the heat current autocorrelation of dSi before (black line) and
after (red line) the lowpass 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 lowpass filter with filtering frequency 0.5Hz.................................... 94
Figure 5.11 log plot of the normalized lowpass (0.5Hz) filtered heat current
autocorrelation of dSi (black line) and Si
136
(red line). ................................................... 95
xiv
Figure 6.1 Fitted effective phonon life time
eff
?
of dSi to its experimental heat
conductivity. Temperature dependence of Si
136
was estimated assuming the equal
eff
?
... 99
Figure 6.2 Atomic model of dSi .................................................................................... 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 bondstretching
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 dSi.......................... 106
Figure 6.7 Pressure dependence of 8a32e bond stretching A
ijk
terms in Si
136
. ............. 106
Figure 6.8 Pressure dependence of 32e96g bond stretching Aijk terms in Si
136
........... 107
Figure 6.9 Pressure dependence of 96g96g type (a) bond stretching A
ijk
terms in Si
136
.
......................................................................................................................................... 107
Figure 6.10 Pressure dependence of 96g96g type (b) bond stretching A
ijk
terms in Si
136
.
......................................................................................................................................... 108
Figure 6.11 Mode Gr?neisen parameters of dSi 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 FirstPrinciples method with a 444? ? qpoint 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 Onetoone correspondence of the index for the displacement to the index of
atom being displaced......................................................................................................... 26
Table 2.2 Onetoone correspondence of index for the pair displacement to that for the
single atom displacement.................................................................................................. 32
Table 4.1 LDA calculated static (T=0K) BirchMurnaghan equation of state of the ground
state diamond phase Si (dSi) and the metastable typeII 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 dSi and Ref. [65] for Si
136
. ........................................................ 58
Table 4.2 Comparison of acoustic velocities (m/s) for dSi (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 clathratestructured 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 dSi 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 supercell of dSi........................................................ 101
Table 6.2 List of coordinates of all 4 first nearest neighbors of one representative atom
for each Wyckoff site in our chosen supercell 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 dSi and Si
136
....................... 110
Table 7.1 Our LDA calculated static equilibrium properties of MgO fitted by 3
rd
order
BirchMurnaghan 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
Atomscale numerical simulation techniques, especially those based on the first
principles electronicstructure theories, have become powerful tools to understand
structures, formations, dynamics, and many other physical and chemical properties of real
and complex materials systems. Firstprinciples 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,
firstprinciples 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 firstprinciples 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 lowcost 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 firstprinciples density functional theory (DFT). We have
successfully developed a realspace supercell 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
principlesmethod 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 quasiharmonic approximation (QHA), and (2) lattice
thermal conductivity based on the kinetic transport theory. An alternative simulation
3
approach based on the GreenKubo theory is also implemented to evaluate the lattice
thermal conductivity as comparison.
Results of firstprinciples calculation of two material systems are presented in
this dissertation. The first is nanoopen crystalline typeII silicon clathrate (Si
136
)
materials. This new Si allotrope has an opencage structure, isostructural with low
density inclusion compounds of H
2
Oice. It has a cubic framework in which each cubic
unitcell contains sixteen 20atom ?cages? (dodecahedra) and eight 28atom ?cages?
(hexakaidecahedra). In addition to the elemental ?guestfree? 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
?glasslike? thermal conductivity due to scattering of acoustic heatcarrying 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 (dGe).
2
In 2000, J. Dong et al.
3
performed a
theoretical calculation for both guestfree and guestencapsulated 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 guestfree 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 guestfree 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 expandedframework semiconducting crystals,
including studies of both fundamental thermal properties and thermal conductivity. As a
first step in this area, we use firstprinciples theoretical methods to predict the measurable
thermal properties (such as heat capacity and thermal expansion) of the guestfree
clathrate Si
136
. Our results are analyzed and then compared with previous data on the
wellknown ground state diamondstructured phase of this element (dSi)
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 clathratestructured polymorph
Si
136
should exhibit a region of negative thermal expansion below 124K, like the
diamondstructured phase. This prediction has been confirmed by an experimental
measurement. Then we calculate the timecorrelation functions of heat current in both d
Si and Si
136
using classical molecular dynamics (MD) with an empirical Tersoff potential
at equilibrium microcanonical (N,E,V) conditions, and derive the nonequilibrium
thermal transport properties based on statistical fluctuationdissipation theory
7
(the
GreenKubo formula). The calculation indicates that thermal conductivity of Si
136
is only
about 10.6% of that of dSi. 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 endmember 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 firstprinciples 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 finitesize 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 firstprinciples method
16
. In
addition to the conventional quasiharmonic 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 firstprinciples 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? ? qpoint 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 firstprinciples 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: dSi and Si
136
. Chapter 6 studies the
3
rd
order lattice anharmonicity and 1
st
order approximation of phononphonon scattering
rates using the same empirical (Tersoff) potentials and reexams the ratio of lattice
thermal conductivity between Si
136
and dSi based on the kinetic transport theory.
7
Chapter 7 presents the firstprinciples 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
FIRSTPRINCIPLES ATOMSCALE SIMULATION AND MODELING OF
CRYSTALLINE SOLIDS
2.1 FirstPrinciples Density Functional Theory (DFT)
An accurate total energy theory, which predicts the energy of an Natom
material system at a given structural configuration, is the foundation for all atomscale
simulation and modeling. Because of the large mass ratios between nuclei and electrons
inside an atom and the fact that electronic structures of coreelectrons are not
environmentsensitive, an atom inside a solid are often considered as a positivelycharged
ion (including nuclei and core electrons) surrounded by a group of negatively charged
valence electrons, and ions and electrons move at different timescales. 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 manybody 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 solidstate materials
systems. Firstprinciples electronicstructure 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 exchangecorrelation 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 singleelectron approach to exactly calculate the total
energy of a manyelectron 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
exchangecorrelation interactions of many electrons. Unless specified otherwise, results
reported in this dissertation are obtained with the LDA. To selfconsistently solve the KS
equation, the electronic wavefunctions 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 jth eigenfunction of the kpoint 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 wavevector. 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 (bondforming)
valence electrons are solved explicitly, and the effects of (environmentinsensitive) core
electrons are approximated with the socalled pseudopotentials (PP). Two widely
adopted pseudopotential types are norm conserving pseudo potentials
22
and ultrasoft
pseudopotentials
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 firstorder
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 FeynmanHellmann (FH) force, and the remaining
12
two terms are called Pulay forces. The FeynmanHellmann theorem
24
states that Pulay
forces vanish when the calculation has reached selfconsistency and the basis set
orthornomality persists and are independent of the atomic positions. Therefore, in the
planewave implemented pseudo potential calculation, the FH 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 AtomScale 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 interatomic 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
threedimensional (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 unitcell, 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 unitcell 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 unitcell 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. basecentered (bac), facecentered (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 facecenter arrangement belongs to the face
centeredcubic (fcc) Bravais lattice type, and its primitive unitcell 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, dSi, Si
136
,
and B1MgO, 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 unitcell. Si
136
crystals adopt the typeII clathrate lattice with 34 atoms per fcc unit
cell (or 136 atoms per cubic cell),
14
The arrangement of atoms in the unitcell can be uniquely specified with the
locations of atoms inside a unitcell 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 unitcell of dSi locate at (0, 0, 0) and (14,14,14) sites, while the one Mg
atom and one O atom in the B1structured 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 ith 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 pointgroup (PG) (rotation) operation and/or a spacegroup (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 dSi and B1MgO 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 pointgroup symmetry operations with possible glide plane and screw axis. We
15
also constructed an atomatom 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 atomatom 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
atomatom 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 Xray 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 wavevector 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 (kspace), of a crystal.
The reciprocal lattice vectors of a crystal are defined based on the corresponding real
space unitcell 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 bodycentered 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 interplane 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 latticeperiodic 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 WignerSeitz
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=?? unitcells (
123
,,NNN are the numbers of unit
cells along lattice vector direction
123
,,aaa
v vv
respectively) satisfying the Bornvon Karman
18
periodic boundary condition, the number of kpoints 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 quasicontinuous.
The evaluation of many physical properties involves integrating the periodic
functions of kvector in the BZ, therefore, the decent dense kpoint sampling is crucial in
order to yield accurate results, usually convergence test is desired in real calculation.
Among all the uniform kpoints sampled in BZ, many of them are dependent on the
others due to the lattice point symmetry. All the independent kpoints 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 MonkhorstPack
28
kpoint
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 kpoint set and then assign the calculated physical quantities to all the other
equivalent kpoints according to symmetry map. If the quantity is a scalar, we can simply
assign the same value at an independent kpoint to all the other equivalent kpoints.
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 kpoints 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 kpoint 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, unitcell models are adopted to describe lattices and atomic basis of
crystals. Based on its unitcell 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 spacegroup symmetry
operators and all the symmetryallowed 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 pseudopotential
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 pseudopotential of
nuclei; And KPOINTS contains the kpoint 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 energyvolume 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 BirchMurnaghan 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 interatomic 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 higherorder 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 unitcell 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 qpoint (kpoint is conventionally used for electronic energy calculation,
and qpoint 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 eigenfrequency (?) and eigenvector e of the ?
th
vibration eigenmode at a
reciprocal space qpoint 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 3n3
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
lowfrequency acoustic vibration correspond to the sound waves in the lattice.
We have implemented and further developed a realspace supercell technique to
calculate the forceconstant matrix ? with firstprinciples 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 zeroforce, 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 Onetoone correspondence of the index for the displacement to the index of
atom being displaced.
1 2 3 4 5 6 L j L 6N5 6N4 6N3 6N2 6N1 6N
atom #1
atom #N
27
The onetoone correspondence between the index of moves for singleatom
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 atomatom 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 FH 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
atomatom mapping under the current Operator , and _f force can be calculated by:
__fforce Operator i force= ? . (2.25)
Since the atomatom 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 FH 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 cancelout effect.
We have further developed a paireddisplacement 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 FH
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 FH 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 Onetoone 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
imove
jmove
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 imove from 1 to 6N, and another loop over
jmove from imove to 6N, note that only those equal to or larger than imove 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 mmove for imove and nmove for jmove 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 (iatom and j
atom) and direction (ibasis and jbasis) for pair move (, )i j ; then find correspondent pair
atom (matom and natom) for (iatom and jatom), and displacement direction (mbasis
and nbasis) for (ibasis and jbasis), and thus get the index move mmove and nmove
(,)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 HF 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 FH 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 dSi 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 qpoints 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 unenforced (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 eigenvector 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 eigenvector 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 FeynmanHellmann theorem :
39
1
(,) (,) () (,)
2
gq
v qi eqi q eqi
?
=?Dv
vv vv v v
(3.7)
3.1.2 QuasiHarmonic 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 quasiharmonic approximation
(QHA), where the harmonic oscillator model is valid, yet the oscillation frequencies
become volumedependent. 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 FeynmanHellmann 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 eigenvectors 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 unitcell,
,q ?
?v and ()
v
qe
r v
are eigenfrequency and eigenvector
of ?th eigenmode 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 phononphonon scattering rate in phonon life time
and thermal conductivity calculation.
3.1.3 Anharmonicity Induced PhononPhonon Scatterings
Any lattice imperfection (such as finitesize 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 threephonon scattering mechanism. Four
phonon scattering may be important at high temperature, which will be explored in future
work.
41
In general, threephonon 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 3phonon
processes are shown in Figure 3.1.
Figure 3.1 Schematic figures for two types of 3phonon 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 quantummechanical timedependent 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 3phonon 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 BoseEinstain 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 firstprinciples total energy theory discussed in previous section.
3.2.2 Macroscopic Thermal Properties within the QHA
To connect the firstprinciples calculated thermodynamic potentials with
experimental measurements, we can further evaluate the thermal properties based on the
macroscopic thermodynamic theories. For example, to predict the isothermal 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 rewrite 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 DulongPetit 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 NonEquilibirum 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 nonequilibrium 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 fullwidth at halfmaximum (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 wellknown BoseEinstein 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 3phonon 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 qvector 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 3phonon 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 3phonon
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 SingleMode 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 GreenKubo Formula
Nonequilibirum transport properties are also associated with the fluctuation
phenomena at thermal equilibirum. According to the FluctuationDissipation Theorem
(FDT) of GreenKubo
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 steadystate flow at a nonequilibrium 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 timecorrelation 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 timedependent heat current is obtained by performing
molecular dynamics (MD) simulations. Empirical Tersoff potential has been used in our
study of silicon materials.
The WienerKhintchine 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 timedependent 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,
diamondstructured silicon (dSi). Several metastable forms, including amorphous (a
Si)
41
and nanostructured (nanoSi) silicon
42
, also have a variety of industrial applications.
Highdensity polymorphic forms of Si, which are synthesized under highpressure
conditions, have been explored in great detail both theoretically and experimentally.
43
A
common feature of the observed highpressure phases is the increase of the coordination
numbers of Si atoms from four (as in the ground state dSi and metastable aSi 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 lowdensity
crystalline form of elemental Si was synthesized by removal of Na atoms from the
clathratestructured framework compounds Na
x
Si
136
44,45,4
prepared by controlled thermal
decomposition of NaSi. In contrast to the dense highpressure phases, the atoms in this
55
new low density Si material are fully tetrahedrally bonded, and it is a widergap
semiconductor (E
g
= 1.9 eV)
4
. The new Si allotrope has the typeII clathrate structure,
isostructural with lowdensity inclusion compounds of H
2
Oice, that has a cubic
framework in which each cubic unitcell contains sixteen 20atom ?cages? (dodecahedra)
and eight 28atom ?cages? (hexakaidecahedra) (For example, Figure 1 in Ref. 44). In
addition to the elemental ?guestfree? 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 typeI
clathrate structure. Although various series of guestencapsulated binary or ternary
compounds involving Si framework atoms have been prepared in many laboratories, no
one has yet reported a guestfree elemental solid with the typeI 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 guestencapsulated 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 widegap
semiconducting
4
to metallic
46
and even superconducting (T
c
=68 K) in (Na,Ba)
containing Si clathrates.
47,48
In recent years, the potential applications of Si, Ge, and Sn
based clathratestructured 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 ?glasslike? thermal conductivity
due to scattering of acoustic heatcarrying 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 ?phononglass electroncrystal? paradigm (PGEC).
53,54
This model then raised the possibility of ?tuning? the hostguest 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 guestfree Ge
46
clathrates could be lowered by at least one order of
magnitude compared with that of the corresponding diamondstructured Ge crystals.
3
Recently, Nolas et al.
52
demonstrated experimentally that the guestfree 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 clathratestructured material, or the flattening in the phonon
dispersion relations associated with formation of the lowdensity structure, that is
correlated with the presence of fivemembered 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 expandedframework semiconducting crystals. As a first
step in this area, we used firstprinciples theoretical methods to predict the measurable
thermal properties (such as heat capacity and thermal expansion) of the guestfree
57
clathrate Si
136
. Our results are analyzed and compared with previous data on the well
known ground state diamondstructured 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 clathratestructured 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 facecenteredcubic (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
typeII 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 dSi and Si
136
crystals are
fitted to a 3
rd
order BirchMurnaghan equation of state (BMEoS), 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 dSi.
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 metastable phases. This result is
consistent with the observed metastability and recovery to ambient conditions of the
guestfree clathrate phase Si
136
, synthesized experimentally.
4
The calculated (static)
equilibrium cubic lattice constants are 0.5397 nm and 1.455 nm for the 8atom unitcell
of dSi and 136atom unitcell of Si
136
respectively. These are in excellent agreement
with experimental measurements of 0.5431 nm and 1.4644 nm for dSi
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 dSi. The mass density of the guestfree
clathrate is 13% less than that of dSi. 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) BirchMurnaghan equation of state of the
ground state diamond phase Si (dSi) and the metastable typeII 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 dSi and Ref. [65] for Si
136
.
dSi
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 dSi 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) 4fold ?folded? phonon dispersion plots of dSi, calculated with unitcell
size four times that of the primitive cell ( 4
cell primitive
aa= ?=2.16 nm), along with (b) the
replotted 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 tetrahedrallycoordinated Si atoms with similar sp
3
SiSi
bonding, we expect the local force constants and vibrational properties to have analogous
behavior. Three regions within the calculated ?
i
( q
r
) relations in dSi (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 dSi (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 dSi in Figure 4.3(a) using a unitcell of four times of the primitive unitcell,
along with the replotted 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 dSi, with the dispersion relations ?folded? back towards the first
Brillouin zone of dSi. This observation is important, because it means that the elastic
properties and phonon propagation relations within dSi and Si
136
at small wavevectors
are similar to each other. To further illustrate this similarity, we listed the calculated
group velocities of acoustic phonon modes in dSi and Si
136
in Table 4.2.
Table 4.2 Comparison of acoustic velocities (m/s) for dSi (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
dSi
EXPT.
dSi
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 dSi 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 dSi: especially in the lowfrequency regime between 150200 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 dSi
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 5membered rings; these cause the expansion of the unit
cell motif and the decrease in density, compared with dSi. 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 dSi, 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 dSi. 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 firstprinciples methods. Also, the nearly
dispersionless phonon branches occurring near 200 cm
1
predicted by the firstprinciples
calculations were more ?spread out? in frequency, in the empirical study.
Figure 4.4 The theoretical mode Gr?neisen parameters ?
i
along three highsymmetry
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 dSi 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 dSi, 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 dSi 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 dSi (Figure 4.1(c)). Although the boundaries between the three
regions are less well defined for Si
136
compared with dSi, 99% of the optic phonon
branches for Si
136
can be characterized as derived from TAlike, LAlike, or optic
branches issued from dSi. The values of the Gr?neisen parameters
i
? for Si
136
for these
groups of modes are similar to those for dSi; 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 dSi.
The calculated zonecenter 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 highfrequency modes that
are associated with SiSi 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 cm1 that corresponds mainly to a trace of dSi impurity in the sample. However,
there is also a calculated mode at 516 cm1 for Si
136
clathrate at this position.
66
Figure 4.6 Selected Raman spectra of Si
136
collected at compression. Asterisk marks
diamondphase silicon.
There are some discrepancies between the theoretical and experimentally
observed pressure shifts, mainly occurring in the midfrequency 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 SiSi
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 dSi.
68
Table 4.3 LDA calculated and experimentally measured Raman frequencies (?) and
mode Gr?neisen parameters (?
i
) for clathratestructured 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
2phonon 920
2phonon 965
69
4.4 Thermodynamic Potentials and TP 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 BMEOS. 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 dSi 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 dSi; 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 dSi (solid line) and Si
136
(dashed
line).
70
This small entropy difference (?S) between the Si
136
and dSi 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
dSi 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 dSi 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 TP 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
clathratestructured
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 dSi.
52
Our quasiharmonic computational results are
clearly valid for dSi; 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 dSi 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 dSi and
Si
136
.
72
We then used our LDA results to predict the ?(T) relations at P=0. It is well
known that dSi exhibits a region of negative thermal expansion (NTE) between T =0
150 K, that is often ascribed to anharmonic phononphonon interactions, but that is well
reproduced by quasiharmonic 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 quasiharmonic model. The predicted ?(T) behavior at low T for
dSi and Si
136
calculated at the same level of theory is shown in Figure 4.11. The result
for dSi 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 0125K. Low temperature Xray 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 Naloss occurring during formation of the guestfree clathrate.
However, the data points collected are sufficient to show at least a flattening in the ?(T)
data between 5200 K, with a minimum apparent near 90120 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 dSi and Si
136
.
Figure 4.12 Highresolution Xray 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 clathratestructured Si
136
obtained by fitting to high angle powder Xray
diffraction data. The line drawn through the data points is a guide to the eye.
The origin of the lowT NTE behavior in dSi 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 guestfree clathrate
polymorph of silicon (Si
136
) based on firstprinciples calculations combined with
experimental Xray and Raman scattering measurements. The Si
136
clathrate is metastable
75
compared with the dSi 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 diamondstructured
and ?expandedvolume? clathrate polymorphs, our current studies reveal that the thermal
properties of the two phases involving longwavelength phonons are similar to each other.
The vibrational properties of Si
136
phonons are similar to those of dSi, and they can be
understood in terms of Brillouin zone reduction following the unit cell expansion
between dSi 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 GREENKUBO 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 solidstate cooling for microelectronic devices. The key
limiting factor of current thermoelectric technologies is the lack of highperformance
thermoelectric materials, whose efficiency is described with a dimensionless 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 ?glasslike? low thermal conductivity due
to the scattering of acoustic heatcarrying 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 ?phononglasselectroncrystal? paradigm
(PGEC).
53,54
This model then raised the possibility of ?tuning? the hostguest 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 typeI 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 twoorder of
magnitude lower than that of dGe.
In addition to the importance to thermoelectric applications, such glasslike 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 GreenKubo theory
and molecular dynamics simulations, they demonstrated that the lattice thermal
conductivity (? ) of guestfree Ge
46
clathrates could be lowered by about one order of
78
magnitude compared with that of the corresponding diamondstructured 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 guestfree Ge clathrates remains unverified since no laboratories have
successfully removed encapsulated guest atoms from the Gebased clathrates. The only
pristine clathrate system that has been synthesized is the typeII silicon clathrate (Si
136
).
Si
136
, which is very similar to dSi in the aspects of local SiSi 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 guestfree 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 guestfree 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 dSi 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 Gebased typeI clathrate. The current study directly examines the typeII 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 dSi 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 SiSi 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 BirchMurnaghan equations of state of the
two Si phases, the diamondstructure Si (dSi) and typeII 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 dSi
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 dSi, and
consequently its calculated V
0
and B
0
are nearly identical to those of room temperature
experimental data of dSi. 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 dSi and Si
136
together with
experimental data.
dSi 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 dSi and Si
136
from this empirical potential
calculation are shown in Figure 5.1 and 5.2 respectively.
Figure 5.1 Phonon dispersion of dSi 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 dSi calculated with the Si Tersoff potential are in
good agreement with those obtained with the firstprinciples DFT method.
5.3 Results and Discussions
In the current study, a 2744atom super cell, i.e. 7 7 7? ? of the 8atom cubic unit
cell, and a 3672atom super cell (a 333? ? of the 136atom cubic cell) have been chosen
dSi and typeII Si
136
clathrate respectively. In the MD simulation, we chose the single
simulation timestep 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 microcanonical (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 dSi 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 timeaveraged correlation function G(t) is timeconsuming.
Instead, we first do Fouriertransformation for the heat current function J(t), then
86
calculate the powerspectrum 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 timedecay 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 dSi 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
dSi 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 dSi is 34.45 W/K/M/ps.
Figure 5.4 shows our calculated normalized heat current autocorrelation in dSi 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 dSi.
(a) an overall look; (b) closeup look at the beginning of the MD run; (c) closeup 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 dSi. Based
on this definition, we obtained the calculated lattice thermal conductivity of dSi 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 dSi 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 underestimation of lattice anharmonicity by the
adopted Si Tersoff potential is likely to be the same in both dSi 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 timedependent correlation function
()g t
at room temperature
for Si
136
. (a) an overall look; (b) closeup 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 dSi. 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 dSi. On the other hand, the
correlation function of the heat current in Si
136
is dramatically different from that in dSi.
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 dSi. 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 timeintegration 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 dSi. As the G(0) are comparable in dSi and Si
136
, our GreenKubo
calculations predicted that the lattice thermal conductivity of Si
136
(? at 300K ~ 51
W
mK
) is only about 10.6% of that dSi (482
W
mK
). Our predicted oneorder of
magnitude reduction in typeII Si clathrates is comparable to the previous findings in
91
typeI Ge clathrates, while the experiment measurement of Nolas et al.
5
reported a much
larger twoorder 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 dSi 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 highfrequency components in dSi are negligibly
small, while there are several strong highfrequency 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 highfrequency 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 dSi.
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 nonzero ?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 nonzero for many
vibration modes, the only nonvanishing 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
phononphonon scattering. To quantitatively estimate the contribution from each term,
we performed the lowpass filtering on the g(t) plot to eliminate the highfrequency
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 dSi 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 dSi which contains similar tetrahedral SiSi
bonds, which is associated with the oscillation feature in Si clathrates.
94
Figure 5.9 Normalized the heat current autocorrelation of dSi before (black line) and
after (red line) the lowpass 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 lowpass 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 dSi 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 lowpass (0.5Hz) filtered heat current
autocorrelation of dSi (black line) and Si
136
(red line).
5.4 Conclusion
We have carried out a series of large supercell MD simulations of dSi 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 GreenKubo 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 dSi.
Our predicted 89.4% reduction of ? in typeII Si clathrate is close to the predicted
reduction in typeI 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 dSi and Si
136
using
GreenKubo formula with empirical Tersoff potential. The calculation revealed that the
thermal conductivity of Si
136
is only about 10.6% of that of dSi 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 clathratestructured material, or the
flattening in the phonon dispersion relations associated with formation of the lowdensity
structure, which is correlated with the presence of fivemembered 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. 128atom dSi 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 dSi 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 dSi 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 lifetime of each phonon mode. As the firstorder 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 dSi. We have found the
temperature dependence of
eff
? (blueplusdashed 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 dSi have the same effective phonon life time:
136Si dSi
eff eff
??= ,
we have calculated heat conductivity for Si
136
which is shown as the reddiamonddashed
square in Figure 6.1. Mode heat capacity and group velocity for both systems were
calculated with kpoint grids 16 16 16? ? and 888? ? for dSi and Si
136
respectively. The
experimental data are based on natural Si crystals, which include isotope disorder. At
T=300K, the isotopeenriched dSi 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 dSi 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 dSi. Figure 6.1 showed
that thermal conductivity of Si
136
is much lower than that of dSi 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 GreenKubo 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 xdirection, ydirection, or zdirection 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 dSi 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 supercell of dSi.
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 dSi.
102
Table 6.2 List of coordinates of all 4 first nearest neighbors of one representative atom
for each Wyckoff site in our chosen supercell 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 dSi. 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 dSi 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 xy 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 bondstretching terms as long as they both are
perpendicular to z?.
Figure 6.4 A schematic show of coordinate transformation to yield a bondstretching 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 dSi (8a8a), while there are four types of bonds in Si
136
, and they are 8a32e,
32e96g, 96g96g(a), 96g96g(b). Pressure dependence of bondstretching
ijk
A for all
these 5 types of bond has been investigated and they are shown in Figure 6.6 for dSi and
Figure 6.76.10 for Si
136
respectively.
106
Figure 6.6 Pressure dependence of bond stretching A
ijk
terms in dSi.
Figure 6.7 Pressure dependence of 8a32e bond stretching A
ijk
terms in Si
136
.
1
65
y
107
Figure 6.8 Pressure dependence of 32e96g bond stretching A
ijk
terms in Si
136
.
Figure 6.9 Pressure dependence of 96g96g type (a) bond stretching A
ijk
terms in Si
136
.
108
Figure 6.10 Pressure dependence of 96g96g type (b) bond stretching A
ijk
terms in Si
136
.
As one can see from Figure 6.76.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 dSi with that of Si
136
, we
first calculated the averaged bondstretching 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 dSi 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 dSi and Si
136
System B(GPa)
333
A ( eV/?
3
)
333
/AB
336
A ( eV/?
3
)
336
/AB
dSi 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 dSi 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 dSi 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 dSi. 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 dSi. Section 6.3 calculated the lattice anharmonicity and
analyzed the bondstretching 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 dSi. And note that this ratio is a ballpark 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 finitesize grain boundaries). Current understanding on lattice
anharmonicity and its pressure dependence is limited. Using the statistical linear response
theory (GreenKubo 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 phononphonon
scattering at compression may behave differently at low and high pressure ranges. The
interatomic potential adopted in that study was a nonempirical 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 interatomic interactions,
(nonequilibrium) 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 phononphonon scattering is explicitly
evaluated in the GreenKubo 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 selfconsistent planewave methods.
Recently, Oganov and Dorogokupets reported a study about the anharmonicity effects on
the thermodynamic potentials of MgO using a firstprinciples method
99
. In addition to
the conventional quasiharmonic 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 firstprinciples calculation of harmonic phonon spectra and
3
rd
order lattice anharmonicity in MgO. Explicit calculation of both anharmonicity and
massdisorder 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 endmember 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 finitedisplacement (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 phononphonon
scattering rates, and such rates are needed for evaluating phonon life time and lattice
thermal conductivity based on nonequilibrium 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? ? qpoint sampling over the BZ. In addition to the
anharmonicity induced phononphonon scattering, we also considered the mass disorder
116
(isotope effect) induced phononphonon 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 firstprinciples 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
Twoatom facecenteredcubic (fcc) unitcell 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? ? MonkhorstPack 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
semicore 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
BirchMurnaghan 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
BirchMurnaghan 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 128atom super cell model with only ?point sampling in the Brillouin zone
was used in the calculation for HF 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 highpressure 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 PTV
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 singleatom 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 bondingstretching
anharmonicity between two neighboring MgO 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 dSi
123
and Si
136
, but this is based on the FH force calculated from firstprinciples, 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
Phononscattering 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 eigenfrequency 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?104 7.398?104
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 phononphonon 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??qgrids.
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)
?
isotopepure
(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 isotopepure MgO (
isotope pure
?
?
). At 500K, isotopepure 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 phononphonon 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 FirstPrinciples method with a 444? ? qpoint 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
anharmonicityinduced phonon scattering with just one single parameter: the effective
phonon lifetime (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 supercell finitedisplacement algorithm, we
studied the harmonic and anharmonic lattice dynamics in MgO at highpressure. 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 phononphonon 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, isotopepure
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 qgrid 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 firstprinciples 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 firstprinciple 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 supercell finitedisplacement
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 SCFD algorithm is currently implemented
with the VASP code, but it can also be easily implemented with any other firstprinciple
codes with any forms of pseudopotentials and electron wave basis sets. This algorithm is
ideal for parallel computation using lowcost 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 firstprinciples 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 quasiharmonic
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 nonequilibrium 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 guestfree clathrate polymorph of silicon (Si
136
) based on first
principles calculations combined with experimental Xray and Raman scattering
measurements. The Si
136
clathrate is metastable compared with the dSi 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 diamondstructured and ?expandedvolume? 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 dSi, and they can be understood in terms of Brillouin zone
reduction following the unit cell expansion between dSi 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 GreenKubo formulism to calculate the thermal conductivity of
dSi 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 dSi.
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 dSi 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 dSi, 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 supercell finitedisplacement algorithm, we then
studied the harmonic and anharmonic lattice dynamics in MgO at highpressure. 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?? qgrids 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 qpoint sampling deserves a much
more careful investigation.
135
REFERENCES
1
G. A. Slack, Thermoelectric MaterialsNew 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.quantumespresso.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:
McGrawHill, pp. 108112 (1999).
88
W.S. Capinski, H.J. Maris, E. Bauser, I. Silier, M. AsenPalmer, 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. AsenPalmer, T. Ruf, M. Cardona,
and E. Gmelin, Appl. Phys. Lett. 71, 2109 (1997).
112
Shiichiro 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
HermannMauguin notation (HM 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 HM 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 2fold rotation
axes; improper rotation, denoted by
n
S in S notation or n in HM 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 righthand 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 facecentered cell and I for bodycentered 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, dSi and its
isotope phase silicon clathrate type II (Si
136
) both belongs to space 227, but they have
different combination of Wyckoff sites.
In dSi 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 4fold 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 Xray 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 firstprinciples 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 34atom 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)