Thermal conductivity in PbTe from first

principles

A. H. Romero, E. K. U. Gross, M. J. Verstraete and Olle Hellman

Linköping University Post Print

N.B.: When citing this work, cite the original article.

Original Publication:
A. H. Romero, E. K. U. Gross, M. J. Verstraete and Olle Hellman, Thermal conductivity in
PbTe from first principles, 2015, Physical Review B. Condensed Matter and Materials Physics,
(91), 21, 214310.
http://dx.doi.org/10.1103/PhysRevB.91.214310
Copyright: American Physical Society
http://www.aps.org/

Postprint available at: Linköping University Electronic Press

http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-120224
 PHYSICAL REVIEW B 91, 214310 (2015)

Thermal conductivity in PbTe from first principles

A. H. Romero,1,2 E. K. U. Gross,2 M. J. Verstraete,3 and Olle Hellman4,5
1
Department of Physics, West Virginia University, 207 White Hall, 26506, West Virginia, USA
2
Max-Planck-Institute für Mikrostrukturphysik, Weinberg 2, D-06120 Halle, Germany
3
Department of Physics and European Theoretical Spectroscopy Facility, Université de Liège, av du 6 août, 17, B-4000 Liège, Belgium
4
Department of Applied Physics and Materials Science, California Institute of Technology, Pasadena, California 91125, USA
5
Department of Physics, Chemistry and Biology (IFM), Linköping University, SE-581 83, Linköping, Sweden
(Received 20 April 2015; published 29 June 2015)
We investigate the harmonic and anharmonic contributions to the phonon spectrum of lead telluride and
perform a complete characterization of how thermal properties of PbTe evolve as temperature increases. We
analyze the thermal resistivity’s variation with temperature and clarify misconceptions about existing experimental
literature. The resistivity initially increases sublinearly because of phase space effects and ultra strong anharmonic
renormalizations of specific bands. This effect is the strongest factor in the favorable thermoelectric properties of
PbTe, and it explains its limitations at higher T . This quantitative prediction opens the prospect of phonon phase
space engineering to tailor the lifetimes of crucial heat carrying phonons by considering different structure or
nanostructure geometries. We analyze the available scattering volume between TO and LA phonons as a function
of temperature and correlate its changes to features in the thermal conductivity.

DOI: 10.1103/PhysRevB.91.214310 PACS number(s): 31.15.E−, 61.48.−c, 71.20.Tx, 78.30.Na

