First-principles Theoretical Study of Lattice Thermal Conductivity of
Crystals and Earth Minerals at High Temperatures and Pressures.
Moses Chou Ntam
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 4, 2012
Dissertation Abstract
First-principles Theoretical Study of Lattice Thermal Conductivity of
Crystals and Earth Minerals at High Temperatures and Pressures.
Moses Chou Ntam
Doctor of Philosophy, August 4, 2012
173 Typed Pages
Directed by Jianjun Dong
The focus of this dissertation is the first-principles calculation of lattice thermal
conductivity, a non-equilibrium thermal transport property, for a wide range of solid
materials including both ideal crystals and iron (Fe)-bearing mineral solid solutions at
high temperatures and pressures. Our computational technique combines first-principles
density functional theory (DFT), quantum scattering theory, and the Peirls-Boltzmann
transport theory (PBTT) within the single mode excitation approximation (SMEA).
Lifetimes of individual phonon modes have been directly evaluated over a wide range
of temperature-pressure conditions, including those relevant to the Earth?s deep interior,
without any empirical extrapolation. Our first-principles calculated lattice thermal
conductivity is directly derived from the lifetimes, group velocities, and heat capacities of
individual phonon modes.
An important input of our computational technique is the DFT-predicted microscopic
inter-atomic potentials. In this study, inter-atomic potentials of solid materials were
calculated up to third order using an efficient real-space super cell finite displacement
(RSSFD) algorithm that we developed recently. This computational technique has predictive
ii
capability over the range of validity of the DFT calculated atomic energy and forces. The
robustness of the predicted density dependence of harmonic force constants and third order
lattice anharmonicity tensors governs the reliability of the derived pressure dependence
of phonon lifetimes, group velocities and heat capacities. To study iron-bearing mineral
solid solutions, the vibrational virtual crystal approximation (vVCA) was implemented to
calculate the configurationally averaged harmonic force constant matrices.
The most computationally intensive step of our calculation involves the evaluation of
phonon scattering rates due to third order lattice anharmonicity. We have optimized our
algorithm for speed, efficiency and massive parallelization on supercomputers. Our algorithm
has been successfully tested on our PRISM computer cluster, the dense memory cluster
(DMC) of the Alabama super computer authority (ASA) and the Ranger super computer
at the Texas super computer center. Based on group symmetry theory, we only directly
evaluate the phonon scattering rates at the irreducible q-points in the first Brillouin zone
and rest are reconstructed from symmetry. For now our computational technique can handle
systems as large as 20-atoms per unit cell.
Our newly developed computational techniques have been successfully adopted to study
the two most abundant lower mantle minerals: silicate perovskite Mg(1?x)FexSiO3 (20-atoms
per orthorhombic unit cell) with x = 0 and x = 12.50%, and ferropericlase Mg(1?x)FexO
(2-atoms per face centered cubic unit cell) with x = 0 and x = 12.50%, as well as corundum
structured Al2O3 (10-atoms per rhomohedral unit cell).
Our calculation shows that the low frequency acoustic modes are more effective carriers
of heat compared to the optical phonon modes. In MgO and Mg(1?x)FexO the acoustic
phonon modes are seen to account for nearly 95% of the overall lattice thermal conductivity.
The effect of Fe-substitution on the thermal conductivity of Fe-free periclase has been
discussed. Fe is observed to lower the thermal conductivity of ferropericlase by significantly
reducing the phonon lifetimes of of the effective heat carrying phonon modes. The behavior
iii
of ?MgO and ?MgSiO3 have been studied. Both follow a T?1 dependence typical of insulating
materials however, the pressure increase of ?MgSiO3 is substantially weaker than ?MgO. For
example, the normalized thermal conductivity ??0 at 300 K from 0 to 135 GPa is seen to
increase from 1 to 3 for MgSiO3, while the same value goes from 1 to 8 for MgO. Implications
for heat flow in the Earth?s lower mantle have been discussed.
iv
Acknowledgments
I would like to begin by acknowledging the blessings of health, strength and a sound
mind from God without which this work would not be possible. I dedicate this dissertation
to my parents Joseph Ntam (of blessed memory) and Comfort Ntam who never had any
formal education but invested their lives in raising me and my seven siblings and laid a solid
foundation which made it possible for me to reach this far in my academic journey.
I hereby express my profound gratitude to Prof. JJ. Dong, who supervised this work, for
his time and dedication and for constantly giving me new challenges throughout this work.
His understanding and flexibility in allowing me to maintain my own schedule is appreciated.
I am grateful for the following persons: Mr Philip Forrest (physics department
IT-specialist) who maintained the PRISM computing cluster that enabled me to carry out
most of the computations reported in this dissertation. I remember times when Phil had to
work overnight and sometimes for days, to bring the computing cluster back online when
there was an outage. Without his dedication and sacrifice, this work might not have been
completed. Dr A.B. Chen who interviewed me and introduced me to the physics program at
Auburn University, Dr J.D. Perez (Physics department chair) for his support, generosity, and
encouragement throughout my studies at Auburn University, Dr Marlin Simon for advice on
teaching and who was always available to share his wisdom on future faculty career related
questions, Dr Chin-Che Tin my teaching supervisor while I taught Engineering Physics II, for
all the useful discussions we had on teaching, Dr Hinata who supervised me while I taught
general physics II for all the advice I received from him. Dr Overtoun Jenda (assistant
provost for diversity and multi-cultural affairs), for professional advice and mentoring, Dr
Florence Holland for useful discussions and career advice while I was pondering what to
v
do after my Ph.D degree. Dr Bin Xu and Dr Tang Xiaoli former members of our research
group for useful research discussions. Dr Eddy Kwessi for useful discussions on latex which
I used to typeset and compile this dissertation. Besides those mentioned here, I gratefully
acknowledge the support and prayers from all my friends during the course of this work.
Special thanks to My siblings (Mary Ewei Buh, Magaret Ndum Kah, John Mai Ntam,
Godlove Nji Ntam, Denis Beng Ntam, Elvis Cheghe Ntam and James Mueghe Ntam), for
their affection, unconditional support and prayers from Cameroon throughout the duration
of my graduate studies at Auburn University.
I am indebted to my secondary and high school physics teachers Mr Chetambong John and
Mr Mbenga Isaac for inspiring me and laying a solid foundation that enabled me to pursue
a career in physics.
Finally I would like to acknowledge financial support for my Ph.D research from the
National Science Foundation (NSF grant # EAR 075784) and the teaching assistantship
program of the Physics Department at Auburn University. I am most grateful to the physics
department at AU for the opportunity to hone my teaching skills while serving as teaching
assistant and instructor.
vi
Style manual or journal used Bibliography follows Institute of Physics (IOP) style
Computer software used The document preparation package TEX (specifically LATEX)
together with the departmental style-file auphd.sty.
vii
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii
List of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix
1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Lower Mantle Heat Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Lower Mantle Minerals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Lattice Thermal Conductivity of Lower Mantle Minerals . . . . . . . . . . . . 4
1.4 Challenges of Measuring ? at High Temperatures and Pressures . . . . . . . . 8
1.5 First-principles Calculations of Lattice Thermal Conductivity: Motivations and
Recent Progress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.6 Outline of Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2 FIRST-PRINCIPLES CALCULATION OF HARMONIC AND ANHARMONIC
LATTICE DYNAMICS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.1 Interatomic Potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2 Harmonic Lattice Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.1 Harmonic Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.2 Second Quantization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3 Third Order Lattice Anharmonicity . . . . . . . . . . . . . . . . . . . . . . . 22
2.4 Real Space Super Cell Calculation of Lattice Dynamics . . . . . . . . . . . . 24
2.5 Quasi-harmonic Approximation (QHA) and Mode Gr?uneisen Parameter . . . 27
2.6 Vibrational Virtual Crystal Approximation(vVCA) . . . . . . . . . . . . . . . 29
3 PEIRLS-BOLTZMANN TRANSPORT THEORY . . . . . . . . . . . . . . . . . 31
3.1 Peierls-Boltzmann Transport Theory (PBTT) . . . . . . . . . . . . . . . . . 31
3.2 Phonon Lifetime (?): Fermi?s Golden Rule, Single Mode Excitation
Approximation (SMEA), and Scattering Rates Due to Mass Disorder . . . . 34
3.3 Numerical Recipe for ?anh Calculation . . . . . . . . . . . . . . . . . . . . . . 37
4 Lattice Thermal Conductivity of ??Al2O3 Crystal . . . . . . . . . . . . . . . . . . 40
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.2 Lattice Dynamics and Equilibrium Thermal Properties . . . . . . . . . . . . . 41
4.2.1 Symmetry Analysis and First-principles Calculation of Total Energies . 41
4.2.2 Quasi-harmonic Phonons . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.3 Phonon Lifetimes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.4 Lattice Thermal Conductivity ? . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
viii
5 LATTICE THERMAL CONDUCTIVITY OF PV-MgSiO3 AT LOWER MANTLE
CONDITIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.2 First-principles Calculation of Total Energies, Vibrational Properties and
Phonon Lifetimes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3.1 Lattice Dynamics and Related Thermal Properties . . . . . . . . . . . . 59
5.3.2 Phonon Lifetimes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.3.3 Pressure Dependence of Thermal Conductivity . . . . . . . . . . . . . . 64
5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.4.1 Comparison of the Behavior of ?MgO and ?MgSiO3 . . . . . . . . . . . . 67
5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
6 LATTICE THERMAL CONDUCTIVITY OF FERROPERICLASE . . . . . . . . 70
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.2 Static Equation of State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.3 Vibrational Virtual Crystal Approximation of Lattice Dynamics of Ferropericlase 72
6.4 Phonon Lifetimes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.5 Lattice Thermal Conductivity of fp-Mg(1?x)FexO . . . . . . . . . . . . . . . . 75
6.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
7 Implications for Heat Flow Across the Core-Mantle Boundary (CMB) . . . . . . . 81
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
7.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
8 CONCLUSIONS AND FUTURE WORK . . . . . . . . . . . . . . . . . . . . . . . 94
Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
A Finite Difference Method For Harmonic Force Constants . . . . . . . . . . . . . . 110
A.1 Symmetry Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
A.2 Calculation of Harmonic Force Constants . . . . . . . . . . . . . . . . . . . . 111
B Interpolation of Irreducible Inter-atomic Potentials (IIIP) . . . . . . . . . . . . . . 113
B.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
B.2 Interpolation of Independent Force Constants . . . . . . . . . . . . . . . . . . 114
B.3 Volume Dependence of Harmonic Force Constants ?ij . . . . . . . . . . . . . 115
C Thermal Equation of State of Al2O3 . . . . . . . . . . . . . . . . . . . . . . . . . 124
D Equations of State (EOS) PV-MgSiO3 and Lattice Dynamics of PPV-MgSiO3 . . 127
D.1 Static Equation of State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
D.2 Thermal EOS of PV-MgSiO3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
D.3 Lattice Dynamics of PPV-MgSiO3 . . . . . . . . . . . . . . . . . . . . . . . . 131
E Density Dependence of Phonon Lifetimes in PV-MgSiO3 . . . . . . . . . . . . . . 134
ix
List of Figures
1.1 Isobaric ?latt of MgO derived based on our LDA calculated thermal equation
of state. Discrete symbols represent experimental measurements. Inset:
Normalized ?latt v.s. pressure. For the seven isotherms (300K, 500K, 1000K,
1500K, 2000K, 2500K, 3000K), ?0 is the lattice thermal conductivity at 0
GPa. Our new calculation shows a slight improvement over the previously
published results thanks to the optimization of our computational technique,
for speed and accuracy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.1 Mode Gr?uneisen Parameter(?i) derived from harmonic force constant matrix
and third order lattice anharmonicity tensor in MgO compared. Blue
circles represent mode Gr?uneisen parameter (? ) from third order lattice
anharmonicity tensor (Aijk ), orange circles represent ? from harmonic force
constants matrix (?ij ). Both ? s follow the same trend, showing that our
calculated Aijk are indeed reliable. . . . . . . . . . . . . . . . . . . . . . . 29
2.2 Vibrational Virtual Crystal Approximation (vVCA) . . . . . . . . . . . . . . 30
3.1 Energy and momentum conservation for 3-phonon processes . . . . . . . . . 35
4.1 LDA calculated density dependent harmonic force constants of Al2O3 based
on a 120-atom super cell model. A total of 1441 independent ?ij were
extracted from the HF forces from VASP calculation. Majority of the ?ij
terms are small and insensitive to volume change. Less than 20 inter-atomic
force constants have magnitudes exceeding 5 eV/?A2 and show an almost
linearly increasing pressure dependence with increasing pressure. . . . . . . 45
4.2 Vibrational phonon frequencies and density of states. Our LDA calculated
phonon spectrum (solid lines) agree well with measurements from experiment
(discrete symbols). Low frequency(acoustic) phonon modes have steeper
slopes (groupvolicities) compared tothehighfrequency (optical) modes which
are mostly flat. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.3 Pressure dependence of anharmonicity tensor elements (Aijk). Aijk elements
are observed to vary linearly with volume/density. From 0 GPa (? 84 ?A3)
to 50 GPa ( ? 75 ?A3), Aijk increases linearly from 40 eV/?A3 to 60 eV/?A3,
representing a 50% increase, ? 1% per GPa. . . . . . . . . . . . . . . . . . 47
x
4.4 Mode Gr?uneisen Parameters calculated from the force constants matrix(blue
circles) and from the third order anharmonicity tensors(orange circles)
compared . We use this information to independently test the reliability of
our calculated third order anharmonicity tensors. The extraction of Gr?uneisen
Parameters from third order lattice anharmonicity is described in section 2.5 48
4.5 Volume dependence of mode heat capacity. Lines represent the heat capacity
of individual phonon modes at the sampled at the irreducible K-points in
Brillouin zone at 6 different densities/unit cell volumes. . . . . . . . . . . . . 49
4.6 Volume/density dependence of square of phonon mode group velocities. Vg2
shows a mild volume dependence. Remember that according to knetic
transport theory, lattice thermal conductivity is given by parenleftbig? = 13summationtextCvvg2?parenrightbig 49
4.7 Density Dependence of phonon lifetimes in ?-Al2O3 the at 300K (a) & (b)
Transverse acoustic modes (c) Longitudinalacoustic phonon mode; (d)Optical
phonon modes. Lifetimes of optical phonon modes are almost pressure
insensitive. Acoustic modes increase almost linearly with decreasing volume
(increasing pressure). The optical phonon mode with the longest lifetime is
only about half the lifetime of the longest acoustic mode. This singles out
acoustic phonons as the dominant carriers of heat in ?-Al2O3. . . . . . . . . 50
4.8 Density dependence of mode-? from acoustic phonon modes at T= 300K.
Acoustic modes thermal conductivity follows a volume dependence identical
to that of the acoustic modes phonon lifetimes shown in figure 4.7 a, b
& c. This is expected as the other determining variables in the thermal
conductivity equation Cv and V2g, figures 4.5 and 4.6 are almost flat(volume
insensitive) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.9 Density dependence of mode- thermal conductivity from optical phonon
modes at T= 300K. ?optical is almost pressure insensitive, consistent with
the pressure dependence of ?optical, figure 4.7 d, and Cv and V2g, (figures 4.5
and 4.6) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.10 Density dependence of ??1. We have directly calculated the lattice thermal
conductivity of ??Al2O3 at six densities ranging from 3.98 gcc?1 to 4.51 gcc?1
and at 10 different temperature points from 300K to 2500K.??1 shows a linear
temperature dependence at constant density. At constant temperature ??1
increases almost linearly with increasing density. . . . . . . . . . . . . . . . 52
4.11 Pressure dependence of isothermal ?. Based on the V(T,T) dependence from
our LDA calculated thermal EOS we have derived the pressure dependence of
?. The numbers 1-6 represent the six unit cell volumes (densities) for which
? was calculated. ? isotherms show a linear pressure dependence. From this
pressure dependence, we can derive the temperature dependence of isobaric
? shown in figure 4.12. From the labels 1-6 one can preview the temperature
dependence of ? at conditions of constant pressure. . . . . . . . . . . . . . . 53
4.12 Isobaric ? of ?-Al2O3 at 0 GPa and 50 GPa. We show a comparison between
our calculated thermal conductivity and experiment at 0GPa. Discrete
symbols denote experiment data while solid continuous line represent our
first-principles calculation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
xi
4.13 Pressure dependence of ?. At 0 GPa, we predict nearly 3 %/GPa and 4.6
%/GPa increase in ? at 300 K and 2500K respectively. . . . . . . . . . . . . 54
5.1 Density dependence of harmonic force constants in PV-MgSiO3. Independent
force constants shown for 160-atom periodic supercells of PV. . . . . . . . . 59
5.2 LDA calculated phonon mode frequencies in PV-MgSiO3 and MgO at 0 GPa. 60
5.3 LDA calculated phonon mode frequencies in PV-MgSiO3 at 120 GPa . . . . 60
5.4 LDA calculated phonon mode group velocities in PV-MgSiO3 at four different
densities shown. Group velocities measured along the Z-direction are more
sensitive to pressure. From a density of 5.465 g/cm3 (? 120 GPa) to
4.141g/cm3 (? 0 GPa)the phonon group velocity in the ??Z direction
decreases from around 15 km/s to 11 Km/s. Velocities measured in the other
high symmetry directions X and Y are almost insensitive to pressure. . . . 61
5.5 Density dependence of 3rd order lattice anharmonicity tensor (Aijk) elements
in PV-MgSiO3. Most Aijk elements vary nearly linearly with volume/density.
Most of the elements are very small. . . . . . . . . . . . . . . . . . . . . . . 61
5.6 Mode Gr?uneisein parameters in PV-MgSiO3. Blue circles represent the
Mode Gr?uneisein parameters calculated from the second order force constants
while the orange circles show the same parameters calculated via our LDA
calculated third order lattice anharmonicity tensors. This step is used to
independently test the reliability of our lattice anharmonicity data since Aijk
is not known to have measured counterparts to which we can compare our
calculated data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.7 Volume/density dependence of phonon lifetimes in PV-MgSiO3 all 60 modes
at two high symmetry points: X (?a,0,0) and Z (0,0, ?c) in the Brillouin Zone 63
5.8 Temperature dependence of ??1 in PV-MgSiO3. ? is calculated at 4 distinct
densities corresponding to unit cell volumes shown on the figure, at 9
different temperatures ranging from 300 K to 4000K, for a total of 36
density-temperature configurations. ? is seen to follow the T?1 dependence
typical of lattice thermal conductivity in insulating materials. . . . . . . . . 64
5.9 Pressure dependence of ??1 in PV-MgSiO3 at isothermal conditions. This
is derived based on the P-V-T relations obtained from our LDA calculated
thermal EOS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.10 Main: Temperature dependence of isobaric ? in PV-MgSiO3. Inset: Pressure
variation of normalized lattice thermal conductivity. Using the coefficients
from ? = a0 + a1P + a2P2 the isobaric lattice thermal conductivity was
derived. Four isobars 0, 50, 100 and 135 GPa are shown. Inset: Pressure
variation of normalized ?. ?0 represents the thermal conductivity at 0 GPa
for each isotherm. From 300 K to 4000K, ?/?0 increases from 3 to 9 at 135
GPa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
6.1 Static EOS of ferropericlase fitted to 3rd order BM EOS . . . . . . . . . . . 71
xii
6.2 Phonon dispersion relations and mode group velocities in Mg0.875Fe0.125O
The black line/symbols shows the result obtained by replacing the mass of
Mg in the MgO calculation by the weighted mass Mg0.875Fe0.125. The blue
lines/symbols show the result obtained from the vVCA approach where the
interatomic force constants are directly calculated as described in section 2.6.
Fe substitution is seen to have the effect of softening phonon frequencies and
mode group velocities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.3 Mode Gr?unesein parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.4 Volume/density dependence of phonon lifetimes of acoustic and optical
phonon modes in fp as predicted by vVCA compared to corresponding modes
in MgO. Panels under a) show the density dependence of the lifetimes for
acoustic modes and b) shows the lifetimes of optical modes. The longest
optical lifetime is nearly a factor of 10 lower than the longest acoustic mode,
confirming the notion that acoustic modes are more effective heat transporters
than the optical modes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.5 Thermal conductivity associated with anharmonicity . . . . . . . . . . . . . 75
6.6 Temperature dependence of ??1 in MgO. ? was calculated at 5 distinct
densities corresponding to unit cell volumes shown on the figure, at 12
different temperatures ranging from 300 K to 3000K. Our calculated ? follows
a T?1 dependence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.7 Temperature dependence of ??1 in fp calculated at the corresponding unit
cell volumes for the Fe-free MgO results presented in figure 6.6. ? does not
quite follow the T?1 dependence seen in MgO. . . . . . . . . . . . . . . . . 77
6.8 Pressure dependence of isothermal ? in MgO derived from our LDA calculated
thermal EOS. For each isotherm, we used our calculated thermal EOS to
determine the pressure corresponding to the unit cell volumes 12, 13, 14, 16
& 18 ?A3 labeled 1, 2, 3, 4 & 5 along the isotherms: 300 K - 3000 K. Only
seven of the 12 isotherms are shown, for clarity. By fitting this to a second
order polynomial in P (? = a0+a1P +a2P2) we used the coefficients a0,a1,a2
from the 12 isotherms to predict ? at isobaric conditions. . . . . . . . . . . . 78
6.9 Pressure dependence of ? at isothermal conditions in fp. This pressure
dependence was derrived using the same approach described under figure
6.8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.10 Main: Temperature dependence of isobaric ? in MgO. Using the coefficients
from ? = a0 + a1P + a2P2 we derived the isobaric thermal conductivity.
Four isobars 0, 50, 100 and 150 GPa are shown. Inset: Pressure variation
of normalized ?. ?0 represents the thermal conductivity at 0 GPa for each
isotherm. From 300 K to 3000K, ?/?0 increases from 8 to 25 at 135 GPa. . . 79
6.11 Main: Temperature dependence of isobaric ? in Mg(1?x)FexO; Inset: Pressure
dependence of ??0 from 0GPa to 135 GPa. ?0 is the thermal conductivity at
0GPa for each isotherm. From 300 K to 3000K, ?/?0 increases from 2.4 to 5
at 135 GPa. This is almost a factor of 5 lower than in MgO. . . . . . . . . . 79
xiii
7.1 Vibrational phonon frequencies of (Mg1?xFex)O (a) and
(Mg1?xFex)SiO3-perovskite (b) with iron concentration x= 12.5% at
deep mantle pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7.2 Calculated lattice thermal conductivity for MgSiO3 perovskite as a function
of pressure and temperature. Calculations are shown for Fe-free crystal
(translucent symbols and lines) and mineral solid solution including 12.5% Fe
(bold color symbols and lines). Also plotted are experimental data fromOsaka
& Ito, 1991 on MgSiO3 and from Manthilake et al., 2011 on (Mg,Fe)SiO3 at
temperatures from 473 to 1073 K. . . . . . . . . . . . . . . . . . . . . . . . . 86
7.3 Isobaric thermal conductivity as a function of temperature. Computed lattice
conductivity (solid purple line) and calculated contribution from radiative
heat flow (red dashed and dotted lines, from optical absorption values
provided by Goncharov et al., 2008 and Keppler et al., 2008 respectively, as a
function of temperature at constant mid-mantle pressure of 80 GPa. The total
thermal conductivity of perovskite is given by the sum of ?rad +?latt (purple
dotted and dashed lines, corresponding to the two estimates for radiative
contributions to the thermal conductivity) . . . . . . . . . . . . . . . . . . . 87
7.4 Thermal conductivity calculated down a sample lower mantle geotherm
(Jeanloz ) and thermal boundary layer. Lattice thermal conductivity
calculated from first principles estimates for MgO (solid blue line, Tang &
Dong, 2010), and MgSiO3 perovskite (solid purple line, this study). Thin red
lines show calculated radiative contribution to thermal conductivity from two
optical absorption estimates (dashed line: Goncharov et al., 2008; . . . . . . 88
7.5 A plot summarizing the observational tradeoffs used to determine the total
heat flow across the core-mantle boundary (contoured, labels in terawatts).
Commonly assumed values from Lay et al., 2008 for the thermal boundary
thickness, the temperature at the bottom of the mantle, and the temperature
at the top of the core are shown on the right axis. Vertical solid green lines
delineate our newly-determined range of estimates for lower mantle thermal
conductivity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
B.1 Volume dependence of one inter-atomic Force constant calculated using
different finite atomic displacements (? = 0.02, 0.04, 0.06, 0.08)?A . . . . . . 115
B.2 Discrete symbols represent directly calculated values of irreducible
interatomic force constants, solid lines are the values of force constants
predicted using the fitting parameters from the calculated values according
to equation B.1. Predicted force constants coincide with those from direct
calculation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
B.3 Main: The dominant irreducible harmonic force constants corresponding
to inter-atomic forces ?Ca+2?Ca+2,?O?2?Ca+2,?O?2?O?2. Inset: Minor
inter-atomic force constants (IFC). The major inter-atomic force constants
are clearly volume dependent implying that the associated phonon frequencies
are volume dependent. The mode Gr?uneisen parameter gives the volume
dependence according to QHA theory. Most minor IFCs are volume
independent. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
xiv
B.4 Volume dependence of Force constants in SrO. Discrete symbols in the main
graph are force constants between Sr-Sr, O-Sr and O-O. The inset shows the
density dependence of the remaining IFCs . . . . . . . . . . . . . . . . . . . 121
B.5 Volume dependence of Force constants in MgO. Main: IFCs between Mg-Mg,
O-Mg and O-O. Inset: Rest of the IFCs. . . . . . . . . . . . . . . . . . . . 122
C.1 Main: Volume-pressure relationship in ??Alumina at five temperature
points: 300K, 500K, 100K, 1500K, 2000K. Inset: Static EOS fitted to 3rd
order Birch Murnaghan EOS . . . . . . . . . . . . . . . . . . . . . . . . . . 125
C.2 Thermal EOS from first-principles theoretical calculation(LDA +PAW)
compared with measurements. Discreet symbols represent measurements
done by Wachtman et al., 1962; Schauer, 1965; Amatuni et al., 1976; Adelber
& Traverse, 1984; Fiquet et al., 1999 and White & Roberts, 1983. Red solid
line represents a reproduction of LDA+PAW calculation previously done by
Bin Xu, 2009(former member of our research team) . . . . . . . . . . . . . . 126
D.1 Main: Static EOS of MgSiO3 fitted to 3rd order Birch-Murnaghan EOS.
Inset: Difference in enthalpies of PV and PPV phases plotted as a function
of pressure. Our calculation predicts that at ambient temperature, the
PV?PPV transition occurs at 101 GPa which is consistent other theoretical
predictions for example Tsuchiya et al. (2004). . . . . . . . . . . . . . . . . 128
D.2 P-V-T isotherms in PV-MgSiO3 derived within the QHA. Solid lines represent
our LDA calculation, dashed lines represent previous calculation by Oganov
et al. (2000), circles denote Experimental measurement by Knittle et al. (1987)130
D.3 Thermal EOS of PV-MgSiO3 derived within the QHA. Solid lines represent
our LDA calculation, dashed lines represent previous calculation by Karki et
al. (2000), Discrete symbols are experimental measurement by Wang et al.
(1994) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
D.4 Density dependence of harmonic force constants (eV/?A2) in PPV-MgSiO3.
Independent force constants shown for 160-atom periodic supercells of PPV. 131
D.5 LDA calculated phonon mode frequencies in PPV-MgSiO3 at 118 GPa . . . 132
D.6 LDA calculated phonon mode group velocities in PPV-MgSiO3 at 118GPa.
Mode group velocities in PPV are of the same order of magnitude as those in
PV at corresponding pressure. . . . . . . . . . . . . . . . . . . . . . . . . . 132
D.7 Density dependence of 3rd order lattice anharmonicity tensor elements
(eV/?A3) in PPV-MgSiO3. Most third order lattice anharmonicity tensors are
insensitive to the change in volume/density. A few terms show a decreasing
trend with compression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
E.1 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #1-6 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
E.2 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #7-12 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
E.3 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #13-18 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
xv
E.4 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #19-24 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
E.5 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #25-30 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
E.6 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #31-36 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
E.7 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #37-42 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
E.8 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #43-48 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
E.9 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #49-54 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
E.10 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #55-60 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
E.11 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #61-66 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
E.12 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #67-72 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
E.13 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #73-78 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
E.14 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #79-84 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
E.15 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #85-90 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
E.16 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #91-96 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
E.17 Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #97-100 in BZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
xvi
List of Tables
1.1 Layered structure of the Earth and composition . . . . . . . . . . . . . . . . 1
4.1 Space groups, formula units per primitive unit cell Z and Wyckoff sites of
??Al2O3: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.2 Irreducible atomic displacements (number of total energy calculations
required), Independent and dependent ?ij and Aijk terms from crystal symmetry. 42
4.3 Calculated lattice thermal conductivity fitted to second order in pressure at
10 different temperatures: ? = a0+ a1P + a2P2. The parameters a0,a1,a2
are then used in conjunction with the calculated thermal EOS to predict the
isobaric lattice thermal conductivity at different temperatures as shown in
figure 4.12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.1 Space groups, formula units per primitive unit cell Z and Wyckoff sites of 2
polymorphs of MgSiO3: perovskite and post-perovskite. . . . . . . . . . . . . 56
5.2 Irreducible representation of the mode symmetries at ?-point for the 2
polymorphs of MgSiO3: perovskite and post-perovskite. . . . . . . . . . . . . 57
5.3 Simulation cell volumes and corresponding pressures for PV-MgSiO3
calculation T=300 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.1 Static Equation of state of fp-Mg(1?x)FexO . . . . . . . . . . . . . . . . . . . 71
8.1 Contributions to overall ?latt from acoustic and optical phonon modes in MgO,
??Al2O3 and PV-MgSiO3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
B.1 Fitting parameters a0, a1, a2, a3 for the 52 independent harmonic force
constants of MgO fitted to third order according to equation B.1. The original
LDA calculated harmonic force constants were extracted based on a 128-atom
supercell and a finite atomic displacement ? =0.04 ?A. ii = (ia?1)?Ns+idir,
jj = (ja?1)?Ns +jdir. Where (ia, ja)= 1 to 128 and Ns=3 . . . . . . . . 123
C.1 Third-order BM-EOS parameters for ??Al2O3 . . . . . . . . . . . . . . . . . 124
C.2 Gibbs Free energy (F0), Volume per unit cell(V0), Bulk modulus (B0 ) and
derivative of bulk modulus (B?0) at 0 GPa . . . . . . . . . . . . . . . . . . . . 124
D.1 Static EOS parameters of MgSiO3 PV and PPV fitted to 3rd order BM EOS
at ambient conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
D.2 LDA calculated thermal EOS parameters for PV-MgSiO3 at 0 GPa. Values in
parentheses are from Tsuchiya et al. (2005) . . . . . . . . . . . . . . . . . . 129
xvii
E.1 Simulation cell volumes and corresponding pressures for PV-MgSiO3 at 300 K. 134
xviii
List of Abbreviations
BM Birch-Murnaghan
BZ Brillouin zone
CMB Core-mantle boundary
cond Conduction
conv Convection
CPU Central processing unit
DFT Density functional theory
e.m electromagnetic
EOS Equation of State
f.c.c. face centered cubic
GGA Generalized gradient approximation
H-F Hellmann-Feynman
IFC Inter-atomic Force Constants
IIIP Interpolation of Irreducible Inter-atomic Potentials
IPAD Irreducible paired atom displacement
ISAD Irreducible single atom displacement
xix
LD lattice dynamics
LDA Local density approximation
LM Lower mantle
LO-TO Longitudinal Optic-Transverse Optic
LTC Lattice Thermal Conductivity
MD molecular dynamics
NEMD Non-equilibrium molecular dynamics
PAW Plane augmented waves
PBTT Peirels-Boltzmann Transport Theory
QHA Quasi-harmonic approximation
rad radiation, radiative
RSSFD Real space supercell finite displacement
RT Relaxation time
RTA Relaxation time approximation
SMEA Single mode excitation approximation
TC Thermal Conductivity
TCBL Thermo-chemical boundary layer
VASP Vienna ab initio simulation package
vVCA vibrational Virtual Crystal Approximation
xx
List of Acronyms
? thermal expansivity
? phonon frequency
?i phonon frequencies
? atomic displacement
? phonon lifetime
? mode Gr?uneisen parameter
? thermal conductivity
? kinematic viscosity
D thermal diffusivity
?latt lattice thermal conductivity
?rad radiative thermal conductivity
?ij harmonic force constants matrix
Aijk third order lattice anharmonicity tensor
J heat flux
PV perovskite
PPV post-perovskite
fp ferropericlase
ModeVg mode group velocity
xxi
ModeCv mode heat capacity
T temperature
P pressure
T ?P temperature-pressure
? density
PBTT Peierls-Boltzmann transport theory
RSSFD real space supercell finite displacement
IIIP Interpolation of irreducible inter-atomic potentials
xxii
Chapter 1
INTRODUCTION
1.1 Lower Mantle Heat Flow
The Earth has a layered internal structure1. Right beneath the surface, is the crust
(about 200-300 km thick), usually categorized into continental and oceanic crust. Under the
crust, is Earth?s mantle which is conventionally divided into three regions associated with the
phase transitions of mantle minerals revealed by the seismic discontinuity. The upper mantle
is mainly made of olivine (Mg,Fe)2SiO4 , pyroxenes (Mg,Fe)SiO3 and garnet (Al-bearing
silicates). The transition zone is mainly made of spinel (MgFeAl2O4) and majorite
(Mg3(Fe,Al,Si)2(SiO4)3) minerals. The lower mantle consists of perovkite (Mg,Fe)SiO3 and
ferropericlase ((Mg,Fe)O). Between the lower mantle and the Earth?s core (innermost layer),
is the D?? layer. This layer is thought to contain the recently discovered post-perovskite2
phase of (Mg,Fe)SiO3. The core has two sub-layers, the liquid outer core made of iron-nickel
alloy, and the solid inner core made of iron. The layers, chemical composition and depth are
summarized3 in table 1.1
Table 1.1: Layered structure of the Earth and composition
Layer Depth (Km) Chemical Composition T(K) P(GPa)
crust 0-10 low density crystalline rocks 300 0
Upper mantle 10-410 olivine + pyroxene + garnet 600 0-13
Transition 410-660 spinels + majorite 13-23
Lower mantle 660-2600 perovskite + ferropericlase 1200 23-125
D?? layer 2600-2900 post-perovskite + ferropericlase 2200 125-135
Outer core 2900-5100 liquid iron alloy 4000 330
Inner core 5100-6400 solid iron 5000 364
1
Both temperature and pressure increase with the depth. At the surface, the Earth?s
pressure is 1 atm and 300 K. It is currently estimated that the pressure and temperature
reach to 364 GPa and 5000K respectively at the center of core. Governed by the second law of
thermodynamics, heat spontaneously flows from Earth?s hot core to its cool surface, through
radiation, convection, and conduction. The total heat loss at the Earths surface is relatively
well-constrained (about 46 Terawatts)4. However, because the rates of radiogenic heat
generation and secular cooling remain poorly constrained for the mantle region, independent
evaluations of the amount of heat flowing from the core into the lower mantle are much
needed to understand the global heat balance. The lower mantle heat flow is one of most
important geophysical processes, because of its close relevance to the outer core convection5
and the generation of the Earths magnetic field6. The knowledge of this process is essential
to our understanding of the origin, evolution, and present-day state of the Earths physical
system.
The amount of heat from the core is largely regulated by thermal conduction at
the thermo-chemical boundary layer (TCBL) above the core-mantle boundary (CMB). As
defined by Fouriers law of heat conduction,
J = ??T (1.1)
determination of the CMB heat flow relies on knowledge of both the temperature profile
around the CMB and thermal conductivity (?) of the lowermost mantle. thermal
conductivity (? ) also plays a crucial role in mantle dynamics. For example, ? dictates
the amount of heat absorbed by the cold subducted slabs from the hot surrounding mantle.
As ? appears in the Rayleigh number (Ra, equation 1.2), a measurement of the convective
vigor of a system, the associated effect of lower mantle viscosity controls the structure,
2
thickness and dynamics of CMB7, the style and structure of mantle convection8.
Ra = g0??Td
3
?v (1.2)
Here g is acceleration due to gravity, ? is the thermal expansion coefficient, ?T is the
temperature difference between the reference mantle temperature and the Core mantle
boundary, v is the kinematic viscosity, ? is the thermal conductivity. Although the effect
of strongly varying viscosity is paramount, recent geodynamic simulations have clearly
demonstrated that (1) the value of assumed average lower mantle conductivity is of close
relevance to the strength of mantle upwellings9, and (2) moreover the temperature derivatives
of thermal conductivity impact even more significantly the development of lower mantle
plumes10,11. Gubbins and willis et al.6 have shown that lower mantle thermal conductivity
is correlated to Earth?s magnetic field generation.
Determination of thermal conductivity of Earth minerals at relevant high pressure and
temperature conditions is clearly a subject of solid state physics research. The output of the
solid state physics study will help earth scientists to develop interlocked deep-Earth evolution
models. For example, the current estimates of QCMB are widely scattered, ranging from 5
to 15 terawatts (TW)5. The lower-bound-estimates suggest that additional sources of heat
in the mantle such as enhanced internal heating or an additional thermal boundary layer
are required to accommodate the observed budget of 47 ? 3 terawatts from the surface of
the Earth1,5. The upper-bound-estimates suggest a recent evolution of the inner core/outer
core system and high internal temperatures for early Earth conditions12.
1.2 Lower Mantle Minerals
The Earth?s lower mantle is composed predominantly of magnesium oxides and silicates.
Understanding and predicting the behavior of these compounds is key to understanding
3
the behavior and evolution of the Earth?s lower mantle. The most abundant mineral in
the Earth?s lower mantle is magnesium silicate MgSiO3, comprising more than 70% of the
lower mantle volume13. Fe2+ and/or Fe3+ distributed in these phases1 however the exact
oxidation and spin state of Fe are subject of current scientific debate. The distribution of
Fe is believed to be around 10%. Perovskite PV-MgSiO3 is the stable phase of magnesium
silicate (MgSiO3). Above 120GPa, PV-MgSiO3 transforms into a post-perovskite phase.
PPV-MgSiO3 The PV ? PPV transition is thought to explain seismic discontinuities in the
Earth?s D? layer.
Ferropericlase-(Mg,Fe)O is known to be the second most abundant mineral in the
Earth?s lower mantle constituting about 20% of lower mantle volume. It is the result of
substitution of Mg in Fe-free MgO by Fe2+. Understanding the behaviour of (Mg,Fe)O is
key to understanding the dynamics of the Earth?s lower mantle. Moreover the high-spin to
low-spin transition in (Mg,Fe)O, has been the focus of recent studies14?16.
1.3 Lattice Thermal Conductivity of Lower Mantle Minerals
Heat is a nonmechanical means of energy transfer between a system and its environment.
In the presence of a temperature difference, heat flows spontaneously from high to low
temperature regions. The heat flux J, defines the rate of heat energy transfer through
a given surface. Heat is transmitted via three mechanisms: conduction(cond), convection
(conv) and radiation (rad). As a result, the total heat flux J , has three components as shown
in equation 1.3, and heat transfer via these three mechanisms often occurs simultaneously.
J = Jconv + Jrad + Jcond (1.3)
Convective heat transfer is important in all situations where there is bulk flow of matter,
carrying heat along with the mass flow. Thus Jconv is the dominant heat transfer mechanism
4
in liquid media. The Earth?s lower mantle is not liquid, therefore Jconv is not important in
this region.
Heat transfer by radiation is ballistic in nature and it occurs via electromagnetic (e.m)
waves and does not require a material medium. According to the Stefan-Boltzmann law for
black body radiation,
Jrad ? T4 (1.4)
Due to this T4 dependence, the ballistic Jrad is a very important component of the
total heat flux in all high temperature laboratory experiments that have steep temperature
gradients and involve direct (or boundary-to-boundary or ballistic) radiative transfer wherein
light from the source warms the thermocouple with little to no participation of the medium.
Corruption of measurements intended to probe Jcond by spurious radiative transfer and mis-
interpretation of this process as diffusive has led to incorrect forms for the diffusive thermal
conductivity. Hofmeister also suggests that radiative transfer exists at lower T because (1)
direct radiative transfer increases gradually with T, i.e., no abrupt onset exists, (2) a small
radiative component is clearly seen by 350 K in laser flash analysis (LFA) experiments,
even when samples are coated with metal and/or graphite. A small amount of radiative
transfer at 298 K probably exists in most experiments17?20 Hofmeister also argues that
the ballistic radiative transfer may be large in DAC experiments21 because thermal losses
at sample interfaces provide an error in the opposite direction. In contrast, it is widely
assumed that the ballistic radiative heat transfer is not important inside Earths mantle
because mantle minerals are grainy22. The temperature gradient inside the Earth is only
about several kelvins over a kilometer. It is reasonable to approximate that minerals reach
thermal equilibrium within its grains, and radiative heat will be absorbed and randomly
scattered from grain-to-grain. At the geological length scale, the radiative heat transfer
contributes to an effective diffusive thermal conductivity. This indirect/diffusive thermal
conductivity due to radiative photons can not be directly measured in laboratory. Instead
5
they are modeled based on the optical properties of minerals, equation 1.5 The Rosseland
mean approximation for an optically thick medium conducting heat by radiative transport
is:
Qrad = ?16?n
2T3
3?R ?T = ??rad?T (1.5)
with n2<?
R>
= ?4?T3
?integraltext
0
n2?
??
dIb,?
dT d?.
Here n is the index of refraction, ? is the Stefan-Boltzmann constant, ?R is the Rosseland
mean extinction coefficient, ? is wavelength, n? is the spectral index of refraction, ?? is the
spectral extinction coefficient, and Ib,? is the Planck black body intensity function.
Heat transfer by conduction is diffusive in nature, without a bulk flow of matter and is
usually driven by a temperature gradient. Jcond is the dominant mechanism of heat transfer
in solids and is achieved via electrons, phonons (lattice vibrations) and other thermal
excitations (quasi-particles). In metallic solids, heat is conducted mainly by electrons.
According to Wiedemann-Franz (W-F) law equation 1.6 the ratio of thermal conductivity
to the electrical conductivity of a metal is proportional to the temperature(T).
?
? = LT (1.6)
Where ? is the thermal conductivity, and ? is the electrical conductivity. L is the Lorentz
number23. We can see from W-F law that as temperature increases, the thermal conductivity
increases while the electrical conductivity decreases.
In non-metallic solids the dominant mechanism of heat transfer is via lattice vibrations
(phonons). The conduction of heat via the propagation of vibrations in a solid in response
to a temperature gradient is often referred to as lattice thermal conductivity (LTC). At low
temperatures and in grainy(optically thick) and insulating media, the dominant mechanism
of heat transfer is conduction by phonons. In such media, lattice thermal conductivity lattice
thermal conductivity (?latt ), is the fundamental property that regulates the flow of heat in
6
response to a temperature difference. According to Fourier?s law ? relates the heat flux J,
through a material to the temperature gradient ?T according to equation 1.1 Even though
heat transfer is often associated with macroscopic objects, the physical properties that
constrain the flow of heat in a material are ultimately tied to the microscopic constituents
of the material. Inside crystal solids, phonons can be scattered by crystalline defects, grain
boundaries and lattice anharmonicity. The phonon lifetime (?), which characterizes the
the length of time that elapses before a phonon mode is scattered, is key to determining
the thermal conductivity of a material. Other factors necessary to determine the thermal
conductivity of a material include the phonon mode group velocities, and heat capacities.
While temperature is known to strongly affect the phonon lifetime (? ) and heat capacity,
high pressures significantly alter the elastic constants, phonon group velocities and densities
of states. As a result, first-principles study of temperature and pressure dependence of
thermal conductivity can advance our understanding of thermal transport in solid media.
The work presented in this thesis is focused on the conduction of heat in non metallic
crystal solids and Earth minerals via lattice waves (phonons). Detail treatment of the
other mechanisms of heat conduction can be found elsewhere24?29. The study of thermal
conductivity is of great technological interest because of the importance of both very high and
low thermal conductivity materials in various heat management technologies. High thermal
conductivity materials have potential applications in heat management in electronics30?32.
Low thermal conductivity materials on the other hand are the focus of recent studies in
the search for high efficient thermoelectric materials33?38 and materials for thermal-barier
coating applications39?42. The study of thermal conductivity provides useful guidelines to
the design of novel materials with unique applications in heat management.
LTC of the lower mantle ?LM regulates the heat flux across the core mantle boundary
(CMB) and impacts convection throughout the mantle. ?LM is crucial in determining the
total heat budget of the Earth. Despite its importance, ?LM is poorly constrained due to the
7
challenges in carrying out measurements at temperature and pressure conditions relevant
the Earths interior.
Understanding and constraining the thermal conductivity of a material, is an example
of a material science problem with importance in Earth science and at the heart of solid state
physics. First-principles calculations offer us new insights to material properties at extreme
conditions, yet are very challenging.
1.4 Challenges of Measuring ? at High Temperatures and Pressures
Conventional measurements of thermal conductivity (?), involve determining a
temperature gradient across a sample subjected to steady heat flow43?46. Other approaches
involve transient heating of the sample and then a study of the temperature evolution with
time. These techniques measure the thermal diffusivity (D), which is related to the ? by
equation 1.7:
D = ??C
p
(1.7)
Where ? is the density and Cp is the isobaric heat capacity. Over recent years, various
new techniques have been developed and used to measure ? or D. Contact methods like
?Angstr?om?s method47,48, Transient hot-wire technique49, and contact-free techniques like
laser flash analysis (LFA)19, Optical thermal grating technique50?52, Optical pulse transient
heating method53?56, and time-domain thermoreflectance (TDTR)57 have been developed
and used with great success to measure ? at various T,P conditions. Contrary to contact
methods where the sample is in direct contact with the heating source or the thermocouple
used to determine heat flow the sample has no direct contact in the contact-free techniques.
Measurements from contact methods are known to include systematic errors arising from
thermal contact resistance and differential thermal expansion, leading to underestimation
of ?latt. Contact-free techniques like the LFA avoid these problems but are limited by
8
requirement of very large sample size. A comprehensive review of the strengths and
weakness of these techniques is reported elsewhere19. Other challenges faced in accurately
constraining ? or D in experiment include the instability of some samples at ambient
conditions, dealing with unwanted direct radiative transfer in optically thin or partially
transparent samples, and eliminating errors associated with surface roughness and presence
of materials with P-T dependent thermophysical properties43. All these challenges often
limit the pressure range for which ? can be accurately determined in the laboratory. For
instance, the lattice thermal conductivity of ferropericlase-Mg(1?x)FexO and Fe-bearing
silicate perovskite-Mg(1?x)FexSiO3 has not been measured at lower mantle temperatures
and pressures58. Despite the success achieved in measuring equilibrium thermodynamic
properties of solids, certain temperature and pressure (T-P) conditions relevant to studies of
the Earth?s deep interior, remain inaccessible to experiments. Current geophysical estimates
are based on long extrapolations from low-pressure data11,17,43.
1.5 First-principles Calculations of Lattice Thermal Conductivity: Motivations
and Recent Progress
First-principles calculations offer an invaluable alternative approach to experimental
techniques. The terminology of first-principles or ab initio simulations is referred to
simulations that are based on quantum and statistical mechanical theory to predict the
behavior of matter at the atomistic scale and how it affects the macroscopic behavior of the
material, without any inputs and/or fitting of experimental data. A fundamental challenge of
material theory and simulation is the fact that each piece of solid material consists of around
1023 strongly interacting positively charged ions and negatively charged valence electrons.
Due to their small masses and rapid speeds, electrons inside a solid have to be described
with quantum mechanics, instead of classical Newtonian mechanics.
9
The density-functional theory (DFT) developed by Hohenberg, Kohn and Sham59,60 in
the 1960s has profound influence in materials science and engineering. This theory adopts
a series of systematical approximations to map a complicated many-electron problem into
a relatively simple independent electron problem. More importantly, this formalism allows
theorists to take advantage of fast computers andlarge storagespace to calculate a wide range
of complex quantum electronic structures (i.e. chemical bonding) using numerical methods,
and then to predict atomic structures and atomic vibration from first-principles, i.e. without
any a priori experimental knowledge. In the past twenty years, many efficient software
packages like VASP61,62, ABINIT63,64, QUANTUM EXPRESSO65 have been developed and
widely adapted to study atomic structure and equilibrium thermodynamic properties of
crystals66.
These first-principles methods have been used to study a wide range of important
material properties such as melting67 electronic structure66,68,69, phase transitions70?79,
elasticity80?84, vibrational phonon spectra85,86, structural stability87?89, Raman spectra90,91
and thermal equation of state29,92?98. Recent advances in parallel computing technologies
and declining cost of hardware and storage have made it possible to apply density
functional theory (DFT) based first-principles methods99 to study more complex material
systems over a wide range of temperature and pressure conditions. Some advantages of
predictive first-principles theoretical models include the selection of materials for important
technological applications, the design of novel materials with promising properties100 and
understanding behavior of materials key to constraining heat flow in the Earth?s deep interior.
Despite the successes recorded in the prediction of equilibrium thermal properties,
first-principles calculations of non-equilibrium properties, such as lattice thermal
conductivity or viscosity are limited due to significant increase in computational intensity.
There are no software packages available and calculations rely on time -consuming in-house
methodology development. To better constrain the heat flow in the Earth?s deep interior,
10
robust theoretical models with predictive capability are much needed. The motivation of
this work is to predict lattice thermal conductivity (?latt ) for a wide range of materials at
high temperatures and pressures relevant to the Earth?s lower mantle. Starting from simple
cubic 2-atom MgO, our goal is to study more complex mineral compositions like 10-atom
trigonal Al2O3, 20-atom orthorhombic MgSiO3, and iron-bearing mineral solid solutions like
(Mg,Fe)O.
In the last five years, several theoretical approaches58,101?108 for predicting lattice
thermal conductivity reported in literature109. These methods can be categorized into
three types: (1)Green-kubo method110,111 based on the fluctuation-dissipation theorem
and molecular dynamics simulations, (2) Non-equilibrium molecular dynamics(NEMD), (3)
Peierls-Boltzmann transport theory.
The Green-Kubo method relates the lattice thermal conductivity of a system to the time
required for temperature induced fluctuations in the instantaneous heat flux to dissipate
equation 1.8.
?ij = Vk
BT2
integraldisplay ?
0
?Ji(t)Jj(0)?dt (1.8)
where i,j are components of lattice thermal conductivity tensor, V is the volume of the
system, kB is the Boltzmann constant, T is the temperature, Ji(t) is the instantaneous heat
flux in the ith direction at time t. One major advantage of the Green-Kubo formalism is that
one can determine the entire lattice thermal conductivity tensor from one simulation. This
is particularly important in calculations of non-isotropic materials. A serious limitation to
this method is that it often requires a long time for the heat flux correlation to decay to zero
and the need to identify the self energy terms of atoms108.
Non-equilibrium molecular dynamics (NEMD) simulations impose a fixed temperature
gradient on a supercell and the heat flux required to maintain it is calculated. NEMD
11
simulations require very large simulation cells and contain noticeable finite size effects, which
sometimes are corrected with unreliable empirical models58.
Peierls-Boltzmann transport theory describes thermal conductivity in terms of specific
heat capacity (Cv), group velocities (Vg), and lifetimes ? of all the phonon modes: ? =
1
3
summationtextC
vvg2?. Both Cv and Vg are calculated based on harmonic lattice dynamics, while the
phonon lifetime ? is estimated either using equilibrium molecular dynamics simulations106,
or quantum scattering theory101.
s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48 s51s48s48s48
s48
s50s48
s52s48
s54s48
s56s48
s32s32s84s104s101s111s114s121s58s32s97s110s104s97s114s109s111s110s105s99s105s116s121s32s43s32s105s115s111s116s111s112s101
s32s32s84s104s101s111s114s121s58s32s97s110s104s97s114s109s111s110s105s99s105s116s121s32
s32s32s69s120s112s116s58s32s72s111s102s109s101s105s115s116s101s114s32s40s50s48s48s55s41
s32s32s69s120s112s116s58s32s84s80s82s67s32s112s111s108s121s32s99s114s121s115s116s97s108
s32s32s69s120s112s116s58s32s84s80s82s67s32s115s105s110s103s108s101s32s99s114s121s115s116s97s108
s32s32s69s120s112s116s58s32s83s108s97s99s107s32s40s49s57s54s50s41
s32s32s69s120s112s116s58s32s67s104s97s114s118s97s116s32s38s32s75s105s110s103s101s114s32s40s49s57s53s55s41
s40
s87
s109
s45
s49
s75
s45
s49
s41
s84s101s109s112s101s114s97s116s117s114s101s40s75s41
s48 s49s53 s51s48 s52s53 s54s48 s55s53 s57s48 s49s48s53 s49s50s48 s49s51s53
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s49s56
s50s48
s50s50
s50s52
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s49s56
s50s48
s50s50
s50s52
s32s51s48s48s75
s32s53s48s48s75
s32s49s48s48s48s75
s32s49s53s48s48s75
s32s50s48s48s48s75
s32s50s53s48s48s75
s32s51s48s48s48s75
s80s114s101s115s115s117s114s101s40s71s80s97s41
Figure 1.1: Isobaric ?latt of MgO derived based on our LDA calculated thermal equation of
state. Discrete symbols represent experimental measurements. Inset: Normalized ?latt v.s.
pressure. For the seven isotherms (300K, 500K, 1000K, 1500K, 2000K, 2500K, 3000K), ?0 is
the lattice thermal conductivity at 0 GPa. Our new calculation shows a slight improvement
over the previously published results thanks to the optimization of our computational
technique, for speed and accuracy.
12
All the first-principles calculations reported in this dissertation are based on our
newly developed computational method that combines Peierls-Boltzmann transport theory,
quantum scattering theory, and firs-principles DFT methods101,112. Figure 1.1 shows our
first-principles predicted thermal conductivity of MgO over a wide temperature-pressure
range. Our predicted values at 0 GPa are in excellent agreement with available
measurements. Our techniques are ideal for massive parallelization. During the course of my
Ph.D study, we have significantly improved the numerical efficiency of our techniques and
adopted them to study of lattice thermal conductivity, for the following minerals: ??Al2O3
(10-atoms/unit cell), MgSiO3 (20-atoms/unit cell) and ferropericlase-(Mg,Fe)O.
1.6 Outline of Dissertation
The rest of this dissertation is organized as follows: Chapter 2 covers the theory of
lattice vibrations. Chapter 3 presents the Peierls-Boltzmann Transport Theory, in chapter 4
we report the lattice thermal conductivity of corundum structured Al2O3 crystal. In chapter
5, the lattice thermal conductivity of silicate perovskite is reported. Chapter 6 presents the
vibrational virtual crystal approximation of lattice dynamics and lattice thermal conductivity
in ferropericlase Mg(1?x)FexO with x=12.5%. In chapter 7, the implications for heat flow
across the Core-Mantle Boundary (CMB) are discussed. Chapter 8 presents a summary of
the dissertation, and recommendations for future work.
13
Chapter 2
FIRST-PRINCIPLES CALCULATION OF HARMONIC AND ANHARMONIC
LATTICE DYNAMICS
2.1 Interatomic Potentials
In solids, the chemical bonds, which hold constituent atoms together, are strong yet not
100% rigid. Consequently the atoms oscillate about their equilibrium positions. When the
amplitudes of these vibrations are small compared to the interatomic spacing, the interaction
potential energy between atoms can be expanded as a function of the atomic displacements
from their equilibrium positions. Under this assumption, the total interatomic potential
energy of the crystal lattice can be expressed as follows:
V = Vstatic +
summationdisplay
?,i,?
parenleftbigg ?V
?x?,i (?)
parenrightbigg
0
x?,i (?)
+12
summationdisplay
?,i,?
summationdisplay
??,j,?
parenleftbigg ?2V
?x?,i (?)?x?,j (??)
parenrightbigg
0
x?,i (?)x?,j (??)
+16
summationdisplay
?,i,?
summationdisplay
??,j,?
summationdisplay
???,k,?
parenleftbigg ?3V
?x?,i (?)?x?,j (??)?x?,k (???)
parenrightbigg
0
x?,i (?)x?,j (??)x?,k (???)
+ 124
summationdisplay
?,i,?
summationdisplay
??,j,?
summationdisplay
???,k,?
summationdisplay
????,k,?
parenleftbigg ?4V
?x?,i (?)?x?,j (??)?x?,k (???)?x?,l (????)
parenrightbigg
0
x?,i (?)x?,j (??)x?,k (???)x?,l (????)
+??? (2.1)
14
Here V represents the total interatomic potential energy, V0 is the binding energy of the
crystal when all the atoms are at their equilibrium positions. x?,i (?) denotes the ? cartesian
component of the displacement of the ith atom in the ?th unit cell. The subscript 0 indicates
that the derivatives are evaluated with atoms at the equilibrium positions. Consequently, the
first order derivatives
parenleftBig
?V
?x?,i(?)
parenrightBig
0
, which represent the net forces on the atom Fi,at equilibrium
vanish and the second term in equation 2.1 is zero. The second order derivatives of the
total inter-atomic potential energy ??i,?j (?,??) =
parenleftBig
?2V
?x?,i(?)?x?,j(??)
parenrightBig
0
can be interpreted as
harmonic spring force constants. The elements of the harmonic force constants matrix have
the translational symmetry:
??i,?j (?,??) = ??i,?j (0,????) = ??i,?j (????,0) (2.2)
Equation 2.1 can be re-written as follows:
V = V0 + 12
summationdisplay
?,i,?
summationdisplay
??,j,?
??i,?j (?,??)x?,i (?)x?,j (??)
+16
summationdisplay
?,i,?
summationdisplay
??,j,?
summationdisplay
???,k,?
A?i,?j,?k (?,??,???)x?,i (?)x?,j (??)x?,k (???)
+ 124
summationdisplay
?,i,?
summationdisplay
??,j,?
summationdisplay
???,k,?
summationdisplay
????,k,?
B?i,?j,?k,?l (?,??,???,????)x?,i (?)x?,j (??)x?,k (???)x?,l (????)
+??? (2.3)
Where:
??i,?j (?,??) =
parenleftbigg ?2V
?x?,i (?)?x?,j (??)
parenrightbigg
0
A?i,?j,?k (?,??,???) =
parenleftbigg ?3V
?x?,i (?)?x?,j (??)?x?,k (???)
parenrightbigg
0
B?i,?j,?k,?l (?,??,???,????) =
parenleftbigg ?4V
?x?,i (?)?x?,j (??)?x?,k (???)?x?,l (????)
parenrightbigg
0
(2.4a)
15
2.2 Harmonic Lattice Dynamics
2.2.1 Harmonic Approximation
Within the harmonic approximation, third and higher order lattice anharmonicity are
neglected, and the Hamiltonian of a vibrating lattice can be expressed as:
H = V0 +KE +Vharmonic (2.5)
Where KE is the kinetic energy defined as:
KE =
summationdisplay
?,i,?
p2?,i (?)
2mi (2.6)
and
Vharmonic = 12
summationdisplay
?,i,?
summationdisplay
??,j,?
??i,?j (?,??)x?,i (?)x?,j (??) (2.7)
Here p?,i (?) and mi are the ? cartesian component of the momentum and mass of the ith
atom in the ?th unit cell respectively. ??i,?j (?,??) is defined above.
Solving the Hamiltonian of an infinitely large periodic lattice in real space is not trivial.
However, we can take advantage of the translational symmetry of the crystal and rewrite
the Hamiltonian in the reciprocal (?q) space using Fourier transformations:
X?,i (q) = 1?N
summationdisplay
?
x?,i (?)e?iq?? (2.8)
P?,i (q) = 1?N
summationdisplay
?
p?,i (?)eiq??
where ? = 1,2,??? ,3Na that labels the index of branch and Na is the number of atoms in a
primitive unit cell. N is the number of unit cells in the supercell model of the crystal, q is
16
the wave vector in the reciprocal space.
Note that X?,i (q) and P?,i (q) also have translational symmetry in reciprocal space:
X?,i (q + g) = X?,i (q)
P?,i (q + g) = P?,i (q)
(2.9)
Here g is a reciprocal lattice point and eiq?? = 1. Since x?,i (l) and p?,i (l) are real, we can
show that X?? (q) = X? (?q) and P?? (q) = P? (?q)
Conversely, in real space,
x?,i (?) = 1?N
summationdisplay
q
X?,i (q)eiq?? (2.10)
p?,i (?) = 1?N
summationdisplay
q
P?,i (q)e?iq??
With these transformations, the equation 2.6 becomes:
KE = 12
summationdisplay
q,q?,i,?
1
miP?,i (q)P?,i (q
?) 1
N
summationdisplay
?
e?i(q+q?)??
= 12
summationdisplay
q,q?,i,?
1
miP?,i (q)P?,i (q
?)?q+q?,0
= 12
summationdisplay
q,i,?
1
miP?,i (q)P?,i (?q)
= 12
summationdisplay
q,i,?
1
miP?,i (q)P
?
?,i (q)
= 12
summationdisplay
q,i,?
|P?,i (q)|2
mi
17
Letting h = ???? i.e ? = ?? + h
Equation 2.7 becomes
Vharmonic = 12
summationdisplay
i,?
summationdisplay
j,?
summationdisplay
h
??i,?j (h,0)
summationdisplay
q,q?
X?,i (q)X?,j (q?)eiq?h 1N
summationdisplay
??
ei(q+q?)???
= 12
summationdisplay
i,?
summationdisplay
j,?
summationdisplay
q
summationdisplay
h
??i,?j (h,0)X?,i (q)X?,j (?q)eiq?h (2.11)
We define
??i,?j (q) =
summationdisplay
h
??i,?j (h,0)eiq?h
Note that,
???i,?j (q) = ??i,?j (?q) = ??j,?i (q) (2.12)
So the equation 2.11 becomes
Vharmonic = 12
summationdisplay
i,?
summationdisplay
j,?
summationdisplay
q
??i,?j (q)X?,i (q)X?,j (?q) (2.13)
The Hamiltonian can now be expressed as:
H = V0 +
summationdisplay
q
H (q) (2.14)
Where
H (q) = 12|P?,i (q)|
2
mi +
1
2
summationdisplay
j,?
??i,?j (q)X?,i (q)X?,j (?q)
(2.15)
18
We can solve the equation of motion of lattice vibration using the canonical equations:
?H
?X?,i (q) = ?
?P?,i (q)
?H
?P?,i (q) =
?X?,i (q) (2.16)
Now,
?H
?X?,i (q) =
1
2
summationdisplay
j,?
??i,?j (q)X??,j (q) + 12
summationdisplay
j,?
??j,?i (?q)X??,j (q)
(2.17)
Using equation 2.12, the first canonical equation becomes:
?H
?X?,i (q) =
summationdisplay
j,?
??i,?j (q)X??,j (q) = ? ?P?,i (q) (2.18)
Similarly, the second canonical equation becomes:
?H
?P?,i (q) =
P??,i (q)
m? =
?X?,i (q) (2.19)
Using equation 2.12 we obtain,
?X?,i (q) + 1
m?
summationdisplay
j,?
???i,?j (q)X?,j (q) = 0 (2.20)
Let
X?i (q) = 1?m
?
ei?(q)t?e?i (q) (2.21)
19
equation 2.20 becomes
??2?
m??e?,i (q) +
1
m?
summationdisplay
j,?
???i,?j (q) 1?m
?
?e?,j = 0 (2.22)
ie
summationdisplay
j,?
bracketleftbigg??
?i,?j (q)?
m?m? ??
2???
bracketrightbigg
?e?,j = 0 (2.23)
Where ?e?i (q) and ?(q) are the eigenvectors and eigenfrequencies of the dynamic matrix
D(q) which is defined as:
D?i,?j (q) = ?
?
?i,?j (q)?
m?m? =
summationdisplay
h
??i,?j (h,0)?
m?m? e
?iq?h (2.24)
For each q-point, a D(q) matrix has 3Na eigenvalues/eigenvectors. The eigenvalues
represent the frequencies (?) of normal modes of vibrations. The dependence of the vibration
frequencies on the wave vectors is often referred to as phonon dispersion relation and the
slope of a dispersion curve (?q?(q)) represents the velocity of propagation of the phonon
through the crystal lattice also known as the phonon group velocity Vg (q). Of the 3Na
modes, three are usually the acoustic modes while the remaining 3Na?3 are optical modes.
2.2.2 Second Quantization
I now introduce normal mode coordinates,
Q? (q) =
summationdisplay
i,?
?m
ie??,i (q,?)?X?,i (q) (2.25)
P? (q) =
summationdisplay
i,?
1?
mie?,i (q,?)?P?,i (q)
20
Conversely,
X? (q) =
summationdisplay
i,?
1?
mie?,i (q,?)?Q? (q) (2.26)
P? (q) =
summationdisplay
i,?
?m
ie??,i (q,?)?P? (q)
Note that Q?i (q) = Qi (?q) and P?i (q) = Pi (?q).
e?,i (q,?) is an eigenvector with the following property:
summationdisplay
i,?
e?,i (q,?)?e??,i (q,?) = ???
After some derivation it can be shown that within the harmonic approximation, the
Hamiltonian takes the form:
H (q) = 12
summationdisplay
?
parenleftbig|P
? (q)|2 +?2? (q)|Q? (q)|2
parenrightbig (2.27)
Under the second quantization, the creation and annihilation operators are defined in terms
of the normal coordinates as follows:
?a? (q) = 1?2planckover2pi1?
? (q)
parenleftBig
?? (q) ?Q? (q) +i?P? (?q)
parenrightBig
(2.28)
?a+? (q) = 1?2planckover2pi1?
? (q)
parenleftBig
?? (q) ?Q? (?q)?i ?P? (q)
parenrightBig
(2.29)
Conversely
?Q? (q) = planckover2pi1?
2planckover2pi1?? (q)
parenleftbig?a+
? (?q) + ?a? (q)
parenrightbig
?P? (q) = i
?planckover2pi1?
? (q)
2
parenleftbig?a+
? (q)??a? (?q)
parenrightbig (2.30)
21
Here
bracketleftBig?
Q? (q), ?P? (q)
bracketrightBig
= iplanckover2pi1?????qq?
bracketleftbig?a
? (q),?a+? (q)
bracketrightbig = ?
???qq? (2.31)
In terms of ?a? (q) and ?a+? (q) the Hamiltonian takes the form:
?H (q) = summationdisplay
?
planckover2pi1?? (q)
parenleftbigg
?a+? (q)?a? (q) + 12
parenrightbigg
(2.32)
=
summationdisplay
?
planckover2pi1?? (q)
parenleftbigg
?n? (q) + 12
parenrightbigg
(2.33)
?n? (q) is the population of phonons with wave vector q in the ?th mode.
2.3 Third Order Lattice Anharmonicity
From the expansion of the total interatomic potential energy of a crystal in terms of
atomic displacements, all terms of order 3 or higher are associated with the anharmonicity
of the lattice. The third order lattice anharmonicity is the leading term that causes a
perturbation to the harmonic lattice vibrations and anharmonicity induced phonon-phonon
scattering.
Recall from 2.3 that we defined the third order lattice anharmonicity tensor as:
A?i,?j,?k (?,??,???) =
parenleftbigg ?3V
?x?,i (?)?x?,j (??)?x?,k (???)
parenrightbigg
(2.34)
22
In reciprocal space the third term in equation 2.3 becomes
H(3)anh = 16
summationdisplay
?,i,?
summationdisplay
??,j,?
summationdisplay
???,k,?
A?i,?j,?k (?,??,???)x?,i (?)x?,j (??)x?,k (???)
= 16
summationdisplay
q,i,?
summationdisplay
q?,j,?
summationdisplay
q??,k,?
A?i,?j,?k (?,??,???)X?,i (q)X?,j (q?)X?,k (q??)ei(q??+q????+q??????)
(2.35)
Let h = ????? ie ? = ??? + h
and h? = ?? ???? ie ?? = ??? + h?
H(3)anh = 16
summationdisplay
q,i,?
summationdisplay
q?,j,?
summationdisplay
q??,k,?
summationdisplay
???
A?i,?j,?k (q,q?)X?,i (q)X?,j (q?)X?,k (q??)ei(q+q?+q??)????
= 16
summationdisplay
q,i,?
summationdisplay
q?,j,?
summationdisplay
q??,k,?
A?i,?j,?k (q,q?)X?,i (q)X?,j (q?)X?,k (q??)?q+q?+q??,g
(2.36)
Where
A?i,?j,?k (q,q?) =
summationdisplay
h,h?
A?i,?j,?k (h,h?,0)ei(q?h+q??h?) (2.37)
We have X? (q)X? (q?)X? (q??)
= 1?m
?m?m?
summationdisplay
?,?,?
Q? (q)Q? (q?)Q? (q??)e?,? (q)?e?,? (q?)?e?,? (q??)
= 1?m
?m?m?
parenleftbiggplanckover2pi1
2
parenrightbigg3
2 summationdisplay
?,?,?
e?,? (q)?e?,? (q?)?e?,? (q??)radicalbig
?(q)?(q?)?(q??) Y??? (q,q
?,q??)
(2.38)
23
Here
Y??? = parenleftbig?a+? (?q) + ?a? (q)parenrightbigparenleftbig?a+? (?q?) + ?a? (q?)parenrightbigparenleftbig?a+? (?q??) + ?a? (q??)parenrightbig
(2.39)
Using equations 2.38 and 2.39, we have implemented an algorithm to calculate the
phonon-phonon scattering rates due to lattice anharmonicity within the single mode
excitation approximation (SMEA). The details of the SMEA are reported in chapter 3.
2.4 Real Space Super Cell Calculation of Lattice Dynamics
We developed a real space supercell finite displacement ( RSSFD) method to evaluate
interatomic potentials up to third order from first-principles. In the RSSFD approach, we
displace atoms in a real space super cell by a series of finite displacements and then evaluate
harmonic force constants and third order lattice anharmonicity tensor elements from forces
obtained via the Hellmann-Feynman (H-F) theorem.
F?,i (?) = ? ?E?u
?,i (?)
= ?
angbracketleftbigg
?
vextendsinglevextendsingle
vextendsinglevextendsingle ?H
?u?,i (?)
vextendsinglevextendsingle
vextendsinglevextendsingle?
angbracketrightbigg
(2.40)
where ?,i is the component of the force and atom index respectively, ? is the unit cell index,
H is the DFT hamiltonian and ? is the electronic ground-state wavefunction.
This technique has been successfully used to predict the full phonon spectrum of many
materials78,79,85,101,112?117. The computational cost of the RSSFC method depends on the
size of supercell and the range of of interatomic force constants. For ionic solids, the effect
of induced Born effective charges due to lattice vibrations is often added as a correction to
the phonon spectrum. Also referred to as LO-TO splitting, this often requires a separate
calculation following a method proposed by Kunc and Martin118.
24
Our RSSFD technique involves two types displacement schemes: The first scheme which
we call irreducible single atom displacement (ISAD) involves displacement of single atoms
based on crystal symmetry, to determine the H-F forces needed to derive the harmonic
force constants. The second algorithm called irreducible paired atom displacement (IPAD)
involves the displacement of atoms in pairs, to determine the H-F forces needed to evaluate
the third order lattice anharmonicity tensors. The details of the ISAD are reported here119.
In this section I will describe the IPAD algorithm for the calculation of the third order lattice
anharmonicity tensor A?i,?j,?k (?,??,???).
In our IPAD scheme, we displace a a pair of atoms simultaneously by a finite amount ?
which is much smaller than the interatomic distance. For example if we displace the jth atom
in ? direction and the kth atom in ? direction by (?,?), (?,??), (??,?), and (??,??),
25
respectively. The ? component of the H-F forces on the ith atom are:
F++?,i (?) = ?(??i,?j (?,??) +??i,?k (?,???))??
?12 (A?i,?j,?j (?,??,??) +A?i,?k,?k (?,???,???) + 2A?i,?j,?k (?,??,???))??2
?16
?
?? B?i,?j,?j,?j (?,??,??,??) +B?i,?k,?k,?k (?,???,???,???)
+2B?i,?j,?j,?k (?,??,???,????) + 2B?i,?j,?k,?k (?,??,???,????)
?
????3
(2.41a)
F+??,i (?) = ?(??i,?j (?,??)???i,?k (?,???))??
?12 (A?i,?j,?j (?,??,??) +A?i,?k,?k (?,???,???)?2A?i,?j,?k (?,??,???))??2
?16
?
?? B?i,?j,?j,?j (?,??,??,??)?B?i,?k,?k,?k (?,???,???,???)
?2B?i,?j,?j,?k (?,??,???,????) + 2B?i,?j,?k,?k (?,??,???,????)
?
????3
(2.41b)
F?+?,i (?) = ?(???i,?j (?,??) +??i,?k (?,???))??
?12 (A?i,?j,?j (?,??,??) +A?i,?k,?k (?,???,???)?2A?i,?j,?k (?,??,???))??2
?16
?
?? ?B?i,?j,?j,?j (?,??,??,??) +B?i,?k,?k,?k (?,???,???,???)
+2B?i,?j,?j,?k (?,??,???,????)?2B?i,?j,?k,?k (?,??,???,????)
?
????3
(2.41c)
F???,i (?) = +(??i,?j (?,??) +??i,?k (?,???))??
?12 (A?i,?j,?j (?,??,??) +A?i,?k,?k (?,???,???) + 2A?i,?j,?k (?,??,???))??2
+16
?
?? B?i,?j,?j,?j (?,??,??,??) +B?i,?k,?k,?k (?,???,???,???)
+2B?i,?j,?j,?k (?,??,???,????) + 2B?i,?j,?k,?k (?,??,???,????)
?
????3
(2.41d)
From equation 2.41, simple derivation yields
A?i,?j,?k (?,??,???) = ?F
++
?,i (?) +F
+?
?,i (?)?F
?+
?,i (?) +F
??
?,i (?)
4?2 (2.42)
26
2.5 Quasi-harmonic Approximation (QHA) and Mode Gr?uneisen Parameter
Within the harmonic approximation, lattice vibrations are treated as a collection of
independent harmonic oscillators. Phonons are quasi-particles that represent the quantized
thermal excitations of harmonic lattice vibrations. According to this approximation, phonons
do not interact with each other, and they have infinitely large lifetimes. Also, phonon
frequencies do not depend on vibration amplitudes and unit cell volumes. i.e. they are
temperature and pressure independent. This approximation however, cannot explain several
material properties associated with lattice anharmonicity such as small yet finite thermal
expansion and finite phonon lifetimes. When anharmonic contributions to the interatomic
potentials is small, anharmonic effects in a crystal be evaluated within the framework of the
quasi-harmonic approximation(QHA), which still describes inter-atomic potentials within
the harmonic approximation, but at the same time, force constant matrices are treated as
volume dependent. Consequently, phonon frequencies are volume dependent. The volume
dependence of the vibrational frequencies is described by a parameter known as the mode
Gr?uneisen parameter (? ). From QHA theory, ? is defined as follows:
?iparenleftbig?qparenrightbig = ?dln?i
parenleftbig?qparenrightbig
dlnV = ?
V
?iparenleftbig?qparenrightbig
d?iparenleftbig?qparenrightbig
dV (2.43)
By invoking Hellman-Feynman theorem, ? can be expressed in terms of the logarithmic
volume derivative of the dynamical matrix
?iparenleftbig?qparenrightbig= ? V
2?iparenleftbig?qparenrightbig2
angbracketleftbige
i(
?q)vextendsinglevextendsingledD
parenleftbig?qparenrightbig
dV
vextendsinglevextendsinglee
i(
?q)angbracketrightbig (2.44)
Where the dynamical matrix D(vectorq) is defined in equation 2.24
This expression can be rewritten as follows:
27
?iparenleftbig?qparenrightbig= ? 1
2?iparenleftbig?qparenrightbig2
summationdisplay
??
summationdisplay
kk?
dD??(?q)
dlnV
e??ki (?q)e?k?i (?q)?
MkMk? (2.45)
The mode Gr?uneisen parameters (?) are used to quantify the degree of lattice
anharmonicity in a crystal because they are determined by the volume dependence of
the dynamical matrix D(vectorq) in reciprocal space or force constants matrix ?(0,vectorh) in real
space. According to perturbation theory120,121, the logarithmic volume derivative of the
force constants matrix ?(0,vectorh) is directly associated with third order lattice anharmonicity
according to equation 2.46
parenleftBigg
d?(0,vectorh)
dlnV
parenrightBigg
??
= 13
summationdisplay
?
summationdisplay
l?l??
A???(0,h,h?)r?(h?) (2.46)
By combining equations 2.46, 2.45 and 2.24, it is straight forward to show that the mode
Gr?uneisen parameter can be derived from the third order lattice anharmonicity tensors
according to equation 2.47
?iparenleftbig?qparenrightbig= ? 16?
i2
(?q)
summationdisplay
k
summationdisplay
k?h
summationdisplay
k??h?
summationdisplay
???
A?k0,?k?h,?k??h?(0,h,h?)e
??k
i (
?q)e?k?
i (
?q)
?m
kmk?
ei?q?
?h
r?k??h? (2.47)
In our first-principle calculation of interatomic potential energy of solids, ? serves as an
important internal check onthe robustness of our calculated third order lattice anharmonicity
tensor (Aijk ). This is because unlike phonon frequencies and group velocities which have
experimental measurements against which our calculated values can be compared, the third
order lattice anharmonicity tensor do not have experimentally measured counterparts. To
independently test the reliability of our calculated Aijk data, we evaluate ? first from the
calculated dynamic matrices which have been verified to be accurate by comparing with
previous measurements and calculations where these exist. Next we evaluate ?i via the
28
calculated third order lattice anharmonicity. By comparing ? independently evaluated from
the two methods, we can establish the reliability of our Aijk data. This is illustrated in
Figure 2.1. Here ? from ?ij and Aijk for MgO are plotted on the same axes. It is observed
that the two ? capture the same features. Generally, the quality of the agreement in the two
sets of ? reflects the reliability of the calculated third order lattice anharmonicity tensors.
s48
s49
s50
s51
s52
s32
s77
s111
s100
s101
s32
s71
s114
s252
s110
s101
s105
s115
s101
s110
s32
s80
s97
s114
s97
s109
s101
s116
s101
s114
s40
s105
s41
s88 s76
s76
s77s103s79 s80s32s61s32s48s32s71s80s97
s32s68s101s114s105s118s101s100s32s102s114s111s109s32s76s97s116s116s105s99s101s32s65s110s104s97s114s109s111s110s105s99s105s116s121s32s84s101s110s115s111s114
s32s68s101s114s105s118s101s100s32s102s114s111s109s32s72s97s114s109s111s110s105s99s32s70s111s114s99s101s32s67s111s110s115s116s97s110s116s115s32s77s97s116s114s105s120
Figure 2.1: Mode Gr?uneisen Parameter(?i) derived from harmonic force constant matrix and
third order lattice anharmonicity tensor in MgO compared. Blue circles represent ? from
Aijk , orange circles represent ? from ?ij . Both ? s follow the same trend, showing that
our calculated Aijk are indeed reliable.
2.6 Vibrational Virtual Crystal Approximation(vVCA)
Lower mantle minerals are not pristine crystals but solid solutions that contain
compositional disorder. For example, about 10% of Fe2+ or Fe3+ exist in perovskite and
ferropericlase. Figure 2.2 shows a sketch of differences between 1-dimensional crystals (AO
and BO) and 1-dimensional pseudo binary solid solutions, (A1?x,Bx)O where both A and
B atoms are randomly distributed at the cation sites. In Figure 2.2, ?AOij and ?BOij represent
the harmonic force constants in perfect crystals of AO and BO respectively. ?(SS)ij represents
29
the force constants in the disordered mineral solid solution containing a random substitution
of cation-A with cation-B. Both ?AOij and ?BOij have translational symmetry. But the A/B
substitution disorder breaks the ?(SS)ij translational symmetry.
Unit Cell
(a) Periodic Perfect Crystal AO, ?(AO)ij
(b) Periodic Perfect Crystal BO, ?(BO)ij
(c) Non-periodic Solid Solution (A,B)O, Compositional Disorder, ?(SS)ij
(d) Periodic Virtual Crystal XO, Ensemble Average: ?(XO)ij ? ??(SS)ij ?
Figure 2.2: Vibrational Virtual Crystal Approximation (vVCA)
An exact solution for the lattice dynamics of mineral solid solutions requires large
supercell models. As an approximation, we implemented a vibrational virtual crystal
approximation (vVCA). Within the vVCA, we first generate a supercell model with a
specific A/B substitution configuration ?(ss?1)ij Then we generate all the symmetrically
equivalent substitution configurations, and derive their force constant matrices with group
theory. ?(ss?2)ij ,?(ss?3)ij ,?(ss?4)ij ??? are derived with ?(A,B)Oij = ?(A(1??)B?)Oij . We perform a
configurational average which represents the force constants in the virtual crystal of XO.
Notice that crystalline periodicity is restored in ??ssij?vV CA. It is worth pointing out that
??ssij?vV CA is not simply the arithmetic average: (1?x)?AOij +x?BOij . It includes a much more
realistic description of A/B substitution.
This substitution is known to modify the interatomic forces especially in the vicinity of
the substitution. In this study we find that the bond modification leads to a softening of
phonon frequencies and group velocities, details reported in chapter 6.
30
Chapter 3
PEIRLS-BOLTZMANN TRANSPORT THEORY
3.1 Peierls-Boltzmann Transport Theory (PBTT)
All thermally excitable particles/quasi-particles contribute to heat conduction. In
optically thick and electrically insulating solids, the main carriers of heat are phonons, i.e.
the quasi-particles that represent quantization of lattice vibrations. At the microscopic level,
the kinetic particle theory describes the heat flux according equation 3.1
vectorJ = 1
N?q
summationdisplay
?q,i
n(?q,i)planckover2pi1?(?q,i)vectorvg(?q,i) (3.1)
Here, n(vectorq,i) represents the phonon occupation (distribution) function of mode (vectorq,i), vectorvg(vectorq,i)
is the group velocity and ?(vectorq,i) is the phonon frequency. At thermal equilibrium, n = n0,
where no is the well-known Bose-Einstein equilibrium occupation (distribution) function for
phonons:
no = 1eplanckover2pi1?/k
BT ?1
(3.2)
The net heat flux at thermal equilibrium vectorJ0 is zero:
vectorJ0 = 1
N?q
summationdisplay
?q,i
n0(?q,i)planckover2pi1?(?q,i)vectorvg(?q,i) = 0 (3.3)
31
According to the Peierls-Boltzmann thermal transport theory, the variation of the
phonon distribution function with time can be expressed as:
parenleftbiggdn
dt
parenrightbigg
=
parenleftbiggdn
dt
parenrightbigg
scattering
(3.4)
The total derivative parenleftbigdndtparenrightbig contains three terms:
parenleftbiggdn
dt
parenrightbigg
= ?n?t +
parenleftbiggdn
dt
parenrightbigg
diffusion
+
parenleftbiggdn
dt
parenrightbigg
field
(3.5)
For heat conduction, an external field is absent, i.e. we have (dndt)field = 0.
At steady state conditions, ?n?t = 0. Therefore,
parenleftbiggdn
dt
parenrightbigg
diffusion
=
parenleftbiggdn
dt
parenrightbigg
scattering
(3.6)
The diffusion term can be expressed as:
(dndt)diffusion = vectorvg ? vector?n(vectorr,vectorq,t)
?
parenleftBig
vectorvg ? vector?T
parenrightBig?n0
?T
= Cvplanckover2pi1?
parenleftBig
vectorvg ? vector?T
parenrightBig
(3.7)
where Cv = ???T = ??T parenleftbign0 + 12parenrightbigplanckover2pi1? = planckover2pi1?
parenleftBig
?n0
?T
parenrightBig
is the phonon mode heat capacity.
Within the relaxation time approximation (RTA), the scattering term in equation 3.6
can be approximated by equation 3.8
parenleftbiggdn
dt
parenrightbigg
scattering
? ?
parenleftbiggn?no
?
parenrightbigg
= ??n? (3.8)
32
Consequently we have,
?n(vectorq,i)
?(vectorq,i) = ?
Cv(vectorq,i)
planckover2pi1?(vectorq,i)
parenleftBig
vectorvg(vectorq,i)? vector?T
parenrightBig
(3.9)
The heat flux in equation 3.1 can be rewritten as:
vectorJ = 1
Nvectorq
summationdisplay
vectorq,i
parenleftbign0(vectorq,i) + ?n(vectorq,i)parenrightbigplanckover2pi1?(vectorq,i)vectorv
g(vectorq,i) (3.10)
Using equation 3.3, (where vectorJ0 = 0) in conjunction with equation 3.9, we can further simplify
equation 3.10:
vectorJ = 1
Nvectorq
summationdisplay
vectorq,i
parenleftbigg
??(vectorq,i)Cv(vectorq,i)planckover2pi1?(vectorq,i)
parenleftBig
vectorvg(vectorq,i)? vector?T
parenrightBigparenrightbigg
planckover2pi1?(vectorq,i)vectorvg(vectorq,i)
= ? 1N
vectorq
summationdisplay
vectorq,i
Cv(vectorq,i)?(vectorq,i)
parenleftBig
vectorvg(vectorq,i)? vector?T
parenrightBig
vectorvg(vectorq,i)
(3.11)
The thermal conductivity tensor ? of a single crystal is defined as:
vectorJ = ??vector?T (3.12)
Comparing equation 3.12 with equation 3.10, we can clearly see that:
??? = 1N
vectorq
summationdisplay
?q,i
Cv(?q,i)vg?(?q,i)vg?(?q,i)?(?q,i) (3.13)
33
We can show that the non-diagonal terms ??? = 0 for ? negationslash= ?.
The averaged thermal conductivity in a randomly oriented polycrystal is defined as:
? = 13 (?11 +?22 +?33) = 13N
vectorq
summationdisplay
vectorq,i
Cv(vectorq,i)vectorVg(vectorq,i)2?(vectorq,i) (3.14)
The heat capacity and group velocity are both determined from the phonon spectrum
and vibrational density of states of the solid. The challenge in determining ? lies with the
evaluation of ?(vectorq,i) from first-principles.
3.2 Phonon Lifetime (?): Fermi?s Golden Rule, Single Mode Excitation
Approximation (SMEA), and Scattering Rates Due to Mass Disorder
Phonon relaxation time (RT) is a quantity associated with non-equilibrium transport
process. However, in practice, RT is often considered as equivalent to the phonon lifetime
which is an equilibrium quantity, because both quantities are related to the rate of phonon
scattering. In this dissertation,the phrases ?phonon lifetime? and ?phonon relaxation time?
are used interchangeably.
Inside an insulating crystal, phonons can be scattered by crystalline defects
(compositional disorder), grain boundaries or lattice anharmonicity. In the work reported in
this dissertation, we considered two types of perturbations to the harmonic lattice dynamics
namely: mass disorder (either isotope or Fe/Mg substitution) and intrinsic third order lattice
anharmonicity. We neglect phonon scattering due to grain boundary since its contribution is
comparably small at room temperature and above. The overall phonon lifetime is obtained
from the sum of the scattering rates due to the individual scattering processes according to
Matthiessen?s rule.
parenleftbig??1parenrightbig
eff =
parenleftbig??1parenrightbig
3?phonon +
parenleftbig??1parenrightbig
mass (3.15)
34
According to Fermi?s golden rule based on quantum time dependent perturbation theory,
the transition rate (?fi ) from an initial quantum state to the final state |f?, under a
perturbation ?H is defined by equation 3.16.
?fi = 2?planckover2pi1 |?f | ?H|i?|2?(?f ??i) (3.16)
|i? and |f? are the initial and final quantum states with energies ?i and ?f respectively. In
the case of third order lattice anharmonicity induced phonon scattering, ?H is the same as
H(3)anh defined in equation 2.36 Considering only scattering due to 3-phonon processes, we
can rewrite scattering rate equation as follows
?ni(?q)
?i(?q) = ?
parenleftbiggdn
i(
?q)
dt
parenrightbigg
3?phonon
(3.17)
Phonon-phonon scatterings due to third order lattice anharmonicity involves 2 phonons
combining into one phonon or one phonon splitting up into two phonons as depicted in
figure 3.1. These are commonly referred to as 3-phonon processes.
q
q??
q?
q
q??
q?
planckover2pi1? = planckover2pi1?? + planckover2pi1???
q +G = q? +q??
planckover2pi1? + planckover2pi1?? = planckover2pi1???
q +G+q? = q??
Figure 3.1: Energy and momentum conservation for 3-phonon processes
35
Based on equation 3.16, 2.36, 2.37, 2.38 and 2.39, we can derive that:
parenleftBig
dni(?q )
dt
parenrightBig
3?phonon
=
summationtext
?q ?
summationtext
?q ??
??
?
??
bracketleftBigparenleftbig
ni(?q) + 1parenrightbig
parenleftBig
nj(?q?) + 1
parenrightBigparenleftBig
nk(?q??)
parenrightBig
?parenleftbigni(?q)parenrightbig
parenleftBig
nj(?q?)
parenrightBigparenleftBig
nk(?q??) + 1
parenrightBigbracketrightBig
S?q
??
?q,?q ?+
1
2
bracketleftBigparenleftbig
ni(?q) + 1parenrightbig
parenleftBig
nj(?q?)
parenrightBigparenleftBig
nk(?q??)
parenrightBig
?parenleftbigni(?q)parenrightbig
parenleftBig
nj(?q?) + 1
parenrightBigparenleftBig
nk(?q??) + 1
parenrightBigbracketrightBig
S?q
?,?q ??
?q
??
?
??
(3.18)
Where S?q
??
?q ,?q ? and S
?q ?,?q ??
?q relate the scattering rates to lattice anharmonicity and phonon
spectrum, also enforcing the conservation of energy and momentum.
S?q
?,?q ??
?q =
h2
242?N?(?q,i)?(?q?,j)?(?q??,k)
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingle
summationdisplay
???
A???(?q,?q?)ei?(q)ej?(?q?)ek?(?q??)
(m?m?m?)1/2
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingle
2
??[planckover2pi1?(?q,i)?planckover2pi1?(?q?,j)?planckover2pi1?(?q??,k)]??q??q ???q ??,?g (3.19)
S?q
?,?q ??
?q =
h2
242?N?(?q,i)?(?q?,j)?(?q??,k)
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingle
summationdisplay
???
A???(?q,?q?)ei?(q)ej?(?q?)ek?(?q??)
(m?m?m?)1/2
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingle
2
??(planckover2pi1?(?q,i) + planckover2pi1?(?q?,j)?planckover2pi1?(?q??,k))??q+?q ???q ??,?g (3.20)
The single mode excitation approximation (SMEA) assumes that only one phonon mode
is displaced out of equilibrium and its phonon occupation number relaxes back to equilibrium
while phonon occupation numbers of all other modes maintain their equilibrium value. In
terms of the phonon distribution functions, this can be illustrated by equation 3.21
ni(?q) = nio(?q) + ?ni(?q)
nj(?q?) ? njo(?q?)
nk(?q??) ? nko(?q??)
(3.21)
36
After enforcing energy conservation and neglecting terms of order 2 and higher, equation
3.17 can be expressed as:
parenleftBig
dni(?q)
dt
parenrightBig
3?phonon
=
?ni(?q)?summationtext
?q ?
summationtext
?q ??
braceleftBigparenleftBig
nko(?q??)?njo(?q?)
parenrightBig
S?q
??
?q ,?q ? ?
1
2
parenleftBig
njo(?q?) +nko(?q??) + 1
parenrightBig
S?q
?,?q ??
?q
bracerightBig (3.22)
The scattering rate (the inverse of phonon lifetime) associated with third order lattice
anharmonicity can be expressed as:
1
?anh(?q,i) = ?
summationdisplay
?q ?
summationdisplay
?q ??
braceleftbiggparenleftBig
nko(?q??)?njo(?q?)
parenrightBig
S?q
??
?q,?q ? ?
1
2
parenleftBig
njo(?q?) +nko(?q??) + 1
parenrightBig
S?q
?,?q ??
?q
bracerightbigg
(3.23)
The temperature independent phonon scattering rate due to mass disorder is given by
equation 3.24
?(?q,i)mass =
2?
planckover2pi1
parenleftBig
planckover2pi1?(?q ,i)
2
parenrightBig
V0
8?3
3Nsummationtext
i=1
integraltext
BZ d
?q?.
bracketleftbigg
?
braceleftBig
planckover2pi1?(?q,i)?planckover2pi1?(?q?,i?)
bracerightBigsummationtext
?
g?|e??(?q,i).e?(?q?,i?)|2
bracketrightbigg (3.24)
Here g? = summationtext
k
fk(?)[1?mk(?)/?m(?)]2, fk(?) is the fraction of kth isotope of atom ?
that has mass mk(?), and ?m(?) is the average mass of atom with index ?. e?(?q,i) and
?(?q?,i?) are respectively the eigenvector and frequency of the phonon mode (?q,i).
3.3 Numerical Recipe for ?anh Calculation
We outline here the steps for calculating the phonon lifetimes (?anh) associated with
anharmonicity induced phonon scattering. The first step in the calculation of ?anh is
the calculation of the harmonic force constants matrix (?ij ) and third order lattice
37
anharmonicity tensor (Aijk ). This is achieved via our recently developed RSSFD algorithm
described in section 2.4. The harmonic force constants matrix is needed to determine
the phonon frequencies which are necessary inputs together with the third order lattice
anharmonicity tensor for the calculation of the 3-phonon coupling (TPC) terms (S?q
??
?q,?q ?
and S?q
?,?q ??
?q ) equations 3.19 and 3.20. A fine sampling grid for the calculation of the
phonon spectrum and TPC terms is adopted based on numerically converged results. Group
symmetry analysis is next performed to determine which of the ?q-points in the sampling
grid are irreducible. We save this in a file called Q?list. The number of irreducible q-points
based on 16?16?16 and 8?8?6 sampling grids for MgO and PV-MgSiO?3 calculation are
145 and 100 respectively. The phonon lifetime is then directly evaluated at these ?q-points.
To take advantage of parallel computing technology, we have divided our calculation of ?anh
into two steps that are embarrassingly parallel. First we calculate the TPC terms at the
irreducible ?q-points in the sampling grid are irreducible. We save this in a file called Qlist.
The phonon lifetime is then directly evaluated at these ?q-points, using our SMEA based
code TPC.f90. The inputs to this code are:
? List of the irreducible ?q-points in the sampling grid.
? third order lattice anharmonicity tensor
? phonon frequencies (?i )
? masses of constituent atoms
The next step is to read the precalculated TPC terms and evaluate the scattering rates
at given temperatures, enforcing the conservation of energy and momentum over a finer
sampling grid. This is accomplished by using a code named TauSMEA.f90. The inputs to
this code are:
? List of the irreducible ?q-points in the sampling grid.
38
? The number of temperature points
? The dimensions of the sampling grid adopted in TPC calculation
? The phonon frequencies calculated over the fine grid
? A finer grid parameter for enforcing energy and momentum conservations
? A file containing the temperature points
The phonon scattering rates in the whole Brillouin zone are then reconstructed using group
theory based on the crystal symmetry. We combine this with the temperature independent
scattering rates from mass disorder to obtain the overall phonon scattering rate.
39
Chapter 4
Lattice Thermal Conductivity of ??Al2O3 Crystal
4.1 Introduction
Alumina (Al2O3) is a well-known ceramic material with a wide range of technological
applications due to its outstanding mechanical and optical properties. Alumina constitutes
nearly 15% of minerals in the Earth?s crust and is thus a material of geophysical interest.
Its unique blend of low thermal expansion and high compressive strength makes it
suitable for use in thermal shock applications. Few chemicals attack alumina making it
a good candidate for chemical barrier coatings Alumina shows good electrical insulation at
high temperatures, good wear resistance and is very hard At ambient conditions, Al2O3
crystallizes in the corundum structure (??Al2O3, space group R?3C). Experimental88,122?124
and theoretical71,73,76,125,126 studies show that ??Al2O3 is stable up to 90 GPa. The space
group, formula units per primitive unit cell Z and Wyckoff sites of ??Al2O3 are summarized
in table 4.1.
Table 4.1: Space groups, formula units per primitive unit cell Z and Wyckoff sites of
??Al2O3:
Phase Space group Z Species Wyckoff site
? R?3c 2 O 6e
Al 4c
Thermal conductivity of Corundum-structured Al2O3 has been studied
experimentally127?129 however experimental and theoretical data on the pressure and
dependence of phonon lifetimes are lacking. We have developed a first-principles
computational technique to calculate the lattice thermal conductivity of solids at extreme
40
pressure and temperature conditions. Our technique has been successfully adopted to
study two-atom cubic MgO101,112. We present here our first-principles calculation of lattice
thermal conductivity of the trigonal corundum phase of Alumina (Al2O3) (Ten atoms per
unit cell). We directly evaluated the phonon-phonon scattering rates using our calculated
third order lattice anharmonicity tensors, within the framework of quantum scattering
theory. Combining this with our calculated vibrational phonon frequencies and group
velocities, within the framework of the Boltzmann transport theory we predict the lattice
thermal conductivity of ?? Al2O3 at six different densities and ten temperature points
ranging from 300K- 2500K, corresponding to a total of 60 T-P configurations without
any empirical fitting. The motivation of the current work is to test the applicability of
our newly developed computational technique on a more complex material system like
?-Al2O3 with 10 atoms per unit cell. To our knowledge this is the first ab initio study of
the volume/density dependence of individual phonon mode lifetimes, mode group velocity
and mode heat capacity. By combining this with the calculated thermal EOS of ?-Al2O3,
we derive the pressure dependence of lattice thermal conductivity in ?-Al2O3 depicted in
figures 4.12 and 4.13.
4.2 Lattice Dynamics and Equilibrium Thermal Properties
4.2.1 Symmetry Analysis and First-principles Calculation of Total Energies
Symmetry Analysis and Computational Time
To calculate the harmonic force constants and third order lattice anharmonicity tensors,
10-atom unit cells of volumes ranging from 75 -90 ?A3 were first relaxed to their ground state.
The independent interatomic force constants (?ij) and third order lattice anharmonicity
tensors (Aijk) for 120-atom supercells were then determined based on group symmetry
41
analysis. The results of the symmetry analysis to determine the independent terms is
summarised in the table 4.2.
Table 4.2: Irreducible atomic displacements (number of total energy calculations required),
Independent and dependent ?ij and Aijk terms from crystal symmetry.
Irreducible moves Independent terms Dependent terms ?
?ij 18 (ISAD) 1,441 63,348 0.06 ?A
Aijk 5,568 (IPAD) 163,267 7,668,135 0.06 ?A
After performing group symmetry analysis, we find that a total of 18 irreducible single
atom displacements (ISAD) and 5568 irreducible paired-atom displacement (IPAD) are
needed to determine the Hellmann-Feynman (HF) forces required to extract ?ij and Aijk
based on 120-atom periodic supercells. These numbers translate directly into number of
VASP114,130 calculations to be executed per atom. As you can see, computing the HF forces
for the extraction of third order lattice anharmonicity is a formidable task, demanding a
total of 5,568N (N = number of atoms in per supercell) calculations, with N computations
requiring ? 2,623s. This is equivalent to 4,058 computational hours or 169 days (5months)
if these computations were to be performed in series one after the other. Fortunately, this
step is highly parallelised. As a result, depending on the number of available compute nodes,
several computations can be performed concurrently. These computations were performed on
the Auburn University physics department?s PRISM research cluster, a distributed memory
computer system (Beowulf cluster). This cluster is comprised of 128 AMD Athlon MP CPUs.
With 128 CPUs running in parallel, the 5568 computations can be distributed across the
cluster with 44computations assigned per CPU. This dramatically reduces the computational
time from 5 months to 31 computer hours (? 1.5 days) and makes the task feasible. HF
forces were computed using the RSSFD method described in section 2.4.
42
Total Energy Calculations
Total energy calculations were performed within the framework of density functional
theory (DFT) with a plane wave basis set and plane augmented waves (PAW) method which
has been implemented in the VASP code114,131. We adopted the local density approximation
(LDA) to handle the exchange and correlation functional. Plane wave basis functions with
energies up to 400.00 eV were used for the PAW method. A totalenergy convergence criterion
of 10?9 eV was specified, for self-consistent iterations. Monkhorst-Park K-point sampling
grids of 6 ? 6 ? 4 were used for Brillouin zone integrations. The calculated static energy
at various unit cell volumes was fitted to 3rd order Birch-Murnaghan EOS. The static EOS
parameters are presented in table C.1 in appendix C, together with previous calculations
and experimental data for comparison.
Harmonic and Anharmonic Inter-atomic Potentials
Both harmonic and anharmonic inter-atomic potentials were obtained based on the HF
forces that were computed for all the irreducible distorted 120-atom supercell models using
a single Brillouin zone center q-point sampling grid. We next interpolated the calculated
independent ?ij and Aijk over the calculation volumes. Figures 4.1 and 4.3 show the
irreducible LDA calculated density dependent harmonic and anharmonic lattice dynamics of
Al2O3 based on a 120-atom super cell model. Harmonic force constants are obtained using
a irreducible single-atom displacement(ISAD) algorithm. Third order lattice anharmonicity
tensors are extracted using a irreducible paired-atom displacement (IPAD) scheme described
in section 2.4. Both results suggest a nearly linear density dependence of harmonic force
constants and third order lattice anharmonicity tensors. The volume dependence of the
inter-atomic potentials enables us to predict the inter-atomic potentials for any arbitrary
volume within our calculation rangewithout the need to performa separate VASP calculation
as explained in appendix B.
43
4.2.2 Quasi-harmonic Phonons
Phonon frequencies are obtained by solving equation 2.23. To account for the correction
to the dynamical matrix resulting from LO-TO splitting, we performed a separate calculation
according to Kunc and Martin?s method118. In a purely harmonic system, the inter-atomic
potentials are volume independent as well as the phonon frequencies. Figure 4.1 showing our
LDA calculated 2nd order inter-atomic force constants, clearly reveals a volume dependence.
in the limit of small deviations from the equilibrium volume, this anharmonic behavior
can be approximated via the quasi-harmonic approximation (QHA). Within the QHA, the
crystal is considered harmonic at fixed volumes, while the phonon frequencies become volume
dependent. This volume dependence is described via the mode Gr?uneisen parameter (? )
defined by equation 2.43 in section 2.5. In this calculation, the static energy was fitted to the
3rd-order Birch-Murnaghan equation of state (BM-EOS) and the thermal free energy was
fitted to the 2nd-order BM-EOS. F (T,V) = Estatic (V) + Fvib (T,V). Our calculated zero
pressure thermal expansion coefficient (TEC) of ?-Al2O3 as a function of temperature using
PAW method are compared with former measurements in figure C.2.
Figure 4.2 shows our calculated vibrational phonon frequencies, density of states and
group velocities at 0 GPa and 300K. Experimental data (discrete symbols) is included
for comparison. The overall agreements between our LDA + PAW calculated vibrational
phonon frequencies and experiment, are good and within the typical accuracy of ab initio
calculations. Among the calculated 30 phonon branches, the low-frequency (acoustic) ones
have steeper slopes implying larger phonon mode group velocities (Vg = d?dK). This suggests
that acoustic phonon modes are more effective carriers of heat compared with the flat optical
modes.
We believe that our LDA calculated phonon dispersions are reliable based on the
excellent agreement with measurements from experiment as shown in figure 4.2. The
calculated lattice anharmonicity tensors on the other hand are not known to have measured
44
s55s53s46s48 s55s55s46s53 s56s48s46s48 s56s50s46s53 s56s53s46s48
s45s49s48s46s48
s45s53s46s48
s48s46s48
s53s46s48
s49s48s46s48
s49s53s46s48
s50s48s46s48
s50s53s46s48
s51s48s46s48
s51s53s46s48
s105
s106
s32
s40
s101
s86
s47
s197
s50
s41
s86s111s108s117s109s101s40s197
s51
s47s117s110s105s116s32s99s101s108s108s41
s86
s48
Figure 4.1: LDA calculated density dependent harmonic force constants of Al2O3 based on
a 120-atom super cell model. A total of 1441 independent ?ij were extracted from the HF
forces from VASP calculation. Majority of the ?ij terms are small and insensitive to volume
change. Less than 20 inter-atomic force constants have magnitudes exceeding 5 eV/?A2 and
show an almost linearly increasing pressure dependence with increasing pressure.
45
s48
s50s48s48
s52s48s48
s54s48s48
s56s48s48
s49s48s48s48
s32s81s72s65s40s84s104s105s115s32s119s111s114s107s41
s69s120s112s116s32s40s83s99s104s111s98s101s114s32s101s116s32s97s108s41
s40s98s41s86s68s79s83
s32
s32
s70
s114
s101
s113
s117
s101
s110
s99
s121
s32
s40
s99
s109
s32
s45
s49
s41
s40s97s41s80s104s111s110s111s110s32s68s105s115s112s101s114s115s105s111s110
s90
s68s65 s65s114s98s105s116s114s97s114s121s32s117s110s105s116s115
s32
s32
s32
s32
s48
s50
s52
s54
s56
s49s48
s49s50
s32
s32
s40s99s41s71s114s111s117s112s32s86s101s108s111s99s105s116s121s40s75s109s47s115s41
s90
s68s65
s48
s50
s52
s54
s56
s49s48
s49s50
Figure 4.2: Vibrational phonon frequencies and density of states. Our LDA calculated
phonon spectrum (solid lines) agree well with measurements from experiment (discrete
symbols). Low frequency(acoustic) phonon modes have steeper slopes (group volicities)
compared to the high frequency (optical) modes which are mostly flat.
46
s55s53s46s48 s55s55s46s53 s56s48s46s48 s56s50s46s53 s56s53s46s48
s45s54s48
s45s52s48
s45s50s48
s48
s50s48
s52s48
s54s48
s65
s105
s106
s107
s40
s101
s86
s47
s197
s51
s41
s86s111s108s117s109s101s47s117s110s105s116s32s99s101s108s108s40s197
s51
s41
s45s54s48
s45s52s48
s45s50s48
s48
s50s48
s52s48
s54s48
Figure 4.3: Pressure dependence of anharmonicity tensor elements (Aijk). Aijk elements are
observed to vary linearly with volume/density. From 0 GPa (? 84 ?A3) to 50 GPa ( ? 75
?A3), Aijk increases linearly from 40 eV/?A3 to 60 eV/?A3, representing a 50% increase, ? 1%
per GPa.
counterparts. To verify their robustness, we evaluate mode Gr?uneisen parameters using the
extracted lattice anharmonicity and compared with the ones evaluated from the already
verified harmonic force constants. The comparison of the mode Gr?uneisen parameters is
shown in figure 4.4. Both mode Gr?uneisen dispersion plots capture the same features,
showing that our calculated lattice anharmonicity tensors are indeed reliable. The process
of extracting the mode Gr?uneisen parameters from ?ij and Aijk is described in section 2.5
47
s48
s49
s50
s32
s32
s77
s111
s100
s101
s32
s71
s114
s252
s110
s101
s105
s115
s101
s105
s110
s32
s80
s97
s114
s97
s109
s101
s116
s101
s114
s90
s68
s65
s32s68s101s114s105s118s101s100s32s102s114s111s109s32s97s110s104s97s114s109s111s110s105s99s105s116s121s32s116s101s110s115s111s114
s32s68s101s114s105s118s101s100s32s102s114s111s109s32s112s104s111s110s111s110s32s115s112s101s99s116s114s117s109
Figure 4.4: Mode Gr?uneisen Parameters calculated from the force constants matrix(blue
circles) and from the third order anharmonicity tensors(orange circles) compared . We
use this information to independently test the reliability of our calculated third order
anharmonicity tensors. The extraction of Gr?uneisen Parameters from third order lattice
anharmonicity is described in section 2.5
48
78 79 80 81 82 83
Volume/Unit Cell
5.0?10-8
1.0?10-7
1.5?10-7
2.0?10-7
Mode heat Capacity
Figure 4.5: Volume dependence of mode heat capacity. Lines represent the heat capacity of
individual phonon modes at the sampled at the irreducible K-points in Brillouin zone at 6
different densities/unit cell volumes.
78 79 80 81 82 83
Volume/Unit Cell
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
Vg
2 (10
8 Km
2 S
-2 )
Figure 4.6: Volume/density dependence of square of phonon mode group velocities. Vg2
shows a mild volume dependence. Remember that according to knetic transport theory,
lattice thermal conductivity is given by parenleftbig? = 13summationtextCvvg2?parenrightbig
49
4.3 Phonon Lifetimes
Phonon lifetimes are evaluated based on an 8?8?8 k-point grid for a total of 512
k-points, 150 of which are irreducible. Using a 10-atom unit cell, the life times of the 30
phonon modes are evaluated at each k-point using Fermis golden rule (equation 3.16). The
full life times are obtained using the symmetry relations between the irreducible k-points
and the dependent ones. Figure 4.7 shows the lifetimes of the phonon modes at all the 150
irreducible k-points in Brillouin zone.
s55s56 s55s57 s56s48 s56s49 s56s50 s56s51 s56s52 s56s53
s48
s53
s49s48
s49s53
s50s48
s50s53
s51s48
s51s53
s52s48
s100s41s32s79s112s116s105s99s97s108s32s109s111s100s101s115
s99s41s32s76s65
s32
s32
s40
s112
s115
s41
s97s41s32s84s65s49
s55s56 s55s57 s56s48 s56s49 s56s50 s56s51 s56s52 s56s53
s48
s49s48
s50s48
s51s48
s52s48
s98s41s32s84s65s50
s32
s32
s55s56 s55s57 s56s48 s56s49 s56s50 s56s51 s56s52 s56s53
s48
s53
s49s48
s49s53
s50s48
s50s53
s51s48
s51s53
s52s48
s32
s32
s40
s112
s115
s41
s86s111s108s117s109s101s47s85s110s105s116s32s67s101s108s108s32s40s197
s51
s41
s55s56 s55s57 s56s48 s56s49 s56s50 s56s51 s56s52 s56s53
s48
s53
s49s48
s49s53
s50s48
s50s53
s51s48
s51s53
s52s48
s32
s32
s86s111s108s117s109s101s47s85s110s105s116s32s67s101s108s108s32s40s197
s51
s41
Figure 4.7: Density Dependence of phonon lifetimes in ?-Al2O3 the at 300K (a) & (b)
Transverse acoustic modes (c) Longitudinal acoustic phonon mode; (d)Optical phonon
modes. Lifetimes of optical phonon modes are almost pressure insensitive. Acoustic modes
increase almost linearly with decreasing volume (increasing pressure). The optical phonon
mode with the longest lifetime is only about half the lifetime of the longest acoustic mode.
This singles out acoustic phonons as the dominant carriers of heat in ?-Al2O3.
50
4.4 Lattice Thermal Conductivity ?
The inverse of thermal conductivity of Al2O3 scales linearly with temperature at
constant density, figure 4.10. The thermal conductivity shown, includes the contributions
from anharmonicity induced phonon scattering and isotope-induced phonon scattering up to
lower mantle temperatures. Our calculated ?(T) at ambient pressure figure 4.12, is in good
agreement with measurements from experiment127?129. At 0 GPa and 300 K we predict ?
as ? 23 Wm?1K?1 at and 56 Wm?1K?1 at 50 GPa and 300 K, representing a factor of 2.4
increase in ? per GPa
78 79 80 81 82 83
Volume/Unit Cell
0.0
50.0
100.0
150.0
200.0
250.0
300.0
Mode Thermal Conductivity
Figure 4.8: Density dependence of mode-? from acoustic phonon modes at T= 300K.
Acoustic modes thermal conductivity follows a volume dependence identical to that of the
acoustic modes phonon lifetimes shown in figure 4.7 a, b & c. This is expected as the other
determining variables in the thermal conductivity equation Cv and V2g, figures 4.5 and 4.6
are almost flat(volume insensitive)
51
78 79 80 81 82 83
Volume /Unit Cell
0.0
5.0
10.0
15.0
20.0
25.0
30.0
Mode Thermal Conductivity
Figure 4.9: Density dependence of mode- thermal conductivity from optical phonon modes
at T= 300K. ?optical is almost pressure insensitive, consistent with the pressure dependence
of ?optical, figure 4.7 d, and Cv and V2g, (figures 4.5 and 4.6)
s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48
s48s46s48s48
s48s46s49s48
s48s46s50s48
s48s46s51s48
s48s46s52s48
s32
s32
s61s32s52s46s53s49s32s103s47s99s109
s51
s61s32s52s46s51s55s32s103s47s99s109
s51
s61s32s52s46s50s51s32s103s47s99s109
s51
s61s32s52s46s49s48s32s103s47s99s109
s51
s61s32s52s46s48s53s32s103s47s99s109
s51
s40
s87
s45
s49
s46s109
s46s75
s41
s84s101s109s112s101s114s97s116s117s114s101s40s107s41
s48s46s48s48
s48s46s49s48
s48s46s50s48
s48s46s51s48
s48s46s52s48
s61s32s51s46s57s56s32s103s47s99s109
s51
Figure 4.10: Density dependence of ??1. We have directly calculated the lattice thermal
conductivity of ??Al2O3 at six densities ranging from 3.98 gcc?1 to 4.51 gcc?1 and at 10
different temperature pointsfrom300Kto 2500K.??1 shows a linear temperature dependence
at constant density. At constant temperature ??1 increases almost linearly with increasing
density.
52
s48 s49s48 s50s48 s51s48 s52s48 s53s48
s48
s49s48
s50s48
s51s48
s52s48
s53s48
s54s48
s55s48
s48
s49s48
s50s48
s51s48
s52s48
s53s48
s54s48
s55s48
s32
s49s32s32s61s32s55s55s46s53s48s32s197
s51
s32s40 s32s61s32s52s46s53s49s32s103s47s99s109
s51
s41
s50s32s32s61s32s32s56s48s46s48s48s32s197
s51
s32s40 s32s61s32s52s46s51s55s32s103s47s99s109
s51
s41
s51s32s32s61s32s32s56s50s46s53s48s32s197
s51
s32s40 s32s61s32s52s46s50s51s32s103s47s99s109
s51
s41
s52s32s32s61s32s32s56s51s46s53s54s32s197
s51
s32s40 s32s61s32s52s46s49s48s32s103s47s99s109
s51
s41
s53s32s32s61s32s32s56s53s46s48s48s32s197
s51s32
s40 s32s32s61s32s52s46s48s53s32s103s47s99s109
s51
s41
s54s32s32s61s32s32s56s55s46s53s48s32s32s197
s51
s40 s61s32s51s46s57s56s32s103s47s99s109
s51
s41
s51s48s48s75
s49s53s48s48s75
s49s48s48s48s75
s50s53s48s48s75
s40
s87
s109
s45
s49
s75
s45
s49
s41
s80s114s101s115s115s117s114s101s32s40s71s80s97s41
s53s48s48s75
s45s65s108
s50
s79
s51
s85s110s105s116s32s67s101s108s108s32s86s111s108s117s109s101s32s40s68s101s110s115s105s116s121s41
s50s48s48s48s75
Figure 4.11: Pressure dependence of isothermal ?. Based onthe V(T,T) dependence fromour
LDAcalculated thermal EOS we have derived the pressure dependence of?. The numbers 1-6
represent the six unit cell volumes (densities) for which ? was calculated. ? isotherms show a
linear pressure dependence. From this pressure dependence, we can derive the temperature
dependence of isobaric ? shown in figure 4.12. From the labels 1-6 one can preview the
temperature dependence of ? at conditions of constant pressure.
s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48
s48
s49s48
s50s48
s51s48
s52s48
s53s48
s54s48
s69s120s112s116s58s32s75s105s110s103s101s114s121s32s40s49s57s53s52s41
s32s69s120s112s116s58s32s75s105s110s103s101s114s121s32s40s49s57s53s55s41s32s80s111s108s121s32s99s114s121s115s116s97s108
s32s69s120s112s116s58s32s75s105s110s103s101s114s121s32s40s49s57s53s55s41s32s83s105s110s103s108s101s32s99s114s121s115s116s97s108
s32s69s120s112s116s58s32s83s108s97s99s107s40s49s57s54s50s41
s32s69s120s112s116s58s32s86s105s115s104s110s101s118s115s107s105s32s40s49s57s54s56s41
s32s84s104s101s111s114s121s58s32s48s32s71s80s97
s32s84s104s101s111s114s121s58s32s53s48s32s71s80s97
s40
s87
s109
s45
s49
s75
s45
s49
s41
s84s101s109s112s101s114s97s116s117s114s101s40s75s41
s48 s49s48 s50s48 s51s48 s52s48 s53s48
s49s46s48
s49s46s53
s50s46s48
s50s46s53
s51s46s48
s51s46s53
s52s46s48
s52s46s53
s49s46s48
s49s46s53
s50s46s48
s50s46s53
s51s46s48
s51s46s53
s52s46s48
s52s46s53
s32s51s48s48s75
s32s49s48s48s48s75
s32s49s53s48s48s75
s32s50s48s48s48s75
s32s50s53s48s48s75
s80s32s40s71s80s97s41
Figure 4.12: Isobaric ? of ?-Al2O3 at 0 GPa and 50 GPa. We show a comparison between
our calculated thermal conductivity and experiment at 0GPa. Discrete symbols denote
experiment data while solid continuous line represent our first-principles calculation.
53
s48 s49s53 s51s48 s52s53 s54s48 s55s53
s49s46s48
s49s46s53
s50s46s48
s50s46s53
s51s46s48
s51s46s53
s52s46s48
s52s46s53
s53s46s48
s32s51s48s48s75
s32s50s53s48s48s75
s32
s32
s100
s108
s110
s47
s100
s80
s32
s40
s37s47
s71s80
s97
s41
s80s114s101s115s115s117s114s101s40s71s80s97s41
Figure 4.13: Pressure dependence of ?. At 0 GPa, we predict nearly 3 %/GPa and 4.6
%/GPa increase in ? at 300 K and 2500K respectively.
Table 4.3: Calculated lattice thermal conductivity fitted to second order in pressure at 10
different temperatures: ? = a0+ a1P + a2P2. The parameters a0,a1,a2 are then used
in conjunction with the calculated thermal EOS to predict the isobaric lattice thermal
conductivity at different temperatures as shown in figure 4.12
T a0 (Wm?1K?1) a1 (Wm?1K?1GPa?1) a2 (Wm?1K?1GPa?2)
300K 22.544 0.74449 0.00004956
500K 12.372 0.40730 0.00034499
750K 7.6949 0.27119 0.00041005
1000K 5.4235 0.20468 0.00035945
1250K 4.0608 0.16568 0.00032107
1500K 3.1491 0.13990 0.00028940
1750K 2.4901 0.12191 0.00026943
2000K 1.9958 0.10818 0.00024826
2250K 1.6059 0.097737 0.00023343
2500K 1.2890 0.089519 0.00022246
54
4.5 Conclusion
The density dependent harmonic and anharmonic inter-atomic potentials in ?-Al2O3
are calculated using our supercell finite-displacement algorithm. Phonon frequencies and
the thermal equation of state derived from the harmonic force constants are in excellent
agreement with experiment. Our LDA calculated 3rd order lattice anharmonicity has been
independently tested to be reliable via the mode Gr?uneisen parameters. Phonon life times
for individual phonon modes are directly evaluated based on quantum phonon scattering
theory. Acoustic phonon modes are found to have noticeably longer lifetimes compared to
the optical modes. Moreover the lifetimes of the acoustic modes increase almost linearly with
pressure while the optical modes are almost insensitive to pressure change. In ?-Al2O3 the
acoustic phonon modes contribute about 60% of the total thermal conductivity at ambient
conditions. Lattice thermal conductivity of 60 T-P conditions are derived without any
empirical extrapolation for a material with 10-atoms per unit cell. Our predicted ? at
ambient conditions is about 23 Wm?1K?1 , in good agreement with available experiment.
dln?
dP , the pressure dependence of ? is around 3 and 5% per GPa at 0GPa for 300K and 2500K
respectively.
55
Chapter 5
LATTICE THERMAL CONDUCTIVITY OF PV-MgSiO3 AT LOWER
MANTLE CONDITIONS
5.1 Introduction
MgSiO3 is known to exist in two phases, the perovskite(PV) and post-perovskite (PPV)
phases. The PV ? PPV transition is reported to occur at 125GPa and 2500K132. The
study of MgSiO3 has attracted much attention since the discovery of the PV ? PPV
transition which is believed to explain the seismic discontinuities observed in the Earth?s D??
layer72. The raman spectrum of PV and PPV-MgSiO3 has been studied theoretically132 and
experimentally133. PV-MgSiO3 has an orthorhombic structure with the symmetry group
#62 Pnma and is believed to be thermodynamically stable up to 110 GPa at 0K72,132.
PPV-MgSiO3 the high pressure phase of MgSiO3 belongs to the Cmcm symmetry group and
has 2 formula units per unit cell. The Space groups, formula units per primitive unit cell Z
and Wyckoff sites for the two phases are summarised in table 5.1
Table 5.1: Space groups, formula units per primitive unit cell Z and Wyckoff sites of 2
polymorphs of MgSiO3: perovskite and post-perovskite.
Phase Space group Z Species Wyckoff site
Perovskite Pnma 4 Mg 4a
Si 4c
O 4c, 8d
Post-perovskite Cmcm 2 Mg 4a
Si 4c
O 4c, 8f
56
Table 5.2: Irreducible representation of the mode symmetries at ?-point for the 2 polymorphs
of MgSiO3: perovskite and post-perovskite.
Phase Raman active IR active Silent
Perovskite 7Ag + 5B1g + 7B2g + 5B3g 10B1u + 8B2u + 10B3u 8Au
Post-perovskite 4Ag + 3B1g + 1B2g 6B1u + 6B2u + 4B3u 2Au
Based on group symmetry, the sixty (3 acoustic + 57 optical) in PV and 30 (3+27)
phonon modes have the following irreducible representation at the ?- point summarised in
table 5.2.
The Earths lower mantle is mostly composed of magnesium silicate in the perovskite
structure (PV-MgSiO3. A recently published measurement of the thermal conductivity of
MgSiO3 perovskite134 disagrees significantly with the only other existent measurement107
and the measured pressure dependence implies values of 9-12 W/m/K when extrapolated
to the core/mantle boundary. The disagreement between the two measured values and the
uncertainty in the long extrapolations required to estimate the value of silicate perovskite?s
thermal conductivity at the pressure and temperature conditions of the lower mantle
motivates our first-principles calculations of the phonon behavior of perovskite which governs
this materials thermal conductivity at the conditions of the lower mantle. The value of
lower mantle thermal conductivity is important to to modeling the total budget inside the
earth. Large uncertainties exist in our current estimates of lower mantle thermal conductivity
because of the long extrapolations in both temperature (T) and pressure (P) of laboratory
measured thermal conductivity of lower mantle minerals at relatively low temperatures and
pressures.
57
5.2 First-principles Calculation of Total Energies, Vibrational Properties and
Phonon Lifetimes
The VASP software package114,130 was used to calculate the relaxed equilibrium atomic
structures, total energies, force constants, and 3rd order lattice anharmonicity tensors using
Bl?och?s projector augmented wave (PAW) pseudopotential135,136 and a plane-wave (PW)
basis set within the local density approximation (LDA) to the density functional theory
(DFT). The details of LDA-PW calculations are as follows: For the PV-MgSiO3 phase,
structures of 20-atom (10-atom for PPV phase) unit cells of different volumes were first
relaxed to their ground state. Independent force constants for 160-atom periodic supercells
were then identified by the crystal symmetry. Both harmonic and anharmonic force constants
were obtained based on the atomic forces that were computed for all the irreducible distorted
supercells using a single Brillouin zone (BZ) center q-point sampling grid. Next, for the PV
phase only, the lifetime of each phonon mode was calculated based on the Fermis golden rule
equation 3.16. This is the most computational intensive step because all the configurations
of three-phonon scattering need to be properly counted. All the phonon modes, acoustic or
optic, are calculated using the same numerical approach. No empirical approximations of
the frequency or q-point dependence are assumed. In the current study of PV-MgSiO3, the
BZ integration was approximated by a summation over a discrete 8 ? 8?6 q-grid. A finer
64?64?48 q-grid was used to interpolate the phonon energies (frequencies) based on the
chosen 8 ? 8 ? 6 q-grid to accurately account for the energy and momentum conservation.
The lattice symmetry of the crystals was analyzed, and only the phonon lifetimes of an
irreducible set of phonon modes are calculated directly and the rest are reconstructed using
the group theory.
58
5.3 Results
5.3.1 Lattice Dynamics and Related Thermal Properties
Both harmonic and anharmonic force constants were obtained based on the HF forces
that were computed for all the irreducible distorted supercells using a single Brillouin zone
(BZ) centered q-point sampling grid. show the independent inter-atomic force constants in
PV and PPV respectively. Independent force constants shown in figures 5.1 and D.4 for
160-atom periodic supercells were identified by the crystal symmetry and obtained based on
the Hellmann-Feynman forces that were computed for all the irreducible distorted supercells
using a single Brillouin zone (BZ) centered q-point sampling grid. An atomic displacement
of ? = 0.06 ?A was adopted in our ISAD scheme
Table 5.3: Simulation cell volumes and corresponding pressures for PV-MgSiO3 calculation
T=300 K.
V0 (?A3) 122 127 137 147 157
P (GPa) 119.684 95.264 57.100 29.45 9.211
120 140 160 180Volume/Unit Cell
-20
0
20
40
60
80
Harmonic Force Constants
Figure 5.1: Density dependence of harmonic force constants in PV-MgSiO3. Independent
force constants shown for 160-atom periodic supercells of PV.
59
s48
s49s53s48
s51s48s48
s52s53s48
s54s48s48
s55s53s48
s57s48s48
s49s48s53s48
s80
s104s111s110s111s110s32s70
s114
s101
s113s117s101
s110s99
s105
s101
s115
s32s40
s99
s109
s45
s49
s41
s32
s32
s80
s104
s111
s110
s111
s110
s32
s70
s114
s101
s113
s117
s101
s110
s99
s105
s101
s115
s32
s40
s99
s109
s45
s49
s41
s80s86s45s77s103s83s105s79
s51
s50s48s45s97s116s111s109s32s117s110s105s116s99s101s108s108
s80s32s61s32s48s32s71s80s97
s48
s49s53s48
s51s48s48
s52s53s48
s54s48s48
s55s53s48
s32
s88
s75 s76
s90 s88 s83 s82 s89
s32s32s32s32s32s77s103s79
s50s45s97s116s111s109s32s117s110s105s116s32s99s101s108s108
Figure 5.2: LDA calculated phonon mode frequencies in PV-MgSiO3 and MgO at 0 GPa.
s48
s50s53s48
s53s48s48
s55s53s48
s49s48s48s48
s49s50s53s48
s49s53s48s48
s80
s104
s111
s110
s111
s110
s32
s32
s70
s114
s101
s113
s117
s101
s110
s99
s105
s101
s115
s32
s40
s99
s109
s45
s49
s41
s32
s90 s88 s83 s82 s89
s49s50s48s32s71s80s97
s32
s80s86s32s45s32s40s77s103s83s105s79
s51
s41s32s76s68s65s32s67s97s108s99s117s108s97s116s105s111s110
Figure 5.3: LDA calculated phonon mode frequencies in PV-MgSiO3 at 120 GPa
60
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s90 s88 s83 s82 s89
s32s32
s71
s114
s111
s117
s112
s32
s86
s101
s108
s111
s99
s105
s116
s121
s32
s40
s75
s109
s47
s115
s41
s32
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s32
s32s32
s71
s114
s111
s117
s112
s32
s86
s101
s108
s111
s99
s105
s116
s121
s32
s40
s75
s109
s47
s115
s41
s32
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s90 s88 s83 s82 s89s90 s88 s83 s82 s89
s90 s88 s83 s82 s89
s32
s32s32
s71
s114
s111
s117
s112
s32
s86
s101
s108
s111
s99
s105
s116
s121
s32
s40
s75
s109
s47
s115
s41
s32
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s100s41s32 s32s61s32s52s46s49s52s49s32s103s47s99s109
s51
s32
s32s32
s71
s114
s111
s117
s112
s32
s86
s101
s108
s111
s99
s105
s116
s121
s32
s40
s75
s109
s47
s115
s41
s32
s97s41s32 s32s61s32s53s46s52s54s53s32s103s47s99s109
s51
s98s41s32 s32s61s32s52s46s56s54s55s32s103s47s99s109
s51
s99s41s32 s32s61s32s52s46s53s51s54s32s103s47s99s109
s51
Figure 5.4: LDA calculated phonon mode group velocities in PV-MgSiO3 at four different
densities shown. Group velocities measured along the Z-direction are more sensitive to
pressure. From a density of 5.465 g/cm3 (? 120 GPa) to 4.141g/cm3 (? 0 GPa)the phonon
group velocity in the ??Z direction decreases from around 15 km/s to 11 Km/s. Velocities
measured in the other high symmetry directions X and Y are almost insensitive to pressure.
120 140 160Volume/Unit Cell-200
-100
0
100
200
A ijk
Figure 5.5: Density dependence of 3rd order lattice anharmonicity tensor (Aijk) elements
in PV-MgSiO3. Most Aijk elements vary nearly linearly with volume/density. Most of the
elements are very small.
61
s48
s49
s50
s51
s52
s77
s111
s100
s101
s32
s71
s114
s252
s110
s101
s105
s115
s101
s110
s32
s80
s97
s114
s97
s109
s101
s116
s101
s114
s40
s103
s41
s32
s90
s88
s83 s82 s89
s32s103s95s102s111s114s99s101s32s99s111s110s115s116s97s110s116s115
s32s103s95s97s110s104s97s114s109s111s110s105s99s105s116s121
s48s71s80s97
s32
s80s86s32s45s32s40s77s103s83s105s79
s51
s41s32s32s76s68s65s32s67s97s108s99s117s108s97s116s105s111s110
Figure 5.6: Mode Gr?uneisein parameters in PV-MgSiO3. Blue circles represent the Mode
Gr?uneisein parameters calculated from the second order force constants while the orange
circles show the same parameters calculated via our LDA calculated third order lattice
anharmonicity tensors. This step is used to independently test the reliability of our lattice
anharmonicity data since Aijk is not known to have measured counterparts to which we can
compare our calculated data.
62
5.3.2 Phonon Lifetimes
The lifetime of each phonon mode was calculated based on the Fermis golden rule
equation 3.16. All the phonon modes, acoustic or optic, are calculated using the same
numerical approach described in section 5.2. No empirical approximations of the frequency
or q-point dependence are assumed. Figures 5.7 show the volume/density dependence of
phonon lifetimes of all 60 phonon modes in PV-MgSiO3 sampled at three high symmetry
points of the Brillouin Zone: X (?a,0,0), Y (0, ?b,0) and Z (0,0, ?c). The volume dependence
of the 60 modes at all the 100 irreducible directions (K-points) of the BZ is reported in
appendix E in figures E.1 to E.17. Our calculation reveals a weak pressure dependence of
phonon lifetimes. We also observe an almost linear volume dependence.
s49s50s53 s49s51s48 s49s51s53 s49s52s48 s49s52s53 s49s53s48 s49s53s53
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s49s56
s50s48
s90s32s32s80s111s105s110s116s58s32s40s48s44s48s44 s99s41 s88s32s32s80s111s105s110s116s58s32s40 s97s44s48s44s48s41
s97s41s32s32s32
s32
s32
s40
s112
s115
s41
s98s41s32
s49s50s53 s49s51s48 s49s51s53 s49s52s48 s49s52s53 s49s53s48 s49s53s53
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s49s56
s50s48
s86s111s108s117s109s101s47s85s110s105s116s32s67s101s108s108s32s40s197
s51
s41
s32
s32
Figure 5.7: Volume/density dependence of phonon lifetimes in PV-MgSiO3 all 60 modes at
two high symmetry points: X (?a,0,0) and Z (0,0, ?c) in the Brillouin Zone
63
5.3.3 Pressure Dependence of Thermal Conductivity
To predict ?Lattice from ambient condition to the pressure-temperature conditions of
CMB without extrapolation, we repeated our calculations for 4 density configurations from
(5.465 gcm?3 to 4.141 gcm?3) at 9 temperature conditions(from 300K to 4000K). The
pressures of the studied 36 density-temperature configurations (as listed in Table below are
then derived based onthe thermal equations of state (EOS) predicted from the first-principles
quasi-harmonic theory (QHA)117. Our calculated ?PV at ambient conditions (300K, 0GPa),
are 5.70 Wm?1K?1and 5.46 Wm?1K?1 respectively for isotopically pure crystals and natural
crystals that contain isotope mass disorder, and are in close agreement with the measurement
of 5.1 Wm?1K?1 by Osako134.
s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48 s51s48s48s48 s51s53s48s48 s52s48s48s48
s48s46s48
s48s46s50
s48s46s52
s48s46s54
s48s46s56
s49s46s48
s49s46s50
s49s46s52
s49s46s54
s49s46s56
s50s46s48
s50s46s50
s50s46s52
s32s49s50s50s32s197
s51
s32
s32
s40 s61s32s53s46s52s54s53s32s103s47s99s109
s51
s41
s32s49s51s55s32s197
s51
s32s40 s32s61s32s52s46s56s54s55s32s103s47s99s109
s51
s41
s32s49s52s55s32s197
s51
s32s40 s32s61s32s52s46s53s51s54s32s103s47s99s109
s51
s41
s32s49s53s55s32s197
s51
s32s40 s32s61s32s52s46s49s52s49s32s103s47s99s109
s51
s41
s40
s87
s45
s49
s109
s75
s41
s84s101s109s112s101s114s97s116s117s114s101s40s75s41
s85s110s105s116s32s99s101s108s108s32s86s111s108s117s109s101s32s40s68s101s110s115s105s116s121s41s32
Figure 5.8: Temperature dependence of ??1 in PV-MgSiO3. ? is calculated at 4 distinct
densities corresponding to unit cell volumes shown on the figure, at 9 different temperatures
ranging from 300 K to 4000K, for a total of 36 density-temperature configurations. ? is seen
to follow the T?1 dependence typical of lattice thermal conductivity in insulating materials.
At the bottom of the mantle close to the core/mantle boundary, a phase transformation
from the PV to the PPV structure may occur137?139 . To assess whether this transformation
64
s48 s49s53 s51s48 s52s53 s54s48 s55s53 s57s48 s49s48s53 s49s50s48 s49s51s53 s49s53s48
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s32s32s61s32s32s49s50s50s32s197
s51
s32
s32
s40 s61s32s53s46s52s54s53s32s103s47s99s109
s51
s41
s50s32s32s61s32s32s49s51s55s32s197
s51
s32s40 s32s61s32s52s46s56s54s55s32s103s47s99s109
s51
s41
s51s32s32s61s32s32s49s52s55s32s197
s51
s32s40 s32s61s32s52s46s53s51s54s32s103s47s99s109
s51
s41
s52s32s32s61s32s32s49s53s55s32s197
s51
s32s40 s32s61s32s52s46s49s52s49s32s103s47s99s109
s51
s41
s52s48s48s48s75
s51s48s48s48s75
s50s48s48s48s75
s49s48s48s48s75
s51s48s48s75
s40
s87
s109
s45
s49
s75
s45
s49
s41
s80s114s101s115s115s117s114s101s32s40s71s80s97s41
s53s48s48s75
s85s110s105s116s32s67s101s108s108s32s86s111s108s117s109s101s32s40s68s101s110s115s105s116s121s41
Figure 5.9: Pressure dependence of ??1 in PV-MgSiO3 at isothermal conditions. This is
derived based on the P-V-T relations obtained from our LDA calculated thermal EOS.
may influence thermal conductivity in the core/mantle boundary layer, we performed a single
calculation of thermal conductivity of a PPV-structured MgSiO3 at T=3000 K, P =130 GPa,
which yields similar values to the MgSiO3 perovskite structure. Based on our calculations
and the materials shared low symmetries, we expect that the difference in the intrinsic ?latt
of PV and PPV will be negligible.
65
s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48 s51s48s48s48 s51s53s48s48 s52s48s48s48
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s32s48s71s80s97
s32s53s48s71s80s97
s32s49s48s48s71s80s97
s32s49s51s53s71s80s97
s84
s104
s101
s114
s109
s97
s108
s32
s67
s111
s110
s100
s117
s99
s116
s105
s118
s105
s116
s121
s40
s87
s109
s45
s49
s75
s45
s49
s41
s84s101s109s112s101s114s97s116s117s114s101s40s75s41
s48 s49s53 s51s48 s52s53 s54s48 s55s53 s57s48 s49s48s53 s49s50s48 s49s51s53
s49
s50
s51
s52
s53
s54
s55
s56
s57
s49s48
s49
s50
s51
s52
s53
s54
s55
s56
s57
s49s48
s32s51s48s48s75
s32s53s48s48s75
s49s48s48s48s75
s32s49s53s48s48s75
s32s50s48s48s48s75
s32s50s53s48s48s75
s32s51s48s48s48s75
s32s51s53s48s48s75
s32s52s48s48s48s75
s48
s80s40s71s80s97s41
Figure 5.10: Main: Temperature dependence of isobaric ? in PV-MgSiO3. Inset: Pressure
variation of normalized lattice thermal conductivity. Using the coefficients from ? = a0 +
a1P +a2P2 the isobaric lattice thermal conductivity was derived. Four isobars 0, 50, 100 and
135 GPa are shown. Inset: Pressure variation of normalized ?. ?0 represents the thermal
conductivity at 0 GPa for each isotherm. From 300 K to 4000K, ?/?0 increases from 3 to 9
at 135 GPa.
66
5.4 Discussion
5.4.1 Comparison of the Behavior of ?MgO and ?MgSiO3
Both ?MgO and ?MgSiO3 follow a T?1 dependence at the isochoric condition as expected
for insulator thermal conductivity, but the pressure increase of ?MgSiO3 is substantially
weaker than that of ?MgO . For example, dln(?)dP at 300 K is 0.25 %GPa?1 for PV, while
the same value is ? 4 % GPa?1 for MgO. This weak pressure-dependence has a strong
effect at core/mantle conditions (3000 K, 135 GPa) where ?MgSiO3 is calculated to be ?
1.64 W/m/K, an order of magnitude smaller than most previous estimates for MgSiO3140,
and smaller than the calculated value for lattice thermal conductivity of MgO at similar
conditions, ? 40 W/m/K101. The very different lattice thermal conductivities calculated
for MgO and MgSiO3 and their different pressure-dependencies can be understood based
on their contrasting phonon behavior interpreted in the context of the kinetic transport
theory. According to the Peierls-Boltzmann equation, the bulk lattice thermal conductivity
is proportional to the phonon heat capacity, group velocity and lifetime. Of these, the
heat capacity plays a lesser role, especially at high temperatures where the heat capacity
of each phonon approaches the Dulong-Petit limit 3kB, making little contribution to the
differences between thermal conductivity of the two systems. The phonon dispersion curves
for the two crystals are also plotted in Figure 5.2. Note that 3 out every 6 phonon modes
(50%) are acoustic modes in MgO, while only 3 out of 60 modes (5%) are acoustic modes
in PV-MgSiO3. Acoustic phonon modes, especially those of q-points near the BZ center,
behave like classical elastic waves and they usually have large phonon group velocities. In
contrast, optic phonon modes can not be approximated with the classical elastic waves
and they usually show a much flatter dispersion. In both crystals, low-frequency acoustic
phonons (e.g. less than 200 cm?1) are effective heat carriers, while high-frequency optic
phonons (e.g. above 700 cm?1) are poor heat carriers. The drastic difference in ?latt of the
67
two crystals is due to the distinct behaviors of mid-frequency (200 to 700 cm?1) phonons,
which are slow optic phonons in PV-MgSiO3, but much faster acoustic phonons in MgO. The
flatter dispersions in PV-MgSiO3 also likely increase the phonon scattering rates even when
the lattice anharmonicity coupling strength remains the same, since more phonon-triplets
can simultaneously satisfy the energy and momentum conservation conditions of phonon
scattering. As a result, the PV-MgSiO3 modes have lifetimes nearly one order of magnitude
shorter than those in MgO. It is worthy to note that in PV-MgSiO3, phonon Gr?uneisen ratios,
a parameter that is commonly adopted to quantify lattice anharmonicity, are less than 30%
lower than those in MgO; Yet the comparison of the phonon lifetimes shows that the lifetimes
of most phonon modes of PV-MgSiO3 are at least one order magnitude shorter than those
of MgO. This result indicates that any empirical models that correlate the Gr?uneisen ratios
and phonon lifetimes are only qualitative at best. Realistic evaluation of phonon lifetimes of
individual modes requires an accurate assessment of the number of the conservation-allowed
phonon triplets, the anharmonic coupling strengths, and the frequencies of the phonons.
5.5 Conclusions
In summary, we have calculated the harmonic and anharmonic lattice dynamics of
PV-MgSiO3 based on 20-atom unit cells. The density dependence of inter-atomic force
constants, phonon group velocities and phonon lifetimes has been studied. Our study
reveals a week pressure dependence of mode group velocities, and the phonon lifetimes of
optical modes. Our calculations further demonstrate that the pressure dependence of ?latt is
clearly material-dependent, and simple elastic theories are likely insufficient. Many phonon
group velocities increase with increasing pressure in both crystals, contributing positively
to the pressure dependence of the lattice thermal conductivity. However, acoustic and
optic phonons have different pressure dependencies. While the group velocities of acoustic
phonon scale approximately as the square root of bulk modulus, there is no simple relation
68
between elastic properties and the group velocities of optic phonons. Group velocities of some
mid-frequency optic modes in PV-MgSiO3 are found to decrease with increase of pressure,
leading to a much complicated pressure dependence of the averaged phonon group velocity.
Unlike the high-frequency optic phonons, the mid-frequency optic phonons in PV-MgSiO3
still have noticeable contribution to the bulk lattice thermal conductivity, and they are
largely responsible for the distinct difference in the pressure dependencies of lattice thermal
conductivity in the two crystals.
69
Chapter 6
LATTICE THERMAL CONDUCTIVITY OF FERROPERICLASE
6.1 Introduction
Ferropericlase (fp)-Mg(1?x)FexO is thought to be the second most abundant lower mantle
(LM) mineral with nearly 20% of the volume fraction141,142. Its thermal conductivity is
important to our understanding of the heat flow across the core-mantle boundary (CMB).
Yet, it remains poorly constrained due to the challenges in carrying out measurements at
LM conditions. We report here, a direct calculation of the lattice thermal conductivity of
fp-Mg(1?x)FexO at lower mantle conditions without any empirical extrapolation. Using our
calculated lattice thermal conductivity of Fe-free MgO crystal101 as the starting point, we
have evaluated the iron effects on lattice thermal conductivity of fp by calculating the phonon
scattering rates due to lattice anharmonicity and Mg/Fe mass disorder within the vibrational
virtual crystal approximation (vVCA). We find that phonon scattering due to Fe/Mg mass
disorder is significant. Our study based on an iron content of 12.5% shows a significant
lowering of lattice thermal conductivity of fp even at the high temperature conditions of the
LM. At 3000K and 135GPa we predict the ? of fp to be 11Wm?1K?1 which is nearly 20%
of the ?MgO.
6.2 Static Equation of State
The pressurevolume relation was derived using static equation of state (EOS) calculated
using the LDA+U technique. Using a unit cell 8 times the size of the 2-atom f.c.c.
MgO unit cell, one of the eight Mg atoms was substituted with Fe. The Brillouin zone
70
integration for electronic energy was calculated using the tetrahedron method over a 8?8?8
Monkhorst-Pack grid. The calculated energy v.s. volume data shown in figure 6.1 were fitted
with a 3rd order BirchMurnaghan EOS. Our calculated EOS parameters are summarised in
table 6.1
Our calculated static equation of state agrees well with experiment.
Table 6.1: Static Equation of state of fp-Mg(1?x)FexO
V0 (?A3 per atom) B0(GPa) B?0
XFe=0.00 (This work: LDA) 9.05 172.7 4.2
XFe=0.125 (This work: LDA+U) 8.85 187.5 4.12
XFe=0.17(Experiment)143 8.92 186 4.6
XFe=0.25(Theory GGA+U)144 9.37 170 4.1
s54 s55 s56 s57 s49s48 s49s49
s45s54s46s55s53
s45s54s46s53s48
s45s54s46s50s53
s45s54s46s48s48
s45s53s46s55s53
s45s53s46s53s48
s32
s69
s40
s101
s86
s47
s97
s116
s111
s109
s41
s86s111s108s117s109s101s47s97s116s111s109s32s40s197
s51
s41
s83s116s97s116s105s99s32s69s113s117s97s116s105s111s110s32s111s102s32s83s116s97s116s101
s77s103
s40s49s45s120s41
s70s101
s120
s79
s88s61s49s50s46s53s37
Figure 6.1: Static EOS of ferropericlase fitted to 3rd order BM EOS
71
6.3 Vibrational Virtual Crystal Approximation of Lattice Dynamics of
Ferropericlase
The harmonic and anharmonic lattice dynamics were calculated using our recently
developed finite displacement supercell technique112 based on a 128-atom supercell. Atomic
forces are derived using the first-principles LDA+U method. The equivalent B1-structured
force constants matrices were reconstructed using the vVCA approach described in section
2.6. The cation mass was replaced by the average mass of Fe and Mg in the phonon
calculations. We assume that 3rd order lattice anharmonicity tensors are not sensitive to
the Fe-content in the dilute limit. Figure 6.2 shows the phonon mode frequencies and
group velocities of Fe-free MgO compared to ferropericlase (fp) with 12.5% Fe content
calculated based on 128-atom supercell. The black line/symbols shows the result obtained
by replacing the mass of Mg in the MgO calculation by the weighted mass Mg0.875Fe0.125.
The blue lines/symbols show the result obtained from the vVCA approach. The overall
trend consistently shows softening of phonon mode frequencies and group velocities with
increasing Fe content. Based on the similar harmonic lattice dynamics, we assume that
the lattice anharmonicity are similar, so we use the Fe-free MgO lattice anharmonicity to
evaluate anharmonicity induced scattering rates for fp
72
s48
s50s48s48
s52s48s48
s54s48s48
s56s48s48
s49s48s48s48
s32
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s70s101s114s114s111s112s101s114s105s99s108s97s115s101s32s45s32s77s103
s40s49s45s120s41
s70s101
s120
s79
s40s98s41s32s80s104s111s110s111s110s32s70s114s101s113s117s101s110s99s105s101s115s40s99s109
s45s49
s41
s88 s75 s76
s40s97s41s32s71s114s111s117s112s32s86s101s108s111s99s105s116s121s40s75s109s47s115s41
s88 s75 s76
s40s77s103s79s41s32s88s32s32s32s61s32s32s48s46s48s48s48
s40s109s97s115s115s41s32s88
s70s101
s32s61s32s32s48s46s49s50s53
s40s118s86s67s65s41s88
s70s101
s61s32s32s48s46s49s50s53
s32s40s77s103s79s41s88
s70s101
s32s61s32s48s46s48s48s48
s32s40s109s97s115s115s41s88
s70s101
s32s32s61s32s48s46s49s50s53
s32s40s118s86s67s65s41s88
s70
s32s61s32s48s46s49s50s53
s32
Figure 6.2: Phonon dispersion relations and mode group velocities in Mg0.875Fe0.125O The
black line/symbols shows the result obtained by replacing the mass of Mg in the MgO
calculation by the weighted mass Mg0.875Fe0.125. The blue lines/symbols show the result
obtained from the vVCA approach where the interatomic force constants are directly
calculated as described in section 2.6. Fe substitution is seen to have the effect of softening
phonon frequencies and mode group velocities.
s48s46s48
s48s46s53
s49s46s48
s49s46s53
s50s46s48
s32
s88 s75
s76
s40s77s103s79s41s32s32s88
s70s101
s32s61s32s32s48s46s48s48s48
s40s109s97s115s115s41s32s32s32s88
s70s101s32
s61s32s32s48s46s49s50s53
s40s118s86s67s65s41s32s88
s70s101
s32s61s32s32s48s46s49s50s53
s77
s111
s100
s101
s32
s71
s114
s252
s110
s101
s105
s115
s101
s110
s32
s80
s97
s114
s97
s109
s101
s116
s101
s114
s86s111s108s117s109s101s47s97s116s111s109s32s61s32s54s46s53s48s32s197
s51
s70s101s114s114s111s112s101s114s105s99s108s97s115s101s32s45s32s77s103
s40s49s45s120s41
s70s101
s120
s79
Figure 6.3: Mode Gr?unesein parameter
73
6.4 Phonon Lifetimes
s49s51 s49s52 s49s53 s49s54 s49s55 s49s56
s48
s53s48
s49s48s48
s49s53s48
s50s48s48
s50s53s48
s51s48s48
s51s53s48
s52s48s48
s32
s32
s40
s112
s115
s41
s97s41s32s65s99s111s117s115s116s105s99s32s32s109s111s100s101s115
s40s77s103
s40s49s45s120s41
s70s101
s120
s41s79
s49s51 s49s52 s49s53 s49s54 s49s55 s49s56
s48
s49s48
s50s48
s51s48
s52s48
s53s48
s54s48
s55s48
s56s48
s57s48
s49s48s48
s40s77s103
s40s49s45s120s41
s70s101
s120
s41s79
s98s41s32s79s112s116s105s99s97s108s32s109s111s100s101s115
s32
s32
s49s51 s49s52 s49s53 s49s54 s49s55 s49s56
s48
s53s48
s49s48s48
s49s53s48
s50s48s48
s50s53s48
s51s48s48
s51s53s48
s52s48s48
s77s103s79
s32
s32
s40
s112
s115
s41
s86s111s108s117s109s101s47s85s110s105s116s32s67s101s108s108s32s40s197
s51
s41
s49s51 s49s52 s49s53 s49s54 s49s55 s49s56
s48
s49s48
s50s48
s51s48
s52s48
s53s48
s54s48
s55s48
s56s48
s57s48
s49s48s48
s77s103s79
s32
s32
s86s111s108s117s109s101s47s85s110s105s116s32s67s101s108s108s32s40s197
s51
s41
Figure 6.4: Volume/density dependence of phonon lifetimes of acoustic and optical phonon
modes in fp as predicted by vVCA compared to corresponding modes in MgO. Panels under
a) show the density dependence of the lifetimes for acoustic modes and b) shows the lifetimes
of optical modes. The longest optical lifetime is nearly a factor of 10 lower than the longest
acoustic mode, confirming the notionthat acoustic modesare moreeffective heat transporters
than the optical modes.
The phonon lifetimes are evaluated based on an 16?16?16 k-point grid for a total
of 4096 k-points, 145 of which are irreducible. The phonon-phonon scattering rates are
directly evaluated from for all the irreducible k-points in the Brillouin zone and life times
in the whole Brillouin Zone (BZ) are obtained using the symmetry relations between the
irreducible k-points and the dependent ones. Figure 6.4 shows the density dependence of all
the acoustic and optical phonon modes at the 145 irreducible k-points in the Brillouin zone
74
for MgO and fp respectively. Quantitatively, we our study reveals that the lower frequency
(acoustic) modes account for nearly 95% of the total conductivity in MgO and fp. It is
observed that the phonon lifetimes of the accoustic modes in fp are almost factor of two
lower than those in MgO. In both cases, the modes having longer lifetimes show increasing
liftimes with increasing density (decreasing volume). The short lived modes are almost
insensitive to pressure. Most of the modes have lifetimes that vary almost linearly with
volume/density.
6.5 Lattice Thermal Conductivity of fp-Mg(1?x)FexO
Figure 6.5 shows that the lattice thermal conductivities associated with anharmonicity
scattering are comparable in MgO and fp (see the orange and blue dashed lines). Including
scattering due to mass disorder results in a reduction in thermal conductivity by a factor of
4-10
s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48 s51s48s48s48
s48
s53s48
s49s48s48
s49s53s48
s50s48s48
s50s53s48
s51s48s48
s51s53s48
s65s116s32s51s48s48s75
s97s110s104
s47
s116s111s116s97s108
s32 s49s48
s86s61s49s52s46s48s48s197
s51
s47s85s110s105s116s32s99s101s108s108
s32
s116s111s116s97s108
s40s102s112s41
s32
s97s110s104s32
s40s102s112s41
s32
s97s110s104
s40s77s103s79s41
s40
s87
s109
s45
s49
s75
s45
s49
s41
s84s101s109s112s101s114s97s116s117s114s101s40s75s41
s65s116s32s51s48s48s48s75
s97s110s104
s47
s116s111s116s97s108
s32 s52
Figure 6.5: Thermal conductivity associated with anharmonicity
75
Figure 6.7 presents the thermal conductivity of fp for xFe = 0.125 in the LS. ? was
directly calculated at 5 density points and from 300K to 3000K without any empirical
extrapolation. ?(T,P) is then derived using our first-principles calculated thermal EOS
V(T,P).
s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48 s51s48s48s48
s48s46s48s48
s48s46s48s50
s48s46s48s52
s48s46s48s54
s48s46s48s56
s48s46s49s48
s48s46s49s50
s48s46s49s52
s48s46s49s54
s32s49s50s46s48s48s32
s32s49s51s46s48s48
s32s49s52s46s48s48
s32s49s54s46s48s48
s32s49s56s46s48s48
s77s103s79
s45
s49
s40
s87
s45
s49
s109
s75
s41
s84s101s109s112s101s114s97s116s117s114s101s40s75s41
s85s110s105s116s32s99s101s108s108s32s118s111s108s117s109s101s40s197
s51
s41
Figure 6.6: Temperature dependence of ??1 in MgO. ? was calculated at 5 distinct densities
corresponding to unit cell volumes shown on the figure, at 12 different temperatures ranging
from 300 K to 3000K. Our calculated ? follows a T?1 dependence.
76
s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48 s51s48s48s48 s51s53s48s48 s52s48s48s48
s48s46s48s53
s48s46s49s48
s48s46s49s53
s48s46s50s48
s48s46s50s53
s48s46s51s48
s48s46s51s53
s48s46s52s48
s70s101s114s114s111s112s101s114s105s99s108s97s115s101s32s45s32s77s103
s40s49s45s120s41
s70s101
s120
s79
s32s49s50s46s48s48s32
s32s49s51s46s48s48
s32s49s52s46s48s48
s32s49s54s46s48s48
s32s49s56s46s48s48
s45
s49
s40
s87
s45
s49
s109
s75
s41
s84s101s109s112s101s114s97s116s117s114s101s40s75s41
s85s110s105s116s32s99s101s108s108s32s118s111s108s117s109s101s40s197
s51
s41
Figure 6.7: Temperature dependence of ??1 in fp calculated at the corresponding unit cell
volumes for the Fe-free MgO results presented in figure 6.6. ? does not quite follow the T?1
dependence seen in MgO.
77
s48 s50s53 s53s48 s55s53 s49s48s48 s49s50s53 s49s53s48 s49s55s53 s50s48s48
s48
s49s48s48
s50s48s48
s51s48s48
s52s48s48
s53s48s48
s50s48s48s48s75
s49s32s32s61s32s32s49s50s32s197
s51
s32
s50s32s32s61s32s32s49s51s32s197
s51
s51s32s32s61s32s32s49s52s32s197
s51
s52s32s32s61s32s32s49s54s32s197
s51
s53s32s32s61s32s32s49s56s32s197
s51
s50s53s48s48s75
s51s48s48s48s75
s49s53s48s48s75
s49s48s48s48s75
s51s48s48s75
s40
s87
s109
s45
s49
s75
s45
s49
s41
s80s114s101s115s115s117s114s101s32s40s71s80s97s41
s53s48s48s75
s67s97s108s99s117s108s97s116s101s100s32 s40s80s44s86s41s32s100s101s112s101s110s100s101s110s99s101s32s98s97s115s101s100s32s111s110s32s116s104s101s114s109s97s108s32s69s79s83s32s111s102s32s77s103s79
s77s103s79
s85s110s105s116s32s67s101s108s108s32s86s111s108s117s109s101
Figure 6.8: Pressure dependence of isothermal ? in MgO derived from our LDA calculated
thermal EOS. For each isotherm, we used our calculated thermal EOS to determine the
pressure corresponding to the unit cell volumes 12, 13, 14, 16 & 18 ?A3 labeled 1, 2, 3, 4
& 5 along the isotherms: 300 K - 3000 K. Only seven of the 12 isotherms are shown, for
clarity. By fitting this to a second order polynomial in P (? = a0 +a1P +a2P2) we used the
coefficients a0,a1,a2 from the 12 isotherms to predict ? at isobaric conditions.
s48 s50s53 s53s48 s55s53 s49s48s48 s49s50s53 s49s53s48 s49s55s53 s50s48s48
s48
s49s48
s50s48
s51s48
s52s48
s50s48s48s48s75
s49s32s32s61s32s32s49s50s32s197
s51
s32
s50s32s32s61s32s32s49s51s32s197
s51
s51s32s32s61s32s32s49s52s32s197
s51
s52s32s32s61s32s32s49s54s32s197
s51
s53s32s32s61s32s32s49s56s32s197
s51
s50s53s48s48s75
s51s48s48s48s75
s49s53s48s48s75
s49s48s48s48s75
s51s48s48s75
s40
s87
s109
s45
s49
s75
s45
s49
s41
s80s114s101s115s115s117s114s101s32s40s71s80s97s41
s53s48s48s75
s67s97s108s99s117s108s97s116s101s100s32 s40s80s44s86s41s32s100s101s112s101s110s100s101s110s99s101s32s98s97s115s101s100s32s111s110s32s118s86s67s65s32s116s104s101s114s109s97s108s32s69s79s83s32s111s102s32s102s112
s102s112s45s40s77s103
s40s49s45s120s41
s70s101
s120
s41s79s32s32s32s32s32s32s32s32s32s32s32s32s32s120s32s61s48s46s49s50s53
s85s110s105s116s32s67s101s108s108s32s86s111s108s117s109s101
Figure 6.9: Pressure dependence of ? at isothermal conditions in fp. This pressure
dependence was derrived using the same approach described under figure 6.8
78
s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48 s51s48s48s48
s48
s53s48
s49s48s48
s49s53s48
s50s48s48
s50s53s48
s51s48s48
s51s53s48
s52s48s48
s32s48s71s80s97
s32s53s48s71s80s97
s32s49s48s48s71s80s97
s32s49s53s48s71s80s97
s40
s87
s109
s45
s49
s75
s45
s49
s41
s84s40s75s41
s77s103s79
s48 s49s53 s51s48 s52s53 s54s48 s55s53 s57s48 s49s48s53 s49s50s48 s49s51s53
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s49s56
s50s48
s50s50
s50s52
s50s54
s50s56
s51s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s49s56
s50s48
s50s50
s50s52
s50s54
s50s56
s51s48
s32s51s48s48s75
s32s53s48s48s75
s32s49s48s48s48s75
s32s49s53s48s48s75
s32s50s48s48s48s75
s32s50s53s48s48s75
s32s51s48s48s48s75
s80s114s101s115s115s117s114s101s40s71s80s97s41
Figure 6.10: Main: Temperature dependence of isobaric ? in MgO. Using the coefficients
from ? = a0 + a1P + a2P2 we derived the isobaric thermal conductivity. Four isobars 0,
50, 100 and 150 GPa are shown. Inset: Pressure variation of normalized ?. ?0 represents
the thermal conductivity at 0 GPa for each isotherm. From 300 K to 3000K, ?/?0 increases
from 8 to 25 at 135 GPa.
s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48 s51s48s48s48
s48
s53
s49s48
s49s53
s50s48
s50s53
s51s48
s51s53
s52s48
s52s53
s53s48
s32s48s71s80s97
s32s53s48s71s80s97
s32s49s48s48s71s80s97
s32s49s53s48s71s80s97
s40
s87
s109
s45
s49
s75
s45
s49
s41
s84s40s75s41
s40s77s103
s40s49s45s120s41
s70s101
s120
s41s79s32s32s32s32s32s32s120s32s61s32s48s46s49s50s53
s49s53 s51s48 s52s53 s54s48 s55s53 s57s48 s49s48s53 s49s50s48 s49s51s53
s49s46s48
s49s46s53
s50s46s48
s50s46s53
s51s46s48
s51s46s53
s52s46s48
s52s46s53
s53s46s48
s53s46s53
s54s46s48
s49s46s48
s49s46s53
s50s46s48
s50s46s53
s51s46s48
s51s46s53
s52s46s48
s52s46s53
s53s46s48
s53s46s53
s54s46s48
s32s51s48s48s75
s32s53s48s48s75
s32s49s48s48s48s75
s32s49s53s48s48s75
s32s50s48s48s48s75
s32s50s53s48s48s75
s32s51s48s48s48s75
s80s114s101s115s115s117s114s101s40s71s80s97s41
Figure 6.11: Main: Temperature dependence of isobaric ? in Mg(1?x)FexO; Inset: Pressure
dependence of ??0 from 0GPa to 135 GPa. ?0 is the thermal conductivity at 0GPa for each
isotherm. From 300 K to 3000K, ?/?0 increases from 2.4 to 5 at 135 GPa. This is almost a
factor of 5 lower than in MgO.
79
6.6 Conclusions
Substituting Mg with Fe lowers lattice thermal conductivity through two mechanisms:
(1) adding additional disorder that scatters phonons and therefore further reduces phonon
lifetimes, and (2) lowering phonon group velocities because of weaker chemical bonds and
larger Fe mass. As a first step to quantitatively estimate the factor of reduction due to Fe
substitution, we neglect the changes due to chemical difference of Mg and Fe, and focus only
on the effects of mass difference, i.e. approximating the iron-bearing mineral solid solution
systems with the same force constants and lattice anharmonicity, but with additional Mg/Fe
mass disorder and heavier average mass (Mg1?xFex). Our approximation provides the upper
limits of ?latt in the solid solution systems, and it is valid for the Mg-rich system. We
also carried out a test based on more computationally intensive vibrational virtual crystal
approximation (vVCA) method for Fe-bearing Mg1?xFexO model containing 12.5% low-spin
ferrous Fe2+. We confirmed that the mass effects are significantly larger than those due to
the change of force constants at the high temperatures.
80
Chapter 7
Implications for Heat Flow Across the Core-Mantle Boundary (CMB)
Parts of this Chapter 7 have been submitted for publication in Nature Geoscience by
Tang, Ntam, Dong, Rainey, and Kavner.
7.1 Introduction
The core/mantle boundary (CMB) is the major thermal boundary layer in the Earth?s
deep interior, and the heat flux cross the CMB (JCMB ) governs the thermochemical
trajectory of the Earth. Yet, estimates of JCMB are widely scattered, ranging from 5 to 15
terawatts5. The lower-bound-estimates suggest that additional sources of heat in the mantle
such as enhanced internal heating or an additional thermal boundary layer are required to
accommodate the observed budget of 47 ? 3 terawatts from the surface of the Earth1,5.
The upper-bound-estimates suggest a recent evolution of the inner core/outer core system
and high internal temperatures for early Earth conditions12. To constrain these interlocked
deep-Earth evolution models, we provide a comprehensive evaluation of lower mantle (LM)
thermal conductivity ?LM, anchored on our first-principles calculations of the pressure- and
temperature-dependent lattice contribution to the heat transport in the most abundant
lower mantle mineral (Mg,Fe)SiO3 perovskite, which demonstrate a heretofore unrecognized
shallow pressure dependence. Including a re-evaluation of conflicting data sets constraining
the radiative contribution to overall ?, we show that values for ?LM are much lower than
previously thought, and approximately constant throughout the lower mantle, with little
sensitivity to lateral variations in temperature.
81
The thermal conductivity ?, a material property relating the heat flux J to the
temperature gradient ?T, via Fouriers law, J =??T, is notoriously difficult to measure,
especially at the extreme conditions relevant to Earth and planetary interiors. Discrepancies
exist among measured and calculated ? values for the important deep Earth materials
magnesium silicate perovskite (PV-MgSiO3 Osako et al134, Manthilake et al107 and periclase
(MgO)48,58,101,106, due to large uncertainties arising from differing models used to extrapolate
those values to high pressures and high temperatures19,112,145,146, and confusion regarding the
relative importance of contributions to overall thermal conductivity from lattice phonons
(?latt) versus photon (?rad)147,148. The goal of this multidisciplinary study is to address all of
these issues resulting in a robust determination of lower mantle thermal conductivity (?LM)
and its dependence on pressure (P), temperature (T), and composition. A major challenge
to constrain ?LM is the lack of reliable physical models for ?latt of complex minerals at
extreme P-T conditions. Based on the density functional theory (DFT) and the kinetic
Peierls-Boltzmann transport theory, we implemented a first-principles method to directly
predict ?latt of pure crystals without any empirical P-T extrapolations details reported in
chapter 5. This method has been recently used to calculate ?latt for MgO101 giving good
agreement with independent experiments48 and theory58. However, it is significantly more
challenging to accurately compute all the phonon-phonon scattering rates in the 20-atom
orthorhombic unit-cell of PV-MgSiO3 crystals than in the 2-atom cubic unit-cell of MgO
crystals. In this study, we further improved the numerical efficiency of our algorithm and
used more than 2 million CPU hours at massive parallelized computer clusters to calculate
?latt of PV-MgSiO3 at 36 P-T conditions, ranging from that of ambient to that at the CMB.
Our calculated value for iron-free PV-MgSiO3 at ambient conditions (gray lines in Figure
7.2) is in good agreement with one set of measurements under corresponding conditions134,
but noticeably lower than a recent set107. Two important outcomes of our calculations
are the significant reduction of ?latt of PV-MgSiO3 in comparison with MgO, and the
82
much lower pressure-dependence of the thermal conductivity of PV-MgSiO3 than that of
MgO. These unexpected results arise from the detailed differences in the phonon behavior
between the high-symmetry, 2-atom unit cell MgO and the lower symmetry, 20-atom unit
cell PV-MgSiO3, and are at the root of our estimate of low thermal conductivity and weak
pressure dependence. Our theoretical method provides a series of testable hypotheses about
how crystalline symmetry, structure, and pressure affect thermal conductivity of mineral
insulators
7.2 Results
The calculated phonon spectra of the iron-bearing solid solutions are shown in Figure
7.1. The scattering rates due to Mg/Fe mass disorder are evaluated using the Fermis
golden rule equation 3.16 based on Fe-bearing phonon spectra. In the dilute limit, the
total phonon scattering rates in the Fe-bearing solid solution systems are approximated as
sum of scattering rates due to both lattice anharmonicity and the Mg/Fe mass disorder.
Considering the predominant effect of Mg/Fe mass disorder in the current calculation, we
assume that lattice anharmonicity induced phonon scattering rates in the mineral solid
solutions are the same as those in the iron-free crystal. Iron (Fe) content affects the ?latt
of lower mantle minerals. However, it is difficult to untangle its contribution from that of
lattice anharmonicity in experimental measurements. In this study, we analyzed the Mg/Fe
mass disorder induced phonon scattering of substituting 12.5% of the Mg for Fe in both the
perovskite and oxide. Including iron further lowers ?latt , although the reduction at high
pressures and temperatures is not as significant as it is closer to ambient conditions (Fig1).
While the expected T?1 dependence of ?latt for iron-free PV-MgSiO3 at constant density is
observed, the calculation for PV-MgSiO3 down a typical lower mantle geotherm shows ?latt
is slightly below ?1 Wm?1K?1 and approximately constant as a function of depth (Fig2),
in contrast to the much higher-valued ?latt and steeper depth dependence for (Mg,Fe)O101.
83
Figure 7.1: Vibrational phonon frequencies of (Mg1?xFex)O (a) and
(Mg1?xFex)SiO3-perovskite (b) with iron concentration x= 12.5% at deep mantle pressure
84
Earth?s mantle is optically thick, and localized radiative heat transport may contribute to
the overall diffusive heat transport. For this case, the photon radiative contribution to ?
can be estimated using the Rosseland mean approximation149,150. We provide two estimates
based on a consistent analysis of the two sets of measured absorption spectra for Fe-bearing
perovskite at high pressures147,148. Choosing between these two estimates is beyond the scope
of this work, but we were comforted to note that the two spectra showed a similar spectral
structure and pressure dependence. The resulting estimates of radiative contribution to heat
flow differ by a factor of 2 to 3; however we note that this is a significant improvement of
the order-of-magnitude difference claimed in the original studies. This partial reconciliation
of two disparate data sets arises from treating the optical absorption data consistently in
the Rosseland calculation of ?rad. Generally the calculated ?rad is affected by scattering
(e.g. grain boundary scattering) as well as absorption. Here we neglect a possible scattering
component due to the short absorption length scales (? 100 m) compared with the expected
grain sizes of the lower mantle (? 1 mm). The presence of additional scattering will only
serve to lower the overall radiative thermal conductivity; thus we consider our estimates
for ?rad to be upper bounds. The two dashed red lines in Figures 7.3 and correspond to
calculations incorporating the two estimates for ?rad at constant pressure as a function of
temperature (Fig.2), and down a lower mantle geotherm (Figure 7.4). Our two calculations
for the ?rad of perovskite are similar in magnitude to the computed ?latt, but with a strong
positive temperature dependence (Figure 7.3) such that at low temperatures the lattice
values dominate and at high temperatures the radiative values dominate. Total thermal
conductivity is the sum of the lattice phonon and radiative contributions. Our resulting two
estimates of total ? are approximately constant with temperature (Figure 7.3) and depth in
the lower mantle (Figure 7.4). This temperature- and pressure- stability of the perovskite
thermal conductivity is a robust result of our analysis. Modest temperature perturbations
about the geotherm will increase the radiative contribution to ? while attenuating the lattice
85
contribution. Changes in iron content or other impurities will likely have the same effect
on both radiative and lattice contributions; for example an increase in iron content will
decrease the value of both components. Finally, we determine the total ?LM based on the
Maxwell-Garnett model, which is an effective medium approximation for examining the total
conductivity of a composite where a minor phase exists as random inclusions within a matrix
phase151,152. As a first-order representation of the lower mantle, we consider a composite
material consisting of a matrix of iron-bearing perovskite with ? 20% iron-bearing oxide
existing as isolated inclusions. Because the length scales for both phonon scattering (? nm)
and typical optical depth for photon scattering (? 10-100 m) are both likely to be less than
the average grain size of individual perovskite and oxide grains (? 1 mm), we perform the
calculation by first summing the lattice and radiative contributions to ? for each phase, and
then combining the two phases in the Maxwell Garnett composite model Figure 3 shows
Figure 7.2: Calculated lattice thermal conductivity for MgSiO3 perovskite as a function of
pressure and temperature. Calculations are shown for Fe-free crystal (translucent symbols
and lines) and mineral solid solution including 12.5% Fe (bold color symbols and lines). Also
plotted are experimental data from Osaka & Ito, 1991 on MgSiO3 and from Manthilake et
al., 2011 on (Mg,Fe)SiO3 at temperatures from 473 to 1073 K.
86
Figure 7.3: Isobaric thermal conductivity as a function of temperature. Computed lattice
conductivity (solid purple line) and calculated contribution from radiative heat flow (red
dashed and dotted lines, from optical absorption values provided by Goncharov et al., 2008
and Keppler et al., 2008 respectively, as a function of temperature at constant mid-mantle
pressure of 80 GPa. The total thermal conductivity of perovskite is given by the sum of
?rad +?latt (purple dotted and dashed lines, corresponding to the two estimates for radiative
contributions to the thermal conductivity)
87
Figure 7.4: Thermal conductivity calculated down a sample lower mantle geotherm (Jeanloz
) and thermal boundary layer. Lattice thermal conductivity calculated from first principles
estimates for MgO (solid blue line, Tang & Dong, 2010), and MgSiO3 perovskite (solid
purple line, this study). Thin red lines show calculated radiative contribution to thermal
conductivity from two optical absorption estimates (dashed line: Goncharov et al., 2008;
88
Figure 7.5: A plot summarizing the observational tradeoffs used to determine the total heat
flow across the core-mantle boundary (contoured, labels in terawatts). Commonly assumed
values from Lay et al., 2008 for the thermal boundary thickness, the temperature at the
bottom of the mantle, and the temperature at the top of the core are shown on the right
axis. Vertical solid green lines delineate our newly-determined range of estimates for lower
mantle thermal conductivity.
89
two estimates for total ?LM down a calculated mantle temperature profile153, including the
effect of a ? 200 km-thick thermal boundary layer between the convecting lower mantle
and the outer core. Our calculations suggest that thermal conductivity will increase across
the steep temperature gradient of the thermal boundary layer, mostly due to increased
radiative transfer. The thermal conductivity is lower if lower temperatures are assumed
at the base of the convecting mantle, but the thermal boundary temperature dependence
likely becomes steeper due to an increased gradient across the core/mantle boundary. When
the uncertainties in the temperatures at the bottom of the convecting lower mantle and
the top of the outer core5 are mapped onto an uncertainty in the thermal conductivity, we
obtain values of values of 3.5 to 5.5 at the bottom of the convecting lower mantle. The
thermal conductivity values rise across the thermal boundary layer if there is an increase
in the radiative contribution to total heat flow in the lowermost boundary layer, which
would occur in the absence of major phase and/or chemistry changes within the thermal
boundary. These values represent a significant adjustment of the ?LM downward from the
typically used value of 10 Wm?1K?1. Models of thermal evolution of the mantle and the
core are coupled via the ?LM. Figure 4 illustrates how tradeoffs in the assumed values for
thickness of the lower mantle thermal boundary, temperature at the top of the outer core,
and temperature at the foot of the convecting mantle influence estimates of the total heat
crossing the core mantle boundary. The contours in Figure 4 show the calculated total heat
flow across the CMB given estimates of the thermal boundary layer conductivity (X-axis) and
the temperature gradient across the boundary layer (left Y-axis). The temperature gradient
across the boundary layer depends on estimates of the temperature in the outer core, the
temperature of the bottom of the convecting mantle, and the thickness of the boundary layer;
all of these parameters have some uncertainty. Often-used values (temperature of 2800 K
at the bottom of the convecting mantle, 3800 K at the top of the core5, and a thermal
boundary layer thickness that varies from 100 km to 200 km154 are shown as a dashed
90
horizontal line. Using these values, the total heat flow across the core mantle boundary
inferred from our new thermal conductivity estimates is calculated to be in the range of
3.5 to 5 TW, much lower than previous estimates of ? 10 TW. Our new estimates of heat
flux across the core/mantle boundary relax some of the constraints on timing of inner core
growth12, but pose difficulties (but not insurmountable ones) in accounting for the global
heat balance of the mantle5,155. Our new thermal conductivity values permit a heat flux
larger than 5 TW across the core/mantle boundary, but require a steeper temperature
gradient to do this. For example, a 10 TW heat flux across the core/mantle boundary
can be most easily accomplished by assuming a much thinner thermal boundary layer, equal
to 75 to 100 km. This is permissible within observational constraints, and is consistent with
seismic observations of the boundary layer width above the core/mantle boundary156). If
the seismological constraints are taken as indicative of the thermal boundary layer thickness,
then one must consider a laterally varying thickness of the thermal boundary layer. Our
results suggest that variations of the thermal conductivity as a result of lateral variations in
temperature or pressure will not be significant. However, local changes in modal abundance
of minerals, chemistry such as iron content, presence of melt, and/or mechanical texturing
may all serve to significantly perturb the value of the thermal conductivity. Therefore, our
results suggest that any nonlinear influence that thermal conductivity has in geodynamical
models11 arises not from pressure- or temperature- dependence, but is solely a function of
compositional heterogeneity at the base of the lowermost mantle.
7.3 Discussion
The thermal evolution of the whole Earth is poorly constrained, as laid out in a review
by Lay et al5. As pointed out in that review, one of the keys to resolving that uncertainty is a
better constraint ofthematerialproperty ofthermal conductivity relevant totheEarths lower
mantle. For ? 20 years, a value of ? 10 W/m/K has been used, and this value was based on
91
a 1991 measurement at ambient conditions, with a long extrapolations in both pressure and
temperature. The uncertainties in that extrapolated value are large (over a factor of two)
and therefore our understanding of Earths thermal-mechanical-chemical evolution of both
the mantle and the core systems are uncertain. The previously used pressure extrapolations
of perovskite were based on models of the pressure-dependent behavior of phonons in simple
high-symmetry iron-free crystals. These new show that previous pressure-dependent models
are insufficient in complex materials systems, where contribution from optic phonons are
non-negligible. Combining the evaluation of Mg/Fe mass disorder effects and our prediction
of lattice thermal conductivity of (Mg,Fe)SiO3 results in a lower estimate for thermal
conductivity throughout the Earths deep mantle. Besides the first-principles calculations
themselves, we provide a synthesis of the bottom line contribution for two issues that
contribute to understanding heat flow in the deep Earth: a radiative contribution and the
composite effect, drawing on an extensive multidisciplinary literature that quantitatively
describes these effects. Particularly, the radiative contribution to the overall heat transport
has caused a great deal of uncertainty in the literature, with two papers published in Nature
and Science in 2008 making opposite claims based on similar measurements. As part of this
study, we return to the optical absorption data sets as originally published, and treat them
both in a self-consistent and physically correct manner. The result is a partial reconciliation
of the two estimates, with predictions for radiative conductivity in the lower mantle that
differ by a factor of 2-3, rather than the order-of-magnitude disagreement originally cited in
the literature.
7.4 Conclusions
The end result is a new estimate of the lower mantle thermal conductivity throughout
the lower mantle. The number is still uncertain, mostly due to lingering uncertainties in
the radiative component of thermal conductivity and the composite effect of combining
92
oxide which has much higher thermal conductivity with perovskite, which has much lower.
Even with this uncertainty, two results of our study are extremely robust. First, our study
requires that much lower values of thermal conductivity, ? 3.55.5 W/m/K rather than
previous estimate of 5-15 W/m/K be adopted in our Earths thermal evolution models. This
new range of values has the advantage of relaxing some of the constraints on timing of
the inner core growth, and allowing for a more reasonable initial temperature profile at
the origin of whole Earth evolution. The second robust result is our finding that thermal
conductivity is likely single-valued as a function of depth throughout the lower mantle and
insensitive to local variations in temperature. This arises because the lattice component and
radiative component of thermal conductivities have approximately similar magnitudes, but
with opposite temperature and weak pressure dependencies. This means that geodynamical
models can rest soundly on the assumption of a constant value of thermal conductivity
throughout the convecting lower mantle. This simplifies calculations significantly, and lessens
the likelihood of nonlinearities due to temperature-dependent thermal conductivity.
93
Chapter 8
CONCLUSIONS AND FUTURE WORK
In summary, I have presented here the results of our recent first-principles calculations
of lattice thermal conductivities of Earth forming minerals using our recently developed
computational techniques based on first-principles DFT. Under methodology development,
we have optimized our recently developed reals space supercell finite displacement algorithm
to extract the harmonic and anharmonic inter-atomic potentials for complex material
systems up to 20-atoms per unit cell and also implemented the vibrational virtual crystal
approximation for the study of mineral solid-solutions.
Using the extracted third order lattice anharmonicity I have evaluated the phonon
scattering rates of individual phonon modes based on quantum mechanical scattering theory
at high temperature and pressure conditions relevant to the Earth?s lower mantle, without
any empirical fitting. Our calculation of phonon-phonon scattering rates was within the
framework of the single mode excitation approximation (SMEA). To better understand the
pressure dependence of thermal conductivity, I have studied the volume/density dependence
of the group velocities and phonon lifetimes of individual phonon modes. My study reveals
that phonon scattering rates for optical modes are almost insensitive to pressure. The
pressure dependence of the lattice thermal conductivity appears to be largely determined by
the pressure behavior of the scattering rates of acoustic phonon modes. Acoustic phonon
modes are found to be the effective carriers of heat in MgO, (Mg,Fe)O (3 acoustic + 3 optic)
and ??Al2O3 (3-acoustic +27 optic) , while the optical modes account for most of the heat
transport in PV-MgSiO3 (3 acoustic + 57 optic). Table 8.1 summarizes the contribution to
the overall lattice thermal conductivity at ambient conditions:
94
Table 8.1: Contributions to overall ?latt from acoustic and optical phonon modes in MgO,
??Al2O3 and PV-MgSiO3
System ?acc (Wm?1K?1) ?opt(Wm?1K?1) ?latt (Wm?1K?1)
MgO 46.30 8.45 54.75
??Al2O3 16.59 9.70 26.29
PV-MgSiO3 1.63 3.08 4.71
We have successfully predicted the lattice thermal conductivity of PV-MgSiO3 at lower
mantle conditions. While discrepancies remain between our predicted values at ambient
conditions and measurements from experiment, we do observe unusual behavior in existing
experimental data. We are currently repeating our calculation to systematically test for:
different Brillouin zone K-point sampling grids, third order lattice anharmonicity cut-offs
and the GGA formalism.
Using our recent implementation of the vibrational virtual crystal approximation we
have studied the effect of Fe on the thermodynamic and thermal transport properties of
ferropericlase (Mg,Fe)O based on and Fe content of 12.5%. Fe is observed to lower the
thermal conductivity of ferropericlase. We have also estimated the effect of Fe on the thermal
conductivity of (Mg,Fe)SiO3 in the dilute limit, assuming same harmonic and anharmonic
lattice dynamics as Fe-free PV-MgSiO3, and discussed the implications for heat flow in the
Earth?s lower mantle.
Lifetimes of phonon modes in MgO are observed to follow the ? ? ??2 behavior
consistent with Klemens? formular157?159(equation 8.1) while MgSiO3 phonon lifetimes do
not follow this rule.
1
??q,i = ?
2
?q ,i
2KBT
Mv2?q,i
?2?q,i
?maxi (8.1)
Here ? is the phonon lifetime, ?, Vg and ? are the mode Gr?uneisen parameter, mode group
velocity and phonon frequency of a phonon mode (?q,i); ?maxi is the largest frequency in the
ith branch. According to Klemens? rule, low-frequency phonons will have longer lifetimes
and thus contribute more to the thermal conductivity.
95
Future Work
The behavior of phonon mode lifetimes in PV-MgSiO3 deserves further investigation. All
DFT calculations reported here were performed within the LDA framework; It is necessary
to repeat these calculations within the generalized gradient approximation (GGA), and test
the convergence of calculated lattice thermal conductivity on the q-point sampling grid for
the Brillouin Zone. The thermal conductivity calculations reported in this dissertation were
performed within the Peierls-Boltzmann transport theory (PBTT) also known as the phonon
gas model. Using a 2-D lattice model, Sun and Allen160 claim that the phonon gas model
breaks down at high temperatures. How this applies to real complex mineral systems is
not yet well known. In future, it would be desirable to investigate the high temperature
behavior of PV-MgSiO3 using the Green-Kubo formalism. Immediate extensions of the
the work presented in this dissertation would be a self consistent solution of the Boltzmann
transport equation, and a careful investigation the effect of different Fe concentrations on the
thermodynamic and thermal transport properties of (Mg,Fe)O and (Mg,Fe)SiO3 by directly
calculating the harmonic and anharmonic lattice dynamics, taking into account the effect of
spin polarization in Fe.
96
Bibliography
1. D. Anderson, New theory of the Earth (Cambridge Univ Pr, 2007).
2. K. Hirose and T. Lay, Elements 4, 183 (June 2008).
3. T. Tsuchiya, J. Tsuchiya, K. Umemoto, and R. Wentzcovitch, Earth and Planetary
Science Letters 224, 241 (2004).
4. S. Labrosse and C. Jaupart, Earth and Planetary Science Letters 260, 465 (2007).
5. T. Lay, J. Hernlund, and B. Buffett, Nature Geoscience 1, 25 (2008).
6. D. Gubbins, A. Willis, and B. Sreenivasan, Physics of the Earth and Planetary Interiors
162, 256 (2007).
7. J. Brown and R. McQueen, Journal of Geophysical Research 91, 7485 (1986).
8. C. Matyska and D. Yuen, Physics of the Earth and Planetary Interiors 154, 196 (2006).
9. J. B. Naliboff and L. H. Kellogg, Physics of the Earth and Planetary Interiors 161, 86
(2007).
10. A. Van den Berg, E. Rainey, and D. Yuen, Physics of the Earth and Planetary Interiors
149, 259 (2005).
11. A. M. Hofmeister and D. A. Yuen, Journal of Geodynamics 44, 186 (2007).
12. B. Buffett, Geophysical research letters 29, 7 (2002).
13. Y. Syono, Pure and Applied Geophysics (PAGEOPH) 143, 759 (1994).
97
14. H. Hsu, K. Umemoto, Z. Wu, and R. Wentzcovitch, Reviews in Mineralogy and
Geochemistry 71, 169 (2010).
15. H. Spiering, E. Meissner, H. K?oppen, E. M?uller, and P. G?utlich, Chemical Physics 68,
65 (1982).
16. J. Lin, G. Vank?o, S. Jacobsen, V. Iota, V. Struzhkin, V. Prakapenka, A. Kuznetsov,
and C. Yoo, Science 317, 1740 (2007).
17. A. M. Hofmeister, Physics of the Earth and Planetary Interiors 170, 201 (2008).
18. A. Hofmeister, Physics and chemistry of minerals 33, 45 (2006).
19. A. Hofmeister, Proceedings of the National Academy of Sciences 104, 9192 (2007).
20. M. Pertermann and A. Hofmeister, American Mineralogist 91, 1747 (2006).
21. A. M. Hofmeister, Physics of the Earth and Planetary Interiors 180, 138 (2010).
22. V. Solomatov and C. Reese, J. Geophys. Res 113, B07408 (2008).
23. L. lorenz, Ann. Phys. 13, 422 (1881).
24. A. van den Berg, D. Yuen, G. Beebe, and M. Christiansen, Physics of the Earth and
Planetary Interiors 178, 136 (2010).
25. Y. Touloukian, R. Powell, C. Ho, and P. Klemens (1970).
26. M. Graf, S. Yip, J. Sauls, and D. Rainer, Physical Review B 53, 15147 (1996).
27. F. Blatt, Physics of electronic conduction in solids, vol. 191 (McGraw-Hill New York,
1968).
28. K. Ohta, K. Hirose, M. Ichiki, K. Shimizu, N. Sata, and Y. Ohishi, Earth and Planetary
Science Letters 289, 497 (2010).
98
29. A. Oganov, J. Brodholt, and G. Price, Earth and Planetary Science Letters 184, 555
(2001).
30. K. Watari and S. Shinde, MRS Bulletin 26, 440 (2001).
31. S. Shind?e and J. Goela, High thermal conductivity materials (Springer Verlag, 2006).
32. J. Klett, R. Hardy, E. Romine, C. Walls, and T. Burchell, Carbon 38, 953 (2000).
33. W. Kim, J. Zide, A. Gossard, D. Klenov, S. Stemmer, A. Shakouri, and A. Majumdar,
Physical review letters 96, 45901 (2006).
34. D. Clarke, Surface and Coatings Technology 163, 67 (2003).
35. K. Kurosaki, A. Kosuga, H. Muta, M. Uno, and S. Yamanaka, Applied Physics Letters
87, 061919 (2005).
36. M. Dresselhaus, G. Chen, M. Tang, R. Yang, H. Lee, D. Wang, Z. Ren, J. Fleurial, and
P. Gogna, Advanced Materials 19, 1043 (2007).
37. G. Nolas, G. Slack, D. Morelli, T. Tritt, and A. Ehrlich, Journal of Applied Physics
79, 4002 (1996).
38. D. Morelli and G. Meisner, Journal of applied physics 77, 3777 (1995).
39. J. Wu, X. Wei, N. Padture, P. Klemens, M. Gell, E. Garc??a, P. Miranzo, and M. Osendi,
Journal of the American Ceramic Society 85, 3031 (2002).
40. X. Cao, R. Vassen, and D. Stoever, Journal of the European Ceramic Society 24, 1
(2004).
41. R. Vassen, X. Cao, F. Tietz, D. Basu, and D. St?over, Journal of the American Ceramic
Society 83, 2023 (2000).
42. M. Winter and D. Clarke, Journal of the American Ceramic Society 90, 533 (2007).
99
43. A. F. Goncharov, P. Beck, V. V. Struzhkin, B. D. Haugen, and S. D. Jacobsen, Physics
of the Earth and Planetary Interiors 174, 24 (2009).
44. P. Andersson and G. Backstrom, Review of Scientific Instruments 47, 205 (1976).
45. R. Barker Jr, R. Chen, and R. Frost, Journal of Polymer Science: Polymer Physics
Edition 15, 1199 (1977).
46. R. Frost, R. Chen, and R. Barker, Journal of Applied Physics 46, 4506 (1975).
47. H. Ftjjisawa, N. Ftjjii, H. Mizutani, H. Kanamoei, and S. Akimoto (1968).
48. T. Katsura, Physics of the earth and planetary interiors 101, 73 (1997).
49. B. Hakansson, P. Andersson, and G. Backstrom, Review of scientific instruments 59,
2269 (1988).
50. E. Abramson, J. Brown, and L. Slutsky, The Journal of Chemical Physics 115, 10461
(2001).
51. E. Abramson, L. Slutsky, M. Harrell, and J. Brown, The Journal of chemical physics
110, 10493 (1999).
52. M. Chai, J. Brown, and L. Slutsky, Physics and chemistry of minerals 23, 470 (1996).
53. P. Beck, A. Goncharov, V. Struzhkin, B. Militzer, H. Mao, and R. Hemley, Applied
Physics Letters 91, 181914 (2007).
54. A. Goncharov, P. Beck, V. Struzhkin, B. Haugen, and S. Jacobsen, Physics of the Earth
and Planetary Interiors 174, 24 (2009).
55. A. Goncharov, V. Struzhkin, J. Montoya, S. Kharlamova, R. Kundargi, J. Siebert,
J. Badro, D. Antonangeli, F. Ryerson, and W. Mao, Physics of the Earth and Planetary
Interiors 180, 148 (2010).
100
56. A. Hofmeister, Applied Physics Letters 95, 096101 (2009).
57. D. Cahill and J. Li, Copyright 2011 Wen-Pin Hsieh (2011).
58. S. Stackhouse, L. Stixrude, and B. Karki, Physical review letters 104, 208501 (2010).
59. P. Hohenberg and W. Kohn, Phys. Rev 136, B864 (1964).
60. W. Kohn, L. Sham, et al., Phys. Rev 140, A1133 (1965).
61. URL http://cms.mpi.univie.ac.at/vasp/vasp/vasp.html.
62. J. Hafner, Journal of computational chemistry 29, 2044 (2008).
63. X. Gonze, J. Beuken, R. Caracas, F. Detraux, M. Fuchs, G. Rignanese, L. Sindic,
M. Verstraete, G. Zerah, F. Jollet, et al., Computational Materials Science 25, 478
(2002).
64. X. Gonze, Zeitschrift f?ur Kristallographie 220, 558 (2005).
65. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli,
G. Chiarotti, M. Cococcioni, I. Dabo, et al., Journal of Physics: Condensed Matter 21,
395502 (2009).
66. R. Cohen, L. Boyer, M. Mehl, W. Pickett, and H. Krakauer, Perovskite: A Structure
of Great Interest to Geophysics and Material Science (1989).
67. R. Cohen and J. Weitz (1996).
68. R. Cohen, NATO ASI SERIES C MATHEMATICAL AND PHYSICAL SCIENCES
543, 201 (1999).
69. Z. Wu, X. Chen, V. Struzhkin, and R. Cohen, Physical Review B 71, 214103 (2005).
101
70. R. M. Wentzcovitch, Y. G. Yu, and Z. Wu, Reviews in Mineralogy and Geochemistry
71, 59 (2010).
71. K. Thomson, R. Wentzcovitch, and M. Bukowinski, Science 274, 1880 (1996).
72. A. Oganov and S. Ono, Nature 430, 445 (2004).
73. A. Oganov and S. Ono, PNAS 102, 10828 (2005).
74. R. Caracas and R. Cohen, AGU Fall Meeting Abstracts 1, 0191 (2004).
75. R. Caracas and R. Cohen, Geophysical research letters 32, L06303 (2005).
76. H. Cynn, D. Isaak, E. Cohen Ronald, M. Nicol, and O. Anderson, Am. Mineral. 75,
439 (1990).
77. L. Stixrude, R. Cohen, R. Yu, and H. Krakauer, American Mineralogist 81, 1293 (1996).
78. J. Dong and O. Sankey, J. Phys: Condens. Matter 11, 6129 (1999).
79. J. Dong, O. Sankey, S. Deb, G. Wolf, and P. McMillan, Phys. Rev. B 61, 11979 (2000).
80. W. Duan, B. Karki, and R. Wentzcovitch, Am. Mineral. 84, 1961 (1999).
81. R. Caracas and R. Cohen, AGU Fall Meeting Abstracts 1, 07 (2005).
82. R. Caracas and R. Cohen, Geophysical research letters 32, L16310 (2005).
83. G. Steinle-Neumann, L. Stixrude, and R. Cohen, Physical Review B 69, 219903 (2004).
84. D. Machon, P. McMillan, B. Xu, and J. Dong, Physical Review B 73, 094125 (2006).
85. J. Dong, O. Sankey, and G. Kern, Phys. Rev. B 60, 950 (1999).
86. B. Karki and R. Wentzcovitch, Physical Review B 68, 224304 (2003).
87. A. Oganov, M. Gillan, and G. Price, Journal of Chemical Physics 118, 10174 (2003).
102
88. S. Ono, A. Oganov, T. Koyama, and H. Shimizu, Earth Planet. Sci. Lett. 246, 326
(2006).
89. L. Stixrude and R. Cohen, Nature 364, 613 (1993).
90. R. Caracas and R. Cohen, Arxiv preprint cond-mat/0603222 (2006).
91. R. Cohen, Geophysical Research Letters 14, 1053 (1987).
92. J. Tsuchiya, T. Tsuchiya, and R. Wentzcovitch, J. Geophys. Res 110, B02204 (2005).
93. Z. Wu, R. Wentzcovitch, K. Umemoto, B. Li, K. Hirose, and J. Zheng, J. Geophys. Res
113, B06204 (2008).
94. I. Inbar and R. Cohen (1995).
95. D. Isaak, R. Cohen, and M. Mehl, Journal of Geophysical Research 95, 7055 (1990).
96. F. Marton, J. Ita, and R. Cohen, Arxiv preprint physics/0010083 (2000).
97. X. Sha and R. Cohen, Arxiv preprint arXiv:0708.0183 (2007).
98. P. Dorogokupets, Physics and Chemistry of Minerals 37, 677 (2010), ISSN 0342-1791,
10.1007/s00269-010-0367-2.
99. R. Wentzcovitch, L. Stixrude, and J. Rosso, Theoretical and computational methods in
mineral physics: geophysical applications (Mineralogical Society of America, 2010).
100. J. Hafner, C. Wolverton, and G. Ceder, MRS bulletin 31, 659 (2006).
101. X. Tang and J. Dong, Proceedings of the National Academy of Sciences 107, 4539
(2010).
102. D. Evans and G. Morriss, Statistical mechanics of nonequilibrium liquids (Cambridge
Univ Pr, 2008).
103
103. B. Daly, H. Maris, K. Imamura, and S. Tamura, PHYSICAL REVIEW-SERIES B- 66,
24301 (2002).
104. F. M?uller-Plathe, Journal of Chemical Physics 106, 6082 (1997).
105. C. Nieto-Draghi and J. Avalos, Molecular Physics 101, 2303 (2003).
106. N. De Koker, Physical review letters 103, 125902 (2009).
107. G. Manthilake, N. de Koker, D. Frost, and C. McCammon, Proceedings of the National
Academy of Sciences 108, 17901 (2011).
108. R. Cohen, Arxiv preprint cond-mat/9708134 (1997).
109. S. Stackhouse and L. Stixrude, Reviews in Mineralogy and Geochemistry 71, 253
(2010).
110. R. Kubo, Reports on Progress in Physics 29, 255 (1966).
111. M. Green, The Journal of Chemical Physics 22, 398 (1954).
112. X. Tang and J. Dong, Physics of the Earth and Planetary Interiors 174, 33 (2009).
113. W. Frank, C. Els?asser, and M. F?ahnle, Phys. Rev. Lett. 74, 1791 (1995).
114. G. Kresse, J. Furthm?uller, and J. Hafner, Europhys. Lett 32, 729 (1995).
115. A. Garc??a and D. Vanderbilt, Phys. Rev. B 54, 3817 (1996).
116. J. Dong, J. Tomfohr, O. Sankey, K. Leinenweber, M. Somayazulu, and P. McMillan,
Phys. Rev. B 62, 14685 (2000).
117. X. Tang and J. Dong, Bulletin of the American Physical Society (2006).
118. K. Kunc and R. Martin, Phys. Rev. Lett. 48, 406 (1982).
104
119. B. Xu, Ph.D. thesis, Auburn University, AL (2009).
120. T. Barron and M. Klein, Dynamical Properties of Solids (North-Holland, Amsterdam,
1974).
121. J. Fabian and P. Allen, Physical review letters 79, 1885 (1997).
122. N. Funamori and R. Jeanloz, Science 278, 1109 (1997).
123. T. Mashimo, K. Tsumoto, K. Nakamura, Y. Noguchi, K. Fukuoka, and Y. Syono,
Geophys. Res. Lett 27, 2021 (2000).
124. J. Lin, O. Degtyareva, C. Prewitt, P. Dera, N. Sata, E. Gregoryanz, H. Mao, and
R. Hemley, Nat. Mater. 3, 389 (2004).
125. F. Marton and R. Cohen, The American mineralogist 79, 789 (1994).
126. R. Caracas and R. Cohen, Geophys. Res. Lett 32, L06303 (2005).
127. W. Kingery, Journal of the American Ceramic Society 37, 88 (1954).
128. W. Kingery, J. Francl, R. Coble, and T. Vasilos, Journal of the American Ceramic
Society 37, 107 (1954).
129. G. Slack, Physical Review 126, 427 (1962).
130. G. Kresse and J. Furthm?uller, Comput. Mater. Sci 6, 15 (1996).
131. G. Kresse and J. Furhmuller, Phys. Rev. B 47, R558 (1993).
132. R. Caracas and R. E. Cohen, Geophys. Res. Lett 33, 12 (2006).
133. D. Andrault, M. Muoz, N. Bolfan-Casanova, N. Guignot, J.-P. Perrillat, G. Aquilanti,
and S. Pascarelli, Earth and Planetary Science Letters 293, 90 (2010).
134. M. Osako and E. Ito, Geophysical Research Letters 18, 239 (1991).
105
135. P. Bl?ochl, Phys. Rev. B 50, 17953 (1994).
136. G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
137. M. Murakami, K. Hirose, K. Kawamura, N. Sata, and Y. Ohishi, Science 304, 855
(2004).
138. K. Catalli, S. Shim, and V. Prakapenka, Nature 462, 782 (2009).
139. B. Grocholski, K. Catalli, S. Shim, and V. Prakapenka, Proceedings of the National
Academy of Sciences 109, 2275 (2012).
140. A.-L. Auzende, J. Badro, F. J. Ryerson, P. K. Weber, S. J. Fallon, A. Addad, J. Siebert,
and G. Fiquet, Earth and Planetary Science Letters 269, 164 (2008).
141. A. Dziewonski and D. Anderson, Physics of the earth and planetary interiors 25, 297
(1981).
142. R. van der Hilst and H. K?arason, Science 283, 1885 (1999).
143. S. Speziale, V. Lee, S. Clark, J. Lin, M. Pasternak, and R. Jeanloz, J. geophys. Res
112, B10212 (2007).
144. K. Persson, A. Bengtson, G. Ceder, and D. Morgan, Geophysical research letters 33,
16306 (2006).
145. A. Hofmeister, Science 283, 1699 (1999).
146. N. de Koker, Earth and Planetary Science Letters 292, 392 (2010).
147. A. Goncharov, B. Haugen, V. Struzhkin, P. Beck, and S. Jacobsen, Nature 456, 231
(2008).
148. H. Keppler, L. Dubrovinsky, O. Narygina, and I. Kantor, Science 322, 1529 (2008).
106
149. M. Alloy and D. Menezes, Brazilian Journal of Physics 37, 1183 (2007).
150. R. Siegel and J. Howell, Thermal radiation heat transfer, vol. 1 (Taylor & Francis,
2002).
151. R. Progelhof, J. Throne, and R. Ruetsch, Polymer Engineering & Science 16, 615
(1976).
152. O. Levy and D. Stroud, Physical Review B 56, 8035 (1997).
153. R. Jeanloz and S. Morris, Annual Review of Earth and Planetary Sciences 14, 377
(1986).
154. D. Turcotte and G. Schubert, Geodynamics (Cambridge Univ Pr, 2002).
155. K. Collaboration, A. Gando, Y. Gando, K. Ichimura, H. Ikeda, K. Inoue, Y. Kibe,
Y. Kishimoto, M. Koga, Y. Minekawa, et al., Nature Geoscience 4, 647 (2011).
156. J. Hernlund, C. Thomas, and P. Tackley, Nature 434, 882 (2005).
157. R. Tye, Thermal conductivity, vol. 1 (Academic press, 1969).
158. P. Klemens, Proceedings of the Royal Society of London. Series A. Mathematical and
Physical Sciences 208, 108 (1951).
159. P. Klemens, Thermochimica acta 218, 247 (1993).
160. T. Sun and P. Allen, Physical Review B 82, 224305 (2010).
161. G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
162. X. Tang, Theoretical study of thermal properties and thermal conductivities of crystals
(ProQuest, 2008).
163. K. Umemoto and R. Wentzcovitch, PNAS 105, 6526 (2008).
107
164. D. Schiferl, W. Denner, H. Schulz, and W. Holzapfel, J. Appl. Phys. 49, 4411 (1978).
165. J. Wachtman, W. Tefft, D. Lam, and R. Stinchfield, J. Am. Ceram. Soc. 43, 334 (1960).
166. A. Schauer, Can. J. Phys. 43, 523 (1965).
167. G. White and R. Roberts, High Temperatures. High Pressures(Print) 15, 321 (1983).
168. P. Aldebert and J. Traverse, High Temp-High Press 16, 127 (1984).
169. A. Amatuni, T. Malyutina, V. Chekhovskoi, and V. Petukhov, High Temp-High Press
8, 565 (1976).
108
Appendices
109
Appendix A
Finite Difference Method For Harmonic Force Constants
A.1 Symmetry Analysis
In crystal solid models, many of the atomic displacements are identical. As a result,
computational time can be greatly reduced by carrying out first principle calculations only
for the truly independent atomic displacements (irreducible moves). After making a choice
of the sizes of the unit and super cell models, a group symmetry analysis is then performed to
determine the atomic displacements that are truly independent. This information is stored in
a file called IrreducibleMoveList. The symmetry analysis also determines the dependency
relations between the dependent partial forces and the independent ones and also reveals
which terms are zero. The symmetry analysis information is stored in the following files:
Phielement depNZ.dat, Phielement indNZ.dat and Phielement zeroes.dat, which contain
the dependency information, indices of independent elements , and the zero terms in the force
constants matrix, respectively(NZ stands for ?non-zero?). The MgO calculation reported in
this thesis, was based on 2-atom face centered cubic (f.c.c.) unit cell and 128-atom cubic
super cell models. After the group symmetry analysis, we find that out of the 128 atoms in
the super cell only 2 are truly independent, the rest are mirror images of these two under
group symmetry operations. According to our finite displacement algorithm, we need to
displace each of the independent atoms by a small but finite amount(?) in two directions
(+?,??). This results in a total of four-independent moves i.e. we need to perform four
total energy calculations to obtain one full dynamic matrix as opposed to 256 calculations
if we had to displace all the 128 atoms in the supercell. Consequently, our computational
110
time is reduced by a factor of 64 due to crystal symmetry alone. The symmetry analysis also
reveals that force constants matrix of MgO has 52 non-zero independent elements, 48524
non-zero elements that are linearly dependent on the independent elements and a total of
25344 elements which are zeroes. We do not need to spend any computational time on the
zero terms. Before performing a symmetry analysis to determine which atoms we need to
displace, it is necessary to relax the ionic positions in the supercell to their ground state
values. In the MgO calculation the relaxation was performed using VASP code with the
following parameters:
? NSW = 999
? ISIF = 2
? IBRION = 2
? EDIFF = 1E-9
? LREAL=.TRUE.
? ISMEAR = 0
For a 128-atom cubic MgO supercell, it takes about one hour to complete the relaxation
step. After relaxation, the One advantage of this step is that the 4 computations can be
done concurrently if there are enough compute nodes, in which case it takes about one hour
to complete the relaxation of the displaced calculate the partial forces on all the atoms due
to the 4 displacements.
A.2 Calculation of Harmonic Force Constants
The next step is to calculate the partial forces resulting from the displacement of the
unique atoms in the super cell according to the information in the IrreducibleMoveList
111
file. We use a code called displacePair.x to displace the desired atoms in the super cell. For
total energy calculations, we use the VASP code130,161. The VASP code requires the following
inputs for a total energy calculation: POTCAR (wave functions, basis set) INCAR (Energy
minimization parameters), POSCAR (displaced super cell) and KPOINTS (for Brillouin
Zone sampling). One distinct advantage of this step is that it is highly parallelised, making
it possible to perform all four total energy calculations concurrently if there are enough
free compute nodes. For more complex systems with a large number of unique atomic
displacements, sometimes we do not have enough compute nodes to run all computations
concurrently. In such situations we submit several computations (jobs) per node and let
them run one after the other. we use the portable batch system (PBS) queue software to
manage job submission deletion from the computing cluster. The H-F forces are determined
from the calculated total energies according to the Hellmann-Feynman theorem (equation
2.40) and used to evaluate the harmonic force constants. The algorithm for calculating the
partial forces and force constants matrix elements when a single atom is displaced by a small,
yet finite displacement (?) is described here162
112
Appendix B
Interpolation of Irreducible Inter-atomic Potentials (IIIP)
B.1 Introduction
We propose a new algorithm to interpolate the calculated irreducible inter-atomic
force constants over the volumes of the simulation unit cells. This technique has several
potential advantages. First, the interpolation of the irreducible harmonic force constants
over the densities/volumes can serve as a useful guide in choosing the optimal value of
atomic displacement (delta) in our real space super cell finite displacement algorithm. A
too small value will lead to higher numerical error while very large delta would lead to more
anharmonicity. Secondly the interpolation could potentially reveal the level of numerical
uncertainties associated with our first-principles calculation thereby helping us control the
quality of our data. We describe the IIIP technique next: The first step is to perform
a first-principles total energy calculation to obtain the forces on different atoms based on
different values of atomic displacement (?). Next, the set of inter-atomic force constants
for each ? are interpolated over the volumes of the simulation unit cells used in the total
energy calculations. By comparing the smoothness of the inter-atomic force constant curves,
it is easy to establish which choice of ? is optimal or converged. Additionally, numerical
uncertainties can be eliminated by fitting the calculated inter-atomic potentials up to 3rd
order Birch-Murnaghan EOS. One distinct benefit of this BM EOS fitting approach is the
predictive capability it has over the entire range of validity of the original calculation.
This process is described in section B.3. Once we establish the trend of volume/density
dependence of the calculated harmonic force constants, we can use this dependence to derive
113
harmonic force constants for any arbitrary density/volume within the range of validity of the
original first principles calculation without the need to perform new total energy calculations.
The interpolation process is described fully in the next section below with MgO as a case
study.
B.2 Interpolation of Independent Force Constants
After extracting the independent elements of the force constants matrix at several unit
cell volumes, we perform an interpolation of the elements over the volumes. In our MgO
study, we extracted the independent force constants matrix elements at 20 different unit cell
volumes between 10 ?A3 and 20 ?A3 using 4 different atomic displacements: (0.02, 0.04, 0.06,
0.08) ?A. Figure B.1 below sows an interpolation of the O?2-O?2 interatomic force constant
from calculations based on four different finite displacements (?). The results show that
an atomic displacement of 0.02 ?A or less would be too small as reflected by the irregular
behavior of the force constant with volume, a clear deviation from the trend. It can be seen
from this figure that we have a converged result for atomic displacements 0.04 , 0.06 & 0.08
?A. For our real space super cell displacement calculations for MgO, we adopt an atomic
displacement value of 0.04 ?A.
114
s49s48 s49s50 s49s52 s49s54 s49s56 s50s48
s53
s49s48
s49s53
s50s48
s50s53
s51s48
s51s53
s52s48
s52s53
s53s48
s53s53
s54s48
s54s53
s32s32 s32s61s32s48s46s48s50s32s197
s32 s32s61s32s48s46s48s52s32s197
s32 s32s61s32s48s46s48s54s32s197
s32 s32s61s32s48s46s48s56s32s197
s40
s101
s86
s47
s197
s50
s41
s86s111s108s117s109s101s47s117s110s105s116s32s99s101s108s108s32s40s197
s51
s41
s77s103s79
s70s105s110s105s116s101s32s68s105s115s112s108s97s99s101s109s101s110s116
Figure B.1: Volume dependence of one inter-atomic Force constant calculated using different
finite atomic displacements (? = 0.02, 0.04, 0.06, 0.08)?A
B.3 Volume Dependence of Harmonic Force Constants ?ij
Fromour study of MgO, SrO and CaO, ?-Al2O3 we have established that the interatomic
force constants show a robust volume dependence according to Birch-Murnaghan equation
of state up to third order as shown in equation B.1
?(V)ij = a+bx +cx2 +dx3 (B.1)
Where x = V?23, V is the unit cell volume; a, b, c and d are parameters to be determined
by fitting the calculated harmonic force constants to the unit cell volumes. For second order
BM EOS fitting, d=0. These coefficients can be used to predict ?ij for any arbitrary volume
within the range of the calculation from which the parameters were determined.
115
From this robust volume dependence of ?ij we observe that besides determining which
atomic displacement (?) is appropriate for our RSSFD method, the interpolation of ?ij
over the unit cell volumes has three distinct advantages. The first distinct advantage is
that by establishing a robust volume dependence of the harmonic force constants, one can
predict the independent force constants for other intermediate unit cell volumes without
need for a new VASP calculation. This is particularly advantageous for systems with lower
crystal symmetry, requiring longer computational times, to completely determine the H-F
forces needed to extract the harmonic force constants for the full dynamic matrix. In such
cases one can perform total energy calculations for a few unit cell volume sizes and then
use the established volume dependence of the force constants (equation B.1) to generate the
independent force constants for any arbitrary intermediate unit cell volumes within the range
of validity of original the first principle calculation. This hypothesis was tested in our MgO
study, by performing total energy calculations to directly determine the force constants for
unit cell volumes 10 ?A- 20 ?A at steps of 2 ?A. From the volume dependence, we predicted the
force constants for the other unit cell volumes 11, 13, 15, 17 and 19 ?A. We then performed
total energy calculations to determine the force constants of MgO at these same unit cell
volumes from first principles. The two results are almost indistinguishable. In figure B.2
we show the results of directly calculated independent force constants compared with those
predicted from the volume dependence. Table B.1 shows the fitting parameters for the
52 independent harmonic force constants for our MgO calculation obtained by fitting our
directly calculated force constants at five unit cell volumes (12 ?A to 20 ?A) to equation B.1.
The calculation was done using a finite displacement of ? = 0.04?A. The first and second
columns represent the indices of the atom pairs (including direction of displacement: x, y or
z) corresponding to the independent interatomic force constants having fitting parameters
a0, a1, a2, a3 (columns 3-6). One can use these fitting parameters in equation B.1 to generate
?ii,jj for any arbitrary volume between 12 ?A and 20 ?A. According to our notation, the 128
116
s54 s56 s49s48 s49s50 s49s52 s49s54 s49s56 s50s48 s50s50 s50s52 s50s54
s45s50s48
s48
s50s48
s52s48
s54s48
s56s48
s49s48s48
s49s50s48
s49s52s48
s49s54s48
s49s56s48
s50s48s48
s50s50s48
s45s50s48
s48
s50s48
s52s48
s54s48
s56s48
s49s48s48
s49s50s48
s49s52s48
s49s54s48
s49s56s48
s50s48s48
s50s50s48
s32
s32
s79
s45s50
s45s79
s45s50
s32s80s114s101s100s105s99s116s101s100s32s102s114s111s109s32s51
s114s100
s32s111s114s100s101s114s32s66s77s32s69s79s83s32s102s105s116s116s105s110s103
s32
s79
s45s50
s45s79
s45s50
s32s102s114s111s109s32s100s105s114s101s99s116s32s86s65s83s80s32s99s97s108s99s117s108s97s116s105s111s110
s79
s45
s50
s45
s79
s45
s50
s40
s101
s86
s47
s197
s50
s41
s86s111s108s117s109s101s47s117s110s105s116s32s99s101s108s108s32s40s197
s51
s41
s77s103s79 s61s48s46s48s52s32s197s32
Figure B.2: Discrete symbols represent directly calculated values of irreducible interatomic
force constants, solid lines are the values of force constants predicted using the fitting
parameters from the calculated values according to equation B.1. Predicted force constants
coincide with those from direct calculation.
117
atoms in our supercell are indexed from 1-128 with 1-64 being O?2 and 65-128 being Mg2+.
Atoms can be displaced in the directions (?x,?y,?z). we label these directions idir and
jdir taking values 1,2,3. ii and jj in table B.1 are related to the atom index ia, ja (1-128)
as shown in equation B.2.
ii = (ia?1)?Ns +idir (B.2)
jj = (ja?1)?Ns +jdir (B.3)
Ns is the number of degrees of freedom (Ns=3 in our MgO calculation) The second advantage
of this volume/density interpolation of harmonic force constants is that from the volume
dependence of the force constants one can directly evaluate the Mode Gr?uneisen parameter
(?i) via the analytic logarithmic volume derivative of the force constant elements without
requiring two different force constant matrices to evaluate (?i) via finite difference.
??q ,i ? ?dln?(
?q,i)
dlnV = ?
1
2?2(?q,i)
angbracketleftbige
i(
?q)vextendsinglevextendsingledD(
?q)
dlnV
vextendsinglevextendsinglee
i(
?q)angbracketrightbig (B.4)
d(D(vectorq))
dlnV =
summationdisplay
l
parenleftbiggd(?
ij)
dlnV
parenrightbigg?exp[?ivectorq ?(vectorR
jl ? vectorRi0)]?
mimj (B.5)
The term d(?ij(V ))dlnV can be obtained by performing a logarithmic volume derivative of
equation B.1 as follows: d(?ij(V ))dlnV = V d(?ij(V ))dV = ?23 (bx + 2cx2 + 3dx3), where x = V?23.
In our previous approach namely the finite difference approach (FDA) for ? , two force
constant matrices corresponding to two different unit cell volumes, were necessary to evaluate
(?i) for any given unit cell volume as illustrated in equation B.6.
?i(V) = ?V?2
i
summationdisplay
??
parenleftbiggD
2 ?D1
V2 ?V1
parenrightbigg
??
?e?i (?)?ei (?) (B.6)
118
Here V1 < V < V 2. Additionally, (?i) obtained this way involves taking the ratios of two
differences. When the difference in the unit cell volumes is very small compared to the
difference in the force constant matrices, this leads to divergence of (?i) which is difficult to
distinguish from numerical error. This problem is does not exist in our new approach, as
we compute ? by directly evaluating the logarithmic derivative of the analytic form of the
inter-atomic force constants. Details of the derivation of (?i) from the third order lattice
anharmonicity are reported in section 2.5.
The third distinct advantage of our interpolation method is that one can use the volume
dependence in equation B.1 to weed out numerical errors from calculated harmonic force
constants. The volume/density dependence of all the irreducible force constants of three
B1-structured ionic oxides: MgO, CaO and SrO is shown in the figures B.3, B.4 and B.5.
In table B.1, we present the fitting parameters a0,a1anda2 for all the 52 IFCs obtained
by fitting our calculated MgO IFCs to third order BM EOS. One can use these parameters
to predict IFCs for any arbitrary unit cell volume as follows: ?(V)ij = a + bx + cx2 + dx3
with x = V?23, V being the unit cell volume in ?A3.
119
s50s48 s50s50 s50s52 s50s54 s50s56 s51s48
s45s53
s48
s53
s49s48
s49s53
s50s48
s50s53
s51s48
s51s53
s52s48
s45s53
s48
s53
s49s48
s49s53
s50s48
s50s53
s51s48
s51s53
s52s48
s32
s32
s70
s111
s114
s99
s101
s32
s67
s111
s110
s115
s116
s97
s110
s116
s115
s40
s101
s86
s47
s197
s50
s41
s86s111s108s117s109s101s47s117s110s105s116s32s99s101s108s108s32s40s197
s51
s41
s32
s79s45s79
s32
s79s45s67s97
s32
s67s97s45s67s97
s86
s101s113
s67s97s79s32s32s32s71s71s65s32s99s97s108s99s117s108s97s116s105s111s110 s32s61s32s48s46s48s52s32s197
s32s32
s50s48 s50s50 s50s52 s50s54 s50s56 s51s48
s45s50
s45s49
s48
s49
s32
s32
Figure B.3: Main: The dominant irreducible harmonic force constants corresponding
to inter-atomic forces ?Ca+2?Ca+2,?O?2?Ca+2,?O?2?O?2. Inset: Minor inter-atomic force
constants (IFC). The major inter-atomic force constants are clearly volume dependent
implying that the associated phonon frequencies are volume dependent. The mode Gr?uneisen
parameter gives the volume dependence according to QHA theory. Most minor IFCs are
volume independent.
120
s50s56 s51s48 s51s50 s51s52 s51s54 s51s56 s52s48
s45s52
s45s50
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s49s56
s50s48
s45s52
s45s50
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s49s56
s50s48
s32
s32
s32
s79s45s79
s32
s79s45s83s114
s32
s83s114s45s83s114
s70
s111
s114
s99
s101
s32
s67
s111
s110
s115
s116
s97
s110
s116
s115
s40
s101
s86
s47
s197
s50
s41
s86s111s108s117s109s101s47s117s110s105s116s32s99s101s108s108s40s197
s51
s41
s83s114s79s32s32s71s71s65s32s32s99s97s108s99s117s108s97s116s105s111s110
s32s61s32s48s46s49s50s32s197s32
s50s56 s51s48 s51s50 s51s52 s51s54 s51s56 s52s48
s45s50
s48
s50
s52
s54
s56
s49s48
s49s50
s32
s32
Figure B.4: Volume dependence of Force constants in SrO. Discrete symbols in the main
graph are force constants between Sr-Sr, O-Sr and O-O. The inset shows the density
dependence of the remaining IFCs
121
s49s48 s49s50 s49s52 s49s54 s49s56 s50s48
s45s50s48
s48
s50s48
s52s48
s54s48
s56s48
s49s48s48
s45s50s48
s48
s50s48
s52s48
s54s48
s56s48
s49s48s48
s32
s32
s79s45s79
s32
s79s45s77s103
s32
s77s103s45s77s103
s70
s111
s114
s99
s101
s32
s67
s111
s110
s115
s116
s97
s110
s116
s115
s40
s101
s86
s47
s197
s50
s41
s86s111s108s117s109s101s47s117s110s105s116s32s99s101s108s108s40s197
s51
s41
s86
s101s113
s77s103s79s32s32s71s71s65s32s32s99s97s108s99s117s108s97s116s105s111s110
s32s61s32s48s46s48s52s32s197s32
s49s50 s49s52 s49s54 s49s56 s50s48 s50s50
s45s51
s45s50
s45s49
s48
s49
s32
s32
Figure B.5: Volume dependence of Force constants in MgO. Main: IFCs between Mg-Mg,
O-Mg and O-O. Inset: Rest of the IFCs.
122
Table B.1: Fitting parameters a0, a1, a2, a3 for the 52 independent harmonic force constants
of MgO fitted to third order according to equation B.1. The original LDA calculated
harmonic force constants were extracted based on a 128-atom supercell and a finite atomic
displacement ? =0.04 ?A. ii = (ia?1)?Ns +idir,
jj = (ja?1)?Ns+jdir. Where (ia, ja)= 1 to 128 and Ns=3
ii jj a0 (eV?A?2) a1 (eV) a2 (eV?A2) a3 (eV?A4)
1 1 92.47391711 - 1718.68492426 9513.01862554 - 10341.25717161
1 4 - 0.40482956 6.49643673 - 50.47378240 - 35.91076121
1 5 - 0.17783410 0.30949676 - 36.67019948 - 43.94750904
1 7 0.53994533 - 10.76994845 53.72063476 - 107.40040433
1 8 0.37080575 - 7.53578988 31.35209757 - 65.40606807
1 16 0.01730697 - 0.49509185 - 7.81719061 21.56255900
1 17 0.13912270 - 2.51110079 6.74094948 - 7.21171552
1 19 - 0.05807758 0.98499042 - 5.47536188 10.67920890
1 21 - 0.05531149 1.11998260 - 2.78379867 2.45922988
1 22 - 1.22095350 22.96027588 - 113.37737647 231.52908879
1 31 0.34524746 - 6.20819702 39.39422950 - 77.29218029
1 64 0.16080048 - 3.04080758 14.23910972 - 24.40120860
1 70 0.99249115 - 17.84161631 114.09655148 - 220.87478657
1 79 - 0.50588086 9.02857233 - 51.75492914 95.66906824
1 82 0.28496069 - 5.23341210 26.09410055 - 48.12514085
1 84 0.14895321 - 2.70672078 11.64119048 - 20.46924927
1 94 1.51639353 - 28.53442043 152.16348598 - 300.70442071
1 127 0.68312471 - 12.33400267 68.40260069 - 125.81115843
1 193 0.00193383 - 0.79754313 7.08808208 - 21.91981947
1 194 - 0.45313432 8.47513272 - 25.01964235 27.05598007
1 196 - 0.11916618 1.87574547 - 4.71566224 6.88718954
1 197 - 0.19734462 3.83688393 - 12.86162054 19.16890443
1 202 0.40561829 - 4.47310995 - 180.15834192 848.42384002
1 208 - 0.06743555 1.10193891 - 6.02138866 10.28276693
1 211 0.17313127 - 3.08862312 13.15758171 - 24.99260162
1 214 - 0.37560994 6.78677966 - 42.39902043 85.27290586
1 216 - 0.33670298 8.30745970 - 33.05358760 52.30120608
1 223 - 0.00007024 0.14532545 3.24479382 - 8.87989198
1 224 - 0.06653759 1.38235819 - 4.92958327 7.32881241
1 226 - 0.61500325 11.89659663 - 50.50666734 85.62592642
1 228 - 0.30876852 5.72822889 - 24.97879821 48.41965778
1 238 - 43.18835776 804.52823771 - 4054.14648323 3410.42072467
1 259 - 0.23835309 4.51025336 - 23.89790863 41.51155741
1 271 - 0.60784654 11.47865189 - 58.30196984 103.42874279
193 193 77.71621195 - 1443.35035960 7783.51846977 - 6646.35810740
193 196 1.28008731 - 21.86394086 124.70053340 - 420.54143451
193 197 1.95716047 - 38.57045061 144.07149558 - 360.73823814
193 199 0.30112093 - 4.43579562 20.44040415 - 42.14864279
193 200 0.23494728 - 5.86283241 21.43306968 - 35.01362786
193 208 0.20159417 - 3.13639796 6.76555682 - 11.54059049
193 209 0.20895303 - 4.18091948 14.54034306 - 21.96807984
193 211 - 0.00253489 0.06808854 0.05610270 - 1.69263644
193 213 - 0.04190171 0.39104627 - 2.42241567 7.81401863
193 214 - 0.72810237 11.62803953 - 34.78894406 103.15658186
193 223 0.00514154 - 0.12667549 4.87967830 - 9.53843389
193 256 0.08221283 - 1.33803224 3.58785005 - 7.55114810
193 262 0.56285035 - 11.29421051 67.71731699 - 112.13662209
193 271 - 0.08423031 1.07591453 - 7.55637038 16.43573280
193 274 0.29945488 - 6.19568117 29.59797576 - 51.39155161
193 276 0.15848508 - 2.79811492 11.14573352 - 18.49262242
193 286 0.18605395 - 16.68427019 120.61658739 - 359.72292314
193 319 0.40464050 - 7.95578245 41.07855666 - 72.20951864
123
Appendix C
Thermal Equation of State of Al2O3
Static EOS
The present calculation is a repeat of previous calculation done by Bin Xu a former
member of our research group119. The previous calculation by Bin, is well reproduced to
within 3 % and our new data is consistent with previous calculations and experiments.
Table C.1: Third-order BM-EOS parameters for ??Al2O3
Source V0 (?A3/atom) K (GPa) K?
Static calculation
LDA+PAW (this work) 8.354 256.5 4.200
Calculation71 8.441 258.9 4.01
Calculation126 8.10 248 4.13
Calculation73 8.486 252.6 4.237
300 K (QHA)
LDA+PAW (this work) 8.526 235.06 4.000
Calculation163 8.498 251.0 4.04
Experiment164 8.484 254.4 4.275
Thermal EOS
Table C.2: Gibbs Free energy (F0), Volume per unit cell(V0), Bulk modulus (B0 ) and
derivative of bulk modulus (B?0) at 0 GPa
T(K) F0(eV/atom) V0(?A3) B0 (GPa) B?0
300 -8.18302100 85.2618400 235.06238000 4.00336200
500 -8.21569000 85.6318600 228.98451400 4.09717000
1000 -8.36735800 86.8357200 211.67306800 4.35688500
1500 -8.58576600 88.2734500 192.45940900 4.64781800
2000 -8.85050700 89.9599700 171.34529700 4.98039700
124
Figure C.2 shows our LDA claculated thermal EOS derived within the quasi-harmonic
approximation (QHA), together with measurements from experiments165?169 at 0 GPa and
300 K
s48 s49s48 s50s48 s51s48 s52s48 s53s48
s55s53s46s48
s55s55s46s53
s56s48s46s48
s56s50s46s53
s56s53s46s48
s56s55s46s53
s57s48s46s48
s32s51s48s48s75
s32s53s48s48s75
s32s49s48s48s48s75
s32s49s53s48s48s75
s32s50s48s48s48s75
s80s114s101s115s115s117s114s101s32s40s71s80s97s41
s85
s110
s105
s116
s32
s99
s101
s108
s108
s32
s86
s111
s108
s117
s109
s101
s32
s40
s197
s51
s41
s55s46s53s48
s55s46s55s53
s56s46s48s48
s56s46s50s53
s56s46s53s48
s54s48 s54s53 s55s48 s55s53 s56s48 s56s53 s57s48
s45s56s53s46s48
s45s56s50s46s53
s45s56s48s46s48
s45s55s55s46s53
s45s55s53s46s48
s45s55s50s46s53
s69
s40
s101
s86
s41
s85s110s105s116s32s67s101s108s108s32s86s111s108s117s109s101s32s40s197
s51
s41
Figure C.1: Main: Volume-pressure relationship in ??Alumina at five temperature points:
300K, 500K, 100K, 1500K, 2000K. Inset: Static EOS fitted to 3rd order Birch Murnaghan
EOS
125
s48 s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48
s48s46s48
s48s46s53
s49s46s48
s49s46s53
s50s46s48
s50s46s53
s51s46s48
s51s46s53
s52s46s48
s32s32s84s104s101s111s114s121s32s40s80s65s87s41s32s66s105s110s32s88s117s32s40s50s48s48s57s41s32
s69s120s112s101s114s105s109s101s110s116
s32s32s87s97s99s104s116s109s97s110s32s101s116s32s97s108s40s49s57s54s50s41
s32s32s83s99s104s97s117s101s114s40s49s57s54s53s41
s32s32s65s109s97s116s117s110s105s32s101s116s32s97s108s32s40s49s57s55s54s41
s32s32s65s108s100s101s114s98s101s114s116s32s38s32s84s114s97s118s101s114s115s101s40s49s57s56s52s41
s32s32s70s105s113s117s101s116s32s101s116s32s97s108s32s40s49s57s57s57s41
s32s32s87s104s105s116s101s32s97s110s100s32s82s111s98s101s114s116s115s32s40s49s57s56s51s41
s32
s40
s49
s48
s45
s53
s75
s45
s49
s41
s84s32s40s75s41
s32
s32
s80s61s48s71s80s97
s45s65s108
s50
s79
s51
s112s97s119
Figure C.2: Thermal EOS from first-principles theoretical calculation(LDA +PAW)
compared with measurements. Discreet symbols represent measurements done by Wachtman
et al., 1962; Schauer, 1965; Amatuni et al., 1976; Adelber & Traverse, 1984; Fiquet et al.,
1999 and White & Roberts, 1983. Red solid line represents a reproduction of LDA+PAW
calculation previously done by Bin Xu, 2009(former member of our research team)
126
Appendix D
Equations of State (EOS) PV-MgSiO3 and Lattice Dynamics of PPV-MgSiO3
D.1 Static Equation of State
Total energies at fixed unit cell volumes ranging from 122 ?A3 to 160 ?A3 were calculated
using the VASP package with PAW method. An 8?8?6 Monkhorst-Park k-points grid was
used to smaple the first brillouin zone. The following energy cutoff parameters were used in
the INCAR file.
? System = PV INCAR static
? NSW = 999
? ISIF = 4
? ENCUT = 600.00
? EDIFF = 1E-9
? LWAVE= .FALSE.
? LCHARG= .FALSE.
127
Table D.1: Static EOS parameters of MgSiO3 PV and PPV fitted to 3rd order BM EOS at
ambient conditions
V0 (cm3 per mol) B0(GPa) B?0
PV(This work: LDA) 24.15 248.56 3.84
PPV(This work: LDA) 24.08 227.75 4.0
PV(Theory92:LDA) 24.71 246.1 4.0
PPV(Theory92: LDA) 24.66 215 4.41
s54s46s48 s54s46s53 s55s46s48 s55s46s53 s56s46s48 s56s46s53
s45s55s46s56
s45s55s46s54
s45s55s46s52
s45s55s46s50
s32s80s80s86s45s77s103s83s105s79
s51
s32s80s86s45s77s103s83s105s79
s51
s32
s69
s40
s101
s86
s47
s97
s116
s111
s109
s41
s86s111s108s117s109s101s32s40s197
s51
s47s97s116s111s109s41
s83s116s97s116s105s99s32s69s113s117s97s116s105s111s110s32s111s102s32s83s116s97s116s101
s53s48 s54s48 s55s48 s56s48 s57s48 s49s48s48 s49s49s48 s49s50s48 s49s51s48 s49s52s48 s49s53s48
s45s48s46s48s51
s45s48s46s48s50
s45s48s46s48s49
s48s46s48s48
s48s46s48s49
s48s46s48s50
s48s46s48s51
s72
s80s86
s60s32s72
s80s80s86
s72
s80s86
s62s32s72
s80s80s86
s72
s112
s118
s45
s72
s112
s112
s118
s40
s101
s86
s41
s80s114s101s115s115s117s114s101s32s40s71s80s97s41
s49s48s49s32s71s80s97
Figure D.1: Main: Static EOS of MgSiO3 fitted to 3rd order Birch-Murnaghan EOS. Inset:
Difference in enthalpies of PV and PPV phases plotted as a function of pressure. Our
calculation predicts that at ambient temperature, the PV?PPV transition occurs at 101
GPa which is consistent other theoretical predictions for example Tsuchiya et al. (2004).
128
D.2 Thermal EOS of PV-MgSiO3
Table D.2: LDA calculated thermal EOS parameters for PV-MgSiO3 at 0 GPa. Values in
parentheses are from Tsuchiya et al. (2005)
T(K) V0(cm3/mol K0(GPa) K?0
1.00 24.31 249.36317700 3.92435100
300.00 24.38 (24.71) 243.63925300 (246.1) 3.97535500(4.00)
1000.00 24.95 216.84633000 4.14579900
1500.00 25.46 196.36412500 4.27543300
2000.00 26.03 175.07605500 4.41773400
2500.00 26.76 152.94810800 4.57829400
3000.00 27.60 129.90054500 4.76629500
3500.00 28.68 105.80279200 4.99892500
Figure D.2 shows our LDA calculated P-V-T relations together with previous
calculations and experiment for comparison. Our results agree well with previous
calculations. In our calculation of lattice thermal conductivity, we need the pressure
dependence of simulation cell volumes, to be able convert our calculated isochoric thermal
conductivity to isobaric conditions for comparison with experiments, since experiments are
usually carried out under isobaric conditions. Figure D.3 presents our calculated thermal
expansivity together with previous calculations and experiment. Our LDA calculated
thermal EOS is consistent with both previous calculations and experiment.
129
s49s51s48 s49s52s48 s49s53s48
s48
s50s48
s52s48
s54s48
s56s48
s49s48s48
s49s50s48
s49s52s48
s32s32s84s104s101s111s114s121s58s32s79s103s97s110s111s118s32s101s116s32s97s108s46s32s40s50s48s48s48s41
s32s69s120s112s116s58s32s75s110s105s116s116s108s101s32s101s116s32s97s108s46s32s40s49s57s56s55s41
s84s104s105s115s32s119s111s114s107s58s32s76s68s65
s32s32s51s48s48s75
s32s32s49s48s48s48s75
s32s32s50s48s48s48s75
s32s32s51s48s48s48s75
s32s32s52s48s48s48s75
s80
s114
s101
s115
s115
s117s114
s101
s40
s71
s80
s97
s41
s85s110s105s116s32s67s101s108s108s32s86s111s108s117s109s101s32s40s197
s51
s41
Figure D.2: P-V-T isotherms in PV-MgSiO3 derived within the QHA. Solid lines represent
our LDA calculation, dashed lines represent previous calculation by Oganov et al. (2000),
circles denote Experimental measurement by Knittle et al. (1987)
s48 s53s48s48 s49s48s48s48 s49s53s48s48 s50s48s48s48 s50s53s48s48 s51s48s48s48
s48
s49
s50
s51
s52
s53
s54
s55
s48
s49
s50
s51
s52
s53
s54
s55
s32
s32
s32s32s84s104s101s111s114s121s32s40s84s104s105s115s32s119s111s114s107s41
s32s32s84s104s101s111s114s121s58s32s75s97s114s107s105s32s101s116s32s97s108s46s32s32s40s50s48s48s48s41
s32s32s68s105s114s101s99s116s32s69s120s112s116s58s32s87s97s110s103s32s101s116s32s97s108s32s40s49s57s57s52s41
s32
s32
s40
s49
s48
s45
s53
s75
s45
s49
s41
s84s40s75s41
s48s32s71s80s97
s53s48s32s71s80s97
s50s53s32s71s80s97
Figure D.3: Thermal EOS of PV-MgSiO3 derived within the QHA. Solid lines represent our
LDA calculation, dashed lines represent previous calculation by Karki et al. (2000), Discrete
symbols are experimental measurement by Wang et al. (1994)
130
D.3 Lattice Dynamics of PPV-MgSiO3
Using standard first-principles total energy and force calculations we have calculated
the harmonic force constant matrices and third order lattice anharmonicity tensors of
PPV-MgSiO3 using an efficient real space super-cell finite-difference algorithm. Phonon
spectra derived from the force constant matrices agree well with experiment and previous
calculations.
60 70 80Volume/Unit Cell
-20
0
20
40
60
Hamornic Force Constants
Figure D.4: Density dependence of harmonic force constants (eV/?A2) in PPV-MgSiO3.
Independent force constants shown for 160-atom periodic supercells of PPV.
131
s48
s50s48s48
s52s48s48
s54s48s48
s56s48s48
s49s48s48s48
s49s50s48s48
s49s52s48s48
s32
s90
s80
s104
s111
s110
s111
s110
s32
s32
s70
s114
s101
s113
s117
s101
s110
s99
s105
s101
s115
s32
s40
s99
s109
s45
s49
s41
s32
s90 s84 s83
s82s89
s32
s80s80s86s32s45s32s40s77s103s83s105s79
s51
s41s32s32s76s68s65s32s99s97s108s99s117s108s97s116s105s111s110
s49s49s56s71s80s97
Figure D.5: LDA calculated phonon mode frequencies in PPV-MgSiO3 at 118 GPa
s48
s50
s52
s54
s56
s49s48
s49s50
s49s52
s49s54
s32
s90
s77
s111
s100
s101
s32
s71
s114
s111
s117
s112
s32
s86
s101
s108
s111
s99
s105
s116
s105
s101
s115
s40
s75
s109
s47
s115
s41
s32
s90 s84 s83
s82s89
s80s80s86s32s45s32s40s77s103s83s105s79
s51
s41s32s32s76s68s65s32s99s97s108s99s117s108s97s116s105s111s110
s49s49s56s71s80s97
FigureD.6: LDAcalculated phononmodegroupvelocities inPPV-MgSiO3 at118GPa. Mode
group velocities in PPV are of the same order of magnitude as those in PV at corresponding
pressure.
132
55 60 65 70Volume/Unit Cell
-50
0
50
100
A ijk
Figure D.7: Density dependence of 3rd order lattice anharmonicity tensor elements (eV/?A3)
in PPV-MgSiO3. Most third order lattice anharmonicity tensors are insensitive to the change
in volume/density. A few terms show a decreasing trend with compression.
133
Appendix E
Density Dependence of Phonon Lifetimes in PV-MgSiO3
The volume dependence of the lifetimes of the 60 phonon modes at all the 100 irreducible
directions (K-points) of the BZ are shown in figures E.1 to E.17. To our knowledge this
is the first explicit ab initio calculation of density dependence of anharmonicity induced
phonon scattering rates. Different modes are observed to have different density/pressure
dependencies. Lifetimes of most phonon modes are seen to increase with decreasing volume
(increasing pressure). Overall our calculation reveals a weak pressure dependence of phonon
lifetimes. The variation of the phonon mode lifetimes with volume is almost linear. Table
E.1 shows the Volumes and thermal EOS-derived pressures at 300 K for which we directly
calculated the phonon lifetimes.
Table E.1: Simulation cell volumes and corresponding pressures for PV-MgSiO3 at 300 K.
V0 (?A3) 122 127 137 147 157
P (GPa) 119.68 95.26 57.10 29.45 9.21
134
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #1
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #2
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #3
120 130 140 150
Unit Cell Volume
0
2
4
6
8
10
12
14
16
18
20
Phonon Lifetimes(ps)
Kpoint #4
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #5
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #6
Figure E.1: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #1-6 in BZ
135
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #7
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #8
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #9
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #10
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #11
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #12
Figure E.2: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #7-12 in BZ
136
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #13
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #14
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #15
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #16
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #17
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #18
Figure E.3: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #13-18 in BZ
137
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #19
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #20
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #21
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #22
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #23
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #24
Figure E.4: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #19-24 in BZ
138
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #25
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #26
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #27
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #28
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #29
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #30
Figure E.5: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #25-30 in BZ
139
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #31
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #32
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #33
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #34
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #35
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #36
Figure E.6: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #31-36 in BZ
140
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #37
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #38
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #39
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #40
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #41
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #42
Figure E.7: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #37-42 in BZ
141
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #43
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #44
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #45
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #46
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #47
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #48
Figure E.8: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #43-48 in BZ
142
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #49
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #50
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #51
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #52
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #53
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #54
Figure E.9: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #49-54 in BZ
143
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #55
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #56
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #57
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #58
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #59
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #60
Figure E.10: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #55-60 in BZ
144
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #61
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #62
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #63
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #64
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #65
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #66
Figure E.11: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #61-66 in BZ
145
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #67
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #68
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #69
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #70
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #71
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #72
Figure E.12: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #67-72 in BZ
146
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #73
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #74
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #75
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #76
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #77
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #78
Figure E.13: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #73-78 in BZ
147
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #79
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #80
120 130 140 150
Unit Cell Volume
0
2
4
6
8
10
12
Phonon Lifetimes(ps)
Kpoint #81
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #82
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #83
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #84
Figure E.14: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #79-84 in BZ
148
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #85
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #86
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #87
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #88
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #89
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #90
Figure E.15: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #85-90 in BZ
149
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #91
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #92
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #93
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #94
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #95
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #96
Figure E.16: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #91-96 in BZ
150
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #97
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #98
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #99
120 130 140 150
Unit Cell Volume
0
2
4
6
8
Phonon Lifetimes(ps)
Kpoint #100
Figure E.17: Density dependence of phonon lifetimes for 60 phonon modes at irreducible
K-points #97-100 in BZ
151