THEORETICAL ANALYSIS OF ELECTRONIC, THERMAL, AND MECHANICAL

PROPERTIES IN GALLIUM OXIDE

By

Marco Domenico Santia

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Electrical Engineering – Doctor of Philosophy

2019

ABSTRACT

THEORETICAL ANALYSIS OF ELECTRONIC, THERMAL, AND MECHANICAL

PROPERTIES IN GALLIUM OXIDE

By

Marco Domenico Santia

In recent years, Ga2O3 has proven to be a promising semiconductor candidate for a wide
array of power electronics and optoelectronics devices due to its wide bandgap, high break-

down voltage, and growth potential. However, the material suffers from a very low thermal

conductivity and subsequent self-heating issues. Additionally, the complexity of the crystal

structure coupled with the lack of empirical data, has restricted the predictive power of mod-

elling material properties using traditional methods. The objective of this dissertation is to

provide a detailed theoretical characterization of material properties in the wide bandgap
semiconductor Ga2O3 using first-principles methods requiring no empirical inputs. Lattice
thermal conductivity of bulk β − Ga2O3 is predicted using a combination of first-principles
determined harmonic and anharmonic force constants within a Boltzmann transport formal-

ism that reveal a distinct anisotropy and strong contribution to thermal conduction from

optical phonon modes. Additionally, the quasiharmonic approximation is utilized to estimate

volumetric effects such as the anisotropic thermal expansion.

To evaluate the efficacy of heat removal from β − Ga2O3 material, the thermal bound-
ary conductance is computed within a variance-reduced Monte-Carlo framework utilizing

first-principles determined phonon-phonon scattering rates for layered structures containing
chromium or titanium as an adhesive layer between a β − Ga2O3 substrate and Au con-
tact. The effect of the adhesive layer improves the overall thermal boundary conductance

significantly with the maximum value found using a 5 nm layer of chromium, exceeding the

more traditional titanium adhesive layers by a factor of 2. This indicates the potential of

heatsink-based thermal management as an effective solution to the self-heating issue.

Additionally, this dissertation provides a detailed characterization of the effect of strain
on fundamental material properties of β−Ga2O3 . Due to the highly anisotropic nature of the
crystal, the effect strain can have on electronic, mechanical, and optical properties is largely

unknown. Using the quasi-static formalism within a DFT framework and the stress-strain

approach, the effect of strain can be evaluated and combined with the anisotropic thermal

expansion to incorporate an accurate temperature dependence. It is found that the elastic

stiffness constants do not vary significantly with temperature. The computed anisotropy

is unique and differs significantly from similar monoclinic crystal structures, indicating the

important role of the polyhedral linkage to the reported anisotropy in material properties.

Lastly, the dependence of the dielectric function with respect to strain is evaluated using a

modified stress-strain approach. This elasto-optic, or photoelastic, effect is found to be sig-

nificant for sheared crystal configurations. This opens up a potential unexplored application
space for Ga2O3 as an acousto-optic modulation device.

Copyright by
MARCO DOMENICO SANTIA
2019

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

CHAPTER 1

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Ultra Wide Bandgap Semiconductors . . . . . . . . . . . . . . . . . . . . . .
1.2 Gallium Oxide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Beta Phase
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Alpha Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Summary of Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 2 ELECTRONIC STRUCTURE . . . . . . . . . . . . . . . . . . . . .
2.1 Atomic Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Band Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Effective Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
4
4
6
8

11
12
16
20

25
27
29
32
32
33
35
36
37
38
39
40
41
46
49
50
50
51
52
55
57
58
63
64
67

3.3 β − Ga2O3

CHAPTER 3 THERMAL TRANSPORT . . . . . . . . . . . . . . . . . . . . . . .
3.1 Boltzmann Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Phonon Boltzmann-Transport Equation . . . . . . . . . . . . . . . . .
3.2 Methods of Computing Interatomic Force-Constants . . . . . . . . . . . . . .
3.2.1 Finite-Displacement Supercell Approach . . . . . . . . . . . . . . . .
3.2.2 Non-analytic Corrections . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 Other Scattering Mechanisms
. . . . . . . . . . . . . . . . . . . . . .
3.2.4 Gaussian Smearing . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.5 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Phonon Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Thermal Conductivity . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.4 Memory Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.5 Choice of Functional
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Phonon Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 Thermal Conductivity . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.4
Isotope Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.1 Variance-Reduced Monte-Carlo . . . . . . . . . . . . . . . . . . . . .
3.5.2 Diffuse Mismatch Model
. . . . . . . . . . . . . . . . . . . . . . . . .
3.5.3
Source-Drain Thermal Conductance . . . . . . . . . . . . . . . . . . .
3.5.4 Au/Ti/Ga2O3 Thermal Contact . . . . . . . . . . . . . . . . . . . . .

3.4 α − Ga2O3

Scattering Rates

Scattering Rates

3.5 Layered Structures

v

3.5.5 Au/Cr/Ga2O3 Thermal Contact . . . . . . . . . . . . . . . . . . . . .
3.5.6 Temperature Dependence
. . . . . . . . . . . . . . . . . . . . . . . .
3.5.7 Role of Electron-Phonon Interaction . . . . . . . . . . . . . . . . . .

68
69
70

CHAPTER 4 THERMODYNAMICS AND MECHANICAL PROPERTIES . . . . .
72
4.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
4.1.1 Elastic stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
4.1.2 Quasi-harmonic approximation . . . . . . . . . . . . . . . . . . . . .
76
4.1.3 Microscopic vs Macroscopic Gruneisen . . . . . . . . . . . . . . . . .
79
4.1.4 Linear Thermal Expansion . . . . . . . . . . . . . . . . . . . . . . . .
80
4.1.5
Isothermal elastic stiffness . . . . . . . . . . . . . . . . . . . . . . . .
82
4.1.6
Isentropic elastic stiffness . . . . . . . . . . . . . . . . . . . . . . . . .
82
4.1.7 Elastic moduli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
4.1.8
Sound velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
84
4.2 Quasiharmonic Approximation . . . . . . . . . . . . . . . . . . . . . . . . . .
85
4.2.1 Mechanical and elastic properties . . . . . . . . . . . . . . . . . . . .
88
4.2.2 Elastic Moduli α − Ga2O3
. . . . . . . . . . . . . . . . . . . . . . .
96
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
98
4.3.1 Dielectric Function . . . . . . . . . . . . . . . . . . . . . . . . . . . .
99
4.3.2 Normal Strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.3.3
Shearing Strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.3.4 Photo-elastic Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

4.3 Optical Properties

CHAPTER 5 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

vi

LIST OF TABLES

Table 1.1: Figure-of-Merit (FOM) numbers for a selection of primary semiconductor
candidates used in electronic devices. All FOMs normalized to Ge. Data
collected from [2].

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 2.1: Lattice parameters for β− Ga2O3 used in DFT calculations for LDA and

GGA compared with experimental values . . . . . . . . . . . . . . . . . .

Table 2.2: Lattice parameters in the hexagonal phase for α − Ga2O3 used in DFT
calculations for LDA and GGA compared with experimental values. All
values provided in Angstroms.

. . . . . . . . . . . . . . . . . . . . . . . .

Table 2.3: Comparison to experimental bandgaps of α−Ga2O3 [4] and β−Ga2O3 [1]

using both GGA and hybrid functional

. . . . . . . . . . . . . . . . . . .

Table 2.4: Computed effective mass values along crystallographic directions using
GGA and hybrid functionals. Kane non-parabolicity parameters are also
computed. Measured data for β − Ga2O3 taken from [19]. . . . . . . . . .

Table 3.1: Born effective charges and high-frequency dielectric tensor determined

from density functional perturbation theory (DFPT) for β − Ga2O3

. . .

2

13

16

19

23

35

Table 3.2: Percentage of modal contributions to thermal conductivity for 300 and
1050 K. There is significant anisotropy in these contributions with the
[010] having a significantly larger dependence on optical modes. The
composition does not change significantly as temperature is increased.

. .

45

Table 3.3: Characteristic MFP, Λ0, indicating the length scale above which each

phonon mode contributes to κ.

. . . . . . . . . . . . . . . . . . . . . . . .

Table 3.4: Born effective charges and high-frequency dielectric tensor determined

from density functional perturbation theory (DFPT) for α − Ga2O3

. . .

46

51

Table 4.1: Calculated lattice parameters (Å) at 0 K and bulk moduli, B (GPa), at
T = 0, 300, 1000 K as calculated within the QHA. Results using LDA
and GGA functionals as well as the measured values of Ref. [73] are listed.

86
Table 4.2: Comparison of calculated elastic stiffness constants Cij (GPa) of β− Ga2O3 . 89

vii

Table 4.3: Calculated equilibrium values for bulk modulus B, shear modulus G,
Young’s modulus (E), Vicker’s hardness (HV
), and Poisson’s ratio (ν).
All values determined by Cij and are reported in GPa except the dimen-
sionless ν. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 4.4: Sound velocities in β − Ga2O3 along different crystallographic directions
for the transverse (vt), fast mixed mode (v+), and slow mixed mode
(v−). Velocities are given in units of ×105 cm/s. Comparisons to the
lattice dynamical results extracted directly from the phonon dispersion
are given in parentheses.

. . . . . . . . . . . . . . . . . . . . . . . . . . .

95

95

Table 4.5: Calculated bulk sound velocity and computed Debye temperature (ΘD
).
) and transverse mode (vt) modes were averaged
(K). Velocities are given in ×105 cm/s.

The longitudinal (vL
(vavg) and used to compute ΘD

.

96

Table 4.6: Calculated elastic stiffness (GPa) constants using LDA compared to

α−Al2O3

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

viii

LIST OF FIGURES

Figure 1.1: Application space for various frequencies of operation for common wide

bandgap materials

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.2: Timeline of crystal growth wafer sizes achieved for primary semiconduc-

tor materials used in device engineering.

. . . . . . . . . . . . . . . . . .

Figure 1.3: (left) Conventional monoclinic unit cell of β − Ga2O3 containing 20
atoms compared to the reduced primitive cell (right) containing 10
atoms. Note there exist two types of gallium atoms corresponding to
the tetrahedral coordination and octahedral coordination which simul-
taneously exist in both unit cells.

. . . . . . . . . . . . . . . . . . . . . .

Figure 1.4: Crystal structure of conventional unit cell of α − Ga2O3
Figure 1.5: MIST Epitaxy of α−Ga2O3 layers. Growth occurs on sapphire substrate

and then is lifted-off [5] . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . .

Figure 2.1: Brillouin zone for monoclinic β − Ga2O3 (spacegroup C2/m) with high

symmetry points labelled . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.2: Rhombohedral primitive unit cell for α− Ga2O3 consisting of 10 atoms.

All gallium cations in this phase are octahedrally coordinated.

. . . . .
Figure 2.3: Rhombohedral Brillouin zone for α − Ga2O3 with high symmetry points
Figure 2.4: Electronic bandstructure for β − Ga2O3 from DFT calculations (GGA) .
Figure 2.5: Electronic bandstructure for α − Ga2O3 from DFT calculations (GGA) .
Figure 2.6: Curvature effective mass (left) is defined by the parabolicity of the band
while the transport effective mass (right) depends on the gradient of the
dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.7: Polynomial fits of energy of the lowest conduction band near CBM com-

pared to DFT data for β − Ga2O3 . . . . . . . . . . . . . . . . . . . . . .

Figure 2.8: Transport mass of β − Ga2O3 as a function of energy from the band

edge along the primary crystallographic directions . . . . . . . . . . . . .

ix

2

3

5

7

8

14

15

16

17

18

22

23

24

Figure 3.1: Umklapp and normal phonon-phonon scattering processes. The latter
do not contribute to thermal resistance directly, and only effect scatter-
ing. They are improperly considered as resistive in the RTA and can
lead to a poor estimation of thermal conductivity in materials such as
Diamond which possess a large number of normal scattering events.

. . .

Figure 3.2: Orientation of crystal axes used in this work for monoclinic Ga2O3 ,

spacegroup C2/m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.3: Phonon dispersion of Ga2O3 using a real-space finite-displacement ap-
proach. Perturbative calculations are used to resolve LO/TO splittings.
Γ →A1, A2, and A3 corresponds to propagation in the directions of the
conventional unit cell basis vectors.

. . . . . . . . . . . . . . . . . . . . .

28

38

40

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.4: (a) Calculated and published experimental thermal conductivities vs.
temperature for various crystal orientations. (b) The fitting function
κ=AT−m is plotted with a dashed line (and parameters) for each κ
tensor component.

42
Figure 3.5: Total anharmonic scattering rate for β − Ga2O3 at series of temperatures. 43
Figure 3.6: Dimensionless phase space volume available for each mode (j) and irre-
ducible q-point (q) to participate in three-phonon scattering processes[38].
Above the acoustic-optical phonon gap energy the available phase space
is small, leading to reduced optical phonon scattering. . . . . . . . . . . .

44

Figure 3.7: Accumulated thermal conductivity in three crystallographic directions
for each mode as Λmax is increased. Fits to a parametric function are
shown by the dashed lines. . . . . . . . . . . . . . . . . . . . . . . . . . .

47

Figure 3.8: Convergence of Iterative and RTA methods with increasing mesh density
Figure 3.9: Calculated phonon dispersion of α − Ga2O3 .
51
Figure 3.10: Total anharmonic scattering rate for α − Ga2O3 at series of temperatures. 52
Figure 3.11: Calculated thermal conductivity of α − Ga2O3 . Unlike the monoclinic

. . . . . . . . . . . . . . .

counterpart, the axes are nearly isotropic . . . . . . . . . . . . . . . . . .

48

53

Figure 3.12: Effective thin-film thermal conductivities for in-plane and cross-plane
transport with respect to the c-axis of α − Ga2O3 . Specularity was set
to p=0.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54
Figure 3.13: Isotopic scattering rates using natural isotopic distributions for β − Ga2O3 55

x

Figure 3.14: Isotopic scattering rates using natural isotopic distributions for α − Ga2O3 56
Figure 3.15: Thermal conductivity of α − Ga2O3 and β − Ga2O3 using both nautral

and pure isotope mixtures. . . . . . . . . . . . . . . . . . . . . . . . . . .

57

Figure 3.16: Overview of the distribution function of a computational particle and
the simplification of the deviational particle approach. The blue shaded
area represents the actual deviation from equilibrium which is simu-
lated in VRMC. Only a very small deviation from the Bose-Einstein is
common in problem where the thermal gradient is small.

. . . . . . . .

Figure 3.17: Overview of lifetime of deviational particle in the variance-reduced MC
scheme. Orange arrows represent the initial deviations imposed by a
source intensity function. Green shows the interface scattering mech-
anisms, pink the reservoir absorption and blue an intrinsic scattering
event where no interface effects are necessary to include.

. . . . . . . . .

Figure 3.18: Ohmic contact of Au/Ti on Si-implanted Ga2O3
Figure 3.19: Debye phonon density of states comparison for β − Ga2O3 and the in-
terface materials. Chromium and Titanium both contain a significantly
improved overlap in the density of states, leading to increased elastic
scattering and phonon transmission at the interface. . . . . . . . . . . . .

. . . . . . . . . . . . .

Figure 3.20: Thermal contact for Au/β − Ga2O3 using a wedge adhesive. This il-
lustrates the experimental setup of [42] this work compares to. In this
work, successive simulations are run for the range of layer thickness
presented in the wedge of adhesive Ti/Cr.

. . . . . . . . . . . . . . . . .

Figure 3.21: VRMC results for Ti/Ga2O3 interface with a Ga2O3 layer thickness
of 400 nm and an increasing layer thickness for titanium. Results are
compared to FDTR measurements of [42]. VRMC results under predict
measured data moderately and follow a similar trend in layer thickness,
with an increasing conductance until around 4 nm Ti. . . . . . . . . . . .

Figure 3.22: VRMC results for Cr/Ga2O3 interface with a Ga2O3 layer thickness of
400 nm and an increasing layer thickness for chromium. Results are com-
pared to FDTR measurements of [42] as well the Cr/Al2O3 measurements
for the same structure [49]. VRMC results under predict measured data
moderately and reach a similar plateau value at 5 nm. The peak in the
measured Cr/Ga2O3 results is predicted to be from amorphous layer
formation[42]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

61

64

65

66

67

69

xi

Figure 3.23: Thermal boundary conductance for fixed layer spacing as a function of

temperature. GaN Measurements taken from [50] and β−Ga2O3 measurements
70
from [42] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.1: Free-energy curves computed within QHA for β−Ga2O3 and minimized

to determine the Gibbs free energy function (red line).

. . . . . . . . . .

Figure 4.2: Bulk modulus computed within QHA for β − Ga2O3 . At 0 K the bulk

modulus value is B0 = 166.5 GPa . . . . . . . . . . . . . . . . . . . . . .

Figure 4.3: Average value of Gruneisen parameter of β − Ga2O3 computed macro-
scopically through the QHA method compared to the microscopic for-
mulation utilizing 3rd order IFCs . . . . . . . . . . . . . . . . . . . . . .
Figure 4.4: Frequency-dependent average mode Gruneisen parameters in β−Ga2O3 for
. . . . . . . . . . . .

a sequence of temperatures computed within QHA.

Figure 4.5: Linear thermal expansion coefficient for LDA (solid lines) and GGA
(dashed lines) calculations corresponding to applied strains εi as a func-
tion of temperature. Due to symmetry, ε4 and ε6 are zero.

. . . . . . . .

Figure 4.6: Heat capacity for β − Ga2O3 at constant pressure (blue) and constant
volume (red) computed within the QHA compared to recently measured
data of [20].

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Figure 4.7: Crystallographic orientations of β − Ga2O3 planes . . . . . . . . . . . . .
Figure 4.8: Tensor projection of Youngs modulus of β − Ga2O3
. . . . . . . . . . .
Figure 4.9: Monoclinic β−Ga2O3 polyhedral linkage illustrating the unique bonding
mechanism between octahedrally coordinated Gallium atoms. Edges of
the octahedrally coordinated Gallium atoms’ polyhedrons are shared.

.

78

79

81

82

86

88

90

91

92

Figure 4.10: Monoclinic BiMnO3 polyhedral linkage contrasts the behavior seen in
β − Ga2O3 regardless of the material possessing the same crystal type.
There are no shared edges among the octahedrally coordinated cations’
polyhedron which is predicted to be the cause of the anisotropy difference. 93

Figure 4.11: Tensor projection of Youngs modulus of Monoclinic BiMnO3 illustrat-
ing the significantly lower anisotropy compared to similar monoclinic
structure β − Ga2O3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.12: Isothermal (solid) and isentropic (dashed) elastic stiffness constants for

LDA (left) and GGA (right) functionals for β − Ga2O3 .

. . . . . . . .

93

94

xii

Figure 4.13: Temperature dependence of bulk moduli determined from computed

elastic stiffness constants up to 1400 K. . . . . . . . . . . . . . . . . . . .

Figure 4.14: Tensor projection of Youngs modulus of Rhombohedral α − Ga2O3 .
Significantly lower anisotropy is seen compared to the monoclinic phase
β − Ga2O3

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

98

Figure 4.15: Computed Real (left) and imaginary (left) components of dielectric func-

tion of β − Ga2O3 without external strains . . . . . . . . . . . . . . . . . 101

Figure 4.16: Dielectric function (real and imaginary) of β − Ga2O3 with normal ex-

ternal strains in the (cid:126)x (top), (cid:126)y (middle) and (cid:126)z (bottom) directions.

. . . 102

Figure 4.17: Dielectric function of β − Ga2O3 with applied shearing strains (solid
lines) along with the unstrained calculation (dashed lines). There is
significantly increased anisotropy when strain is applied, increasing the
brirefringence. The strongly fluctuating anisotropy with respect to fre-
quency also indicates dichroism in the refractive index.

. . . . . . . . . . 104

xiii

CHAPTER 1

INTRODUCTION

In recent decades, there has been a tremendous surge of interest in the development of

semiconducting material for use in a wide range of applications where silicon is not ideal. The

majority of these applications consist of optoelectronics and high power electronics, where

nanoscale devices are often implemented in integrated environments. Due to the limited

performance of traditional semiconductors in these environments, the focus has shifted to

wider band gap materials. The intrinsic electrical, thermal, and optical properties of wide

bandgap materials provide several advantages for the engineering and minituarization of

cutting edge electronic devices.

1.1 Ultra Wide Bandgap Semiconductors

The most popular materials for use in power electronic devices have been wide band

gap (WBG) semiconductors such as the nitrides (GaN, AlN, InN, etc.) as well as SiC and

ZnO. Wide band gaps are generally preferred for electronics applications as the large energy

separation between valence and conduction band allow for higher temperatures of operation

at higher voltages. Additionally, they can be used for many optical applications due to

their large band gap which renders them transparent to visible light and often near-UV

as well. While the traditional WBG semiconductors exhibit excellent properties in many

instances they have several limitations in practice as well as in cost-effective fabrication [1].

This coupled with the increasing desire to find maximally large bandgap semiconductors has

facilitated the development of ultra-wide bandgap materials which surpass the maximum

breakdown voltages of even the nitrides. Leading this effort has been diamond, AlGaN,
cubic BN, and Ga2O3 , the latter being the primary focus of this work. These materials
surpass in many cases a 5 eV bandgap and can potentially allow for significantly improved

performance for high power devices. An overview of the power and frequency ranges each

1

Figure 1.1: Application space for various frequencies of operation for common wide bandgap
materials

Figure-of-Merit
Johnson
Keyes
Baliga
Huang (chip-area)
Huang (high-temp)
Huang (material)

Ge
1
1
1
1
1
1

Si GaAs
7
5
1
3
7
106
18
4
0.07
0.35
2
6

4H-SiC GaN
83
5
6220
340
0.01
18

63
7
2429
192
0.01
13

AlN Diamond Ga2O3
200
375
17
1
24615
459736
1109
6053
0.00
0.00
76
22

167
20
176282
2462
0.00
72

Table 1.1: Figure-of-Merit (FOM) numbers for a selection of primary semiconductor candi-
dates used in electronic devices. All FOMs normalized to Ge. Data collected from [2].

material may provide optimal performance in, is presented in Figure 1.1

This trend is further characterized by one of many possible ’figure-of-merit’ (FOM) pa-

rameters that can be used to estimate a semiconductors performance for a specific applica-

tion. These FOM numbers are used to predict the future performance of a material in a

new class of device technologies. The FOM numbers for the most prevalent candidates for

semiconductor materials in future device designs are presented in Table 1.1. The focus of
the work presented in this dissertation is on the last column, or Ga2O3 based devices. This

2

Figure 1.2: Timeline of crystal growth wafer sizes achieved for primary semiconductor ma-
terials used in device engineering.

material excels especially in the Baliga FOM, which characterizes the ability to minimize

conduction losses and is expressed by,

BFOM = 
µE3
G

(1.1)

The material claims a BFOM that is a factor of 3443 improved from silicon for example and
a factor of 4 improvement over GaN. The primary advantage to Ga2O3 as a semiconductor
over the other high-performing materials in the table such as Diamond and AlN, is the

potential for commercial scalable growth beyond what any other wide bandgap material
has achieved. Although the use of Ga2O3 in electronic devices only began within the last
couple of decades, the progress made in developing larger wafers is unmatched. This trend

is visualized in Figure 1.2. This compares favorably to other UWBG semiconductors such as

Diamond which is generally considered to possess the most desirable material properties for

electronics application. However, diamond still suffers from very limited growth potential

with the largest wafers available only spanning a matter of millimeters even after many

decades of research.

3

1.2 Gallium Oxide

Ga2O3 has long been studied due to its close relationship with sapphire and other poly-
morphs of aluminum oxide. Interest in the material among the electronic device community

has increased significantly in recent years, due to the wide array of applications. Some

of these include power electronic devices, solar-blind UV photodetectors, gas sensors, and

transparent conducting films of optoelectronic devices. Many of the applications stem from

the wide bandgap and large window of optical transparency. The successful development
of Ga2O3 transistors, such as MOSFETs and MESFETs [2], coupled with the low-cost and
scalable growth using melt growth techniques has led to a surge of interest in this material.

The primary phase of interest has been the thermodynamically most stable phase, which
is the monoclinic β phase. However there are actually five commonly identified polymorphs,
labeled α, β, γ, δ and 
. In this dissertation, I will focus primarily on the β phase but will also
begin to explore the α phase as it has recently had success as an electronic device material
It also has exhibited some growth advantages over β − Ga2O3 ,
in a metastable state.
primarily owing to its higher degree of symmetry as will be discussed later. In both phases,

the potential application space is large and the material has found success in nearly all the

cases outlined in the previous section. Generally speaking, most of the critically important

semiconductor device properties scale approximately exponentially with band gap and thus
the 4.9eV band gap of Ga2O3 can find useful applications nearly anywhere that GaN or
GaAs could. For example, Ga2O3 is optically transparent up to 250nm while remaining
electrically conductive so it has found use as a window on optical devices [2]. In this section,
I will briefly review the material properties of the α and β phases.

1.2.1 Beta Phase

Comparing to other commonly used wide bandgap (WBG) semiconductor materials, β −
Ga2O3 has a relatively complex crystal structure with low symmetry. The crystal structure

4

Figure 1.3: (left) Conventional monoclinic unit cell of β − Ga2O3 containing 20 atoms com-
pared to the reduced primitive cell (right) containing 10 atoms. Note there exist two types
of gallium atoms corresponding to the tetrahedral coordination and octahedral coordination
which simultaneously exist in both unit cells.

for both the conventional unit cell and the reduced primitive cell are given in Figure 1.3. The
β phase is base centered monoclinic (space group C2/m) with 20 atoms in the conventional
cell and 10 in the primitive cell.
In this symmetry, there are 3 unique crystallographic
directions ((cid:126)a,(cid:126)b, and (cid:126)c) as well as a monoclinic angle β between the (cid:126)a and (cid:126)c axes..

Interestingly, the β phase contains two unique species of gallium possessing different
coordination numbers. There are both tetrahedral and octahedral coordinations labeled

Ga(I) and Ga(II) respectively. This, in addition to the unique crystallographic directions,

leads to anisotropy in the material properties. This anisotropy is important to characterize

and will be considered throughout this work. Because of the reduced coordination of half
the metal ions comparing to the α phase for example, which contains only octahedrally
coordinated cations, the interatomic distances are relatively large. The closest nearest-

neighbor length between two gallium cations is 3.04 Å.

To achieve electronic conduction, numerous types of dopants have been successfully used
in β − Ga2O3 to achieve N-type conduction. Most commonly Si and Sn have been good
shallow donors with small activation energies. Electron densities can be controlled between
1015 − 1019cm−3 [2]. There has not however been any reported success with p-type doping,
so hole conduction has not been realized which remains to be one of the materials largest

bottlenecks. It is not generally expected that hole conduction will be improved by doping

5

LatticeThermalConductivityinβ-Ga2O3M.D.Santia,N.Tandon,andJ.D.AlbrechtDepartmentofElectricalandComputerEngineering,MichiganStateUniversity,EastLansing,MI48824,USA(Dated:April24,2015)Latticethermalconductivityiscomputedfromfirst-principlesandBoltzmannTransportEquationcalculationsforbulkβ-Ga2O3.Forceconstantsforboththesecondandthirdorderpotentialinteractionsarecomputedwithareal-spacefinite-displacementapproach.Phonondispersionaswellasanharmonicpropertiesarethencomputedandusedtocomputebulkthermalconductivitytensorfortempertaturesrangingfrom25Kto1050K.Calculatedvaluesat300Kareκ[100]=16.06WK−1m−1,κ[010]=21.54WK−1m−1andκ[001]=21.15WK−1m−1whichareingoodagreementwithexperimentallyobservedvalues.Decomposingκintomodecontributionsrevealthatopticalphononmodescontributesignificantlytotheoverallthermalconductivity,asmuchas44%at300Kinthe[010]direction.β-Ga2O3hasreceivedsignificantattentionasasemi-conductorinrecentyearsduetoitslargebandgapandthepotentialcrystalgrowthviaamodifiedCzochralskimethod1.Thematerialisidealforhigh-powerelectronicdevicesduetoitshighelectricalbreakdownstrength.AtthenanoscaleofthesedevicesacompleteunderstandingofthethermalpropertiesofthebulkmaterialbecomesessentialformanagingtheaccumulationofJouleheat-ingincriticalregionssuchasthechannelofthedevice.β-Ga2O3hasseenrecentsuccessforremoving”hotspots”inGaNhigh-electron-mobilitytransistors(HEMT)asathin-filmgateoxide2.Understandingthermaltransportpropertiesatthenanoscaleisnecessaryfortheengineer-ingofthesedevices.Therearefewexperimentallymeasuredvaluesforκ1,3–5whichhaveobservedκinthe[100]and[010]crys-tallographicdirections,yieldingroomtemperaturevaluesofaround13WK−1m−1and21WK−1m−1respectivelyanddemonstratingsignificantanisotropy.Thislowκhaslargelybeenattributedtothephonon-phononUmk-lappdominatedscatteringhowevertheexactmechanismshaveyettobeexplored.Recentadvancesinabinitiomodellinghaveledtothedevelopmentofamethodofcomputinganharmonicforceconstantsfromfirst-principles6,7.Knowledgeofanhar-monicforceconstantsallowsforphonon-phononinter-actionstobedeterminedwithoutsignificantapproxima-tions.Itiswellknownthatthermalconductivityisadirectresultofthesethree-phononprocesses8.ThishasledtothedevelopmentofaccuratenumericalmethodsofsolvingthephononBoltzmannTransportEquation(BTE)7,9.Inthiswork,weutilizetheShengBTEpackageof9forsolvingtheBTEforlatticethermalconductivity.Thesemethodshaveproventobeveryaccurateatrepro-ducingexperimentallyobservedthermalproperties7,9,10.ForadetaileddescriptionofthetheorybehindtheBTEmethodutilizedinthiswork,wereferto9.DensityFunctionalTheory(DFT)total-energycalcu-lationsarerequiredforthecomputationofinteratomicforceconstants(IFC).Foranharmonicaswellashar-monicforceconstants,singleatomsaredisplacedinasu-percellforvariousconfigurationsbasedonthesymmetry(a) (b) FIG.1:Monoclinicβ-Ga2O3crystalstructure.VisualizationsdonewithVESTApackage15ofthecrystal.Totalenergycalculationsfromthevariousdisplacedsupercellsallowsforces,andsubsequentlyIFCs,tobecomputed.Inthiswork,DFTcalculationsareper-formedinQuantumESPRESSO11,whilePhonopy12isusedforharmonicphononpropertiesandShengBTE9isusedforsolvingtheBTEaswellasforcomputingthird-orderIFCs.β-Ga2O3isamonoclinicstructuredcrystalwithspace-groupC2/mand20atomsintheconventionalunitcell.LatticeparameterswereoriginallydeterminedbyGeller13tobe:a=12.214˚A,b=3.0731˚A,c=5.7981˚Aandβ=103.83◦.Thestructurecontainstwonon-equivalenttypesofGalliumatomsclassfiedbythenatureoftheirbonding.Thefirst,Ga(I),hasatetrahedralco-ordinationwiththesurroundingOxygenatomswhilethesecond,Ga(II),hasanoctahedralcoordination.Thecon-ventionalcellcanthenbetransformedintoaprimitivecellcontaining10atomsperunitcell14asshowninFig.1Insuchamonoclinicstructure,thereisanon-trivialconstructionoftheBrillouinzoneandthishasledtoam-biguityinpreviousliterature16–19.Weadopttheorienta-tionusedby14.DetailsonthesymmetrypathandBril-louinZoneorientationcanbefoundthesupplementarymaterial.Itisimportanttonotethattoappropriatelycomparetoexperimentalresults,thesymmetrypathin14ismodifiedtocontainreal-spacecrystallographicdirec-tionsoftheconventionalunitcell,ratherthanthatofthedue to the difficulty in finding acceptors with small activation energy given the very large
bandgap of β − Ga2O3 . This is a relatively common trend in WBG oxides such as ZnO,
MgO, Ga2O3 and Al2O3 where the conduction band states come from metal atoms and the
valence band states from O 2p orbitals. These orbitals present very large effective masses

and small dispersion. Holes have been found to prefer localization in the form of polarons

rather than in the form of free holes in the valence band [3].

The largest bottleneck β − Ga2O3 faces is the well known low thermal conductivity, or a
high thermal resistance. This large thermal resistance coupled with the high temperatures

of device operation in power electronics will lead to a localized thermal effects known as

"self-heating". While the phenomena of self-heating is not a new concept, it has become

an increasingly impassable obstacle in recent years due to the minituarization of devices.
Unfortunately, β − Ga2O3 has a thermal conductivity an order of magnitude lower some of
the other wide bandgap materials such as GaN, AlN and SiC. Comparing to another ultra

wide bandgap material like diamond, which claims the highest known thermal conductivity,
the issue can make the use of Ga2O3 very difficult in practice. This problem has gained much
attention in the device engineering community as of recently and characterizing the details of
the thermal behavior has become a crucial task for further improvements of Ga2O3 devices.
The thermal conductivity has been a difficult quantity to measure in full due to the anisotropy

of the crystal structure and subsequent material property tensors. Many precisely cut crystal

orientations of high quality must be available for measurement to resolve a highly anisotropic

tensor. This can be prohibitive due to the infancy of the crystals growth processes for
Ga2O3 and the likely existence of residual strains that may effect the material properties.

1.2.2 Alpha Phase

The α phase of Ga2O3 is of rhombohedral, corundum, symmetry and is iso-structural to
sapphire. The conventional unit cell presented in Figure 1.4. In this configuration, the oxygen

ions are hexagonally closed-packed while gallium occupies two-thirds of the octahedral sites.

6

Figure 1.4: Crystal structure of conventional unit cell of α − Ga2O3

There are two independent lattice parameters a = 4.98 Å and c = 13.43 Å. Recently, the α
phase of Ga2O3 has begun to be explored as a possible solution for limitations β−Ga2O3 faces
in device engineering. As mentioned previously, there has not been success in developing
bipolar β−Ga2O3 devices that can be manipulated from n-type to p-type. However, there has
been success in other similar oxide semiconductors for use in pn junction heterostructures.
Namely, α−Ir2O3 and α−Rh2O3, which both have exhibited p-type conductivity by the
Seebeck effect [4] and possess the same symmetry as the α phase of Ga2O3 . Additionally,
the bandgap of this phase is even larger than the β phase which leads to a larger breakdown
voltage. This has motivated and led to preliminary success at device fabrication of bipolar
α − Ga2O3 transistors [3].

Material properties of α− Ga2O3 are generally better than its β− Ga2O3 counterpart for
electronics applications. The bandgap of α − Ga2O3 is the largest of all the polymorphs of
Ga2O3 with a value of 5 eV increasing to over 5.16 eV for tin doped thin films [1]. Unlike the β
phase, α−Ga2O3 is only metastable, however it can be maintained under ambient conditions
unlike the other metastable phases of Ga2O3 . This phase can be produced by a transition of

7

Figure 1.5: MIST Epitaxy of α − Ga2O3 layers. Growth occurs on sapphire substrate and
then is lifted-off [5]

the β phase under hydrostatic pressures at higher temperature. Experimentally, the phase
transition has been observed at P = 4.4 GPa at high temperatures and at P = 19.2 GPa
at cold compression [1]. Because of this, the ability to fabricate bulk crystals is prohibitive.
All applications of α − Ga2O3 are in thin-film form, grown in epitaxial layers on a sapphire
substrate which then can be lifted off and isolated, or layered with another substrate, as

illustrated in Figure 1.5. This offers a unique possibility to control thermal heat flux by means

of more efficient heatsink thermal energy extraction. The thinner the layer and the lower

the mismatch between layer and substrate, the more likely it becomes that a high thermal
conductivity substrate will more rapidly dissipate heat. The sapphire Al2O3 substrate has
the same crystal structure with only 4% lattice mismatch.

Overall, this phase is less understood than the β phase. Although the high bandgap
is appealing, the difficult growth process can be prohibitive especially for measurements of

bulk crystalline properties. Theoretical predictions of the material properties are the only

option when the crystal cannot be easily isolated. Thus far, very few theoretical works have
been performed for any optical, thermal or electronic properties of α − Ga2O3 .

1.3 Summary of Objectives

As mentioned previously, perhaps the most appealing benefits of several polymorphs of
Ga2O3 is the potential for readily scalable bulk crystal growth by several possible methods.
This has been the primary factor leading to the massive increase in Ga2O3 research in recent
years [6]. However there is still not a clear picture of what the limitations of the material are
and how they can be overcome. This leaves many open questions about how Ga2O3 compares

8

to other WBG semiconductors and how effective it could be in practice.

Experimentally, it is difficult to characterize thermal properties of ternary or doped struc-

tures due to the need to fabricate specific structures meant solely for thermal measurements.

Additionally, a highly anisotropic material also can prove to be a challenge to accurately an-

alyze due to the necessity of several high quality crystals grown in specific orientations to re-

solve the entire thermal conductivity tensor. Fortunately, there has been rapid developments

in computational modelling of thermal transport properties in recent years. Reliable theo-

retical models that can extract spatial and temporal distributions of thermal gradients are

beginning to become computationally feasible without any empirical inputs. First-principles

methods of determining harmonic and anharmonic phonon properties in crystals systems of

all symmetry classes is now possible [7]. This dissertation will seek to apply these new meth-
ods to Ga2O3 to further understand the thermal conductivity behavior and breakdown the
origins of the high resistance. Additionally, the topic of thermal boundary conductance will

be analyzed using recently developed Monte-Carlo methods to assess the viability of forming
effective thermal contacts with Ga2O3 using higher thermal conductivity metals to reduce
the self-heating effects. At the time of writing, these methods have not been rigorously

tested for realistic device models, so this work will additionally serve to set a benchmark of

the robustness of the variance reduced Monte-Carlo algorithm.

Understanding the origins of the reportedly low thermal conductivity is an essential

tool for further refinement and management of the thermal properties. This has been a
challenge for β − Ga2O3 due to the lack of fundamental knowledge of material mechanical
and thermodynamic properties. For example, the elastic properties are mostly unknown
for all phases of Ga2O3 .
of information about the material and can be linked to heat capacity, thermal expansion,

In general, elastic properties can reveal a significant amount

phonon dispersion, Gruneisen parameters, melting points, etc. When engineering devices

operating in a high temperature environment, the characterization of thermal expansion

and overall thermodynamic stability is important to ensure reliability and efficiency. In this

9

dissertation, first-principles methods are used to evaluate the thermodynamic properties

within a quasiharmonic approximation for phonons for a series of volumes of the unit cell.

This methodology allows for a robust and accurate determination of fully anisotropic material

properties which would otherwise be difficult to measure. Additionally, the effect that strain

has on these fundamental properties can be more directly evaluated. Similarly, a significant

portion of this work will focus on the predicted optical response properties such as dielectric
function to supplement the understanding of the anisotropic properties of Ga2O3 -based opto-
electronic devices such as solar-blind photo-detectors. Using this, the effect of mechanical
strain on the system can be included and used to determine photo-elastic behavior of β −
Ga2O3 . The photo-elasticity is a difficult to measure quantity that provides a useful estimate
of the optical behavior in a material known to undergo significant strain. These strains can

appear as residual strains during crystal growth for example, and can lead to a high degree

of birefringence even in the ideal bulk crystal.

10

CHAPTER 2

ELECTRONIC STRUCTURE

Understanding and predicting the properties of semiconductors through electronic transport

requires an accurate description of the electronic band structure of the material. Histori-

cally, computation has relied on empirical models such as tight-binding and empirical pseu-

dopotential methods which utilize experimentally determined parameters. These methods

become infeasible as crystal complexity increases, as the necessary parameters are generally

unavailable and difficult to extract from measurements. The alternative, and more predic-

tive, approach is to utilize first-principles atomistic calculations to determine band struc-

ture and other transport properties. Specifically, density functional theory (DFT) using a

planewave pseudopotential basis has demonstrated high accuracy with reasonable computa-

tional requirements. Although these methods have vastly increased computational overhead

to empirical methods, they have proven themselves in recent years to be highly accurate

without any need for measured data for considerably complex systems. This is especially
important for a material as new and unrefined as Ga2O3 , where an improved understand-
ing of the electronic properties at the atomic level can aid in the engineering of transport

properties which directly influence device performance. The amount of information available

on the electronic structure of this material is surprisingly limited and the exact nature of

the bandgap is still not well understood. This is in part due to the complexity of the Bril-
louin zone which contains a very complex high symmetry path. Furthermore, β − Ga2O3 is
known to possess a significant anisotropy in both electrical and optical properties, with the

former largely attributed to an anisotropic effective mass and the latter due to band-to-band

selection rules.

Throughout this work the QuantumEspresso [8] package will be utilized for the electronic

structure calculations. The majority of the DFT calculations utilize the generalized gradient

approximation (GGA) with projector augmented wave (PAW) potentials. At times however

11

the local density approximation (LDA) is also utilized and compared to as the former is

known to overestimate the lattice constant and subsequently the cell volume. For specifi-

cally electronic structure, GGA is generally considered more accurate [8] and will be used

for the exchange-correlation functional with the PAW potentials used to represent the inter-

action between ionic cores and valence electrons. The kinetic energy cutoff for Kohn-Sham

planewave basis is 125 Ry and typical Brillouin zone integrations are performed on a 8x8x4

grid. Structural relaxation is performed on a more dense grid with 16x16x8 to ensure the

forces between atoms are properly minimized. This step is important when extending the

DFT calculation to phonons, where there is a greater sensitivity to residual forces.

2.1 Atomic Structure

In the case of β − Ga2O3 there is a mixed bonding character between the gallium ions as
mentioned previously. They occupy both tetrahedral and octahedral sites while the oxygen

ions form a distorted cubic packing structure with the unique sites comprising two threefold

coordinated and one fourfold coordinated [9]. This mixed bonding character leads to chem-

ical bonds which exhibit both covalent and ionic bonding characteristics and are heavily

directionally-dependent (anisotropic). This anisotopy has been realized geometrically and in

measurements of the electrical conductivity tensor [2].

The geometry of the β phase is more complex and can be represented by its monoclinic
conventional cell or reduced primitive cell with lattice parameters a, b, c and β which is the
angle between a and c. The lattice vectors for the conventional unit cell are expressed in the
traditional form [9] by

aconv
1
aconv
2
aconv
3

= aˆx + 0ˆy + 0ˆz

= 0ˆx + bˆy + 0ˆz

= c cos β ˆx + 0ˆy + c sin βˆz

12

and the primitive cell transformation is mapped in the following way

aprim
1
aprim
2
aprim
3

3

2

1 + aconv
)/2
1 + aconv