Heat conversion by using thermoelectric power generation ambient conditions and makes it such a good starting point
has received a huge amount of interest in the last few years: for nanostructuring and doping. For the crystalline system, the
transforming a temperature gradient to a voltage difference low thermal conductivity has been correlated to the presence
promises to recover waste heat in thermal engines, transform- of large anharmonic effects at the  point and zone boundary,
ing it into electrical energy. The thermoelectric efficiency of a large change in the Gruneissen parameters (in particular of a
a material is captured by the figure of merit, ZT = T S 2 σ/κ, soft TO mode) and a very low speed of sound [5,10–17]. The
where T is the temperature, S is the Seebeck coefficient, and σ electronic structure is intimately linked to the resonant bonding
and κ are the electrical and thermal conductivities. The search and vibrational properties, as shown by Lee et al. [18]. It is
for good thermoelectrics is centered on finding materials with a also important to note that recent champion thermoelectric
high figure of merit, which implies large electric and small ther- materials, such as SnSe, are related to PbTe, and it is expected
mal conductivities. Since the 1990s, a sequence of new materi- that they share the same strong anharmonicity [13]. An
als which offer new paradigms in this field and a large number understanding of the anharmonic effects in PbTe will clarify
of energy harvesting applications have been proposed [1–4]. many of the issues raised in other similar thermoelectrics.
Lead telluride is an industrial standard and a reference high In the following we spotlight the fundamental reason for
performance thermoelectric, considered to be one of the mate- which PbTe has a large thermal resistivity, how its anharmonic-
rials to beat with alternative paradigms such as nanostructuring ity goes beyond the second order response in the interatomic
or nanoalloying. The properties of PbTe in its halite structure force constants, and how the optical-acoustic branch crossing
have been well documented, and it has been the basis of a is removed with increasing temperature. The importance of
large set of investigations, both in its pristine and alloyed full anharmonicities, to be able to describe neutron scattering
configurations. The value of ZT for the crystalline phase of data up to room temperature, has been underscored by Li et al.
PbTe is close to 1.4 [5]. Recent investigations have reported [13] using the same methods. Additionally Lee et al. [18]
values above 1.5 which can be obtained by n or p doping, have pointed out that resonant bonding creates an anisotropy
and the possibility of band gap engineering has been recently in the interatomic force constants, in particular along the
explored [2,5]. This value can be further increased by alloying (100) direction, which corresponds to the TO mode and is
or nanostructuring. For example, Hsu et al. recently reported a responsible for a large part of the anharmonicity. They have
large figure of merit, close to 2.2 at 800 K, for nanostructured also correlated the phonon lifetimes with the available volume
PbTe [6]. Alloying and nanostructuring minimizes the lattice for three phonon scattering. We also observe these effects, now
thermal conductivity, and one must try to keep the electrical with infinite order anharmonic renormalization of frequencies
conductivity unaffected. A different approach to increasing the (essential for PbTe), and we quantify the available phonon
figure of merit is to enhance the electrical conductivity and S scattering volume.
by engineering the electronic band structure [7–9]. The use We demonstrate here the strong dependence of the thermal
of PbTe as a starting material is the common factor in a long resistivity with respect to volume, which converts to a depen-
series of such improvements of ZT . dence with respect to thermal expansion. For this quantity we
In order to understand the behavior of composite nanos- also clarify differences with previously reported results, where
tructured or alloyed materials based on PbTe, it is important good comparison was claimed, but with the wrong experimen-
to be able to explain the thermoelectric properties of pristine tal data. Additionally, we also observe strong differences in
PbTe, and, in particular, its low thermal conductivity between the Grüneisen parameters obtained from the quasiharmonic
450 and 800 K, which is the basis for its high efficiency at approximation compared with fully anharmonic methods. This

1098-0121/2015/91(21)/214310(7) 214310-1 ©2015 American Physical Society

ROMERO, GROSS, VERSTRAETE, AND HELLMAN PHYSICAL REVIEW B 91, 214310 (2015)

is a difference that has been previously discussed in the vector q. Through perturbation theory [22], one can include
references above: We show here that a fair comparison of the phonon self energy qs + iqs where:
these two methodologies, under the same conditions, exposes 
the pathologies of the full QHA self energy (its use is often 18   qq q 2
qs ( ) = − 2 dq dq  ss  s  
restricted to the imaginary part).  s  s  BZ
There have been numerous theoretical studies on the 
lattice dynamics of PbTe. Recent papers have focused on nq s  + nq s  + 1 nq s  + nq s  + 1
× +
demonstrating the sensitivity of the optical phonons to changes ωq s  + ωq s  + ωq s  + ωq s  −
in volume [12] and temperature [13,14,19,20]. Experimentally, 
nq s  − nq s  nq s  − nq s 
the phonon spectrum shows a strong dependence on tem- + +
perature. The coupling between the acoustic and transverse ωq s  − ωq s  + ωq s  − ωq s  −

optical modes is identified to strongly correlate with the 12  qq̄q q̄
thermal conductivity [19,20]. A clear theoretical explanation + dq sss  s  (2nq s  + 1) (2)
 s  BZ
of the anomalous stiffening of the LO mode with increasing  π    qq q 2
temperature was only given last year by Li et al. [13]. At qs ( ) = dq dq  ss  s  
the same time, there is a very strong anharmonicity between s  s 
16 BZ
the ferroelectric transverse mode and the longitudinal acoustic
×[(nq s  + nq s  + 1)δ( − ωq s  − ωq s  )
mode, which was studied in part by Li et al. in a followup
paper [21]. +2(nq s  − nq s  )δ( − ωq s  + ωq s  )], (3)
The anharmonic influence issue had already been raised, in
particular in Refs. [10,12], and the large difference between
predicted quasiharmonic approximation (QHA) results and 20

