ELECTRONIC STRUCTURE AND THERMOELECTRIC PROPERTIES OF NARROW BAND
GAP SEMICONDUCTORS AND PSEUDO-GAP SYSTEMS
By
Dat Thanh Do

A DISSERTATION
Submitted
to Michigan State University
in partial fulfillment of the requirements
for the degree of
Condensed Matters Physics - Doctor Of Philosophy
2013

ABSTRACT
ELECTRONIC STRUCTURE AND THERMOELECTRIC PROPERTIES OF NARROW BAND
GAP SEMICONDUCTORS AND PSEUDO-GAP SYSTEMS
By
Dat Thanh Do
The direct energy conversion from heat to electricity without any moving part using the thermoelectric effects is attractive and has many promising applications in power generation and heat
pumping (refrigeration) devices. However, the wide use of thermoelectric materials in daily life is
still not practical due to its low efficiency. Presently, there has been renewed interest in thermoelectrics with new strategies for improving the efficiency, mostly by controlling the morphology
and dimensionality, and manipulating the electronic structure of novel complex systems. There are
new questions arise and several old questions have not yet been answered. My group at Michigan State University has been actively working on narrow band gap semiconductors and pseudogap systems including BiTe, BiSe, SbTe, AgPbm SbTem+2 (LAST), III-VI compounds, Mg2 Si,
Cu3 SbSe4 , Fe2 VAl, ZrNiSn, etc. , focusing on their thermoelectrics related properties. This dissertation tries to address some of the fundamental questions on the electronic structure and related
thermoelectric properties of narrow band gap semiconductors and pseudo-gap systems, employing the state-of-art density functional theory (DFT) and Boltzmanns transport equation. I focuses
on some well-known materials including Heusler compound, Fe2 VAl, tetrahedrally bonded compounds (Cu3 SbSe4 and related systems) and a newly reported novel nanocomposite system of
Half-Heusler–Heusler (ZrNiSn–ZrNi2 Sn). For the first two systems (Heusler and tetrahedrally coordinated systems), the effects of correlated electrons in the system containing d-electrons, valency
of sp-elements, and the role of lone-pair electrons in band gap formation and effects on thermoelectric properties are fully investigated. For the third system, Half-Heusler–Heusler nanocomposite, I
report a detailed energetics analysis of the nanostructure formation and its effects on the electronic
structures of the interested systems.

Copyright by
DAT THANH DO
2013

To my family

iv

ACKNOWLEDGEMENTS

I am extremely grateful to those who in anyway contribute to my work at Michigan State University
(MSU). First of all, I would like to express my sincere thanks to my supervisor Prof. S. D. Mahanti
for helping me become a theoretical physicist, for his patient, understanding and inspiration. Over
last five year, he has dedicated countless amount of time and effort to support me in my research
as well as in writing this dissertation.
I would like to thank Prof. Donald Morelli,Prof. Carlo Piermarocchi,Prof. Wayne Repko, and
Prof. Chong-Yu Ruan for serving on my guidance committee and for providing useful advices
when I needed.
I would also like to thank my collaborations, Prof. Donald Morelli’s Group (MSU), Prof.
Christopher Wolverton’s Group (Northwestern University – NWU), Prof. Vidvuds Ozoline’s (University of California at Los Angeles – UCLA), Prof. Mercouri G. Kanatzidis’s (NWU), especially
to those who I have been directly working with, Dr. Eric Skoug, Dr. Chang Liu (MSU), Dr. Yongsheng Zhang, Dr. Daniel Shoesmaker, Dr. Thomas Chasapi, Dr. Utpal Chatterjee (NWU), for their
valuable contributions.
My gratitude also goes to the past members of my group at MSU, Dr. Khang Hoang, Dr. Rak
Zsolt, Dr. Mal-Soon Lee, for stimulus and useful discussions.
Special thanks goes to my family, my friends for always being by my side and supporting me
when I needed the most.
I would like to acknowledge the Energy Frontier Research Center (EFRC) of at MSU, Revolutionary Materials for Solid State Energy Conversion, a research center of the U.S. Department
of Energy, for funding my research and providing an extremely active networking environment,
giving me chance to work with and learn from people from different disciplines.
Last but not least, I would like to thank Vietnam Education Foundation (VEF) for giving me
scholarship to study at MSU.

v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

LIST OF SYMBOLS AND ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . xiv
CHAPTER 1
INTRODUCTION . . . . . . . . . . . . . . .
1.1 Thermoelectrics . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Brief history of thermoelectrics . . . . . . . . . . .
1.1.2 Mechanisms of thermoelectric effects . . . . . . .
1.1.3 Thermoelectric applications . . . . . . . . . . . .
1.1.4 Figure of merit and device efficiency . . . . . . . .
1.2 Theoretical aid in pursuing better figure of merit (ZT ) . . .
1.3 Narrow band gap semiconductors and pseudo-gap systems
1.4 Outline of the dissertation . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

. 1
. 1
. 1
. 4
. 5
. 8
. 11
. 13
. 16

CHAPTER 2
THEORETICAL METHODOLOGY . . . . . . . . . . . . . . .
2.1 Basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Independent-electron – Hartree and Hartree-Fock approximation . . . . . . . .
2.2.1 Hartree approximation . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Hartree-Fock approximation (HFA) . . . . . . . . . . . . . . . . . . .
2.3 Exchange and correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Fundamentals of density functional theory (DFT) . . . . . . . . . . . . . . . .
2.5 Functionals for exchange and correlation . . . . . . . . . . . . . . . . . . . . .
2.5.1 Local (spin) density approximation (LSDA) and the gradient correction
(GGA) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.2 Local treatment of strongly correlated electron systems: LDA+U . . . .
2.5.3 Non-local treatment of exchange-correlation in correlated electron systems: Hybrid functional theory . . . . . . . . . . . . . . . . . . . . . .
2.6 Boltzmann’s equation for transport . . . . . . . . . . . . . . . . . . . . . . . .
2.7 Computational details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

. . 33
. . 34
. . 37
. . 39
. . 42

CHAPTER 3
3.1
3.2
3.3
3.4

HEUSLER COMPOUND – Fe2 VAl: EFFECTS OF ON-SITE
COULOMB REPULSION . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Computational details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Electronic structure: Pseudo-gap versus real gap . . . . . . . . . . . . . . . . .
Thermoelectric properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 GGA calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 GGA+U calculation : the effect of finite band gap . . . . . . . . . . . .
3.4.3 Comparison with experiment and other theoretical results . . . . . . . .

17
17
20
20
22
25
26
30

vi

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

44
44
48
49
56
56
59
60

3.5

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

CHAPTER 4
DIAMOND-LIKE TERNARY COMPOUNDS . . . . . . . . . . .
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Crystal structures of Famatinite and Enargite . . . . . . . . . . . . . . . . . . . .
4.3 Electronic structure of Cu3 SbSe4 : band gap formation and the role of Sb lone
pairs and non-local exchange effect . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 GGA calculation and indication of the existence of Sb lone pair electrons
4.3.2 Failure of local exchange approximation – Effect of non-local exchange
interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.3 Unphyscal value of U and its meaning . . . . . . . . . . . . . . . . . . .
4.4 Bond and band gap in the a class of I3 -V-VI4 Famatinite and Enargite compounds
4.4.1 Atomic energy levels of constituents and some schematic studies. . . . .
4.4.2 Results obtained using local approximations . . . . . . . . . . . . . . . .
4.4.3 Effect of non-local exchange . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Thermopower calculation of Cu3 SbSe4 . . . . . . . . . . . . . . . . . . . . . . .
4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

77
81
83
84
85
92
92
96

CHAPTER 5
HALF-HEUSLER–HEUSLER COMPOSITE . . . . . . .
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Methods of electronic structure calculation and modeling nanostructures
5.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 Formation energy and energetics of clustering . . . . . . . . . .
5.3.1.1 HH with excess Ni . . . . . . . . . . . . . . . . . . .
5.3.1.2 FH with deficient Ni . . . . . . . . . . . . . . . . . .
5.3.2 Electronic structure . . . . . . . . . . . . . . . . . . . . . . . .
5.4 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

98
98
103
105
105
105
108
112
117

CHAPTER 6

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

. 66
. 66
. 69
. 71
. 74

CONCLUSION AND OUTLOOK . . . . . . . . . . . . . . . . . . . 119

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

vii

LIST OF TABLES

Table 2.1: Summary of approximations used in DFT . . . . . . . . . . . . . . . . . . . . . 32
Table 3.1: Summary of experimental studies of energy gap in Fe2 VAl . . . . . . . . . . . . 44
Table 3.2: Summary of previous theoretical studies on Fe2 VAl electronic structure . . . . . 45
Table 3.3: Summary of Fe2 VAl band gap using different methods . . . . . . . . . . . . . . 54
Table 4.1: Summary of crystal structures of Cu3 AB4 . . . . . . . . . . . . . . . . . . . . . 71
Table 4.2: Nearest neighbor distances (Å) and band gaps (eV) with different exact exchange ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
Table 4.3: Atomic energies of the s and p valence electrons for different constituents and
for the monovalent Cu, the energy of Cu d-level is included . . . . . . . . . . . 84
Table 4.4: Effect of different V element substitutions on the band gap of Cu3 SbSe4 with
and without ionic relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
Table 5.1: Formation energies (eV) of HH-FH systems (for the lowest energy configurations)105
Table 5.2: Ground state energy as a function of Ni3 ’s relative position with respect to
Ni1 -Ni2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

viii

LIST OF FIGURES

Figure 1.1: Schematic diagram of a simple thermo-couple made of two conductors A and
B with junction 1 and 2 kept in a heat sink and heat source respectively. (For
interpretation of the references to color in this and all other figures, the reader
is referred to the electronic version of this dissertation.) . . . . . . . . . . . . .

2

Figure 1.2: Schematics explain (a) the mechanism of Seebeck, temperature gradient drives
the movement of carriers, creating a current, and (b) Peltier effects, the current of carriers transfer the heat from cool side to the hot side. . . . . . . . . . .

5

Figure 1.3: Power generation efficiency versus temperature of the hot side, plotted for
different energy conversion technologies. The cold side is assumed to be
at room temperature. The efficiency range of some of the other renewable
energy technologies is marked in a bar on the right-hand side of the graph.
In this graph PV is photovoltaic; CSP denotes concentrated solar power; SJ
and MJ denotes single- and multiple-junction; and Org, TE, and TPP denote organic, thermoelectric and thermionic devices, and thermal power plan
respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

Figure 1.4: ZT of several TE materials and year of their discovery, showing how ZT or
materials have been improved. . . . . . . . . . . . . . . . . . . . . . . . . . .

8

Figure 1.5: Competing relation between S, σ , and κ . . . . . . . . . . . . . . . . . . . . . 10
Figure 1.6: Schematic diagram of the density of state of a material with a delta function
defect state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Figure 1.7: Band gap of several typical semiconductors . . . . . . . . . . . . . . . . . . . 14
Figure 1.8: Diagram explaining different scheme of band gap in a semiconductor. . . . . . 15
Figure 2.1: The factor f (x) in Hartree-Fock approximation for homogeneous electron gas. . 24
Figure 2.2: The ladder of density functional theory. Notations on the left side represent
the terms added to the theory, i.e. only density n functional is considered
in LSDA, GGA adds gradient correction ∇n, meta-GGA adds kinetic energy
density τ , hybrid GGA takes into account parts of exact non-local exchange
exact , RPA includes unoccupied orbital. . . . . . . . . . . . . . . . . . 31
energy Ex
Figure 3.1: Band structure of Fe2 VAl showing d-characters of Fe and V (left); and DOS
(right). EF denotes the Fermi level set to be zero. The size of the symbols
represents the strength of orbital characters. . . . . . . . . . . . . . . . . . . . 49

ix