= (aconv
= (−aconv
= aconv

)/2

2

The detailed parameters of the crystal structures after geometric relaxation for both LDA
and GGA functionals for the β-phase is provided in Table 2.1.

Expt[10] GGA
12.41
3.07
5.86
103.72

12.23
3.04
5.80
103.82

LDA
12.26
3.04
5.81
103.75

a
b
c

β

Table 2.1: Lattice parameters for β − Ga2O3 used in DFT calculations for LDA and GGA
compared with experimental values

It is crucial to provide clarity on the shape of the Brillouin zone (BZ) as there has not been

a consistent convention adopted in past literature [9]. This is attributed to the flexibility

in how the axes of a monoclinic crystal are assigned. This topic was recently clarified in

detail in [9] where it was determined that some inconsistency has arisen due to the non

interchangeable BZs of the primitive and conventional cells. Essentially, it was found that

the smallest parallelpiped of the reciprocal lattice vectors of the conventional cell actually

contain a smaller volume than the analog in the primitive cell. This is a non-intuitive result

that leads to an irreducible BZ of the primitive cell which does not fully span reciprocal

space when performing calculations on the conventional cell. This is contrary to higher

symmetry crystal structures and here, the BZ of the primitive and conventional cell must be

considered independently of one another for consistency. In general the monoclinic system

can be represented by one of five unique BZ shapes depending on the choice of base-centering
(c-base centered or b-base centered). In this work, the b-axis is used as the dyadic axis so
it is b-base centered, leading to the spacegroup symbol C2/m. The full BZ is presented in

13

Figure 2.1: Brillouin zone for monoclinic β−Ga2O3 (spacegroup C2/m) with high symmetry
points labelled

Figure 2.1 with high symmetry points and path labelled for clarity. The high-symmetry path

used in this BZ was constructed in [9] intentionally to ensure a high symmetry path which

avoids discontinuities.

The α phase of Gallium Oxide is a simpler crystal structure possessing a trigonal sym-
metry in the form of a rhombohedral lattice, spacegroup R3c which can also be expressed
in hexagonal coordinates. The conventional unit cell can be composed as a series of six for-

mula units totalling 30 atoms of nearly closed packed ABAB stacking of the oxygen anions.

The gallium ions meanwhile occupy two-thirds of the octahedral interstitial along the [0001]

direction. The primitive unit cell can be simplified to two formula units totalling ten atoms
as shown in Figure 2.2. The crystal has the same structure as sapphire, which is the α
phase of Aluminum Oxide (Al2O3 ). This structure is well known and has been researched
extensively so there is significantly more literature to compare to than β − Ga2O3 . The
form of the lattice which is typically adopted requires the z axis to be aligned parallel to
the rhombohedral (cid:104)111(cid:105) direction which is the c-axis if using the hexagonal representation

14

Figure 2.2: Rhombohedral primitive unit cell for α − Ga2O3 consisting of 10 atoms. All
gallium cations in this phase are octahedrally coordinated.

of the rhombohdral lattice. The lattice vectors can be expressed as

√

3ahex/3, chex/3(cid:17)

a1 =(cid:16)0,
a2 =(cid:16)−ahex/2,−
a3 =(cid:16)ahex/2,−

√
√

3ahex/6, chex/3(cid:17)
3ahex/6, chex/3(cid:17)

(2.1)

In this definition, the lattice parameters can be defined in terms of the rhombohedral lattice
constant a and subsequent angle α,

ahex = a(cid:112)2(1 − cos α)
chex = a(cid:112)3(1 + 2 cos α)

The lattice parameters for this phase are provided in Table 2.2 with a corresponding Brillouin

zone given by the shape in Figure 2.3.

15

Expt.[10] GGA LDA
5.31
4.99
13.41
55.93

5.321
4.97
13.43
55.79

5.4
5.10
13.64
55.87

a
ahex
chex
α(◦)

Table 2.2: Lattice parameters in the hexagonal phase for α−Ga2O3 used in DFT calculations
for LDA and GGA compared with experimental values. All values provided in Angstroms.

Figure 2.3: Rhombohedral Brillouin zone for α − Ga2O3 with high symmetry points

2.2 Band Structure

It is well known that planewave DFT computations of electronic structure fail to re-

solve an accurate bandgap between the lowest conduction band and highest valence band.

This is due to the inability of traditional DFT to incorporate the electron number density

dependence in the exchange-correlation functional term. The details of this restriction are

beyond the scope of this work and are explained in detail in [11]. However, as the materials

of interest are semiconducting, the bandgap is a key parameter and must be kept in mind

when comparing measurements to DFT calculations. The experimental room-temperature

16

Figure 2.4: Electronic bandstructure for β − Ga2O3 from DFT calculations (GGA)

band-gap has been reported to be around 4.9 eV [12]. In general, oxides of post-transition
metals such as Ga2O3 and Al2O3 have an (n−1)d10ns2 electronic configuration with densely
packed metal ions. The bonding and non-bonding O 2p states give rise to the valence band
while the conduction band arises from the Gas-Op interactions [13].
It is these Gas-Op
interactions that lead to the band gap as well as the dispersion of the conduction band. As

a result, these are the orbitals that will participate in charge transport upon doping.

The computed bandstructure is given in Figure 2.4 along the high-symmetry path which

coincides with the BZ given previously. Notably, there is very little dispersion in the va-

lence bands which indicates a low hole mobility. Additionally, it is interesting to note the

conduction band minimum is located at Γ, but the valence band maximum is not located
at this position as was previously speculated thus causing the semiconductor to technically

be an indirect bandgap. In the work of [9] it was realized that the valence band maximum

17

-202468101214X1YΓNXMIFZZI1ΓLYΓF1Energy (eV)Figure 2.5: Electronic bandstructure for α − Ga2O3 from DFT calculations (GGA)

is actually located at the I1 point of the BZ which was not previously reported on. It does
not in fact have a significant impact on the overall behavior however due to the dispersion

of the valence bands being very small. Other valleys appear at the Z, Y, N and X points.

Given the significantly anisotropy in within the crystal axes, it is relatively surprising to note

the lack of corresponding anisotropy in this electronic structure. For this reason, it could

reasonably be predicted that the electrical transport properties express a three-dimensional

character to them and are not largely anisotropic. This property will be analyzed in further

detail using effective mass calculations.

For the α phase, similar behavior is seen to the β phase and is presented in Figure 2.5,
both show small dispersion in the valence bands and possess a conduction band minimum

(CBM) at the BZ center (Γ-point). The computed bandgap is about 1eV larger than the
β. The uppermost valence bands from the O 2p orbital are mostly flat and display little

18

dispersion. This leads to a valence band minimum which is more sensitive to the precise
arrangement of atoms. Interestingly, in the α phase, the existence of an indirect bandgap is
more pronounced with the indirect gap being off the high symmetry axes somewhere along Γ
to P a VBM 0.3 eV higher than at Γ. In both cases, the bandgap is somewhere in the range
of 2.5-3.5 eV which is lower than the experimentally determined values by a large margin.

This is expected and is a limitation of DFT without using hybrid functionals that can correct

the exchange-correlation term.

To correct for this inaccuracy, the use of hybrid functionals has been adopted in several

works [9, 14, 15] using the Heyd, Scuseria, and Ernzerhof (HSE) [16] potential however this

has consistently shown a bandgap which still underestimates the bandgap to around 4eV.

To correct the bandgap in this work, the relatively new Gau-PBE functional [17] functional

is adopted for its accuracy and computational efficiency. This functional approximates the

Coulomb term in the Hartree-Fock potential to a Gaussian function. The results for the

bandgaps compared to the GGA and experimental values is provided in Table 2.3. Excellent

exp
∆
g

5.3-5.6
4.7-4.9

∆GGA

g
3.29
2.54

gauPBE
∆
g
5.66
4.69

α − Ga2O3
β − Ga2O3

Table 2.3: Comparison to experimental bandgaps of α − Ga2O3 [4] and β − Ga2O3 [1] using
both GGA and hybrid functional

agreement is found with the experimental bandgaps in both phases when using the hybrid

gauPBE functional. It is important to note however, due to the computational complexity

of hybrid functional calculations, the ability to perform the calculation for a dense mesh in

reciprocal space is limited. Computing the full band structure along the high symmetry path

would require significant computational resources.

19

2.3 Effective Mass

Traditional electron transport computations are most commonly performed using a drift-

diffusion or hydrodynamic model. Due to the complexity of these calculations, the ability to

incorporate a full electronic band structure without any loss of detail is limited. Instead, an

analytic approximation to the conduction band minimum (CBM) and valence band maximum

(VBM) are usually used in conjunction with an effective mass approximation to analytically

represent the energy dispersion for the important valleys. This approximation allows for a

model of how an electron in a periodic potential will respond to external perturbation by

representing the electron as a free-electron with a renormalized mass, i.e the effective mass.
There is an increasing need to investigate these approximations in β − Ga2O3 where little
data is available both the electronic structure and the effective masses. This type of analysis

is most commonly performed simply by evaluating the curvature of the CBM/VBM within

a simple harmonic oscillator expression,

1
mc

=

1
¯h2

∂2E
∂k2

(2.2)

The most common model of band structure used in transport solvers is the parabolic band

model. This method assumes a fully parabolic representation of energy dispersion where the

effective mass parameter generates the entirety of the band structure about a given valley,

E(k) =

¯h2k2
2mc

(2.3)

which is often valid for semiconductors where the carrier concentration is low and the elec-

tron energy is near the band minimum. The effective mass can be measured experimentally

or determined from first-principles calculations where a full band-structure is used. This

approximation is an efficient method of computing electronic properties, but makes the fun-

damental assumption that bands are parabolic in nature. As semiconductors used in elec-

tronics becomes more sophisticated and less silicon-like, the need to consider non-parabolicity

becomes increasingly important. In a material at elevated temperatures, or with a higher

20

carrier concentrations for example, the electrons are excited to higher up in the band where

there is a distinct non-parabolic character typically. This non-parabolicity can be modelled

in a variety of different methods but generally introduces another representation of the ar-
tificial mass known as the "transport mass" mt. The transport mass will be equivalent to
the effective mass for completely parabolic bands. The simplest approach to characterizing

non-parabolicity is derived [18] with the relationship between momentum and velocity of an
electron wavepacket mv = ¯hk and gives,

1
mt

=

1
¯h2k

∂E
∂k

(2.4)

One of the more robust and common approaches however is to use an expansion in energy

which introduces non-linearity,

¯h2k2
2mt,0

= γ(E) = E + αE2 + βE3 + . . .

(2.5)

where the transport mass is expressed at the band edge, mt,0. Taking this expression to the
first order results in the Kane quasi-linear dispersion relation [18]

where the transport effective mass then can be determined through differentiation,

¯h2k2
2mt,0

= E(1 + αE)

mt(E) = mt,0(1 + 2αE)

(2.6)

(2.7)

Using this expression, the non-parabolicity of the bands can be characterized through the
parameter α. It corresponds to the flattening of the bands as the energy moves away from
the band edge. The flattening results in a larger overall transport mass. This is shown
in Figure 2.6 The non-parabolicity α will be α = 0 for a perfectly parabolic band, and
will be positive or negative for conduction and valence bands respectively. Computing the
curvature and Kane dispersion parameter is done in this work for β − Ga2O3 and compared
to the recent effective mass measurements [19]. It is important to note the effect that the

exchange-correlation term will have on the overall curvature. As previously shown, the

21

Figure 2.6: Curvature effective mass (left) is defined by the parabolicity of the band while
the transport effective mass (right) depends on the gradient of the dispersion

band gap is significantly underestimated with the LDA/GGA functionals. Here the hybrid

Gau-PBE functional is once again used to correct not just for the width of the bandgap,

but also the curvature and anisotropy. The same bandstructure previously given is utilized

and the curvature is computed along 1D slices of the BZ, i.e high-symmetry paths along

the crystallographic directions. The curvature can be computed using a finite-difference

approach or a quadratic-least squares approach. The non-parabolicity in the Kane-dispersion

is computed using a 6th order polynomial to the eigenvalues computed in DFT for up to

0.25eV from the band edge. The transport mass is then computed through differentiation

as previously noted.

It is important to sample very densely near the CBM as the finite-

difference method will amplify any numerical noise that comes with the energy eigenvalues

from DFT. For this reason the least-squares approach is generally considered more robust.
For β − Ga2O3 , energy band reconstruction was very accurate with respect to the DFT
results, this is shown for the [100] direction of the lowest conduction band in Figure 2.7. For

the non-parabolic fits, the transport masses corresponding to these fits are given in Figure 2.8

as a function of energy from the band-edge. The transport masses show relatively little

anisotropy in the conduction band which is surprising considering the geometric anisotropy

22

Figure 2.7: Polynomial fits of energy of the lowest conduction band near CBM compared to
DFT data for β − Ga2O3 .

of the crystallographic axes. The effective mass parameters computed for both β−Ga2O3 and
α − Ga2O3 are presented in Table 2.4 and compared to the measured data available.

α − Ga2O3

[100]
0.133
0.261
0.33

[010]
0.133
0.261
0.33

[001]
0.141
0.272
0.34

m∗
GGA
m∗
hybrid
α
Exp

β − Ga2O3

[100]
0.126
0.223
0.38
0.288

[010]
0.140
0.301
0.41
0.283

[001]
0.131
0.292
0.40
0.286

Table 2.4: Computed effective mass values along crystallographic directions using GGA and
hybrid functionals. Kane non-parabolicity parameters are also computed. Measured data
for β − Ga2O3 taken from [19].

It is important to note the little variance in effective mass between the two phases, this

implies the effect of the oxygen coordination on the mass is minimal. Additionally, the

anisotropy is relatively low, under 5%, implying the geometric anisotropy does not have a

significant impact on effective mass.

23

0.000.020.040.060.080.100.12k (Å1)0.000.050.100.150.20energy (eV)polynomial order 6finite diff parabolicfive point parabolicweighted parabolicKane quasi-linearDFTFigure 2.8: Transport mass of β − Ga2O3 as a function of energy from the band edge along
the primary crystallographic directions

24

CHAPTER 3

THERMAL TRANSPORT

The most prevalent bottleneck in power electronics development is often the ability to manage

the very large thermal accumulation. This issue becomes more severe as the device scales

down in size, resulting in a great challenge for engineering on-chip electronics. This self-

heating of devices can drastically reduce lifetime of the device. For example a 50 K increase

in temperature can reduce the time to failure of a device by an entire order of magnitude.

Additionally, there can often be a significant influence of the thermal build-up on electronic

properties of a system. Hence, a proper assessment of the thermal transport becomes an

essential component to predicting a materials response properties. The effect of heat on

the system can play an indirect but important role as well. For example, the complexity of
the crystal structure in Ga2O3 naturally results in anharmonicity in thermal expansion, and
this phenomena can play a significant role in the evolution of thermal stresses and growth
defects such as cracks. Growth behavior is still a rapidly evolving field for Ga2O3 and
hence thermal characteristics are of equal or greater interest in this environment, so the

need for characterization is great. This need further increases when considering the known
low thermal conductivity of Ga2O3 is considered the greatest bottleneck of the material for
engineering purposes.

There are few reports of experimentally measured values for the thermal conductivity
κ [20–22]. The majority of the data has been measured along the [100] and [010] crys-
tallographic directions, yielding room temperature values of around 13 and 21 W/mK, re-

spectively and exhibiting significant anisotropy. These low overall conductivity values have

largely been attributed to scattering dominated by phonon-phonon Umklapp processes, how-

ever the exact mechanisms have yet to be explored. In this chapter, a first principles cal-
culation of the thermal conductivity tensor (κ) is presented with a detailed analysis of the
scattering contributions by phonon mode.

25

The need for sophisticated methods of computing thermal transport properties is rela-

tively new and stems from the length-scales of modern power electronics devices and tran-

sistors. Macroscopic systems can appropriately apply the well-known Fouriers Law of heat

and set up the system as a partial differential equation in a more traditional finite-element

type framework. This style of modelling heat breaks down at the nanoscale where simulta-

neous ballistic and diffusive phonon transport can be occurring. Additionally these phonons

can interact and scatter significantly with dislocations, vacancies, interfaces, etc. The only

reasonable solution to the problem is to perform more complete solutions of the boltzmann

transport equation with detailed attention to the force-interactions of atoms and their sur-

roundings, i.e vibrational frequencies at a microscopic level.

Historically, the predictive power of computational tools for thermal transport proper-

ties involving complicated underlying mechanisms such as thermal conductivity, have been

limited. Computational models generally have relied on significant empirical data to pa-

rameterize to, and even when empirical data has been available, the underlying models

assume a degree of simplicity in the crystal structure. The low symmetry space group of
β − Ga2O3 does not adhere to the required simplicity of such models. For example the
"gray model" of phonon dispersion that assumes one acoustic phonon mode with a single

group velocity and one optical mode of zero velocity. This approximation is often utilized

and is effective for simple systems in which optical modes are known to contribute very little

thermal transport. Generally, this model proves to be surprisingly effective as most thermal

transport is within the acoustic phonon modes which possess the highest group velocities.

For a complex crystal with many atoms in the unit cell and a large number of corresponding

phonon modes however, this assumption is not guaranteed to suffice. This theoretical limi-

tation, coupled with the complexity of taking detailed measurements of anisotropic thermal
properties of ideal bulk material has hindered the predictive capabilities of Ga2O3 device
design. Using recently developed first-principles methods however, it is possible to determine

not only macroscopic thermal properties of materials but also prove a deeper insight into

26

the underlying mechanisms.

3.1 Boltzmann Equation

Thermal conduction at the nanoscale is dictated by the interaction and diffusive motion

of the phonons of the system. Expressing this transport phenomena is handled in a variety

of different ways depending on the type of carrier, the conservation laws it must adhere to,

and the interplay between the mean-free path and life-time of the particle. For the case

of particles in which the characteristic length-scale of relevant phenomena far exceeds the

mean-free path, the transport is called diffusive. In other words, the particles on average

will collide with other particles before interacting with a boundary. Quantitatively this is

characterized by the Knudsen parameter

Kn =

Λ
L

(3.1)

where Λ represents the mean-free path (MFP). Transport is considered ballistic when Kn
far exceeds 1 and carrier scattering events are negligible compared to boundary effects. For

realistic systems of phonon transport in an electronic device, the length scales tend to lead
to Knudsen parameters between Kn = 0.1 and Kn = 10. This is considered a transitional
transport regime which is the most complex to model and requires a more general transport

equation to be solved which incorporates ballistic and diffusive transport. The evolution

of such a system is described by a Boltzmann-type equation describing the diffusion and

scattering of phonons:

∂f
∂t

+ ∇kω(k, p) · ∇xf = [

∂f
∂t

] collisions

(3.2)

where the frequency of the phonon is represented by ω, the polarization (mode index) by
p, the wave-vector by k and the phonon distribution function by f. The key to any solver
chosen to model this formula is in how the right-hand side is taken into account, i.e the scat-

tering terms. Unlike the electron gas models, phonons obey Bose-Einstein statistics rather

27

Figure 3.1: Umklapp and normal phonon-phonon scattering processes. The latter do not
contribute to thermal resistance directly, and only effect scattering. They are improperly
considered as resistive in the RTA and can lead to a poor estimation of thermal conductivity
in materials such as Diamond which possess a large number of normal scattering events.

than Fermi-Dirac. This means they do not need to adhere to particle number conservation,

which leads to many forms of interactions which can occur within the phonon population in

isolation. This is also neglecting the many external scattering mechanisms which often take

place in reality, which will have varying degrees of importance. Although higher-order inter-

actions are possible, in most cases the primary interaction giving rise to thermal resistance is

the three-phonon interaction, where two phonons combine to a single phonon or one phonon

decays into two. Both of these processes conserve energy and total momentum accordingly.

These scattering events are illustrated in Figure 3.1.

These two are classified more succinctly as:

Type I (absorption) :

Type II (emission) :

(cid:48) = k
(cid:48) = ω

(cid:48)(cid:48) + H
(cid:48)(cid:48)

k + k

ω + ω

k − k
ω − ω

(cid:48) = k
(cid:48) = ω

(cid:48)(cid:48) + H
(cid:48)(cid:48)

28

where H = (cid:126)G is the reciprocal lattice vector for Umklapp processes and H = 0 for normal
processes. Umklapp processes are those which do not conserve local momentum due to the

transference of momentum to the lattice, and hence slow the phonons down. This loss of

momentum from the reference frame of the phonon gives rise to a thermal resistance, as the

slowing down of the heat-carrying phonons will lead to less rapid heat dissipation. Normal

processes on the other hand fully conserve momentum and thus do not directly provide

resistance, but they must be considered in the BTE regardless as they do change the overall

phonon occupation.

3.1.1 Phonon Boltzmann-Transport Equation

The difficulty in solving the BTE in general lies in the form of the collision operator, or the

right-hand side of the BTE. In the most general form, this could be expressed as [23],

{fk(cid:48)p(cid:48)(fkp + 1) − fkp(fk(cid:48)p(cid:48) + 1)}Qk(cid:48)p(cid:48)

kp

(cid:48)(cid:48)
{(fkp + 1)(fk(cid:48)p(cid:48) + 1)fk(cid:48)(cid:48)p(cid:48)(cid:48) − fkpfk(cid:48)p(cid:48)(fk(cid:48)(cid:48)p(cid:48)(cid:48) + 1)}Q
kp,k(cid:48)p
{(fkp + 1)fk(cid:48)p(cid:48)fk(cid:48)(cid:48)p(cid:48)(cid:48) − fkp(fk(cid:48)p(cid:48) + 1)(fk(cid:48)(cid:48)p(cid:48) + 1)}Qk(cid:48)p(cid:48),k(cid:48)(cid:48)p

kp

[

∂f
∂t

] coll = (cid:88)k(cid:48),p(cid:48)
+ (cid:88)k(cid:48)p(cid:48),k(cid:48)(cid:48)p(cid:48)(cid:48)
2 (cid:88)k(cid:48)p(cid:48),k(cid:48)p(cid:48)(cid:48)

+

1

(3.3)

(3.4)

(3.5)

where Q is defined as the transition probability matrix which is determined by the Hamil-
tonian of interaction given relevant physical conservation laws. This form leads to a signif-

icantly non-linear integro-differential equation for even the simplest of geometries and is in

most cases unrealistic to solve in full. Traditionally, the problem can be solved by assum-

ing a relaxation time approximation (RTA) is valid [23]. This assumes that all momentum

destroying processes tend toward an equilibrium Planck distribution and in many cases is

sufficient for determining thermal conductivity. However, this assumption does not accu-

rately incorporate normal scattering events and treats them as momentum non-conserving.

In situations where normal scattering events dominate the thermal conduction, such as in

low temperature environments, this approximation is not sufficient.

29

Recently, a more robust solution has been introduced which allows for a more realistic

treatment of phonon interactions [24] which utilizes a linearized form of the BTE and has

been implemented in the ShengBTE software package [24]. This approximation assumes
the temperature gradient is small enough such that the distribution function fkp
can be
expanded to first order in ∆T . The resulting form of the BTE then can be determined
from the phonon dispersion, group velocities, and relaxation time determined in an iterative

procedure. In this linearized BTE form, the difficult part of the calculation then becomes
the determination of relaxation times τ 0
λ

for a phonon mode λ,

1
τ 0
λ

=

1
N

(

+(cid:88)λ(cid:48)λ(cid:48)(cid:48)

Γ+

λλ(cid:48)λ(cid:48)(cid:48) +(cid:88)λ(cid:48)λ(cid:48)(cid:48)

1
2

λλ(cid:48)λ(cid:48)(cid:48) +(cid:88)λ(cid:48)
Γ−

Γλλ(cid:48))

(3.6)

where N represents the grid density of sampled reciprocal space in the BZ and λ has been
used to simplify notation and represent a (kp) pair with the scattering rates Γ+ and Γ− for
absorption and emission being expressed as,

Γ+
λλ(cid:48)λ(cid:48)(cid:48) =

¯hπ
4

|V +
λλ(cid:48)λ(cid:48)(cid:48)|2δ(ωλ + ωλ(cid:48) − ωλ(cid:48)(cid:48))

Γ−
λλ(cid:48)λ(cid:48)(cid:48) =

¯hπ
4

−
|V
λλ(cid:48)λ(cid:48)(cid:48)|2δ(ωλ − ωλ(cid:48) − ωλ(cid:48)(cid:48))

0

f(cid:48)
0 − f(cid:48)(cid:48)
ωλωλ(cid:48)ωλ(cid:48)(cid:48)
f(cid:48)
0 + f(cid:48)(cid:48)
0 + 1
ωλωλ(ωλ(cid:48)(cid:48)

(3.7)

(3.8)

Here, conservation laws are considered by the dirac delta functions representing Fermi’s

golden rule combined with the collision matrix elements V

±
λλ(cid:48)λ(cid:48)(cid:48)[23],

V

±

λλ(cid:48)λ(cid:48)(cid:48) = (cid:88)i∈u.c.(cid:88)j,k (cid:88)αβγ

Φαβγ
jik

p(cid:48)(cid:48)(k)

λ(i)eβ
eα

p(cid:48),±q(cid:48)(j)eγ

(cid:112)MiMjMk

Given this form of collision matrix, only the phonon eigenvectors eα
λ

interatomic force constants (IFC) are needed,

Φαβγ

ijk =

∂3E
i ∂rβ

j ∂rγ

k

∂rα

30

(3.9)

and anharmonic

(3.10)

Anharmonic IFCs describe the atomic force interactions between atoms within a crystal

lattice to a third order approximation to the interatomic potential. Analogously, a second

order approximation to the interatomic potential would yield harmonic forces such as those

found in phonons where restoring forces among atomic bonds behave harmonic oscillator

like. Computing the third-order anharmonic IFCs is non-trivial. Primarily this owes to the

fact that most empirical potential models found via force-matching methods are not fitted

to parameters which capture the relevant anharmonic behaviors. This in part is due to the

limited physical parameters which actually results from the anharmonicity in interatomic

potential. Gruneisen parameters can make a good fitting parameter however they are diffi-

cult to characterize experimentally [23]. For simple systems such as silicon there has been

some success with Stillinger-Weber potentials but not in a robust fashion [25]. There is no

theoretical restriction for the use of an empirical force-field potential in this computation,

if one could construct one that captures anharmonic force-interactions. But this remains

a challenging task where the predictive power of such models becomes questionable and

requires extensive validation. Alternatively, the ability to perform these calculations using

first-principles methods has become feasible for most systems given enough computational

resources [26].

To account for this inaccurate representation of normal modes, a full linearized BTE

can be solved for bulk material subjected to a temperature gradient which accounts for all

scattering terms. The details of this formulation are described in detail in [24]. Essentially,

the linear system takes the form

AF = B

(3.11)

where F is a solution variable for a phonon mode k, Fk = (Fk,x, Fk,y, Fk,z). They
represent a cartesian vector of generalized mean free paths which are used in the expression

31

for thermal conductivity,

This is solved starting from the simpler RTA solution Fk,α = τkvk,α

κα,β =(cid:88)k

Ckvk,αFk,β

(3.12)

as the initial guess

and then an iterative procedure is performed until the solution to this linear system is con-

verged. This method is considerably more expensive but necessary for predicting thermal

conductivity in systems possessing a large number of normal scattering events such as Dia-

mond.

3.2 Methods of Computing Interatomic Force-Constants

Recent advances in ab-initio modeling have led to rigorous methods of computing har-

monic and anharmonic interatomic force constants from first-principles energy calculations

[26–28] on a series of displaced supercell configurations. Knowledge of the anharmonic real-

space force constants allows for the phonon-phonon interactions to be computed for a given

reciprocal space mesh. Significant prior work has focused on the three-phonon scattering

processes and this has led to the development of accurate numerical methods of solving

the phonon Boltzmann Transport Equation (BTE) to compute lattice thermal conductivity

through phonon scattering.[24, 27]

3.2.1 Finite-Displacement Supercell Approach

Density Functional Theory (DFT) total energy calculations are required for the computation

of interatomic force constants (IFCs). For anharmonic as well as harmonic force constants,

single atoms are displaced in a supercell for various configurations based on the symmetry of

the crystal. Total energy calculations from the various displaced supercells allows forces, and

subsequently IFCs, to be computed. Upon displacement of the individual atoms, applied

forces on the remaining atoms in the unit cell are computed and divided by the displacement

distance to yield the force-constant. For polar semiconductors, the Born charges are also

32

added into the formulation to account for non-analytic long-range Coulombic forces. The

born charges can be determined experimentally or by means of density-functional perturba-

tion theory (DFPT). This method is sometimes known as the "direct method" [26] and is

considered an extension of the "frozen-phonon" method that allows for phonon dispersion

to be computed in arbitrary directions and points in the BZ. In this supercell approach,
the atomic displacement of the sth atom in the unit cell is displaced in the α direction by
a ±δα
, which can be tuned slightly but often is sufficient at a value of 0.015 Angstroms.
The larger the displacement, the smaller the reciprocal space mesh is required to resolve

s

the force-constants accurately, however too large of a value can lead to inaccuracies due to

higher order force-interactions becoming dominant. The matrix element of the IFC can then

be written as

su = −(cid:34) F β

φαβ

u (+δα

u (−δα
s )

s ) − F β
2δα
s

(cid:35)

(3.13)

This is then used to generate the dynamical matrix at arbitrary points in reciprocal space

[26],

Dαβ

st (q) =

1√

MsMt(cid:88)R

st (R)eiq·R
φαβ

(3.14)

this is now diagonalized to compute the square of the phonon frequencies. This method

generalizes to arbitrary order force-constants assuming proper numbers of displacements are

made. In this work, only up to third-order interactions are considered as they will contain

the majority of the relevant phenomena to thermal conductivity. The third-order force-

constants require two atomic displacements to be made at a time and higher-order would

require more. In general, fourth-order IFCs can be computed using this methodology as well,

but it becomes computational prohibitive due to the number of displacements required.

3.2.2 Non-analytic Corrections

It is important to take note of the effect that the ionicity of the crystal plays on the overall

lattice dynamical properties.

In general, the displacement of an atom in an ionic crystal

33

generates a dipole moment which can itself couple with other dipoles formed within the

crystal by means of a macroscopic electrostatic field, and subsequently effect the longitudinal

optical (LO) phonon mode energies. This electrostatic field does not however impact the

transverse optical (TO) modes. These dipole-dipole interactions have an impact on the

harmonic and anharmonic force constants in such a way that a splitting of the LO/TO

modes occurs at the Γ point which is dependent on the direction of approach to the Γ
point. These "non-analytical" contributions can be characterized by the non-vanishing Born
effective charge (BEC) tensor Z∗ in the form of a non-analytical correction to the dynamical
matrix, ˜Dαβ
st

which is written in the limit as it approaches the BZ center,

st (q → 0) =
˜Dαβ

=

∗νβ
t

MsMt(cid:80)γ Z
(q · Z∗

4πe2
√
Ω
4πe2
√
Ω
MsMt

∗γα
s

qγ(cid:80)ν Z
(cid:80)γ,ν qγ
ν∞qν
s)α(cid:0)q · Z∗
t(cid:1)β

q · 
∞ · q

qν

(3.15)

Here 
∞ is the high-frequency dielectric tensor and Z∗
is born effective charge tensor for
the sth atom. Both of these parameters are calculated in this work by density functional
perturbation theroy (DFPT) in the QuantumEspresso software package. The force-constants

s

are corrected by this non-analytical term in the phonon dispersion by the following form,

Φαβ
st (R) = φαβ

st (R) + ˜ϕαβ

st

where the correction is now determined by the condition

lim
q→0

1√

MsMt(cid:88)R

st eiq·R = ˜Dαβ
˜ϕαβ

st (q → 0)

Which results in the final expression for non-analytical correction

˜ϕαβ
st =

1
N

4πe2
Ω

(q · Z∗

s)α(cid:0)q · Z∗
t(cid:1)β

q · 
∞ · q

(3.16)

(3.17)

(3.18)

Born effective charges represent the first derivative of the macroscopic polarization. The

results determined from DFPT calculations is presented in Table 3.1

34


∞
Z∗ (Ga1)
Z∗ (Ga2)
Z∗ (O1)
Z∗ (O2)
Z∗ (O3)

xx
3.8
2.88
3.44
-2.10
-2.26
-1.96

yy
4.1
2.74
3.23
-1.45
-2.26
-2.25

zz
3.9
3.06
3.12
-2.50
-1.38
-2.29

xz
0.114
0.285
-0.345
0.067
0.205
-0.217

zx
0.114
0.155
-0.180
0.037
0.228
-0.240

Table 3.1: Born effective charges and high-frequency dielectric tensor determined from
density functional perturbation theory (DFPT) for β − Ga2O3

3.2.3 Other Scattering Mechanisms

There are potentially numerous other sources of scattering which can be present in the

crystal that can lower thermal conductivity. All additional scattering will inevitably lead to

increased thermal resistance so it is important to consider potentially important phenomena.

One such mechanism is the effect of isotopic scattering in which the impurities present due

to naturally occurring isotope composition leads to mass differences which effect the overall

thermal conductivity. This scattering rate is given by the expression,

1/τ

iso
λ =

π
2N0

ω2

λ(cid:88)κ

gκ(cid:88)λ(cid:48) (cid:12)(cid:12)(cid:12)ˆeλ

κ · ˆeλ(cid:48)∗

κ (cid:12)(cid:12)(cid:12)

2

δ(cid:0)ωλ − ωλ(cid:48)(cid:1)

(3.19)

represents the phonon eigenvector for the κ atom in mode λ. The mass variance
where ˆeλ
κ
term gκ is the Pearson deviation coefficient of the masses miκ for isotope i of atom κ and is
given by

gκ =

1
m2

k(cid:88)i

fiκ (miκ − mκ)2

(3.20)

with fiκ representing the concentration and mκ the average mass. The stable isotopes of
gallium are found in concentrations of 60.11% for 69Ga and 39.89% for 71Ga which gives a
deviation parameter of gGa = 1.97× 10−4. The oxygen isotope variance is very low and will
not contribute. The effect of isotopic scattering is important to consider because it is one

of the potentially tunable parameters for controlling thermal conductivity. The ability to

increase thermal conductivity by isotopic enrichment has recently been demonstrated [29].

35

Additionally, because this parameter is the only empirical input to the first principles ther-

mal conductivity workflow, it can be adjusted freely to attempt to characterize a measured

samples isotopic concentration. The last scattering mechanism which can be easily included

is scattering from an ideal boundary. This is done by means of a relaxation time,

τ B
λ =

Leff
vave

(3.21)

where Leff is the effective phonon mean free path which will be on the order of the sample
size when boundary scattering is relevant. The average velocity of acoustic phonons, vave is
defined by a sum of the transverse, vT

, and longitudinal velocities vL

,

vave =(cid:20)1

3(cid:18) 2

vT

+

1

vL(cid:19)(cid:21)−1

(3.22)

In both cases, the scattering rates are incorporated into the full solver by combining them

with the anharmonic scattering rate for a given point by Matthiessen’s rule,

anharm
1/τλ = 1/τ
λ

+ 1/τ

iso
λ + 1/τ

boundary
λ

(3.23)

3.2.4 Gaussian Smearing

One of the more computationally demanding and complex aspects of this workflow lies in

the Brillouin zone integration scheme when integrating properties of interacting phonons.

The reason for this complication is that generally, the three-phonon scattering process for

a point on the sampling mesh will almost always be composed of modes which are not on

the N1 × N2 × N3 q-point mesh to remain energy conserving, ωλ −(cid:0)±ωλ(cid:48) + ωλ(cid:48)(cid:48)(cid:1) = 0.

Several approaches have been utilized to attempt to solve this issue and allow for the delta

functions in the lifetime calculations to be integrated over. One of the more advanced

methods is known as the tetrahedron method [30]. This method entails decomposing a BZ

into tetrahedral regions and using a piecewise linear interpolation on the matrix elements

for integration. The method is generally very accurate however can also be expensive and

can be more difficult to achieve a full convergence. The second method, which is adopted

36

in this work, is a broadening technique which replaces the dirac delta functions with some

smearing function with a finite broadening width that can be tuned. The smearing function

is generally chosen to be a gaussian of the form:

g (ωλ − W ) =

1√
2πσ

e

−(ωλ−W )2

2σ2

(3.24)

where the width of the gaussian is controlled by the parameter σ and will approach a dirac
delta function as σ tends to zero. The quality of this approach then becomes a delicate
balance between keeping σ small enough to approximate a delta function, and keeping it large
enough such that it diminishes the probability of not completing all the scattering channels

needed for fully integrating the points on the sampling grid. This method is extended into
an adaptive smearing framework in the work of [24] where the authors introduce how the σ
parameter can be determined and adjusted automatically. Readers are referred to work of

[24] for complete details on the adaptive smearing technique.

3.2.5 Convergence

Due to the infancy of the first-principles thermal transport work-flow, it is important to

establish convergence criteria for the calculations to provide validation. Especially in the

finite-displacement super-cell approach utilized in this work, there are many parameters that

can be adjusted in the first-principles methodology. Ideally, the largest possible supercell

would be used as this is more likely to assure that all atoms in the unit-cell are considering

all force-interactions with neighboring cells until they have reached a radius where the mag-

nitude is negligible. If this were not the case, there is a risk of contributions from images in

neighboring supercells to be included in the FCs. In practice however, it is not realistic to

increase the supercell beyond what is feasible for first-principles calculations. On the other

hand, using the smallest possible supercell that would sufficiently capture force interactions

would require a higher density k-space mesh to be used for sampling reciprocal space. For

crystals with 80 atoms or greater, only sampling at the Γ point in reciprocal space is suffi-

37

Figure 3.2: Orientation of crystal axes used in this work for monoclinic Ga2O3 , spacegroup
C2/m

cient There is no single solutions for these problems and there are trade-offs with whatever

convergence criteria is chosen. This convergence criteria becomes even more sensitive for the

case of third-order force-constants where more nearest neighbor interactions also must be
converged. For monoclinic β − Ga2O3 the largest supercell that was feasible for harmonic
phonon (second-order force-constants) calculations was 3x3x3 (270 atoms) while anharmonic

(third-order force-constants) was 2x2x2. In both cases, Γ-point sampling was sufficient.
3.3 β − Ga2O3

In this section, the previously outlined methodologies are applied to β − Ga2O3 to deter-
mine a theoretical prediction of thermal conductivity from first-principles[31]. To compare

the anisotropic material properties with the available measured results, the crystallographic

orientation must be established and is illustrated in Figure 3.2.

38

For a monoclinic crystal, the second-rank material property tensor takes the form:

κij =

κxx

0

κxz

0

κyy

0

κxz

0

κzz

.



(3.25)

3.3.1 Phonon Dispersion

Firstly, structural optimization is performed using the plane-wave pseudopotential technique

implemented in [8] using the projector augmented waves method with the generalized gra-

dient approximation [32, 33]. Published lattice parameters were used as the initial structure
for performing structural optimization [34]. A plane-wave cut-off of 125 Ryd with a 16×16×8
k-point grid was used to relax the structure such that the forces on all atoms were below

1 mRyd/a.u. The Born effective charges and dielectric tensor required to resolve the LO-TO
splitting at the Γ point were evaluated using density functional perturbation theory [35].
The structural parameters are used for the calculations of the IFCs, both harmonic and
third order, which were evaluated using a 2×2×2 supercell containing 80 atoms. Due to the
10 atom primitive unit cell, there will be 30 phonon modes overall which will likely result in

significant scattering to occur as many energy and momentum conserving transitions become

possible. The computed phonon dispersion along the high symmetry path is presented in

Figure 3.3. This dispersion shows a significant anisotropy with acousto-optical gaps being

significant (approximately 15 meV) along the [001] direction which limits the possible scat-

tering channels for high energy optical modes along this direction. Meanwhile, the [100] and

[010] directions do not possess this gap but do however experience optical modes with signif-

icant group velocities. In general, optical modes are not expected to have significant group

velocities and are assumed to be effectively flat and without dispersion but this is not the

case here. This implies the optical modes may potentially carry a non-negligible amount of

heat and could contribute to the thermal conduction. In all directions, there exists relatively

steep acoustic modes with group velocities similar to GaN [20]. One distinct feature which

39

Figure 3.3: Phonon dispersion of Ga2O3 using a real-space finite-displacement approach.
Perturbative calculations are used to resolve LO/TO splittings. Γ →A1, A2, and A3 corre-
sponds to propagation in the directions of the conventional unit cell basis vectors.

can be seen near the Γ point, is the mixing of low lying optical modes with the acoustic
modes near the BZ center.

3.3.2 Thermal Conductivity

After phonon dispersion is established, the computation can be extended to the anharmonic

force-constant calculation in a similar manner outlined previously with scattering rates being

calculated subsequently for use in the phonon BTE. In the calculation of thermal conductiv-

ity, the RTA and linearized form of the BTE is solved in an iterative manner. Anharmonic

IFCs are used to compute scattering rates and the mean free path (MFP) parameters for each
mode which are then used to compute κ. For this calculation, a 11×11×11 k-space sampling
mesh is utilized for solving the BTE over temperatures ranging from 25 K to 1050 K.

The calculated values at 300K are κ[100]=16.06, κ[010]=21.54 and κ[001]=21.15 W/mK.
These values compare well with the measured values of κ[100]=13 and κ[010]=21 W/mK.[12,
20–22] The theoretical results systematically overestimate κ owing to the various minority

40

2priatelycomparetoexperimentalresults,thesymmetrypathin[14]ismodifiedtocontainreal-spacecrystallo-graphicdirectionsoftheconventionalunitcell,ratherthanthatoftheprimitiveunitcellusedtosimplifycalcu-lations.ThesevectorsaredenotedasA1,A2,andA3cor-respondingtothe[100],[010],and[001]crystallographicdirectionsrespectively.DensityFunctionalTheorycalculationsareperformedusingtheplane-wavepseudopotentialpackageQuantumESPRESSO,wherewehaveadoptedthegeneralizedgradientapproximation(GGA)inthePerdew-Burke-Ernzerhofparameterization.Forpseudopotentials,wehaveusedprojectoraugmentedwave(PAW)potentialswithanenergycutoffvalueof125Ryanda16x16x8k-pointgridwithanenergyconvergenceof1x10−12Ryfortheself-consistentDFTcalculations.Structuralparame-tersareoptimizedwithavolume-cellrelaxationmethodtoensureequilibriumconditions.Properequilibriumstructureparametersarecriticalforcalculationsinvolv-ingdisplacedsupercellsthusourforceconvergencecrite-riamustbesetsufficientlysmall.Wehavefoundthataforceconvergencecriteriaof1x10−4a.uproperlyrelaxesthestructuralparameterstoequilibrium.Convergedlatticeparametersare:a=12.276˚A,b=3.054˚A,c=5.647˚Awhichisinagreementwithexperimentalresults[13].Weadditionallyuseaperturbativecalcula-tionusingthePHononpackageofQuantumEspressotodeterminetheBornEffectiveCharges(BEC)anddielec-trictensorrequiredtoresolveLO/TOsplittings.Phonondispersioniscalculatedwitha2x2x2supercelland2x2x2k-spacesampling.Thisrequirestotal-energyDFTcalculationsfor20displacedsupercellscontaining80atomseach.Phonopythenusedtoreconstructhar-monicIFCsandgenerateharmonicpropertiessuchasdispersionwhichisshownalongthepreviouslydiscussedBrillouinzonepathinFig.2.Additionallytheprojecteddensityofstates(PDOS)iscomputedandshownalong-side2.Wehavefoundthedispersiontobeingoodagree-mentwithpreviousDFTcalculations[19,20].Anharmonicforceconstantswerethencomputedwithasimilarreal-spacefinite-displacementmethodasimple-mentedin[9].Forthisweusea2x2x2supercellwithΓpointk-spacesamplinganduptothirdnearest-neighborinteractions.Thisrequiresatotal-energyDFTcalcula-tiontobeperformedfor392displacedsupercellsof80atoms.Thecalculatedthirdorderforceconstantsfilealongwithvariousotherinputfilesusedinthisworkisincludedinthesupplementarymaterialsalongwithade-taileddescriptionofprocedureused.Finally,theBTEissolvedintheShengBTEframework[9]inwhichalinearizedformoftheBTEissolvedinaniterativemanner.AnharmonicforceconstantsareusedtocomputescatteringratesandMean-Free-Path(MFP)parametersforeachmodewhicharethenusedtocom-putethethermalconductivitytensorofthebulkcrystal.Forthiscalculation,wehaveemployeda10x10x10k-10080604020Energy (meV)PDOSFIG.2.Phonondispersionofβ-Ga2O3usingareal-spacefinite-displacementapproach.PerturbativecalculationsareusedtoresolveLO/TOsplittings.Thesymmetrypathusedcanbefoundin[14]aswellasthesupplementarymaterial.ItisimportanttonotewehaveintroducedpointsA1,A2,andA3tocorrespondtothecrystallographicdirectionsoftheconventionalunitcell.spacesamplingmeshovertemperaturesrangingfrom25Kto1050K.Aftercalculationsareperformed,asym-metrizationofthecartesianκtensorisperformedbasedonthepointgroupsymmetryofthecrystaltoassureconsistencywhenrotatingthetensorintothecrystallo-graphicdirectionsoftheconventionalcell.Theformofthetensorhasbeengivenin[21]:κij=κ110κ130κ220κ310κ33(1)wherebysymmetryκ13=κ31.ShengBTEalsoallowsforsolutionstobeobtainedintherelaxationtimeapproximation(RTA).Thissignifi-cantlyreducesthecomputationalrequirementshoweverisnotguaranteedtoresultinanaccuratesolution.Westriclyadheretotheiterativesolutioninthisworktomaximizeaccuracy.Additionally,itmustbenotedthatShengBTEac-countsforboundaryscatteringbymeansofthecommonapproximationgivenin[8]aswellastheapproximationofscatteringduetonaturallyoccuringisotopes[7].Inthetheconventionalunitcell,thea([100])andb([010])directioncoincidewiththexandydirectionsre-spectively.Thevalueofκisfoundat300Ktobeκ[100]=16.06WK−1m−1,κ[010]=21.54WK−1m−1andκ[001]=21.15WK−1m−1.Whichisingoodagreementwiththeexperimentalvaluesofκ[100]=13WK−1m−1andκ[010]=21WK−1m−1[1,3–5].Generally,thetheoreti-calcalculationisexpectedtoovershoottheexperimentalvaluesduetovariousminoritysourcesofscatteringfoundinexperimentalsamplesofthecrystal.Thefulltemper-aturedependentκisshownforthethreecrystallographicdirectionsinFig.3alongwithfitsforlowandhightem-peraturebehavior.Eachregionisfittedtoκ=AT−m[8]whereAandmsources of scattering found in experimental samples, e.g., boundary effects at low tempera-
tures. The temperature dependence of κ is plotted in Figure 3.4 along with fits for low and
high temperature behavior.

Each range is fitted to κ = AT−m where A and m are fit parameters shown in the
Figure 3.4b inset. For high temperatures, the system is expected to be Umklapp dominated
and to exhibit a T−1 dependence [36], and that is the trend of the results presented here.
There is good agreement with the experimentally observed trends, especially the degree

of anisotropy in the [010] and [100] directions. At low temperature (approaching 25 K) those

results were not in agreement with the experimentally observed value for the [100] and [001]

directions from Guo et al. This is not surprising given the importance of boundary scattering

in that range and the lack of explicitly considered boundary effects considered in this BTE
description. At very low T, optical modes contribute less than 1 % to κ and acoustic modes
near Γ dominate transport. To capture the relevant behavior near the center, very dense
k-meshes are required. Unfortunately, very dense k-meshes are computationally unfeasible

due to memory scaling issues.[37]

3.3.3 Scattering Rates

The scattering rates are presented in Figure 3.5 with a temperature dependence. Overall,

there is a significant dependence on temperature which is an expected trend. This is ex-

plained by analogy of the phonon gas model, in which a very strong temperature dependence

exists in the Umklapp processes. Generally, the phonon scattering rate is dependent on the

phonon density in the gas. Using a bose-einstein description of the phonon number density

as has been done throughout this work, there is an inherent temperature dependence in the

distribution function. This leads to a phonon gas density which will be low for temperatures

smaller than the Debye temperature in which phonon-phonon interactions become effectively

zero. This is seen in Figure 3.5a when compared to Figure 3.5d which show nearly an or-

der of magnitude difference in scattering rate. As the temperature increases, the Umklapp

41

Figure 3.4: (a) Calculated and published experimental thermal conductivities vs. temper-
ature for various crystal orientations. (b) The fitting function κ=AT−m is plotted with a
dashed line (and parameters) for each κ tensor component.

42

3T ≤ 200KT > 200KA (x105) mA (x105) mκxx4.731.880.1991.13κyy8.811.930.2471.27κzz0.231.210.3501.26κxz7.222.260.1221.39(b)(a)FIG.3.(a)Calculatedandpublishedexperimentalthermalconductivitiesvs.temperatureforvariouscrystalorienta-tions.(b)Thefittingfunctionκ=AT−misplottedwithadashedline(andparameters)foreachκtensorcomponent.phononsabovethegap.ThiscanbeseeninFig.4,whichrepresentsthephasespacevolumesatisfyingenergyandmomentumconservationinthethreephononprocess[7].Thelackofagapinthe[010]directioninadditiontothelowergroupvelocityalong[010]then[001]wouldleadtolowervaluesofκifassumingacousticmodescarryalltheheatinthesystem.ReferringtoTable??,wecanseetheaveragegroupvelocityofopticalphononsisrelativelyhigh.Itislikelythatgiventhenumberofopticalmodesinβ-Ga2O3inadditiontothenon-negligiblegroupve-locities,theircontributiontoκmustbeexplored.Toinvestigatefurther,welookatthepercentcompositionofeachmodefromourBTEsolutioninTableI.Atroomtemperature,opticalphononmodesinthe[010]directionoverallcontributeover44%ofthetotalthermalconduc-tivity.Athightemperaturestheircontributionincreasestonearly50%.ThesecompositionsareconsistentthebehaviorpredictedfromthecomputedgroupvelocitiesgiveninTable??.302004006008001000Temperature (K)101102Thermal Conductivity (Wm−1K−1)[100] This Work[010] This Work[001] This Work[100] Guo et al. (Sn-doped)[010] Guo et al. (Sn-doped)[001] Guo et al. (Sn-doped)[010] Galazka et al.[100] Handwerg et al.FIG.3.Calculatedthermalconductivityfrom75Kto1050K.Thefittingfunctionκ=AT−misshownwithadashedline.Thefitisperformedforthelow-temperatureregionbelow200Kandthehightemperatureregionabove200K.TABLEI.FitparametersforthermalconductivityatlowandhightemperatureLowT(≤200K)HighT(>200K)A(x105)mA(x105)m[100]4.711.880.181.11[010]6.561.900.241.27[001]0.231.210.291.25arefitparametersshowninTableI.Forhightempera-tures,thesystemisexpectedtobeUmklappdominatedandistheoreticallyexpectedtofollowaT−1dependence[8].Thereisgoodagreementwiththeexperimentaltrendsincludingthedegreeofanisotropyinthe[010]and[100]directions.Toexplaintheoriginoftheanisotropyineachcrystallographicdirection,wethencomparecalculatedgroupvelocityaverages.Valuesweredeterminedas¯v=pv2LA+v2TA1+v2TA2foracousticmodes.CalculationswereperformedanalogouslyforopticalmodesandresultsaresummarizedinTableII.Atlowtemperature(25K)ourresultswasnotinagree-TABLEII.Averagegroupvelocityofacousticmodesalongcrystallographicdirections.κ(WK−1m−1)¯v(meV·˚A)300K1050KAcousticOptical[100]16.064.05606.973.4[010]21.545.73651.7212.1[001]21.155.42690.0120.70102030405060708090Energy (meV)0.40.60.81.01.21.4P3 (1x10−7)OpticalLATA1TA2FIG.4.Phasespacevolumecorrespondingtoallowedthree-phononprocesses.Abovethegapfrequency,thevolumequicklydiesoff.mentwiththeexperimentallyobservedvalueforthe[100]and[001]directionsfrom[3].Atsuchlowtemperaturesopticalmodescontributelessthan1%toκandacousticmodesstaynearthezonecenterduetotheirlowenergies.Tocapturetherelevantbehaviornearthecenter,verydensek-meshesarerequired.Unfortunately,verydensek-meshesarecomputationallyinfeasibleduetomemoryscalingissues[22].Additonally,theavailablemeasurementsonsuchlowtemperaturesislimitedandinconsistent[3,4].Thelowtemperatureregionisnotexplicitlyinterestingforthecommunityofinterestingforβ-Ga2O3asitsstandardusesareinelectronicdevicesoperatinginhightempera-turesregimes.Previousworks[3,4]onβ-Ga2O3thermalconductivityhaveassumedtheacousticmodeswerethedominantcon-tributionandopticalmodewerenegligibleduetotheirlowgroupvelocity.Howeverwehavefoundthatthisas-sumptioncannotexplainthe[010]direction’shighvalueofκ.FromFig.2,thereisnoopticalgapinthedisper-sionalongthe[010]directionwhileagapdoesexistinthe[100]and[001].Gapsareoftenattributedtoatomicmassdifferences[10],anditsexistencepreventslowfre-quencyacousticphononsfrominteractingwithopticalphononsabovethegap.ThiscanbeseeninFig.4,whichrepresentsthephasespacevolumesatisfyingen-ergyandmomentumconservationinthethreephononprocess.[7].Thelackofagapinthe[010]directioninadditiontothelowergroupvelocityalong[010]then[001]wouldleadtolowervaluesofκifassumingacousticmodescarryalltheheatinthesystem.ReferringtoTableII,wecanseetheaveragegroupvelocityofopticalphononsisrelativelyhigh.Itislikelythatgiventhenumberofopticalmodesinβ-Ga2O3inadditiontothenon-negligiblegroupve-FIG.4.Phasespacevolumecorrespondingtoallowedthree-phononprocesses.Abovethegapfrequency,thevolumequicklydiesoff.TABLEI.ModecompositionofThermalConductivity.Totalpercentcompositionofopticalmodesisshownfollowedbyadecompositionintolow(≤30meV)andhighenergy(>30meV)Contribution(%)toκat300KLATA1TA2Optical(Low/High)[100]21.041.720.516.8(12.3/4.5)[010]19.520.016.244.3(37.8/6.5)[001]31.022.821.025.2(18.5/6.7)Contribution(%)toκat1050K[100]20.540.419.919.2(12.7/6.5)[010]18.518.615.547.4(37.4/10.0)[001]29.821.920.128.2(19.2/9.0)Onesourceofthisstrongopticalcontributionmaybefromthe”mixing”oflowfrequencyopticalmodeswiththeLAmodeinFig.2.Thecontributionfromopticalmodesisparticularlyimportanttobeawareofwhenim-puritiesarepresent[9].Deformationstothechemicalbondinginthecrystalcansignificantlyaltertheabilityoftheopticalphononstocarryheat[23].ReferringtothePDOSinFig.2wecanseethatatopticalfrequencies,Oatomshaveasignificantlydominantcontribution.Stan-dardimpuritiesintroducedinβ-Ga2O3havebeenviaOvacancies[24,25],sowecanexpectsubstitutionorva-canciesofOsitestoimpactthemodecompositionandpossiblythermalconductivity.Thescaleofthiseffectisunknownforβ-Ga2O3andwillbeexploredinfuturework.Todeterminelengthscalesofrelevantcontributions,weanalzyethecumulativeκinFig.5foreachmodeasafunctionofamaximumMFPvalueweimposed.Thisisinformativefordeterminingwhatlengthscaleseachprocesses dominate and the normal processes become negligible.

(a) 100K

(b) 300K

(c) 500K

(d) 800K

Figure 3.5: Total anharmonic scattering rate for β − Ga2O3 at series of temperatures.

In the analysis of measured data, some have concluded that the acoustic modes domi-

nate thermal transport and that the optical modes contribute negligibly due to low group

velocities.[20, 21] It is found that this assumption cannot explain the comparatively large
value of κ[010]
propagating along [010] while gaps exist in the dispersions for [100] and [001]. Gaps are

. From Figure 3.3, there is no optical phonon gap in the dispersion for states

often attributed to atomic mass differences [28] and can suppress the interaction of acoustic

phonons with optical phonons above the gap. This can be seen in Figure 3.6, which plots each
mode and irreducible q-point’s dimensionless scattering phase space volume, P j
3 (q), satisfy-
ing energy and momentum conservation in the three phonon processes[38]. This parameter

43

Figure 3.6: Dimensionless phase space volume available for each mode (j) and irreducible q-
point (q) to participate in three-phonon scattering processes[38]. Above the acoustic-optical
phonon gap energy the available phase space is small, leading to reduced optical phonon
scattering.

can be expressed as

where,

and,

(±)
D
j

(−)

3 (cid:19)

(±)
j

(q)

P3 =

P

1
2

P

2

(±)

(+)
3 +

3Ω(cid:18)P
3 =(cid:88)j (cid:90) dqD
δ(cid:16)ωj(q) ± ωj(cid:48)(cid:0)q

(3.26)

(3.27)

(3.28)

(q) = (cid:88)j(cid:48),j(cid:48)(cid:48)(cid:90) dq

(cid:48)

(cid:48)(cid:1) − ωj(cid:48)(cid:48)(cid:0)q ± q

(cid:48) − G(cid:1)(cid:17)

In general, the smaller the phase space, the larger the phonon lifetimes and the higher

the thermal conductivity. One important feature seen here is the phonon energy gap in

the optical modes at around 65 meV. The lack of a gap in the [010] direction in addition
along [010] than [001] would lead to lower values of κ under the
to the lower TA mode vG

44

302004006008001000Temperature (K)101102Thermal Conductivity (Wm−1K−1)[100] This Work[010] This Work[001] This Work[100] Guo et al. (Sn-doped)[010] Guo et al. (Sn-doped)[001] Guo et al. (Sn-doped)[010] Galazka et al.[100] Handwerg et al.FIG.3.Calculatedthermalconductivityfrom75Kto1050K.Thefittingfunctionκ=AT−misshownwithadashedline.Thefitisperformedforthelow-temperatureregionbelow200Kandthehightemperatureregionabove200K.TABLEI.FitparametersforthermalconductivityatlowandhightemperatureLowT(≤200K)HighT(>200K)A(x105)mA(x105)m[100]4.711.880.181.11[010]6.561.900.241.27[001]0.231.210.291.25arefitparametersshowninTableI.Forhightempera-tures,thesystemisexpectedtobeUmklappdominatedandistheoreticallyexpectedtofollowaT−1dependence[8].Thereisgoodagreementwiththeexperimentaltrendsincludingthedegreeofanisotropyinthe[010]and[100]directions.Toexplaintheoriginoftheanisotropyineachcrystallographicdirection,wethencomparecalculatedgroupvelocityaverages.Valuesweredeterminedas¯v=pv2LA+v2TA1+v2TA2foracousticmodes.CalculationswereperformedanalogouslyforopticalmodesandresultsaresummarizedinTableII.Atlowtemperature(25K)ourresultswasnotinagree-TABLEII.Averagegroupvelocityofacousticmodesalongcrystallographicdirections.κ(WK−1m−1)¯v(meV·˚A)300K1050KAcousticOptical[100]16.064.05606.973.4[010]21.545.73651.7212.1[001]21.155.42690.0120.70102030405060708090Energy (meV)0.40.60.81.01.21.4P3 (1x10−7)OpticalLATA1TA2FIG.4.Phasespacevolumecorrespondingtoallowedthree-phononprocesses.Abovethegapfrequency,thevolumequicklydiesoff.mentwiththeexperimentallyobservedvalueforthe[100]and[001]directionsfrom[3].Atsuchlowtemperaturesopticalmodescontributelessthan1%toκandacousticmodesstaynearthezonecenterduetotheirlowenergies.Tocapturetherelevantbehaviornearthecenter,verydensek-meshesarerequired.Unfortunately,verydensek-meshesarecomputationallyinfeasibleduetomemoryscalingissues[22].Additonally,theavailablemeasurementsonsuchlowtemperaturesislimitedandinconsistent[3,4].Thelowtemperatureregionisnotexplicitlyinterestingforthecommunityofinterestingforβ-Ga2O3asitsstandardusesareinelectronicdevicesoperatinginhightempera-turesregimes.Previousworks[3,4]onβ-Ga2O3thermalconductivityhaveassumedtheacousticmodeswerethedominantcon-tributionandopticalmodewerenegligibleduetotheirlowgroupvelocity.Howeverwehavefoundthatthisas-sumptioncannotexplainthe[010]direction’shighvalueofκ.FromFig.2,thereisnoopticalgapinthedisper-sionalongthe[010]directionwhileagapdoesexistinthe[100]and[001].Gapsareoftenattributedtoatomicmassdifferences[10],anditsexistencepreventslowfre-quencyacousticphononsfrominteractingwithopticalphononsabovethegap.ThiscanbeseeninFig.4,whichrepresentsthephasespacevolumesatisfyingen-ergyandmomentumconservationinthethreephononprocess.[7].Thelackofagapinthe[010]directioninadditiontothelowergroupvelocityalong[010]then[001]wouldleadtolowervaluesofκifassumingacousticmodescarryalltheheatinthesystem.ReferringtoTableII,wecanseetheaveragegroupvelocityofopticalphononsisrelativelyhigh.Itislikelythatgiventhenumberofopticalmodesinβ-Ga2O3inadditiontothenon-negligiblegroupve-300 K
[010]
20
20
16
44

[001]
31
23
21
25

1050 K
[010]
19
19
15
47

[001]
30
22
20
28

[100]
21
40
20
19

[100]
21
42
20
17

LA (%)
TA1 (%)
TA2 (%)
Optical (%)

Table 3.2: Percentage of modal contributions to thermal conductivity for 300 and 1050 K.
There is significant anisotropy in these contributions with the [010] having a significantly
larger dependence on optical modes. The composition does not change significantly as tem-
perature is increased.

assumption that acoustic modes dominate. To investigate further, an examination of the
contributions to κ of each mode is performed using our BTE solution in Table 3.2.

At 300K, optical phonon modes contribute over 44% of the total κ010. At high temper-
atures their contribution increases to nearly 50%. These compositions are consistent with
the behavior predicted from the computed vG
higher κ without other contributions.

for acoustic modes, which would not support

Contributions from optical modes are not completely surprising and have been found in

many complex materials [28]. One source of the strong optical phonon contribution may

arise from the hybridization of low energy optical modes with the longitudinal acoustic (LA)

mode. This causes a flattening of the LA dispersion for states away from the zone center
and therefore a reduction in vG
this case, the contribution from optical modes may be strongly influenced by the existence

leading to reduced thermal transport from LA phonons. In

of impurities[28], where deformations to the chemical bonding in the crystal can significantly

alter the ability of the optical phonons to carry heat.[39] Referring to the PDOS in Figure 3.3

it is seen that at optical phonon energies, O atoms are the dominant contribution. Standard
point defects in Ga2O3 have been via O vacancies or substitutional impurities at O sites (like
Sn), so it can be expected that these impact the mode composition and thereby κ [40, 41].
This potentially opens up an opportunity for engineering of thermal conductivity through

impurities and defects.

To determine length scales of relevant phonon transport, the cumulative contributions are

45

parameterized using κ/κaccum = 1 + Λ0/Λmax where κaccum is the conductivity associated
with accumulated contributions up to a MFP cutoff of Λmax. The full value is κ and Λ0
can be interpreted as representative of the MFP of relevant heat-carrying phonons.[24] The
parameterized contributions are shown in Figure 3.7 for each mode as a function of Λmax.
This is informative for determining the length scales at which each phonon mode becomes

relevant for carrying heat.[28] Results for each mode are summarized in Table 3.3.

LA
187
62
380

[100]
[010]
[001]

Λ300K
0
TA1
320
63
174

(nm)

TA2
331
82
384

Optical LA
8
44
15
11
10
92

Λ1050K
0
TA1
79
16
45

(nm)

TA2
75
21
88

Optical
2
2
2

Table 3.3: Characteristic MFP, Λ0, indicating the length scale above which each phonon
mode contributes to κ.

There is a significant directional dependence on the MFP. This indicates an anisotropy

in the influence boundary effects will have as a function of length. This effect has been
experimentally observed for Ga2O3 thin films were used as gate contacts for high frequency
GaN devices. It was shown that κ spanned 1.5 orders of magnitude when an interface was
introduced [42].

3.3.4 Memory Scaling

There are several parameters that can be tuned in the workflow presented. It is important to

verify the set of parameters used is converged. Firstly, the sampling mesh for the DFT calcu-

lation in this work was chosen to be a 2x2x2 k-mesh in reciprocal space for the harmonic force

constants and gamma point only for anharmonic force constants. The mesh was confirmed

to converge for the 2x2x2 case by comparing with a gamma only mesh for harmonic force

constants. The force constants experience a change of at most .1% and thus are sufficiently

converged even with just a gamma point sampling. This is expected considering the large

46

Figure 3.7: Accumulated thermal conductivity in three crystallographic directions for each
mode as Λmax is increased. Fits to a parametric function are shown by the dashed lines.

size of the unit cell for β-Ga2O3 which grows considerably larger when using a supercell.
This large volume in real space shrinks the volume of reciprocal space considerably. When

reciprocal space is sufficiently small, gamma point only sampling will sufficiently sample the

brillouin zone for our purposes. The same principle applies to the anharmonic force constant

calculation. For the anharmonic calculation a comparison was performed analyzing the effect

of increasing the nearest neighbor interactions from 3 to 6 and have not found significant
changes to the final thermal conductivity (under 1 Wm−1K−1).

The algorithm for solving the BTE is not easily parallelizable and scales poorly with

respect to the size of the sampling mesh which will be discussed in more detail in a later

section. For this reason it becomes prohibitively expensive to use a dense sampling mesh. To

47

02468100246810Cumulative Thermal Conductivity (Wm−1K−1)10-1100101102103104Maximum Mean-Free-Path0246810LATA1TA2Optical[100] [010] [001] Figure 3.8: Convergence of Iterative and RTA methods with increasing mesh density

validate the convergence of the result with respect to the sampling mesh, an extrapolation

procedure can be used as explained in [24]. Results are computed for several mesh densities

(6x6x6, 7x7x7, 8x8x8, 9x9x9, 10x10x10, 12x12x12 and 14x14x14) and fitted using:

κnq = κ∞(1 − c
nq

)

(3.29)

where nq is the number of q-points sampled in each direction, κnq
is the corresponding
thermal conductivity, c is fitting parameter and κ∞ is the final thermal conductivity ap-
proximated to an infinitely dense mesh. Results for this extrapolation at 300K in both the

iterative and relaxation time approximation (RTA) are shown Figure 3.8

The final convergence parameter that must be addressed is the size of the supercell used

in both force constant calculations. Upon testing, the chosen supercell size is a 2x2x2 su-

percell based on published literature [43]. This cell consists of 80 atoms. Expanding this

cell to a 3x3x3 supercell is computationally very demanding due to the number of atoms

increasing to 270. Additionally, the size of the 3x3x3 cell will lead to configurations contain-

48

10152025Mesh Density nq20.521.021.522.022.5Average Thermal Conductivity (Wm−1K−1)IterativeRTAing displacements too small to be resolved without increasing the convergence threshold to
smaller than 10−5. Which becomes prohibitively expensive with such a large cell. Increas-
ing the displacement size can be a solution but can also potentially induce a anharmonic

content which will not be accounted for. For supercells larger than 2x2x2, it is likely that a

perturbation theory approach would be more computationally feasible and accurate.

The low symmetry of the monoclinic structure in addition to the large number of atoms

in the supercell introduces very large memory requirements. Additionally, the algorithm

requires a full set of phonon eigenvectors to be loaded into memory on each node for compu-

tation. So the number of processors cannot be very large for this structure. This particularly

input file required approximately 64 GB of memory per node. Additionally, the memory scal-

ing behaves as N5 where N is the number of k-mesh points in a single direction. So increasing

the mesh density further is unrealistic without access to very large memory nodes or running

on very few processors. The latter introduces issues with computational time. After testing

a sequence of mesh densities, It was found the 11x11x11 mesh to give the best solution

with reasonable resource requirements. The temperature was swept for a sequence of 20

temperatures, so it was important to conserve resources when possible. The majority of

the simulations performed in this work were performed on two 256GB memory nodes and

utilized 4 Xeon E5-2670 2.5GHz cores.

3.3.5 Choice of Functional

One topic that has not been explored in detail thus far is the effect of the exchange-correlation

term on the anharmonic force-constants and subsequently the thermal conductivity. This

effect is largely unstudied in even simple materials but its effect can potentially be important.

The effect of the LDA vs GGA functionals has proven to be significant in the determination of

mechanical properties such elastic constants due to the underestimation and overestimation

of the lattice constants respectively. Recent studies investigating the effect have found the

predictions of LDA and GGA to be nearly equivalent for thermal conductivity [44]. This

49

result is non-intuitive considering the LDA formulation over-binds while the GGA under-

binds and softens the IFCs. It was speculated in [44] that this phenomena arises from the

cancellation of errors in the respective approximation methods for the exchange-correlation

functionals. The LDA case generates larger third-order IFCs due to the overbinding but is

cancelled by the increased acousto-optical gap in the phonon dispersion, decreasing phase-

space for phonon-phonon scattering events and vice-versa for the GGA predictions.
3.4 α − Ga2O3

There does not currently exist any literature on thermal properties of single crystal α
phase due to the difficulty in maintaining the metastable phase for measurement. The
recent surge of interest in α − Ga2O3 as an alternative to β − Ga2O3 has increased the need
for a theoretical understanding analogous to this chapters treatment of β − Ga2O3 . The
complexity of the alpha phase is less than that of the beta phase due to the higher symmetry

the crystal possesses so the computational requirements are less restrictive.

3.4.1 Phonon Dispersion
The phonon dispersion for α − Ga2O3 is computed using a 3x3x3 supercell and 2x2x2 mesh
in reciprocal space to determine force-constants. The results along high-symmetry path

are presented in Figure 3.9 The non-analytical correction were again computed within a

DFPT framework in QuantumEspresso and give the born effective charges as well as the

high-frequency dielectric function summarized in Table 3.4 Overall the phonon dispersion
behaves similarly to the β phase, with comparable group velocities. The most distinct feature
is the lack of a sizable acousto-optical gap for any of the high-symmetry directions. The lack

of a gap implies smaller IFCs and softer bonds with an increased phase-space for anharmonic

scattering. The reduction in IFCs likely stems from the lack of an tetrahedrally coordinated

gallium cation present in the corundum structure. The lack of a gap is likely to cause a
reduction in thermal conductivity compared to the β phase.

50

Figure 3.9: Calculated phonon dispersion of α − Ga2O3 .


∞
Z∗ (Ga1)
Z∗ (Ga2)
Z∗ (O1)
Z∗ (O2)
Z∗ (O3)

xx
3.9
3.20
3.20
-2.48
-1.96
-1.96

yy
3.9
3.23
3.23
-1.88
-2.40
-2.18

zz
4.0
3.23
3.23
-2.08
-2.08
-2.30

Table 3.4: Born effective charges and high-frequency dielectric tensor determined from
density functional perturbation theory (DFPT) for α − Ga2O3

3.4.2 Scattering Rates

Anharmonic force-constants were computed and used in determining the phonon-phonon

scattering rates which are presented in Figure 3.10. The lack of a pronounced gap is evi-

denced further by these plots, as well as the increased overall scattering, with the magnitudes
increasing as much as 60% over the β phase. As temperature increases, the magnitude grows
more rapidly than the β phase as well.

51

10080604020Energy (meV)(a) 100K

(b) 300K

(c) 500K

(d) 800K

Figure 3.10: Total anharmonic scattering rate for α − Ga2O3 at series of temperatures.

3.4.3 Thermal Conductivity

Solving the BTE in the RTA framework gives the thermal conductivity tensor presented
in Figure 3.11. Unlike the β phase, the thermal conductivity is relatively isotropic with
a negligible anisotropy along crystallographic directions. Most importantly however, the
thermal conductivity is approximately half of the magnitude of the β phase. This is not
surprising given the larger scattering rates and lower phonon lifetimes.

The α − Ga2O3 polytype is most often grown and used in thin-film form, and it is
important to characterize what thermal conductivity suppression may arise from the size of

the film thickness. This also can be done within the first-principles framework presented
here following the methodology of [45]. Firstly, the in-plane (κ(cid:107)) and cross-plane (κ⊥) is

52

Figure 3.11: Calculated thermal conductivity of α − Ga2O3 . Unlike the monoclinic coun-
terpart, the axes are nearly isotropic

expressed within the RTA as an effective thermal conductivity of the form, eb

κeff (L) =(cid:88)k

Sk(L)Ck (cid:107)vk(cid:107) Λk cos2 ϑk

(3.30)

where the film thickness is L, and a suppression function S(L) is constructed to incorporate
the additional phonon scattering that a film boundary will provide. Additionally, ϑ represents
the angle between group-velocity vector and the direction of transport and Λk
is the MFP.
The in-plane suppression function is formulated as,

1 − p exp(cid:16)− 1

S(cid:107) =

K(cid:17) − (1 − p)K(cid:104)1 − exp(cid:16)− 1
K(cid:17)(cid:105)
1 − p exp(cid:16)− 1
K(cid:17)

(3.31)

where the parameter p indicates the specularity of the film, ranging from 0 to 1. For p = 0,
the film boundaries are perfect absorbers and for p = 1, perfect reflectors. The Knudsen
parameter is given by K and is computed within a transformed coordinate system [45] as,

K =

(ˆvˆz/(cid:107)v(cid:107)) · Λ

L

53

=

ˆΛˆz
L

(3.32)

Figure 3.12: Effective thin-film thermal conductivities for in-plane and cross-plane transport
with respect to the c-axis of α − Ga2O3 . Specularity was set to p=0.2

In the cross-plane direction, the transport vector u is along the normal direction n and hence
the film boundaries act as perfect absorbers and no specularity appears in the suppression

function expression,

with a Knudsen number that is simpler to calculate and can be written as,

S⊥ =

1

1 + 2K

K = Λ| cos ϑ|/L

(3.33)

(3.34)

Using the IFCs and subsequent scattering rates from the bulk thermal conductivity calcula-

tion, the effect of the thin-films are assessed and presented in Figure 3.12. It is seen that the

effect of the thin-film scattering is relatively strong and the thermal conductivity is nearly

half of its bulk value for 100 nm film thickness for cross-plane thermal conductivity. The

in-plane reduction is slightly less due to the nature of the suppression function formula-

tion. This reduction in thermal conductivity is at length-scales well in excess of what has

already begun to see engineering applications. For example, the recent success of a metal-

54

 0 2 4 6 8 10 12 100 nm 1 m 10 m 100 mThermal Conductivity (W/mK)Film ThicknessFigure 3.13: Isotopic scattering rates using natural isotopic distributions for β − Ga2O3

semiconductor-metal photodetector in [46] utilizes an α− Ga2O3 thin layer thickness of only
10 nm, an order of magnitude lower than where the thermal conductivity was found to be

halved.

3.4.4

Isotope Effect

The effect of the naturally occurring isotope scattering for various Gallium based semicon-

ductors has recently been investigated in [29] where it was determined that in the case of

GaN for example, the thermal conductivity can be increased as much as 65% through iso-

topic enrichment from a naturally occurring isotope concentration to pure isotope. This

isotopic scattering mechanism plays the role of a mass-disorder scattering term similar to

the effect an impurity would have. The effect of isotope concentration has been analyzed and
is presented here. The isotopic scattering rates are also presented in Figure 3.13 for the β
phase, And for the α phase in Figure 3.14 to be notably different than what was found in the

55

Figure 3.14: Isotopic scattering rates using natural isotopic distributions for α − Ga2O3

β phase with lower magnitude of scattering rates for the whole range of frequencies. Overall,
the effect the isotopic scattering has on the thermal conductivity can be seen in Figure 3.15

and appears to be a minor effect. The disparity between pure and natural isotope mixtures

is more prevalent at lower temperatures, where the effect of impurities is in general more
relevant. At most there is a 20% increase in thermal conductivity in the β phase with the
notably less in the α phase which is consistent with the isotopic scattering rates presented.
Considering the difference between the two crystal structures, it can be inferred that the

effect of isotope variation throughout the structure has more impact on the tetrahedrally
coordinated Gallium cations found in β − Ga2O3 which are not present in α − Ga2O3 .

56

Figure 3.15: Thermal conductivity of α − Ga2O3 and β − Ga2O3 using both nautral and
pure isotope mixtures.

3.5 Layered Structures

Use of a high thermal conductivity heat-sink for thermal management has been specu-
lated to be the only practical solution of remedying self-heating effects in Ga2O3 devices.
One could assume using as highly conductive a substrate/heat-sink as possible would yield

the greatest results but this has not been the case in practice [1].

In the case of nitride

semiconductors, CVD diamond has been used for heat removal and it was found that the

quality of the interface between the materials was limiting. There has been some preliminary
success with Ga2O3 self-heating reduction experimentally by using higher thermal conduc-
tivity heatsink materials, however it is difficult to characterize and predict what materials

will prove to be effective at heat removal due to the difficulty in achieving a strong adhesion

of metal contacts to the semiconducting material. Considering this limitation, it becomes

important to produce a more complete phonon model that incorporates boundary/interface

scattering effects in addition to simply phonon-phonon scattering. Realistic electronic de-

57

vices are typically made from layering of various semiconductor materials as well as metallic

ohmic contacts and this leads to a complicated temperature distribution which is controlled

primarily by how the energy carrying particles, whether it be phonons or electrons, behave

at an interface. Traditionally the multilayer problem could be simulated for basic material

layers using a very simplified Fourier model with some approximate model for the effect of an

interface. This however does not consider the atomistic properties of the underlying mate-

rial and assumes the layers to be isolated from one another in terms of quantum mechanical

considerations. The more accurate method to represent a realistic structure would be using

a molecular dynamics style simulation with a full description of atomic interactions. This is

often unrealistic however due to the very large number of atoms which would be required and

the corresponding computational expense. The solution then becomes to find an accurate

way to utilize the BTE approach outlined previously while still considering the interface.

As a remedy of the curse of dimensionality, there has been recent success at combining

the linearized BTE with a deviational monte-carlo solver that performs calculations orders

of magnitude faster than the traditional MC counterpart.

In this case, the equilibrium

component of the system is determined using the established linearized phonon BTE and the

non-equilibrium components split into particle-like wave packets with stochastic trajectories

that can be determined statistically, with MC techniques [45].

3.5.1 Variance-Reduced Monte-Carlo

Firstly, it is important to understand why traditional MC techniques are inefficient for com-

putations of properties such as thermal conductance by introducing a simple toy problem.

Suppose there is a spherical domain in which the center is fixed at a temperature of 301

K and the rest of the domain is fixed at 300 K. As with many particle MC techniques, the

problem can be discretized and phonons can be lumped into "superparticles" which are com-

putational constructs to represent a larger population of phonons with similar properties.

These particles are assigned an occupation given by the Bose-Einstein distribution and the

58

local temperature. The heat flux of this system is then determined by allowing particles

to move and scatter until they have reached a steady-state condition and the subsequent

heat-flux can be determined by summing the net particle number gain/loss for each "cell"

of the system. This is a very inefficient method of simulation due to how close the particles

already are to the equilibrium distribution function with respect to the very small tempera-

ture gradient they experience. This then requires a very fine resolution for simulation with a

very large number of computational particles to simulate. The problem becomes significantly

more complex when considering extrinsic events such as an electron-phonon scattering event

in which the resolutions required for the electron scale differ significantly from the phonon

scale.

Recently, a solution to one-dimensional thermal transport of multi-layered structures

with isothermal boundary conditions was introduced in [47]. The methodology implements

a variance-reduced Monte-Carlo scheme to solve the boltzmann equation for phonons. The

main idea behind these methods is to solve the phonon BTE using an approximation of the
nonequilibrium component, gk
approximation in situations such as the previously outlined example, where the phonon

, of the distribution of individual modes. This is a feasible

distribution functions are very nearly the equilibrium Bose-Einstein condition, which is of

course trivially obtained analytically. This concept is illustrated in Figure 3.16

This approximate gk

consists of a sum of a number of delta functions (energy packets)

all with the same energy and each with a well defined position. These packets are then

tracked using a stochastic monte-carlo method and their trajectories can be analyzed. Rather

than being actual phonons or super-particles of phonons, these packets are a mathematical

construct and labelled as “deviational particles”. They are simply abstract representations of

corrections to the Bose-Einstein distribution function which can be represented as positive or

negative depending on whether they are emitted from the hot or cold thermal reservoir. It is

also important to note that the solutions of trajectories of deviational particles is performed

in the linear approximation, where the temperature deviations are small. This assumption

59

Figure 3.16: Overview of the distribution function of a computational particle and the
simplification of the deviational particle approach. The blue shaded area represents the
actual deviation from equilibrium which is simulated in VRMC. Only a very small deviation
from the Bose-Einstein is common in problem where the thermal gradient is small.

is necessary to prevent a location dependence of phonon properties within the respective

materials and the deviational particles are independent of each other.

The lifetime of the deviational particles is illustrated, including the possible scattering

events that a particle experiences, in Figure 3.17.

This amounts to the following workflow,

1. Launch Event - This is the initial source function essentially, where we impose a

deviational intensity function from the reservoirs given by

Iiso =

(vk · n) H (vk · n) ωk [fBE (ωk, Tiso) − fBE (ωk, Tref )]

¯h

V (cid:88)k

Particles are randomly assigned to reservoirs with a corresponding probability

phot = |Ihot| / (|Ihot| + |Icold|)

pcold = 1 − phot

60

Figure 3.17: Overview of lifetime of deviational particle in the variance-reduced MC scheme.
Orange arrows represent the initial deviations imposed by a source intensity function. Green
shows the interface scattering mechanisms, pink the reservoir absorption and blue an intrinsic
scattering event where no interface effects are necessary to include.

2. Advection steps - Advancing the particle in time is the primary step of the MC
process. Using the group velocities of the particle vx, a random time step, χ, is
generated from an exponential distribution and the particle advances its position ∆x =
χvx, from here there are several possible final states depending on what the particle
experiences,

a) Layer/interface - the particle is reflected/transmitted according to a diffuse mis-

match model and the particle is assigned a new phonon mode while conserving
energy. For two materials M1 and M2 the incident particle is assigned a new
phonon mode from the distribution[45],

where G is the gaussian regularization of the delta function and the total number

k · n(cid:17) GM (ωk) /(cid:16)V M N M
tot(cid:17)

(3.35)

=(cid:16)vM

k · n(cid:17) H(cid:16)vM

61

of modes is given by

Ntot = NA × NB × NC × Nbranches

(3.36)

b) Boundary - when the particle reaches a boundary, it has completed the trajectory