experimental values was pointed out. Delaire et al. concluded (a)

that the anharmonicity comes mostly from higher order terms
in the interatomic force constants, which affect the longitudinal 16
acoustic phonons and the heat they carry. One of the most
important features reported in PbTe in its B1 phase is the
Frequency (meV)

anharmonic behavior of the LA-TO coupling along the -X 12

direction, and in particular the crossing occurring around 1/3
along that path. In Ref. [10] and then in Ref. [13], inelastic neu-
tron scattering data has been reported, where the temperature
8
effects are demonstrated to strongly affect the zone boundary
by hardening the TO modes, such that the crossing between
LA and TO is lifted as a function of temperature. In these
same references the quasiharmonic then fully anharmonic 4

band structures demonstrated the importance of the latter.

Below we also use the TDEP method to analyze in detail the
pathologies of the quasiharmonic approximation, the higher 0
Γ X Γ L K
temperature phonons, and in particular the dramatic effect the
renormalization has on the thermal resistivity. It is important to 7

note that the method used to describe the anharmonic thermal (b)
behavior in this paper goes beyond the first few terms in 6

the expansion of the interatomic force constants, and that it

Mode Grünesien parameter

captures the complete anharmonic behavior to all orders. 5

At zero temperature, it is common practice to describe

the system with a model Hamiltonian as a Taylor expansion 4
of the full ion-electron system as a function of the ionic
displacements: 3

 p2 1   αβ α β 2
H = i
+  u u
i
2mi 2 ij αβ ij i j
1
1  αβγ α β γ 1  αβγ δ α β γ δ
+ ij k ui uj uk +  u u u u . (1)
3! ij k 4! ij kl ij kl i j k l 0

Here p indicates momentum and ui the displacement of atom Γ X Γ L K

i. The indices ij kl run over all atoms, and αβγ δ over Cartesian
components. The tensors  are the interatomic force constants FIG. 1. (Color online) Phonon dispersion relation and Mode
of increasing order. From the second order expansion, the Grünesien parameters for T = 0 K for a chosen path within the
frequencies ωqs are obtained as a function of mode s and wave Brillouin zone.

214310-2
THERMAL CONDUCTIVITY IN PbTe FROM FIRST . . . PHYSICAL REVIEW B 91, 214310 (2015)

10 (a) Potential at 0K (b) Potential at finite T

Harmonic Delaire et al
9
Quasiharmonic Cochran et al
8
Quasiharmonic+Δ TDEP
Energy of TO(Γ) (meV)

3 Displacement Displacement
2

