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