and is considered terminated. The iterative process moves on to the launch of the

next particle and stores the heat flux contributions for post-processing

c) Pure advection - if the particle does not experience any boundary or interfaces

along the trajectory, it is assigned a new phonon mode solely based on the intrinsic

phonon lifetimes (determined from the phonon-phonon scattering rates)

3. Summing heat flux - To analyze the trajectories, the system is partitioned into bins.
At each bin, the total number of deviational particles which successfully crossed can be

summed to determine the total heat flux. From this, an effective thermal conductivity

for the entire structure can be expressed:

κeff =

qnet

(Th − Tc)/Ltot

(3.37)

Traditionally, the effectiveness of a thermal contact is more readily expressed by the

thermal conductance.

reff =

Thot − Tcold

qnet

=

Ltot
κeff

(3.38)

This parameter characterizes the amount of heat which crosses an interface in a given tem-

perature gradient.

In other terms, it can be described as the ratio of the temperature

discontinuity at an interface to the heat flux flowing across the interface. Often it is solely

this parameter which determines the effectiveness of a thermal coating and the ability to

tailor the TBC is crucial for thermal management engineering. On the atomic scale, TBC is

determined by the mismatch of phonon density of states of the two materials of the interface.

However in practice, the TBC can be significantly dependent on properties which describe