Figure 3.2: Band structure of Fe2 VAl in absolute scale for UFe of (a) 0 eV, (b) 1 eV, (c)
2 eV and (d) 4 eV. The Fe-eg and Fe-t2g bands are pointed out by arrows. . . . . 51
Figure 3.3: Band structure of Fe2 VAl in absolute scale for UV of (a) 0 eV, (b) 1 eV, (c)
2 eV and (d) 4 eV.The V-eg bands are pointed out by an arrow. . . . . . . . . . 53
Figure 3.4: Band structure of Fe2 VAl obtained in different methods: (a) GGA-PBE , (b)
GGA+U , (c) mBJ and (d) PBE0. . . . . . . . . . . . . . . . . . . . . . . . . . 54
Figure 3.5: Thermopower (S) and chemical potential (µ ) versus carrier concentration
(ne , nh ) in the absence of U (the pseudo-gap case) for n- and p-type doping
in Fe2 VAl. The Fermi level is chosen to be 0. . . . . . . . . . . . . . . . . . . 56
Figure 3.6: S as functions of temperature obtained using the GGA band structure (U=0 eV,
the pseudo-gap case) at different concentrations: for n-type (closed symbols)
and p-type (open symbols) doping. . . . . . . . . . . . . . . . . . . . . . . . . 58
Figure 3.7: S as functions of T for GGA+U calculation (UF e=4 eV and UV =1 eV) at different concentration: for n-type (closed symbols) and p-type (open symbols)
doping. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
Figure 3.8: Calculated p-type S as a function of n for the cases of: pseudo-gap (GGA calculation) (solid line) and gapped (GGA+U calculation), Eg =0.55 eV, (dashed
line) at 300 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
Figure 3.9: Experimental thermopower obtained by doping Si on Al site (symbols) and
theoretical curve (lines) with U=1 eV (Eg = 0.07 eV ). (carrier concentrations
are given in unit of cm−3 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
Figure 4.1: Crystal structure of Famatinite (Fa) and Enargite (En) structure are represented as double unit cell of the Zinc blend (top panel) and the Wurtzite
(bottom panel). Green represents S (Se, Te), gray represents Zn which are
replaced by Copper/Silver (blue) and red represents Zn which are replaced
by Sb (P, As, Bi). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
Figure 4.2: Crystal structure of Cu3 SbSe4 (left) and Cu3 SbSe3 (right), Cu is blue ball,
Sb is brown, and Se is green. . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
Figure 4.3: Band structures of Cu3 SbSe4 (a) and Cu3 SbSe3 (b) using different methods.
In (a), the bold (red) line is “BOI”, the doted lines present the HSE band
structure using GGA position. UCu−d = 10 eV was used in GGA+U calculation. The energy is given in absolute scale and the Fermi energy is denoted
by EF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

x

Figure 4.4: (a) Total DOS of Se4 and Se3 (top panel) and the projected DOS of Se4
(bottom panel), along with isosurface plot of charge density associated with
(b) the lowest band at -14.5 eV, (c) the Sb-s band at -9.5 eV, and (d) the
“BOI”. For clarity, we only show SbSe4 tetrahedron with Sb at the center.
The DOS are given in unit of number of states per eV per formula unit. . . . . . 75
Figure 4.5: Simple bonding-antibonding scheme in Cu3 SbSe4 showing energy levels in
the atomic limit (left column) intermediate bonding-antibonding states (middle column) and bonding-antibonding states in crystal (right column). (Cu d
and Sb-p are not shown for simplicity). The number of states associated with
each level is given in the bracket. . . . . . . . . . . . . . . . . . . . . . . . . . 77
Figure 4.6: Band structure of Cu3 SbSe4 in absolute energy scale for the ratios of exact
exchange α = (a) 0, (b) 0.1, (c) 0.2 and (d) 0.3 with fat band representation
of orbital character, showing that the bands with dominating d characters
affected strongly; there is an exchange of band character between BOI and
BOI-3 at a critical value of α , αc ≈ 0.2. . . . . . . . . . . . . . . . . . . . . . 78
Figure 4.7: (a) Effect of U on band structure and (b) comparison between GGA+U and
HSE06 band structures. The energy is given in the absolute scale. . . . . . . . 82
Figure 4.8: Visualization of the atomic levels of constituent atoms. . . . . . . . . . . . . . 84
Figure 4.9: GGA bandstructure of (a) Cu3 PS4 (b) Cu3 PSe4 (c) Cu3 AsS4 (d) Cu3 AsSe4
(e) Cu3 SbS4 (f) Cu3 SbSe4 (g) Cu3 BiS4 and (h) Cu3 BiSe4 . Where the first
two compounds are En and the rest compounds are Fa. . . . . . . . . . . . . . . 86
Figure 4.10: GGA+U bandstructure of (a) Cu3 PS4 (b) Cu3 PSe4 (c) Cu3 AsS4 (d) Cu3 AsSe4
(e) Cu3 SbS4 (f) Cu3 SbSe4 (g) Cu3 BiS4 and (h) Cu3 BiSe4 . Where the first
two compounds are En and the rest compounds are Fa. . . . . . . . . . . . . . . 87
Figure 4.11: Relation between band gap and the distance between V and VI elements. . . . . 89
Figure 4.12: Band structure of Na3SbSe4 (Fa structure) obtained using (a) GGA and (b)
HSE06 with fatband representation showing atomic orbital associated with
energy levels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
Figure 4.13: HSE06 bandstructure of (a) Cu3 PS4 (b) Cu3 PSe4 (c) Cu3 AsS4 (d) Cu3 AsSe4
(e) Cu3 SbS4 (f) Cu3 SbSe4 (g) Cu3 BiS4 and (h) Cu3 BiSe4 . Where the first
two compounds are En and the rest compounds are Fa. . . . . . . . . . . . . . . 93
Figure 4.14: Band gap as a functions of V-VI distance for Fa and En: a comparison between GGA+U and HSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

xi

Figure 4.15: Theoretical seebeck coefficient (negative means n-type, positive means ptype) of Se4 using Boltzmann’s equation for different carrier concentration
(1018 cm−3 ) together with experimental data in the similar range of doping
levels. x is the excessed holes per unit cell, x=0.01 corresponds to hole concentration of 50 × 1018 cm−3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
Figure 5.1: Improvement of powerfactor (PF) (a), and TEM image showing FH-nanoinclusions in HH(b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
Figure 5.2: (a) Phase diagram of HH-FH mixture ZrNi1+x Sn, (b) Comparison between
measured susceptibility and density of state calculated using Korringa-KohnRostocker method with coherent potential approximation and local density
functional (KKR-CPALDA) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Figure 5.3: Crystal structure of (a) Half-Heusler (HH) and (b) Full-Heusler (FH) compounds, where Zr (orange circles) and Sn atoms (gray circles) form a NaCl
sub-lattice, inside which there are 8 small cubes. HH is formed by filling
every other cube with Ni atoms (blue circles), while FH is formed by filling
all the cubes. The empty-cube sites (quasi-particle of Ni-vacancy) in HH is
presented by a dashed-line circle. . . . . . . . . . . . . . . . . . . . . . . . . . 101
Figure 5.4: Left column presents energetics of two extra Ni in HH: energy vs. distance.
Right column presents the preferable configuration of the excess Ni (orange
circle) (only Ni matrix (blue circles) is shown). . . . . . . . . . . . . . . . . . . 106
Figure 5.5: Left column presents energetics of three extra Ni in HH: position of the third
Ni (Ni3 ) is given in polar coordinate, with respect to Ni1 -Ni2 , whose positions were previously determined by studying the energetics of a pair of excess Ni (FIG. 5.4). Ground state energy of each configuration is represented
by color code from low (purple) to high (yellow) energy. Right column shows
the preferred configuration of the Ni3 (only Ni matrix (blue circle) is shown).
The excess Ni are orange circles. . . . . . . . . . . . . . . . . . . . . . . . . . 107
Figure 5.6: Conventional notations for Ni-sites: “f” for Ni at (1/4,1/4,1/4) and equivalent
atoms (dark blue), “h” for Ni at (3/4,3/4,3/4) and equivalent atoms (bright yellow).108
Figure 5.7: The top figure shows energetics of FH with two Ni-vacant sites (white circles); the bottom shows structures of the three most preferable configurations,
which have ground state energies close to each other: (a) f-f1, (b) f-f2, and
(c) f-h4. Only Ni matrix (blue circles) is shown. . . . . . . . . . . . . . . . . . 109
Figure 5.8: Band structures of (a) HH and (b) FH in the Brillouin zone (BZ) of 2 × 2 ×
2 cubic supercell; and (c) the density of states (DOS) within energy range
[-2,2]; in the onset, we zoom in the range [-1,1] to show the fundamental
difference between HH and FH: the former is metal whereas the latter is
semiconductor, with a band gap ∼ 0.5 eV. . . . . . . . . . . . . . . . . . . . . 111
xii

Figure 5.9: Charge difference between (a) HH+1Ni and (b) FH-1Ni and their parent compounds (HH and HF), where blue color indicates charge depletion whereas
red color indicates charge accumulation. The figure show the contour plot on
(-110) plane. Charge difference is defined as ρ X (r) − ρ 0 (r), where ρ X (r),
and ρ 0 (r) are charge density of defective and pure system respectively. Host
Ni are blue circles, excess Ni is red circle, Zr are brown circles, Sn are gray
circles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
Figure 5.10: Evolution of electronic structure of HH with excess Ni: band structures (a)
HH+1Ni, (b) HH+2Ni, and (c) HH+3Ni, and (d) density of states of HH with
1, 2, 3 and 4 additional Ni. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
Figure 5.11: The isosurface of the charge density associated with the in-gap defect states
in HH+1Ni. Host Ni are blue circles, excess Ni is the red circle, Zr are brown
circles, Sn are gray circles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
Figure 5.12: Evolution of electronic structure of FH with deficient Ni: band structures (a)
FH, (b) FH-1Ni, and (c) FH-2Ni , and density of states of those systems. The
inset show the DOS, indicated by the arrow, around -1 eV in small range
[-1.2,-0.9]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

xiii

LIST OF SYMBOLS AND ABBREVIATIONS

α

Ratio of exact Hartree-Fock exchange in hybrid functional HSE06

ε0

Electrical susceptibility

εi

Eigenvalue of a single-electron wavefunction

ηc

Carnot efficiency

h
¯

Planc’s constant

κ

Thermal conductivity

κi

Thermal conductivity of element i

κelec Thermal conductivity due to electrons
κlatt

Thermal conductivity due to lattice vibration

ω

Screening parameter in hybrid functional HSE06

Ψ

Many-body wavefunction

ψi

Single-electron wavefunction

σ

Electrical conductivity

τk

Relaxation time of carrier with momentum k

σ ,σi

Spin index

a0

The first Bohr radius, a0 =

Exc

Exchange-correlation energy

4πε0 h2
¯
2 ≈ 0.529Å
me e

1
f (ε ) Fermi-Dirac distribution, f (ε ) = β (ε −µ )
+1
e

U

On-site Coulomb repulsion parameter in LDA+U

COP Coefficient of Performance, percent (%)
DFT Density functional theory
FH

(Full) Heusler compounds, in this thesis it is usually referred to ZrNi2 Sn

GGA Generalized gradient approximation
HFA Hartree-Fock approximation

xiv

HH

Half-Heusler compounds, in this thesis it is usually referred to ZrNiSn

HSE06 Hybrid functional with screened exact exchange potential proposed by Heyd, Scuseria and
Ernzerhof
IEA

International Energy Agency

KS

Kohn-Sham

L(S)DA Local (spin) density approximation
LDA+U, GGA+U LDA, GGA with on-site Coulomb repulsion correction
NBGS Narrow-band-gap semiconductor
RPA

Random phase approximation

TE

Thermoelectric(s)

TEG Thermoelectric generator
TI

Topological insulator

E

Electrical field

E

Total energy of a system

e

Electron’s charge

ˆ
H

Hamiltonian operator

kB

Boltzmann’s constant, kB = 8.6173324 × 105 eV/K

me

Electron’s mass

N

Number of electrons

n
ˆ

n(r) = ∑i=1,N δ (r − ri ) is the particle density operator
ˆ

ˆ
O

Operator corresponding to an observable

S

Seebeck coefficient or thermopower

Si

Thermopower of element i

T

Operating temperature (T = (Th + Tc )/2)

Tc

Temperature of heat sink

Th

Temperature of heat source

ˆ
T

Kinetic energy operator

ˆ
Vext

External potential operator
xv

ˆ
Vint

Electron-electron interaction potential operator

Z

Figure of merit

ZT

Dimensionless figure of merit, where T is temperature

xvi

CHAPTER 1
INTRODUCTION

1.1 Thermoelectrics
In this era of technology, energy is a key factor in any economy. In fact, our world runs on energy.
According to International Energy Agency (IEA) "World energy outlook 2012"[1], 81% of total
primary energy consumption is fossil fuel preserved in the form of coal, petroleum, and natural
gas. The energy consumption has increased about 35% from 1990 to 2010 and the same rate is
expected for the next tens years. This brings the world’s energy demand up to about 15 billions
tonne of oil equivalent (toe) (or about six hundred millions Watt-second) per year by 2035. The
reserve sources of fossil fuel are limited and will run out in the future. There is therefore an urgent
need to use energy more efficiently and to find other sustainable forms of energy such as solar,
wind power, bio fuel cell, etc. Among the ongoing directions, thermoelectric converters, with high
enough efficiency, can serve well for both purposes.

1.1.1

Brief history of thermoelectrics

Thermoelectric converter directly converts heat gradient into electrical current or vice versa. The
mechanism behind such devices is the thermoelectric effect which is usually referred to as SeebeckPeltier effect to honor the two physicist who discovered it almost two centuries ago.
In 1821, the Baltic German physicist Thomas Johann Seebeck discovered that a circuit of two
dissimilar metals could deflect the compass needle when exposed to a temperature gradient.[2]
Seebeck, however, first thought it was "Magnetism induced by the temperature difference". It
was quickly realized that the temperature difference produces a potential (voltage V) which drives
an electrical current through a closed circuit, and that current deflects the compass needle. It
was Hans Christian Ørsted, a Danish physicist, who first used the term "thermoelectricity". The

1

SINK

τA I∇T
A

.

˙
Q2 = Π2 I

.

˙
Q1 = Π1 I

SOURCE

1

2

I
B
.

T0

τB I∇T

T0 + ∆T

Figure 1.1: Schematic diagram of a simple thermo-couple made of two conductors A and B with
junction 1 and 2 kept in a heat sink and heat source respectively. (For interpretation of the
references to color in this and all other figures, the reader is referred to the electronic version of
this dissertation.)

potential difference generated was found to be independent of the distribution of temperature but
proportional to the temperature difference (∆T ) between two junctions[3]
VAB = SAB ∆T,

(1.1)

where the proportional coefficient SAB is called thermopower or Seebeck coefficient. VAB and SAB
are positive if the current flows from A to B at the hot junction. This is the mechanism of a simple
thermo-couple as depicted in figure 1.1.
About a decade after Seebeck’s discovery, in 1834, the French physicist Jean Charles Athanase
Peltier discovered the reverse process when observing that an electrical current could heat or cool
the junction of two dissimilar metals.[4] In 1838, Lenz showed that depending on the direction of
the electrical current, heat could be removed or added to the junction.[5] The heat being transferred
in a unit of time is proportional to the current
˙
QAB = ΠAB I,

(1.2)

˙
where ΠAB is called the Peltier coefficient. Similar to the equation of Seebeck effect, QAB and
ΠAB are positive if heat is absorbed when the current flows from A to B.

2

Seebeck-Peltier effects are additive, meaning if there are two junction AB and BC the effective
proportionality coefficients of AC is the summation of the former
VAC = VAB +VBC = (SAB + SBC )∆T,

(1.3)

˙
˙
˙
QAC = QAB + QBC = (ΠAB + ΠBC )I.

(1.4)

This additive property implies that the the Seebeck(Peltier) coefficient of any pair AB is the difference between absolute Seebeck(Peltier) coefficients, S(Π),
SAB = SA − SB and ΠAB = ΠA − ΠB .

(1.5)

These absolute quantities were formulated later by William Thomson (also known as Lord Kelvin).[4]
Thomson studied the Seebeck-Peltier effects and discovered a new effect called Thomson effect
˙
QA = τA I∇T

(1.6)

predicting that the rate of heat absorption per unit length of a conductor is proportional to the electrical current I and temperature gradient ∇T . The proportional coefficient τ is called Thomson
coefficient. He then derived the relation between the three thermoelectric coefficients by applying
the first and second laws of thermodynamics on the thermoelectric circuit. The first law of thermodynamics requires that the work done by the Seebeck potential must equal to the thermal energy
absorbed from the surounding
V = VAB = VB −VA = Π1 − Π2 +

∫ T +∆T
0
T0

(τB − τA )dT

(1.7)

or
dV = dΠ + (τB − τA )dT.

(1.8)

This equation is called the first Thomson (or Kelvin) relation. The second law of thermodynamics
requires that the total change in entropy must be zero which give us
∫

0=

T0 +∆T τ − τ
Π1
Π2
B
A
−
+
dT.
T0 T0 + ∆T
T
T0

3

(1.9)

From (1.8) and (1.9) one gets
dV
dΠ
=
+ τB − τA
dT
dT
dΠ Π
0=
− + τB − τA
dT T

S=

(1.10)
(1.11)

and obtains the second Thomson relation
Π = ST

(1.12)

Differentiating (1.12) with respect to T and using (1.10), one gets
dS
τ − τB
= A
dT
T

(1.13)

which give the expression for absolute Seebeck coefficient of A at given T
SA =

∫ T
τA
0

T

dT

(1.14)

and the absolute Peltier coefficient of A
ΠA = T SA .

1.1.2

(1.15)

Mechanisms of thermoelectric effects

Seebeck-Peltier effect involves two physical phenomena: (i) charge-carrier diffusion and (ii) phonon
drag. The former occurs when a conductor exposed to a temperature gradient. There is unbalance
in hot and cold carrier concentration at the two ends of a conductor, the hot carriers tend to diffuse
from the hot end to the cooler one to reach a steady state. The charge-carrier diffusion builds up a
net charge at the cold end, giving rise to an electrostatic potential counteracting the chemical potential of the diffusion and brings the system to equilibrium. In the beginning, thermocouple was made
of metals, however, more practical approach is using alternative n-type and p-type semiconductors
(figure 1.2) connected by metals. In these devices, charges flow through the n-type element, cross
the metallic connector and pass into the p-type one, creating a current. The diffusion of charges
(and so the mobility and the effective mass of carriers) is affected by scattering. There are two
4

(b)

(a)

Heat source
.
J
+
−
− −
+
+ ..
.. .
p
n − − . +. .+
+
−−
+
.
.

Cooled side
.
J +
− −
+ +.
.− . .
. . +p
.
.−
n−
+
−
+ +
− .
.
.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.
.

.

.

.

.

.

.
.

.

.

.

.

.

.

.
.

.

.

.

.

.

.

.

.

.
.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

Disipated Heat

Heat sink

Figure 1.2: Schematics explain (a) the mechanism of Seebeck, temperature gradient drives the
movement of carriers, creating a current, and (b) Peltier effects, the current of carriers transfer the
heat from cool side to the hot side.
categories of scattering,[6] namely lattice vibrations or phonon and lattice imperfection, hence, the
thermopower of a material depends greatly on impurities, imperfections and structural changes.
In addition to charge-carrier diffusion, the latter effect, phonon drag, is found to produce large
thermoelectric effects in some semiconductors at low temperature (but not very low temperature).[6]
The physics is that in pure system, at low temperature, electron-phonon scattering is predominant.
Electrons have a very long relaxation time, and transfer their momentum to phonons in the direction of the electric current. This energy transport contributes to the Peltier heat which may be much
larger than the former effect. This happens for
1
T ≈ θD ,
5

(1.16)

where θD is the Debye temperature. At higher temperature the phonon-phonon scattering suppresses the phonon-drag effect.

1.1.3

Thermoelectric applications

About 90 per cent of the world’s energy is generated by burning fossil fuel, and in the best scenario
the ratio of fossil fuel is still about 80 percent by 2035[1]. The efficiency of typical heat engines is
about 30 – 40 percent, leaving more than half of energy lost in form of waste heat. Thermoelectric
5

devices can convert part of that waste into electricity and greatly improve the operating efficiency
of the engines. Moreover, the ability of operating without mechanically moving parts and with
very low maintaining requirement makes thermoelectric device greatly promising for a wide range
applications. They can fall into two categories: (i) power generation (using Seebeck effect), and
(ii) heat pumping or refrigeration (using Peltier effect).
The first commercial thermoelectric generator (TEG), made of ZnSb and constantan in USSR,
was dated back to 1948, produced electricity strong enough to power a radio.[7] Since then, the
field has been moving much further. Its applications ranges from big space missions[8–10] to
micropower harvesting.[11]
The thermoelectric module called Radioisotope Thermoelectric Generator (RTG) has played
a key role in many missions of the United States National Aeronautics and Space Administration
(NASA). The first use of RTGs was in the spacecraft "Transit 4A" launched in 1961 and then
in many other missions after. Recently an eighth RTG configuration, called the Multi-Mission
Radioisotope Thermoelectric Generator (MMRTG), has recently been used on the Mars Science
Laboratory rover, Curiosity. RTGs use radioactive isotope plutonium-238 decay as the heat source
and are preferable over other types of power generators because of its (i) long life time which
suitable a long mission in the deep space, (ii) ability to operate in an extreme environment such
as the high-radiation belts of Jupiter, extreme temperature on moon, or dust storm on Mars, (iii)
operational independence which guarantees to providing needed power regardless of the spacecraft
position or orientation and (iv) its reliability which has been proven to run for decades.
The most active area of thermoelectric application is in the automotive industry. Both major automotive companies have thermoeletric development programs.[12] It is estimated that in
a vehicle using internal combustion engine, only 25 percent of provided energy is used to move
it, 75 percent is just wasted, in which about 40 percent can be used for a thermoelectric power
generator.[12] TEG can use the wasted heat to recharge the battery. This can help to reduce or
even eliminate the need for alternator, and as a result reduce the load on the engine, therefore,
improve fuel efficiency by the order of 10 percent.[11] Similar applications can be used in any

6

Figure 1.3: (Source: Zebarjadi et al. [14]) Power generation efficiency versus temperature of the
hot side, plotted for different energy conversion technologies. The cold side is assumed to be at
room temperature. The efficiency range of some of the other renewable energy technologies is
marked in a bar on the right-hand side of the graph. In this graph PV is photovoltaic; CSP denotes
concentrated solar power; SJ and MJ denotes single- and multiple-junction; and Org, TE, and
TPP denote organic, thermoelectric and thermionic devices, and thermal power plan respectively.

thermal power plan to recover wasted heat.
Along with the development of science and technology, small and portable applications need
power sources smaller and lighter than conventional batteries. A TEG using combustible fuel could
be a promising substitution. There is, however, one big challenge that it needs at least 500◦ C in
order to reach desired efficiency.[11] Going to much smaller in size, TEG proves its strength since
other types are not scalable or at lease have much lower efficiency.[13] At small scale, TEG can
havest very small amount of heat for low power applications, such as wireless sensor networks,
mobile devices or medical applications.[11]
Despite the low efficiency in current thermoelectrics (< 10%), thanks to its unique properties
such as having no mechanical moving parts, fluid or refrigerants, flexible shape (they can be very
small), long life time and controllable, TE cooler and heater (heat pump) become common in daily

7

3

Bulk
Thin film

Bi2Te3/Sb2Te3-SL
LAST

2

PbTe-QDSL

ZT

Pb0.98Tl0.02Te
RECoSb3
Zn4Sb3
1
Bi2Te3

n-SiGe

CsBi2Te3

PbT
0
1950

1960

1970

1980
Year

1990

2000

2010

Figure 1.4: (Source: Tritt [15]) ZT of several TE materials and year of their discovery, showing
how ZT or materials have been improved.

life.[14] Such applications include portable mini cooler, refrigerator, coffee warmer. Thermoelectric heat pumps are also found in temperature/weather control devices such as cooled seat, weather
vest, temperature regulators of semiconductor lasers, and medical.[14] However, for large scale applications involving higher powers, one need to develop (or discover) thermoelectrics with higher
efficiency (∼ 20%). This is the driving force behind current research in novel thermoelectrics.

1.1.4

Figure of merit and device efficiency

The maximum efficiency that any thermal process can have is the Carnot efficiency ηc of an adiabatic circle. The efficiency of a TE device, usually referred as coefficient of performance (COP),
is given by
COP = ηc × M,

8

(1.17)

where M is the merit factor which is positive and smaller than unity. COPs are different for the
refrigerator and the power generator. This thesis is concerned mainly with the latter for which the
COP is given by

√
∆T
1 + ZT − 1
COPG = ηc × M =
×√
,
Th
1 + ZT + Tc /Th

(1.18)

where Th , Tc , and ∆T are temperatures of the hot (heat source) and the cold (heat sink) ends, and
the difference between them respectively, T = 1/2(Th + Tc ) is the operating temperature. The
figure of merit Z is given by
Z=

1 S2 σ
,
T κ

(1.19)

in which σ and κ are the electrical and thermal conductivity respectively and S is thermopower
(Seebeck coefficient). One often refers to ZT as the dimensionless figure of merit. Note that ZT
depends only on properties of materials such as S, σ , and κ . The larger the ZT , the higher is
the efficiency. The relation between efficiency and temperature for different ZT is given in figure
1.3. For comparison, the efficiency of other energy sources and their temperature range is also
given in the same figure. In automotive industry, the state of thermoelectric is currently at the first
generation of thermoelectric module with ZT =1 and the second generation with ZT =2 is under
progress.[16] To make thermoelectric applications economically viable, we need ZT of ∼3. Even
though, there is, in principle, no upper limit of ZT , obtaining high ZT , however, is not an easy job.
The history of thermoelectrics reveals that for almost a century after the Seebeck’s discovery of
thermoelectric effect, ZT barely reached unity (figure 1.4). It is during the last few decades that
the thermoelectric research has come with some break through ideas and new directions which can
improve ZT up to ∼ 2.5. According to equation 1.19, a good thermoelectric material needs (i)
large Seebeck coefficient, (ii) high electrical conductivity, and (iii) low thermal conductivity. The
difficulty in designing high ZT materials is the conflicting requirements of properties. An example
of this competing relation is given for a homogeneous conductor as followings.
(i) In order to get large Seebeck coefficient, which, for homogeneous conductor, is given by[17]
2
8π 2 kB ∗ ( π )2/3
m T
,
S=
3n
2eh2

9

(1.20)

1
S

σ

κelec

ZT

ZT

PF=S2σ

0
1018

1019
1020
Carrier concentrations (cm-3)

1021

Figure 1.5: Competing relation between S, σ , and κ
where kB and h are the Boltzmann, Planck constants respectively, one needs a low carrier
concentration n with large effective mass m∗ , or low carrier mobility µ .
(ii) On the other hand, a high electrical conductivity requires large n and highµ since

σ = neµ ,

(1.21)

where e is the charge of the carriers.
(iii) The electronic part of the thermal conductivity can be related to electrical conductivity by
Wiedemann-Franz law

κ = κlatt + κelec
= κlatt + LT σ = κlatt + LT neµ ,

10

(1.22)
(1.23)

ε
.

g(ε )

εF
δ (ε )

. .
Host

Figure 1.6: Schematic diagram of the density of state of a material with a delta function defect
state.
where κlatt and κelec are thermal conductivities due to lattice vibrations (phonons) and electrons, and L is Lorentz number. It means in order to get low thermal conductivity one needs
small n and µ .
This conflicting relation is demonstrated in figure 1.5. Hence, a large ZT could be only found
with optimized value of S, σ , and κ , which is a challenging problem. It turns out that semiconductors with small band gaps (∼0.2–0.5 eV) or pseudo-gap systems are the best candidates for
thermoelectric applications at room temperature.[18] For a brief discussion of semiconductors and
pseudo-gap systems see section 1.3.

1.2 Theoretical aid in pursuing better figure of merit (ZT )
Extensive works have been done in trying to maximize the magnitude of ZT . Approaches to improve the performance of TE materials can be categorized into two types: (i) engineering electronic
structures to optimize S,σ and κelec , and (ii) engineering phonons to minimize κlatt . Following
the first approach (engineering electronic structures), Mahan and Sofo [17] on their quest to “The

11

best thermoelectric” predicted that the enhancement of S can be achieved in a system where exist
a rapid change in the density of states g(ε ) near the Fermi level. This prediction can be easily
understood using the Mott’s formula for thermopower which is given by

π 2 kB
S=
k T
3 q B
π 2 kB
=
k T
3 q B

{
{

}
d[ln(σ (ε ))]
dε
ε =εF

}
1 dn(ε ) 1 d µ (ε )
,
+
n dε
µ dε
ε =εF

(1.24)

where n(ε ) = g(ε ) f (ε ) is the carrier density, f (ε ) = 1/(Exp[(ε − µ )/kB T ] is the Fermi function,

µ (ε ) is carrier’s mobility at energy ε , and εF is the Fermi level. The Mott’s equation predicts
that large S can be obtained when DOS (the first term in the bracket of Eqn. 1.24) or carrier’s
mobility (the second term in the bracket of Eqn. 1.24) varies rapidly at the Fermi level. One
approach on this direction is to find a material with a rapidly-varying region in DOS which is
accessible by doping. An example for this approach is the Heusler compoud Fe2 VAl[19], which
has a steep DOS near the band edge. Chapter 3 will be dedicated for discussing the electronic and
thermoelectric properties of this system. Another approach to achieve a large change of DOS in
a small range of energy is to introduce a delta (resonant) states into the electronics structure of
the host materials (Fig. 1.6). This approach has been quite successful in improvement of ZT . For
example, in the classic thermoelectric material PbTe, theoretical calculations by Ahmad et al. [20]
(2006) predicted that doping of group III elements (Ga, In, Tl) creates resonant levels. Heremans
et al. [21] later reported a large enhancement of ZT in Tl-doped PbTe whose ZT = 1.5 was found
at 700 K compared to ZT = 0.7 for the host material PbTe.
Another effect which possibly improves thermoelectric performance and can be achieved by
engineering electronic structure is energy filtering. Detailed discussions of this effect can be found
in reference [22] and [14], and references therein. This effect can enhance the performance of thermoelectric materials by blocking low-energy electron, and thus increase the average heat transfer
by an electron. This effect can be achieved by nanostructure inclusion in the host materials. An
exciting example of this approach is the work on Half-Heusler–(Full) Heusler nanocomposites initiated by Makongo et al. [23] who reported a significant improvement of the powerfactor by ∼25%.
12

I will discuss more about this system in chapter 5 of this thesis.
Engineering phonons approach to lower the thermal conductivity of thermoelectric materials
was considered as early as 1974 by Rowe [24], and later by Slack and Hussain [25]. Following
this path, one looks for materials with complex structures and/or with large number of atoms per
unitcell since the thermal conductivity is predicted to be inversely proportional to the number of
atoms per unitcell[26]. This is the rationale for studying some of materials discussed in chapter
4 of this thesis. Recently, Skoug et al. [27], Zhang et al. [28] brought back the important relation
between lone-pair electrons and the aharmonicity of phonon spectrum of materials. Zhang et
al. calculated the phonon spectrum and suggested that lone-pair electrons created arharmonicity
of and softened phonon modes, producing large Grüneisen parameter, and thus lowing thermal
conductivity due to phonon-phonon scattering. Chapter 4 will discuss more on this mater for
tetrahedrally bonded compounds.

1.3 Narrow band gap semiconductors and pseudo-gap systems
The area of semiconductor is one of the fastest developing fields in both science, technology and
industry which totally change the life of the mankind. The first production of semiconductor transistor was announced just sixty years ago. Nowadays, it involves virtually every thing around us
from simple daily task to the way we communicate and work. It is hard to imagine how the world
would turn out to be if not for the development of semiconductor. The first recorded observation
of semiconductor dates back to 1833 when Michael Faraday reported that the electrical conductivity of sulphurette of silver (silver sulfide) increases with increasing of temperature[29] which he
thought to be related to the tension. The properties of semiconductor came to be understood about
a century later, thanks to the quantum mechanics study by Wilson in 1931.[30] Semiconductor
is now known to have a forbidden gap in the energy spectrum, in which the electronic states are
not allowed to exist. The increase in the conductivity of a semiconductor with increasing temperature is known to be due to the excited electrical carriers when electrons jump from the valence
bands to the conduction bands. There is a thin line that distinguishes semiconductor from insulator
13

3
AlP

Band gap Eg (eV)

2.5

AlAs
GaP

2

AlSb
GaAs

1.5

CdTe

InP

Si

1

GaSb

Ge

PbS InAs

0.5

PbTe

PbSe

0

InSb
HgTe

-0.5
5.2

5.4

5.6
5.8
6
6.2
Lattice constant (Angstroms)

6.4

6.6

Figure 1.7: (Source: Lin and Kuo [31]) Band gap of several typical semiconductors. The shade
indicates the range of narrow-band-gap semiconductors.

since both have a forbidden band gap in the electronic structure. The gap of semiconductor is usually smaller than that of insulator (∼5.5 eV for diamond) so that the excitation occurs at ordinary
temperatures. Figure 1.7 shows values of band gap for several well-known semiconductors.
A vast amount of work has been devoted to the understanding of semiconductor or related
effects, covering from fundamental physics to technological application. Many aspects of semiconductor have been revealed and many more are being explored. Examples of the applications
and interesting phenomena in semiconductors can be found in overwhelming numbers of publications in the literature, including, but not limit to, thermoelectrics, spintronics, electronics, etc.
. Recently, a new phenomenon that has attracted a lot of attention lately is in the field of topological insulator (TI) discovered in 2005.[32]. Except for its misleading name, TIs are small band
gap semiconductors with topologically protected surface states. In this thesis, I only deal with
semiconductors for thermoelectric applications.
Semiconductors can be classified into two categories, direct and indirect band gap. In direct-

14

Indirect

DOS

Direct
Narrow-gap

.

EF Eg

DOS
.

.

EF Eg
.

.

.

Metal

Indirect

.

Pseudo-gap

Wide band gap

Direct

EF

.

.

EF

.
.

Figure 1.8: Diagram explaining different scheme of band gap in a semiconductor.

band-gap semiconductors (for example, GaAS, InP, GaN, etc), the lowest of the conduction bands
and the highest of valence bands occur at the same symmetric point in the first Brillouin Zone
(IBZ), whereas in indirect-band-gap semiconductors, they occur at different symmetric points
(see Fig. 1.8). Another way to classify semiconductors is based on the magnitudes of their band
gaps, one has either large-band-gap, narrow-band-gap or pseudo-gap semiconductors. A diagram
demonstrate the different types of band gap in semiconductor is given in figure 1.8. The definition
of narrow-gap semiconductor is not quite strict. Since the history of narrow-band-gap semiconductors (NBGS) came along with its application in infrared emitting and detecting, usually NBGS
is limited to those with the band gap (Eg ) corresponding to an infrared absorption cut-off wavelength over 2µ m (∼0.6 eV).[33] However, one could also refer to NBGS as those having band
gap smaller than that of silicon (∼1.1 eV). The renewed interest in NBGS for its potential thermoelectric application was inspired by the works of Mahan and Sofo in from 1989–1994[18, 34]
who suggested that NBGS are the best candidates for TE application. Beside NGBS, pseudo-gap

15

materials are also interesting systems which can be used in thermoelectric applications. Pseudogap materials have an overlap between conduction bands and valence bands, giving rise to a small
but finite density of states at the Fermi level. Our group at Michigan State University has been
working extensively on the branch of narrow-band-gap (NBGS) and pseudo-gap semiconductors
with the focus on their thermoelectric (TE) application (an introduction on thermoelectrics is given
in section 1.1). We have been working on a wide range of materials including well known NBGS
such as BiTe, BiSe, SbTe, AgPbm SbTem+2 (LAST), III-VI compounds, Mg2 Si, Cu3 SbSe4 , etc.

1.4 Outline of the dissertation
In this dissertation, I have used first principle calculations using density functional theory (DFT),
which is discussed in chapter 2, and reported the results of my works done in the last several years
on the electronic and thermoelectric properties of selected NBGS including Heusler compounds
(mainly Fe2 VAl), tetrahedrally bonded I3 -V-VI4 chalcogenides (I = Cu, V = P, As, Sb, and Bi, and
VI = S, Se, and Te), and ZrNiSn-ZrNi2 Sn mixture.
The results can be divided into two parts: the first part includes chapter 3 and 4, and the second
part includes chapter 5. The first part focuses on understanding the fundamental properties of
Heusler compound Fe2 VAl and the tetrahedrally coordinated compounds. I address the questions
about the band gap of these materials using different levels of approximation in DFT and study
their thermoelectric properties using Boltzmann’s transport equation. Most of the works in this
part have been published in references [35], [36], and [37]. The second part dedicates for studying
the formation of the Heusler phase nanostructures in the Half-Heusler matrix or vice versa, and
analyzing their electronic structures with a discussion on how they can affect the thermoelectric
properties of these materials.
In chapter 6, I will summarize my findings and propose some further works.

16

CHAPTER 2
THEORETICAL METHODOLOGY

Electronic structure plays a key role in understanding physical properties of materials, in particular,
their thermoelectric properties. The methods of calculating electronic structure, especially the
Density functional theory (DFT), have developed rapidly during the last quarter of the last century,
due to the development of computing techniques and the power of modern computers. There is a
large selection of available DFT-based computational codes (programs) which can be easily used
to calculate electronic structure of materials. DFT has become the most successful tool in studying
physical properties materials. Even though DFT is an exact theory, all the calculations are done
using, at least, some approximations. Each approximation is suitable for some specific cases,
thus, it is very important to understand the theoretical basis behind these calculations and their
limitations. In this chapter I shall review some of the fundamentals of electronic structure theory,
focusing on DFT, and discuss several approximations to it. The Boltzmann’s transport equation
will also be discussed as a bridge connecting calculated electronic structure to the thermoelectric
(transport) properties of a material.
Some of the contents of this chapter follow the spirit of the "Electronic Structure: Basic Theory
and Practical Methods" by Martin [38].

2.1 Basic equations
Electronic structure, following the discovery of electrons in the late 1800s, is the concept used for
the electron configuration, or the arrangement of electrons in atoms, molecules, and solids. Theory
of electronic structure is a combination of quantum mechanics and statistical physics.
The problem of dealing with a realistic electron system is that it is a many-body system. In order to understand the behavior of such electrons, one has to start with the Hamiltonian for electrons

17

and nuclei,
2
e2
h2 2 1
Z I Z J e2
¯2
¯
ˆ = − h ∑ ∇2 − ∑ ZI e + 1 ∑
H
−∑
∇I + ∑
,
2me i i i,I |ri − RI | 2 i̸= j ri − r j
2 I̸=J |RI − RJ |
I 2MI

(2.1)

where electrons, with charge e and mass me , are denoted by lower case subscripts, and nuclei, with
charge ZI and mass MI , are denoted by upper case subscripts. For simplicity, we adopt the Hartree
atomic units where h = me = 4π /ε0 = 1. Using the Born-Oppenheimer approximation[38], one
¯
can neglect the kinetic energy of nuclei and rewrites the Hamiltonian as
ˆ
ˆ ˆ
ˆ
H = T + Vext + Vint + EII ,

(2.2)

1
ˆ
T = − ∑ ∇2 ,
2 i i

(2.3)

where the electron kinetic energy is

ˆ
Vext is the potential due to the nuclei, given by
ZI
,
i,I |ri − Ri |

ˆ
Vext = ∑ VI (|ri − Ri |) = ∑
i,I

(2.4)

ˆ
Vint is the electron-electron interaction,
1
1
ˆ
Vint = ∑
,
2 i̸= j ri − r j

(2.5)

and EII is the interaction between static nuclei. In general, EII can include other terms that contribute to the total energy but does not directly affect the electronic behavior. Also, Eext can have
other "external potentials" acting on electrons, for instance, electric fields, magnetic fields, etc.
The wavefuntions describing electrons is the solution of the time-dependent Schrödinger equation
i¯
h

dΨ({ri };t)
ˆ
= HΨ({ri };t),
dt

(2.6)

where the many-body state Ψ is a function of electron position ri and time t. It is noted that the
spin subscript is omitted and assumed to be included in the coordinate ri . The wavefunction can
be written as Ψ({ ri } ;t) = Ψ({ ri })e−iEt , and it also satisfies the time-independent Schrödinger
equation
ˆ
H |Ψ⟩ = E |Ψ⟩ ,
18

(2.7)

where E is the energy associated with an eigenstate |Ψ⟩. The ground state wavefunction Ψ0 is the
state with the lowest energy and can be determined by minimizing the total energy with respect to
parameters appearing in |Ψ⟩.
The time-independent expression for an observable is given by the expectation value of the
ˆ
corresponding operator O,
ˆ
⟨O⟩ =

ˆ
⟨Ψ O Ψ⟩
.
⟨Ψ|Ψ⟩

(2.8)

Applying this equation to the density operator n(r) = ∑i=1,N δ (r − ri ), one gets an important
ˆ
quantity in the electronic structure theory, the density of particles,
∫

d 3 r2 · · · d 3 rN ∑σ1 |Ψ(r, r2 , r3 , ..., rN )|2
⟨Ψ |n| Ψ⟩
ˆ
n(r) =
= N∫
,
⟨Ψ|Ψ⟩
d 3 r1 , d 3 r2 · · · d 3 rN |Ψ(r1 , r2 , r3 , ..., rN )|2

(2.9)

where N is the total number of electrons, σ1 is the spin index associated with r. The integral
form of the density is obtained by tracing over all variables, taking into account the symmetry of
wavefunction over electron positions. If the σ1 index is omitted, one gets the density of one spincomponent. Similarly, the total energy of the system is obtained by calculating the expectation
value of the Hamiltonian,
ˆ
⟨Ψ H Ψ⟩
ˆ
ˆ
E=
= ⟨T ⟩ + ⟨Vint ⟩ +
⟨Ψ|Ψ⟩

∫

d 3 rVext (r)n(r) + EII ,

(2.10)

It is convenient to organize the terms in the energy into charge-neutral groups,
(
)
ˆ
ˆ
E = ⟨T ⟩ + ⟨Vint ⟩ − EHartree + E CC ,
where
1
EHartree =
2

∫

d 3 rd 3 r′

n(r)n(r′ )
,
|r − r′ |

(2.11)

(2.12)

is the self-interaction energy of the density n(r) treated as a classical charge density, and
E CC = EHartree +

∫

d 3 rVext (r)n(r) + EII ,

(2.13)

ˆ
is the classical Coulomb energy. The middle term in 2.11, ⟨Vint ⟩ − EHartree , is the difference
between the Coulomb energies of interacting, correlated electrons with density n(r) and that of a
19

continuous classical charge distribution having the same density. This term defines the exchangecorrelation energy, Exc , which is an essential and perhaps the most problematic part in DFT. A
discussion of this term is the central focus of this chapter. It is noted that all the Coulomb and
Hartree terms are long-range interactions, whereas the exchange-correlation term is short-range,
since the long-range part is canceled in the difference.

2.2 Independent-electron – Hartree and Hartree-Fock approximation
It is impossible to solve the Schrödinger equation exactly, especially for a large many-electron
system. It is important to find an effective approximation. One approximation which is used
extensively till today is the independent-electron approximation (mean-field theory). This approximation assumes that the electrons are uncorrelated except that they follow the exclusion principle
and an electron moves in an average potential provided by the others. There are two basic approaches which characterize independent-electron approximation, they are Hartree, and HartreeFock approximations.

2.2.1

Hartree approximation

In the Hartree approximation, the effective potential includes the effect of real interactions in a
averaging way but does not include the Pauli principle explicitly in the wavefunction (Fermion
antisymmetry). The many-electron wavefunction can be considered as a multiple of the singleparticle wavefunctions,
N

Ψ({ ri }) = ∏ ψi (ri ),

(2.14)

i=1

where N is the number of electrons. The many-body problem, thus, becomes an effective singleelectron equation,

[

]

1
ˆ
He f f ψiσ (r) = − ∇2 +Veσf f (r) ψiσ (r) = εiσ ψiσ ,
2

(2.15)

where Veσf f (r) is an effective potential acting on electrons of spin σ and position r. The ground
state for many non-interacting electrons is obtained by filling the lowest eigenstates obeying the

20

exclusion principle. The use of this approximation is the core of DFT which will be discussed later
in this chapter.
At finite temperature, it is straight forward to apply statistical mechanics of non-interacting
fermions, with the Fermi-Dirac equilibrium distribution of electron,
f (εiσ ) =

1

e

β (εiσ −µ )

,

(2.16)

+1

where µ is the chemical potential of the electrons. One can get the total energy by weighted sum
of the eigenvalues of 2.15
ˆ
E(T ) = ⟨H⟩ = ∑ f (εiσ )εiσ

(2.17)

i,σ

And the one-electron density operator can be defined as
ˆ
ρ = ∑ |ψiσ ⟩ fiσ ⟨ψiσ | ,

(2.18)

i

and one can get the density matrix for explicit spin and position representation,

ρ (r, σ ; r′ , σ ′ ) = δσ ,σ ′ ∑ ψiσ ∗ (r) f (εiσ )ψiσ (r′ ),

(2.19)

i

the diagonal part gives the electron density,
nσ (r) = ρ (r, σ ; r, σ ) = ∑ f (εiσ ) ψiσ (r) .
2

(2.20)

i

An example of the simplest model system in condensed matter is the homogeneous electron
gas. In this model the nuclei are replaced by a uniform positive charge. The system is specified by
its density n = N/Ω, which can be characterized by the parameter rs such that
4π 3
1
rs = Ω/N = .
3
n
The Hamiltonian for this system is given by
[
]
∫
1
1
1
n2
ˆ
H = − ∑ ∇2 +
− d 3 rd 3 r′
2 i i 2 i̸∑j |ri − r j |
|r − r′ |
=
˜
which can be rewritten in terms of r = r/rs as
(
)]
( )2 [
∫
a0
1 ˜ 2 1 rs
1
3
2r 1
ˆ
H=
∑ − 2 ∇i + 2 a0 ∑ |˜ i − r j | − 4π d ˜ |˜ | ,
rs
r ˜
r
i
j̸=i
21

(2.21)

(2.22)

(2.23)

where a0 is the first Bohr radius. This expression shows that the properties of a system with
the density rs /a0 are equivalent to a system at fixed density but with scaled electron-electron
interaction e2 → (rs /a0 )e2 .
Applying the non-interacting approximation, the solutions of the Schrödinger equation are nor1
malized plane waves ψk = ( 1/2 )eik·r with energy εk = 1 k2 , where k is the wavevector. The
2
Ω
density for each spin can be found

ρ (r, r′ ) = ρ (|r − r′ |) =

1
(2π )3

∫

′

dk f (εk )eik·(r−r )

(2.24)

using r = |r − r′ | one finds
]
3 [
kF
sin(kF r) − kF r cos(kF r)
,
ρ (r) = 2 3
eπ
(kF r)3

(2.25)

where kF is the Fermi wavevector given by
4π
(2π )3
(kF )3 =
N.
3
Ω
2.2.2

(2.26)

Hartree-Fock approximation (HFA)

In the Hartree-Fock approach, one starts from a properly antisymmetrized determinant wavefunction (satisfying Fermi statistics) with fixed number of electrons, N, and finds the single determinant
that minimizes the total energy for the full interacting Hamiltonian. In the absence of spin-orbit
interaction, the determinant wavefunction can be written as

ϕ1 (r1 , σ1 )
Ψ({ ri , σi }) =

ϕ2 (r1 , σ1 )
.
.
(N!)1/2
.
1

ϕ1 (r2 , σ2 ) · · · ϕ1 (rN , σN )
ϕ2 (r2 , σ2 ) · · · ϕ2 (rN , σN )
,
.
.
..
.
.
.
.
.

ϕN (r1 , σ1 ) ϕN (r2 , σ2 ) · · · ϕN (rN , σN )

22

(2.27)

where ϕi (ri , σi ) = ψiσ (r j )αi (σ j ), αi (σi ) is the spin eigenstate. The expectation value of the Hamiltonian is given by
ˆ
⟨Ψ|H|Ψ⟩ = ∑

[

∫

i,σ

1
+
2 i, j,∑,σ
σ
i

−

]

1
drψiσ ∗ (r) − ∇2 +Vext (r) ψiσ (r) + EII
2

1
2 i,∑
j,σ

∫

∫

σ j∗
σj ′
1
σ∗
σi
drdr′ ψi i (r)ψ j (r′ )
′ | ψi (r)ψ j (r )
|r − r

j

drdr′ ψiσ ∗ (r)ψ σ ∗ (r′ )
j

1
ψ σ (r)ψiσ (r′ ),
|r − r′ | j

(2.28)

where the first term is the sum over the one-body expectation values, the second term is the interaction energy between ions which can be neglected as discussed earlier, the third term is the
direct or Hartree energy, the forth term is the exchange or the Fock part. The Hamiltonian includes
unphysical self-interaction energy (i = j, σi = σ j ), however, they are canceled between the Hartree
and the Fock terms. The exchange potential exists only between electrons with same spins and is
non-local, making it more difficult to calculate. Variation of ψiσ ∗ (r) then leads to the Hartree-Fock
equation

− 1 ∇2 +Vext (r) + ∑
2
j,σ



∫

σ j∗

dr′ ψ j

j

−∑

∫

σj

(r′ )ψ j (r′ )

1  σ
ψi (r)
|r − r′ |

dr′ ψ σ ∗ (r′ )ψiσ (r′ )
j

j

1
ψ σ (r) = εiσ ψiσ (r).
|r − r′ | j

(2.29)

This equation can be written as
[

]

1
ˆi
ˆ i,σ
He f f ψiσ (r) = − ∇2 + Ve f f (r)
2

ψiσ (r) = εiσ (r),

(2.30)

where
ˆ i,σ
ˆ i,σ
Ve f f (r) = Vext (r) +VHartree (r) + Vx (r),

(2.31)

ˆ i,σ
Vx (r) = −

(2.32)

∑
j

∫

ψ σ (r)
j
′ ψ σ ∗ (r′ )ψ σ (r′ ) 1
dr j
i
|r − r′ | ψiσ (r)

For the homogeneous electron gas, the Hartree-Fock equation can be solved analytically, leading to
k
1
εk = k2 + F f (x),
2
π
23

(2.33)

Figure 2.1: The factor f (x) in Hartree-Fock approximation for homogeneous electron gas.

where

(

1 − x2
1+x
f (x) = − 1 +
ln
2x
1−x

)
,

(2.34)

with x = k/kF .
Figure 2.1 shows that f (x) is always negative with f (0) = −2 and f (x → ∞) = 0. There is a
singularity at the Fermi level (x = 1), however the value of f (x) is well defined, f (x → 1) = −1.
Thus, the band width is boarder (by ∆W = kF /π for homogeneous electron gas) in Hartree-Fock
than in Hartree approximation. The singularity at the Fermi level is due to the long-range Coulomb
interaction. This problem leads to zero density of state at the Fermi level. This failure of HartreeFock approximation persists even for electron in a solid if the system is metallic. It can be avoided
if there is a finite band gap (i.e. insulator or semiconductor) or if the Coulomb interaction is
screened to be an effective short-range interaction.

24

2.3 Exchange and correlation
As mentioned above, the key problem in electronic structure theory is the exchange-correlation
part. The interactions between electrons always involve pairs of electron, thus, the two-body
correlation function is sufficient to describe many properties. If we denote n(r, σ ; r′ , σ ′ ) as the
probability of finding two electrons with spins σ and σ ′ , at positions r and r′ respectively, and
n(r, σ ) as the probability of finding an individual electron with spin σ at position r, the measure
of correlation function for a pair of electrons is given by
∆n(r, σ ; r′ , σ ′ ) = n(r, σ ; r′ , σ ′ ) − n(r, σ )n(r′ , σ ′ )

(2.35)

The electron correlation is neglected in Hartree-Fock approximation, (except that required by the
exclusion principle) and must be distinguished from the exchange hole in HFA, which is given by
2

∆nHFA (r, σ ; r′ , σ ′ ) = −δσ ,σ ′

∑

ψiσ ∗ (r)ψiσ (r′ )

with
nHFA (r, σ ; r′ , σ ′ ) =

,

(2.36)

i

ψi (r, σ ) ψi (r′ , σ ′ )

2

1
2! ∑ ψ (r, σ ) ψ (r′ , σ ′ )
i, j
j
j

(2.37)

On the other hand, the exchange term explicitly appears in the Hartree-Fock approximation, it
includes two effects: Pauli exclusion and the self-term (i = j in eqn. 2.28) which must cancel the
unphysical self-term in the direct Coulomb Hartree energy. The exchange term lowers the total
energy, which may be interpreted as the interaction of each electron with a positive "exchange
hole" surrounding it. The exchange energy in HFA can be written in terms of the exchange hole as
∫
∫
[
]
1
∆n
(r, σ ; r′ , σ ′ )
ˆ
Ex = ⟨Vint ⟩ − EHartree (n) HFA = ∑ d 3 r d 3 r′ HFA
2 σ
|r − r′ |

(2.38)

As one can see from Eqn. 2.38, calculation of Ex involves computation of a non-local quantity
∆nHFA (r, σ ; r′ , σ ′ ) which is computationally demanding. As we will see later, DFT dramatically
reduces this computational complexity.

25

2.4 Fundamentals of density functional theory (DFT)
The original form of density functional theory was proposed by Thomas[39] and Fermi[40] in
1927, shortly after Schrödinger [41] (1926) published his paper marking the beginning of wave
mechanics. However it did not attract much attention due to the rough approximation which
missed essential physics and chemistry, such as shell structures of atoms and binding energy of
molecules. The modern formulation of DFT was formulated about 40 years later by Hohenberg
and Kohn (1964)[42], who considered the electron density a basic (fundamental) variable from
which other properties can be calculated. Later, Kohn and Sham (1965)[43] developed the theory
in terms of solving a set of single-particle equations, also known as Kohn-Sham equations. One
of the reasons which makes DFT a powerful tool in condensed matter theory, as Kohn put it in
his Nobel lecture[44], is that it breaks the exponential wall in many-electron wavefunction method
which makes it extremely difficult to calculate the electronic structure for a system with number
of electrons in the neighborhood of N0 ≈ 10 or larger. Thanks to the development in numerical
techniques as well as the improvement of modern computers, DFT can handle much larger number
of electrons.
The DFT is based upon two Hohenberg-Kohn theorems[42]. Here, only the theorems are
presented for the sake of completeness. The proofs of these theorems are fairly simple and can be
found in Hohenberg and Kohn’s original paper or in the reference [38] or [44].
Theorem 1. For any system of interacting particles in an external potential Vext (r), the potential
Vext (r) is determined uniquely, except for a constant, by the ground state particle density n0 (r).
Since the Hamiltonian is thus fully determined, except for a constant shift of the energy, it
follows that the many-body wavefunctions for all states (ground and excited) are also determined.
It naturally follows that
Corollary 1. All properties of the system are completely determined given only the ground state
density n0 (r).

26

This theorem, however, only states that in order to understand properties of a system, only the
ground state density is required but does not give any guidance on how to find that density, which
is given in the second theorem.
Theorem 2. A universal functional for the energy E[n] in terms of the density n(r) can be defined,
valid for any external potential Vext (r). For any particular Vext (r), the exact ground state energy
of the system is the global minimum value of this functional, and the density n(r) that minimizes
the functional is the ground state density n0 (r).
This leads to a lemma
Corollary 2. The functional E[n(r)] alone is sufficient to determine the exact ground state energy
and density.
One can easily generalize these theorems to a system with several types of particles since they
only change the term

∫ 3
d rVext (r)n(r) in the Halmiltonian. Any such pair of external potential

and density will obey Hohenberg-Kohn theorems. For example, one can formulate the spin density
functional theory for a system under external magnetic field where spin up and spin down electrons
are treated as different types of particles. Within this theory, the particle density is the sum of upand down-spin density,
n(r) = n↑ (r) + n↓ (r),

(2.39)

s(r) = n↑ (r) − n↓ (r).

(2.40)

E = EHK [n, s] = E ′ [n],

(2.41)

and the spin density is

The energy functional now becomes

where n in E ′ [n] is a function of both position r and spin σ .
Even with two Hohenberg-Kohn theorems and their lemmas, it is still difficult to solve the
Schrödinger equation to get the energy and the density of electrons. Another essential part in DFT

27

is the Kohn-Sham ansatz[43] which, in principle, leads to exact calculation of electronic structure
for many-body system using independent-particle methods. Kohn and Sham proposed that the
ground state density of an interacting system is equal to that of a suitably chosen non-interacting
system, leading to soluble independent-particle equations, with all the many-body terms included
in a exchange-correlation functional of the density. This approach has proven to be very powerful
and has become the basis of most “first-principles” or “ab-initio” calculations.
The Kohn-Sham ansatz is based on two assumptions:
1. The exact ground state density can be represented by the ground state density of an auxiliary
non-interacting system.
2. The auxiliary Hamiltonian is chosen to have the usual kinetic term and an effective local
potential Veσf f (r).
Then, the density of particles and the ground state energy functional in Kohn-Sham approximation can be written as
Nσ

n(r) = ∑ n(r, σ ) = ∑ ∑ |ψiσ (r)|2 ,
σ

∫

EKS = Ts [n] +

(2.42)

σ i=1

drVext (r)n(r) + EHartree [n] + EII + Exc [n],

(2.43)

where N σ is number of electrons with spin σ , Ts is the kinetic energy of the independent-electron
system, given by
σ

σ∫

1 N
1 N
Ts = − ∑ ∑ ⟨ψiσ |∇2 |ψiσ ⟩ = ∑ ∑
2 σ i=1
2 σ i=1

d 3 r|∇ψiσ (r)|2 ,

(2.44)

the Hartree energy is given by
1
EHartree [n] =
2

∫

d 3 rd 3 r′

n(r)n(r′ )
.
|r − r′ |

(2.45)

All the effects of many-body exchange and correlation are incorporated into Exc which can be
written as
ˆ
ˆ
Exc [n] = ⟨T ⟩ − Ts [n] + ⟨Vint ⟩ − EHartree [n].
28

(2.46)

So, in principle, if Exc is known the the exact ground state energy and density of the particles can
be found.
Using variational method, one gets the Kohn-Sham equations
σ
(HKS − εiσ )ψiσ (r) = 0,

(2.47)

where
1
σ
σ
HKS (r) = − ∇2 +VKS (r),
2
1
σ
= − ∇2 +Vext (r) +VHartree (r) +Vxc (r).
2

(2.48)
(2.49)

The eigenvalues of the Kohn-Sham equations have, in general, no physical meaning with only
one exception – the highest eigenvalue in a finite system, which is minus of the ionization energy
and is exact.[45]
σ
In order to find the form of Exc and Vxc , one can approximate them as a local or nearly local

functional of the density. The exchange-correlation energy can be written as
∫

Exc [n] =

dr n(r)εxc ([n], r),

(2.50)

where εxc is an energy per electron. Note that for simplicity only non-polarized case is considered.

εxc is related to the exchange-correlation hole and their relation can be found using the "adiabatic
connection" formula with the scaled electronic charge λ e2 , where λ varies from 0 (non-interacting)
to 1 (actual value),
∫ 1

dV
1
Exc [n] =
d λ ⟨Ψλ | int |Ψλ ⟩ − EHartree [n] =
dλ
2
0

∫

d 3 rn(r)

∫

d 3 r′

nxc (r, r′
¯
|r − r′ |

(2.51)

in which, nxc (r, r′ ) is called "coupling-constant-averaged hole",
¯
nxc (r, r′ ) =
¯

∫ 1
0

d λ nλ (r, r′ ).
xc

(2.52)

Here nxc (r, r′ ) is the exchange hole, given in eqn. 2.37, summed over parallel and anti-parallel
spin. Finally one gets the expression for the exchange-correlation density
1
εxc ([n], r) =
2

∫

29

d 3 r′

nxc (r, r′
¯
.
|r − r′ |

(2.53)

This equation shows that the exact exchange-correlation energy can be considered as the average of the exchange-correlation hole over the interaction from independent-particle (e2 = 0) to
fully interacting electrons (e2 = 1). This becomes an important rationale for improvements of the
approximation for Exc , i.e. the hybrid functional approach which will be discussed later in this
chapter.
σ
The exchange-correlation potential Vxc now can be found from the functional derivative of Exc ,
σ
Vxc (r) = εxc ([n], r) + n(r)

δ εxc ([n], r)
.
δ n(r, σ )

(2.54)

It is important to note that the second term in 2.54 (sometimes called the "response potential") is
discontinuous at a band gap. This leads to a "derivative discontinuity" that changes all the electrons’ potential by a constant when a single electron is added to the system. Thus, even with exact
Kohn-Sham theory, the obtained band gap is not the same as the actual one. This discontinuity can
be understood by analyzing the kinetic energy of the independent-electron system Ts in Eqn. 2.44.
Since, ψiσ are different for different bands, Ts changes discontinuously when going from an occupied to an empty band. Thus it has discontinuous derivative with respect to density for filled
bands. It can also be seen that the exact exchange-correlation functional must have discontinuous
derivative. This discontinuity, however, cannot be incorporated into simple functionals such as
local density (LDA) and generalized gradient (GGA) approximations (which are discussed later in
this chapter), leading to what has been known as "Band-gap problem" in semiconductors.

2.5 Functionals for exchange and correlation
The central problem in DFT is how to treat the exchange-correlation energy Exc or the exchangeσ
correlation potential Vxc (r). This has been an active field in the DFT in the last several decades.

All treatments of exchange-correlation functional are approximated and can be classified into two
main categories: (1) local and semi-local approximation, i.e. local (spin) density approximation
(LDA/LSDA), generalized gradient approximation (GGA), meta-GGA and (2) non-local approximation such as hybrid functional (hyper-GGA), random phase approximation, etc.
30

Heaven of exact DFT
unoccupied

RPA+...
exact
Ex

τ
∇n
n

Hybrid GGA

.

Meta-GGA
GGA
LSDA
Hartree world

.
Figure 2.2: The ladder of density functional theory. Notations on the left side represent the terms
added to the theory, i.e. only density n functional is considered in LSDA, GGA adds gradient
correction ∇n, meta-GGA adds kinetic energy density τ , hybrid GGA takes into account parts of
exact , RPA includes unoccupied orbital.
exact non-local exchange energy Ex
The improvement of exchange-correlation approximations can be visualized as climbing a ladder (sometimes called Jacob ladder)[46] – which can be sometimes slippery – of which the ground
is mean field Hartree system and the top is the exact Kohn-Sham exchange-correlation functional
(figure 2.2). As one climbs up the ladder, the exchange-correlation functionals become more complicated, more sophisticated and, in general, more accurate (but not always). (list of some approximations are given in Table 2.1) The first three rungs of the ladder are local/semi-local approximations where the first rung of the ladder is L(S)DA in which exchange-correlation is functional of
only the local spin density nσ (r), the second rung is GGA which adds the gradient correction of the

31

Table 2.1: Summary of approximations used in DFT
Rung
LDA
GGA

Meta-GGA
Hybrid

Name
Exchange Correlation
a
PZ
PZ
PZ
b
PBE
PBE[47]
PBE[47]
c
PW91
PW91[48] PW91[48]
EV93d
EV93[49]
LSDA
e -PBE
WC
WC[50]
PBE[47]
f
TPSS
TPSS[51] TPSS[51]
g
mBJ
mBJ[52]
PBE[47]
B3LYP
Bh [53]
LYPi [54]
B3PW91
B[53]
PW91[48]
B1LYP
B[53]
LYP[54]
B1WC
WC[50]
PW91[48]
PBE0
PBE[47]
PBE[47]
j
HSE
PBE[47]
PBE[47]

a Perdew and Zunger

%HF/screening
–
–/–
–/–
–/–
–/–
–/–
–/–
16-20/No[55]
16-20/No
/No[56]
16-20/No[57]
25/No[58]
25/Yes[59–61]

f Tao, Perdew, Staroverov and Scuseria
g modified Becke-Johnson potential

b Perdew, Burke and Ernzerhof
c Perdew and Wang

h Becke

d Engel and Vosko
e Wu and Cohen

i Lee, Yang and Parr
j Heyd, Scuseria and Ernzerhof

density, and the third rung is meta-GGA, where the kinetic energy dependence is included. These
semi-local approximations work well for many systems where the density varies slowly. They fail
where one has a non-uniform, high density and rapidly varying regions, or has an open systems
where the number of electrons fluctuates. The fourth and the fifth rungs of the ladder are non-local
approximations where the non-local exchange is included. Computationally, the difficulty of calculations changes slightly at the first three rungs of approximation but then increases steeply by
the fourth rung. It is one of the reasons why the non-local approximation was not widely used until
last several years when the computing technologies became affordable for such calculations.
It is important to point out that while for the first three and the fifth rungs, the functionals do
not involve any fitting, for the forth one, empirical fitting is unavoidable.
The next section gives a brief introduction to some of the approximations including LSDA,
GGA, LDA+U and hybrid functional. These approximations have been used in my calculations of
electronic structures of complex multi-component systems. Descriptions of other approximation

32

can be found in literature, such as in reference [38].

2.5.1

Local (spin) density approximation (LSDA) and the gradient correction (GGA)

The first local density approximation was proposed by Dirac [62] to be used with Thomas-Fermi
density functional[39, 40]. Later, in 1951, Slater [63], starting from Hartree-Fock approximation,
proposed the local spin density approximation (LSDA) which included spin-dependent functional.
Within LSDA the exchange-correlation can be written in general as
LSDA
Exc [n↑ , n↓ ] =

∫

d 3 rn(r)εxc (n↑ (r), n↓ (r))
∫
[
]
= d 3 rn(r) εx (n↑ (r), n↓ (r)) + εc (n↑ (r), n↓ (r))(r) ,

(2.55)

where εxc is assumed to be the same as in a homogeneous electron gas. In this equation the spinquantization axis is assumed to be the same everywhere (collinear spin). For the non-collinear
case, one can replace the density by the local spin density matrix
β

ρ αβ (r) = ∑ fi ψiα ∗ (r)ψi (r).

(2.56)

i

Now the Kohn-Sham hamiltonian becomes a 2 × 2 matrix. The only complication now is that one
has to find the local axis of spin quantization. The form of εxc is, however, the same.
LSDA works the best for systems close to a homogeneous gas (nearly free-electron metal) and
the worst for the inhomogeneous ones (atoms, molecules or localized-electron systems). The reason is that, unlike in Hartree-Fock approximation, the unphysical self-term in Hartree interaction
is only approximately canceled by the exchange term. This error is negligible in a homogeneous
(delocalized) system. Furthermore, LSDA also misses the derivative discontinuity of energy functional which generally underestimates the band gap in semiconductors and insulators.
The first step to improve LSDA is by taking into account the effect of spatially varying density (inhomogeneous) by including the gradient of density in the functional. There are various
approaches to achieve this goal, among which generalized gradient approximation (GGA) is the

33

most successful and widely used. In GGA, the Exc can be written as
GGA
Exc [n↑ , n↓ ] =

=

∫
∫

d 3 rn(r)εxc (n↑ , n↓ , |∇n↑ |, |∇n↓ |, ...),
hom
d 3 rn(r)εx (n)Fxc (n↑ , n↓ , |∇n↑ |, |∇n↓ |, ...),

(2.57)

hom is the exchange energy of the unpolarized homogeneous electron gas and F is a
where εx
xc

dimensionless function.
For the exchange part, Fx , it is easy to show that only spin-unpolarized consideration is sufficient. The expansion of Fx is given, to some lowest orders, by
Fx = 1 +

10 2 146 2
s +
s +··· ,
81 1 2025 2

(2.58)

where sm are dimensionless reduced density gradients defined as
sm =

1 |∇m n|
(2kF )m n

(2.59)

The correlation part is more difficult and can be expanded in sm as
LDA
εc
Fc = LDA (1 − 0.219.51s2 + · · · ).
1
εx

(2.60)

There are many forms for Fx (n, s), where only terms with order of s = s1 or lower are considered, higher terms (m ≥ 2) are neglected. Among those forms, there are three commonly
used forms proposed by Becke [53], Perdew and Wang [48] and Perdew et al. [47]. For small
s, (0 < s ≤ 3) – slowly varying density – which is relevant for most physical systems, different
GGAs give similar improvements. Even though GGA has been successfully used to predict properties of many system, its accuracy is still not enough and the improvement over LDA is modest.

2.5.2

Local treatment of strongly correlated electron systems: LDA+U

As discussed earlier, in LSDA, the self-interaction energy is only partly canceled, and the discontinuity in derivative of energy with respect to density at integral values of n is not included. These
problems are severe for systems with localized electrons such as d of f electrons. There exist
34

corrections to this problem within the local approximation including self interaction correction
(SIC)[64] and LDA+U[65]. While SIC can produce very well the nature of localized electrons, its
one-electron energies are usually in strong disagreement with experiments. Here we only discuss
the LDA+U approximation, which can be regarded as a first-principles method since U can be
calculated from ab-initio calculations, making it parameter-free.
LDA+U is the method that involve LDA- or GGA-type (more correctly GGA+U) calculations
with an additional orbital-dependent Hubbard-like interaction, 1 U ∑i̸= j ni n j , where N = ∑i ni is
2
the total number of localized (says d) electrons. Here, U is an effective Coulomb repulsion energy
needed to place two electrons on the same orbital, and is given by
U = E(d n+1 ) + E(d n−1 ) − 2E(d n ),

(2.61)

where d denotes localized orbitals such as those occupied by d and f electrons. One also has
to subtract the Coulomb energy which is already included in LDA, E = UN(N − 1)/2, to avoid
double counting. The energy in LDA+U, then , is given by
E = ELDA −

U
1
N(N − 1) + U ∑ ni n j .
2
2 i̸= j

(2.62)

The orbital energies εi are derivatives of the total energy with respect to orbital occupations ni :
1
εi = εLDA +U( − ni ).
2

(2.63)

This equation means that, if one neglects the hybridization between the localized and extended
states, the effect of U is to shift the LDA orbital energy by −U/2 for occupied bands (ni = 1) and
U/2 for unoccupied bands (ni = 0). Taking into account the exchange effect between the same
spin electrons, the energy becomes
E = ELDA −

J
1
U
U −J
N(N − 1) − N(N − 2) + U ∑ nσ n−σ +
nσ nσ ,
i i
2
4
2 i
2 i̸∑j i j
=

35

(2.64)

where U and J can be obtained by using the Slater integrals F k [66]
U = F 0,

(2.65)

U − J = F 0 − (F 2 + F 4 ),

(2.66)

F2 + F4
14

(2.67)

J=

One can calculate U and J by from constrained DFT calculations using supercell model. In a
method proposed by Gunnarsson et al. [67] one removes the transfer integrals between the d orbitals and the other orbitals in a system and obtains U by using the equation U = δ 2 E/δ n2 .[66, 68]
d
Various examples of LDA+U calculations are given in the reference [65]. A well-known example is the case of CuO-based superconductors[65] which are found to be non-magnetic metals in LDA and GGA calculations whereas LDA+U calculations give the correct antiferromagnetic insulators. Another example showing the severe errors in LDA/GGA which can be corrected by LDA+U is the case of CoO Mott-insulator. CoO has a large band gap of 2.4 eV[38]
which is predicted to be a metal within LDA/GGA, with Co-d state located at the Fermi level[69].
LDA+U[65, 70] correctly opens a band gap in good agreement with experiment.
Anisimov et al. [65] brought up an important relationship between LDA+U and the Green
function (GF) method; GF correctly describes the frequency-dependent quasiparticle energies. A
simple approximation to the GF method is GW approximation (GWA)[71], where G and W stand
for Green and response functions respectively. The authors note that both GW approximation and
LDA+U are Hartree-Fock-like theories which give the corrections to the electron self-energy. The
only difference between the two methods is the procedure for calculating the screened Coulomb
interaction. In LDA+U, it can be done by using constrained DFT, whereas GWA requires the
computation of a much more difficult response function W. Thus, LDA+U can be considered, at
least for localized orbitals, an approximation of GWA.

36

2.5.3

Non-local treatment of exchange-correlation in correlated electron systems: Hybrid
functional theory

In spite of the success of local DFT approximations, they fails in many cases where it is believed
that non-local effect of exchange-correlation are important. In chapter 4 I will discuss one of such
systems. To tackle this problem, hybrid functionals have been introduced. The rationale for the
hybrid functional theory arises from the equation 2.51 which can be written as
Exc =

∫ λ
0

d λ Exc,λ ,

(2.68)

where
Exc,λ = ⟨Ψλ |

dVint
|Ψλ ⟩ − EHartree [n].
dλ

(2.69)

At λ = 0 the non-interacting Kohn-Sham system is recovered and at λ = 1 it is the fully interacting
system. The first formulation of hybrid function was proposed in 1993 by Becke [72] who assumed
a linear dependence on λ and approximated the integral by
1 HF
h
DFA
Exc = (Ex + Exc ),
2

(2.70)

HF
DFA
where Ex = Exc,λ =0 is the exact exchange energy of Kohn-Sham orbital, Exc is the exchange-

correlation energy at λ = 1 in the density functional approximation (DFA), such as LDA or GGA.
Later Becke [55] parameterized the formulation including three fitting parameters (α0 , αx , αc ),
h
LDA
HF
DFA
Becke + α E ,
Exc = Exc + α0 (Ex − Ex ) + αx Ex
c c

(2.71)

Becke is the exchange energy in the form proposed by Becke himself[53]. For convenience,
where Ex

each formulation of hybrid functional method is usually named in such a way that it indicates which
exact exchange and which correlation function used. For example, two popular formulations of
such three-parameter hybrid functional are B3PW91 and B3LYP (some of other formulations are
given in table 2.1), both of which use the exact exchange proposed by [55], the former uses the
GGA exchange-correlation by Perdew and Wang [48] while the latter uses that by Lee et al. [54].
This three-parameter formulation later was simplified to one-parameter with αx = 1 − α0 = 1 − α

37

and αc = 1,[56] i.e.,
DFA
h
DFA
Exc = Exc + α (Ex − Exc ).

(2.72)

Perdew, Ernzerhof and Burke[58] suggest the use of α = 1/4. Combining with PBE[47] correlation one ends up with the popular un-screened hybrid functional PBE0[58].
In all the approximation discussed above, the bare Coulomb interactions are considered which
include short-range and long-range parts. This greatly increases the demand of computational resources for hybrid functional calculations. It is, however, known that the short-range exchange interaction decays exponentially as a function of the distance. Furthermore, the long-range Coulomb
interaction causes a divergence near the Fermi level as discussed earlier. This problem can be
avoided by screening the interaction. The formulation for hybrid functionals based on a screened
Coulomb potential was proposed by Heyd, Scuseria and Ernzerhof (HSE) in 2003[59] which was
then amended in 2006[61] (HSE06). In this formalism the Coulomb operator is split into shortrange (SR) and long-range (LR) components:
1 er f c(ω r) er f (ω r)
=
+
,
r
r
r

(2.73)

where er f c(ω r) = 1 − er f (ω r) and ω is an adjustable parameter. The short- and long-range parts
correspond to the er f c and er f terms respectively. At ω = 0, the long-range term vanishes and
a full range operator recovered, which is PBE0-type. At ω = ∞, the interaction is completely
screened so that it falls back to (semi)local approximation. In HSE06, the PBE correlation is used
and the energy functional can be written as
HF,SR

PBEh = α E
Exc
x

PBE,SR

(ω ) + (1 − α )Ex

PBE,LR

(ω ) + Ex

PBE
(ω ) + Ec .

(2.74)

Thus, in this formulation, the non-local HF exchange is included only in the SR term whereas the
LR term is that of PBE. This approach not only helps to reduce the computing demand but also
removes the Fermi surface singularity. HSE showed that screened-Coulomb-interaction hybrid
calculations greatly improved the DFT results[59–61] especially in obtaining the band gaps which
were in much better agreements with experiments. HSE-type calculations have proved to be suc-

38

cessful in many systems[59–61] with "standard" values of parameters, α = 1/4 and ω = 0.2. In
this thesis, most of the HSE calculations were done using these values.

2.6 Boltzmann’s equation for transport
As discussed in chapter 1, Boltzmann’s equation of transport is the most successful and practical
approach to study thermoelectric properties of materials. Applying this approach, one can calculate
transport quantities, particularly thermoelectric properties, from the electronics structure obtained
from DFT or similar calculations. In this section, I give a brief introduction about Boltzmann’s
transport equation in the presence of temperature gradient, more comprehensive discussions can
be found in other textbooks on solid states physics, such as “Principles of the Theory of Solids” by
Ziman [73].
For a steady state, the Boltzmann’s equation is given by
[

]
[
]
[
]
∂ fk (r,t)
∂ fk (r,t)
∂ fk (r,t)
+
+
= 0,
∂t
∂t
∂t
di f f
f ield
scatt

(2.75)

where the subscripts di f f , f ield and scatt denote the terms corresponding to diffusion, external
field and scattering respectively, fk (r,t) is the probability of having an electron (hole) with wave
vector k at position r in space at a given time t. Note that fk (r,t) is not for the equilibrium but for
0
the steady state. At equilibrium, fk (r) is the single particle Fermi distribution,
0
fk =

(

1

)

ε −µ
exp k kT (r) + 1
B

,

(2.76)

where εk is the energy corresponding to the wavevector k, µ is the chemical potential and kB is the
Boltzmann constant. (The band index has been omitted for simplicity.) We assume that each state
k is two-fold spin degenerate and ignore the spin index since we do not consider magneto transport
in this thesis. The Boltzmann’s equation states that the net rate of change in fk (r,t) is zero for any
k at given point in space.

39

In order to find the explicit solution to the Boltzmann’s equation we apply Liouville’s theorem
for r,
fk (r,t) = fk (r − tvk , 0),

(2.77)

where vk = ∂ εk /∂ k. The diffusion term in 2.75 can be written as
[

]
∂ fk (r)
∂ f ∂T ∂r
∂f
=− k
·
= − k vk · ∇T.
∂t
∂T ∂r ∂t
∂T
di f f

(2.78)

The derivative of fk with respect to T can be evaluated by using the definition of Fermi distribution,
(
)
εk −µ
(
)
exp k T
∂ fk
εk − µ
∂ fk εk − µ
B
=[
= −
.
(2.79)
(
)
]2
∂T
∂ εk
T
εk − µ
kB T 2
exp k T + 1
B

One, then, gets the rate of change of the distribution function due to diffusion:
[

]
(
(
)
)
∂ fk (r)
∂ fk
vk · ∇T
= −
(εk − µ ) −
.
∂t
∂ εk
T
di f f

(2.80)

Similarly, one can apply Liouville theorem in k-space and consider only external electric field E,
where
hk = eE,
¯˙

(2.81)

to get the formulation for the external field term,
[

)
]
(
∂ fk
∂ fk
v · E.
=e −
∂ t f ield
∂ εk k

(2.82)

The scattering term in 2.75 is more complicated and in general is energy-dependent. Within the
scope of this thesis, we make relaxation time approximation which have been successfully applied
in many system.[74] Within this approximation the scattering term can be given by
[

]
0
fk − fk
∂ fk
g
=−
= − k,
∂ t scatt
τk
τk

(2.83)

0
where gk = fk − fk , τk is the carrier’s relaxation time. Equation 2.75 can be rewritten as

(
(
)
)
(
)
∂ fk
∂ fk
vk · ∇T
g k = eτk −
v · E + τk −
(εk − µ ) −
∂ εk k
∂ εk
T

40

(2.84)

In order to make use of this equation to get other transport properties of materials one can start
with the microscopic and macroscopic Ohm’s law which relate the current density to the applied
electric field and temperature gradient. The microscopic Ohm’s law is given by
1
1
∑ e fkvk = V ∑ egkvk,
V k
k

J=

(2.85)

since there is no current at equilibrium,
0
∑ e fk vk = 0,

(2.86)

k

where V is the volume in real space.
Using the expression of gk and grouping out the common terms, one can rewrite J as


(
)
∂ fk
)
(
1
∂f
1 ∑k τk − ∂ εk (εk − µ )vk vk


(
)
J = ∑ e2 τk − k vk vk · E +
· (−∇T ) .
∂ fk
V k
∂ εk
eT
∑ k τk − ∂ ε v k v k

(2.87)

k

Comparing with the macroscopic Ohm’s law in the presences of an external electric field and
temperature gradient, given by
(
)
←
→
→
→
→
→←
J = ← E + ← S (−∇T ) = ← E + S (−∇T ) ,
σ
σ
σ

(2.88)

←
→
→
one get the expressions for conductivity, ← , and thermopower (Seebeck coefficient), S , tensors,
σ
(
)
2
∂ fk
← =e
→
σ
τ −
v v ,
V ∑ k
∂ εk k k
k
←
→
→
→ ←
S = (← )−1 A ,
σ

(2.89)
(2.90)

(
)
←
→
∂ fk
e
where A = V T ∑k τk − ∂ ε (εk − µ )vk vk .
k

For isotropic system electrical conductivity σ and thermopower S can be calculated using

σ

= e2

(
)
∂ f0
dε −
Σ(ε ),
∂ε

+∞
∫

(2.91)

−∞
+∞
∫

e
S=
Tσ

−∞

(
)
∂ f0
Σ(ε )(ε − µ ),
dε −
∂ε

41

(2.92)

where µ is the chemical potential, e is the electron charge, f0 is the Fermi-Dirac distribution
function, and Σ(ε ), called the transport distribution function, is given by
Σ(ε ) =

1
2
k)
k)
k)).
∑ v(n,⃗ τ (n,⃗ δ (ε − ε (n,⃗
3 ⃗

(2.93)

n,k

In Eqn. (2.93), n is the band index and the summation is over the first Brillouin zone (BZ), v(n,⃗
k)
is the carrier velocity of an electron in state (n,⃗ and τ (n,⃗ is the relaxation time.
k),
k)
By analogy, one can follow a similar procedure for electronic heat current and obtain the expression for zero-field electronic thermal conductivity due to electrons
(
)
1
∂ fk
←
→
κ elec =
∑ τk − ∂ εk (εk − µ )2vkvk.
VT2 k

(2.94)

From the expression of conductivity and thermopower, one can easily see the inverse relation
between them, i.e. S decreases with increasing σ , whereas the opposite can be said for κelec . These
competing relations, which are discusses in Chapter 1, make the job of optimizing thermoelectric
properties a difficult and challenging one.

2.7 Computational details
In this thesis, most of the calculations were done using the projector-augmented wave (PAW)
method[75, 76] as implemented in the VASP code.[77–79] (exceptions will be noted). In all these
calculations we used a plane-wave energy cutoff of 400 eV and an energy convergence criterion
(between two successive self-consistent cycles) of 10−4 eV (total energy/unit cell). Perdew-BurkeErnzerhof (PBE)[47] formulation was used for the exchange-correlation functional.
For comparison, the meta-GGA called modified Becke-Johnson potential (m-BJ) is also considered for some systems. These calculations were carried out using linearized augmented plane
wave (LAPW) method as implemented in Wien2K[80] package, with a plane-wave cutoff parameter RKMAX=7.
For GGA+U calculations, we used the formalism proposed by Dudarev et al.[81], where the
on-site repulsion and exchange are incorporated through a parameter Ue f f = U − J, where U and
42

J are the Coulomb and exchange terms (screened) respectively. In this approach, total energy is
rewritten as:
Ue f f ,i
EGGA+U = EGGA + ∑
2 ∑
σ
i

)

[(

∑ nσ 1,m1
i,m

m1

(
−

)]

∑

m1 ,m2

nσ ,m nσ ,m
i,m1 2 i,m2 1

(2.95)

where EGGA and EGGA+U are the energies in GGA and GGA+U approximations, respectively;
m1 and m2 are the orbital quantum numbers (m = −2, −1, 0, 1, 2 for the d states) and nσ ,m are
i,m1 2
the matrix elements of the density operator nσ associated with spin σ and site i in this basis. For
ˆi
simplicity, hereinafter, we will refer to Ue f f as U.
Boltzmann transport theory was used to calculate thermopower S as a function of carrier concentration (both n and p types) and temperature. For simplicity we used rigid band and constant
relaxation time approximations (for a justification see a recent paper by Lee and Mahanti, Phys.
Rev. B 2012) to calculate S. These calculations were done using the BoltzTrap pacakage of Madsen and Singh[68]. For accurate computation of transport velocities (and hence S) we used dense
k-meshes (ranging from 21 × 21 × 21–41 × 41 × 41 for different systems).

43

CHAPTER 3
HEUSLER COMPOUND – Fe2 VAl: EFFECTS OF ON-SITE COULOMB REPULSION

3.1 Introduction
(Full) Heusler compounds are a class of promising thermoelectric materials. They have a stoichiometric composition with the general formula X2Y Z, where X and Y are two different transition
metals (TMs) and Z is a metalloid. They crystallize in a cubic structure corresponding to the space
group L21 . After their first discovery by Friedrich Heusler in 1903, more and more Heusler-type
alloys have been found and studied extensively.[85] Among them, Fe2 VAl based compounds came
to the research community’s attention in 1997 when Nishino et al. [82] pointed out a possible
existence of d-electron heavy-fermion behaviour (commonly seen in f -electron system) in this
compound. Their heat capacity measurement at low temperature (1.6-6 K) showed a linear term
with a γ value of 14 mJ mol−1 K−2 from which a large electron effective mass (10me , where me
is the free electron mass) was estimated. Photoemission spectrum showed a clear Fermi edge indicating metal-like characteristics. In contrast, by fitting the resistivity to an exponential function
of Eg /T in the temperature range 400 to 800 K, Nishino et al. [82] found that Fe2 VAl behaved
like a semiconductor with Eg ≈ 0.1 eV. These experimental results are quite intriguing. This dual
nature was also seen in NMR experiments by Lue and Ross [83] By measuring the T -dependence
of Knight shift (Ks ) and the nuclear spin relaxation time (T1 ) and fitting these two quantities to
functions of Eg /T , they found a band gap of 0.21-0.22 eV for Ks and 0.27 eV for T1 . Surprisingly,
Table 3.1: Summary of experimental studies of energy gap in Fe2 VAl
Group - year
Nishino et al. - 1997[82]
Lue et al. - 1998[83]

Okamura et al. - 2000[84]

Measured quantity
Resistivity (400–800 K)
NMR (250–550K)
• Knight shifts (Ks )
• Spin relaxation rate (T1 )
Photo-conductivity (9–295K)

44

Eg (eV)
0.1
0.21–0.22
0.27
0.1 - 0.2

Table 3.2: Summary of previous theoretical studies on Fe2 VAl electronic structure
Group - Year

Method

a(Å)

Eg (eV)

Singh & Mazin - 1998[87]
Weinert & Watson 1998[88]
Weht & Pickett - 1998[89]
Guo et al. - 1998[90]
Kumar et al. - 2009[91]

LSDA, LAPW (LMTO)
LDA, (FLASTO)

5.76

-0.2
-0.2

g(EF )
(eV−1 / f .u.)
0.3
–

GGA, PBE
LSDA
GGA, PBE

5.76
5.68
5.712

-0.2∼-0.1
-0.2
–

0.1
0.54
0.19

in addition to the activated behaviour, they also found metallic characteristics through Korringa
law. [86] Using the Korringa relation, they found a finite DOS g(EF ) at the Fermi level (EF ), and
estimated its V-d component to be about 1.7 × 10−2 eV−1 per V atom. The actual g(EF ) should
be somewhat larger. Okamura et al. [84] using photoconductivity measurements also reported similar results for Fe2 VAl with a band gap of 0.1-0.2 eV and a finite DOS at the Fermi level. The
experimental results of the band gap and the temperature range in which they were measured are
summarized in Table 3.1. Based on the transport measurements and the finite DOS at EF , it was
argued that Fe2 VAl was a semimetal with a pseudo-gap. This implies that there is a small overlap
between the valence band and the conduction band. This overlap gives rise to a negative band
gap Eg = Ec,min − Ev,max < 0, where Ec,min and Ev,max are the minimum of the conduction and
the maximum of the valence band, respectively. However the value of this negative band gap is
not known from the experiments. Note that the notation of negative band gap Eg and the positive
pseudo-gap can be distinctly different.
The existence of a negative band gap Eg in the electronic structure of Fe2 VAl has been found
in several density functional theory (DFT) calculations within either LDA or GGA; [87–91] a
summary of these theoretical calculations is given in Table 3.2. All the LDA/GGA calculations
indicate that Fe2 VAl has a well-developed negative band gap, with the valence band maximum
(VBM) at the Γ point and the conduction band minimum (CBM) at the X point (Γ and X are two
symmetry points in the first Brillouin Zone of a FCC lattice), coming mainly from the Fe-t2g and
V-eg states, respectively. The negative band gap was found to be −0.1 ∼ −0.2 eV, and there is

45

a small but finite DOS at the Fermi level, qualitatively consistent with the Korringa law and the
linear T -dependence of the heat capacity at low T . However, as we will discuss later, the calculated
DOS at EF is an order of magnitude larger than that estimated from Korringa’s law and an order
of magnitude smaller to explain the observed large low-T heat capacity. Nevertheless, these band
structure calculations also revealed that Fe2 VAl had sharp edges in the DOS near the Fermi level
which is a desirable feature for a good thermoelectric, as suggested by Mahan and Sofo. [17]
Encouraged by the sharp edges in the calculated DOS of Fe2 VAl, thermoelectric properties
have been studied experimentally, including our own EFRC1 group at Michigan State University[19].
In addition, several attempts have been made to improve its thermopower S, for example by incorporating anti-site defects (increasing ratio of V/Fe) [92] or by doping at the Al site with Si, Ge, and
Sn. [19, 93, 94] In the present work we will not be concerned with the V/Fe anti-site defects since
they give rise to localized states near the Fermi energy and their effects on the transport properties
are difficult to calculate accurately. Our focus here is to see how well one can understand the experimental data without involving the effect of defects. Experimentally it is found that nominally
pure Fe2 VAl is a p-type thermoelectric (in agreement with Hall measurement. [95]) In the temperature range T = 100 − 400 K, S was found to be 10–50 µ V/K. When Al was partially replaced
by Si, Ge, or Sn[19, 93, 94], more electrons were put in the system and turned it into an n-type
thermoelectric. These measurements show that the highest S (in magnitude) is about -150 µ V/K,
obtained for 5–6 % of Al substitution corresponding to an electron doping of 1.1–1.3×1020 cm−3
(assuming each dopant replacing an Al donates one electron).
In spite of extensive experimental works, there is no direct spectral evidence supporting the
existence of either a real or a pseudo gap in Fe2 VAl. From a theoretical prospective, there is also no
systematic study of the relationship between the LDA/GGA electronic structure and thermoelectric
properties in the pseudo-gap regime and how does theory compare with experiment. Furthermore,
we know that LDA and GGA do not do well vis-à-vis band gaps in semiconductors due to the
self-interaction error inherent in the LDA/GGA potential and vanishing of the discontinuity of the
1 Energy Frontier Research Center, a project funded by the US Department of Energy

46

exchange-correlation potential as a function of the level occupation at EF as discussed in chapter 2.
The errors in LDA and GGA are severe in describing accurately the localized electrons in d or f
states, i.e. transition metal and rare earth compounds. Whether there will be similar problems in
pseudo-gap or zero-gap systems (when the overlap between conduction band and valence band
is small) is not known. Recently, Bilc and Ghosez [96] used B1WC hybrid functional [57] (see
table 2.1) to study Fe2 VAl. They used ratio of α = 0.16 exact exchange and found it to be a
semiconductor with a band gap of 0.34 eV. Later in this chapter I will compare this value of band
gap to that obtained using LDA+U and other method as well as hybrid functional with different
parameter α .
In this chapter, I explore the effect of Coulomb repulsion associated with the d-electrons of
Fe and V electrons Fe2 VAl using GGA+U method (see Sec. 2.5.2), which gives quite reasonable
results for relatively low computational cost. This method has been found to be quite successful
in several other Heusler compounds containing d-electrons such as Co2 FeSi, Co2 MnSi. [97, 98] I
systematically investigate the effect of the on-site Coulomb repulsion U at both Fe and V sites on
the electronic structure and thermoelectric properties of Fe2 VAl. U was first treated as a parameter
to understand how the band structure is affected by the values of U at different transition metal
sites. Then the values of U were calculated using constrained DFT method (see section 2.7). The
calculated band structure for these values of U then is used to investigate the carrier concentration
and temperature dependence of thermopower, look for their optimum values and compare with the
experiment.
The chapter is arranged as follows. In Sec. 3.2, I briefly describe the computational procedure
beyond which already discussed in section 2.7. Results and discussions on electronic structure
and thermoelectric properties are presented in Sec. 3.3 and Sec. 3.4. A brief summary is given in
Sec. 3.5. Most of materials presented in this chapter are published in the reference [35].

47

3.2 Computational details
The on-site Coulomb interaction was added using GGA+U method within the scheme suggested
by Dudarev et al. [81] (see section 2.5.2). The values of U for Fe and V are calculated using the
constrained DFT method, where the effective Coulomb repulsion energy is given by
)
)
(
(
n 1 n
n 1 n
Ue f f =εd↑
+ ,
− εd↑
+ , −1
2 2 2
2 2 2
(
(
)
)
n 1 n
n 1 n
− εF
+ εF
+ ,
+ , −1
2 2 2
2 2 2

(3.1)

where εd↑ (n↑ , n↓ ) and εF (n↑ , n↓ ) are respectively the spin-up d-eigenvalue and the Fermi energy
for the configuration of n↑ up-spins and n↓ down-spins, n is the total number of d-electrons.
Since this method of calculating U involves varying the occupation number, which is not allowed in VASP, The Wien2k package [80] was used, following the procedure suggested by Madsen
and Novák. [68] A 2 × 2 × 2 fcc supercell was used within PBE-GGA [47] with a plane wave cutoff
RKmax =7 and a Fourier expansion cutoff Gmax =9. Values of Ue f f were found to be 4.0 and 1.5 eV
for Fe and V, respectively. The value of Ue f f for V is smaller than that of Fe as expected. Note
that in the atomic limit the value of U very high. In extended systems, since the effective Coulomb
interaction is screened by the cloud of electrons, the value of U is much smaller. For example,
U ≈ 20 eV for atomic Fe, whereas it is ∼ 6.2 − 6.8 eV for Fe metal [66] and ∼ 4.8 − 7.4 eV for Fe
oxide. [68] Thus the small value of Ue f f of Fe that we obtained suggests a strong screening effect
in Fe2 VAl.
In this chapter, the effect of Coulomb repulsion was first investigated with various values of U
on Fe and V sites separately to understand how different bands are affected by U. For this purpose
three situations were consider for each atom: weak (Ue f f =1 eV), intermediate (Ue f f =2 eV), and
strong (Ue f f =4 eV) repulsion limits (comparing to typical d bandwidth of 1.5–2 eV). GGA+U
calculations were then done with U values obtained using the constrained DFT method. In the rest
of this chapter the suffix “e f f ” will be omitted from Ue f f .
Thermopower S was calculated using Boltzmann transport equation in constant relaxation time
and rigid band approximations. [74] Since Fe2 VAl is a cubic system, its transport tensors are
48

Figure 3.1: Band structure of Fe2 VAl showing d-characters of Fe and V (left); and DOS (right).
EF denotes the Fermi level set to be zero. The size of the symbols represents the strength of
orbital characters.
diagonal and all the diagonal elements are the same and electrical conductivity σ and thermopower
(Seebeck coefficient) S are given by Eqn. 2.91 and 2.92 respectively. Here τ (n,⃗ is assumed to be
k)
constant. The energy dispersion was first obtained with Monkhorst-Pack mesh of 41 × 41 × 41 kpoints, then the transport coefficients were calculated from the resulting eigenvalues employing the
Boltztrap code developed by Madsen and Singh. [74] As mentioned earlier, this chapter discusses
only the thermopower S.

3.3 Electronic structure: Pseudo-gap versus real gap
In Fig. 3.1 (left panel) the energy bands were shown together with Fe- and V-d characters along
different symmetry directions of the first BZ. The right panel gives the DOS of Fe2 VAl which
clearly shows a pseudo-gap structure near EF . The zero of energy has been chosen to be at the
Fermi level. The present GGA band structure results agree with the previous GGA calculations of
other groups [87–91] in which Fe2 VAl has a negative band gap of Eg ∼ −0.17 eV. Certain general
49

features of the band structure are the following: It has a deep and rather broad Al s-band from -10
to -6 eV. Note that here 0 eV is assumed at Fermi level. The bands from -5 to 2 eV result from
the hybridization of Fe and V d-orbitals with a small admixture of Al p-orbitals (only Fe-d and
V-d characters are shown in Fig. 3.1). The pseudo-gap is formed by a nondegenerate V-eg band
minimum at the X point and a three-fold degenerate Fe-t2g band at the Γ point. Near 0.3 eV there is
a very narrow band originating from Fe-eg states which gives rise to a sharp edge in the conduction
band DOS. These Fe- and V-d bands are important features which are relevant to the low-energy
physics in general and the thermoelectric properties in particular. Weht and Pickett [89] suggested
that the excitonic effect may be important in this region.
The band structure of Fe2 VAl near the Fermi-level shows that there are three-fold degenerate
hole pockets centered around the Γ point with nearly same effective mass and a nondegenerate
electron pockets centered around the X points. Since there are three inequivalent X points in the BZ
this compensates for the 3-fold degeneracy of valence band when one compares the contributions
from holes and electrons to the thermopower. The calculated DOS at the Fermi energy for the
undoped system is 0.27 eV−1 per formula unit (f.u.). This gives a γ value for the linear heat
capacity of 0.63 mJ mol−1 K−2 , nearly the same as 0.69 mJ mol−1 K−2 reported by Guo et al.
[90] in their LSDA calculation. The small difference between the two theoretical values may be
due to the different lattice constants and different methods which have been used. However, both
these values are an order of magnitude smaller than the experimentally measured low-T γ value of
14 mJ mol−1 K−2 .[82]
As seen in Fig. 3.1 (right panel), there is a rapid rise in the DOS near 0.3 eV (in the conduction
band) and a similar one but less rapid at -0.2 eV (in the valence bands). As we will discuss later,
one obtains large values of S (magnitude) when Fe2 VAl is heavily doped, either n-type or p-type,
so that the chemical potential approaches these rapidly increasing parts of the DOS. This result is
expected, as already discussed in section 1.2.
When the on-site Coulomb repulsion U is turned on, as one would expect (see Chap. 2) the
occupied d-states are pushed down in energy whereas the unoccupied states are pushed up. This

50

Figure 3.2: Band structure of Fe2 VAl in absolute scale for UFe of (a) 0 eV, (b) 1 eV, (c) 2 eV and
(d) 4 eV. The Fe-eg and Fe-t2g bands are pointed out by arrows.

is clearly seen in Figs. 3.2 and 3.3. In order to see how different bands change due to the on site
Coulomb repulsion, the bands structure was plotted in the absolute scale focusing on the bands near
the Fermi level in the range from -11 to -7 eV. Fig. 3.2 shows the bands structures with UFe =0,1,
2, and 4 eV at the Fe site only. In this case V bands are relatively unaffected while the unoccupied
Fe-eg bands are pushed up and the occupied Fe-t2g bands are pushed down. For a quantitative
measure of the effect of UFe , we define a quantity ∆eg ,t2g which is the splitting between the Fe-eg
conduction band and the Fe-t2g topmost valence band at the Γ point. The value of ∆eg ,t2g changes
by 0.22, 0.48 and 1.18 eV for UFe =1, 2, and 4 eV, respectively. As expected, the larger the value
of UFe , the larger is the change in the energy difference. However, the changes in the energies
are about 22∼30 % of the values of the parameter UFe . This is due to the hybridization of Fe-d,
V-d and some Al-p orbitals, which weakens the effect of U. In the weak repulsion regime, the
parameter UFe has very small effect on the band gap. In the intermediate and strong repulsion
cases the effects are significant. The relative shift of the bands results in changing the pseudo
gap to a real gap. For UFe = 1 eV the pseudo-gap remains while the values of the band gap are

51

0.002 eV and 0.27 eV corresponding to UFe = 2 eV and 4 eV, respectively.
Looking more carefully at the changes in the band structures (see Fig. 3.2), one can point out an
interesting change in Fe-t2g bands at the X point. In the absence of UFe , the t2g bands of Fe splits
into two levels, where a nondegenerate band is located at ≈0.15 eV above the two-fold degenerate
one. Note that the nondegenerate band consists of Fe-t2g only while the two-fold degenerate band
shows strong hybridization, with a small contribution of Fe-t2g and large contribution of p-like
plane-wave states. On the other hand, at the Γ point, the VBM is 3-fold degenerate coming from
hybridization of Fe- and V-t2g orbitals. In the presence of UFe , the bands without hybridization
are strongly affected by UFe , whereas the bands with stronger hybridization, specially when p
states are involved, are not affected much by UFe . This results in the change in the ordering of
these levels (at X and Γ) when UFe =4 eV. The nondegenerate band is pushed down more strongly
than the two-fold degenerate one. As a consequence, when UFe =4 eV, reversed ordering of these
levels is observed at the X point, i.e. the two-fold degenerate Fe-t2g band becomes the highest
occupied band. This makes the system a direct band gap semiconductor and also brings more hole
pockets into the system even with small doping levels. This would make Fe2 VAl a better p-type
thermoelectric (which will be discussed in detail later in this chapter). Similar effect is also seen
in Fe-eg conduction bands. Without UFe , at the X point, the second and fourth lowest-energy
bands have pure Fe-eg character while the third one is a hybridized band of Fe-eg and V-t2g . By
turning on the Fe-site Coulomb repulsion, the energy difference between the second and fourth
bands remain nearly the same, whereas that between the second and third bands is decreased.
Analogous band shifting is seen when we turn on the Coulomb repulsion at the V site only:
The Fe bands remain nearly unchanged whereas the V-eg bands, since they are mainly unoccupied,
are shifted upwards (shown in Fig. 3.3). There is however a difference in the way V d-bands react
when we include U effect on the V site. The V-eg bands, most of which lie above the Fermi level,
is affected strongly by UV since it is predominantly of V-eg character. The energy of this band is
changed about 13–20 % of the UV value. On the other hand, the occupied bands which hybridize
with Fe-t2g are shifted by less than 1 % of the UV value (Note that UFe =0). In the case when one

52

Figure 3.3: Band structure of Fe2 VAl in absolute scale for UV of (a) 0 eV, (b) 1 eV, (c) 2 eV and
(d) 4 eV.The V-eg bands are pointed out by an arrow.

considers UV at the V site only, even a weak repulsion causes a significant change in the band gap.
For UV =1 eV, Fe2 VAl becomes an almost zero-gap system; for UV =2 eV there is a real gap of
0.21 eV; when UV =4 eV the V-eg band is pushed up so much that the flat Fe-eg band becomes the
lowest conduction band and the gap is now 0.32 eV.
When we turn on U at both Fe and V sites the net effect appears to be roughly a combination
of individual Fe and V effects. For the U values obtained from constrained DFT as described in
Sec. 3.2 (UFe =4 eV and UV =2 eV), the calculated band structure is shown in Fig. 3.4(b). A band
gap of 0.55 eV is obtained and it is both direct and indirect since the highest occupied states at the
Γ and the X points are nearly degenerate (as discussed in the case UFe =4 eV above).
As mentioned in Sec. 3.1, several improvements beyond LDA/GGA have been done to take
better account of localized d-electron systems. Here I compare the results obtained using different
approximations including improved GGA by Engel and Vosko [49], hybrid functional: PBE0, [58],
HSE06[59–61] and B1WC [57] as well as the recently propose modified Becke-Johnson (mBJ)[52]
which is claimed to give good band gap values which agree well with those obtained using the GW

53

Figure 3.4: Band structure of Fe2 VAl obtained in different methods: (a) GGA-PBE , (b) GGA+U
, (c) mBJ and (d) PBE0.
mBJ and PBE0 bands structures were calculated with Wien2K

Table 3.3: Summary of Fe2 VAl band gap using different methods
Method∗
Parameter
GGA–PBE –
GGA–EV
–
mBJ
–
GGA+U
UFe = 4 eV,UV = 1.5 eV
PBE0
α = 0.25
HSE06
B1WC

α = 0.25, ω = 0.2
α = 0.16
α = 0.25

∗

‡

For detailed discussion on methods of
calculations see Sec. 2.5.3 and table 2.1
†

⊥

Calculated with VASP

54

Band gap (eV)
-0.17†
-0.06†
0.22‡
0.55†
0.58‡
1.8†
1.1†
0.34⊥‡
0.62‡

Calculated with Wien2K
Reference [96]

method. [71] The summary of methods and corresponding band gaps are given in Table 3.3 and
the band structures obtained from some of these methods are shown in Fig. 3.4. Note that some of
the results were calculated using Wien2k [80] package (see table 3.3).
We found that both PBE0 and B1WC which are implemented in Wien2K with α = 0.25 give
quite large band gaps, 0.58 and 0.62 eV, which are very close to our GGA+U result. On the other
hand, PBE0 implemented in VASP gives a really large band gap of ∼ 1.8 eV. With VASP, when
the exact exchange is screened, i.e. as in HSE06, the obtained band gap is reduced to ∼1.1 eV,
which, however is still twice of that calculated with Wien2K hybrid functional. Note that GGA
results agree extremely well between the two DFT codes. This brings up a questions of how nonlocal exchange is implemented in those programs and how to critically validate the large number
of electronic structure codes which are handily available.
Reducing the value of α to 0.16 in B1WC method results in the band gap of 0.34 eV. These
values are, however, larger than that obtained from mBJ method which gives 0.22 eV (in reasonable
agreement with experimental band gap). It is important to mention that mBJ method, among all
the methods we have tried, gives the energy level ordering at the X point very similar to that of
seen in the GGA+U calculation. However the change in the band structure seen in mBJ (compare
to GGA calculation) is not as large as that seen in GGA+U. Hence the top of the valence band is
still at the Γ point and the band gap is indirect in the mBJ band structure.
If mBJ method indeed gives the right band gap for Fe2 VAl and this compound is an intrinsic
narrow band gap semiconductor then the large gap (0.55 eV) obtained from GGA+U calculation
suggests that the constrained DFT method proposed by Anisimov and Gunnaarsson overestimates
the value of U, at least in these ternary compounds with small band gaps. In their method, Anisimov and Gunnarsson only considered the screening from s and p electrons and neglected the effect
of d electrons themselves. This might underestimate the screening effect on U because small band
gap in this system involves excitations of d-electrons. Suppression of these excitations will tend
to give large U values. Having noted this shortcoming of the present U calculations smaller values of U have been checked. Results show that mBJ result could be reproduced by GGA+U with

55

Figure 3.5: Thermopower (S) and chemical potential (µ ) versus carrier concentration (ne , nh ) in
the absence of U (the pseudo-gap case) for n- and p-type doping in Fe2 VAl. The Fermi level is
chosen to be 0.
UFe =3 eV and UV =1 eV. The similarity in the results obtained from mBJ and GGA+U calculation
is worth a further study.

3.4 Thermoelectric properties
Thermopower S was calculated using Eqns. (2.91) and (2.92) for different doping levels (different
concentrations) and for different temperatures (T ) in the range 100–700 K.

3.4.1

GGA calculation

First let us consider the case when U=0 eV (GGA calculation). This corresponds to the negative
gap or pseudo-gap picture in which Fe2 VAl is a semi-metal with a small but finite DOS at EF .
In the absence of any doping the electron concentration (ne ) and hole concentration (nh ) are the
∫ +∞
1
same. We define the carrier concentration as n = V (N − −∞ D(ε ) f0 (ε )d ε ), where N is the total

56

number of electrons in the system, V is the volume, D(ε ) is the density of states and f0 (ε ) is
the Fermi-Dirac distribution function. The sign of n corresponds to the sign of the charge of
the carriers. ne = −n when n is negative and nh = n when n is positive. With this definition,
carrier concentration is equivalent to the doping level in the experiment. We change n by changing
the chemical potential and see how S changes with n, and T . Fig. 3.5 shows thermopower (S)
and chemical potential (µ ) as a function of n for two different temperatures, 300 K and 700 K.
Although the magnitude of S increases with T , the higher the T , the smoother is the variation of S
with n, caused by enhanced thermal broadening of the Fermi distribution function (see Eqn. 2.92).
The magnitude of S has two peaks for n-type doping, the first one at ne ≈ 5.1 × 1020 cm−3 and
the second at ne ≈ 5.0 × 1021 cm−3 . For p-type doping there is only one peak at and nh ≈ 2.3 ×
1021 cm−3 . These three concentrations correspond to the chemical potential of 0.13, 0.43, and
-0.18 eV, respectively.
We can understand these peaks in magnitude of S by analysing the band structure of Fe2 VAl.
Let’s first understand the n-type thermopower as a function of n. Without doping the chemical
potential lies in the overlap region between conduction band and valence band (see Fig. 3.4-(a)).
In this region the hole and electron contribution in S cancel each other resulting in very small magS σ +S σ
nitude of S following the equation S = e σe +σh h , where Se , σe and Sh , σh are thermopowers and
e h

electrical conductivities associated with the electrons and holes, respectively. When ne increase,
the chemical potential increases and eventually gets out of the overlap region at ∼0.13 eV. As a
result, the hole-electron cancellation is suppressed. The suppression of hole-electron cancellation
gives rise to the first peak in the n-type S. The second peak in n-type S is obtained when n is high
enough for chemical potential to reach the flat band of Fe at 0.43 eV. At this point the magnitude
of S gets a large value due to the rapid change in DOS.
Similarly, one can see that the peak in p-type S is due to both suppression of hole-electron
cancellation effect and the rapid rise in DOS near -0.18 eV. The coincidence of the two events:
getting out of the overlap region and reaching the rapid rise in DOS is the reason why there is only
one peak in the magnitude of p-type S.

57

Figure 3.6: S as functions of temperature obtained using the GGA band structure (U=0 eV, the
pseudo-gap case) at different concentrations: for n-type (closed symbols) and p-type (open
symbols) doping.

It appears from these results that if the band structure of Fe2 VAl is adequately described by
LDA or GGA, i.e. it is a pseudo-gap system, then to have large thermopower values one has to go
to rather high doping levels when the cancellation effects are nearly absent. However, the carrier
concentrations should not be so large that Pauli suppression for degenerate carriers start to reduces
|S|.
Now let us look at thermopower as a function of temperature. Fig. 3.6 plots S versus T for
several carrier concentrations. We will consider n- and p-type dopings separately.
n-type doping: In the temperature range 100–700 K, S has the highest value (magnitude) for
ne ≈ 4 − 6 × 1021 cm−3 . The magnitude of S increases with T , but the increase is rapid at low T
and then slows down at high T . At high T , large number of electron and hole excitations start to
dominate, resulting in a cancellation and reduction in S, as seen in figure 3.6 where T > 700 K).
This electron-hole cancellation is also the reason for small magnitude of S values seen for small n.
p-type doping: One sees behaviour similar to the n-type doping. Again S increases with T and

58

has small values for low concentrations due to electron-hole cancellation. S is largest (100 µ V/K
) for nh ≈ 2×1021 cm−3 . Interestingly, GGA calculation shows that Fe2 VAl is a slightly better
p-type thermoelectric than n-type. For example, at 700 K, the highest value of the magnitude of
S is 110 µ V/K for p-doping whereas it is only 80 µ V/K for n-doping, about 30 % smaller. To
understand this difference we look at the band structure. The electron pockets center around the X
points and are non-degenerate (ignoring spin degeneracy), whereas the hole pockets center around
the Γ point but are three-fold degenerate. Since there are three inequivalent X points in the BZ
and only one Γ point, the difference in the S value cannot be ascribed to the band degeneracy. We
believe it is due to the difference in the effective masses of electrons and holes. The Fe dominated
holes have larger effective mass and lead to larger S values. This property of the Fe2 VAl is different
from that of ZrNiSn,[99] which belongs to the class of so called Half-Heusler compounds. In
the latter compound the electron pockets centered around the X points are heavier than the hole
pockets centered around the Γ point. Consequently, the magnitude of the thermopower is larger
for the electrons in ZrNiSn.[99] For more discussions on this compounds see chapter 5.

3.4.2

GGA+U calculation : the effect of finite band gap

As we have discussed earlier, the effect of incorporating U is to decrease the overlap between the
valence and the conduction bands and eventually open up the gap. For UFe =4 eV and UV =1.5 eV
(values calculated using constrained DFT), the band gap is 0.55 eV. For the calculation of the
transport coefficients we choose the zero of energy at the middle of the gap. In Fig. 3.7 we give
the results of S as a function of T for different values of the concentration. The T -dependence of S
for the case of real gap is much simpler compared to the pseudo-gap case. Since there is very little
electron-hole cancellation, S follows the Pizarenko relation, i.e. the magnitude of S increase with
decreasing carrier concentration. Thus, increasing carrier concentration reduces the magnitude of
S. And for the same reason S does not saturate up to 700 K as the thermal excitations of electrons
and holes are suppressed due to the large band gap.
As mentioned above, in the presence of U, valence band maxima near Γ and X points become

59

Figure 3.7: S as functions of T for GGA+U calculation (UF e=4 eV and UV =1 eV) at different
concentration: for n-type (closed symbols) and p-type (open symbols) doping.
nearly degenerate. Since more hole pockets contribute to transport (both at the Γ and X points)
(see Fig. 3.4(b)) we found that opening up of the gap by U leads to a better p-type thermoelectric.
The value of S for the p-type at concentration of 5 × 1019 cm−3 at 700 K is about 3/2 times larger
than that for the n-type at the same carrier concentration and temperature.

3.4.3

Comparison with experiment and other theoretical results

For comparison, the theoretical results are compared to the experimental results on Fe2 VAl1−x Gex
by Nishino et al. [93] and Fe2 VAl1−x Six by Skoug et al. [19] for x values in the range 0.0–0.2.
Although these two systems are similar and show qualitatively similar behaviours, there are quantitative differences which might help us in testing our theoretical calculations, particularly the
adequacy of the pseudo-gap (or semimetal) picture.
Let us first look at the nominally undoped system (x=0). Both experiments show p-type behaviour in the temperature range 100–400 K. S of nominally undoped system increases slightly
with T and saturates between 300–400 K. The values of S at 300 K, however, differ by nearly a
60

Figure 3.8: Calculated p-type S as a function of n for the cases of: pseudo-gap (GGA calculation)
(solid line) and gapped (GGA+U calculation), Eg =0.55 eV, (dashed line) at 300 K.
factor of 2 between the two measurements (25 and 50 µ V/K obtained by Nishino et al. and by Skoug et al. respectively). One can argue that this difference can be due to the difference in the hole
concentrations and the small values of S can be understood within the pseudo-gap model when
the hole concentration lies between 1020 –1021 cm−3 , because in this model, the electron–hole
cancellation effect makes S to be not only small in magnitude but also a weak function of carrier
concentration over a wide range. (see Fig. 3.8). In contrast, for a finite gap case (especially for a
large gap 0.55 eV as obtained in the GGA+U calculation) one can not understand the peak in S because S follows the Pizarenko relation and changes rapidly as the function of carrier concentration
(see Fig. 3.8). From Hall measurement done by Nishino et al. [95], the hole concentration for the
nominally undoped sample is found to be ≈ 5 × 1020 cm−3 . For this concentration the GGA+U
band structure gives S ≈ 140µ V /K at 300 K, which is about two times larger than the experimental
values (25 and 50 µ V/K), while GGA band structure gives S ≈ 30µ V /K in quite good agreement
with experiment. This result would then favour the pseudo-gap or semimetal picture. Physically,

61

this means that when the hole concentration is not too large, small S values can be obtained near
300 K only when there is significant electron-hole cancellation.
When one replaces Al by Ge, Si, or Sn, the dopants give an extra electron/dopant to the network and when x > 0.03, the alloys show n-type behavior. One significant feature in the n-type
compounds studied by both Nishino and Skoug’s groups is that for the optimum doping (x ≈0.05–
0.06), S shows a minimum between 200–250 K and the value of S near the minimum is about 150 µ V/K[19, 93]. For larger x (presumably for larger carrier concentration), the minimum moves
towards higher T .
Studying S as the function of T corresponding to the pseudo-gap case (U=0) (Fig. 3.6), we
found that that for the entire range of electron concentration ne ≈ 5×1019 −1022 cm−3 , there is no
discernible minimum in S at temperature lower than 500 K and also the magnitude of the S values
are less than 60 µ V/K (two times smaller than the experimental values) in the temperature range
100–700 K. Both these observations are in disagreement with experiments. So neither the pseudogap model nor the model with large band gap (∼0.55 eV) can explain the observed minimum in S
at T ∼ 250K and its magnitude.
In an attempt to see if a finite gap picture could explain the experimental data I have tried using
smaller values of U than that obtained from the constrained DFT calculation. We found that with
very small U (≈1 eV) (and almost zero band gap,Eg ≈ 0.07eV ) we could qualitatively reproduce
the experimental data of S for the n-type system, particularly the minimum at ∼ 250 K (see Fig.
3.9). However in order to get such an agreement with experiment one has to use n ≈ 5×1019 cm−3
which is about one order of magnitude smaller than that obtained from the Hall measurement.[93]
For the p-type (nominally undoped) system with hole concentration of 5 × 1020 cm−3 (obtained
from Hall measurement[93]) the magnitude of S is roughly of the order of the experiment done by
Skoug et al. However, S does not saturate in the range 0–700 K.
The above analysis shows that it is difficult to understand the experimental data for both the
nominally undoped and n-doped cases using either the pseudo-gap or the real gap model. This
conclusion agrees with the work of Bilc and Ghosez [96]who used B1WC hybrid functional

62

Figure 3.9: Experimental thermopower obtained by doping Si on Al site (symbols) (by Skoug
et al. [19]) and theoretical curve (lines) with U=1 eV (Eg = 0.07 eV ). (carrier concentrations are
given in unit of cm−3 )
method [57] and found a band gap of about 0.34 eV which lies between our value of 0.55 eV
and the experimental value of about 0.2 eV. With this value of band gap they could not reproduce
the experimental data. So just scaling the gap to a small value cannot explain the n,T-dependence
of S. In order to understand the difference between experiment and theory, one has to overcome
several weaknesses of the present calculations.
First, the effect of inter-site Coulomb repulsion between d-electrons at neighbouring V- and
Fe-sites has not been included within the present GGA+U approach. Also, GGA+U is essentially
a mean-field approximation, and does not include electron-hole excitonic correlations which may
alter the electronic states near the Fermi energy in the pseudo-gap or small gap regime. Second,
we have used a rigid band model to calculate S. It is possible that the dopants distort the electronic
structure near the Fermi energy and hence the rigid band model is not valid any more. [100] Third,
there may be defects (Fe-V antisite defects, vacancies, and interstitials) which alter the states near
the Fermi energy, giving rise to not only new (and perhaps smaller than the band gap) excitation

63

energy scales but also to localized states that act as carrier traps.

3.5 Summary
In this chapter I have discussed how the band structure of Fe2 VAl changes when one takes into
consideration the effect of intra-site Coulomb repulsion between d-electrons at the Fe and V sites
beyond GGA using GGA+U approximation. In Heusler compounds, which contain two different
types of transition metal atoms, the effect of U on the band gap depends sensitively on which one
of the two transition metal sites one is dealing with. For example, even small values of UV (∼1
eV) at the V site can open up the gap because the lowest conduction band is predominantly V
d-character and gets pushed upward in energy strongly by UV . In contrast, small values of UFe
(∼1 eV) at the Fe site do not open up the gap, it rather changes the band structure slightly. This
is because the top of the valence band is strongly hybridized and is shifted downwards in energy
only slightly. Thus one needs a larger value of UFe (∼2 eV) at the Fe site to open the band gap
when UV = 0.
The values of U calculated using constrained density functional theory (UFe =4 eV for Fe and
UV =1.5 eV for V) are much smaller than the atomic values of U (∼20 eV) indicating a strong
screening effect. This is not surprising and is seen in many transition metal compounds. What
surprising is that the values which have been calculated for Fe is smaller than that one finds in Fe
metal (∼6 eV) since the screening should be stronger in the metal.
Using the values UFe =4.0 eV and UV =1.5 eV give a gap of 0.55 eV in Fe2 VAl, considerably
larger than that obtained from transport measurements and NMR (eg ≈ 0.2 eV ). In this case, S
shows characteristics of usual narrow band-gap semiconductor. Both GGA and GGA+U calculations suggest that Fe2 VAl is a better p-type thermoelectric either due to larger hole effective mass
(in GGA) or due to large degeneracy (in GGA+U).
Calculations, however, point out that the present models (both pseudo and real intrinsic gap
cases), cannot explain all the experimental data on S consistently for both n- and p-doped cases.
Band structure calculated using modified Becke-Johnson potential tend to give small band gaps
64

(0.22 eV). But these results are still not enough to explain the thermopower data. Several limitations of the present models, such as the electron-hole (excitonic) correlation effects, defect induced
changes in band structure, etc. , were pointed out and further work should be done.

65

CHAPTER 4
DIAMOND-LIKE TERNARY COMPOUNDS

4.1 Introduction
Among approaches to improve ZT , reducing lattice thermal conductivity is considerably less complicated (discussed in Sec. 1.1), in the sense that it does not involve other parameters (namely, S,

σ , and κelect ). This direction has been quite successful, bringing ZT to the current state of 2–
3.[14, 101] One way to reduce thermal conductivity is to increase the number of atoms in a unit
cell, following the work of Slack and his colleagues[26] which suggested that thermal conductivity
is inversely proportional to the number of ions per unit cell

κlatt =

¯ 3
BM δ θD
n2/3 T γ 2

,

(4.1)

¯
where B is a constant, M is the average atomic mass, δ is the average volume per atom, θD is the
average Debye temperature, T is the temperature, n is the number of atoms per unit cell and γ is
the Grüneisen constant.
For diamond-like (tetrahedrally bonded) compounds, Grimm and Sommerfeld’s rule[102] provides a simple way to search for compounds with greater number of atoms per unit cell. The rule
predicts that systems with the average of four valence electrons per atom will have diamond-like
structure where all atoms are tetrahedrally coordinated. Applying this rule, one can start from the
multiple unit cells of a simple system (say, Silicon) and replace one constituent with two or more
of different types, keeping the same average number of valence electrons. For example, we can
double the unit cell of Si and then replace two Si by one Zn and one Se to get sphalerite ZnSe
(zincblende/wurtzite structure). Similarly, starting from ZnSe with zincblend structure, one can
double the unit cell to (ZnSe)2 , keep the Se sublattice intact, and replaces two Zn (divalent – II)
atoms by Cu (monovalent – I) and In (trivalent – III) to get the chalcopyrite structure. Another
example is to triple the ZnSe unit cell and replace the three divalent anions by two monovalent and
66

one tetravalent cations, Cu2 GeSe3 (in general I2 -IV-VI3 ).
While most of diamond-like compounds have been extensively studied for their photovoltaic
application due to this wide band gaps, such as CuInSe2 [103] which has a band gap of 1.032 eV,[104]
recently several Cu-Sb-Se(S) based ternary compounds have been found to exhibit promising thermoelectric properties. [27, 28, 105–109] Among these, Cu3 SbSe4 (Se4) is a narrow band gap
semiconductors [109–111] which gives a large p-type thermopower. [106] Even though several
theoretical studies have been done on this compound,[109–111] its electronic and thermoelectric
properties are not well understood. In my work, I will be concerned with Se4 and the related systems in the class of I3 -V-VI4 compounds, where the unit cell is quadruple of ZnSe and 4 divalent
cations are replaced by 3 monovalent (Cu,Ag) cations and one pentavalent (P, As, Sb, Bi) cation.
The I3 -V-VI4 compounds are classified into the family of famatinite (Fa) and enargite (En) which
have tetragonal and orthorhombic structures respectively.
Tetrahedrally coordinated semiconductors have been studied extensively and have played an
important role in our understanding of the relationship between coordination, bonding, and band
gap. This relationship is easy to understand in monoatomic covalent solids C, Si, Ge, Sn. In these
solids spatially directed sp3 hybrids form bonding (B) and antibonding (AB) bands.[112, 113]
Electrons fill up the bonding bands (valence bands) following the simple Lewis octet rule[114,
115], the antibonding bands (condution bands) are empty, there exists a band gap (Eg ) between the
valence band maximum and the conduction band minimum. There is a decrease in the splitting
between the bonding and anti-bonding bands due to increase in the lattice constant as one goes
from C to Sn, with a concurrent reduction in the band widths. The net result is a decrease in Eg
which is nearly zero for Sn. Although this qualitative relationship between bonds and band gap is
useful, the actual value of Eg depends on other subtleties like the details of band dispersion, spin
orbit interaction etc.
The inter-relationship between bonds, bands and band gaps discussed above can be extended
to tetrahedrally bonded III-V and II-VI binary compounds.[112, 113] Here again the octet rule
plays the dominant role. However the splitting between the valence and the conduction bands and

67

therefore the band gap depends on other parameters, namely the differences in electron affinity
(EA) and ionization energy (IE) of the two components. The bonds now pick up both covalent and
ionic characters. The classic works of Pauling[112] and Phillips[113] in relating the ionicity and
covalency of a particular bond to the above parameters (Pauling) and to the dielectric properties
(Phillips) have given a fundamental understanding of the bonding and structure of these binary
compounds. However, these theories do not explicitly address the relationship between band gaps
and the nature of bonds. (Phillips ionicity scale does depend indirectly on the band gap through the
dielectric constant). The general trend seen in covalent semiconductors, i.e. decrease in Eg with increasing lattice constant (bond length) is also seen in binary tetrahedrally bonded semiconductors.
For example Eg = 6.2 eV in BN (lattice constant a=0.361 nm) and 0.23 eV in InSb (a=0.648 nm).
Although this tendency is quite general the actual values of the band gap depends on the details
of band dispersion. Progress in the in electronic structure theories (ab-initio DFT both local and
non-local, GW, and others) have helped us in addressing the relationship between bonding, band
structure and band gaps in these compounds.[38]
Extending the above ideas further to tetrahedrally coordinated ternary compounds, in this case,
Fa and En compounds, becomes challenging because of the competition between the natures of
different bonds. Both these compounds contain group V ions which can exhibit multiple valency.
According to a simple valence count, with Cu (I) being 1+ and Se (VI) 2-, for the charge to be
compensated V should be pentavalent (5+) to explain the semiconducting properties of these compounds. However, GGA calculations for Se4 (Wang et al. [116]) suggests that Se4 is metallic
while experiments reported a band gap of 0.1–0.4 eV.[106, 110, 111] A recent calculation employing the hybrid functional scheme developed by Heyd, Scuseria and Ernzerhof (HSE), which
includes non-local exchange interaction, gives a band gap of 0.4 eV[109]. Nevertheless, in all the
previous studies of Se4, the questions of Sb valency, the precise nature of the bands near the gap
region and the physical origin of the gap formation in Se4 have not been addressed.
This chapter will try to answer some of the questions mentioned above. It is arranged as following. In section 4.2 the crystal structures of Fa and En are briefly reviewed. I will discuss the

68

band gap formation in Se4 in 4.3 by carrying out extensive ab-initio electronic structure calculations using several local/semilocal approximations such as local density approximation (LDA),
GGA,[47] GGA+U[66] and recently developed modified Becke-Johnson (m-BJ)[52] exchange potential. In addition, I have also used non-local approximations to exchange correlation potential
(HSE06)[59–61] to calculate the electronic structure. Then in section 4.4, I show that the physics
of band gap formation is generic for the whole family of I3 -V-VI4 compounds and study how the
band gap varies for different compounds. Then the temperature and carrier concentration dependence of thermopower in Se4 are presented and how they depend on the details of the electronic
structure is discussed. Thermoelectric properties of Se4 and summary of the chapter are given in
section 4.5 and 4.6 respectively. The results presented in this chapter are partially published in
references [36] and [37].

4.2 Crystal structures of Famatinite and Enargite
Systems are considered in this chapter belong to the family of I3 -V-VI4 compounds where (I) is
Cu or Ag, (V) is P, As or Sb, and (VI) is Se or Te. While most Ag-compounds are not found in
the literature, Cu-compounds usually crystallize in the Famatinite structure (Fa) and some in the
Enargite structure (En). As discussed in the previous section, Famatinite and Enargite structures
may be considered as derivatives of the well known Zinc blende and Wurtzite structures respectively. As shown in Fig. 4.1, Fa and En can be obtained by doubling the ZnS unit cell (Zincblende
or Wurtzite respectively) and replacing four Zn atoms by three group I element (Ag or Cu) and one
group V element (P, As, Sb). This keeps the valence electron count intact (4 divalent Zn vis-à-vis
3 monovalent Ag/Cu and 1 nominally pentavalent P/As/Sb). In our calculations, due to different
space groups (body centered tetragonal for Fa and simple orthorhombic for En), the Fa unit cell
contains one formula unit while En unit cell has two formula units. This difference gives rise to
different Brillouin zones and a difference in the number of bands for a given wave vector k, i.e.
En has twice as many bands as Fa.
In both the structures, each atom is tetrahedrally coordinated to four nearest neighbors, in which
69

Figure 4.1: Crystal structure of Famatinite (Fa) and Enargite (En) structure are represented as
double unit cell of the Zinc blend (top panel) and the Wurtzite (bottom panel). Green represents S
(Se, Te), gray represents Zn which are replaced by Copper/Silver (blue) and red represents Zn
which are replaced by Sb (P, As, Bi).

70

Table 4.1: Summary of crystal structures of Cu3 AB4
B\A
S

Se

Te

P
En
d = 2.06exp. [117,
118]
2.11GGA
En
d
= 2.29 −
exp. [117]
2.36
2.29
−
2.32GGA
×

As
Fa/En(Tc
=
580)[111]
d = 2.20exp. [119]
2.31GGA
Fa
d = 2.50GGA

Sb
Fa
d = 2.38exp. [118]
2.48GGA

Bi
Fa
d = 2.65GGA

Fa
d = 2.54exp. [120]
2.65GGA

Fa
d = 2.78GGA

?

?

×

En = Enargite, Fa = Famatinite, × = does not exist, ? = unknown, d = nearest neighbor
distance between atoms A and B (unit Angstrom). GGA indicates that parameters were relaxed
using GGA approximation in the present work.
I and V atoms bond to four VI atoms while VI atoms bond to three I and one V. Table 4.1 presents
a summary of crystal types of I3 -V-VI4 compounds and the bond lengths between V and VI. As
the constituent atoms go from smaller to larger radii, the bond lengths increase accordingly. We
will show later that the change in bond lengths and the band gaps are intimately related.

4.3 Electronic structure of Cu3 SbSe4 : band gap formation and the role of
Sb lone pairs and non-local exchange effect
This section will try to resolve the puzzle about the band gap in Se4 mentioned in section 4.1, by
comparing electronic structure obtained using different approximations including GGA, GGA+U
and HSE06. In order to understand how the structural property affects the electronic structure of
Se4, we also carry out calculations with Cu3 SbSe3 (hereinafter, it will be referred to as Se3) which
has one Se less than Se4, and therefore, does not have the tetrahedrally coordinated structure.
Let us first briefly review the experimentally known atomic structures of these two compounds.
Below 700 K Se4 has the Famatinite lattice structure, crystallizing in base centered tetragonal
¯
structure (space group I42m)[111] (Fig. 4.2(left)). At room temperature, Se4 has lattice constants
of a = 5.661 Å and c = 11.280 Å. [120] Each Cu and Sb is coordinated tetrahedrally to Se while

71

Figure 4.2: Crystal structure of Cu3 SbSe4 (left) and Cu3 SbSe3 (right), Cu is blue ball, Sb is
brown, and Se is green.

each Se is coordinated to three Cu atoms and one Sb atom. Unlike Se4, Se3 has a more complicated
orthorhombic structure (space group Pnma) with lattice constants a = 7.986 Å, b = 10.614 Å, and
c = 6.837 Å at room temperature[121] (Fig 4.2(right)). The local coordination and bond lengths in
this compound are very different from those in Se4. Sb has three Se and one Cu nearest neighbors
(nns), whereas Cu has four Se and one Cu nns. The fact that Sb has an asymmetric network
facilitates the existence of a dynamically active 5s2 lone pair of electrons in Se3. This lone pair of
electrons, however, do not actively involve in band gap formation of Se3. This structural distinction
is believed to be the origin of the difference in their physical properties (the physics of band gap
formation and the unusual lattice anharmonicity of Se3[28]) in spite of similar compositions with
identical elements.

4.3.1

GGA calculation and indication of the existence of Sb lone pair electrons

The GGA-PBE band structures for Se4 and Se3 are shown in Fig. 4.3(a) and 4.3(b) (left panels),
respectively. The present GGA calculation for Se3 agrees reasonably well with the recently pub72

Figure 4.3: Band structures of Cu3 SbSe4 (a) and Cu3 SbSe3 (b) using different methods. In (a),
the bold (red) line is “BOI”, the doted lines present the HSE band structure using GGA position.
UCu−d = 10 eV was used in GGA+U calculation. The energy is given in absolute scale and the
Fermi energy is denoted by EF .

73

lished LDA calculations by Sevik and Ça˘ in (SC) [108], showing that Se3 is a narrow band gap
g
semiconductor with an indirect band gap. However a larger band gap was found, 0.49 eV compared to 0.24 eV by SC [108]. There are several reasons that may account for this difference. First,
SC used LDA which usually gives smaller lattice constant (a=7.924 Å) compared to our GGA
calculation (a=8.075 Å). Second, SC define the band gap using the DOS which is very sensitive to
the smearing method used in computing the DOS, whereas we obtain the band gap using the band
dispersion.
In contrast to Se3, GGA calculation reveals that Se4 is not a semiconductor, disagreeing with
the experimental measurements. There is a narrow band, with a width of ∼1 eV (marked as bold
and red in Fig. 4.3, near the top of the valence band, which is partially filled. We refer to this band
as the “band of interest” or “BOI”. At the Γ point, “BOI” is degenerate with the top two valence
bands. At the X point, it gets lower and overlaps with several valence bands. By exploring the
physical origin of “BOI” and its sensitivity to different approximations, it reveals that the origin of
the “BOI” is intimately related to the question of Sb valency in Se4.
It is quite easy to understand the semiconducting ground state of Se3 using the simple ionic
picture and the octet rule. In Se3, Cu takes the 1+ ionized state, Se is 2-, and Sb is 3+ (trivalent).
With this configuration, all the octets are filled, making Se3 a semiconductor. Using the similar
argument, but with pentavalent Sb (5+) for Se4, one would come to the same conclusion. However,
GGA band structure does not seem to support this simple picture, since there is a nondegenerate
conduction band (BOI) overlapping with the valence bands, making Cu3 SbSe4 a semimetal or a
pseudo-gap system.
In order to see the similarities and the differences between Se4 and Se3 I plot the total density
of states (DOS) for both the compounds (Fig. 4.4(a), top panel) and the projected DOS (PDOS) of
Se4 (Figure 4.4(a), bottom panel). It can be seen that the only significant difference between them
is the "BOI" state near the Fermi level (EF ) in Se4. In both compounds, the Sb 5s bands are quite
narrow (∼0.5 eV) and located approximately at -9.5 eV, while the Sb p bands are above the Fermi
level. This suggests that Sb is nominally 3+ in both the systems. In addition, Cu s is empty. As a

74

Figure 4.4: (a) Total DOS of Se4 and Se3 (top panel) and the projected DOS of Se4 (bottom
panel), along with isosurface plot of charge density associated with (b) the lowest band at
-14.5 eV, (c) the Sb-s band at -9.5 eV, and (d) the “BOI”. For clarity, we only show SbSe4
tetrahedron with Sb at the center. The DOS are given in unit of number of states per eV per
formula unit.

75

consequence, Se-p is fully filled in Se3 making it a conventional semiconductor. In contrast, for
Se4, Se p bands are two electrons deficient from being completely filled (corresponding to nearly
unoccupied “BOI”), leading to a semi-metal. The DOS of Se4 near the Fermi energy indeed shows
a pseudogap structure (Fig. 4.4 (a) inset).
To further understand the nature of the bands in Se4, we look at the partial charge densities
associated with three narrow bands, one at -14.5 eV, one at -9.5 eV, and the “BOI” as shown in
Figs. 4.4(b), 4.4(c), and 4.4(d), respectively. At -14.5 eV, there is a bonding state of Se s and Sb s
which splits off from the other Se s states by about 1 eV. The state at -9.5 eV is mostly Sb s with
some Se p, resulting in a charge density with asymmetric lobes around Se directed away from Sb
which is different than usual s-p hybridization. The “BOI” band is a mixture of Se p, Sb s and Cu
d.
To explain the origin of the “BOI” a simple scheme is proposed for the mixing of orbitals in Se4
as shown in Fig. 4.5, where it ignores the interaction between Se and its neighboring Cu orbitals
which gives rise to the p-d hybridization with energy in the range of -6 eV to 0. In the final picture
this hybridization plays an important role, particularly when the effect of Coulomb correlation is
included through the parameter U. At -14.5 eV the mixing between Sb 5s and the symmetric
combination of four nearest neighboring Se 4s (ψSe,s = ∑µ =1,4 ϕSeµ ,s ) leads to a bonding (B)
combination ψs,B = ϕSb,s + AψSe,s whose energy (-14.6 eV) is about 1 eV lower than the rest
three Se 4s bands. The charge density associated with ψs,B indeed shows the bond charge between
Sb and neighboring Se atoms. The charge density for the state at -9.5 eV looks rather peculiar
with asymmetric lobes around Se atoms (Figure 4.4(c)). It can then be thought of as a bonding
combination of an antibonding (AB) state of ϕSb,s and ψSe,s , i.e. ψs,AB = ϕSb,s − AψSe,s and the
neighboring Se 4p states (ψSe,p ). This results in a cancellation of Se 4s and 4p amplitudes in the
region between Sb and Se, leading to the peculiar shape of the charge density, seen in Figure 4.4c.
The corresponding antibonding combination has a finite density in the region between Sb and Se,
leading to the “BOI”. The charge density associated with “BOI” shows a mixing of Sb-5s, Se-4p,
Se-4s and Cu-d characters. Thus it is a rather subtle and complex state. As we will see, its position

76

Figure 4.5: Simple bonding-antibonding scheme in Cu3 SbSe4 showing energy levels in the
atomic limit (left column) intermediate bonding-antibonding states (middle column) and
bonding-antibonding states in crystal (right column). (Cu d and Sb-p are not shown for
simplicity). The number of states associated with each level is given in the bracket.

and dispersion are quite sensitive to the approximations made in treating the exchange-correlation
potential. However, the general mixing picture discussed above (over energy range ∼10 eV) is
valid for all the approximations used.

4.3.2

Failure of local exchange approximation – Effect of non-local exchange interaction

In order to see if the semiconducting state can be obtained in Se4 within local approximation,
I have tried a recently proposed semi-local model called modified Becke-Johnson (m-BJ) model
which introduces a correction which depends on the kinetic energy density (∑ |∇ψiσ |2 ) to GGA.
This model has been found to give reasonable gap values in a large number of semiconductors
[52]. Within this model, Se3 is found to be a semiconductor with a gap of about 1 eV compared to
0.49 eV obtained in GGA. However, m-BJ does not open up a gap in Se4.
For Se4 to be a semiconductor, the “BOI” must split off from the rest of the Se p bands. Since

77

(a)

Energy(eV)

-2

(b)

(c)

(d)

-4
BOI

EF
-6

21
3
P
Γ
Sb s
Sb p

X|P

Γ X|P
Cu t2g
Cu eg

Γ

X|P

Γ
Se p

X

Figure 4.6: Band structure of Cu3 SbSe4 in absolute energy scale for the ratios of exact exchange
α = (a) 0, (b) 0.1, (c) 0.2 and (d) 0.3 with fat band representation of orbital character, showing
that the bands with dominating d characters affected strongly; there is an exchange of band
character between BOI and BOI-3 at a critical value of α , αc ≈ 0.2.
it includes a small admixture of copper d-states the effect of intrasite d-d Coulomb repulsion can
have important effects on the position and hybridization of “BOI”. For this, I carried out GGA+U
[65] calculations with different U values (0≤ U ≤20 eV) at the copper site. In Fig. 4.3 (middle
panel) we give the band structure with UCu =10 eV. We find that in Se3, by increasing U from 0 to
10 eV, band gap changes from 0.49 eV to 1.10 eV (similar to what is seen in m-BJ). There is a large
rearrangement of most of the valence bands but the dispersion of the top most valence band and the
lowest conduction band does not change much (which is relevant for thermoelectric properties).
However in Se4, the dispersion of the “BOI” changes drastically and there is a level-reordering in
the band structure near the Fermi level. For example, at the X points, the “BOI” is pushed higher
than the valence bands and there is also a dramatic change in the top four bands near the Fermi
level at the Γ point. However, a gap does not open up unless U is larger than 12 eV. This value is
however unphysically large compared to typical values of U for copper, 5–8 eV.[68, 122–125]
In order to circumvent the problem of using large U within local and semi-local theories, the
effect non-local exchange is studied using HSE06 model (α = 0.25). [59–61] Since Yang et al.
78

[109] have already reported the existence of a band gap in Se4 using this model, we wanted to
investigate in detail how non-local exchange affects different bands near the Fermi level and how
the gap opens up. In Se3, HSE06 increases the band gap from 0.49 in GGA to ∼1.5 eV as expected
(Fig. 4.3(b) right panel). In Se4, even for this model, if we use the atomic positions optimized using
GGA, there is no gap (presented as dotted lines in Fig. 4.3(a)). After carrying out further internal
relaxations within HSE06, we find that, compared to GGA values, the distances between Sb and
its four Se neighbors decrease by 3.4%, the distances between Cu and its neighboring Se increase
by ∼1% (detail is shown in table 4.2)and a direct gap of 0.26 eV is formed at the Γ point between
“BOI” and the rest of the occupied bands (Fig. 4.3(a) right panel). This value of band gap is in
good agreement with experimental values (0.1–0.4 eV)[106, 110, 111] but smaller than the value
of 0.4 eV found by Yang et al. [109]. We believe that the difference between our result and that
of Yang et al. [109] is due to the difference in the lattice parameters used, and also due to the
estimation of the gap from DOS calculation by these authors, rather than using the actual band
dispersion.
To further understand the band gap formation in Se4 the effect of non-local exchange on the
band structure is studied using different values of α and look not only at the “BOI” but also three
other bands right below it. These three bands are denoted as BOI-i, where i=1–3 (labelled "1", "2"
and "3" in Figure 4.6). Fig. 4.6 shows the band structure of Se4 with different values of α showing
how different bands are affected. For each value of α , internal relaxation was performed.
For α =0 (i.e. GGA), at the Γ point, three bands “BOI”, BOI-1 and BOI-2 are degenerate and
have mainly Cu t2g and Se p characters. When one moves away from the Γ point the character of
“BOI” changes to a mixture of Sb s, Se p and Cu eg , whereas the characters of BOI-1 and BOI-2
remain the same as at the Γ point. The character of BOI-3 is the opposite of that of BOI, i.e. a
mixture of Sb s, Se p, and Cu eg at the Γ point which changes to a combination of Se p and Cu t2g
by moving away from it. This inverted-band feature was discussed by Wang et al. [116] to point
out the nontrivial topological property of Se4 which, as shown later, is not valid when HSE06 is
used. There is a band-crossing at the X point where the energy of the “BOI” is lower than the other

79

Table 4.2: Nearest neighbor distances (Å) and band gaps (eV) with different exact exchange
ratios

α
Sb-Se Cu1-Se Cu2-Se Eg (eV)
0
2.65
2.44
2.43
–
0.1
2.59
2.46
2.44
–
0.2
2.55
2.47
2.45
0.08
0.3
2.53
2.48
2.45
0.6

three. With increasing α , the bands containing Cu t2g are pushed down faster. There is a removal
of band-crossing at the X point. When α ≥ 0.2, at the Γ point, the “BOI” switches its character
with BOI-3 and remains nearly at the same position while the other bands keep going down. After
the band gap opens up, “BOI” has consistent character throughout the k points sampled, namely,
it is primarily a mixture of Sb 5s, Se 4p and a small admixture of Cu eg . At this point, a gap
formation is observed at the Γ point. Then the “BOI” becomes the lowest conduction band (LCB)
which is predominantly a mixture of Sb 5s and Se 4p. As mentioned above, the band structure
using HSE06 does not support the inverted-band picture suggested by Wang et al. [116] because
of the switching of characters between BOI and BOI-3 when the gap opens up. It is interesting to
note that for α = αc , one sees a linear dispersion (Dirac cone) involving BOI and BOI-3.
Studying the nearest neighbor distances in Se4 as a function of the strength of non-local exchange reveal that by increasing α the Sb-Se distance gets shorter while Cu-Se distance gets
longer.(Table 4.2) This suggests that the non-local exchange interaction is stronger between Sb
and Se due to large overlap between Sb 5s and Se 4p.
It is interesting that, in Se4, local and semi-local approximations to the exchange interaction
fail to open up the gap. Furthermore, employing non-local exchange only is not sufficient to open
a gap, one has to relax the bonds. These observations strongly suggest that the nearest neighbor
distances between Se and Sb, and between Se and Cu play an important role in opening the gap.
Thus, it is believed that non-local exchange interaction between the 4s and 4p electrons of Se
and the neighboring Sb 5s and Cu 3d electrons are responsible for the gap formation. The band
structure of Se4 obtained using HSE06 can be directly tested experimentally by (1) probing the

80

excitations from occupied Cu d bands to the conduction bands, i.e. by determining the peaks in
the excitation spectrum one can confirm the position of Cu d states which are at about 3–4 eV
below the lowest conduction band; and (2) doing NMR measurement in n-doped sample, which
should show a Knight-shift of s-spin density at the Sb nucleus.

4.3.3

Unphyscal value of U and its meaning

Even though GGA+U fail to predict the semiconducting ground state of Se4 within typical (physical) range of U value (unlike the case of Fe2 VAl discussed in chapter 3), one can get results almost
identical to that obtained in HSE06 within GGA+U (local theory) but using an unphysically large
value of U (∼15 eV) for the Cu-d orbitals. One can see that changing U gives similar effect as
changing α . As seen in Fig. 4.7a, increasing U does two things, it lowers the occupied d bands of
Cu which reduces the hybridization between the Cu d and the rest of the s, p bands (the amount
of d-character in BOI at the Γ point decreases) and eventually for a sufficiently large value of
U (> 10eV), a gap opens up between the BOI and the rest of the occupied valence bands. one
also observes a similar effect in HSE06 calculations as one increases the strength of the non-local
exchange by increasing the value of the mixing parameter α as shown in Fig. 4.6.
One important difference between GGA+U and HSE06 is that the Cu-d bands are pushed down
much further in the former, which can be easily observed in Fig. 4.7. In Fig. 4.7a, as U increases,
the mixing between Cu d and Se p decreases, and as a result, the d character in the conduction
bands near Fermi level gets weaker and eventually disappear at large U. This effect is clearly seen
in Fig. 4.7b which compares band structures obtained using GGA, GGA+U and HSE06. If HSE06
indeed gives correct band structure, one should detect a Cu-d peak ∼ 3 eV bellow the Fermi level
in XPS measurements.
Despite the differences in contributing orbitals, the band structure near the Fermi level are
similar in GGA+U (U=15 eV) and HSE06 (α =0.25). Since physics of materials depend mostly on
the electronic structure in this range, GGA+U (even with unphysical U) calculations can give some
good insight about physical systems of interest without consuming huge computational resources

81

(a)

U=0

Energy(eV)

-2

U=4

U=15

U=20

-4
EF
‘

-6

-8 P G X|P
Sb s
Sb p
(b)

G X|P G X|P
Cu t2g
Cu eg

GGA

-2

Energy(eV)

U=10

G X|P G X
Se p

U=15

HSE06

-4
EF
-6

-8 P
Sb s
Sb p

Γ

X N|P

Γ
Cu t2g
Cu eg

X N|P

Γ

X
Se p

N

Figure 4.7: (a) Effect of U on band structure and (b) comparison between GGA+U and HSE06
band structures. The energy is given in the absolute scale.

82

as required by HSE06.

4.4 Bond and band gap in the a class of I3 -V-VI4 Famatinite and Enargite
compounds
In section 4.3.1, a simplified model is proposed for the bonding and anti-bonding states between
the states of different atoms as shown in Fig. 4.5. In this picture, s-level of Cu and p-state of group
Sb have highest energy and donate all their electrons to the others. The s-states of the Se can be
treated as filled core states and the s-states of the monovalent cation as empty states to start. The
p-states of Se anions and s-state (lone pairs) of Sb strongly interact and form the bonding and
anti-bonding bands, giving rise to an unoccupied band (the so called "BOI") near the Fermi level.
The position of the BOI is very sensitive to the approximation schemes used in calculations of the
electronic structures and the distance between the Se atom and its four Se neighbors. The lone
pair states of Sb are somewhat schizophrenic, on the one hand they bond with the s-states of the
surrounding Se s-cores but also like to bond with one properly symmetrized combination of the p
states of the surrounding Se4 cluster. This latter bonding is so strong that it splits off one band from
the top of the Se p-bands giving rise to a band gap. The d-bands of Cu are quite narrow and lie
below but near the Fermi level. Their precise role in the band structure depends on the particular
type of atoms forming the compound. The aim in this section is to extract some of the generic
(universal) roles of the V lone pairs and to show that the proposed bonding-antibonding picture for
Cu3 SbSe4 is rather universal and applicable to the other tetrahedral compounds containing V-VI4
(SbSe4 ) clusters.

4.4.1

Atomic energy levels of constituents and some schematic studies.

In analyzing bonding and antibonding nature of different bonds and bands, it is helpful to know
the energies of atomic levels of the constituent atoms. Table 4.3 gives the atomic energies of the s
and p valence electrons for different constituents and for Cu the energy of the d-level is also given.
Positions of these atomic levels are shown in Fig. 4.8. It is clear that Cu-d and S/Se-s levels are
83

Table 4.3: Atomic energies of the s and p valence electrons for different constituents and for the
monovalent Cu, the energy of Cu d-level is included
Energy(eV)
Es
Ep
Ed

P
-19.22
-9.54
–

As
Sb
Bi
Cu
S
Se
Te
-18.92 -16.03 -15.19 -7.7
-24.02 -22.86 -19.12
-8.98
-8.14
-7.79 –
-11.6 -10.68 -9.54
–
–
–
-20.26 –
–
–

Figure 4.8: Visualization of the atomic levels of constituent atoms.

the lowest which makes S/Se-s states stable and occupied while Cu-d are near the Fermi level as
usual. The energies of the s and p levels of group V atoms increase in going from P to Bi.
In Se4, the band structure was found to be affected significantly by the bond length between
Sb and Se. Since nearest neighbor distances increase from P to Bi (table 4.1), one should expect
competing effect of the shrinking of bonding-antibonding separation and the narrowing of the
bandwidth in forming the band gaps. This competing relation is more complicated due to the
differences in energies of different atomic levels. Furthermore, the role of d electrons of Cu also
have to be considered. In order to isolate the effect of d electrons, one can replace Cu by an alkali
atom and study the difference in electronic structure.
84

4.4.2

Results obtained using local approximations

To understand the crucial role played by the lone pairs (ns2 ) of the group V atoms in the bondingantibonding scheme and whether the underlying mechanism for the formation of band gap proposed in Cu3 SbSe4 can be applied to other members of the I3 -V-VI4 family, firstly, GGA calculations were carried out for systems in this family, Cu3 (P,As,Sb,Bi)(S,Se)4 . Fig. 4.9 shows the band
structures of these compounds, focusing near the Fermi level, in which the BOI are marked as
red thick lines. The GGA calculations show that all the compounds have similar features in their
band structures which resemble the simple picture of the bonding-antibonding scheme discussed
earlier. The difference between En and Fa band structures is that for each k the former has twice
as many number of bands as the latter. As a result En has two BOIs. All the compounds have
Cu-d states occupied right below the Fermi level and the lowest partially unoccupied band is the
BOI. The Cu-d states mix with the Se-p states. These systems like Cu3 SbSe4 are semimetals with
pseudogaps near the Fermi energy. For these systems, the Cu-d states mix with the BOI at the Γ
point of the Brillouin zone, giving rise to a three-fold degenerate band, while the next band has
Sb-s and Se-p characters as the BOI. The only exception is Cu3 PS4 where there is a band gap (of
∼ 0.5 eV) between BOI and the rest of the bands below, which are occupied; it is predicted to be a
semiconductor even with GGA level of approximation. The reason which makes Cu3 PS4 unique
is that it has the shortest nearest neighbor distance as listed in Table 4.1. The overlap between BOI
and valence bands increases in going from P to Bi as well as from S to Se.
As discussed in section 4.3.3, since we are dealing with a large number of systems, and due to
limited computational resources, GGA+U (with effective U=15 eV), instead of HSE06, is a good
start to understand the band structures near the Fermi level. Also GGA+U will be much easier to
implement when one studies the nature of defect states in these compounds because in this case
one will have to use large supercells with number of atoms ∼64 or larger. Of course whenever we
have some questions for a particular system one can go back to HSE06. In Fig. 4.10, GGA+U band
structures of (a) Cu3 PS4 , (b) Cu3 PSe4 , (c) Cu3 AsS4 , (d) Cu3 AsSe4 , (e) Cu3 SbS4 , (f) Cu3 SbSe4 ,
(g) Cu3 BiS4 and (h) Cu3 BiSe4 are given. For most of compounds, a band gap between BOI and
85

Energy (eV)

3

(a)

3

0

0

Energy (eV)

-3 Z Γ X S Y Γ
3 (c)

T R

Γ

U

0

-3 M

(b)

-3 Z Γ X S Y Γ
3(d)

Energy (eV)

Γ

U

0

Γ

N

P

Γ

X

N

-3 M

3 (e)

0

Γ

N

P

Γ

X

N

Γ

N

P

Γ

X

N

Γ

N

P

Γ

X

N

3 (f)

0

-3 M

Γ

N

P

Γ

X

N

-3 M

3(g)
Energy (eV)

T R

3(h)

0

0

-3 M

Γ

N

P

Γ

X

N

-3 M

Figure 4.9: GGA bandstructure of (a) Cu3 PS4 (b) Cu3 PSe4 (c) Cu3 AsS4 (d) Cu3 AsSe4 (e)
Cu3 SbS4 (f) Cu3 SbSe4 (g) Cu3 BiS4 and (h) Cu3 BiSe4 . Where the first two compounds are En
and the rest compounds are Fa.

86

Energy (eV)

3

(a)

3

0

0

Energy (eV)

-3 Z Γ X S Y
3 (c)

Γ T R

Γ

U

0

-3 M

(b)

-3 Z Γ X S Y
3(d)

Energy (eV)

Γ

U

0

Γ

N

P

Γ

X

N

-3 M

3 (e)

0

Γ

N

P

Γ

X

N

Γ

N

P

Γ

X

N

Γ

N

P

Γ

X

N

3 (f)

0

-3 M

Γ

N

P

Γ

X

N

-3 M

3(g)
Energy (eV)

Γ T R

3(h)

0

0

-3 M

Γ

N

P

Γ

X

N

-3 M

Figure 4.10: GGA+U bandstructure of (a) Cu3 PS4 (b) Cu3 PSe4 (c) Cu3 AsS4 (d) Cu3 AsSe4 (e)
Cu3 SbS4 (f) Cu3 SbSe4 (g) Cu3 BiS4 and (h) Cu3 BiSe4 . Where the first two compounds are En
and the rest compounds are Fa.

87

2.5
Cu3PS4(En)

2

Cu2CdGeS4
Cu3PSe4(En)

1.5
Cu3AsS4(Fa)

Eg (eV)

1

Cu3AsS4(En)
Cu3AsSe4(Fa)

0.5

CuInS2
Cu3SbS4(Fa)
Cu2SnS3
CuInSe2

CuInTe2

Cu2SnSe3
Cu3SbSe4(Fa)

0
Cu3BiS4(Fa-En)

Cu2SnTe3

Cu3BiSe4(Fa-En)

linear fit
S compounds
Se compounds
Cu3SbS4 with different a
Cu3SbSe4 with different a
other tetrahedral compounds

-0.5
-1
-1.5
2

2.2

2.4
2.6
V-VI distance (Angstrom)

2.8

3

Figure 4.11: Relation between band gap and the distance between V and VI elements.

the occupied valence bands opens up except in Bi compounds. For those compounds where band
gap does open up, the highest occupied bands are multiply degenerate, which makes these systems
excellent p-type thermoelectrics, just like Cu3 SbSe4 which will be discussed later in this chapter.
We also explored how the band gap and the width of the BOI (now the lowest conduction band)
depends on the U parameter. The larger is the U value, the larger is the band gap, and the narrower
is the BOI.
To see how the band structures depend on the crystal structure, for some of the compounds,
calculations were done for both Fa and En structures. Fig. 4.11 plots the magnitude of the band
gap as a function of average nearest neighbor (nn) distance between V and VI elements for Fa and
En. For comparison data of the other ternary compounds involving lone-pair atoms is also plotted
in Fig. 4.11.

88

Table 4.4: Effect of different V element substitutions on the band gap of Cu3 SbSe4 with and
without ionic relaxation
Element
P
As
Sb
Bi

w/ relaxation
dV −V I (Å) Eg (eV)
2.29
0.92
2.42
0.29
2.56
0.28
2.68
0

w/o relaxation
0
0
0.28
0

Let us first analyze the trend in I3 -V-VI4 compounds. One can see that irrespective of the crystal
structure (En or Fa), the band gap decreases with the increase of the V-VI distance (dV −V I ). The
band gap decreases almost linearly with the dV −V I up to d ≈ 2.6Å as represented by the red solid
line in Fig. 4.11. With larger dV −V I , the gap disappears. With other constituents being the same,
compounds with S have smaller dV −V I than those with Se, and thus the band gaps are larger in the
former. However, for the same dV −V I , Se-compounds have larger band gap. For instance, the band
gap of Cu3 SbSe4 is ∼0.4 eV while Cu3 BiS4 is gapless. One exception is the case of Cu3 SbS4 and
Cu3 AsSe4 . Both of them have dV −V I ≈ 2.4Å but the former has a larger band gap.
To explain these features, one may attempt to look at the atomic energy levels of the constituents (given in table 4.3 and Fig. 4.8). The difference between S and Se is that the atomic levels
of the latter is higher than those of the former, the energy difference for s level is ∼1.16 eV and for
p is ∼0.92 eV. These energy differences may account for the band gap difference between S- and
Se-compounds. On the other hand, when going from As to Sb, there is a big jump in the atomic
levels (Fig. 4.8). This change in the atomic levels between As and Sb may compensate for the
difference in the atomic levels between S and Se, making the band gap of Cu3 AsSe4 smaller than
that of Cu3 SbS4 .
To validate the argument above, dV −V I was varied in Cu3 SbS4 and Cu3 SbSe4 by changing
the lattice parameter, keeping the relative ionic positions intact (no ionic relaxation). The results
are presented in Fig. 4.11 as the purple dotted (Cu3 SbS4 ) and cyan dot-dashed (Cu3 SbSe4 ) curves.
Interestingly, for dV −V I > 2.3Å, two curves are almost identical and follow the general trend which

89

we have discussed above; band gap decreases when the dV −V I increases. This result, however,
does not appear to support the relation between atomic levels and band gaps which is argued above
in which one would expect the difference between two curves due to the difference in atomic level
between S and Se. The two curves only differ from each other for small distances where they split
apart by about 1 eV and the trend turns over; the band gap now decreases with decreasing distance.
This change indicates that for short distances, the increasing of the band width overcomes the
splitting of bonding-antibonding band, reducing the band gap. The turning points are different in
two curves, it occurs at larger distance in Se-compound.
In an attempt to resolve this puzzle, the band structures of Cu3 XSe4 , where X is P, As, Sb or
Bi, are calculated using the same lattice parameter as that of Cu3 SbSe4 for two cases: (1) with the
same positions of atoms (hence, the same dV −V I ) and (2) with the ionic positions relaxed. The
results are shown in table 4.4. It is interesting to note that, without relaxation, only Sb-compound
has a positive band gap whereas other compounds have overlap between BOI and the valence
bands. However, when the atomic positions are allowed to relax, we go back to the general trend
where the dV −V I increases as one goes from P to Bi, and as a result, the band gap decreases
accordingly. Thus, different constituents affect the local geometry of tetrahedrally coordinated
compounds differently, giving different values of the band gaps.
Similar effects are also observed in other tetrahedrally coordinated compounds such as I-IIIVI2 , I2 -IV-VI3 . The values of the band gaps (within GGA+U, with U=15 eV) for some systems
are shown as green triangle in Fig. 4.11. These materials, in general, have larger band gaps than
the I3 -V-VI4 compounds.
To verify that the underlying physics of the band-gap formation is dominated by the interaction
between the lone pairs of V and the p-states of VI and less so by the hybridization with the d-states
of Cu, calculations were done for an artificial compound Na3 SbSe4 . The structural parameters,
since it is, to the best of my knowledge, not known in the literature, are initially chosen to be
the same as Cu3 SbSe4 but then allowed to fully relax. Fig. 4.12 gives the band structures of
Na3 SbSe4 obtained using GGA (U=0) and HSE06 calculations. GGA already opens up a band

90

(a)
5

(b)

Energy (eV)

0
-5

-10
-15 P

Γ
Sb s

X

N P

Sb p

Se s

Γ

X

N

Se p

Figure 4.12: Band structure of Na3SbSe4 (Fa structure) obtained using (a) GGA and (b) HSE06
with fatband representation showing atomic orbital associated with energy levels.
gap of ∼1.5 eV similar to Cu3 PS4 and HSE06 increases the band gap to ∼2.5 eV. The bands are
rather narrow compared to the Cu compounds because the lattice constant for the Na compound
is 6.29 Å compared to 5.73 Å for the Cu compound. Also the distance between Sb and nn Se
decreases from 2.65 Å in the Cu compound to 2.49 Å in the Na compound. This leads to a stronger
interaction between V and VI atoms leading to a larger band gap in Na3 SbSe4 . Another interesting
feature shows up, clearly due to the absence of d-levels, other bands split off from the Se-p bands
but pushed below. This band has strong Sb-p character indicating a strong mixing between the
Sb-p conduction bands and the Se-p valence bands. A closer examination of the band structure of
the Cu compounds also shows this feature.

91

Energy (eV)

3

(a)

3

EF

0

Energy (eV)

-3 Z Γ X S Y Γ
(c)
3

Energy (eV)

Γ

U

0

-3 P
(e)
3

Γ

X

N

-3 Z Γ X S Y Γ
(d)
3

T R

Γ

U

-3 P
(f)
3

Γ

X

N

Γ

X

N

Γ

X

N

0

Γ

X

N

0

-3 P

0

0

0

-3 P
(g)
3
Energy (eV)

T R

(b)

-3 P
(h)
3

0

Γ

X

N

-3 P

Figure 4.13: HSE06 bandstructure of (a) Cu3 PS4 (b) Cu3 PSe4 (c) Cu3 AsS4 (d) Cu3 AsSe4 (e)
Cu3 SbS4 (f) Cu3 SbSe4 (g) Cu3 BiS4 and (h) Cu3 BiSe4 . Where the first two compounds are En
and the rest compounds are Fa.

92

2.5
2

Se compounds

1.5

Eg (eV)

1

HSE06
S compounds

Cu3PS4(En)

Cu3PSe4(En)
Cu3AsS4(Fa)
Cu3SbS4(Fa)
Cu3AsS4(En)
Cu3AsSe4(Fa)
Cu3BiS4(Fa)

0.5
0

Cu3SbSe4(Fa)
Cu3BiSe4(Fa/En)

Cu3BiS4(En)
LDA+U

-0.5

linear fit
S compounds
Se compounds
Cu3SbS4 with different a
Cu3SbSe4 with different a

-1
-1.5
2

2.2

2.4
2.6
V-VI distance (Angstrom)

2.8

3

Figure 4.14: Band gap as a functions of V-VI distance for Fa and En: a comparison between
GGA+U and HSE
4.4.3

Effect of non-local exchange

As discussed above, even though GGA+U can give good understanding of the electronic structure
of Fa and En in the vicinity of the Fermi level, the non-local exchange interaction is important, and
it is always a good idea to check GGA+U result using hybrid functional calculation. Especially,
I would like to address the band structure of the two compounds containing Bi, Cu3 BiS4 and
Cu3 BiSe4 . Both GGA (Fig. 4.9g,h) and GGA+U (Fig.4.10g,h) did not give a gap. The HSE band
structures of Fa and En are given in Fig. 4.13. Indeed, the results obtained from HSE calculations
confirm the features of band structures produced by GGA+U. For the two Bi-compounds, there is
no gap opening and both these compounds are semimetals and the band structures are similar to
those obtained from GGA+U.

93

600

Skoug et al
x=0
x=0.01
x=0.02
x=0.03
Theory
(1018cm-3)
0.5
2
10
50
100
150
500

S (µV/K)

400
200
0
-200
100

300
500
Temperature (K)

700

Figure 4.15: Theoretical seebeck coefficient (negative means n-type, positive means p-type) of
Se4 using Boltzmann’s equation for different carrier concentration (1018 cm−3 ) together with
experimental data reported by Skoug et al. [106] in the similar range of doping levels. x is the
excessed holes per unit cell, x = 0.01 corresponds to hole concentration of 50 × 1018 cm−3 .
Indeed, figure 4.14 show that the relation between band gap and V-VI bond length from HSE
calculation is almost identical to that which is discussed for GGA+U previously. This agreement,
again, confirms that one can use GGA+U with an effectively chosen U to obtain reasonable results
while using much less computational resources than that needed for HSE calculations.

4.5 Thermopower calculation of Cu3 SbSe4
Once the electronic structure is understood, the thermoelectric properties, particularly thermopower,
of materials can be easily obtained using Boltzman’s transport equation. And since the band structures are similar between members of Fa and En, with the exception of the band gaps, one can
study one system and extrapolate the results to the other. For example, a compound with a larger
band gap than that of the studied one would be expected to have thermopower saturating at a higher
temperature since it requires higher higher energy for the cancellation effect of excited electron94

hole to occur. By looking at the band structure, one can also deduce that Fa and En would be better
p-type thermoelectrics since there is a three-fold degenerate band with different effective masses
at the top of the valence band. It is indeed the case as will be shown latter on.
We choose to study Se4 compound which is extensively studied in experiment with which we
will compare our calculated results. With the band structure obtained from HSE06 calculations
using α = 0.25 and 21 × 21 × 21 k-mesh, we calculate thermopower (S) as a function of temperature (T ) for different electron (n) and hole (p) concentrations in the range 0.5-500 (in the unit of
1018 /cm3 ) (Fig. 4.15).
In the p-doped Se4, at 300 K, S is about 600 µ V/K for hole concentration p=0.5 and decreases
to ∼100 for p=150, where concentration is given in unit of 1018 cm−3 . This behavior agrees with
experiment by Skoug et al. [106] For fixed p, the present calculation shows reasonable agreement
with experiment by Skoug et al. [106] in terms of both the magnitude and the temperature dependence but not with Yang et al. [109] (Figure 4.15). For a similar doping level Yang et al. get
a much lower thermopower than Skoug et al. and our calculated values. This suggests that their
samples have, perhaps, larger hole concentrations due to defects (most likely Se vacancies).
S increases with T up to 700 K for high concentrations (p=50, 100 and 150) which agree with
doping level of x=0.01, 0.02 and 0.03 of Skoug et al. [106]. However, for p≈0.5, which coresponds to nominally undoped sample, Skoug et al. [106] found that S peaked around 500 µ V/K at
300 K and decreased slowly with T reaching a value of ∼300 µ V/K at 500 K, while one finds in
theoretical results a dramatic drop from ∼600 µ V/K at 300 K to ∼80 µ V/K at 500 K as a result of
electron-hole excitations. This discrepancy might be caused by various reasons. One of them may
be the pinning of the chemical potential around its 300-K value due to the creation or stabilization
of defects. This assumption is consistent with the hole-concentration measurement showing an
increase of p with T . [106]
Comparing Se4 and Se3 (done by Sevik and Ça˘ in [108]), from a cursory look at their band
g
structures (Figure 4.3), one expects that Se4 should have large S values in the p-doped regime due
to heavy holes and large degeneracy. On the other hand, Se3 should be a better n-type thermo-

95

electric because of flatter conduction band with multiple minima. These expectations are indeed
observed in our calculations. For p=10, the maximum values of S are obtained at the same temperature (∼500 K) for both Se4 and Se3. However the value is larger in Se4 (400 µ V/K) than in Se3
(300 µ V/K). Thus Se4 is a better p-type thermoelectric. In contrast, for n-doping the maximum
|S| in the concentration range 1 < n < 100, is ∼200 µ V/K in Se4 whereas it is around 450 µ V/K
in Se3, making the later a much better n-type thermoelectric. Unfortunately, to the best of my
knowledge, there are no published thermopower measurements in n-doped Se3 and Se4 to confirm
the calculated results.

4.6 Summary
In this chapter, I have discussed the electronic structures of valence compensated ternary compounds given by the formula I3 -V-VI4 which belong to the family of Famatinite and Enargite.
Particular attention is given to Cu3 SbSe4 compound for its interesting properties reported experimentally. The calculations reveal that even though one can look at these compounds as derivatives
of Zinc blende or Wurtzite strucutre, which have been extensively studied and are quite well understood, the electronic properties of Famatinite and Enargite compounds are quite different.
First of all, unlike the monoatomic and binary tetrahedrally bonded compounds, the simple
octet rule cannot be applied to Famatinite and Enargite family, or other diamond-like ternary compounds. It is shown that V ions are primarily in the 3+ valence state, not 5+ as one would infer
from naive valence counting, verifying that there exist lone pair of electrons on the V site. The
effect of the local bonds as well as the role of V lone-pair electrons on the band structure are discussed. V lone pair electrons strongly interact with the surrounding VI-p orbitals, giving rise to
the special BOI (single band in Famatinite and two bands in Enargite) which play the role of the
lowest conduction band. V lone pairs are directly involved in the band gap formation. On the other
hand, the filled d-shells of Cu contribute indirectly in the band gap, only through changing the
bonding nature and bond length between V and VI atoms.
By comparing the electronic structure of a large number of ternary I3 -V-VI4 systems using dif96

ferent methods (local and non-local), we bring out the importance of the non-local exchange in gap
opening. My calculations, however, also show that GGA+U (with an appropriately chosen value of
U, about 15 eV) can reproduce extremely well the electronic structure near the gap. This method
is still very useful in surveying a wide range of materials, giving useful insight to understand the
physics of materials.
Results also suggests that one can easily tune the band gap to a desirable value by choosing
proper constituent or doping elements.
Famatinite and Enargite are expected to be better p-type thermoelectrics due to large degeneracy of the valence bands. Study of Cu3 SbSe4 shows that the band gap and thermopower values in
p-doped Cu3 SbSe4 agree very well with experiment,[106, 110, 111] showing that Se4 should be a
better p-type thermoelectric thermoelectric.
This work gives a foundation for further study of tetrahedrally bonded ternary materials, for
example to investigate the electronic structure of defects in these compounds.

97

CHAPTER 5
HALF-HEUSLER–HEUSLER COMPOSITE

5.1 Introduction
Half-Heusler and Heusler compounds have been of great interest for several decades for thermoelectric, magnetic, half-metallic and many other interesting properties. Development in the study
of this family of compounds can be found in several review articles,[85, 126–133] among which
Graf et al. [85] gave a quite comprehensive overview of the field. An interesting example of
such compounds is Zr-Ni-Sn, a promising thermoelectric material.[134–150] They can go from
semiconducting half-Heusler (HH)[139] limit, ZrNiSn with cubic MgAgAs structure, to metallic
Heusler limit (some times referred to as full-Heusler, FH), ZrNi2 Sn[139] with cubic MCu2 Al-type
structure, where M=Zr,Ti,Hf. Because of their excellent thermoelectric properties, ZrNiSn and
related HH systems, TiNiSn and HfNiSn, and their solid solutions have attracted considerable experimental and theoretical interest over last several years.[23, 140, 142, 143, 148–166] The idea
of going from HH to FH limit in other systems has also been tried. In the Ti-Co-Sn system for
example these two phases represent the edge members of the solid solution TiCo2−x Sn , where
the electrical conductivity changes from semiconducting to metallic as one decreases x from 1 to
0.[167]
Recently Makongo et al. [23, 157] have come up with an exciting observation, which shows
that large enhancements of both thermopower and electrical conductivity (hence the powerfactor
PF) of half-Heusler (HH) phases of (Zr,Hf)NiSn can be achieved at high temperatures through insertion of nanometer scale FH coherent inclusions within the HH matrix, caused by adding a small
concentration of extra Ni atoms. They denoted these systems as HH(1 − x)/FH(x) nanocomposites
which can be alternatively characterized as ZrNi1+x Sn, where 0 < x < 1. As shown in Fig. 5.1,
in Bi-doped HH(0.94)-HF(0.06) they found a PF of 3.3 mW/mK2 compared to 1.6 mW/mk2 for

98

(a)

Powerfactor (W/K2.m)

3x10-3
2x10-3

100% improvement

1x10-3

HH
Bi:HH
Bi:HH(97%)/FH(3%)
Bi:HH(94%)/FH(6%)
Bi:HH(90%)/FH(10%)

0

400

600
Temperature (K)

800

(b)

Figure 5.1: (Source: Makongo et al. [23]) Improvement of powerfactor (PF) (a), and TEM image
showing FH-nano-inclusions in HH(b)

simply Bi-Doped HH around 700 K. Chai and Kimura [160] also saw a large density of nanosized
FH precipitates within HF matrix in slightly Ni rich TiNiSn alloys and found an increase in the
power factor at high temperatures (S2 σ ∼ 3.5 mW/mK2 at 700 K compared to ∼2.5 mW/mK2 at
the same temperature for the parent compound TiNiSn[153]). They ascribed an energy filtering
mechanism proposed by Faleev and Léonard [22] to explain this enhancement. The basic idea
behind this filtering mechanism is to reduce the commonly observed electron-hole cancellation effect in thermopower by selectively allowing only one type of carrier to transport energy and charge.
There is however no direct evidence of such energy filtering in the presence of FH nano particles;
99

(a)
50
45
Ni (HH)
Sn (HH)
Zr (HH)
Zr (FH)
Ni (FH)
Sn (FH)

at%

40
35
30
25
HH

20
0

0.2

Two phase
0.4
0.6
x (Ni)

FH
0.8

1

DOS(EF) (states/eV)

1.1
1
0.9
0.8
0.7
0.6
DOS(EF)
Magenet.Suscept.

0.5
0.4
0

0.02

0.04 0.06
x(Ni)

0.08

12
10
8
6
4
2
0
-2
-4
-6
-8

χ (10-8cm3/g)

(b)

0.1

Figure 5.2: (Source: Romaka et al. [149]) (a) Phase diagram of HH-FH mixture, (b) Comparison
between measured susceptibility and density of state calculated using Korringa-Kohn-Rostocker
method with coherent potential approximation and local density functional
(KKR-CPALDA)[168].

more careful studies of electronic structure are needed to support this idea.
More recently, Romaka et al. [149] have carried out a thorough investigation of structural
phase transitions in half-Heusler - Heusler stannides (Zr,Hf)Ni1+x Sn and have shown how different physical properties change as one goes from HH to FH limit. They found that in ZrNi1+x Sn,

100

(a)

(b)

Figure 5.3: Crystal structure of (a) Half-Heusler (HH) and (b) Full-Heusler (FH) compounds,
where Zr (orange circles) and Sn atoms (gray circles) form a NaCl sub-lattice, inside which there
are 8 small cubes. HH is formed by filling every other cube with Ni atoms (blue circles), while
FH is formed by filling all the cubes. The empty-cube sites (quasi-particle of Ni-vacancy) in HH
is presented by a dashed-line circle.
macroscopic phase separation of HH and FH phases takes place in the region 0.3 ≤ x ≤ 0.7 (see
Fig. 5.2a). Their average XRD measurements matched very well to the cubic MgAgAs structure
for x ≤ 0.3 and to the cubic MCu2 Al-type structure for x ≥ 0.7. However the presence or absence of nanostructures of one phase in the matrix of the other phase in these concentration ranges
(x ≤ 0.3 or x ≥ 0.7) could not be established. One intersting result, relevant to this theis, is that
for 0 < x < 0.1, Romaka et al. [149] found and increase in the density of states (inferred from
susceptibility measurement – Fig. 5.2b) suggesting that new states appear at the Fermi-energy with
increasing x. They also reported a theoretical DOS, calculated using Korringa-Kohn-Rostocker
method with coherent potential approximation and local density functional (KKR-CPALDA)[168],
in good agreement with the susceptibility measurement. However, one has to note that KKR-CPA
method considers the excess atoms to be distributed randomly in the host matrix, thus, it cannot
address the effect of nanostructure on the electronic structure.
One can visualize the structure of these HH and FH compounds as follows (Fig. 5.3). One
starts from a NaCl-type lattice of (Zr,Hf)Sn consisting of two interpenetrating FCC sub-lattices.
Inside the large cubic unit cell of length a, there are 8 small cubes of length a/2. In (Zr,Hf)NiSn,
the centers of 4 out of these 8 cubes are occupied by Ni atoms so that each (Zr,Hf) or Sn is

101

tetrahedrally bonded to 4 Ni atoms. The other 4 small cube centers are empty. In (Zr,Hf)Ni2 Sn all
the 8 small cube centers are occupied by Ni atoms. Thus one can tune from HH to FH or vice versa
by varying the amount of Ni. Due to the additional Ni, FH compounds have an additional inversion
symmetry compared to HH. The lack of mutual solubility of the HH and FH phases with similar
crystal structures (as seen in Fig. 5.3) is quite remarkable in view of the fact that one can start from
one structure and either add or remove Ni atoms to reach the other structure without changing the
basic structure of the ZrSn matrix. Note that the lattice constants are however slightly different
between HH and FH (6.15 Å and 6.32 Å respectively).
There are several fundamental questions that arise in the study of these mixed HH-FH systems,
for example in ZrNi1+x Sn, such as: what is the nature of the phase diagram in the ZrNi1+x Sn in
the T vs x plane? what are the short and long range structural features and how they change with
x? what are the electronic and lattice properties of the mixed system? There are several recent
works on related issues: Romaka et al. [149] have addressed the electronic structure issue using
LDA-single site CPA approximation[168]; Kirievsky et al. [164] have studied the thermodynamics
of Ni-rich TiNiSn compound, i.e. the phase separation and energetics of antisite defects in the
simple cubic supercell of the relevant system. However, as pointed out earlier, CPA approximation
used by Romaka et al. [149] cannot address the clustering of nanostructure and the simple cubic
supercell used by Kirievsky et al. [164] in their electronic structure calculations is too small. I
will discuss some of their results later in this thesis. In this work I address two of these basic
issues. First, do the additional Ni atoms (or Ni vacancies) go randomly to the empty (or occupied)
cube sites or form some sort of local ZrNi2 Sn (or ZrNiSn) nanoclusters? Second, what is the
effect of these nanoclusters, if they form, on the electronic structure and vice-versa. In order to
address these questions I have carried out electronic structure calculations using ab-initio density
functional methods and quasi-periodic model[169] to describe the nanostructures.
The structure of this chapter is as follows. In Sec. 5.2, I discuss briefly some of the details
of electronic structure calculations and the model to describe the nanostructures, beyond which
already discussed in chapter 2. In Sec. 5.3, I discuss the energetics of nanoclustering and the

102

electronic structure, focusing on the nature of the gap states as one adds Ni atoms to ZrNiSn or the
modification in the density of states near the Fermi energy as one adds Ni-vacancies to ZrNi2 Sn.
Finally, in Sec. 5.4, I present a summary of my main findings and discuss further work in this area.

5.2 Methods of electronic structure calculation and modeling nanostructures
It is well known that LDA/GGA calculations usually underestimate the band gaps in semiconductors (as discussed in Chapter 2 as well as seen in Fe2 VAl (Chapter 3) and Cu3 SbSe4 (Chapter 4)),
however, the problem is opposite in the case of MNiSn (M=Ti, Zr, Fh)[139, 170] where the LDA
calculations predicted a band gap of ∼ 0.5 eV, two times larger than that of found from resistivity
measurements[135] (Eg = 0.12, 0.18 and 0.22 for M=Ti, Zr and Fh respectively). I have tested
the case of ZrNiSn with approximations on higher rung of DFT (Fig. 2.2), i.e. hybrid functional
HSE06, and have found the problem worsen (Eg =0.58 eV). Thus, in this chapter all total energy
and electronic structure calculations were done using density functional theory (DFT) within generalized gradient approximation (GGA). For the exchange-correlation potential we used the model
suggested by Perdew, Burke, and Ernzerhof (PBE).[47] We employed projector-augmented wave
(PAW) method[75, 76] as implemented in the VASP code.[77–79] with a plane-wave energy cutoff
of 400 eV and an energy convergence criterion (between two successive self-consistent cycles) of
10−4 eV (total energy/unit cell).
To study the evolution of HH-HF we have done calculations for ZrNi1+x Sn, where x = 0 − 1,
using the supercell method.[169] At x = 0 and x = 1 we have pure HH and FH compounds respectively. The supercell method has been used by Kirievsky et al. [164] for TiNi1 + xSn, however,
with simple cubic 1 × 1 × 1 cell (maximum 8 Ni atoms per unitcell) which gives a very large concentration of defects. They focused primarily on the thermodynamics of the HH-HF mixture and
the quasi-binary-system phase diagram. Here, I use a larger unitcell, 2 × 2 × 2 (maximum 64 Ni
atoms per unitcell) to reduce the interaction between defects in neighboring cells and focus on
the formation of nanostructure of FH (HH) in HH (FH). I study the dilute limit of the defect (i)

103

starting from the HH structure and then adding more Ni atoms to create the FH inclusion and (ii)
starting from the FH structure and then removing some Ni atoms to create HH inclusion. I use the
notation HH+nNi (FH-nNi) to refer to a 2 × 2 × 2 super cell of HH with n additional Ni (FH with
n deficient Ni). For the HH end, to study how the FH phase forms, I tried all possible configurations of excess Ni-pair and looked for the lowest-energy (the most favorable) configuration; once
the preferred pair-configuration was found, I added the third Ni and, again, looked for the lowest
energy configuration. An analogous procedure was applied for the FH end with Ni vacancies.
In all the calculations, structures were allowed to relax. The lattice parameters of HH and FH
are found to be 6.15 Å and 6.321 Å respectively. The defective structures, then, were derived
from the host structures and ionic positions were allowed to relax. The charge densities and total
energies were calculated self-consistently with an 8 × 8 × 8 Monkhorst-Pack[171] k-grid. The
electronic density of states (DOS) was then obtained using a denser k-mesh (12 × 12 × 12), from
the self-consistent charge density and wavefunction.
In studying properties of defects in a material, it is useful to calculate the formation energy of
the defect. There are extensive theoretical studies which have been dedicated to this problem. Here
I follow the work of Zhang et al. [172]. Applying to HH-FH, the formation energy of a defect can
be calculated by
∆H f (X) = ∆E f (X) + nµNi ,

(5.1)

∆E f (X) = E(X) − E(0) + nENi .

(5.2)

where

In these equations, E(X) and E(0) are the ground state energies of a system with defect X and
that of the host system respectively, n is the number of Ni added (n < 0) or removed (n > 0) from
the host-system (HH or FH), ENi is the ground state energy of Ni solid, and µNi is the chemical
potential of Ni. µNi is usually considered as a parameter which can vary from Ni-rich to Ni-poor
environment. The limit of µ can be defined from thermodynamic equilibrium conditions which
was carefully discussed by Zhang et al. [172].

104

Table 5.1: Formation energies (eV) of HH-FH systems (for the lowest energy configurations)
HH+1Ni
∆E f (X)

HH+2Ni

HH+3Ni

HH+4Ni

FH-1Ni

FH-2Ni

0.423

0.758

1.083

1.372

0.417

0.789

5.3 Results and Discussions
5.3.1
5.3.1.1

Formation energy and energetics of clustering
HH with excess Ni

Table 5.1 gives the energies needed to add (or remove) Ni to (from) HH or FH. For simplicity, only

µ = 0 is considered because the differences (not the absolute value) in energies are the studied. All
the energies are positive, which means it costs energy whenever a Ni is added or removed from HH
or HF systems. The energy required to add one Ni to HH host (0.423 eV) is slightly higher than
that needed to remove one Ni from FH matrix (0.417 eV). When the second Ni is added to HH,
the needed energy is smaller than that for the first one (0.758 − 0.423 = 0.335 eV), and so on for
the third Ni. The similarity is seen at the FH end. When more than one Ni is added to HH or taken
out of FH, the minimum energy configuration was found by varying the positions of the extra Ni
atoms or vacancies.
With the first excess Ni’s (Ni1 ) is fixed, when the second Ni (Ni2 is added to the HH matrix,
there are several possible empty cube sites that it can occupy. In Figure 5.4(left column) we plot
the ground state energy of HH+2Ni as a function of the distance between the added Ni-pair (Ni1
and Ni2 ). It is clear that two Ni are attracted, since the energy is lowest at the smallest distance
(4.22Å). The structure of the preferred positions of Ni is shown in figure 5.4, right column. The
ground state energy increases as the distance increases and then saturates at large distances, when
the distribution of excess Ni is more homogeneous over the space. The energy difference between
the lowest- and the highest-energy states is about 0.1 eV.
The saturation of the energy at large distance is due to the finite size of the periodic super-

105

Ground state energy (eV)

-648.84
-648.86
-648.88
-648.9
-648.92
-648.94
-648.96
4
5
6
7
8
9
10 11
Distance (Angstrom) between two excess Ni
Figure 5.4: Left column presents energetics of two extra Ni in HH: energy vs. distance. Right
column presents the preferable configuration of the excess Ni (orange circle) (only Ni matrix
(blue circles) is shown).
cell. Even with a 2 × 2 × 2 unit cell, the size is only 12.3Å×12.3Å×12.3Å, there is an artificial
periodicity of the defects. At large distance, the interaction between Ni2 (Ni1 ) with Ni1 (Ni2 ) in
the neighboring unit cells compensates the interaction between those in the same unit cell. For a
bigger unit cell, one should expect the saturation at a larger distance.
When the third Ni (Ni3 ) is added to HH, the problem is a little bit more complicated. For
convenience, we determine the position of Ni3 using the polar coordination (r,θ ) with respect to
the Ni1 -Ni2 , where the origin is at Ni1 and the x axis is along the Ni1 -Ni2 bond.
Fig. 5.5(left column) plots the relative positions of Ni3 in the Ni1 -Ni2 reference with the ground
state energy presented in color code varying from dark purple (low value) to yellow (high value).
Table 5.2: Ground state energy as a function of Ni3 ’s relative position with respect to Ni1 -Ni2
r(Å)

θ (◦ )

4.349
60

90

6.151
120

180

90

7.533
73.22

8.699
90

10.654

120

90

E (eV) -654.14 -654.10 -654.08 -654.06 -654.03 -654.03 -653.99 -653.99 -654.00

106

Figure 5.5: Left column presents energetics of three extra Ni in HH: position of the third Ni
(Ni3 ) is given in polar coordinate, with respect to Ni1 -Ni2 , whose positions were previously
determined by studying the energetics of a pair of excess Ni (FIG. 5.4). Ground state energy of
each configuration is represented by color code from low (purple) to high (yellow) energy. Right
column shows the preferred configuration of the Ni3 (only Ni matrix (blue circle) is shown). The
excess Ni are orange circles.

The value of the ground state energy of each configuration is given in table 5.2. The results confirm
the attractive interaction between 3 excess Ni atoms, similar to which found for a pair of excess Ni.
For the same r, the energy, in general, increases with increasing θ . For instance, for the smallest
distance, r = 4.349, the ground state energy is −654.14 eV at θ = 60◦ , whereas it is −654.06 eV
at θ = 180◦ . Similarly, for the same θ the ground state energy increases with increasing distance.
The maximum energy difference between two configurations is about 0.14 eV.
Similar to the case of Ni-pair, with 3 excess Ni, the ground state energy also saturates at large
distance, and for the same reason (periodic supercell).
The lowest ground state is found for Ni3 at (4.34Å, 60◦ ), where Ni3 is closest to both Ni1 and
Ni2 ; the crystal structure of this configuration is given in Fig. 5.5(right column). Excess Ni atoms
tend to come close to each other and form a nano-cluster of FH in the HH matrix.

107

Figure 5.6: Conventional notations for Ni-sites: “f” for Ni at (1/4,1/4,1/4) and equivalent atoms
(dark blue), “h” for Ni at (3/4,3/4,3/4) and equivalent atoms (bright yellow).

The results strongly support the HH-FH phase separation observed in experiments[23, 149,
158–160, 163] In experiment, however, the phase separation is only observed at a high concentrations of Ni, for example, Romaka et al. [149] showed that the two macroscopic phases of
ZrNi1+x Sn can be seen only within the range of 0.30 < x < 0.65 (Fig. 5.2a). The present results, on the other hand, predict that the nanostructures of FH form in a HH matrix, at least at
low temperature, even at small concentrations of excess Ni(x ∼ 0.1, or 3 excess Ni in 32 host Ni),
thereby supporting the ideas put forward by Makongo et al. [23]

5.3.1.2

FH with deficient Ni

At the FH end, one should note that there are two different Ni-sites: “h”-site (h) at (1/4,1/4,1/4) and
its equivalent, which is occupied by Ni in HH, and “f”-site (f) at (3/4,3/4,3/4) and its equivalent,
which is empty site in HH (Fig. 5.6). When removing Ni, the Ni-vacancy can be located at either
of these two sites.
If one Ni is removed, due to symmetry, two sites are equivalent. The energy cost for removing
1 Ni out of 2 × 2 × 2 unit cell is ∼ 0.417 eV. Similar to the case of HH+nNi, the energy cost to
remove the second Ni (to add the second vacancy) is smaller than that needed to remove only one
Ni (Table 5.1). When 2 Ni are removed, for the convenience, the first vacancy (V1) is assumed
to be at a f site, the second one can be at either h or f site. We denote the pair of vacancies by

108

(a)

-804.8
f-f
f-h

Ground state energy (eV)

f-h1
-804.85
f-h2
-804.9

f-f3

-804.95
-805

f-f4

f-h3
f-f1
f-f2

-805.05

f-h4

-805.1
3
4
5
6
7
8
9
10
11
Distance (Angstrom) between two Ni-vacancies in FH
(b) f-f1

(c) f-f2

(d) f-h4

Figure 5.7: The top figure shows energetics of FH with two Ni-vacant sites (white circles); the
bottom shows structures of the three most preferable configurations, which have ground state
energies close to each other: (a) f-f1, (b) f-f2, and (c) f-h4. Only Ni matrix (blue circles) is shown.

109

the name of the sites which they occupy,i.e. f-f and f-h, note that f-f and h-h are equivalent. In
Fig. 5.7(a), we plot the ground state energies as a function of the distance between 2 Ni-vacancies,
considering f-f and f-h separately. Fig. 5.7(a) shows that f-f and f-h are energetically different.
For the f-h case, the interaction between two vacancies is repulsive,i.e. the total energy is the
lowest for the largest distance. The relation between energy and distance is almost linear. The
difference between the lowest and the highest energies is about ∼ 0.22 eV. For the f-f case, in
contrast, the interaction between the two vacancies is, in general, attractive, the smaller is the
distance, the lower is the energy. There is, however, a competing effect between repulsion and
attraction at small distance, which makes the ground state energy higher at the smallest distance,
4.47Å (f-f1), than at the larger one, 6.32Å (f-f2) (-805.00 eV for f-f1 compared to -805.055 eV for
f-f2, as shown in Fig. 5.7(a)). The energy saturates at large f-f distance.
The structure of f-f2 (Fig. 5.7(c)) has 2 Ni-vacancies distributing evenly at the every other site
along lattice axis (say x), forming a quasi one dimensional chain. This structure is similar to the
configuration (i) in reference [164], except that they put one vacancy in 1×1×1 simple cubic super
cell, thus it forms a 3-dimensional ordered structure. In their work, Kirievsky et al. [164] calculated
electronic structures of 11 configurations with different fraction of Ni over the eight cubic site in
TiSn matrix and studied their energetics using a quasi binary model (TiNi2 Sn)c (TiSn)1−c , where
the formation energy is given by
mixture − (c · E TiNi2 Sn + (1 − c) · E TiSn ),
∆U = Etot
tot
tot

(5.3)

where Etot is the total energy of corresponding system. They found that the configuration (i),
c = 0.875, had negative free-energy of formation whereas most of the other configurations (except
c = 0.5 which corresponds to Half-Heusler phase TiNiSn) have positive free-energy of formation.
One should note that the definition of formation energy use by in reference [164] is different from
the definition used in this thesis. Furthermore, they have used a 1 × 1 × 1 supercell where the
periodic cell constraints are much significant. Thus, it is not practical to directly compare values
of formation energies between two works.

110

Energy (eV)

(a)
2

(b)
2

1

1

0

0

-1

-1

DOS

-2 R
(c)
1000

Γ

X

M

Γ -2 R

EF

Γ

M

Γ

0

X

1

200
HH
FH

500
0 -1
0 -2

-1

0
Energy (eV)

1

2

Figure 5.8: Band structures of (a) HH and (b) FH in the Brillouin zone (BZ) of 2 × 2 × 2 cubic
supercell; and (c) the density of states (DOS) within energy range [-2,2]; in the onset, we zoom in
the range [-1,1] to show the fundamental difference between HH and FH: the former is metal
whereas the latter is semiconductor, with a band gap ∼ 0.5 eV.
It is interesting that the ground state energy of f-f2 (the lowest energy among f-f) is about
∼ 0.012 eV lower than that of f-h4 (the lowest among f-h) (see Fig. 5.7). Thus, f-f2 and f-h4 can
coexist at x = 0.9375. This suggest a frustration of excess Ni-vacancies in the FH matrix.
To summarize, supercell calculations show that the energetics of FH-nNi is different from that
of HH+nNi. In the former case, the Ni-vacancy tend to distribute almost evenly in spacealthough
there is some suggestion for vacancy ordering. Thus, the nanostructure of the HH phase is not
formed in FH at low vacancy concentrations, whereas in the latter case, excess Ni tend to form
nano-clusters of FH in HH.

111

5.3.2

Electronic structure

Before we discuss the electronic structure of HH-FH mixture, let us first review the electronic structures of the end compounds ZrNiSn and ZrNi2 Sn. Electronic structures of ZrNiSn and ZrNi2 Sn are
well understood and were extensively studied both theoretically[139, 170] and experimentally[135,
136]. ZrNiSn is a semiconductor with an indirect band gap of ∼ 0.5 eV[139, 170, 173] between Γ
and X points in the FCC Brillouin Zone (BZ). The band gap is formed from the hybridization of
Zr-d and Sn-p orbitals. The band gap formation was carefully discussed by Larson et al. [173]. Top
of the valence bands are three-fold degenerate, one of which has a light effective mass while the
others have heavier ones. On the other hand, bottom of the conduction bands is a non-degenerate
band with anisotopic effective mass. All Ni-d and Sn-s bands are occupied, where the latter is
located at about -15 eV below the valence band maximum. ZrNi2 Sn, on the other hand, is a good
metal.
Fig. 5.8 gives the band structures of ZrNiSn (Fig. 5.8a) and ZrNi2 Sn (Fig. 5.8b) in the Brillouin
zone (BZ) of 2×2×2 simple-cubic supercell, and the corresponding density of state (DOS)(Fig. 5.8c).
Going from a fcc primitive unit cell to 2 × 2 × 2 simple-cubic supercell, due to band folding – since
the BZ is smaller in the larger cell – the band structure of ZrNiSn shows a direct (instead of indirect
as in fcc) band gap or ∼ 0.5 at the Γ point, in good agreement with other works[139, 170, 173].
The band structure of both ZrNiSn and ZrNi2 Sn appear more complicated in the 2 × 2 × 2 supercell since we have more bands in the first BZ. Fig. 5.8c shows comparison of the DOS of the 2
compounds, near the Fermi level. The inset shows the DOS close to the Fermi energy.
When extra Ni (Ni1 ) is introduced into HH matrix (or removed from FH), there is a change
in the electronic structure and a redistribution of the charge density with respect to the host system (HH or FH). To analyze this charge redistribution, we calculate the charge-density difference
(charge-difference) ∆ρ = ρ X (r) − ρ 0 (r), where ρ X (r) and ρ 0 (r) are the charge densities of the
defective system and corresponding host system. Fig. 5.9 depicts the charge difference contour
on the (-110) plane for HH+1Ni (Fig. (5.9a) and FH-1Ni (Fig. (5.9b), where blue indicates charge
depletion and red indicates charge accumulation.

112

(a)

(b)

Figure 5.9: Charge difference between (a) HH+1Ni and (b) FH-1Ni and their parent compounds
(HH and HF), where blue color indicates charge depletion whereas red color indicates charge
accumulation. The figure show the contour plot on (-110) plane. Charge difference is defined as
ρ X (r) − ρ 0 (r), where ρ X (r), and ρ 0 (r) are charge density of defective and pure system
respectively. Host Ni are blue circles, excess Ni is red circle, Zr are brown circles, Sn are gray
circles.
For HH+1Ni, the additional Ni introduces more electrons to the HH system, these electrons
localize around Ni1 (Fig. 5.9a). The additional electrons at the excess Ni site cause the rearrangements of charge in the neighboring space. However, this rearrangement occurs mostly at the
Ni-sites. Electrons tend to move further away from Ni1 . In Fig. 5.9a one can clearly see the blue
lobes at the Ni-sites are directed towards Ni1 , whereas the red lobes at the same sites are directed
away from Ni1 . In Fig. 5.9 one sees that the Ni atoms of the host matrix (HH) connect to N1 in
directly through Zr or Sn, let us call these connections Zr- and Sn-channels. With those channels
defined, one can realize that the charge rearrangements are not the same on all Ni, they are larger on
Ni which interacts with the Ni1 through only one channel, either Sn-channel or Zr-channel, where
Ni-Zr(Sn)-Ni1 form a straight line. The charge rearrangements on other Ni-sites, which connect
to Ni1 through two channels of the same type (either two Sn-channels or two Zr-channels), are
smaller even though they may be at shorter distances with respect to Ni1 .
A similar picture is also seen in HF-1Ni (Fig. 5.9b), except that the removal of Ni (V1) takes
away electrons. A subtle difference between FH-1Ni and HH+1Ni is that there is an interesting
charge rearrangement on Ni indicated by arrows in Fig. 5.9b (hereafter it is called Nix ). Nix indirectly connect to V1 through two channels, one Zr-channel and one Sn-channel, the Nix -Sn(Zr)-V1

113

(a)

(b)

(c)

Energy (eV)

1.5
0
-1.5
R
(d)

+4Ni
+3Ni
+2Ni
+1Ni
HH -2

Γ X M

-1

ΓR

Γ X M

0
Energy (eV)

ΓR

Γ X M

1

Γ

2

Figure 5.10: Evolution of electronic structure of HH with excess Ni: band structures (a) HH+1Ni,
(b) HH+2Ni, and (c) HH+3Ni, and (d) density of states of HH with 1, 2, 3 and 4 additional Ni.
angle is about ∼ 110◦ . The charge rearrangement on Nix site does not occur on the Ni-V1 direction
but on the Ni-Sn direction, as shown in Fig. 5.9b. The blue lobe around Nix is closer to V1 than
the red lobe, opposite to those around the other Ni atoms with linear link Ni-Sn(Zr)-V1, whose red
lobes are closer to V1. This difference in charge rearrangement is most likely the reason for the
difference in energetics of FH-nNi from HH+nNi and is a reflection of different local bonding in
FH and HH systems.
In Fig. 5.10, we give the band structures and DOS of HH+nNi systems, note that hereafter we
consider only the lowest-energy state for each n. As one can see in Fig. 5.10(a,b,c), adding Ni
to HH matrix has two effects: splitting of the host HH bands due to impurity induced symmetry
lowering, and introducing defect states inside the gap. There is still some local residual symmetry
which preserve the degeneracy of the valence band maximum at the Γ point. The splitting of the
HH conduction bands reduce the band gap from ∼ 0.5 eV to ∼0.37 eV. The dispersion of the host

114

Figure 5.11: The isosurface of the charge density associated with the in-gap defect states in
HH+1Ni. Host Ni are blue circles, excess Ni is the red circle, Zr are brown circles, Sn are gray
circles.
bands near the gap changes slightly. There are two defect states per additional Ni located just
below the HH conduction band edge. Projected density of states for HH+nNi reveal that the defect
states in the gap are associated primarily with Ni-eg orbitals. When one excess Ni is added to the
HH matrix, it brings 10 more electrons which occupy 5 Ni1 -d bands, 3 t2 g and 2 eg. All the Ni1 -d
bands are below the Fermi level, the t2 g bands are located at ∼ −1.5 eV, whereas the eg bands are
located right below the HH conduction bands, giving rise to the in-gap states. The defect states in
the gap are localized, its bandwidth is ∼ 0.14 eV. The charge density associated with the in-gap
defect states is given in Fig. 5.11, clearly showing the Ni-eg character. The dispersion of the defect
bands is small in x, y, or z directions, but larger in other directions, such as Γ-M (110), Γ-R (111),
which are directions of Ni1 -Zr(Sn) bonds.
When more Ni are added to the HH matrix, more impurity bands are introduced into the gap
region, and eventually fill up the gap, converting the semiconducting state of HH into metallic
state of FH. Comparing the DOS of HH+nNi (Fig. 5.10d), one can clearly see the evolution of
the electronic structure. With one extra Ni, DOS starts to show a peak right below the conduction
bands. As the number of excess Ni increases the peak gets higher and shifts downwards. These

115

(a)

(b)

(c)

Energy (eV)

1.5
0
-1.5
R
(d)

Γ X M

ΓR

Γ X M

ΓR

Γ X M
FH
FH-1Ni
FH-2Ni

60
DOS

Γ

400
0
-1.2 -1.1 -1 -0.9

200
0
-2

-1

0
Energy (eV)

1

2

Figure 5.12: Evolution of electronic structure of FH with deficient Ni: band structures (a) FH, (b)
FH-1Ni, and (c) FH-2Ni , and density of states of those systems. The inset show the DOS,
indicated by the arrow, around -1 eV in small range [-1.2,-0.9].

localized defect states are occupied and act as donor-states which can inject carriers into the conduction bands and the system behaves like a n-type semiconductor. These states can affect the
thermoelectric properties of HH-FH composite systems.
The effect of Ni-vacancy in FH is less subtle than that of Ni excess in HH. Fig. 5.12 gives
the electronic structure of FH-nNi. The dominant effect of Ni-vacancy on the band structure of
FH is the band splitting, which makes it appears to be more complicated. One can see some
extra bands appearing around ∼-1 eV. These extra bands give rise to a small peak in the DOS at
∼-1 eV (Fig. 5.12d inset). These effects, however, just introduce small perturbation in the DOS
(figure 5.12d), which do not change the metallic properties of FH, at least for small concentration
of Ni-vacancy.

116

5.4 Summary and Conclusion
To summarize this chapter, I have discussed the energetics and electronic structures of the HHFH mixture, ZrNi1+x Sn. The evolution from HH to FH and from FH to HH was studied using a
2 × 2 × 2 supercell model. The focus is on the small concentration (< 10%) limit of defects (excess
Ni in HH and Ni-vacancy in FH).
The calculations predict an attractive interaction between excess Ni atoms in the HH matrix,
i.e. the excess Ni atoms tend to come close to each other and form nano-clusters within the HH
matrix. Vacancies in the FH matrix, on the other hand, are found to possibly occupy either of
two crystallographic sites, h(1/4,1/4,1/4) or f(3/4,3/4,3/4), and their energetics are distinctively
different from that of excess Ni in the HH matrix. If two Ni-vacancies in FH occupy the equivalent
sites(either type h or f), the interaction between them is attractive at the distance larger or equal
6.32 Å and repulsive at smaller distance. If two Ni-vacancies in FH occupy the sites of different
types (one at h-site and another at f-site), then the interaction is repulsive, the two vacancies tend
to stay away from each other. Thus, at low concentration of vacancies, the nano-structure of HH
in FH matrix is energetically not preferred.
Analysis of the change charge density shows that excess Ni (or Ni-vacancy) in HH (FH) causes
large charge rearrangement in the host matrix. The charge rearrangement occurs mostly on Nisites and is larger at Ni sites which are connected to the excess Ni (Ni-vacancy) indirectly through
one channel mediated by either a Zr or a Sn atom. The charge rearrangement is modest at the Nisites which are connected to the defect site through channels mediated by two atoms of the same
type. FH with one Ni-vacancy has an interesting charge rearrangement at the Ni-sites connected
to the vacancy through two channels mediated by one Zr and one Sn. This charge rearrangement
is in different directions compared to that at other Ni-sites and is most likely responsible for the
difference in energetics between FH-nNi and HH+nNi.
While the effect of Ni-vacancies in FH on the electronic structure is modest, the excess Ni
in HH changes the HH’s electronic structure dramatically. The excess Ni in HH introduces new
impurity states in the gap of HH band structure, giving rise to finite density of states below the
117

conduction bands. Due to the in-gap impurity states, HH with excess Ni can behave like a n-type
semiconductor.
The present study of ZrNi1+x Sn agrees well with experimental and other theoretical studies
on crystal and electronic structures of this mixture. The calculated results however cannot confirm
or deny the energy filtering effect of the FH nano-inclusion in the HH matrix, which has been
suggested by several authors to be the reason for the improvement in the material’s thermoelectric
properties. In order to understand that problem, larger-supercell calculations and transport studies
are needed.

118

CHAPTER 6
CONCLUSION AND OUTLOOK

Thermoelectrics are an important part in the solution of the energy problems. They can contribute
to both improving energy conversion efficiency (saving waste heat), and power generation. In
the last several years the field of thermoelectrics has developed quickly; however, there are still
emerging needs of new discoveries in material for thermoelectric applications, including finding
new compounds and finding new ways of improving the existing ones.
The job of increasing the figure of merit ZT of thermoelectrics is challenging since there is
a competing relation between quantities that are involved in ZT . There are several ideas of how
to improve ZT including, but not limited to, introducing resonant states, energy filtering, creating
aharmonicity, etc. , all of which, now, can involve the inclusion of low-dimensional structures, and
require a deeper understanding of electronic structure and transport properties in complex multicomponent systems.
The advances in theories and computational methods, density functional theory in particular,
of electronic structure have been extremely helpful in studying physical properties of materials.
Density functional theory calculations have improved greatly in the last decade, in both accuracy
and productivity. The improvement in density functional theory calculations is due to the improvement of approximations used in treating the effect of electron-electron interaction, more interaction
and correction terms are added, going from LDA/GGA to hybrid functional and to random phase
approximation.
In this thesis, I have employed density functional theory first principles cacluations, using several approximations including LDA/GGA, LDA+U, meta-GGA and hybrid functional approximation, to study electronic structures and transport properties of several novel thermoelectric materials, including the Heusler compound Fe2 VAl, Cu3 SbSe4 and its related tetrahedrally coordinated
systems, and the nanocomposites of Half-Heusler and Heusler compounds.

119

Chapter 3 discusses the electronic structure and thermopower of the Heusler compound Fe2 VAl.
Using LDA+U, I have shown that Fe2 VAl is a narrow band gap semiconductor (Eg =0.5 eV),
rather than a pseudo-gap system as predicted in LDA/GGA (negative gap). The values of intra-site
Coulomb repulsion, UFe = 4 eV and UV = 1.5 eV, were calculated using constrained density functional theory. Using electronic structure obtained from LDA+U, I have studied the thermopower of
Fe2 VAl and have shown that the results are in better agreement with experiment, compared to the
results calculated with GGA band structure. Fe2 VAl is shown to be a better p-type thermoelectric
due to large degeneracy at the top of the valence bands. I have also noted that the better agreement
with experimental data can be achieved if a smaller value of U (1 eV) is used. Furthermore, the
pseudo-gap behavior observed in certain experiments is still not understood. This may be due to
new physics, such as defects, dynamic interaction or excitonic effect and further works need to be
done in order to have a better understanding of this material.
In chapter 4, I discussed the band gap problem in Cu3 SbSe4 and related tetrahedrally coordinated compounds I3 -V-VI4 (belonging to the Famatinite and the Enargite families), where I=Cu,
V=P, As, Sb, or Bi and VI=S, or Se, and have pointed out the importance of lone-pare electron of
group V elements in the band gap formation in this family of materials. Both GGA, LDA+U and
hybrid functional calculations were done and discussed. A simple bonding-antibonding scheme
explaining the band-gap formation has been proposed. In this simple picture, the lone-pair electrons of V stereochemically interact with the surrounding VI-p orbital, giving rise to a special
band (BOI) which is the lowest of the conduction bands, forming a band gap. Exceptions are those
systems which contain Bi, all of them are found to be (semi)metal in all the approximations. I
have also studied the relation between the value of band gaps and the bond length between V and
VI elements and have shown that the relation is almost linear, band gap decreases with increasing
bong length. The results suggested that one can tune the band gap in these compounds by choosing
proper composition of dopants.
Chapter 5 has been dedicated to the study of the problem of nanostructure inclusion of Heusler
(ZrNi2 Sn) compounds in the Half-Heusler (ZrNiSn) matrix and vice versa, focusing on the low

120

density limit using 2 × 2 × 2 periodic supercell model. The electronic structures of the nanocomposites have been discussed. It is important to note that all approximations to density functional
theory tend to overestimate the band gap of ZrNiSn. The results predict that extra Ni atoms in
ZrNiSn tend to form nano clusters even at low concentrations, but excess Ni-vacancies in ZrNi2 Sn
do not. The latter tend to distribute almost evenly in space. These findings are in agreement with
experiments. Inclusion of ZrNi2 Sn phase in ZrNiSn introduces in-gap impurity states corresponding to a finite density of state right below the conduction band edge, the more excess Ni are put,
the more states are introduced, in agreement with susceptibility measurement. These in-gap states
change the electronic structure, and consequently the thermoelectric properties of ZrNiSn. The
nanocomposite of ZrNiSn-ZnNi2 Sn acts like a n-type semiconductors. The detail effects of the
nanostructure inclusion on thermoelectric properties of ZrNiSn, however, have not been discussed
in this thesis. Further works are under progress.
In conclusion, I have studied several systems of thermoelectric significance and have contributed to a better understanding of some of their fundamental issues. The results have shown that
even though climbing up the density functional theory ladder of approximation, in general, improves the accuracy of electronic structure calculations, there is no universal approximation which
is good enough to be used for all the materials – LDA+U gave reasonable results for Fe2 VAl but
failed in predicting band gaps in tetrahedrally coordinated systems, and both GGA, LDA+U and
hybrid functional calculation overestimate the band gap in ZrNiSn. It is clear that there are still
some outstanding questions which have not been addressed in this thesis. For examples: What
are dynamic correlation effects in pseudo-gap systems involving d-electrons? How can one make
critical evaluations of non-local exchange (hybrid functional) theories vis-a-vis mBJ, and other approximations? What defects can occur in Heusler and tetrahedrally coordinated systems, and their
roles in electronic structures and thermoelectric properties of these compounds? What are the effects of nano-inclusion on thermoelectric properties of Half-Heusler–Heusler composite? Whether
transport and energy filtering effects are taken place in such the system? Further work addressing
some of those questions is under progress.

121

BIBLIOGRAPHY

122

BIBLIOGRAPHY

[1]

World energy outlook 2012 (International Energy Agency, Organisation for Economic Cooperation and Development, 2012), retrieved from http://www.oecd-ilibrary.org.

[2]

T. Seebeck, Proc. Prussian Acad. Sci. (1825-273).

[3]

G. Magnus, Poggendorf´ Annalen der Physik 83, 469 (1851).
s

[4]

I. B. Cadoff and E. Miller, eds., Thermoelectric materials and devices, Materials Technology
Series (Reinhold Publishing Corporation, New York, 1960).

[5]

G. Mahan, B. Sales, and J. Sharp, Physics Today 50, 42 (1997).

[6]

R. R. Heikes and J. Roland W. Ure, Thermoelectricity: science and engineering (Interscience publishers, New York - London, 1961).

[7]

M. V. Vedernikov and E. K. Iordanishvili, in XVII International Conference on Thermoelectrics, 1998. Proceedings ICT 98 (1998), pp. 37–42.

[8]

NASA, Radioisotope power systems for space exploration (2011), nASA facts.

[9]

New horizons mission powered by space radioisotope power systems, URL http://
energy.gov/ne.

[10] R. R. Furlong and E. J. Wahlquist, Nuclear News 42, 2635 (1999).
[11] G. J. Snyder, The Electrochemical Society Interface p. 54 (2008), and therein references.
[12] J. Yang and F. R. Stabler, Journal of Elec Materi 38, 1245 (2009).
[13] C. B. Vining, Nat Mater 8, 83 (2009).
[14] M. Zebarjadi, K. Esfarjani, M. S. Dresselhaus, Z. F. Ren, and G. Chen, Energy Environ. Sci.
(2011), and therein references.
[15] T. M. Tritt, Annual Review of Materials Research 41, 433 (2011).
[16] J. W. Fairbanks, Automotive thermoelectric generators and hvac (2012), presented at Directions in Engine-Efficiency and Emissions Research (DEER) Conference 2012, URL http:
//www1.eere.energy.gov/vehiclesandfuels/resources/conferences/deer/.
[17] G. D. Mahan and G. O. Sofo, in Proc. Natl. Acad. Sci. (1996), vol. 93, p. 7435.
[18] J. O. Sofo and G. D. Mahan, Phys. Rev. B 49, 4565 (1994).
[19] E. J. Skoug, C. Zhou, Y. Pei, and D. T. Morelli, Jr. Electronic. Matter 38, 1221 (2009).
[20] S. Ahmad, K. Hoang, and S. D. Mahanti, Phys. Rev. Lett. 96, 056403 (2006).

123

[21] J. P. Heremans, V. Jovovic, E. S. Toberer, A. Saramat, K. Kurosaki, A. Charoenphakdee,
S. Yamanaka, and G. J. Snyder, Science 321, 554 (2008).
[22] S. V. Faleev and F. Léonard, Phys. Rev. B 77, 214304 (2008).
[23] J. P. A. Makongo, D. K. Misra, X. Zhou, A. Pant, M. R. Shabetai, X. Su, C. Uher, K. L.
Stokes, and P. F. P. Poudeu, J. Am. Chem. Soc. 133, 18843 (2011).
[24] D. M. Rowe, Journal of Physics D: Applied Physics 7, 1843 (1974).
[25] G. A. Slack and M. A. Hussain, Journal of Applied Physics 70, 2694 (1991).
[26] D. T. Morelli and G. A. Slack, High Thermal Conductivity Materials (Springer, New York,
2006), chap. 2, p. 37.
[27] E. J. Skoug, J. D. Cain, and D. T. Morelli, Appl. Phys. Lett. 96, 181905 (2010).
[28] Y. Zhang, E. Skoug, J. Cain, V. Ozolins, D. Morelli, and C. Wolverton, Phys. Rev. B 85,
054306 (2012).
[29] M. Faraday, Experimental researches in electricity (printers And Publishers To The University Of London, London, 1949), vol. I, 2nd ed., sect. 432-6.
[30] A. H. Wilson, Proc. R. Soc. Lond. A 133, 458 (1931).
[31] H. H. Lin and J. B. Kuo, in Wiley Encyclopedia of Electrical and Electronics Engineering
(John Wiley & Sons, Inc., 2001), ISBN 9780471346081.
[32] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
[33] J. Chu and A. Sher, Introduction (Springer New York, 2008), pp. 1–17.
[34] G. D. Mahan, 65, 1578 (1989).
[35] D. Do, M.-S. Lee, and S. D. Mahanti, Phys. Rev. B 84, 125104 (2011).
[36] D. Do, V. Ozolins, S. D. Mahanti, M.-S. Lee, Y. Zhang, and C. Wolverton, J. Phys.: Condens. Matter 24, 415502 (2012).
[37] D. Do and S. D. Mahanti, arXiv e-print 1306.0503 (2013).
[38] R. M. Martin, Electronic Structure, Basic Theory and Practical Method (Cambridge University Press, 2004).
[39] L. H. Thomas, Proc. Cambridge Philos. Soc. 23, 542 (1927).
[40] E. Fermi, Atti. Accad. Naz. Lincei, Cl. Sci. Fis. Mat. Nat. Rend. 6, 602 (1927).
[41] E. Schrodinger, Phys. Rev. 28, 1049 (1926).
[42] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).

124

[43] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[44] W. Kohn, Rev. Mod. Phys. 71, 1253 (1999).
[45] M. Levy, J. P. Perdew, and V. Sahni, Phys. Rev. A 30, 2745 (1984).
[46] J. P. Perdew, A. Ruzsinszky, L. A. Constantin, J. Sun, and G. I. Csonka, J. Chem. Theory
Comput. 5, 902 (2009).
[47] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
[48] J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
[49] E. Engel and S. H. Vosko, Phys. Rev. B 47, 13164 (1993).
[50] Z. Wu and R. E. Cohen, Phys. Rev. B 73, 235116 (2006).
[51] J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys. Rev. Lett. 91, 146401
(2003).
[52] F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009).
[53] A. D. Becke, Phys. Rev. A 38, 3098 (1988).
[54] C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
[55] A. D. Becke, J. Chem. Phys. 96, 5648 (1993).
[56] A. D. Becke, J. Chem. Phys. 104, 1040 (1996).
[57] D. I. Bilc, R. Orlando, R. Shaltaf, G.-M. Rignanese, J. Íñiguez, and P. Ghosez, Phys. Rev. B
77, 165107 (2008).
[58] C. Adamo and V. Barone, The Journal of Chemical Physics 110, 6158 (1999).
[59] J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 118, 8207
(2003).
[60] J. Heyd and G. E. Scuseria, The Journal of Chemical Physics 121, 1187 (2004).
[61] J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 124, 219906
(2006).
[62] P. A. M. Dirac, Proc. Cambridge Philos. Soc. 26, 376 (1930).
[63] J. C. Slater, Phys. Rev. 81, 385 (1951).
[64] A. Svane and O. Gunnarsson, Pys. Rev. Lett. 65, 1148 (1990).
[65] V. I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein, Journal of Physics: Condensed
Matter 9, 767 (1997).
[66] V. I. Anisimov and O. Gunnarsson, Phys. Rev. B 43, 7570 (1991).
125

[67] O. Gunnarsson, O. K. Andersen, O. Jepsen, and J. Zaanen, Phys. Rev. B 39, 1708 (1989).
[68] G. K. H. Madsen and P. Novák, Europhys. Lett. 69, 777 (2005).
[69] V. I. Anisimov, M. A. Korotin, and E. Z. Kormaev, J Phys.: Condens. Matter 2, 3973 (1990).
[70] M. J. Han, T. Ozaki, and J. yu, Phys. Rev. B 73 (2006).
[71] L. Hedin, Phys. Rev. 139, A796 (1965).
[72] A. D. Becke, The Journal of Chemical Physics 98, 1372 (1993).
[73] J. M. Ziman, Principles of the Theory of Solids (Cambridge University Press, 1964).
[74] G. K. Madsen and D. J. Singh, Computer Physics Communications 175, 67 (2006).
[75] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
[76] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
[77] G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
[78] G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996).
[79] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
[80] K. Schwarz and P. Blaha, Computational Materials Science 28, 259 (2003), proceedings of
the Symposium on Software Development for Process and Materials Design.
[81] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev.
B 57, 1505 (1998).
[82] Y. Nishino, M. Kato, S. Asano, K. Soda, M. Hayasaki, and U. Mizutani, Phys. Rev. Lett. 79,
1909 (1997).
[83] C.-S. Lue and J. H. Ross, Phys. Rev. B 58, 9763 (1998).
[84] H. Okamura, J. Kawahara, T. Nanba, S. Kimura, K. Soda, U. Mizutani, Y. Nishino, M. Kato,
I. Shimoyama, H. Miura, et al., Phys. Rev. Lett. 84, 3674 (2000).
[85] T. Graf, C. Felser, and S. S. Parkin, Progress in Solid State Chemistry 39, 1 (2011).
[86] J. Korringa, Physica 16, 601 (1950).
[87] D. J. Singh and I. I. Mazin, Phys. Rev. B 57, 14352 (1998).
[88] M. Weinert and R. E. Watson, Phys. Rev. B 58, 9732 (1998).
[89] R. Weht and W. E. Pickett, Phys. Rev. B 58, 6855 (1998).
[90] G. Y. Guo, G. A. Botton, and Y. Nishino, J. Phys.: Condens. Matter 10, L119 (1998).

126

[91] M. Kumar, T. Nautiyal, and S. Auluck, Journal of Physics: Condensed Matter 21, 446001
(2009).
[92] C. S. Lue and Y.-K. Kuo, Phys. Rev. B 66, 085121 (2002).
[93] Y. Nishino, S. Deguchi, and U. Mizutani, Phys. Rev. B 74, 115115 (2006).
[94] M. Vasundhara, V. Srinivas, and V. V. Rao, Phys. Rev. B 77, 224415 (2008).
[95] Y. Nishino, Intermetallics 8, 1233 (2000).
[96] D. I. Bilc and P. Ghosez, Phys. Rev. B 83, 205204 (2011).
[97] H. C. Kandpal, G. H. Fecher, C. Felser, and G. Schönhense, Phys. Rev. B 73, 094422 (2006).
[98] B. Balke, G. H. Fecher, H. C. Kandpal, C. Felser, K. Kobayashi, E. Ikenaga, J.-J. Kim, and
S. Ueda, Phys. Rev. B 74, 104405 (2006).
[99] P. Larson, S. D. Mahanti, S. Sportouch, and M. G. Kanatzidis, Phys. Rev. B 59, 15660
(1999).
[100] M.-S. Lee and S. D. Mahanti (to be published), to be published.
[101] G. Snyder, Nature materials 7, 105 (2008), and therein references.
[102] H. Grimm and A. Sommerfeld, Zeitschrift für Physik 36, 36 (1926).
[103] A. Rockett and R. W. Birkmire, Journal of Applied Physics 70, R81 (1991).
[104] S. Chichibu, T. Mizutani, K. Murakami, T. Shioda, T. Kurafuji, H. Nakanishi, S. Niki, P. J.
Fons, and A. Yamada, Journal of Applied Physics 83, 3678 (1998).
[105] E. J. Skoug, J. D. Cain, and D. T. Morelli, Journal of Alloys and Compounds 506, 18 (2010).
[106] E. Skoug, J. Cain, P. Majsztrik, M. Kirkham, E. Lara-Curzio, and D. Morelli, Science of
Advanced Materials 3, 602 (2011).
[107] E. J. Skoug and D. T. Morelli, Phys. Rev. Lett. 107, 235901 (2011).
[108] C. Sevik and T. Ça˘ in, Journal of Applied Physics 109, 123712 (2011).
g
[109] C. Yang, F. Huang, L. Wu, and K. Xu, Journal of Physics D: Applied Physics 44, 295404
(2011).
[110] H. Nakanishi, S. Endo, and T. Irie, Jpn. J. Appl. Phys. 8, 443 (1969).
[111] O. Madelung, Semiconductors: Data Handbook (Springer, 2004).
[112] L. Pauling, The Nature of Chemical Bond and the Structure of Molecules and Crystals (Cornell University Press, New York, 1939).
[113] J. C. Phillips, Bonds and Bands in Semiconductors (Academic Press, New York and London,
1973).
127

[114] G. N. Lewis, J. Am. Chem. Soc. 38, 762 (1916).
[115] I. Langmuir, J. Am. Chem. Soc. 41, 868 (1919).
[116] Y. J. Wang, H. Lin, T. Das, M. Z. Hasan, and a. Bansil, New Journal of Physics 13, 085017
(2011).
[117] J. Garin and E. Parthé, 28, 3672 (1972).
[118] A. Pfitzner and S. Reiser, 217, 51 (2002).
[119] A. Pfitzner and T. Bernert, 219, 20 (2004).
[120] A. Pfitzner, Z. Kristallogr. 209, 685 (1994).
[121] A. Pfitzner, Z. anorg. allg. Chem. 621, 685 (1995).
[122] M. S. Hybertsen, M. Schlüter, and N. E. Christensen, Phys. Rev. B 39, 9028 (1989).
[123] R. Laskowski, P. Blaha, and K. Schwarz, Phys. Rev. B 67, 075102 (2003).
[124] M. Papagno, D. Pacilé, G. Caimi, H. Berger, L. Degiorgi, and M. Grioni, Phys. Rev. B 73,
115120 (2006).
[125] A. Önsten, M. Månsson, T. Claesson, T. Muro, T. Matsushita, T. Nakamura, T. Kinoshita,
U. O. Karlsson, and O. Tjernberg, Phys. Rev. B 76, 115127 (2007).
[126] I. Galanakis, P. Mavropoulos, and P. H. Dederichs, J. Phys. D: Appl. Phys. 39, 765 (2006).
[127] I. Galanakis and P. Mavropoulos, J. Phys.: Condens. Matter 19, 315213 (2007).
[128] M. I. Katsnelson, V. Y. Irkhin, L. Chioncel, A. I. Lichtenstein, and R. A. de Groot, Rev.
Mod. Phys. 80, 315 (2008).
[129] S. Picozzi, in Advances in Solid State Physics, edited by P. D. R. Haug (Springer Berlin
Heidelberg, 2008), no. 47 in Advances in Solid State Physics, pp. 129–141, ISBN 978-3540-74324-8, 978-3-540-74325-5.
[130] T. Graf, S. S. P. Parkin, and C. Felser, IEEE Transactions on Magnetics 47, 367 (2011).
[131] F. Casper, T. Graf, S. Chadov, B. Balke, and C. Felser, Semiconductor Science and Technology 27, 063001 (2012).
[132] Z. Bai, L. Shen, G. Han, and Y. p. Feng, arXiv:1301.5455 (2013).
[133] I. Galanakis, arXiv:1302.4699 (2013).
[134] F. G. Aliev, N. Brandt, V. V. Kozyr’kov, V. V. Moshchalkov, R. V. Skolozdra, Y. V. Stadnyk,
and V. K. Pecharski, Pis’ma Zhr. Eksp. Teor. Fiz. 45, 535 (1987).
[135] F. G. Aliev, V. V. Kozyrkov, V. V. Moshchalkov, R. V. Scolozdra, and K. Durczewski,
Zeitschrift für Physik B Condensed Matter 80, 353 (1990).
128

[136] F. Aliev, Physica B: Condens. Matter 171, 199 (1991).
[137] F. G. Aliev, V. V. Pryadun, R. Villar, S. Vieira, J. M. Calleja, N. Mestres, and R. V. Scolozdra,
Int. J. Mod. Phys. B 07, 383 (1993).
[138] N. Mestres, J. Calleja, F. Aliev, and A. Belogorokhov, Solid State Commun. 91, 779 (1994).
[139] A. Slebarski, A. Jezierski, S. Lütkehoff, and M. Neumann, Phys. Rev. B 57, 6408 (1998).
[140] H. Hohl, A. P. Ramirez, C. Goldmann, G. Ernst, B. Wölfing, and E. Bucher, J. Phys.: Condens. Matter 11, 1697 (1999).
[141] P. Larson, S. D. Mahanti, and M. G. Kanatzidis, Phys. Rev. B 62, 12754 (2000).
[142] Q. Shen, L. Zhang, L. Chen, T. Goto, and T. Hirai, J. Mater. Sci. Lett. 20, 21972199 (2001).
[143] J. D. Germond, P. J. Schilling, N. J. Takas, and P. F. P. Poudeu, MRS Proc. 1267 (2010).
[144] Y. Kimura, T. Tanoguchi, and T. Kita, Acta Materialia 58, 4354 (2010).
[145] P. Qiu, J. Yang, X. Huang, X. Chen, and L. Chen, Appl. Phys. Lett. 96, 152105 (2010).
[146] Y. Kimura, T. Tanoguchi, Y. Sakai, Y.-W. Chai, and Y. Mishima, MRS Proc. 1295 (2011).
[147] M.-S. Lee and S. D. Mahanti, Phys. Rev. B 85, 165149 (2012).
[148] S. Chen, K. C. Lukas, W. Liu, C. P. Opeil, G. Chen, and Z. Ren, Adv. Energy Mater. (2013).
[149] V. Romaka, P. Rogl, L. Romaka, Y. Stadnyk, A. Grytsiv, O. Lakh, and V. Krayovskii, Intermetallics 35, 45 (2013).
[150] D. F. Zou, S. H. Xie, Y. Y. Liu, J. G. Lin, and J. Y. Li, J. Appl. Phys. 113, 193705 (2013).
[151] S. Sakurada and N. Shutoh, Appl. Phys. Lett. 86, 082105 (2005).
[152] L. Chaput, J. Tobola, P. Pécheur, and H. Scherrer, Phys. Rev. B 73, 045121 (2006).
[153] S.-W. Kim, Y. Kimura, and Y. Mishima, Intermetallics 15, 349 (2007).
[154] L. L. Wang, L. Miao, Z. Y. Wang, W. Wei, R. Xiong, H. J. Liu, J. Shi, and X. F. Tang, J.
Appl. Phys. 105, 013709 (2009).
[155] C. Yu, T.-J. Zhu, R.-Z. Shi, Y. Zhang, X.-B. Zhao, and J. He, Acta Materialia 57, 2757
(2009).
[156] H. Hazama, M. Matsubara, R. Asahi, and T. Takeuchi, Journal of Applied Physics 110,
063710 (2011).
[157] J. P. Makongo, D. K. Misra, J. R. Salvador, N. J. Takas, G. Wang, M. R. Shabetai, A. Pant,
P. Paudel, C. Uher, K. L. Stokes, et al., Journal of Solid State Chemistry 184, 2948 (2011).
[158] S. Poon, D. Wu, S. Zhu, W. Xie, T. Tritt, P. Thomas, and R. Venkatasubramanian, Journal
of Materials Research 26, 2795 (2011).
129

[159] J. E. Douglas, C. S. Birkel, M.-S. Miao, C. J. Torbet, G. D. Stucky, T. M. Pollock, and
R. Seshadri, App. Phys. Lett. 101, 183902 (2012).
[160] Y. W. Chai and Y. Kimura, Appl. Phys. Lett. 100, 033114 (2012).
[161] G. Joshi, T. Dahal, S. Chen, H. Wang, J. Shiomi, G. Chen, and Z. Ren, Nano Energy 2, 82
(2013).
[162] O. Appel, M. Schwall, D. Mogilyansky, M. Köhne, B. Balke, and Y. Gelbstein, Journal of
Electronic Materials 42, 1340 (2013).
[163] C. S. Birkel, J. E. Douglas, B. R. Lettiere, G. Seward, N. Verma, Y. Zhang, T. M. Pollock,
R. Seshadri, and G. D. Stucky, Physical Chemistry Chemical Physics 15, 6990 (2013).
[164] K. Kirievsky, Y. Gelbstein, and D. Fuks, Journal of Solid State Chemistry 203, 247 (2013).
[165] X. Yan, W. Liu, S. Chen, H. Wang, Q. Zhang, G. Chen, and Z. Ren, Adv. Energy Mater.
(2013).
[166] R. A. Downie, D. A. MacLaren, R. I. Smith, and J. W. G. Bos, Chem. Commun. 49, 4184
(2013).
[167] Y. Stadnyk, L. Romaka, A. Horyn, A. Tkachuk, Y. Gorelenko, and P. Rogl, Journal of Alloys
and Compounds 387, 251 (2005).
[168] H. Akai, Journal of Physics: Condensed Matter 1, 8045 (1989).
[169] R. M. Nieminen, Topics in Applied Physics: Theory of defects in semiconductors (Springer,
2006), vol. 104, pp. 36–40.
[170] S. Ö˘ üt and K. M. Rabe, Phys. Rev. B 51, 10443 (1995).
g
[171] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
[172] S. B. Zhang, S.-H. Wei, A. Zunger, and H. Katayama-Yoshida, Phys. Rev. B 57, 9642
(1998).
[173] P. Larson, S. D. Mahanti, S. Sportouch, and M. G. Kanatzidis, Phys. Rev. B 59, 15660
(1999).

130