1
FIG. 3. (Color online) Cartoon of the potential felt by a given
0 100 200 300 400 500 600 700 800 atom using the two different approaches to describe the materials
Temperature (K) thermal behavior. Left: QHA—the potential at a given volume is
fixed for all T . Right: TDEP—the ions see a different average
FIG. 2. (Color online) The energy of the TO() mode as a potential energy surface depending on the temperature and the present
function of temperature. The harmonic line has no temperature anharmonicity.
dependence at all, the quasiharmonic has temperature dependence
through the volume, the quasiharmonic+ is with the anharmonic
shifts, which increases the disagreement. The experimental points are observed in Fig. 4. While QHA predicts a fast growth of this
from Refs. [10] (circles) and [23] (triangles).
factor with temperature, as can be observed for most of the
X and  modes, this enhancement does not reproduce the
which depends on the interatomic force constants to third and observed experimental phonon hardening, and the TO mode
qq q qq q q
fourth order ( ss  s  and ss  s  s  ) and describe the lifetime evolution is qualitatively wrong. Our anharmonic calculations
broadening and shifts of the phonon frequencies. predict a slow decay of the Grünesien intensities as temperature
In Fig. 1 we show the phonon dispersion relations and increases, reproducing the saturation of all modes, in particular
mode Grüneisen parameters at 0 K for PbTe. The calculated the TO in Fig. 2.
Grüneisen parameters in PbTe are all positive except for As discussed, when the system is strongly anharmonic, a
the TA modes. When the temperature is increased, and the model Hamiltonian, at T = 0, obtained up to second order
crystal expands, the optical modes should all soften in the is not sufficient, and the Born-Oppenheimer energy surface
quasiharmonic approximation, as reported in An et al. [12]. has a nontrivial temperature dependence. To incorporate this
In Fig. 2 we show the energy of the TO modes at  as a in a lattice dynamical model, we seek explicitly temperature-
function of temperature. In the harmonic approximation, they dependent force constants. We have performed extensive ab
are naturally constant. In the quasiharmonic approximation, initio molecular dynamics (AIMD) to sample the potential
where (T ) = (V (T )) and V (T ) is the volume as a function energy surface. Using the temperature-dependent effective
of temperature, the TO() mode softens rapidly, as indicated potential technique (TDEP), we extract the temperature-
by the Grünesien parameters. If we add the anharmonic dependent interatomic force constants [24–26]. Since AIMD
frequency shift, qs (ωqs ), the mode becomes unstable around incorporates phonon-phonon and electron-phonon coupling
650 K. This is all contrary to experimental trends: the TO() implicitly to all orders, the effective force constants will
mode should stiffen with temperature, not soften. At the  contain that information. They represent the best possible fit
point several frequencies are seen in the neutron experiments, of the potential energy surface, for the current volume and
each corresponding to a peak in the scattering intensity. This temperature, using a finite order of force constants.
is not reproduced at all in the QHA. We employ Born-Oppenheimer molecular dynamics in
In case it was necessary, it is obvious by now that conven- density functional theory (DFT) with the projector-augmented
tional quasiharmonic theory can not adequately describe PbTe wave (PAW) method as implemented in the VASP code [27–30].
at elevated temperature. Including the T = 0 interatomic force To converge the force constants fully, a 5 × 5 × 5 supercell
constants up to fourth order does not help and even worsens (250 atoms) is employed. For the BZ integration we use the 
the agreement, as seen in Fig. 2. The observed potential point and run the simulations in the canonical ensemble, for a
energy well for an atom in the crystal is “hardened out” grid of temperatures and volumes. Temperature was controlled
as temperature increases, because it is now an average over using a Nosé thermostat [31]. Exchange-correlation effects
the potential energy created by the fluctuating atoms close were treated using the AM05 functional [32,33], and we use a
by, as is represented in Fig. 3. This behavior is not taken plane wave cutoff of 250 eV. The simulations run for about 30
into account by the QHA, which leads to its failure, and ps after equilibration with a time step of 1 fs, ensuring proper
is probably not accurately described using few anharmonic ergodicity. A subset of uncorrelated samples was then chosen,
terms in the energy expansion. A clear difference can be and for each of the samples the electronic structure and total

214310-3
ROMERO, GROSS, VERSTRAETE, AND HELLMAN PHYSICAL REVIEW B 91, 214310 (2015)

TDEP, T=42K TDEP, T=168K TDEP, T=294K TDEP, T=463K TDEP, T=631K
10
Mode Grüneisen parameter

-2
Γ X Γ L Γ X Γ L Γ X Γ L Γ X Γ L Γ X Γ L
QH, T=42K QH, T=168K QH, T=294K QH, T=463K QH, T=631K
20
Mode Grüneisen parameter

15

10

Γ X Γ L Γ X Γ L Γ X Γ L Γ X Γ L Γ X Γ L

FIG. 4. (Color online) Mode Grünesien parameters as a function of temperature in the TDEP and quasiharmonic (QH) formalisms. The
trends with increasing temperature are opposite, with the TDEP being the physically reasonable ones.