the quality of the interface mechanically. For example, the specularity, interfacial roughness,

62

interface chemistry and chemical bonding all can play a significant role potentially. In the

work performed here, an idealized model is used in which none of these chemical and me-

chanical factors are considered explicitly. With respect to these factors, any idealized model

adopted here would theoretically serve as an upper-bound estimate, in practice however this

is very dependent on the type of mismatch model which is adopted to represent the effect of

the interface.

3.5.2 Diffuse Mismatch Model

Traditionally, the effect of an interface can be modelled in a couple of different ways. The

method described previously, which is utilized in the VRMC, is essentially a form of diffuse

mistmatch model (DMM). In this method, the scattering at the interface is assumed to be

completely diffuse and elastic. The transmission probabilities are determined solely from the

phonon dispersion available on each side of the interface and the phonons have no "mem-

ory" of the mode it represented before scattering with the interface. The alternative model

commonly adopted for representing interfaces is the acoustic mismatch model (AMM) which

treats phonons in their planewave representation in a continuum model. Then, an acoustic

impedance can be assigned for each side of the interface. This method has shown success for

very low temperatures but is not an effective solution for the temperatures of realistic device

operation or for instances where optical phonons are also contributing to thermal transport.
As was previously shown throughout this chapter, the optical phonons in β − Ga2O3 do
contribute a non-negligible amount of group velocity to allow for heat to be carried by them

