Table Of ContentFirst principles search for n-type oxide, nitride, and sulfide thermoelectrics
Kevin F. Garrity∗
Material Measurement Laboratory, National Institute of Standards and Technology, Gaithersburg MD, 20899
(Dated: June 24, 2016)
Oxides have many potentially desirable characteristics for thermoelectric applications, including
low cost and stability at high temperatures, but thus far there are few known high zT n-type
oxide thermoelectrics. In this work, we use high-throughput first principles calculations to screen
transitionmetaloxides,nitrides,andsulfidesforcandidatematerialswithhighpowerfactorsandlow
thermal conductivity. We find a variety of promising materials, and we investigate these materials
indetailinordertounderstandthemechanismsthatcausethemtohavehighpowerfactors. These
6 materialsallcombineahighdensityofstatesneartheFermilevelwithdispersivebands,reducingthe
1 trade-off between the Seebeck coefficient and the electrical conductivity, but they do so for several
0 different reasons. In addition, our calculations indicate that many of our candidate materials have
2 low thermal conductivity.
n
u
I. INTRODUCTION oxide thermoelectrics has focused on a relatively small
J
number of oxides, mostly binaries and perovskites, leav-
2
ingopenthepossibilitythatbetteroxidethermoelectrics
2 The need for clean efficient power generation has led
to a renewed interest in thermoelectric materials, which exist.
] In this work, we use high-throughput density func-
can directly convert a temperature gradient into electri-
i
c cal power. Thermoelectrics can take advantage of a va- tional theory (DFT) calculations[12–23] to identify
s promising n-type thermoelectric oxides and related ma-
- riety of heat sources, including solar or waste heat, to
l cleanly generate electricity[1–5]. Conversely, they could terials from the Inorganic Crystal Structure Database
r
t be used in cooling applications via the Peltier effect. (ICSD)[24]. The large amount of work required to syn-
m
thesize, optimize, and measure thermoelectrics experi-
There has been an extensive effort over recent years to
. mentallymakethistypeoftheoreticalscreeningofcandi-
t discover and optimize materials with high zT, a dimen-
a date materials particularly desirable. Similar techniques
sionless thermoelectric figure of merit. While there has
m have been used successfully to study the thermoelectric
been significant progress in this area, existing materials
- have not yet managed to provide a combination of high behavior of a variety of materials, including oxides[25–
d
36]. While a fully first principles theoretical calculation
n zT, low materials cost, and high durability that would
of zT remains challenging, especially for oxides, which
o result in widespread adoption. Much of the research on
often have partially localized carriers, we can neverthe-
c thermoelectricshasfocusedonhighmobilitysemiconduc-
[ tors with small band gaps. Unfortunately, many of the less screen materials for both electronic and vibrational
2 most promising candidate materials have practical con- propertiesthatarenecessaryforgoodthermoelectricper-
formance. Inthiswork,weperformsuchascreeningpro-
v cerns(cost,toxicity,stability)whichhavethusfarlimited
cedure,identifyingmanycandidatematerialswithcalcu-
2 their use in applications[6].
2 latedthermoelectricpropertiesthataresimilartoorsur-
In this work, we focus on the less explored group of
6 passexperimentallystudiedn-typeoxides. Furthermore,
wide band gap transition metal oxides, as well as re-
1 we analyze the mechanisms behind the high thermoelec-
lated nitrides and sulfides. While oxides are not usually
0 tricperformanceofthesematerials,findingthattheyfall
. thought of as promising for thermoelectric applications,
1 into a small number of groups with similar properties.
due to their typically low mobilities, the discovery of
0
goodthermoelectricperformanceinp-typeNa CoO and
6 x 2
1 otherlayeredCo-basedmaterialsresultedinanincreased
II. METHODS
: interest in this class of materials[1–3, 7]. n-type mate-
v
rials such as ZnO and SrTiO have also displayed high
i 3
X powerfactors,buttheirzT valueshavethusfarbeenonly A. Calculating Thermoelectric Performance
r moderate, due to high thermal conductivities[8–10]. De-
a spitelimitedsuccessthusfar,oxidesprovidemanypoten- The dimensionless figure of merit for thermoelectrics
tial advantages as thermoelectrics: 1) high thermal and can be written as
chemical stability in air, 2) chemical versatility, allowing
for extensive substitutions and doping, 3) low thermal zT =σS2T/κ, (1)
conductivity,and4)lowcostmaterialsandprocessing[1–
3, 6, 11]. Thus far, much of the experimental work on where σ is the electrical conductivity, S is the Seebeck
coefficient, κ is the total thermal conductivity (electrical
pluslattice),andT isthetemperature. Thepowerfactor,
which determines the electrical response of a material to
∗Electronicaddress: [email protected] a temperature gradient, is S2σ.
2
Unfortunately, the components of zT are not all easy calculating band structure related quantities. We find
to calculate using first principles techniques. Within the that for most materials this correction has a relatively
constant relaxation time approximation, which is used minor effect beyond increasing the band gap, and larger
in this work, S can be calculated from a band struc- gapshavenodirecteffecton thermoelectric performance
ture calculation without any adjustable parameters[37]. as long as the gap is already large enough to avoid
Within the same approximation, it is possible to cal- significant thermal carrier excitation. We perform
culate σ/τ , where τ is the electronic relaxation time. phonon calculations using DFT perturbation theory[51]
e e
Unfortunately, calculating τ from first principles re- without the +U correction.
e
mains challenging[38, 39]. This problem is especially Our main results are done on fully relaxed struc-
severe for oxides, which often display complicated con- tures with initial coordinates from the ICSD. We use
duction mechanisms and polaronic effects at low dop- PYMATGEN[52] to manipulate files from the ICSD to
ing and low temperatures. In this work, we are con- setuptheinitialstructuresforrelaxation. Becausecalcu-
cerned primarily with the opposite regime of high tem- lationsarerunatafixednumberofplanewaves,changes
peratures and high doping, where the carrier mobilities in the unit cell during relaxation can effectively modify
of oxides are typically larger[1–3]. Because we are com- the basis set. To ensure consistency between the ba-
paring materials which are chemically similar, we expect sis set and the final structure, we run each relaxation
them to have broadly similar electron scattering mech- three times, with a force convergence tolerance of 0.001
anisms. Therefore, we will use the quantity S2σ/τe to Ry/Bohr, an energy tolerance of 1 × 10−4 Ry, and a
rank our candidate materials for suitability as thermo- stress tolerance of 0.5 Kbar. For phonon calculations,
electrics. This estimate, which has been used in many we decrease the force tolerance to 5 × 10−5 Ry/Bohr.
previousworks[26,34,40],shouldbesufficienttoatleast The BFGS algorithm as implemented in QUANTUM
screen materials for those with band structures that are ESPRESSO was used for relaxations.
promising for thermoelectric applications, even if deter-
We use maximally localized Wannier functions as im-
mining the final ranking of materials will require experi-
plemented in WANNIER90[53–55] to interpolate band
mental input.
structures and BOLTZWANN, the WANNIER90 trans-
Forreference,firstprinciplestechniquescanreproduce
portmodule,tocalculatetheSeebeckcoefficientandcon-
the thermoelectric properties of SrTiO with τ ≈ 4 fs
3 e ductivity under the relaxation time approximation[56,
at room temperature[33], a typical value for oxides, but
57]. The use of Wannier interpolation allows us to per-
some high mobility oxides like ZnO have much longer
form accurate calculations of thermoelectric quantities
scattering times[41]. All wide band materials have to be
starting from relatively sparse first principles k-point
doped in order to be used as thermoelectrics. In this
grids, which we then interpolate to a k-point spacing of
work, we use the rigid band filling model to estimate the 0.02 ˚A−1. This density is about ten times as dense as
effects of doping, and we rank materials by S2σ/τ at
e the first principles calculation along each direction in k-
their optimum doping.
space. The use of Wannier functions also allows us to
Afteridentifyingmaterialswithpromisingbandstruc-
calculate band structure derivatives analytically, which
tures, we perform more computationally expensive
accuratelytreatsdegeneratepointsintheBrillouinzone.
phonon calculations for a limited number of candidate
InordertouseWannierinterpolationforthiswork,we
materials to estimate the lattice thermal conductiv-
had to develop a procedure for automating the construc-
ity, which is the dominant contribution to the thermal
tion of localized Wannier functions. Because we are in-
conductivity for most thermoelectrics, as described in
terestedinthepropertiesofbothvalenceandconduction
Sec. IIC.
states, we normally include all possible orbitals which
could contribute to states near the Fermi level (see sup-
plemental materials for list). This is in contrast to many
B. Band Structure Calculations applications of Wannier functions, which are concerned
with either only the occupied bands or only a localized
All of our calculations are based on DFT subspace of bands (e.g. d-orbitals). In these cases, the
calculations[42, 43], as implemented in QUANTUM Wannier functions extend over several atoms and may
ESPRESSO[44] and using the GBRV high-throughput be sensitive to the details of the localization procedure.
ultrasoft pseudopotential library[45]. We use a plane For our application, we include all the relevant orbitals,
wave cutoff of 40 Ryd for band structure calculations which results in Wannier functions that are atomic-like
and 45-50 Ryd for phonon calculations. For Brillouin and strongly localized, even before the iterative localiza-
zoneintegration, weuseaΓ-centeredgridwithadensity tion procedure, making the final result more robust. To
of 1500 k-points per atom. calculate the Wannier functions, we use an inner ’frozen’
We use the PBEsol exchange-correlation window default of 4.5 eV above and below the conduc-
functional[46], which provides more accurate lattice tion/valence band edges in order to ensure an accurate
constants and phonon frequencies than other GGA interpolation of the band structure. In testing, our cal-
functionals. We use the DFT+U technique[47–49], with culatedthermoelectricpropertiesareinsensitivetominor
a U value of 3 eV for transition metal d-states[50], when variations in this window.
3
WhiletheWannierizationprocedureweoutlinedabove tothemodelinBjerget. al.,wheretheGruneisenparam-
is relatively robust, there are a few situations that can eter (γ) and the Debye temperature (Θ )are calculated
D
resultinfailuresintheWannierization, whichareidenti- from the first principles phonon dispersion
fied by monitoring the spread of the Wannier functions.
First,iftherearesemicorestatesthatwereexcludedfrom (cid:115)5h¯2 (cid:82)∞ω2g(ω)dω
Θ = n−1/3 0 (2)
the Wannierization that overlap in energy with the va- D 3k2 (cid:82)∞g(ω)dω
lence states, it will be necessary to include those states B 0
in the valence. Second, sometimes there are problems γ2 = (cid:80)i(cid:82) 8dπq3γi2qCiq (3)
including orbitals with high energy (e.g. Sr d-states in (cid:80) (cid:82) dq C
i 8π3 iq
SrTiO ),asthesestatescanbecomedifficulttodisentan-
3 V ∂ω
gle from the free electron-like bands when their energy γ = − iq (4)
iq ω ∂V
becomes too high. In both of these cases, we simply ad- iq
just the orbitals that we include in the Wannierization
where n is the number of atoms per unit cell, ω is the
iq
procedure by hand to fix these problems. Another is-
angularfrequencyofphononmodeiatq-pointq,g(ω)is
sue can arise if the ’frozen’ window overlaps with free
the phonon density of states, γ is the mode Gruneisen
iq
electron-like bands. This can be fixed by adjusting this
parameter, C is the mode specific heat, and V is the
iq
window downward to avoid overlap. We encountered all
volume. The sum for the Gruneisen parameter is only
of these problematic cases only rarely, and we adjust for
performed over modes with ¯hω < k Θ . As per the
iq B D
them when necessary.
discussion in Refs. 63, 69, we square γ to avoid cancel-
iq
One potential drawback of the Wannierization ap-
lationbetweenpositiveandnegativeanharmonicitywhen
proach is the necessity of including a large number of
calculating γ.
emptybandsinanon-self-consistentDFTcalculation,in
Using the Debye temperature and Gruneisen parame-
order to construct well-localized conduction band Wan-
ter calculated above, we then insert them into the Slack
nier functions. However, these extra bands are only re-
model [63, 67, 68, 70, 71]:
quiredonthethesparsek-pointgrid,andinpracticethe
√
computationalcostofthisstepissmallerthantheinitial 0.849×334
κ (T)=
structural relaxation. l 20π3(1−0.514γ−1+0.228γ−2)
×(cid:16)kBΘD(cid:17)2 kBMV 31 ΘD (5)
¯h ¯hγ2 T
C. Thermal conductivity
whereM istheaverageatomicmass. Wefindthatthe
For typical thermoelectrics, the thermal conductivity DebyetemperatureandGruneisenparameterusedinthis
is dominated by the lattice thermal conductivity (κl). waycontainalmostalloftheinformationofthefullBjerg
First principles calculations of the thermal conductiv- model, as shown in table II, which presents correlations
ity have been shown to be accurate for a wide variety of various models and quantities with our reference set
ofmaterials[58–62]. Unfortunately,thesecalculationsre- of thermal conductivities. In fact, as shown in table II,
quiretheanharmonicforceconstants,whicharetoocom- the quantity Θ /γ also has a high rank correlation with
D
putationally expensive to use as an initial screening tool thereferencethermalconductivities,althoughthereisno
for high-throughput calculations, especially as many of computational advantage in using Θ /γ instead of Eq.
D
the materials we consider have large unit cells with rela- 5.
tively low symmetry. The reason we do not use the full Bjerg model is that
There have been various recent attempts to model the in some cases we found difficulty in fully converging the
lattice thermal conductivity without performing a full acoustic modes for large unit cells. These modes depend
calculationoftheanharmonicforceconstants[31,62–66]. on careful cancellation between all of the force constants
Yan et. al. use a Debye-Callaway model with a con- to produce modes with zero eigenvalues at q =Γ, which
stant Gruneisen parameter[31, 66]. Toher et. al.[65] is challenging to achieve numerically. This cancellation
demonstrated that a modeled the Debye temperature can be enforced at Γ by using the acoustic sum rule to
and the Gruneisen parameter, combined using the Slack modify the force constant matrix in various ways, but
model[67, 68] for thermal conductivity, is useful as a we sometimes found results which depend on how the
screening method for thermal conductivity. Another rule was enforced. Therefore, we opted to use a more
screening method by Bjerg et. al.[63, 69] incorporates computationallyrobustprocedureappropriateforahigh-
aspects of the first principles phonon band structure to throughputstudybyusingtheSlackmodelinsteadofthe
approximate the lattice thermal conductivity. Bjerg model.
In this work, we want a method which is both com- We find empirically in testing that using this combi-
putationally feasible to apply to a few dozen compounds nation of the Slack model with the Bjerg definition of
to use as a secondary screening procedure, and accurate theDebyetemperatureandGruneisenparameteroveres-
enoughtoprovideareasonableorderingofcompoundsto timatesthethermalconductivity,sowereport70%ofthe
consider for further study. We employ a method similar modelvalue,whichimprovesthequantitativeaccuracyin
4
principles thermal conductivities for a variety of half-
Heusler compounds [62, 65, 72]. In addition, we com-
pare with experimental thermal conductivities for a few
) additional oxides[73]. The results are shown graphically
K
0 in Fig. 1 and correlations are given in table II, see the
0
3 supplementary materials for details. We find that our
at chosen method is sufficient for screening materials for
-1K those likely to have low thermal conductivity. As shown
1 in table II, the Spearman rank correlation[65] between
-m
W thereferenceandmodeledthermalconductivitiesforthe
entire test set is 0.91, indicating we are able to identify
(
ck promising materials. If we limit the dataset to materi-
a
sl als with κ < 50 W m−1K−1 at 300 K, a more realistic
κ l
rangeforcomplexoxides,therankcorrelationforthefull
model drops to 0.83, which is still reasonable for select-
ing materials to study further. We also note that when
κ (W m-1 K-1 at 300 K)
expt considering the entire dataset, the Debye temperature
(see Eq. 2) alone has a rank correlation of 0.82 with the
FIG.1: Comparisonofreferenceexperimentalandfirstprin- thermalconductivity,andwemakeitasaninitialscreen-
ciplesthermalconductivities(x-axis)andtheSlackmodel(y- ing tool, as it is less computationally expensive than the
axis). full model. However, directly calculating the Gruneisen
parameter, rather than estimating it[31, 65, 66], signifi-
cantly increases the accuracy of the model. Finally, ex-
TABLE I: Correlations of various quantities with the ref-
perimental thermal conductivities are sensitive to many
erence thermal conductivities at 300 K (see Supplementary
factors beyond the scope of this work, including defects
Materials for list). κ is the model used in this work (see
Slack and grain boundaries, which both makes comparisons
Eq. 5) ; Bjerg refers to the full model of Ref. [63]. Pear-
with experiments difficult but increases the possibility
son and Spearman refer to the standard Pearson correlation
of engineering materials to have lower thermal conduc-
and the Spearman rank correlation, respectively. The first
twocolumnsincludeallmaterials,thenexttwoarelimitedto tivities.
materials with κ <50 W m−1K−1 at 300 K.
l
Quantity Pearson Spearman Pearson Spearman
D. Effective Masses
Low κ Low κ
l l
κ 0.83 0.91 0.65 0.83
Slack
Bjerg[63] 0.93 0.92 0.69 0.88 In order to understand the conductivity and Seebeck
Θ 0.74 0.82 0.43 0.71 coefficient, we consider several definitions of the effec-
D
1/γ 0.39 0.75 0.66 0.60 tivemass. UsingourWannierinterpolation,wecalculate
ΘD/γ 0.71 0.89 0.60 0.83 derivatives of the band structure analytically[56] at the
conductionbandminimatofindtheeffectivemasstensor
(m )−1 = 1 d2E . We will sometimes concentrate on
ij h¯2dkidkj
our testing but makes no difference in a ranking of com- mmin, the smallest eigenvalue of mij, which helps deter-
pounds for those with the lowest thermal conductivity. mine the largest value of the conductivity tensor at low
For materials with unstable phonon modes, we cannot temperature. We will also consider miso =(m1m2m3)31,
calculateaGruneisenparameterinameaningfulwayus- the isotropic effective mass, where m1, m2, and m3 are
ing purely harmonic calculations, as the unstable modes the eigenvalues of mij.
must be stabilized anharmonically at finite temperature. A related band structure descriptor we calculate is a
Therefore, we do not estimate a thermal conductivity versionoftheeffectivemassbasedonthedensityofstates
for those compounds. In practice, we expect many ma- (DOS):[31, 74, 75]
terials with unstable modes at zero temperature to have m (E) = ¯h2(cid:112)3 π4g(E)g(cid:48)(E) (6)
lowthermalconductivityduetoanharmonicinteractions, DOS
so the observation of unstable modes is already a useful (cid:82) dEg(E)m (E)(−df )
m (T,n ) = DOS dE . (7)
indicator of anharmonicity. We present our calculated DOS d dEg(E)(−df )
thermal conductivities at 300 K, even though we expect dE
thesematerialstobeusedathighertemperatures,where In this expression, g(E) is the DOS at energy E, f(E)
the thermal conductivity will be lower. is the Fermi function, and T and n are the temper-
d
We establish the validity of this method for screen- ature and doping. This definition of m matches
DOS
ing the thermal conductivity by comparing the model m for a single parabolic band, but it is higher for
eff
with the experimental thermal conductivities of a vari- non-parabolic bands or when multiple bands contribute
ety simple binary semiconductors, as well as the first to the conduction. These features allow m to give
DOS
5
a good description of the Seebeck coefficient for many S2σ/τ
materials[31, 74, 75].
E. Materials selection and screening procedure s)
f
m
We are interested in discovering new n-type thermo- /
S
electric oxides, nitrides, or sulfides. As there are over
3
80,000 entries with oxygen in the ICSD, there is a need -0
1
to significantly limit our search space before proceeding. (
Previous experimental and theoretical work on thermo- τ
/
electricshassuggestedthatgoodthermoelectricstendto σ
have anisotropic and non-parabolic bands and high den-
sities of states, all of which can be created by empty
d-orbitals[33, 40, 76–78]. Furthermore, materials with
emptydorbitalscanusuallybedopedn-type,withsome
|S| (μV/K)
carriersoccurringnaturallyduetooxygenvacancies[1–3].
Therefore,inthiswork,wefocusonmaterialscontaining
FIG. 2: |S|(µV/K) versus σ/τ (103 S ) for the entire
at least one of Y, Sc, Ti, Zr, Hf, Nb, Ta, Mo, or W as e mfs
dataset, at T=700K and doping at 1021 cm−3. The color
well as at least one of O, N, or S. In order to limit com-
scale indicates size of power factor, S2σ/τ .
putational time, we restrict our search set to structures e
with primitive unit cell volumes of less than 300 ˚A3.
Our screening procedure proceeds in several steps, a) b)
with each step becoming more computationally expen-
sinivgem. Fatiresrti,alwse(7c6a6lcucolamtepothuendbsa)nadngdaepliomfianlaltoeutrhoeumrsettaarlts-. V/K) S / m fs)
FSeoerbmecaktecroiaelffiswciietnht,aSg,apan(d59e2lemctartoenriicalcs)o,nwdeucctailvciutyla,tσe/tτhee. S (μ -3 τ (10
For materials with a high power factor (191 materials), σ /
we then calculate the Debye temperature using Eq. 2.
Ftoirnaalnlyd,/foorr alowsubDseebtyoef 1te9mmpaetreartiuarlse,wwitehchailgchulpaoteweκrSflaacck- c) m D O S d) 1/mmin
l
withEq. 5. ResultsarepresentedintableIIwithfurther
details in the supplementary materials.
τ
σ/ τ
2S σ/
2S
III. RESULTS AND DISCUSSION
We performed our screening procedure starting with
766 compounds from the ICSD as discussed in Sec. IIE, DOS at E (m )2/m
F DOS min
consistingof661oxides,60nitrides,and71sulfides(some
compoundscontainmultipleanions). Ofthatlist,wefind FIG. 3: Thermoelectric properties of the entire dataset, at
592 materials with band gaps according to DFT+U (551 fixedT=700Knd=1021cm−3. a)mDOSversusS,b)1/mmin
oxides, 53 nitrides, and 25 sulfides). For these materials, versus σ/τe, where mmin is the smallest component of mij,
c) DOS at E versus the power factor, and d) m2 /m
wecalculateS andσ/τe foravarietyoftemperaturesand F DOS min
versus the power factor.
dopings.
If we sort this list of candidate materials by esti-
mated power factor, S2σ/τ , at 700 K and with opti-
e
mizeddoping,werediscoverseveralcompoundsthathave perior to existing materials. In table II, we list some of
previously been measured to have good n-type thermo- our most promising candidate materials, including those
electric properties. For example, doped TiO2, SrTiO3, with high S2σ/τe and a few with moderate S2σ/τe and
KTaO3,andTiS2haveallbeenmeasuredtohavepromis- low κl. We remove structures that are minor distortions
ing power factors and show up highly in our list[1– of other structures on the list, that have missing atoms
3,9,10,33,35,79–82]. Thisgivesusconfidencethatour in the ICSD, or that are only theoretically proposaled
screening procedure is useful. In addition to these pre- structures; full results are presented in the Supplemen-
viously measured materials, there are a variety of com- tary Materials.
pounds which have not been studied for thermoelectric We begin our analysis by looking for patterns in the
applications and which may have properties that are su- entire dataset. First, we note that under the rigid band
6
TABLE II: Thermoelectric properties of the most promising compounds. The first five columns consist of the compound
name, its space group number, the DFT+U band gap (eV), the isotropic effective mass. The next several columns are the
DOS,m (seeEq. 6-7),S,andσ/τ ,allevaluatedatat700Kandfixed1021cm−3doping. ThenextcolumnisS2σ/τ atthe
DOS e e
optimizeddopingandat700K,andfollowingcolumnisthatoptimizeddoping. Thelastcolumnisthemodeledlatticethermal
conductivity at 300 K. Materials with − as κ have unstable phonon modes or we were unable to converge the Gruneisen
l
parameter.
Material Space Grp. Band Gap m DOS m S σ/τ S2σ/τ Opt. Doping κ
iso DOS e e l
(eV) (103 eV−1˚A−3) (µV/K) 10−3 S 103 W (1021cm−3) W
mfs mK2fs mK
CaTaAlO 15 4.0 2.8 4.08 5.0 -189 42.9 1.7 3 −
5
TiO 225 1.7 1.1 1.53 17.0 -305 9.6 1.7 8 −
2
LiNbO 161 3.5 2.1 3.92 7.0 -245 23.4 1.6 2 40
3
TiO 136 2.1 2.2 3.63 7.5 -250 20.7 1.6 3 <1
2
HfS 164 1.6 1.6 5.51 2.7 -165 56.9 1.5 1 3.4
2
NaNbO 63 1.8 3.8 1.89 0.8 -268 75.6 1.4 4 −
3
Ba TaInO 225 4.3 2.3 2.75 10.0 -285 12 1.4 5 −
2 6
YClO 129 5.1 1.1 3.20 6.6 -182 32.2 1.4 3 6.0
LiTaO 161 3.8 2.9 5.50 3.4 -204 36.3 1.4 1 34
3
Li ZrN 164 1.9 0.6 5.10 3.5 -178 42.5 1.4 2 19
2 2
CaTiSiO 15 3.2 4.0 4.28 6.4 -229 25.5 1.3 2 1.1
5
HgWO 15 2.4 1.5 4.22 3.5 -171 46.5 1.3 2 −
4
P WO 12 2.3 11.8 6.01 3.5 -150 56.1 1.3 0.7 6.5
2 8
ZrS 164 1.1 1.6 4.88 3.3 -172 45.1 1.3 2 22
2
TaPO 85 3.7 2.6 4.38 4.1 -160 49.5 1.3 2 −
5
LaTaO 36 3.4 2.9 5.13 3.7 -183 38.1 1.3 2 26
4
NbTl S 217 2.3 0.9 3.59 8.6 -246 19.2 1.3 2 −
3 4
SrTaNO 140 0.9 0.7 3.31 1.8 -109 122.6 1.3 0.3 −
2
TiS 164 0.4 0.4 4.73 3.7 -171 42 1.3 2 5.3
2
PbTiO 99 2.0 0.6 3.68 1.7 -489 35.3 1.2 7 16
3
Sr TaInO 225 4.3 2.1 2.65 10.0 -266 10.4 1.2 6 −
2 6
SrTiO 140 1.2 1.3 3.56 1.6 -165 39.5 1.2 6 −
3
NaNbN 166 1.1 0.8 3.70 3.3 -142 42 1.1 4 −
2
HfTaNO 25 2.0 0.8 4.22 2.4 -122 73 1.1 0.6 −
3
KNbO 99 1.6 3.0 1.82 1.3 -367 75.5 1.1 3 −
3
HfSiO 141 5.8 2.7 4.36 6.3 -276 18.8 1.1 2 190
4
Y O 164 4.2 1.0 3.44 3.7 -122 41.8 1.1 4 5.8
2 3
ZrO 225 3.7 1.1 3.78 4.6 -145 18.7 1.1 10 −
2
CdTiO 62 2.5 1.0 5.2 4.4 -205 20.7 0.9 2 <1
3
CaTiO 62 2.7 1.1 5.23 3.5 -170 26.8 0.8 1 <1
3
Y Ti O 227 3.1 1.1 2.21 10.7 -224 11.4 0.6 2 <1
2 2 7
model used in this work, most materials have optimal the following sections, we explore in more detail how
dopings of about 1021 cm−3, which corresponds to dop- some of these individual materials achieve this higher
ings on the order to 10%. While this is much higher than expected combination of S and σ/τ .
e
than typical semiconductor thermoelectrics, it is consis- Thetrade-offbetweenS andσ/τ makesfindingasim-
e
tent with the behavior of oxides like SrTiO3[1–3, 9, 10]. ple descriptor of the power factor in terms of features of
Reaching such high doping values may be difficult in the band structure very difficult, even though we can
practice, and will require further experimental and theo- relateσ andS individuallytofeaturesinthebandstruc-
retical work (see for instance [83]). In this work we con- ture. IntableIII,wepresentSpearmanrankcorrelations
centrate on identifying promising materials for further betweenseveralthermoelectricpropertiesandvariousde-
optimization. scriptors ofthe band structure, and several ofthese rela-
In Fig. 2, we plot the values of S versus σ/τ , at 700K tionships are plotted in Fig. 3.
e
and for a fixed doping of 1021 cm−3. The color scale For example, we find that the smallest value of the ef-
shows the value of S2σ/τ [86]. There is a clear trade-off fectivemasstensor,m ,ishighlycorrelatedwithσ/τ ,
e min e
between S and σ/τ , which is consistent with the be- and we plot this relationship in Fig. 3b. Unsurprisingly,
e
havior of simple parabolic bands where S ∝ m and materials with small effective masses usually have high
eff
σ ∝ 1/m , where m is the effective mass[4]. The conductivities, although this relationship can be compli-
eff eff
best materials do not maximize either S or σ/τ , but in- cated by anisotropy in the effective mass tensor or by
e
stead have S and σ/τ values in the center of observed many bands contributing to the conduction. Similarly,
e
range, but with a larger combination than is typical. In as shown in Fig. 3a, we find that we can model the See-
7
3
TABLE III: Spearman rank correlation matrix of various
band structure descriptors and thermoelectric quantities at
2
700Kandatafixeddopingof1021cm−3. S2σ/τ isthepower
e
factor,|S|istheabsolutevalueoftheSeebeckcoefficient,σ/τ
e ) 1
is the electrical conductivity, mmin is the minimum value of eV
mDOijS,tihsetheffeedcetinvseitmyaosfsstteantesosra,tmtDheOSFeirsmdiefilenveedl.inCoErqr.el7a,tiaonnds gy ( 0
r
e
with absolute value above 0.65 are in bold. n
E 1
Quantity S2σ/τ |S| σ/τ m m DOS
e e min DOS
S2σ/τ − 0.44 0.07 -0.03 0.22 0.67
e 2
|S| 0.44 − -0.81 0.66 0.84 0.43
σ/τ 0.07 -0.81 − -0.81 -0.83 -0.21
e 3
m -0.03 0.66 -0.81 − 0.74 0.28 R Γ X M Γ
min
m 0.22 0.84 -0.83 0.74 − 0.32
DOS FIG. 4: Band structure of cubic SrTiO . Energies are rela-
DOS 0.67 0.43 -0.21 0.28 0.32 − 3
tive to the Fermi level.
beck coefficient as S ∝m (T,n ) (see Eq. 7). three classes of materials, although some materials fall
DOS d
Despite these relatively strong relationships for S and into several classes. In general, the mechanisms for high
σ in terms of certain definitions of the effective mass, powerfactorsconsistofcombiningalargenumberofflat
we find that the combination of m2 /m has only a bands near the conduction band minimum with at least
DOS min
weakcorrelationwithS2σ/τ ,asplottedinFig. 3d. The somehighlydispersivebands. Thiscombinationallowsa
e
problem is that as shown in table III, the two definitions large number of carriers, some of which are in dispersive
of the effective mass are strongly correlated with each bands, increasing σ, while at the same time keeping the
other, and dividing one by another does not produce a Fermi level near the conduction band minimum, where
usefuldescriptor. Otherquantitieslikem havesimilar |S| is largest, mitigating the typical trade-off.
iso
problems. Furthermore,alloftheeffectivemasses,aswell
as S and σ/τ , are strongly correlated or anti-correlated
e
witheachother,butnonearebythemselvesstronglycor- A. Symmetry driven degeneracy
related with the power factor. This can be understood
in part by looking at Fig. 2, which shows that the best This group of promising thermoelectrics contains ma-
materials do not lie at the extreme of either S or σ/τ , terials that have symmetries (or near symmetries) which
e
buttheyinsteadhaveanatypicalrelationshipbetweenS cause degeneracies in the conduction band minimum.
and σ/τ . Finding a simple descriptor of that relation- This degeneracy increases the DOS for any given Fermi
e
ship is difficult when the stronger trend is the trade-off level, relative to a material without degeneracies. These
between S and σ/τ . In addition, in many anisotropic degeneracies can be due to a single degenerate minimum
e
materials, the Seebeck coefficient and conductivity are inBrillouinzone, ortheycanbedueto abandstructure
not maximized in the same direction, which makes find- with a conduction band minimum which is repeated due
ing a simple descriptor for the maximum power factor to symmetry.
more difficult. Finally, as we will see below, all of our In addition to having degeneracies, the conduction
best materials have unusual band structures with some bands of these materials all consist of empty tran-
combination of high anisotropy, non-parabolic behavior, sition metal d-orbitals that have highly anisotropic
and multiple bands contributing to the conduction, and dispersions[33,35,40,78]. Theseanisotropicbandstruc-
identifyingthesequalitiesrequiresgoingbeyondasimple tures allow the material to have both low m and
eff
effective mass description of parabolic bands. high m bands at the same minimum, combining large
eff
The best simple descriptor we found for S2σ/τ does Seebeck coefficients with the high conductivity. Simi-
e
not include any effective mass, but instead is just the lar degeneracies and anisotropic bands are behind the
DOS evaluated at the relevant doping and temperature, high power factors of several semiconducting materi-
as shown in Fig. 3c. While there is a clear relationship als which rely of empty p-orbitals instead of empty d-
between the DOS and the power factor, many high DOS orbitals[27, 31, 40, 78].
materials have low power factors, making a high DOS a All of these features are present in the band struc-
useful design criterion but not a sufficient condition for ture of cubic SrTiO , as shown in see Fig. 4, which is
3
good thermoelectric performance. knowntobeagoodn-typethermoelectric. SrTiO hasa
3
In the following sub-sections, we will investigate the single triply degenerate conduction band minimum at Γ
band structures of some of the materials in table II, in due to the t states originating from the Ti-d orbitals.
2g
ordertoevaluatethemechanismsthatallowthesepartic- These bands have highly anisotropic dispersions, with
ular materials to minimize the trade-off between S and one nearly flat band (m = 6.3) and two highly dis-
eff
σ. We find that these materials separate roughly into persive bands (m = 0.4) going from Γ to X. The
eff
8
a)
a)
b) c) Al
b) c)
Hf O
S Ta
Ca
FIG. 5: a) Band structure of HfS . b-c) Side and top views FIG. 6: a) Band structure of CaTaAlO . b-c) Top and
2 5
of HfS . Larger gray atoms are Hf, smaller yellow atoms are sideviewsofCaTaAlO . LargeyellowatomsareCa,medium
S 5
S. gray atoms are Ta, smaller magenta atoms are Al, smallest
red atoms are O.
combination of high degeneracy and high effective mass
bands with very low effective mass bands, which allow by reducing its dimensionality and therefore increasing
for high conductivity, is what allows SrTiO to escape its DOS is well-known, and has been shown in SrTiO
3 3
the normal trade-off between S and σ/τ . superlattices[76, 85]. We note that here we are consider-
e
Similar features are present in many of the other per- ing thermodynamically stable materials, rather than ar-
ovskite variants which we find to be candidate thermo- tificial superlattices, nanowires, or quantum dots, which
electrics (SrTiO , PbTiO , NaNbO , LiNbO , KNbO , should reduce manufacturing costs and increase thermo-
3 3 3 3 3
LiTaO ,Ba TaInO ,CaTiO ,Sr TaInO ,SrTaNO ). In dynamic stability.
3 2 6 3 2 6 2
addition, various phases of TiO and ZrO have similar Wepresenttwoexamplesofeffectivelylowdimensional
2 2
featureswhichleadtohighpowerfactors. Manyofthese materials which we predict have high power factors.
materials have been studied as thermoelectrics before, First, in Fig. 5, we show the band structure and atomic
and the mechanisms leading to their power factors are structure of HfS , which consists of weakly bound two-
2
relatively well-known[27, 31, 33, 35], so we will proceed dimensional hexagonal trilayers. The conduction bands
with a discussion of the next two groups. areveryflatfromM totheminimumatL(m =4.5),
eff
characteristic of two-dimensional materials, but they are
much more dispersive in other directions (m =0.3).
eff
B. Low dimensional conductors Second, in Fig. 6, we show the band structure and
atomic structure of CaTaAlO , which consists of TaO
5 6
While all of the structures studied in this work are octahedra arranged into one-dimensional columns that
three dimensional, in many cases the atoms which dom- are separated from each other by Ca ions and AlO
4
inate the conduction band minima are arranged in tetrahedra. This arrangement of Ta atoms leads to an
two-dimensional layers, one-dimensional lines, or zero- anisotropic band structure with very flat bands from Γ
dimensional dots, which leads to effectively low dimen- to Y (m = 2.4) but stronger dispersion from Γ to A
eff
sional conduction. In some cases, the material itself con- (m = 0.5). There are additional nearly degenerate
eff
sists of weakly bound layers, while in others there are conduction band minima at Y which also contribute to
strong bonds in all three directions, but the transition the DOS. In HfS and CaTaAlO , the high DOS and the
5
metals are arranged in a low-dimensional way. strong anisotropy, which are caused by the low dimen-
Reducing the effective dimensionality of a material sionality, create the conditions for a high power factor.
results in highly anisotropic conduction bands and an Within our set of candidate thermoelectrics, ZrS ,
2
increased DOS at the bottom of bands, which can in- TiS , HfS , YClO, CaTiSiO , WP O , TaPO , NaNbN
2 2 5 2 8 5 2
crease the power factor[77, 78, 84]. The idea of im- have quasi-two-dimensional structures, CaTaAlO ,
5
proving the power factor of a candidate thermoelectric HgWO , LaTaO , and HfSiO have quasi-one-
4 4 4
9
3 metald-bands. Inaddition,inbothYClOandY O ,the
2 3
emptysandd-statesoftheYatomsarebothlocatednear
2
the conduction band minimum, which results in similar
1 behavior to the case where the orbitals come from dif-
)
V 0 ferent atoms. Depending on the crystal structure and
e
gy ( 1 tfrhoemanCioun,sZ,nit, mAga,yCbde,pAouss,iborleHtogeonrgpin-beaenrdesmfprotymsI-bna,nSdns,
r
e
n 2 Tl,Pb,orBitobecomedegeneratewithtransitionmetal
E
bands in this fashion. This type of engineering could al-
3
low one material to take advantage of the high Seebeck
4 coefficientsoftransitionmetaloxideswhileincorporating
the higher mobility of semiconductors, which often have
5
L Γ X K Γ empty s or p orbitals from main group elements. The
exact alignments of empty states from different atoms
FIG. 7: Band structure of Sr TaInO . Bands with greater
2 6 is difficult to predict using DFT+U, so further study of
then 35% In content are colored red, others are black.
these materials to determine the band alignments more
precisely may be necessary.
dimensional structures, and NbTl S , Ba TaInO ,
3 4 2 6
and Sr TaInO have quasi-zero-dimensional structures,
2 6
as their transition metals are separated from each other. D. Thermal Conductivity
There are other possible advantages in using low-
dimensional materials as thermoelectrics besides the in- Due to the high computational cost, we were not able
creased DOS, including potentially lower thermal con- to calculate the thermal conductivity of our full dataset.
ductivity, due to phonon scattering from the atomic For 191 compounds, we calculated the Debye tempera-
layers, as well as the ability to physically separate ture (see Eq. 2), which is fairly strongly correlated with
dopantsfromconductingchannels,whichcanreduceelec- thermal conductivity (see table II and supplementary
tron scattering. One disadvantage is that the thermo- materials). Calculations of the Debye temperature are
electric properties of low-dimensional materials will be bothlesssensitivetotheq-pointsamplingofthephonon
anisotropic, resulting in reduced efficiency in polycrys- bandstructurethantheGruneisenparameterandrequire
talline samples. the phonons at only one volume, making the computa-
tions much faster.
We find that in our set of transition metal oxides, ni-
C. Accidental Degeneracies trides, and sulfides, there is relatively little variation in
the Debye temperature (mean of 342 K, standard devia-
Onefinalmechanismforincreasingthepowerfactorof tion of 66 K), as compared to our test dataset of simple
an n-type oxide is to find or engineer a material with ac- binary and ternary semiconductors (mean 319 K, stan-
cidental degeneracies of the conduction band minimum. dard deviation 234 K). This is likely due to the fact that
Whilethiscanhappenforphysicallysimilarbandswhich all of our compounds contain ionic bonds between light
happen to be degenerate at different points in the Bril- anionsandmediumtoheavytransitionmetals,whilethe
louin zone, here we consider cases where the bands come test dataset contains a range of bonds, from covalent to
fromdifferentorbitalsandhavedifferenteffectivemasses. ionic, and a range of atom masses. In both datasets,
For example, in the double perovskite Sr TaInO , the there is a significant correlation between V−1 and the
2 6
conduction band consists of both Ta d-states and In p- Debye temperature, with a correlation coefficient of 0.68
states,whichhappentobeatsimilarenergies(seeFig.7, in the oxides, and 0.88 in the test set. This suggests
whichhighlightstheInstatesinred). TheInstateshave that looking at oxides with larger unit cells could be
loweffectivemasses(m =0.2−0.3)whiletheTastates beneficial[31, 66].
eff
have much higher effective masses (m =13−62), al- DuetotherelativelyweakvariationintheDebyetem-
eff
lowing the material to take advantage of both types of peraturethroughoutoursetofoxides, theGruneisenpa-
bands, inadditiontotheincreasedDOSprovidedbythe rameter becomes more important to identify the materi-
near degeneracy. The large effective masses of the Ta- als likely to have low thermal conductivity. Many of the
d bands are caused by small overlap between them and oxides we consider have soft or unstable phonon modes,
many of the neighboring In-p orbitals, which results in whichlikelyresultsinstronganharmonicityandlowther-
very flat bands, which contribute to a high DOS and malconductivity, butthisisdifficulttoquantifywithout
high power factor. more involved calculations. Due to the large computa-
Similarmaterialswithtwodifferentatomscontributing tional cost, we are only able to calculate the Gruneisen
totheconductionareHgWO ,Ba TaInO andNbTl S , parameter of the materials in table II. We do not have a
4 2 6 3 4
where the Hg(+2), In(+3) and Tl(+1) ions, contribute large enough database of oxide Gruneisen parameters to
empty s/p-bands at similar energies to the transition identify any trends which would predict which materials
10
will have soft modes without doing phonon calculations. effective mass, respectively, of the conduction band elec-
As shown in table II, many of the materials we have trons. However, to find materials with high power fac-
identifiedashavingpromisingpowerfactorsalsohavelow tors, it is necessary to look for materials which are not
thermalconductivityaccordingtoourmodel. Mostofthe well described by a single parabolic band, but instead
perovskite materials we study have strongly anharmonic have degeneracy, anisotropy, or other features which re-
modes, which leads to relatively low thermal conductiv- sult in a high density of states combined with disper-
ities, both in our calculations and in experiment[73]. In sive bands at the Fermi level. These materials achieve
addition, we find that many of the materials with one or their high power factors due to some combination of
twodimensionalbondingalsohavesoftmodes,likelydue symmetry-enforced degeneracies, low dimensionality, or
to the fact that many of the atoms are relatively free to accidental degeneracies. In addition, we use phonon cal-
vibrate in at least one direction. culations to model the thermal conductivity of our best
candidates, and we find many that have low lattice ther-
malconductivityorthatrequireanharmonicstabilization
IV. CONCLUSIONS of the harmonic modes. We hope further work on these
materials, as well as the understanding gained by exam-
We have used high throughput first principles calcu- iningthemechanismswhichleadtohighpowerfactorsin
lations to search for n-type transition metal oxides, ni- oxides,willleadtoimprovedthermoelectricperformance
trides, and sulfides which are promising for thermoelec- in oxides.
tricapplications. Wefindmanymaterialswithestimated
power factors which are comparable to or surpass previ-
ously studied oxide thermoelectrics. Acknowledgments
Across the entire sample of compounds, we find the
expected correlations between the Seebeck coeffient and We wish to acknowledge discussions with Igor Levin
electricalconductivitywiththeeffectivemassandinverse andhelpwiththeICSDfromVickyKarenandXiangLi.
[1] S.Walia,S.Balendhran,H.Nili,S.Zhuiykov,G.Rosen- cher,T.Mueller,K.A.Persson,andG.Ceder,Comput.
garten, Q. Wang, M. Bhaskaran, S. Sriram, M. Strano, Mater. Sci. 50, 2295 (2011), ISSN 0927-0256.
andK.Kalantar-zadeh,Prog.Mat.Sci.58,1443(2013). [13] A. Jain, G. Hautier, C. J. Moore, S. P. Ong, C. C. Fis-
[2] M. Backhaus-Ricoult, J. Rustad, L. Moore, C. Smith, cher,T.Mueller,K.A.Persson,andG.Ceder,Comput.
and J. Brown, Appl. Phys. A 116, 443 (2014). Mater. Sci. 50, 2295 (2011).
[3] J.He,Y.Liu,andR.Funahashi,J.Mater.Res.26,1762 [14] G. Ceder, Y. Chiang, D. Sadoway, M. Aydinol, Y. Jang,
(2011). and B. Huang, Nature 392, 694 (1998).
[4] M.Zebarjadi,K.Esfarjani,M.S.Dresselhaus,Z.F.Ren, [15] A. R. Akbarzadeh, V. Ozolins, and C. Wolverton, Adv.
and G. Chen, Energy Environ. Sci. 5, 5147 (2012), URL Mater. 19, 3233 (2007), ISSN 1521-4095.
http://dx.doi.org/10.1039/C1EE02497C. [16] D. Morgan, G. Ceder, and S. Curtarolo, Meas. Sci. &
[5] S. Riffat and X. Ma, Applied Thermal Engi- Tech. 16, 296 (2005), ISSN 0957-0233.
neering 23, 913 (2003), ISSN 1359-4311, URL [17] O.Levy,G.L.W.Hart,andS.Curtarolo,J.Am.Chem.
http://www.sciencedirect.com/science/article/ Soc. 132, 4830 (2010).
pii/S1359431103000127. [18] G. H. Johannesson, T. Bligaard, A. V. Ruban, H. L.
[6] M. W. Gaultois, T. D. Sparks, C. K. H. Skriver, K. W. Jacobsen, and J. K. Norskov, Phys. Rev.
Borg, R. Seshadri, W. D. Bonificio, and D. R. Lett. 88, 255506 (2002).
Clarke, Chemistry of Materials 25, 2911 (2013), [19] S.Curtarolo,D.Morgan,andG.Ceder,Calphad29,163
http://dx.doi.org/10.1021/cm400893e, URL (2005), ISSN 0364-5916.
http://dx.doi.org/10.1021/cm400893e. [20] W.SetyawanaandS.Curtarolo,Comput.Mater.Sci.49,
[7] I.Terasaki,Y.Sasago,andK.Uchinokura,Phys.Rev.B 299 (2010).
56, R12685 (1997). [21] X. Zhang, L. Yu, A. Zakutayev, and A. Zunger, Adv.
[8] M.Ohtaki,T.Tsubota,K.Eguchi,andH.Arai,Journal Funct. Mater. 22, 1425 (2012).
of Applied Physics 79 (1996). [22] J. W. Bennett and K. M. Rabe, J. Solid State Chem.
[9] S. Ohta, T. Nomura, H. Ohta, M. Hirano, H. Hosono, 195, 21 (2012).
and K. Koumoto, Applied Physics Letters 87, 092108 [23] S.Curtarolo,G.Hart,M.BuongiornoNardelli,N.Mingo,
(2005),URLhttp://scitation.aip.org/content/aip/ S.Sanvito,andO.Levy,NatureMaterials12,191(2013).
journal/apl/87/9/10.1063/1.2035889. [24] A.Belsky,M.Hellenbrandt,V.L.Karen,andP.Luksch,
[10] T. Okuda, K. Nakanishi, S. Miyasaka, and Y. Tokura, Acta Crystallographica Section B 58, 364 (2002), URL
Phys.Rev.B63,113104(2001),URLhttp://link.aps. http://dx.doi.org/10.1107/S0108768102006948.
org/doi/10.1103/PhysRevB.63.113104. [25] S.Wang,Z.Wang,W.Setyawan,N.Mingo,andS.Cur-
[11] W. Pan, S. Phillpot, C. Wan, A. Chernatynskiy, and tarolo, Phys. Rev. X 1, 021012 (2011).
Z. Qu, MRS Bulletin 37, 917 (2012). [26] I. Opahle, A. Parma, E. J. McEniry, R. Drautz, and
[12] A. Jain, G. Hautier, C. J. Moore, S. P. Ong, C. C. Fis- G. K. Madsen, New Journal of Physics 15, 105010