energies are recalculated using a 3 × 3 × 3 k-point grid and a A crossing of the experimental acoustic and optical
cutoff of 300 eV to obtain more accurate forces. branches in the -X direction is present at 300 K and
Recent experimental results [10,13] provide inelastic neu- disappears by 600 K. This is also reproduced by our calcu-
tron scattering spectra for PbTe as a function of temperature. lations, and is central to the good thermal performance of
In Ref. [13] Li et al. also show TDEP calculations at 50 K PbTe, as will be shown below. For higher temperatures PbTe
and room temperature, which are fully compatible with our loses this topological advantage when the crossing is lifted.
calculations at 300 K. Observe that in our figure we have used The quasiharmonic results do not reproduce any of these
a different path along the Brillouin zone, to illustrate more features, and at 600 K the self energy shifts have rendered the
clearly how the anharmonicity affects the LA/TO crossing. To system completely unstable. Previous perturbation theoretical
compare with these experiments we convolve our data with studies avoid this problem by disregarding either thermal
the inelastic neutron scattering cross section [34]: expansion or [19]. From the quantitative agreement with
experimental data, we conclude that we have built an accurate
2ωqs qs ( ) lattice dynamical model for PbTe at finite temperature with
σqs ( ) ∝  2 . (4)
2 − ω2 − 2ω
qs qs ( ) + 4ωqs qs ( )
2 2 a good description of thermal properties and temperature
qs
dependencies. We now apply it to the calculation of the lattice
Figure 5 shows the theoretical and experimental inelastic thermal resistivity ρ.
neutron scattering cross section. The Debye-Waller factor and Experimental data for the thermal resistivity of PbTe has
the experimental transfer functions have been disregarded. For been measured in Refs. [35] and [36] at low and high temper-
the calculations of the self energy, a 31 × 31 × 31 q-point ature, respectively. In the literature, Ref. [37] is often referred
grid and a Gaussian smearing of 0.1 meV have been used to for the bulk lattice thermal conductivity of PbTe. This is
for the numerical evaluation of  and . The calculated and not correct: The paper does not contain any values of lattice
experimental spectra show remarkable agreement in absolute thermal conductivity. The data attributed to Ref. [37] is actually
and relative shifts of the phonon frequencies. The TDEP that contained in Ref. [35]. More importantly, close inspection
accurately describes the double peak structure of the TO of Ref. [35] shows that what is usually reported for the lattice
mode at , as shown by Li et al. [13] at 300 K. Traces of thermal conductivity at high T is not a measurement but a
this double peak remain at 600 K, which also agrees with linear extrapolation from values obtained between 50 and 280
the experimental picture. The energies of the two peaks as a K. The good agreement with experimental ρ at high T in Refs.
function of temperature are shown with squares in Fig. 2, and [19,20] is thus probably the result of the analytically linear
they agree well with the experimental values. behavior of ρ in perturbation theory and the inappropriate

214310-4
THERMAL CONDUCTIVITY IN PbTe FROM FIRST . . . PHYSICAL REVIEW B 91, 214310 (2015)

Quasiharmonic, 300K TDEP, 300K Experiment, 300K

10 10 10

8 8 8

Energy (meV)

Energy (meV)
Energy (meV)

6 6 6

4 4 4

2 2 2

X Γ X X Γ X X Γ X

Quasiharmonic, 600K TDEP, 600K Experiment, 600K

10 10 10

8 8 8
Energy (meV)
Energy (meV)

Energy (meV)
6 6 6

4 4 4

2 2 2

X Γ X X Γ X X Γ X

FIG. 5. (Color online) PbTe Neutron scattering cross section. Notice the double peaks close to  and the lifting of the optical branches
above the acoustic captured within the TDEP method as compared to experiment. The quasiharmonic are included only for reference.