and thus cannot be easily ignored. One important restriction of the DMM is the mentioned

requirement of elastic scattering mechanisms only. The method does not consider inelastic

sources of scattering such as an electron-phonon scattering event which are in general difficult

to characterize [48]. In this work, the AlmaBTE [45] package is utilized which is restricted to

the assumption of fully elastic scattering such the an incident phonon can feasibly undergo

elastic mode conversion in a timely manner. In other words, if the discrepancy in phonon

63

Figure 3.18: Ohmic contact of Au/Ti on Si-implanted Ga2O3

density of states between the two materials of the interface are significantly different, it is

unlikely to experience primarily elastic phonon scattering.

3.5.3 Source-Drain Thermal Conductance
Using the previously outlined methodology for the case of β − Ga2O3 it is now possible to
estimate the thermal conductance for a real thermal contact in a transistor between source

and drain. The most basic interface in this scenario that has been experimentally considered
[42] has been a direct Au/β − Ga2O3 contact achieved by e-beam evaporated and sputtered
Au layers. Knowing the low thermal conductivity of β − Ga2O3 , this direct ohmic contact
will make a very poor thermal conductor. To improve this ohmic contact, an adhesive layer

is generally added as shown in Figure 3.18 with titanium as an example. In general, reactive

metals such as chromium and titanium can make for good adhesion layers.

In this setup, the layer thickness is an important variable to control the magnitude of

the thermal conductance. As the layer increases in thickness, it can potentially add more

resistance by increasing the amount of intrinsic phonon scattering which occurs before a

phonon reaches an interface. In this work, the VRMC algorithm is used with first-principles
determined force-constants and scattering rates for β − Ga2O3 with the contact materials
chromium and titanium to assess how well β − Ga2O3 can form thermal contacts and if
the thermal boundary conductance (TBC) is a limiting factor in thermal management. The

64

reason for these materials is two-fold, firstly they are known to make good adhesive materials

for metal-dielectric interfaces. Titanium in particular is perhaps the most common metal-

dielectric adhesive layer material. It has recently been shown however that chromium can
significantly improve the TBC of Ga2O3 interfaces due to the overlap it has in the phonon
density of states.
In general, the metal-dielectric interface will be dominated by elastic

phonon scattering incidents, as discussed previously. In these events, phonons of the same

frequency will couple at the interface, hence the higher overlap between phonon density of

states functions, the more phonon transmission will occur and the higher the TBC will be.

In the debye model the phonon DOS functions can be roughly sketched as Figure 3.19. The

Figure 3.19: Debye phonon density of states comparison for β − Ga2O3 and the interface
materials. Chromium and Titanium both contain a significantly improved overlap in the
density of states, leading to increased elastic scattering and phonon transmission at the
interface.

VRMC calculation is run for a sweep of layer thickness for the chromium and titanium, with
the β − Ga2O3 thickness being held at 400nm. This setup was analyzed experimentally
in [42] where a structure similar to that shown in Figure 3.20 was used. The number of
computational particles used is in the VRMC solver N = 10000 with the thermal boundaries

65

Figure 3.20: Thermal contact for Au/β − Ga2O3 using a wedge adhesive. This illustrates
the experimental setup of [42] this work compares to. In this work, successive simulations
are run for the range of layer thickness presented in the wedge of adhesive Ti/Cr.

being held at 300 K for the cold side and 300 K for the warm side. The scattering rates
for β − Ga2O3 were computed from a 12x12x12 grid in reciprocal space. These simulations
represent a perfect boundary which in practice is unrealistic. Surface roughness and surface

chemistry effects can often change the TBC significantly. An initial calculation for the
TBC using a direct Au/Ga2O3 interface is computed with a resulting phonon TBC value
of Gp = 61 MW/m2K which is comparable to the experimentally determined value of of
exp
p = 80 MW/m2K [42]. This value is generally low for a metal-dielectric interface, which
G
is expected considering the small density of states overlap between these two materials.

Maximizing this value is crucial to device engineering as the heat dissipation controls the

stability of integrated circuits. To further investigate, the Kapitza length can be defined by

Lk =

κavg
Gp

(3.39)

This length determines the required thickness of the material with average thermal conduc-
tivity κavg to have an equivalent resistance with the interface. The Kapitza length For the
direct Au/Ga2O3 interface is found to be Lk = 344 nm, which can be restrictive for device
engineering depending on the applications.

66

Au (70nm Layer)Ga2O3 (400nm Substrate)Ti/Cr Adhesion Layer (1-14 nm) Figure 3.21: VRMC results for Ti/Ga2O3 interface with a Ga2O3 layer thickness of 400 nm
and an increasing layer thickness for titanium. Results are compared to FDTR measurements
of [42]. VRMC results under predict measured data moderately and follow a similar trend
in layer thickness, with an increasing conductance until around 4 nm Ti.

3.5.4 Au/Ti/Ga2O3 Thermal Contact

To evaluate the effect of titanium as an adhesion layer, the VRMC method is performed
sequentially for the layer thickness in increments of ∆L = 0.5 nm. Results are presented
in Figure 3.21 and compared to the frequency-domain thermoreflectance (FDTR) results of

[42] Overall, there is a relatively good agreement with the measured trends. The magnitude

of the TBC differs in the VRMC prediction by about 60 MW/m2K, or approximately 30%.
The predicted and measured values reach their peak values of Gp = 190 MW/m2K and
Gp = 254 MW/m2K respectively, at around 4 nm in both cases. For the VRMC results, this
corresponds to a Kapitza length of Lk = 110 nm which is over a factor of 3 improvement
to the direct Au/Ga2O3 contact. This is a very significant improvement to the thermal

67

Adhesion Layer Thickness (nm)TBC (MW m-2K-1)100200300400500600123456789Measured Ti -Ga2O3 (Aller et al.)VRMC Ti -  Ga2O3boundary conductance for engineering applications. Discrepancies between measured and

predicted results are most likely due to the lack of inelastic scattering mechanisms considered

at the boundary. It is also possible that the simulation would require a more dense reciprocal

space mesh, beyond 12x12x12, for computing the anharmonic scattering rates however this

is computational prohibitive.

3.5.5 Au/Cr/Ga2O3 Thermal Contact

For the chromium layer example, it is expected that the TBC will improve if the density

of states model is sufficiently predictive and indeed the primary method of phonon trans-

mission at the interface is due to elastic phonon-phonon scattering. The results for the

chromium case, using the same simulation cell as the titanium example, are presented in
Figure 3.22. Compared to both the Au/Ga2O3 contact and the Au/Ti/Ga2O3 contact,
the effect chromium has on the TBC is significant with the plateau value at 5 nm being
Gp = 395 MW/m2K. This compares well with the measured value of Gp = 454 MW/m2K
at the same width, differing only by 13%. This corresponds to over a factor of 6 improve-
ment to the Au/Ga2O3 contact and a corresponding Kapitza length of only 53 nm. This
is a very large value for boundary conductance which significantly exceeds even many GaN
interfaces [50]. This implies that in thermal engineering of Ga2O3 devices, the interface will
not prove to be a significant obstacle. The trend seen in Figure 3.22 is consistent with the
Al2O3 measurements but differ slightly with the Ga2O3 measurements. Notably, the peak
at around 2 nm in the measured results is unusual. It was speculated in [42] that due to the
high heat of formation of Ga2O3 , it is likely that a thin amorphous layer is forming at the
interface. This phenomena is not seen in the Al2O3 measurements due to the lower heat of
formation found in Al2O3 . Surface chemistry phenomena such as this will not be considered
in simulation, and hence the peak is not present in the computed result.

68

Figure 3.22: VRMC results for Cr/Ga2O3 interface with a Ga2O3 layer thickness of 400 nm
and an increasing layer thickness for chromium. Results are compared to FDTR measure-
ments of [42] as well the Cr/Al2O3 measurements for the same structure [49]. VRMC results
under predict measured data moderately and reach a similar plateau value at 5 nm. The peak
in the measured Cr/Ga2O3 results is predicted to be from amorphous layer formation[42]

3.5.6 Temperature Dependence

To investigate the temperature dependence of the TBC, the VRMC computation is fixed

at the 4 nm plateau value for titanium and 5 nm plateau value for chromium and repeated

for a series of temperatures between RT and 600 K. From the bulk thermal conductivity
results in β − Ga2O3 , it is not expected that there will be a sizable temperature dependence
in the anharmonic scattering rates and hence the TBC will not be strongly temperature

dependent. The results for both Ti and Cr interfaces as well as measured results is presented

in Figure 3.23. Overall, the temperature dependence is small when considering only up to

third-order phonon-phonon scattering mechanism. In general, it is possible that at higher

temperatures (above 700 K), the effect of fourth-order phonon scattering may become more

69

Adhesion Layer Thickness (nm)TBC (MW m-2K-1)100200300400500600123456789Measured Cr -Ga2O3 (Aller et al.)Measured Cr - Al2O3 (Jeong et al.)VRMC Cr -  Ga2O3Figure 3.23: Thermal boundary conductance for fixed layer spacing as a function of temper-
ature. GaN Measurements taken from [50] and β − Ga2O3 measurements from [42]

influential which has been reported to be the case for GaN [51].

3.5.7 Role of Electron-Phonon Interaction

Electronic contributions to the thermal conductivity were not considered so depending on

the importance of the electron-phonon scattering channel at the interface, the transmission

can be underrepresented compared to experiment. If the thermal conductance is expressed
as Gtot, the contributions sum in parallel,

Gtot = (Ge−pGp)/(Ge−p + Gp)

(3.40)

The role of electron-phonon scattering at the interface is not fully understood however.

Many authors have speculated that electron-phonon scattering at interface is negligible and

70

Measured Au-GaN (Donovan et al.)Measured Au/Ti-GaN (Donovan et al.)Au/Cr-Ga2O3Au/Ti-Ga2O3Au-Ga2O3Measured Au/Cr -Ga2O3 (Aller et al.)Measured Au/Ti -Ga2O3 (Aller et al.)the interaction only becomes relevant for transport solely within the metal, and only phonon-

phonon scattering is relevant at the boundary. Simultaneously however, comparisons between

theory and experiment such as the results of this work, demonstrate an underestimation of

TBC when electron-phonon interactions are neglected. Although the metals used in this

study are known to possess a low electron-phonon coupling constant, the large range of
phonon frequencies available in Ga2O3 lead to more numerous possible scattering channels
involving an electron-phonon scattering event. However, it is difficult to draw a conclusion

from this trend alone considering the complexity of the transport problem and the difficulty

in characterizing the exact surface chemistry of a metal-dielectric interface in an experimental

setup. As was mentioned, the possible formation of amorphous layers, can be difficult to

observe but can influence the TBC and this phenomena will not appear in a computational

model without explicit inclusion. Additionally, in a realistic setup, the surfaces will never

truly make a perfect contact, and there will certainly be some degree of roughness at the

interface. The surface roughness will contribute an additional scattering term in theory

which will further increase the thermal resistance. In general, the problem is very difficult to

solve without a fully coupled electron-phonon Monte-Carlo solver using first-principles inputs

which is a considerably complex computational problem that has not yet been solved.

71

CHAPTER 4

THERMODYNAMICS AND MECHANICAL PROPERTIES

An analysis of the elastic and thermodynamic properties is critical for a complete under-
standing of the mechanical stability of Ga2O3 . The elastic stiffness tensor characterizes the
linear response to applied stresses and accurate values allow for the analysis of structural

and lattice dynamical properties in an elastic continuum model. One common application

of elastic constants for electronic device design is to predict the acoustic deformation po-

tential. Recently, deformation potential estimation was used to predict the electron impact

ionization coefficients.[52] However, due to the lack of available data on the elastic prop-
erties of Ga2O3 the elastic properties of gallium nitride were used as a substitute. That
approximation was justified under the assumption that under high electric field conditions

the non-polar optical phonon contribution was negligible.[52] However, it has recently been

shown that this non-polar optical phonon scattering is the dominant contribution to electron

mobility,[53] and a full understanding of elastic properties is needed for such an analysis.

Unfortunately, measuring these quantities is quite challenging due to the variability

of growth quality as well as the anisotropy of monoclinic crystals.[12, 54] Only recently,

elastic stiffness constants were measured at room temperature[55] and also theoretically
predicted.[14, 56] Since the proposed applications for Ga2O3 lead to a wide range of oper-
ating temperatures and the scarcity of measured data, the temperature dependencies of the

elastic properties are also considered in the present investigation.

Advances in first-principles DFT calculations have made robust computations for elastic

constants in any crystal symmetry possible with minimal computational overhead.[57] In this

work calculations were performed for the linear thermal expansion coefficient and elastic stiff-
ness constants (Cij) which are used to evaluate the complete set of elastic moduli including
the bulk modulus, shear modulus Young’s modulus, Vicker’s hardness, and Poisson’s ratio.

Additionally, temperature effects are considered by use of the quasi-static approximation[58]

72

and further evaluation of the isothermal and isentropic Cij.

More recently, Ga2O3 has become a material of great interest for its various applications
it has been used in solar-blind UV photodetector applications[59]

in optoelectronics, e.g.

due to the wide optical bandgap. Characterizing the optical properties in practice is a

challenge for low-symmetry crystals as they exhibit biaxial optical responses and possess

significant birefringence. There have been several recent investigations detailing this intricate

geometry dependence [60] as well as the incident wavelength dependence in the context

of absorption spectra. The effects of mechanical strain have not been considered in the

literature thus far but have been theorized to potentially play a significant role in optical

characteristics [60] and are generally an area of great interest for the development of crystal

growth methodologies[55].

The response of optical properties to strain is commonly known as the elasto-optic effect,

or photoelastic effect, and is characterized by strain dependence of the dielectric function.

Measurements of photoelasticity are notoriously difficult for even the simplest of crystals

due to the many sets of strains that must be individually applied to resolve the tensor. For

low-symmetry crystals with significant anisotropy it is challenging to control strain at the

precision required for resolving measurements and it is challenging to distinguish the response

in individual crystallographic directions. Optical interferometry can possibly be used, but

requires a large number of measurements on precisely cut crystals for all orientations which
is unrealistic for Ga2O3 at present.
principles for the dielectric function spectrum under strained conditions to resolve the elasto-

In this work calculations were performed from first-

optic properties.

4.1 Methodology

In this work, a first principles based approach as implemented in the Quantum Espresso

package using the projector augmented wave method is utilized. Comparisons are made be-

tween the effect of local density approximation (LDA) and the modified generalized gradient

73

approximation of Perdew-Burke-Ernzerof (GGA-PBEsol) on the elastic stiffness constants.

A plane wave basis with a kinetic energy cut-off of 125 Ryd was used for all calculations.
Brillouin zone integration was performed with a 16 × 16 × 8 k-point grid for structural re-
laxation and 8 × 8 × 8 grid for calculating the stress in the elastic constants calculations.
Figure 3.2 shows the orientation of the monoclinic structure.

The method adopted here computes the stress tensor (σij) directly from first principles
after applying a strain (εij) and uses the stress-strain form of Hooke’s law, σij = Cijklεkl
to compute the elastic stiffness constants.[61, 62] This approach was chosen over the “stress

approach” [57] because it reduces the number of calculations required for a low symmetry

structure such as the monoclinic phase being studied in this work. The ElaStic[57] software

package is used to identify the required strained configurations for this crystal. Total en-

ergy calculations were then performed on these strained configurations and post-processed

using ElaStic. Phonon calculations for determining the free-energy in the quasi-harmonic

approximation (QHA) are done in the real-space finite-displacement supercell method as

implemented in Phonopy[7]. Using this methodology, force-constants were computed using a
2× 2× 2 supercell containing 80 atoms with a 4× 4× 3 k-point grid. The phonon density of
states was determined for seven volumes and used to compute the vibrational contribution

to the Helmholtz free energy.

In the Voigt notation there are 13 independent elastic constants[63] for this structure and

the elastic stiffness constants relate the stress and strain by



σ1

σ2

σ3

σ4

σ5

σ6



=



C11 C12 C13

0 C15

C21 C22 C23

C31 C32 C33

0 C25

0 C35

0

0

0 C44

0 C46

C51 C52 C53

0 C55

0

0

0

0 C46

0 C66

0

0

0





ε1

ε2

ε3

ε4

ε5

ε6



(4.1)

where σ1−3 represent normal stresses and σ4−6 shear stresses. The notation is consistent

74

with the axes definitions in Figure 3.2. This choice is not unique and some literature uses a

second viable notation based on a different labeling of lattice vectors.

To minimize the set of strains required to fully determine Cij, the universal linear-
independent coupling strains are used which results in 5 strains for the monoclinic crystal
for a given magnitude of Lagrangian strain. [64] The strain values chosen are η = ±0.007
and η = ±0.013 for all calculations as in the work by Shang et al..[61, 62] Finite temperature
effects were considered by using the QHA which has been discussed in detail[65].

Elastic moduli were computed from the Cij using the Voigt-Reuss-Hill method[66][67].
) was computed using the simplified phenomenological model[68]

Vicker’s Hardness (HV
which depends solely on the shear modulus. Sound velocities along different crystallographic
directions were evaluated from Cij and compared with those obtained directly from the
phonon dispersion of Ga2O3 . [69]

4.1.1 Elastic stiffness

The methodology adopted in this work requires several steps which are briefly described.

Firstly, the elastic stiffness constants at 0 K can be determined using one of two methods,

both of which require DFT calculations to be performed on a series of strained structures.

The first method[57] uses the strain energy density vs strain relationship. This method only

requires total energy to be computed for the strained structures. The method requires many

energy calculations to be performed and is computationally expensive for low-symmetry

structures. For this reason this work adopts the second method, which directly computes
the stress tensor (σij) from first principles after applying a strain (εij) and uses the stress-
strain form of Hooke’s law, σij = Cijklεkl
). This
reduces the number of calculations required for low symmetry structures[57]. Choosing a

, to compute elastic stiffness constants (Cijkl

minimal set of strains to apply is heavily dependent on crystal symmetry. Hence, in this work

the “universal linear-independent coupling strains” (ULICS) [64] were utilized. The ULICS

basis is a linearly independent set of strains that span the “strain space” and consequently

75

represent the minimal set of strains required to fully determine Cij. For monoclinic this
results in 5 strains to apply to the crystal for a given magnitude. Applying several values

of strain is required in practice to assure the result is converged. The values of Lagrangian
strain used range from η = ±0.007 and η = ±0.013 for all calculations performed here.

4.1.2 Quasi-harmonic approximation

To consider finite temperature properties using DFT, this work utilizes the quasi-harmonic

approximation. This approximation assumes the dominant anharmonic contribution to the

phonon disperion be contained within the volumetric expansion rather than the explicit

temperature dependence of anharmonic force constants. This allows the Helmholtz energy

to be written as:

F (V, T ) ≈ E(V ) + Fphonon(V, T )

(4.2)

Where E(V ) is the equation of state function and Fphonon(V, T ) is the free energy con-
tribution due to vibrations about the ground state. Both contributions will require DFT

calculations for several unit-cell volumes. The former is determined by fitting total-energy

of the various volume configurations to an equations of state (EOS), while the latter contri-

bution requires the phonon density of states to be computed at these volumes. Bulk modulus
(B0), volume-temperature relationship V (T ), and heat capacity at constant pressure (CP
) or
) can be determined from the computed F (V, T ). Once the free-energy
constant volume (CV
function is determined, minimizing results in the Gibbs free-energy function which can be

used within several thermodynamical relations for bulk thermodynamic quantities, such as

bulk modulus:

heat capacity at constant pressure,

∂V 2(cid:19)T
B = V (cid:18) ∂2F
CP = −T(cid:18) ∂2G
∂T 2(cid:19)P

76

(4.3)

(4.4)

and heat capacity at constant volume,

CV = −T(cid:18) ∂2F
∂T 2(cid:19)V

These quantities are used to compute the bulk thermal expansion coefficient αV

CP (T ) − CV (T ) = α2

V V (T )T B0(T )

:

(4.5)

(4.6)

QHA analysis is necessary for determining the volumetric thermal expansion relationship and

subsequent thermodynamic properties such as heat capacity. It is also capable of predicting

bulk modulus as a function of temperature which can be compared to a similar prediction

made by direct computation of elastic stiffness constants. With an accurate representa-

tion of the volumetric thermal expansion, it is then possible to estimate the temperature

dependence of mechanical properties such as elasticity withing the so-called quasi-static ap-

proximation. Additionally, the QHA enables the computation of a macroscopic Gruneisen

parameter without needing to consider microscopic behaviors of interatomic force-constants.

Hence, a direct comparison between the anharmonic force-constant approach presented in

the previous chapter can be made with the macroscopic thermodynamic properties presented

here. This type of calculation requires successive phonon calculations to be performed within

DFT and the reader is referred to the previous chapter for the phonon calculation details.

Generating the strained-unit cells required to resolve a volumetric dependence entails ap-

plying hydrostatic strain followed by geometric relaxation within DFT. It is important that

a hydrostatic strain is used to adjust the volume as this ensures a minimal change lattice

parameters while retaining the desired volume change. Significant changes to any one lattice

parameter in an anisotropic crystal structure can lead to improper geometry optimization

that can more readily result in an unintended, and possibly artificial, phase transformation.

The determined free energy curves from a series of 20 different volumes in which phonon

properties are computed is presented in Figure 4.1.

The corresponding bulk modulus as a function of temperature in the QHA is now pre-
sented in Figure 4.2. The 0 K value of B0 = 166.5 GPa compares relatively well to the

77

Figure 4.1: Free-energy curves computed within QHA for β − Ga2O3 and minimized to
determine the Gibbs free energy function (red line).

recently measured data of Ref.[70] which gives a value of B0 = 182.6 GPa. This value is
notably lower than θ − Al2O3 which possesses the same crystal structure and spacegroup,
where B0 = 202.7 GPa [70]. This is not particularly surprising considering the volume of
θ − Al2O3 is as much as 20% lower than its gallium counterpart. The combination of a
smaller volume and lighter atoms leaves the crystal with less overall compressibility and
stiffer bonds which is seen in the significantly larger debye temperature of θ − Al2O3 [61].

78

Figure 4.2: Bulk modulus computed within QHA for β − Ga2O3 . At 0 K the bulk modulus
value is B0 = 166.5 GPa

4.1.3 Microscopic vs Macroscopic Gruneisen

To characterize the anharmonicity of the crystal in both the QHA and the phonon-phonon

scattering frameworks, the Gruneisen parameters can be computed. Large Gruneisen pa-

rameters are associated with strong anharmonicity where atomic motion breaks periodicity

of normal modes. In general they are a material property which can be defined in several

different ways which measure how the phonon frequencies change with volume. They can be

computed microscopically using the 3rd order IFCs,
ej∗
ακej
(cid:112)MκM(cid:48)

j (q)(cid:88) Φ0κ;l(cid:48)κ(cid:48);l(cid:48)(cid:48)κ(cid:48)(cid:48)

γj(q) = − 1
6ω2

αβγ

βκ(cid:48)

κ

× rl(cid:48)(cid:48)κ(cid:48)(cid:48)γeiq·Rl(cid:48)

(4.7)

where in this context they act as a measure of how the restoring forces of harmonic motion

change as the crystal expands. For direct comparison with measured data, it is more easy

79

to express this in an averaged scalar form.

V

(4.8)

γ = (cid:80) γjCj
(cid:80) Cj

V

Additionally, the effect of a changing volume can be computed more directly in the QHA

framework,

with the averaged scalar form being determined by,

γj(q) = −∂ ln ωj(q)

∂ ln V

γ = B

αV
CV

(4.9)

(4.10)

Computing the Gruneisen parameters using both of these methods allow for a more direct

validation and comparison to be performed between the microscopic and macroscopic calcu-

lations. This should serve as a validation of convergence of the 3rd order IFCs as well as a

check that the volumetric dependence within the QHA is consistent. The direct comparison

of the averaged scalar quantities is presented in Figure 4.3. Overall the parameters were

relatively small which is expected in stiff crystals. More importantly, there is a very con-

sistent match between the two methods for a modest range of temperatures. The growing

disparity for increasing temperatures is likely due to the onset of higher order anharmonic

effects in the microscopic formulation arising from fourth-order phonon scattering. This ef-

fect although small can be important to consider for high temperatures. However in this

case, the temperatures where the higher-order scattering become relevant would approach

the phase transition temperatures. The spectrum averaged mode Gruneisen parameters is

provided in Figure 4.4 where a minimal temperature dependence is seen.

4.1.4 Linear Thermal Expansion

Finite temperatures introduce anharmonic effects due to thermal expansion as has been

demonstrated in the previous section. This is especially the case in a monoclinic unit-cell,

where a high degree of anisotropy in the thermal expansion is expected along crystallographic

80

Figure 4.3: Average value of Gruneisen parameter of β − Ga2O3 computed macroscopically
through the QHA method compared to the microscopic formulation utilizing 3rd order IFCs

directions. Determining the linear thermal expansion coefficients (TEC) can be done by the

quasi-static approximation[58]. In this approximation, the thermal expansion is determined

by the optimized lattice vectors for several strained configurations of the unit-cell at 0 K com-

bined with the volumetric thermal expansion determine with the QHA. The strain definition

of thermal expansion at constant pressure[63] (using Voigt notation),

∂T (cid:19)P
αi =(cid:18) ∂εi(T )

(4.11)

Depending on the symmetry of the crystal, there will be applied strains εj which are uniquely
parallel to an individual direction j. Along these directions, the linear thermal expansion
(εj) relation and their corresponding temperatures defined by the QHA V (T ). Using the
computed 0 K values for V (εj) with the inverse relationship T (V ) produces a form of the tem-
perature dependence T (V (εj)) that relates to a specific direction through its corresponding
strain.

81

Figure 4.4: Frequency-dependent average mode Gruneisen parameters in β − Ga2O3 for a
sequence of temperatures computed within QHA.

4.1.5

Isothermal elastic stiffness

Calculating the temperature dependence quasi-statically is done analogously to the linear
thermal expansion. Elastic stiffness constants, Cij(V ), is computed for several volumes and
combined the quasi-harmonic V (T ). Then, Cij(T ) can be computed by interpolation.

4.1.6

Isentropic elastic stiffness

Often it is the case that experimental measurements for Cij at high temperatures are not
isothermal but isentropic. For example, the commonly used resonance method[61] results in

higher elastic wave speed relative to the heat diffusion, which results in adiabatic (constant

entropy for a reversible system) behavior. Therefore, this work uses an expression to compute

82

the isentropic stiffness constants, CS
heat capacity CV

and the linear thermal expansion αj[71]:

ij(T ) in terms of the isothermal constants CT

ij(T ), the

ij = CT
CS

ij +

T V λiλj

CV

λi = − 6(cid:88)k=1

αkCT
ik

(4.12)

(4.13)

4.1.7 Elastic moduli

Using Cij and the corresponding elastic compliance Sij, the bulk modulus (B) and shear
modulus (G) can be computed the Voigt-Reuss-Hill method[66]. Voigt’s approximation con-

siders a uniform strain as opposed to a uniform stress as in Reuss’ approximation. The

approximations serve as a lower and upper bound respectively. Hence, the average (Hill ap-

proach) of these two approximation will be used. These approximations are shown below[67]

BVoigt =

GVoigt =

+

2
9

1
((C11 + C22 + C33) +
(C12 + C13 + C23))
9
((C11 + C22 + C33) − (C12 + C13 + C23))
1
15
1
5

(C44 + C55 + C66)

BReuss = [(S11 + S22 + S33) + 2(S12 + S13 + S23)]−1
GReuss = 15[4(S11 + S22 + S33) − 4(S12 + S13 + S23)

+ 3(S44 + S55 + S66)]−1
(BReuss + BVoigt)

BHill =

GHill =

(GReuss + GVoigt)

1
2
1
2

Poisson’s ratio (ν) and Young’s modulus (E) are then computed with

ν =

1

2(cid:34) B − 2
3 G(cid:35)

B + 1

3 G

83

(4.14)

E =

9BG

G + 3B

(4.15)

Lastly, Vicker’s Hardness (HV

) can be computed using the simplified phenomenological

model [68] which depends solely on the shear modulus,

This model has been demonstrated to be effective in covalent materials such as oxides[68].

G = 6.78 HV

(4.16)

4.1.8 Sound velocities

Thermodynamic quantities can often be accurately estimated from sound velocities and the
Debye temperature (ΘD
validity of the computed elastic constants. In the bulk case, there are simple relationships

). Additionally, they are easily measured and thus can support the

depending only on previous calculations for bulk and shear modulus

ΘD =

where the mean velocity is written as

h

k(cid:20) 3n
4π(cid:18)NAρ

vm

M (cid:19)(cid:21)1/3
t(cid:33)(cid:35)−1/3

1
v3

+

3(cid:32) 2
with the longitudinal velocity, vl =(cid:113)(B + 4

vm =(cid:34)1

v3
l

taining the sound velocity along a specified direction is more complicated and requires the

3)/ρ, and transverse velocity vt =(cid:112)G/ρ. Ob-

solution of the Christoffel equation [72]:

(Γij − ρv2δij)uj = 0

(4.17)

(4.18)

(4.19)

(4.20)

where lk
and lm represent the direction cosines, uj is the displacement vector, and v is
the sound velocity eigenvalue that is solved for by diagonalization. The direction cosines in

Γij =

Cikjmlklm

3(cid:88)k=1

3(cid:88)m=1

84

in our orientation are with respect to the x-axis, l1 = cos α, y-axis, l2 = cos β and z-axis,
l3 = cos γ. For a monoclinic crystal, the solution in an arbitrary direction will result in a two
mixed modes (labelled v− and v+ for the slower and faster respectively), and a transverse
mode (vt)[69]. Along the crystallographic directions, the solutions can be expressed as:

[100] :

vt =(cid:112)C66/ρ
v+ =(cid:115)(C11 + C55) + (C2
v− =(cid:115)(C11 + C55) − (C2

15 − 2C11C55 + C2
11 + 4C2
2ρ
15 − 2C11C55 + C2
11 + 4C2
2ρ

55)1/2

55)1/2

[010] :

vt =(cid:112)C22/ρ
v+ =(cid:115)(C44 + C66) + (C2
v− =(cid:115)(C44 + C66) − (C2

46 − 2C44C66 + C2
44 + 4C2
2ρ
46 − 2C44C66 + C2
44 + 4C2
2ρ

66)1/2

66)1/2

[001] :

vt =(cid:112)C44/ρ
v+ =(cid:115)(C33 + C55) + (C2
v− =(cid:115)(C33 + C55) − (C2

35 − 2C33C55 + C2
33 + 4C2
2ρ
35 − 2C33C55 + C2
33 + 4C2
2ρ

55)1/2

55)1/2

4.2 Quasiharmonic Approximation

Calculated lattice parameters for both LDA and GGA as well as their respective bulk

moduli determined from our QHA calculations, are given in Table 4.1.

85

a

LDA 12.181

b
3.029

c
5.772

β
103.73

GGA 12.276

3.054

5.813

103.71

B

0 K:
176
300 K: 172
950 K: 157

0 K:
161
300 K: 157
950 K: 136

Expt.

12.233

3.038

5.807

103.82

199

Table 4.1: Calculated lattice parameters (Å) at 0 K and bulk moduli, B (GPa), at T =
0, 300, 1000 K as calculated within the QHA. Results using LDA and GGA functionals as
well as the measured values of Ref. [73] are listed.

Unsurprisingly, lattice parameters calculated using GGA are over-estimated while the

LDA parameters are under-estimated. GGA over-prediction is a common phenomena,[33]

and in this case the LDA resolves a more accurate value for the bulk modulus just as it does
in the Al2O3 isostructure.[61] Overestimation of lattice constants is important to note as it

Figure 4.5: Linear thermal expansion coefficient for LDA (solid lines) and GGA (dashed
lines) calculations corresponding to applied strains εi as a function of temperature. Due to
symmetry, ε4 and ε6 are zero.

86

can lead to larger changes in the thermal expansion.

Linear thermal expansion for the applied strains was calculated for both GGA and LDA.

The results are summarized in Figure 4.5. As per our previously defined monoclinic orienta-
tion, [31] an applied strain ε2 is parallel to the b-direction and ε3 to the c-direction. Using
these definitions, the linear TEC in the b and c direction can be written as αb = αε2
and
. Similar values from experimental measurements of (αb = 3.37 × 10−6 K−1 and
αc = αε3
αc = 3.15 × 10−6 K−1) and (αb = 4.2 × 10−6 K−1 and αc = 4.2 × 10−6 K−1) have been
reported previously.[22, 54] In both measurements the lattice constants were fit as linear

functions of temperature, yielding single temperature-independent TEC values that are con-
sidered to be valid over approximately 100− 700K[54]. The LDA functional is found to more
accurately predicts measured values. This can be attributed to the GGA over-estimating

the lattice constants[33]. The origin of the anisotropy in the TEC is likely to be purely

geometric[54].

the heat capacity at constant pressure (Cp) and volume (CV

) are presented in Figure 4.6

and compared to [20] where good agreement is found.

Heat capacity in general is an important indicator of thermal performance in a semi-

conductor for use in power electronics. It is a parameter which represents the amount of

thermal energy which can be stored in the material, hence it directly effects the thermal

management and subsequent performance of the device. Comparing to Gallium Nitride, a

commonly used power electronic material with a thermal conductivity exceeding 200 W/mK,

the heat capacity is actually very similar. This is somewhat surprising when only considering

the phonon dispersion and crystal complexity. There is a significantly more complex crystal
structure in Ga2O3 which contains 30 phonon modes compared to GaN which only has 12
phonon modes due to the 4 atom hexagonal unit cell. The number of modes has a direct

impact on the total volumetric heat capacity, as the more phonon modes available, the more

energy can be stored. In fact, the heat capacities only differ by less than 20% which does not

explain the significant disparity (an entire order of magnitude) in thermal conductivities. It

87

Figure 4.6: Heat capacity for β − Ga2O3 at constant pressure (blue) and constant volume
(red) computed within the QHA compared to recently measured data of [20].

can thus be inferred that the average phonon lifetime in GaN is likely to be significantly
higher in GaN than Ga2O3 and this higher lifetime leads to more thermal energy storage
before phonon-phonon scattering occurs and heat is dissipated or transported through low

energy acoustic modes. This explanation is consistent with the phonon dispersion of GaN

which shows a significant bandgap. These gaps limit the anharmonic phase space for phonon-

phonon scattering events and thus create less resistance from phonon scattering and a higher

overall thermal conductivity regardless of the heat capacity [20].

4.2.1 Mechanical and elastic properties

The 13 independent elastic stiffness constants for Ga2O3 are listed in Table 4.2 for the
DFT equilibrium value (0 K) using both LDA and GGA exchange correlation functionals.
Our results are compared to the recent Cij calculations that used LDA[56] and GGA.[14]
GGA tends to underestimate Cij as compared to LDA since LDA also predicts smaller cell
sizes. [61, 62] The three principle elastic constants C11, C22 and C33 compare closely with

88

Ref. [14] Ref. [56] Ref. [55] Ref. [70]

This Work
LDA GGA GGA
223.1
237.9
347.9
333.2
330.0
356.7
50.3
66.7
79.5
68.6
94.2
112.4
116.5
139.1
118.8
125.3
75.0
92.6
-17.4
-1.8
5.6
12.2
7.3
7.8
9.1
17.4

222.0
327.3
333.6
58.8
78.3
108.4
125.2
112.9
79.9
-1.7
5.2
7.0
9.0

C11
C22
C33
C44
C55
C66
C12
C13
C23
C15
C25
C35
C46

LDA
237
354
357
54
67
95
125
147
95
-18
11
6
19

Expt.
238
359
346
49
91
107
130
152
78
-4
2
19
6

Expt
240
349
345
48
87
103
128
160
72
-1.6
0.3
1
5

θ − Al2O3 [61]
Expt
283
420
435
86
104
124
119
159
83
-30
12
16
23

Table 4.2: Comparison of calculated elastic stiffness constants Cij (GPa) of β − Ga2O3 .

the corresponding references. The remaining diagonal constants C44, and C66 are larger
within our approach while C55 is slightly smaller. The off-diagonal elastic constants show
significant variation among the reports.

There is a notable difference between the LDA and GGA functionals predicted results

which is an expected disparity. This can be attributed to the predicted lattice constants us-

ing each functional, where the LDA resolves a closer approximation to the measured lattice

constant. As such, it can be seen that the LDA approximation more accurately matches

the experimentally measured results. It is worth pointing out however that there is a rel-

ative spread in Table 4.2 within the off-diagonal portions of the tensor where the shearing

strains modulate the anisotropy of the structure. To resolve these tensor components, many

measurements must be performed on several different samples cut in precise configurations

to encapsulate the proper crystallographic orientation. This very complex measurement has

only very recently been achieved in the work of [70].

From these computed stiffness constants, the effect of anisotropy can be analyzed in

89

(a) (001)-plane
(b) (010)-plane
Figure 4.7: Crystallographic orientations of β − Ga2O3 planes

more detail. One notable feature which stands out is the significant disparity between the
C11 and C22 or C33 values. There is a nearly 30% decrease in the value of the stiffness
constant in the [100] direction implying a higher compressibility along the x-direction and

weaker chemical bonding along that direction. Referring to Figure 4.7a, it is seen that the

[100] direction aligns with the largest number of octahedrally coordinated Gallium atoms

(as opposed to the alternative tetrahedral sites). These octahedral polyhedra are more

stable and will resist deformation. Application of a strain along this direction will result

in less deformation of these polyhedra than along for example the [010] direction as seen

in Figure 4.7b. This will lead to a higher resistance to strain deformation along the [010]

which in turn results in a lower compressibility. This concept applies similarly to explaining

the anisotropy in thermal conductivity, where the [100] direction also exhibits the lowest

thermal conductance. This is consistent with the idea that a higher compressibility along this

direction will allow for a lower frequency oscillation to occur in atomic vibrations due to the

lower rigidity of the bonds. The octahedrally coordinated Gallium atoms can vibrate more

freely with less returning forces due to the increased ability to compress without perturbing

the coordination polyhedra. To illustrate the significant anisotropy, the Young’s modulus

ellipsoid is projected along each plane and presented in Figure 4.8. This anisotropy is

90

(a) (yz)-plane

(b) (xz)-plane

(c) (xy)-plane

Figure 4.8: Tensor projection of Youngs modulus of β − Ga2O3

important to consider for engineering purposes where directional phenomena can be utilized

to control material properties. Conversely, it is also important to be aware of for unintended

effects as may be the case when forming thin-film surfaces along certain orientations. It is
seen in the figure that the [100] or x direction contains a significantly lower Young’s modulus
value then the [010] and [001] directions. This behavior is consistent with the anisotropy

also calculated within the elastic stiffness constants. According to the results of [70], this
behavior is also found in θ−Al2O3 , which possess the same spacegroup. However, for similar
monoclinic structures such as BiMnO3 this degree of anisotropy is not seen. This alludes to
the anisotropy not being due to geometry alone and likely has more to due with the nature

of the chemical bonding. This discrepancy is investigated in further detail by considering the

linkage of the coordination polyhedron of each material. In Figure 4.9 the polyhedral linkage

is illustrated and it can be seen that among the two distinct species of Gallium atoms, the

octahedrally coordinated variety has a unique linkage network along the [100] direction which

is not seen along the other primary crystallographic directions. Specifically, the octahedral

polyhedron share edge with one another along [100] which will likely lead to higher binding
energies. Comparing to the BiMnO3 system with similar crystallography, the linkage is vastly
different as shown in Figure 4.10. None of the octahedrally coordinated cations’ polyhedron

91

Figure 4.9: Monoclinic β − Ga2O3 polyhedral linkage illustrating the unique bonding mech-
anism between octahedrally coordinated Gallium atoms. Edges of the octahedrally coordi-
nated Gallium atoms’ polyhedrons are shared.

have any shared edges or faces. This phenomena is analyzed in more detail in [70], where the

authors speculate that the [100] direction is able to produce a longitudinal deformation with

less bond length changes among cation-anion due to the flexibility of the linkage to bend

between tetrahedron and octahedron. It was also revealed in this work that the bondlength

change along this direction to achieve a set value of strain is notably smaller than that of
the [010] and [001] directions. This behavior is not seen in BiMnO3 which is further verified
by the projections of the Youngs’ Modulus ellipsoid presented in Figure 4.11 which shows

very little relative anisotropy.

It is also interesting to note the anisotropy reported on here matches the anisotropy

seen within the thermal conductivity tensor, where the [100] direction demonstrates a no-

tably smaller thermal conductivity (nearly half of the [010]). This anisotropy is visible in

the phonon dispersion and leads to less acousto-optical gaps in the [100] which leads to an

increased likelihood of anharmonic scattering events which provide thermal resistance. To

consider the finite temperature effects that may be present in the elastic moduli, the isother-
mal and isentropic Cij for β − Ga2O3 were predicted and are plotted in Figure 4.12 as a

92

Figure 4.10: Monoclinic BiMnO3 polyhedral linkage contrasts the behavior seen in β −
Ga2O3 regardless of the material possessing the same crystal type. There are no shared
edges among the octahedrally coordinated cations’ polyhedron which is predicted to be the
cause of the anisotropy difference.

(a) (yz)-plane

(b) (xz)-plane

(c) (xy)-plane

Figure 4.11: Tensor projection of Youngs modulus of Monoclinic BiMnO3 illustrating the
significantly lower anisotropy compared to similar monoclinic structure β − Ga2O3

93

Figure 4.12:
(left) and GGA (right) functionals for β − Ga2O3 .

Isothermal (solid) and isentropic (dashed) elastic stiffness constants for LDA

function of temperature for both GGA and LDA functionals. The GGA calculated values are
consistently lower when compared to LDA, as expected. The predicted isentropic Cijs are
almost always greater than the isothermal values. The differences in magnitude of isothermal
and isentropic Cij becomes significant above 400 K. The isentropic and isothermal stiffness
constants C44, C55, and C66 are identical for the entire temperature range considered here.
Using Cij and the corresponding elastic compliance Sij, the bulk modulus (B), shear
) were computed and are shown