linear extrapolation in Ref. [35]. The data of El Sharkawy [36] more reasonable and are smaller than experiment: additional
demonstrate the highly nonlinear behavior of ρ(T ). possible scattering mechanisms (defects, grain boundaries)
We calculate the resistivity from the phonon Boltzmann must add to the resistivity. In particular the sublinear behavior
equation for fixed volume, beyond the relaxation time ap- observed up to 600 K is reproduced correctly. At higher
proximation and including isotope scattering as described in temperatures the expected linear behavior is restored, due to
Ref. [38]. The results are summarized in Fig. 7, using the TDEP the high T linear limit of the Bose Einstein distribution for
and quasiharmonic treatments, together with experimental each given frequency.
results and considering different volumes for the theoretical Additionally, 0 K perturbation theory becomes questionable
calculations. in the case of PbTe: For a few modes, the linewidth and
As temperature increases, the thermal resistivity increases frequency can be comparable and qs  ωqs definitely does
slower than linearly. This is strongly correlated with the not hold. What our results show is that the TDEP renor-
position of the TO mode with respect to the acoustic branches: malized quantities do obey the relation qs (T )  ωqs (T ).
The nested scattering between the LA and TO branches was Finally, the large Grüneisen parameters imply a strong volume
already pointed out by Ref. [19]. When the TO mode is lifted up dependence. In Fig. 7 we show the thermal resistivity for
with increasing temperature, the available scattering volume lattice parameters compressed and expanded by 1%, which
shrinks dramatically, the acoustic phonons are allowed to carry leads to a factor 2 difference in resistivity. The intrinsic
heat unperturbed, and the resistivity is sublinear with T. This exchange-correlation error in current DFT functionals is
process is illustrated in Fig. 6, where at low temperature there certainly of this order of magnitude for the volume, and more
are scattering channels available between the TO and LA. accurate calculations will necessitate further improvements in
When the TO mode stiffens with temperature, these channels the fundamentals of DFT.
no longer satisfy momentum and energy conservation and the Our results clearly demonstrate that a careful examination
lifetimes of acoustic phonons increase. of the phonon response is crucial when dealing with systems as
At 100 K both methods start in decent agreement with strongly anharmonic as PbTe. The TDEP method contains the
experiment. As T increases, however, the QHA results increase anharmonic terms necessary to describe phonon frequencies at
much faster, and the thermal expansion effect pushes the high T . We have noted that the crossing present in the phonon
results away from the experimental curves. TDEP values are dispersion relation is lifted by the temperature renormalization

214310-5
ROMERO, GROSS, VERSTRAETE, AND HELLMAN PHYSICAL REVIEW B 91, 214310 (2015)

(a) T=0K (b) T=150K (c) T=300K

FIG. 6. (Color online) The energy and momentum conserving surfaces for scattering between the TO mode and LA and TA phonons. As
temperature increases and the TO mode stiffens, the available scattering volume for this channel of decay disappears, which is the origin of the
kink in thermal resistivity. The isosurface represents the nesting factor for scattering between bands 2, 3, and 4, while conserving energy and
momentum, integrated over q  and q  .

of the interatomic force constants, with no need for additional can now be taken. Exploiting the volume (and hence stress)
perturbation theory. Quantitative prediction of the exception- dependence of the phonon frequencies, the phonon resistivity
ally low thermal resistivity is possible but still very challenging will be maximized by preserving the optical-acoustic cross-
in this delicate system, where large anharmonicities determine ings. The opposite effect—very different masses lead to a large
the physical behavior. An important effect is also observed phonon gap and impeded phonon scattering—was shown by
in the changes in the electronic structure as a function of Broido and co-workers, in BAs [39] among other compounds,
temperature, Fig. 8, where we have taken the average of the and with anomalous pressure effects [40]. Many other high
electron band structure for 20 uncorrelated configurations. It performance thermoelectrics present low-lying optical modes
is clear that, as the temperature changes, the electron potential and could be manipulated in similar ways: using epitaxy to start
surface also changes, and electron bands both shift and from a lower lattice constant at ambient or uniaxial stress to
broaden. This can have a double effect, in bringing the valence create a preferred direction for the heat flow, as could certainly
band valleys closer to degeneracy (see maxima between  and be done with SnSe [41,42], which is the bulk material with the
X or K, and at L), which will affect the electronic part of ZT . highest reported figure of merit.
Many successful but heuristic attempts have been made It is common to consider phonon band engineering, invoked
to further reduce the low thermal conductivity of PbTe (by in many nanostructuring or nanoalloying strategies, to reduce
doping, alloying, and nanostructuring). We show here that a the speed of sound and break the group velocity of acoustic
direct rational approach, based on first principles calculations, modes, e.g., using disorder or rattlers. Much less is quantified
about the reduction of phonon lifetimes by manipulating the
available phase space for decay, as we have demonstrated

1
Energy (eV)

-1

-2