modulus (G), Poisson’s ratio (ν) and Vicker’s hardness (HV
in Table 4.3 at 0 K for Ga2O3 using both the LDA and PBEsol (GGA) functionals.

Figure 4.13 shows the temperature dependencies of the elastic moduli computed qua-

sistatically. None of the reported values exhibit significant temperature dependencies. For
example, from 0 K to 950 K the bulk modulus only decreases about 10 GPa for LDA and

94

LDA

GGA

Measured (Ref. [74])

B
182
179
180

168
165
167

Voigt:
Reuss:
Hill:

Voigt:
Reuss:
Hill:

ν
0.29
0.30
0.29

0.28
0.30
0.29

HV
13.4
12.2
12.8

12.7
11.5
12.0

G E
91
83
87

234
216
225

86
78
82

222
203
212

233

Table 4.3: Calculated equilibrium values for bulk modulus B, shear modulus G, Young’s
modulus (E), Vicker’s hardness (HV
), and Poisson’s ratio (ν). All values determined by Cij
and are reported in GPa except the dimensionless ν.

LDA

GGA

[100]
[010]
[001]

vt
3.61 (2.87)
7.54 (7.51)
3.23 (2.86)

v+
6.23 (6.48)
3.71 (3.64)
7.63 (7.15)

v−
4.28 (3.88)
3.12 (2.02)
4.28 (3.61)

vt
3.63 (2.47)
7.39 (7.27)
3.04 (2.72)

v+
6.08 (6.32)
3.70 (3.59)
7.49 (7.09)

v−
4.25 (3.84)
2.96 (1.96)
4.25 (3.39)

Table 4.4: Sound velocities in β − Ga2O3 along different crystallographic directions for the
transverse (vt), fast mixed mode (v+), and slow mixed mode (v−). Velocities are given in
units of ×105 cm/s. Comparisons to the lattice dynamical results extracted directly from
the phonon dispersion are given in parentheses.

18 GPa for GGA. Comparing to Table 4.1, the temperature dependence of the bulk modulus
using the QHA is noticeably larger, decreasing 15 GPa from 300K to 900K as opposed to

the 7 GPa decline in the elastic stiffness approach.

Table 4.4 shows the sound velocities in β − Ga2O3 evaluated from the elastic constants
along the crystallographic directions.[69] In all the three directions there is a pure transverse
mode (vt) and two mixed modes (v+ and v−).
In the [100] and [001] directions vt has
the lowest velocity and the mixed v+ mode has the highest. Only in the [010] direction
does the transverse mode have the highest energy. A comparison with the sound velocities

evaluated directly from lattice dynamical studies (phonon band structure) is also included

for completeness.

The Debye temperature in Andersons’s model [75] is evaluated from the average sound

95

Figure 4.13: Temperature dependence of bulk moduli determined from computed elastic
stiffness constants up to 1400 K.

vL

LDA 6.95
GGA 6.79

vt
3.77
3.70

vavg ΘD
4.21
4.13

723 K
557 K

Table 4.5: Calculated bulk sound velocity and computed Debye temperature (ΘD
longitudinal (vL
ΘD

). The
) and transverse mode (vt) modes were averaged (vavg) and used to compute

(K). Velocities are given in ×105 cm/s.

velocity for both GGA and LDA calculations and is shown in Table 4.5. The average sound

velocity was obtained from the longitudinal and transverse sound velocities which were eval-

uated from the elastic moduli.

4.2.2 Elastic Moduli α − Ga2O3

The α phase of Ga2O3 has a simpler corundum symmetry and the mechanical properties are
generally easier to compute. However very little is known about the mechanical properties

96

of this phase, with no estimates of the elastic properties at the time of this writing. It can
be deduced however that the properties should resemble that of sapphire (α−Al2O3 ) which
shares a crystal type. In addition to the shared crystal geometry, sapphire is the material
used as a substrate for epitaxial layer growth of α − Ga2O3 . The lattice mismatch between
these structures is relatively low (4%) so it is natural to expect geometric properties such as

anisotropy to be very similar between these structures. The elastic stiffness tensor of this

structure is expressed generally as

rho
ij =

C



c11

c12

c11

c13
c14
c13 −c14
0
c33

c44

0

0

0

0

c44

0

0

0

0

c14

c11−c12

2



(4.21)

This tensor has only six independent components and the form is heavily dependent on the

orientation of the crystallographic axes. The most common definition is that the z-axis is

parallel to the rhombohedral [111] direction which is also equivalent to the c-axis in the

hexagonal representation. In the LDA approximation, the elastic constants are computed

and compared to sapphire in Table 4.6

α − Ga2O3

364
325
93
128
101
8
180
113

C11
C33
C44
C12
C13
C14
B
G

α−Al2O3 [61]

497
493
154
165
129
19
249
155

Table 4.6: Calculated elastic stiffness (GPa) constants using LDA compared to α−Al2O3

Comparing to the monoclinic phase, the elastic constants are considerably larger. This

is consistent with the decreased volume of the unit cell in the rhombohedral setting, leading

97

(a) (yz)-plane

(b) (xz)-plane

(c) (xy)-plane

Figure 4.14: Tensor projection of Youngs modulus of Rhombohedral α − Ga2O3 . Signifi-
cantly lower anisotropy is seen compared to the monoclinic phase β − Ga2O3

to a lower compressibility and subsequently higher bulk modulus. Additionally, there is
a tighter packing of oxygen cations into the lattice for the α phase due to the hexagonal
symmetry. Similar to the β − Ga2O3 case, the Aluminum isostructure for the α phase shares
the overall form however contains larger magnitudes and a larger bulk modulus. This is

once again likely due to the smaller occupied volume of the unit cell. One interesting feature
that is distinct is the value of C11, which has the highest magnitude for the Al2O3 case
while α − Ga2O3 has a maximum value in the C33 constant. Additionally, the anisotropy
is slightly different, with α − Ga2O3 containing a slightly larger anisotropy in the normal
directions, i.e C11/C33 = 1.12 for α − Ga2O3 and C11/C33 = 1.12 for Al2O3 . Compared to
the monoclinic phase however, the overall anisotropy is also reduced significantly in the 3D

representation of the Youngs modulus, which can be seen in Figure 4.14

4.3 Optical Properties

Due to the large bandgap of β − Ga2O3 , there is significant interest in the optical prop-
erties which can be expected of the material. Characterizing these properties is considerably

more complex for a material with large anisotropy as is the case for monoclinic crystals. The

anisotropy of the crystal relates the optical properties to more complex phenomena such

98

as birefringence and double refraction, wherein electromagnetic waves of different character

and velocity can simultaneously propagate through the same crystallographic direction. This

leads to a range of possible refractive indices which depend on the direction of propagation.

This relationship is often expressed as an ellipsoid which is oriented along the principal axes

of the second-rank tensor (in this case, dielectric function) with the refractive indices being

the semi-axes of the ellipsoid. This common representation of the index of refraction is re-

ferred to as the optical indicatrix and is generally a complicated construct to represent the

anisotropy of any material property tensor [63]. This complexity is furthered extensively in

situations where the already anistropic crystal is subjected to an external influence such as

electric field or strain. However, the effect of geometry also has the potential to simplify the

anisotropy of a material property in many cases which has been observed experimentally for
thin-film β−Ga2O3 [60], where the tensor character of the dielectric function was reduced to
an orthorhombic nature as the film width decreased. This shows the complexity of the indi-

catrix and the need to understand how the material properties will modulate under external

influence. Strain in particular is an important parameter to consider due to the infancy of

the growth methods and the likelihood of residual strain existing in the crystal post-growth.

The effect of these residual strains can significantly influence the optical properties of a semi-

conductor which has been demonstrated in GaAs for example [76]. Experimentally this is a

considerable challenge to characterize and measure.

4.3.1 Dielectric Function

The dielectric function is most commonly described by the sum of the real and imaginary

parts

ε(ω) = ε1(ω) + iε2(ω)

(4.22)

where the imaginary part ε2(ω) relates to the energy dissipation as a wave travels through
the medium and the real part ε1(ω) relates to the energy stored by the medium. As
Ga2O3 belongs to the biaxial class of crystals, the dielectric tensor 
(ω) has several principal

99

optical axes leading to birefringence. This has recently been reported on extensively[60] with

particular emphasis placed on how the orientation of optical axes changes with respect to

incident photon frequency. In the orientation used here the second rank dielectric tensor

takes the form

ε(ω) =

εxx

0

εxz

0

εyy

0

εxz

0

εzz



(4.23)

Several recent works have focused on how influential excitonic contributions to 
(ω) may
be in the visible frequency range. For the case of Ga2O3 there has been a sizable correction
to the onset frequency when excitonic corrections are added to calculations.[77] These calcu-

lations require extensive computational resources and the influence of strain on the method’s

validity has not been investigated for even simple cases thus far, so they were not included in

the current work. In the computations in this work the Drude-Lorentz model[63] was solved
for the imaginary component of the dielectric function, 
2 which entails the summing of tran-
sitions occurring between the unoccupied and occupied electronic states and is expressed as

ε2(ω) =(cid:18)4π2e2

mω2 (cid:19)(cid:88)ij (cid:90) (cid:104)i|M|j(cid:105)2fi (1 − fi) × δ(cid:0)Ej,k − Ei,k − ω(cid:1) d3k

where i is the initial state and j is the final state, fi is the Fermi-Dirac distribution, and M
is the dipole matrix. This is followed by the Kramers-Kronig transformation to determine
the real part, 
1:

(4.24)

(4.25)

ε1(ω) = 1 +

2
(π)

P(cid:90) ∞

0

ω(cid:48)ε2(cid:0)ω(cid:48)(cid:1) dω(cid:48)
(cid:0)ω(cid:48)2 − ω2(cid:1)

The imaginary part of dielectric function is evaluated within the LDA. The effect of temper-

ature is not considered while analyzing the optical properties.

100

Figure 4.15: Computed Real (left) and imaginary (left) components of dielectric function
of β − Ga2O3 without external strains

4.3.2 Normal Strain

The application of a normal strain does not reduce the overall symmetry of the crystal. The
normal strain in the (cid:126)x for example is of the form,


 =

η11 0 0

0

0

0 0

0 0



which transforms the lattice vectors as follows,

(4.26)

(4.27)

(4.28)

(4.29)

(4.30)

a = am(1 + η11),

b = bm

c = cm(cid:113)1 + 2η11 cos2 βm + η2
(cid:113)1 + 2η11 cos2 βm + η2

(1 + η11) cos βm

cos β =

11 cos2 βm,

11 cos2 βm

The effect of these strains along the cartesian directions does not show a significant effect

on the anisotropy or magnitude of the dielectric function even for relatively large strains of

0.75%. The results are summarized in Figure 4.16

101

Figure 4.16: Dielectric function (real and imaginary) of β − Ga2O3 with normal external
strains in the (cid:126)x (top), (cid:126)y (middle) and (cid:126)z (bottom) directions.

4.3.3 Shearing Strain

Applying strain in certain directions breaks the symmetry of the crystal and reduces the
lattice to triclinic symmetry. For example a shearing strain in the xy-plane is of the form,

(4.31)

η =

0

η12 0

η12

0

0

0

0

0



102

which transform the lattice parameters as,

cos α =

1

12,

b =

m + b2

a = am(cid:113)1 + η2
2(cid:113)(a2
c = cm(cid:113)1 + η2
(cid:113)1 + η2
cos βm(cid:113)1 + η2
(cid:113)1 + η2
(cid:113)1 + η2

12(cid:113)(a2

12
12 cos2 βm

cos β =

cos γ =

m)(1 + η2

12) + 4η12ambm,

12 cos2 βm,
cos βm(am(1 + η2
m + b2

12 cos2 βm(cid:113)(a2

12) + 2η12bm)
m)(1 + η2

12) + 4η12ambm

am(1 + η2
m + b2

12) + 2η12bm
m)(1 + η2

12) + 4η12ambm

Specifically, the shearing strains produce non-negligible changes to the form of the dielectric

tensor. This results in significant changes in anisotropy as seen in Figure 4.17.

In both the unstrained and strained cases the results are highly anisotropic and exhibit

dichroism. However, unlike the unstrained situation, the strained configuration exhibits

sizable dichroism at nearly steady-state frequencies.

4.3.4 Photo-elastic Effect

The strain-dependence of the anisotropy in the dielectric tensor is an important phenomena

for a number of applications. For example, one of the more robust methods of characterizing

detrimental extended defects which occur during crystal growth is by means of birefringence

optical microscopy [78]. To quantitatively analyze using this method, knowledge of the

photoelastic constants is necessary. The large strain dependence on optical properties also
potentially opens up an unexplored application space for Ga2O3 in the field of acousto-
optic modulator devices. These devices can be used to control power, frequency, and spatial

directions of lasers by modification of the refractive index using acoustic waves. Recently,

acousto-optic modulators have become more sophisticated and are included on-chip using
integrated optics on materials such as lithium niobate (LiNbO3). Devices of this nature

103

Figure 4.17: Dielectric function of β−Ga2O3 with applied shearing strains (solid lines) along
with the unstrained calculation (dashed lines). There is significantly increased anisotropy
when strain is applied, increasing the brirefringence. The strongly fluctuating anisotropy
with respect to frequency also indicates dichroism in the refractive index.

can be used in numerous ways such as tunable optical filters and optical switches. The

response of the change in inverse dielectric function to strain is referred to as the elasto-

optic, or photoelastic effect. It can be derived from the optical indicatrix ellipsoid mentioned
previously by defining the impermeability Bij and index of refraction nij,

where index ellipsoid is written in the standard notation [63],

1
Bij

=

1
εij

=

1
n2
ij

The change in this impermeability gives rise to the photoelastic effect,

Bijxixj = 1

∆Bij =(cid:16)ε

−1

ηkl(cid:17)ij

−(cid:16)ε

−1

0 (cid:17)ij

= pijklηkl

104

(4.32)

(4.33)

(4.34)

The change in the inverse dielectric function can also be expressed by

ε(cid:19)ij
∆(cid:18)1

= − ∆εij
εiiεjj

=(cid:88)kl

The full strained ellipsoid can now be expressed by

pijklηkl

(cid:0)Bij + pijklηkl(cid:1) xixj = 1

(4.35)

(4.36)

is used to represent the strain vectors. Using this definition, the process of comput-

where ηkl
ing the photoelastic tensor pijkl
is the stress response with the change in inverse dielectric function that is calculated. Due to

is identical to that of the elastic tensor calculation except it

is not easily represented using the more compact voigt notation. Here the most

the lack of symmetry in the photoelastic tensor in the monoclinic case, the full fourth-rank
tensor pijkl
pronounced tensor elements of the unitless pijkl
tion of dielectric properties is given. The response due to normal applied strain is overall the
largest: p1111 = 8 × 10−2 , p2222 = −31 × 10−2 , p3333 = −56 × 10−2 , p1122 = −30 × 10−2
, p1133 = −35 × 10−2 and p3322 = −27 × 10−2 .

quantity which result in the largest modula-

Interestingly, the response due to shearing strains can be significant as well. The non-
negligible values being p1113 = 24 × 10−2, p2213 = 4 × 10−2, p3313 = −10 × 10−2, and
p1333 = 7 × 10−2. These sizable off-diagonal elements, especially p1113, are expected as
these strains correspond to symmetry breaking where the dielectric response should be more

pronounced as was shown in Figure 4.17.

A material’s performance for use in acousto-optical devices can be roughly estimated as
a figure-of-merit (M ) value.[79] This is essentially a measure of the fraction of light deflected
at Bragg incidence and can be represented by

M =

p2maxn6
ρv3ac

(4.37)

where n is the index of refraction, pmax is the maximum photoelastic tensor component, ρ
is the density, and vac is the acoustic wave velocity. This value is often normalized to the
value M0 = 1.51× 10−15 [79]. Utilizing these parameters as calculated and presented earlier,

105

the figure of merit is determined for Ga2O3 to be M/M0 = 6.72. This compares well to the
common acousto-optic material LiNbO3 (M/M0 = 1.8)[79] and demonstrates a promising
and unexplored application space of β − Ga2O3 materials.

Another method of characterizing the effect of strain is by estimating the birefringence

parameter. In general, birefringence stems from the polarization dependence of the refractive

index and subsequently the optical axes of the index ellipsoid. The birefringence can be

characterized by the expression for an angular rotation of the optical axis defined by

tan 2φ =

2εxz

εxx − εzz

(4.38)

Using this expression, at DC (ω = 0) the angle of optical axes from the principal axes is
φ unstrained ∼ 7◦ while the application of a shearing strain increases this value to θ sheared ∼
20◦. This demonstrates a significant dependence on strain even for DC frequencies. The
orientation of the dielectric axes can play a significant role in polarization sensitive applica-

tions.

106

CHAPTER 5

CONCLUSION

This dissertation successfully characterized an array of thermal, elastic, and optical material
properties of both β − Ga2O3 and α − Ga2O3 using only first-principles techniques. This
serves to validate the available measured data on these materials as well as predict prop-

erties which are not currently well understood experimentally. Additionally, this work sets

a benchmark for the robustness of first-principles computation techniques in thermal trans-

port of bulk and layered structures for materials of very low symmetry. In this chapter, the

results of this work will be briefly summarized.

5.1 Summary

Both α − Ga2O3 and β − Ga2O3 are ultra-wide bandgap semiconductors with very
high breakdown fields that has seen applications in high power electronics as well as UV-
transparent conductive films. The application space continues to grow especially for β −
Ga2O3 devices, which has seen significant improvement in crystal growth quality and tech-
nique in recent years. Although the material claims many beneficial material properties and a

wide bandgap, the low symmetry nature has led to a difficult to characterize anisotropy that

is not well understood experimentally. This dissertation has performed a detailed analysis of

the electronic structure using first-principles DFT calculations to resolve the bandstructure
for both α − Ga2O3 and β − Ga2O3 in the LDA, GGA, and hybrid functional frameworks.
This leads to an accurate determination of bandgap and effective mass which agree well

with measured data in both magnitude and anisotropy. The accurate electronic structure

calculations then enabled the investigation of thermal properties, specifically thermal con-

ductivity. It has been speculated for some time that the thermal conductivity will be the
largest bottleneck in the development of more sophisticated Ga2O3 based devices. The value
is very low (around 21 W/mK) compared to other wide bandgap semiconductors such as

107

GaN, diamond, or SiC and creates a thermal management issue in realistic applications such

as on-chip power converters. Using the recently developed finite-displacement supercell ap-

proach, this work computes harmonic and anharmonic phonon properties such as phonon

dispersion, phonon-phonon scattering rates, and thermal conductivity with no empirical in-

puts. The magnitude of thermal conductivity tensor was confirmed to be low as predicted

from the limited measurements available ranging from 13-22 W/mK. The anisotropy is sig-

nificant with a factor of 2 difference between the [100] and [010] crystallographic directions.

The anisotropy was determined primarily to stem from a few different factors, primarily the

importance of a non-negligible group velocity of optical phonon modes in the [010] direction,

leading to an overall contribution of 44% of the total thermal conductivity in this direction.

Additionally, the existence of an acousto-optical gap in the phonon dispersion along the [010]

also limits the scattering phase space and thus increases thermal conductance. Additionally,
α − Ga2O3 was investigated briefly using the same first-principles methodology and was
found to have an even lower thermal conductivity (around 10 W/mK) but more isotropic.
The effect of the film thickness in the α phase was found to be relevant to thermal conduc-
tivity with a significant increase in thermal resistance due to boundary effects for films less

than 100 nm.

Managing the poor thermal conductivity of Ga2O3 is a crucial issue to facilitate further
refinement of Ga2O3 devices. One approach, which was explored extensively in this work,
is by means of forming an effective thermal contact with some higher thermal conductivity

material to act as a heat-sink. By using recently developed first-principles techniques com-

bined with a variance-reduced Monte-Carlo solver, the thermal boundary conductance at a
Au/β − Ga2O3 interface was computed and compared with the setups containing an adhe-
sion layer composed of chromium or titanium to increase phonon transmission. The results

matched very well with recent measurements and both chromium and titanium were deter-

mined to be effective thermal contact layers to increase the thermal boundary conductance

significantly. Chromium had the maximum effect on thermal boundary conductance and

108

this is attributed to the high overlap in the phonon density of states with β − Ga2O3 , which
will lead to a higher phonon transmission within the diffuse mismatch model where only

elastic mode conversions are allowed. More generally, this calculation serves as a benchmark

for the VRMC method, and demonstrates the feasibility of the workflow for realistic device

geometries.

Lastly, the effect of strain on the mechanical and optical properties was characterized.

While the material continues to mature and growth methods become more refined, the effect

that residual strains may have on the crystal is important to consider. In this work, the elastic
stiffness tensor for both α and β phases were computed within a first-principles framework
and the anisotropy of the tensor was evaluated. The results matched very well with recently
measured results and the distinct anisotropy in the β phase was explained through the unique
polyhedral linkage that is present from the existence of multiple coordination complexes

among the cations that is not present in other crystals of similar symmetry. Furthermore, the

quasiharmonic approximation is utilized which incorporated the phonon free-energy functions

of multiple geometric configurations of varying volume. This allowed for an approximation for

the temperature dependent elastic moduli to be determined and extended to the temperature

dependence of the elastic stiffness constants. Finally, the dielectric function is computed as a

function of frequency under a fully spanned set of strains for the crystal. The shearing strains

produce a sizable modulation of anisotropy within the dielectric tensor due to the crystal

symmetry breaking. This modulation was characterized by means of the photo-elastic tensor,
which is then used to compute a figure-of-merit to determine the feasibility of β − Ga2O3 as
an acousto-optic modulator device material. It was determined that the figure of merit of
β − Ga2O3 is sizable and exceed that of lithium niobate, which is a workhorse material for
acousto-optial materials.

109

5.2 Future Work

There are still many obstacles that will need to be overcome to facilitate further improve-

ments in device design. From the results presented in this work, the most immediately useful

topics to continue research on are outlined in this section. Firstly, while the thermal bound-
ary conductance was determined to be large and overall the potential for β − Ga2O3 to
form good thermal contacts is high, this does not alleviate the fact that within the bulk
β − Ga2O3 material itself the thermal conductivity is quite low compared to the other WBG
semiconducting materials. The self-heating issues will not be inherently solved solely by the

inclusion of a good heatsink. It does however imply that by maximizing the area exposed to

a heatsink and minimizing the overall layer thickness, there is likely to be a vastly improved
over thermal performance of a β − Ga2O3 device structure. In the future, developing an
improved Monte-Carlo algorithm allowing for a transient evaluation of temperature profiles

would provide a more complete picture of the device geometries self-heating behavior. The

next topic would be a more complete evaluation of the strain dependent index of refraction.

Including not just the strain-dependent dielectric tensor in the form of the photo-elastic ten-

sor. But by also computing the electro-optic tensor, which characterizes the response of the

dielectric tensor to a static or low-frequency electric field. By computing these properties

within a DFT framework, the electronic and phononic contributions can be isolated and

the influence of phonons, the crystal structure, temperature, and bonding character between

individual atoms can be weighted by their contributions to the optical properties of the

crystal. Lastly, a topic which was not covered in this work, would be an evaluation of the

electron-phonon scattering rates from first-principles. There has been preliminary evidence

that polar optical phonon scattering is the dominant mechanism limiting electron mobility
in β − Ga2O3 [80]. For example, the electron mobility is 10x lower than GaN while the
effective mass is very similar. This was attributed to a large Frohlich interaction stemming

from the low energy phonons. Detailing this interaction using a more complete theoreti-

cal model, using the electron-phonon wannier interoplation method, would be help isolate

110

the phonon modes most responsible for this limited mobility. Overall, significant progress
is being made on the development of Ga2O3 technologies and the work presented in this
dissertation serves to provide a more fundamental understanding of the underlying material

properties and what bottlenecks must still be overcome in the future.

111

BIBLIOGRAPHY

112

BIBLIOGRAPHY

[1] JY Tsao, S Chowdhury, MA Hollis, D Jena, NM Johnson, KA Jones, RJ Kaplar, S Ra-
jan, CG Van de Walle, E Bellotti, et al. Ultrawide-bandgap semiconductors: Research
opportunities and challenges. Advanced Electronic Materials, 2017.

[2] Masataka Higashiwaki, Kohei Sasaki, Hisashi Murakami, Yoshinao Kumagai, Akinori
Koukitu, Akito Kuramata, Takekazu Masui, and Shigenobu Yamakoshi. Recent progress
in Ga2O3 power devices. Semiconductor Science and Technology, 31(3):034001, 2016.
[3] JB Varley, Anderson Janotti, Cesare Franchini, and Chris G Van de Walle. Role of self-
trapping in luminescence and p-type conductivity of wide-band-gap oxides. Physical
Review B, 85(8):081109, 2012.

[4] A Segura, Lluís Artús, Ramón Cuscó, R Goldhahn, and M Feneberg. Band gap of
corundumlike α- ga 2 o 3 determined by absorption and ellipsometry. Physical Review
Materials, 1(2):024604, 2017.

[5] Kentaro Kaneko, Shizuo Fujita, and Toshimi Hitora. A power device material of
corundum-structured α-ga2o3 fabricated by mist epitaxy ® technique. Japanese Jour-
nal of Applied Physics, 57(2S2):02CB18, 2018.

[6] SJ Pearton, Jiancheng Yang, Patrick H Cary IV, F Ren, Jihyun Kim, Marko J Tadjer,
and Michael A Mastro. A review of ga2o3 materials, processing, and devices. Applied
Physics Reviews, 5(1):011301, 2018.

[7] A. Togo, F. Oba, and I. Tanaka. First-principles calculations of the ferroelastic transition
between rutile-type and CaCl2-type SiO2 at high pressures. Phys. Rev. B, 78:134106,
Oct 2008.

[8] Paolo Giannozzi et al. Quantum espresso: a modular and open-source software project

for quantum simulations of materials. J. Phys.: Condensed Matter, 21(39), 2009.

[9] Hartwin Peelaers and Chris G. Van de Walle. Brillouin zone and band structure of

β − Ga2O3. physica status solidi (b), 252(4):828–832, 2015.

[10] Kristina E Lipinska-Kalita, Patricia E Kalita, Oliver A Hemmers, and Thomas Hart-
mann. Equation of state of gallium oxide to 70 gpa: Comparison of quasihydrostatic
and nonhydrostatic compression. Physical Review B, 77(9):094123, 2008.

[11] Richard M Martin. Electronic structure: basic theory and practical methods. Cambridge

university press, 2004.

[12] Zbigniew Galazka et al. On the bulk β−Ga2O3 single crystals grown by the czochralski

method. J. Cryst. Growth, 404(0):184 – 191, 2014.

113

[13] Julia E Medvedeva and Chaminda L Hettiarachchi. Tuning the properties of complex
transparent conducting oxides: Role of crystal symmetry, chemical composition, and
carrier generation. Physical Review B, 81(12):125116, 2010.

[14] J Furthmüller and F Bechstedt. Quasiparticle bands and spectra of ga 2 o 3 polymorphs.

Physical Review B, 93(11):115204, 2016.

[15] Haiying He, Roberto Orlando, Miguel A Blanco, Ravindra Pandey, Emilie Amzallag,
Isabelle Baraille, and Michel Rérat. First-principles study of the structural, electronic,
and optical properties of ga 2 o 3 in its monoclinic and hexagonal phases. Physical
Review B, 74(19):195123, 2006.

[16] Joachim Paier, Martijn Marsman, K Hummer, Georg Kresse, Iann C Gerber, and
János G Ángyán. Screened hybrid density functionals applied to solids. The Jour-
nal of chemical physics, 124(15):154709, 2006.

[17] Jong-Won Song, Koichi Yamashita, and Kimihiko Hirao. Communication: A new hybrid
exchange correlation functional for band-gap calculations using a short-range gaussian
attenuation (gaussian-perdue–burke–ernzerhof), 2011.

[18] Lucy D Whalley, Jarvist M Frost, Benjamin J Morgan, and Aron Walsh. Impact of
non-parabolic electronic band structure on the optical and transport properties of pho-
tovoltaic materials. arXiv preprint arXiv:1811.02281, 2018.

[19] Sean Knight, Alyssa Mock, Rafał Korlacki, Vanya Darakchieva, Bo Monemar, Yoshinao
Kumagai, Ken Goto, Masataka Higashiwaki, and Mathias Schubert. Electron effective
mass in sn-doped monoclinic single crystal β-gallium oxide determined by mid-infrared
optical hall effect. Applied Physics Letters, 112(1):012103, 2018.

[20] Zhi Guo, Amit Verma, Xufei Wu, Fangyuan Sun, Austin Hickman, Takekazu Masui,
Akito Kuramata, Masataka Higashiwaki, Debdeep Jena, and Tengfei Luo. Anisotropic
thermal conductivity in single crystal β-gallium oxide. Applied Physics Letters,
106(11):111909, 2015.

[21] Z. Galazka M. Handwerg, R. Mitdank and S.F. Fischer. Temperature-dependent thermal
conductivity in mg-doped and undoped ga2o3 bulk-crystals. Semi. Sci. and Tech., 30(2),
2015.

[22] Encarnación G. Víllora, Kiyoshi Shimamura, Takekazu Ujiie, and Kazuo Aoki. Electrical
conductivity and lattice expansion of ga2o3 below room temperature. Applied Physics
Letters, 92(20), 2008.

[23] G. P. Srivastava H. M. Tutuncu, Subhash L. Shinde. Length-Scale Dependent Phonon
Interactions. Topics in Applied Physics 128. Springer-Verlag New York, 1 edition, 2014.
[24] Wu Li, Jesús Carrete, Nebil A. Katcho, and Natalio Mingo. Shengbte: A solver of the
boltzmann transport equation for phonons. Comp. Phys. Comm., 185(6):1747 – 1758,
2014.

114

[25] Alistair Ward. First principles theory of the lattice thermal conductivity of semiconduc-

tors. PhD thesis, Boston College, 2009.

[26] Keivan Esfarjani and Harold T. Stokes. Method to extract anharmonic force constants

from first principles calculations. Phys. Rev. B, 77:144112, Apr 2008.

[27] A. Ward and D. A. Broido. Intrinsic phonon relaxation times from first-principles studies

of the thermal conductivities of si and ge. Phys. Rev. B, 81:085205, Feb 2010.

[28] Ruiqiang Guo, Xinjiang Wang, and Baoling Huang. Thermal conductivity of skutteru-
dite cosb3 from first principles: Substitution and nanoengineering effects. Sci. Rep., 5,
Jan 2015.

[29] L Lindsay, DA Broido, and TL Reinecke. Thermal conductivity and large isotope effect

in gan from first principles. Physical review letters, 109(9):095901, 2012.

[30] Mitsuaki Kawamura, Yoshihiro Gohda, and Shinji Tsuneyuki. Improved tetrahedron
method for the brillouin-zone integration applicable to response functions. Physical
Review B, 89(9):094515, 2014.

[31] Marco D. Santia, Nandan Tandon, and J. D. Albrecht. Lattice thermal conductivity in

β-ga2o3 from first principles. Applied Physics Letters, 107:041907, 2015.

[32] P. E. Blöchl. Projector augmented-wave method. Phys. Rev. B, 50(0):17953–17979,

1994.

[33] J. P. Perdew, K. Burke, and M. Ernzerhof. Generalized gradient approximation made

simple. Phys. Rev. Lett., 77(0):3865–3868, 1996.

[34] S. Geller. Crystal Structure of Ga2O3. J. Chem. Phys., 33:676–684, September 1960.
[35] S. Baroni, S. de Gironcoli, and A. D. Corso. Phonons and related crystal properties

from density-functional perturbation theory. Rev. Mod. Phys., 73:515–562, 2001.

[36] J. M. Ziman. Electrons and phonons : the theory of transport phenomena in solids.

Clarendon Press Oxford University Press, Oxford New York, 2001.

[37] Aleksandr Chernatynskiy and Simon R. Phillpot. Phonon transport simulator (phonts).

Comp. Phys. Comm., 192(0):196 – 204, 2015.

[38] L. Lindsay and T. L. Reinecke. Three-phonon phase space and lattice thermal conduc-

tivity in semiconductors. J. Phys.: Condens. Matter, 20(0):165209, 2008.

[39] Tianli Feng and Xiulin Ruan. Prediction of spectral phonon mean free path and thermal
conductivity with applications to thermoelectrics and thermal management: A review.
J. Nanomaterials, 2014:1–25, 2014.

[40] J. B. Varley, J. R. Weber, A. Janotti, and C. G. Van de Walle. Oxygen vacancies and

donor impurities in ga2o3. Appl. Phys. Lett., 97(14), 2010.

115

[41] Encarnación G. Víllora, Kiyoshi Shimamura, Yukio Yoshikawa, Takekazu Ujiie, and
Kazuo Aoki. Electrical conductivity and carrier concentration control in ga2o3 by si
doping. Appl. Phys. Lett., 92(20), 2008.

[42] H. Aller, X. Yu, A. J. Gellman, J. A. Malen, and A. J. H. McGaughey. Thermal
conductance of ga2o3/metal interfaces. In 2018 17th IEEE Intersociety Conference on
Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm), pages
567–571, May 2018.

[43] Adam J. Jackson and Aron Walsh. Oxidation of gan: An ab initio thermodynamic

approach. Phys. Rev. B, 88:165201, Oct 2013.

[44] Marco Arrigoni and Georg KH Madsen. Comparing the performance of lda and gga
functionals in predicting the lattice thermal conductivity of iii-v semiconductor materials
in the zincblende structure: The cases of alas and bas. Computational Materials Science,
156:354–360, 2019.

[45] Jesús Carrete, Bjorn Vermeersch, Ankita Katre, Ambroise van Roekeghem, Tao Wang,
Georg KH Madsen, and Natalio Mingo. almabte: A solver of the space–time dependent
boltzmann transport equation for phonons in structured materials. Computer Physics
Communications, 220:351–362, 2017.

[46] Seung Hyun Lee, Kang Min Lee, Young-Bin Kim, Yoon-Jong Moon, Soo Bin Kim,
Dukkyu Bae, Tae Jung Kim, Young Dong Kim, Sun-Kyung Kim, and Sang Woon Lee.
Sub-microsecond response time deep-ultraviolet photodetectors using α-ga2o3 thin films
grown via low-temperature atomic layer deposition. Journal of Alloys and Compounds,
2018.

[47] Jean-Philippe M Péraud, Colin D Landon, and Nicolas G Hadjiconstantinou. Monte

carlo methods for solving the boltzmann transport equation.

[48] Patrick E Hopkins, Pamela M Norris, and Robert J Stevens.

Influence of inelastic
scattering at metal-dielectric interfaces. Journal of Heat Transfer, 130(2):022401, 2008.
[49] Minyoung Jeong, Justin P Freedman, Hongliang Joe Liang, Cheng-Ming Chow, Vin-
cent M Sokalski, James A Bain, and Jonathan A Malen. Enhancement of thermal
conductance at metal-dielectric interfaces using subnanometer metal adhesion layers.
Physical Review Applied, 5(1):014009, 2016.

[50] Brian F Donovan, Chester J Szwejkowski, John C Duda, Ramez Cheaito, John T Gask-
ins, C-Y Peter Yang, Costel Constantin, Reese E Jones, and Patrick E Hopkins. Thermal
boundary conductance across metal-gallium nitride interfaces from 80 to 450 k. Applied
Physics Letters, 105(20):203502, 2014.

[51] Tianli Feng, Lucas Lindsay, and Xiulin Ruan. Four-phonon scattering significantly
reduces intrinsic thermal conductivity of solids. Physical Review B, 96(16):161201, 2017.

116

[52] Krishnendu Ghosh and Uttam Singisetti. Calculation of electron impact ionization co-
In 72nd Device Research Conference. Institute of Electrical &

efficient in β-Ga2O3.
Electronics Engineers (IEEE), June 2014.

[53] Antonella Parisini and Roberto Fornari. Analysis of the scattering mechanisms con-
trolling electron mobility in β-Ga2O3 crystals. Semiconductor Science and Technology,
31(3):035023, 2016.

[54] Fabio Orlandi, Francesco Mezzadri, Gianluca Calestani, Francesco Boschi, and Roberto
Fornari. Thermal expansion coefficients of β-Ga2O3 single crystals. Applied Physics
Express, 8(11):111101, 2015.

[55] Wolfram Miller, Klaus Böttcher, Zbigniew Galazka, and Jürgen Schreuer. Numerical

modelling of the czochralski growth of ga2o3. Crystals, 7(1), 2017.

[56] Yuichi Oshima, Elaheh Ahmadi, Stefan C. Badescu, Feng Wu, and James S. Speck.

Applied Physics Express, 9(6):061102, 2016.

[57] Rostam Golesorkhtabar, Pasquale Pavone, Jürgen Spitaler, Peter Puschnig, and Claudia
Draxl. Elastic: A tool for calculating second-order elastic constants from first principles.
Computer Physics Communications, 184(8):1861 – 1873, 2013.

[58] Y Wang, J J Wang, H Zhang, V R Manga, S L Shang, L-Q Chen, and Z-K Liu. A
first-principles approach to finite temperature elastic constants. Journal of Physics:
Condensed Matter, 22(22):225404, 2010.

[59] Daoyou Guo, Zhenping Wu, Peigang Li, Yuehua An, Han Liu, Xuncai Guo, Hui Yan,
Guofeng Wang, Changlong Sun, Linghong Li, and Weihua Tang. Fabrication of β-ga2o3
thin films and solar-blind photodetectors by laser mbe technology. Opt. Mater. Express,
4(5):1067–1076, May 2014.

[60] Chris Sturm and Marius Grundmann. Singular optical axes in biaxial crystals and
analysis of their spectral dispersion effects in β − ga2o3. Phys. Rev. A, 93:053839, May
2016.

[61] Shunli Shang, Yi Wang, and Zi-Kui Liu. First-principles elastic constants of α- and

θ-al2o3. Applied Physics Letters, 90(10), 2007.

[62] Shun-Li Shang, Hui Zhang, Yi Wang, and Zi-Kui Liu. Temperature-dependent elastic
stiffness constants of α- and θ-al2o3 from first-principles calculations. Journal of Physics:
Condensed Matter, 22(37):375403, 2010.

[63] J. F. Nye. Physical properties of crystals : their representation by tensors and matrices.

Clarendon Press Oxford University Press, 1985.

[64] R. Yu, J. Zhu, and H.Q. Ye. Calculations of single-crystal elastic constants made simple.

Computer Physics Communications, 181(3):671 – 675, 2010.

[65] Marco D. Santia, Nandan Tandon, and J. D. Albrecht. Thermal properties of β-ga2o3

from first principles. MRS Advances, 1:109–114, 1 2016.

117

[66] Gene Simmons. Single Crystal Elastic Constants and Calculated Aggregate Properties:

A Handbook. M.I.T. Press, Cambridge, Mass, 1971.

[67] Per Söderlind and John E. Klepeis. First-principles elastic properties of α-Pu. Phys.

Rev. B, 79:104110, Mar 2009.

[68] Xue Jiang, Jijun Zhao, and Xin Jiang. Correlation between hardness and elastic moduli

of the covalent crystals. Computational Materials Science, 50(7):2287 – 2290, 2011.

[69] J. Feng, B. Xiao, R. Zhou, and W. Pan. Anisotropy in elasticity and thermal conduc-
tivity of monazite-type {REPO4} (re = la, ce, nd, sm, eu and gd) from first-principles
calculations. Acta Materialia, 61(19):7364 – 7383, 2013.

[70] K Adachi, H Ogi, N Takeuchi, N Nakamura, H Watanabe, T Ito, and Y Ozaki. Unusual
elasticity of monoclinic β- ga 2 o 3. Journal of Applied Physics, 124(8):085102, 2018.
[71] G.F. Davies. Effective elastic moduli under hydrostatic stress—i. quasi-harmonic theory.

Journal of Physics and Chemistry of Solids, 35(11):1513 – 1520, 1974.

[72] Yingchun Ding and Bing Xiao. Anisotropic elasticity, sound velocity and thermal con-
ductivity of tio2 polymorphs from first principles calculations. Computational Materials
Science, 82:202 – 218, 2014.

[73] Kristina E. Lipinska-Kalita, Patricia E. Kalita, Oliver A. Hemmers, and Thomas Hart-
mann. Equation of state of gallium oxide to 70 GPa: Comparison of quasihydrostatic
and nonhydrostatic compression. Phys. Rev. B, 77:094123, Mar 2008.

[74] VI Nikolaev, V Maslov, SI Stepanov, AI Pechnikov, V Krymov, IP Nikitina, LI Guzilova,
VE Bougrov, and AE Romanov. Growth and characterization of β-ga2o3 crystals.
Journal of Crystal Growth, 457:132–136, 2017.

[75] Orson L. Anderson. A simplified method for calculating the debye temperature from

elastic constants. Journal of Physics and Chemistry of Solids, 24(7):909 – 917, 1963.

[76] M Herms, M Fukuzawa, M Yamada, J Klöber, G Zychowitz, and JR Niklas. Photoe-
lastic characterization of residual strain in gaas wafers annealed in holders of different
geometry. Materials Science and Engineering: B, 66(1-3):7–10, 1999.

[77] Kelsey A. Mengle, Guangsha Shi, Dylan Bayerl, and Emmanouil Kioupakis. First-
principles calculations of the near-edge optical properties of β-ga2o3. Applied Physics
Letters, 109(21):212104, 2016.

[78] S Mendelson. Dislocations and plastic flow in nacl single crystals. i. Journal of Applied

Physics, 33(7):2175–2181, 1962.

[79] D Royer. Elastic Waves in Solids. Springer, Berlin New York, 2000.
[80] Nan Ma, Nicholas Tanen, Amit Verma, Zhi Guo, Tengfei Luo, Huili Xing, and Deb-
Intrinsic electron mobility limits in β-ga2o3. Applied Physics Letters,

deep Jena.
109(21):212101, 2016.

118