FIG. 7. (Color online) Thermal resistivity of PbTe versus temper-

ature. The red lines are QHA results with T = 0 force constants and Γ X K Γ L
volumes corresponding to the experimental thermal expansion. The
green lines are calculated using temperature-dependent force con- FIG. 8. (Color online) Electronic band structure average over a
stants, with volumes following the experimental thermal expansion. molecular dynamics trajectory.

214310-6
THERMAL CONDUCTIVITY IN PbTe FROM FIRST . . . PHYSICAL REVIEW B 91, 214310 (2015)

happens in PbTe, confirming the work of Shiga et al. [19]. The coming fellowships (Grant No. PIIFR-GA-2011-911070) and
presence of anharmonicity is intimately linked to the mass of the Donors of the American Chemical Society Petroleum
the atoms and overall softness of the phonon modes, which are Research Fund for partial support of this research under
important features of most successful thermoelectrics. contract 54075-ND10. This material is also based upon work
Summarizing, we describe the relevance of high orders supported by the National Science Foundation under Grant
of the interatomic force constants in the proper description No. 1434897. M.J.V. acknowledges an Action de Recherches
of the thermoelectric properties of PbTe. This result should Concertées (ARC) grant (TheMoTherm # 10/15-03) from the
guide future work on alloy-based thermoelectrics as well as Communauté Française de Belgique, and computer time from
other high temperature applications. The dependence should CECI, SEGI, and PRACE-2IP and 3IP (EU FP7 Grant Nos.
also affect similar compounds, where strong anharmonicities RI-283493 and RI-312763) on Huygens, Hector, and Archer.
have been reported, and can not be accounted for by normal Support from the Swedish Research Council (VR) program
perturbation theories. 637-2013-7296 is gratefully acknowledged. Supercomputer
resources were provided by the Swedish National Infrastruc-
A.H.R. acknowledges the support of the Marie Curie
ture for Computing (SNIC).
Actions from the European Union in the international in-

[1] M. Zebarjadi, K. Esfarjani, M. S. Dresselhaus, Z. Ren, and G. [20] Z. Tian, J. Garg, K. Esfarjani, T. Shiga, J. Shiomi, and G. Chen,
Chen, Energy Environ. Sci. 5, 5147 (2012). Phys. Rev. B 85, 184303 (2012).
[2] A. D. LaLonde, Y. Pei, H. Wang, and G. Jeffrey Snyder, Mater. [21] C. W. Li, J. Ma, H. B. Cao, A. F. May, D. L. Abernathy, G.
Today 14, 526 (2011). Ehlers, C. Hoffmann, X. Wang, T. Hong, A. Huq, O. Gourdon,
[3] A. J. Minnich, M. S. Dresselhaus, Z. F. Ren, and G. Chen, Energy and O. Delaire, Phys. Rev. B 90, 214303 (2014).
Environ. Sci. 2, 466 (2009). [22] D. C. Wallace, Thermodynamics of Crystals, Dover Books on
[4] M. S. Dresselhaus, G. Chen, Z. Ren, G. Dresselhaus, A. Henry, Physics (Dover Publications Inc., Mineola, NY, 1998).
and J. P. Fleurial, JOM 61, 86 (2009). [23] W. Cochran, R. A. Cowley, G. Dolling, and M. M. Elcombe,
[5] A. D. LaLonde, Y. Pei, and G. J. Snyder, Energy Environ. Sci. Proc. R. Soc. London, Ser. A 293, 433 (1966).
4, 2090 (2011). [24] O. Hellman, I. A. Abrikosov, and S. I. Simak, Phys. Rev. B 84,
[6] K. F. Hsu, S. Loo, F. Guo, W. Chen, J. S. Dyck, C. Uher, 180301 (2011).
T. Hogan, E. K. Polychroniadis, and M. G. Kanatzidis, Science [25] O. Hellman, P. Steneteg, I. A. Abrikosov, and S. I. Simak, Phys.
(New York, NY) 303, 818 (2004). Rev. B 87, 104111 (2013).
[7] C. M. Jaworski, M. D. Nielsen, H. Wang, S. N. Girard, W. Cai, [26] O. Hellman and I. A. Abrikosov, Phys. Rev. B 88, 144301
W. D. Porter, M. G. Kanatzidis, and J. P. Heremans, Phys. Rev. (2013).
B 87, 045203 (2013). [27] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
[8] Y. Pei, H. Wang, and G. J. Snyder, Advanced materials (Deerfield [28] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169
Beach, FL) 24, 6125 (2012). (1996).
[9] Q. Zhang, F. Cao, W. Liu, K. Lukas, B. Yu, S. Chen, C. Opeil, [29] G. Kresse and J. Hafner, Phys. Rev. B 48, 13115 (1993).
D. A. Broido, G. Chen, and Z. Ren, J. Am. Chem. Soc. 134, [30] G. Kresse, Comput. Mater. Sci. 6, 15 (1996).
10031 (2012). [31] S. Nosé, Mol. Phys. 52, 255 (1984).
[10] O. Delaire, J. Ma, K. Marty, A. May, M. McGuire, M.-H. Du, [32] R. Armiento and A. E. Mattsson, Phys. Rev. B 72, 085108
D. Singh, A. Podlesnyak, G. Ehlers, M. Lumsden, and B. Sales, (2005).
Nat. Mater. 10, 614 (2011). [33] A. E. Mattsson and R. Armiento, Phys. Rev. B 79, 155101
[11] D. J. Singh, Sci. Adv. Mater. 3, 561 (2011). (2009).
[12] J. An, A. Subedi, and D. J. Singh, Solid State Commun. 148, [34] R. A. Cowley, Rep. Prog. Phys. 31, 123 (1968).
417 (2008). [35] E. D. Devyatkova and I. A. Smirnov, Sov. Phys. Solid State,
[13] C. W. Li, O. Hellman, J. Ma, A. F. May, H. B. Cao, X. Chen, USSR 4, 2507 (1962).
A. D. Christianson, G. Ehlers, D. J. Singh, B. C. Sales, and O. [36] A. A. El-Sharkawy, A. M. Abou El-Azm, M. I. Kenawy, A. S.
Delaire, Phys. Rev. Lett. 112, 175501 (2014). Hillal, and H. M. Abu-Basha, Int. J. Thermophys. 4, 261 (1983).
[14] C.-L. Chen, H. Wang, Y.-Y. Chen, T. Day, and G. J. Snyder, J. [37] V. I. Fedorov and V. I. Machuev, Sov. Phys. Solid State, USSR
Mater. Chem. A 2, 11171 (2014). 11, 1116 (1969).
[15] J. Al-Otaibi and G. P. Srivastava, J. Appl. Phys. 116, 043702 [38] D. A. Broido, M. Malorny, G. Birner, N. Mingo, and D. A.
(2014). Stewart, Appl. Phys. Lett. 91, 231922 (2007).
[16] T. Shiga, T. Murakami, T. Hori, O. Delaire, and J. Shiomi, Appl. [39] D. A. Broido, L. Lindsay, and T. L. Reinecke, Phys. Rev. B 88,
Phys. Express 7, 041801 (2014). 214303 (2013).
[17] Y. Chen, X. Ai, and C.A. Marianetti, Phys. Rev. Lett. 113, [40] L. Lindsay, D. A. Broido, J. Carrete, N. Mingo, and T. L.
105501 (2014). Reinecke, Phys. Rev. B 91, 121202 (2015).
[18] S. Lee, K. Esfarjani, T. Luo, J. Zhou, Z. Tian, and G. Chen, Nat. [41] L.-D. Zhao, S.-H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher,
Commun. 5, 3525 (2014). C. Wolverton, V. P. Dravid, and M. G. Kanatzidis, Nature
[19] T. Shiga, J. Shiomi, J. Ma, O. Delaire, T. Radzynski, A. (London) 508, 373 (2014).
Lusakowski, K. Esfarjani, and G. Chen, Phys. Rev. B 85, 155203 [42] J. Carrete, N. Mingo, and S. Curtarolo, Appl. Phys. Lett. 105,
(2012). 101907 (2014).

214310-7