Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles calculation of thermophysical properties ofsolids with strong phonon anharmonicity
Terumasa TADANO
Research Center for Magnetic and Spintronic Materials,
National Institute for Materials Science (NIMS),
Sengen, Tsukuba, Ibaraki 305-0047

Abstract
Optimizing thermophysical properties of solids, including thermal expansion, thermal conductivity,
and thermodynamic stability is crucial for developing energy- or light-harvesting material that showsa high device performance at operating temperatures. To achieve a simulation-guided design ofsuch materials, lattice vibration (phonon excitation)needs to be modeled accurately beyond the quasiharmonic level. In recent years, several important theoretical and technical improvements havebeen made in the field of ab initio phonon calculation. These methods have overcome the limitation of the conventional phonon calculation method
and further strengthened the predictive power offirst-principles calculations based on density functional theory. In this report, we review some of
the new computational approaches and demonstrate
their validity through applications to a hydrogenrich superconductor and thermoelectric materials
in which phonons are strongly anharmonic.

1

Introduction

Predicting and understanding thermophysical properties of solids only from the crystal structure is
crucial for achieving a simulation-guided design ofhigh-performance energy- and light-harvesting materials including but not limited to thermoelectric,
photovoltaic, and superconducting materials. For
example, lattice thermal conductivity, κL , needs tobe reduced to achieve a high thermoelectric figure-of-merit ZT . Also, the critical temperature (Tc ) of
ferroelectric phase transition plays a crucial role inrealizing a giant dielectric constant of polar semiconductors at ambient temperature, as exemplified
by the commercially used BaTiO3 -based ceramics. Since these finite-temperature thermophysical properties are characterized by the thermal andquantum fluctuation of nuclei, it is necessary toconsider the excitation of phonons for achieving ameaningful prediction.
First-principles phonon calculation is performed
actively and routinely based on density functional
theory (DFT). In most cases, the harmonic approximation (HA) is adapted, assuming that theatomic displacements are much smaller than interatomic distances. While the HA gives reasonable phonon dispersion curves and thermodynamic
quantities for many semiconductors and metals, it
often becomes inaccurate and sometimes breaks
down completely for high-performance thermoelectric and dielectric materials where the atomic
displacements can be as large as 10% of an interatomic distance. Moreover, the HA obviously
fails to describe physical properties related to theanharmonicity of the potential, including phononphonon scattering and structural phase transition.
To solve these limitations of the HA, a more advanced phonon calculation method that can include
phonon anharmonicity is required.
In the last few decades, several new ab initio phonon calculation methods have been proposed as an efficient way to include anharmonic
effects [1–7]. Among them, the self-consistent15

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyophonon (SCP) based approaches [3–7] are particularly important because they enable us to access
finite-temperature phonon frequencies and polarization vectors as well as Tc values of structuralphase transition [8, 9]. Also, the recent advancesin thermal transport theory are noteworthy [10,11].
The advent of these new methods, as well as theirefficient implementations in open-source software,
have opened up a way to study complicated lattice dynamics in emergent materials and to predict
finite-temperature thermophysical properties thatcannot be reached by the conventional method.
In this short report, we introduce the basics of the
new phonon calculation methods and demonstrate
their validity through applications to a hydrogenrich superconductor, thermoelectric clathrate, andthermoelectric tetrahedrite.

2

Recent progress of ab initio
phonon calculation

Self-consistent phonon theory
The self-consistent phonon (SCP) theory, which
was originally developed by Hooton [12], aims tocalculate phonon frequencies and polarization vectors that are renormalized by anharmonicity of thepotential. In this theory, an existence of a welldefined effective one-body Hamiltonian


X
1
†
H0 =
~Ωqν bqν bqν +(1)2qν
is assumed. Here, q and ν are the momentum
and branch index of phonons, respectively, andb†qν (bqν ) is creation (annihilation) operator of thephonon qν. The frequency Ωqν is the effective
frequency that is renormalized by anharmonic effects; therefore, its value differs from the harmonicfrequency ωqν . In addition, since the effective polarization vectors {qν } are also different from theharmonic ones {eqν }, b† and b in Eq. (1) differ
from the creation and annihilation operators of theharmonic Hamiltonian H0 .
Now, let us consider how to determine Ωqν and{qν }. To this end, we introduce the exact Hamiltonian H and associated density operator ρH =e−βH /Tr(e−βH ), with β = 1/kT . The Hamiltonian H includes all anharmonic terms; therefore, it
may be written as H = H0 + U3 + U4 + . . . where
Un is the potential energy of the nth-order anharmonicity. By using H and ρH , the exact Helmholtz
free-energy is given as
FH = Tr(ρH H) +1
Tr(ρH ln ρH ).
β(2)Next, let us substitute a trial density matrix defined as ρH0 = e−βH0 /Tr(e−βH0 ) for ρH in
Eq. (2). We then obtain FH (H0 ) = Tr(ρH0 H) +1
β Tr(ρH0 ln ρH0 ). By using the Jensen’s inequality ehXi ≤ heX i, it is straightforward to show thefollowing Feynman–Gibbs–Bogoliubov inequality
holds:
FH ≤ FH (H0 ).
(3)
In the SCP theory, Ωqν and {qν } are determinedso that the right-hand side of Eq. (3) is minimized.
This is nothing but a mean-field approximation.
Indeed, the SCP theory is a phonon version of the
Hartree–Fock theory.
The minimization of FH (H0 ) can be performedby repeatedly calculating atomic forces in supercells with stochastically sampled atomic configurations and updating the variational parameters in
H0 [5]. An alternative approach we have developed
is to derive the SCP equation via the condition of∂FH (H0 )/∂Xi = 0, where Xi comprises all variational parameters including Ωqν and {qν }. For
the brevity of the explanation, let us assume that thepolarization vectors do not change by anharmonic
effects. Also, we consider the anharmonic terms
up to the fourth-order because fifth- and higherorder terms are less significant. We then obtain thefollowing SCP equation [6, 13]:2
Ω2qν = ωqν+1X
Φ(qν; −qν; q1 ν1 ; −q1 ν1 )
2q ν
1 1

×

~[1 + 2n(Ωq1 ν1 )]. (4)2Ωq1 ν1

The second term on the right-hand side is the renormalization term associated with the quartic anharmonicity, whose magnitude is proportional to themean-square displacement of normal coordinate

16

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of TokyohQ∗qν Qqν i = ~[1 + 2n(Ωqν )]/2Ωqν with n(ω) being the Bose–Einstein distribution function. Sincen(ω) is temperature dependent, the effective frequency Ωqν also changes as a function of temperature.
We have shown that Eq. (4) gives phonon dispersion curves of cubic SrTiO3 that agree wellwith the results of inelastic neutron scattering
measurements. Once we compute all quartic
coefficients Φ(qν; −qν; q1 ν1 ; −q1 ν1 ) in Eq. (4),the effective frequencies at various temperatures
can be obtained quickly. To compute the quartic coefficients, we first calculate the real-space
fourth-order interatomic force constants (IFCs)
in a supercell. Then, the real-space IFCs aretransformed into the normal coordinate representation Φ(qν; −qν; q1 ν1 ; −q1 ν1 ). The computational complexity of this transformation isO(Nqirred. Nq1 Nν2 ), where Nqirred. is the numberof irreducible sets of q points, Nq1 is the number of q1 points, and Nν is the number of phonon
branches. When we consider the anharmonic mixing of polarization vectors, we also need to compute the off-diagonal elements of the quartic coefficients. In that case, the computational complexity becomes O(Nqirred. Nq1 Nν4 ). In our implementation of Eq. (4) in alamode [14], we use ahybrid parallelization of MPI and OpenMP. With
this implementation, the SCP calculation has beensuccessfully performed even for complex structures. For example, a SCP calculation of thermoelectric clathrate Ba8 Ga16 Ge30 (Nν = 162,
Nqirred. Nq1 = 1) with the polarization mixing finishes within 10 hours when we use 432 CPU coresof ISSP system B (sekirei).

Phonon scattering
While the SCP method can compute anharmonicphonon frequencies and polarization vectors efficiently, it is still based on a non-interacting picture [Eq. (1)]. Therefore, the phonon linewidth iszero, which leads to the unphysical result that thephonon lifetime and lattice thermal conductivity
are infinite. To solve this problem, phonon scattering processes must be considered. In bulk semi-

conductors, the most dominant source of phonon
scatterings is the phonon-phonon scattering associated with the lattice anharmonicity. In the lowtemperature region, other scattering processes, including phonon-impurity scattering and phononboundary scattering, also become relevant. However, ab initio treatment of these extrinsic phonon
scatterings is still challenging. Therefore, we will
focus on the intrinsic phonons scattering below.
The intrinsic scattering rate of phonons can beobtained by applying many-body perturbation theory to the anharmonic Hamiltonian. Here, thelowest-order perturbation term is the cubic anharmonicity U3 . The unperturbed Hamiltonian may
be the true harmonic Hamiltonian H0 , only if all
2 > 0), or thephonons are dynamically stable (ωqν
effective one-body Hamiltonian H0 obtained by an
SCP calculation. When we use H0 as the unperturbed Hamiltonian, the phonon self-energy ofthe bubble, which gives a first-order correction, isgiven as [15](B)Σqν (iωm ) =1 X X ~|Φ(−qν; q1 ν1 ; q2 ν2 )|2
16 q q ν νωqν ω1 ω2
1 2

1 2

× ∆(−q + q1 + q2 )f (1, 2, iωm ),(5)
X
1 + n1 + n2
f (1, 2, iωm ) =σiωm + σ(ω1 + ω2 )σ=±1
n1 − n2−, (6)iωm + σ(ω1 − ω2 )where ωm = 2πm/β~. Here, for the brevity of thenotation, we represent ωq1 ν1 as ω1 and n(ω1 ) as n1 .
The function ∆(q) in Eq. (5) becomes 1 when q is
an integral multiple of the reciprocal lattice vector
G and 0 otherwise. Therefore, the double loop
over q points in Eq. (5) can be reduced to a single
loop. After performing an analytic continuation to
the real axis, we obtain the phonon linewidth as(B)(B)Γqν = ImΣqν (ωqν ) and the phonon frequency(B)(B)
shift as ∆qν = −ReΣqν (ωqν ).
In many semiconductors, the linewidth calcu(B)lated from Γqν shows reasonable agreement withexperimental data, including its linear temperaturedependence in the high-temperature range. However, in some cases, the theoretical value may un-17

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

derestimate the experimental linewidth most likely
due to the presence of additional scattering channels of phonons, including phonon-impurity scattering and higher-order phonon-phonon scatterings. Recently, the self-energies of these additional
scattering processes have been evaluated based on
DFT for relatively simple systems and reported to
be important in correctly understanding the thermal transport behavior [16, 17], which is the main
topic of the next subsection.

Thermal transport
First-principles calculation of κL has been an active research topic in recent years because its active
control would lead to efficient energy conversion
and thermal management devices. In particular,significant effort has been placed on elucidating
the atomistic-level origin of ultralow thermal conductivity of efficient thermoelectric materials suchas intermetallic clathrates and SnSe. To computeκL from first principles, the semi-classical Peierls–
Boltzmann theory (PBT) is usually adopted. In the
PBT, we first assume that the heat flux is carriedby phonon quasiparticle; therefore, the heat flux
operator is approximated asXjph =~ωqν vqν b†qν bqν .
(7)qνHere, vqν = ∂ωqν /∂q is the group velocityof phonon, which is readily obtained from thephonon dispersion curves. Then, after introducing
the single-mode relaxation-time approximation, we
obtain the thermal conductivity tensorXκP =cqν vqν ⊗ vqν τqν ,(8)
qν

with cqν being the mode specific heat. The phonon
lifetime τqν can be obtained from the phonon
linewidth as τqν = ~/(2Γqν ). When we consider the dominant bubble self-energy for Γqν , thethermal conductivity shows the temperature dependence of κL ∝ T −1 in the high-temperature
range, in accord with many experimental observation. So far, the PBT-based calculations have been
performed for many thermoelectric materials andexplained their thermal transport behavior successfully.
While the PBT has been successful for many
crystalline solids, it fails to describe the temperature profile of κL observed in amorphous solids
and other disordered solids. In these systems, thethermal conductivity shows little temperature dependence and does not decrease even in the hightemperature region. To understand thermal transport physics in disordered solids, a beyond PBTtreatment would be necessary. One of the mostproblematic assumptions made within the PBT formalism is the approximation of Eq. (7). Morespecifically, when we express the harmonic heat
flux operator j (2) with the creation and annihilation
operators of phonons, we obtain the band diagonal
and off-diagonal terms as follows:j (2) =X ~(ωqν + ωqν 0 )qνν 0= jph +

2vqνν 0 b†qν bqν 0

X ~(ωqν + ωqν 0 )qν6=ν 02vqνν 0 b†qν bqν 0 .
(9)Here, vqνν 0 is the band off-diagonal extension of
the group velocity. In the PBT, the band offdiagonal terms of the heat flux operator is omitted. In the pioneering work by Allen and Feldman [18, 19], however, they showed that the offdiagonal contribution is dominant in the disorderedsolids and successfully explained the temperaturedependence of κL observed in amorphous silicon.
In the Allen–Feldman (AF) theory, however, the effect of phonon scattering is not considered. Therefore, their theory is still not suited for studying thermal transport in disordered solids that also show
strong anharmonicity.
Recently, a new thermal transport theory thatunifies the PBT and the AF theory has been developed [10]. According to their result, the thermal
conductivity is given asκL = κP + κC ,(10)

where κP is the Peierls term defined by Eq. (8), andκC is the coherent term associated with the interband components of the heat flux operator, which

18

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

So far, we have introduced the state-of-the-art
ab initio phonon calculation methods. To study
phonon-related properties by using these methods,
the harmonic, cubic, and quartic IFCs are necessary
as inputs. These IFCs can be estimated by repeatedly calculating atomic forces in a supercell withsuitably chosen atomic configurations. In the conventional supercell approach, only a few atoms in asupercell are displaced. For instance, the harmonic
IFC between atom i and j, Φij , is estimated as
Φij = −[Fj (ui = +h)−Fj (ui = −h)]/2h, where
Fj (ui = ±h) represents the atomic force acting on
the atom j when the atom i is displaced slightly(h ≈ 0.01 Å) from its equilibrium position. This
procedure needs to be repeated for all irreducible
ij pairs. To estimate the cubic IFC, Φijk , in thesame fashion, we need to consider, at most, four
displacement patterns for each element. Therefore,
the number of displacement patterns necessary toobtain all relevant IFCs increases rapidly for thehigher-order terms, and the calculation of quartic
IFCs would be infeasible or extremely expensive if
possible.
To mitigate this issue, we employ the compressive sensing approach [20, 21]. First, to increase
the information density obtained from a single DFT
calculation, we displace all atoms in a supercell in
random directions. For these sampled atomic configurations, we compute the atomic forces; these
displacement-force data sets will be served as training data sets of a linear model. Second, instead of105

conventional
LASSO

104

103

102

8G
a1
Cu 6 Ge
3
12
Sb 0
4S
1
Sn 3
Se
cSr
Ti
O
t-S 3
rT
iO
c3
Cs
Pb
B
t-C
r
sP 3
bB
oCs r3PbBr
3

Efficient calculation of force constants

Ba

While this new theory has been tested only for a
few systems so far, it is expected to improve theprediction accuracy of κL and help us deepen our
understanding of anomalous thermal transport instrongly anharmonic and/or disordered materials.

3Γqν + Γqν 0. (11)(ωqν − ωqν 0 )2 + (Γqν + Γqν 0 )2Tevqνν 0 ⊗ vqν 0 νF×ωqν + ωqν 0

Pb

q ν6=ν 0Si

X X cqν ωqν 0 + cqν 0 ωqν

Sc

κC =

calculating the IFC one by one by the central difference, we estimate all IFCs simultaneously by performing linear regression. Since the atomic forces
F can be expressed as a linear function of IFCs Φ
as F = AΦ with A being a sensing matrix comprising atomic displacements, we can estimate theparameter vector Φ by using ordinary least-squares
or other regression methods. In particular, the L1
penalized regression models such as LASSO areuseful because they select and compute physicallyimportant IFCs somewhat automatically and set theother irrelevant parameters exactly zero.
We have applied LASSO to various systems andtested its performance. As shown in Fig. 1, thenumber of static DFT calculations could be reduced dramatically from what would be required
if the conventional approach were used. In these
LASSO calculations, we gradually increased thetraining data sets and the cutoff radii of interaction until we obtained an accurate model that reproduces the DFT forces within ∼5% error. The
numbers shown for the conventional method are estimated by symmetry argument, where anharmonic
IFCs up to the fourth-order are considered within
the cutoff radii as the LASSO calculation. The efficient approach, which can be 100 times as efficientas the conventional one, enables us to perform the
SCP and PBT calculations of complex systems.

No. of static DFT calcs.

is given as

Figure 1: Number of displacement patterns thatwas used to extract IFCs by LASSO (cyan) compared with the estimated number of patterns required by the conventional approach (red).

19

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo
Harmonic

SCP (T = 0 K)

300

Frequency (meV)

250
200
150
100
50
0
−50
Γ

X U|K

Γ

L

W

X

Figure 2: Phonon dispersion curves of F m3̄mLaH10 at 150 GPa calculated within the HA andthe SCP theory at T = 0 K.

3Applications

Quantum fluctuation in record superconductor LaH10
Recently, record superconductivity at ∼250 K has
been reported for LaH10 under the high pressure of
137–218 GPa [22, 23]. The possibility of the highTc superconductivity in LaH10 was first predicted
by DFT calculations [24,25] and subsequently confirmed experimentally. In the previous DFT predictions, a sodalite-like structure of LaH10 withspace group F m3̄m has been suggested as a crystal structure that realizes a high-Tc . However, if
the phase stability is evaluated based on enthalpy,
LaH10 is predicted to be thermodynamically unstable against decomposition below 200 GPa. Moreover, F m3̄m-LaH10 is predicted to be dynamically unstable over the whole pressure range wherea ∼250 K Tc has been observed experimentally.
These contradictions likely indicate the limitationof the conventional computational approaches.
To solve this puzzle, we have performed phonon
calculations of LaH10 with including anharmonic
effects [26]. Figure 2 shows the phonon dispersion of F m3̄m-LaH10 calculated within the HA
and the SCP theory at 150 GPa. The F m3̄mstructure is dynamically unstable when the HA isemployed, indicating the F m3̄m structure distorts
into a low-symmetry structure. Indeed, we haveconfirmed that the enthalpies of R3̄m- and C2LaH10 are lower than that of F m3̄m-LaH10 by∼50 meV per formula unit (f.u.) at 150 GPa. However, this conclusion changes when we consider thezero-point motion of hydrogen atoms. As shownin Fig. 2, the quartic anharmonicity renormalizes
phonon frequencies significantly and stabilizes thephonons in the entire Brillouin zone even at T = 0
K, which is possible due to the large zero-point
motion of hydrogen atoms. After obtaining the stable phonons, we have also estimated the zero-point
energy (ZPE) and evaluated the relative stabilityof the F m3̄m structure over the distorted phases
based on enthalpy+ZPE. As a result, the F m3̄mphase is found to be stable than the R3̄m structure by ∼ 110 meV / f.u., which agrees better withthe experimental observation. The inclusion of the
ZPE also changes the thermodynamic stability of
F m3̄m-LaH10 . For example, the formation enthalpy defined as∆H[LaH10 ] = H[LaH10 ]71− H[LaH11 ] − H[LaH3 ] (12)88is about 140 meV/f.u. at 200 GPa, indicating that
LaH10 decomposes into LaH3 (Cmcm) and LaH11(P 4/nmm). However, this conclusion changes
when we consider the ZPE additionally; the formation enthalpy becomes negative as ∆H +∆ZPE =−280 meV/f.u.
To summarize, our calculations have clearlyshown that the large quantum fluctuation of hydrogen nuclei plays a crucial role in stabilizing
the highly-symmetric F m3̄m-LaH10 over the wide
pressure range of 137–218 GPa. We have also
calculated the superconducting Tc based on the
Migdal–Eliashberg theory and obtained an excellent agreement with the experimental values obtained for LaH10 and LaD10 [26]. This agreement
further confirms that the F m3̄m phase of LaH10 isindeed responsible for the record superconductivity
at ∼250 K.

20

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

(a)

(b)

Harmonic(c)

SCP (300 K)

Lattice thermal conductivity (W/mK)

35

Frequency (meV)

30
25
20
15
10
5
0

°

X

M

°

R

3.0
∑P (based on HA)
2.5∑P (based on SCP)

2.0

∑C (based on SCP)
∑P + ∑C
Expt. (Sales et al.)Expt. (May et al.)

1.5
1.0
0.5

200

400

600

800

1000

(d)

Harmonic(e) 50SCP (300 K)(f)

Frequency (meV)

40
30
20
10
0
°10

°

H

P

N

°

Lattice thermal conductivity (W/mK)

Temperature (K)
1.0
∑P
∑C

∑P + ∑C
Expt. (Suekuni et al.)

0.8

0.6

0.4

0.2

0.0

200

400

600

800

1000

Temperature (K)Figure 3: Crystal structures and phonon properties of thermoelectric clathrate Ba8 Ga16 Ge30 (BGG) andtetrahedrite Cu12 Sb4 S13 . (a) Crystal structures of (a) BGG, adopted from Ref. [27], and (d) tetrahedrite.
(b,e) Phonon dispersion curves of BGG and tetrahedrite, respectively, calculated within the HA (dottedlines) and the SCP method at 300 K (solid lines). (c,f) Thermal conductivities of BGG and tetrahedrite,
respectively, calculated by different methods. For comparison, the experimental values of Refs. [28–30]
are also shown.

Thermal conductivity in complex thermoelectric materials
The thermal transport in complex thermoelectricmaterials has attracted special attention not only
because their thermal conductivity is unusually
low but because they sometimes show a glasslike
temperature dependence even though they possess periodic crystal structures. In this report, we
will focus on the origin of the unusual temperature dependence of κL observed in type-I clathrate
Ba8 Ga16 Ge30 (BGG) and tetrahedrite Cu12 Sb4 S13 ,whose crystal structures are shown in Figs. 3 (a)and (d). In our previous study [31], we have shownthat the low-κL value of BGG can be attributed
to the strong phonon-phonon scattering inducedby the presence of low-frequency optical phonons.
We could successfully reproduce the experimental κL values up to ∼ 100 K by using Eq. (8).
However, the temperature dependence weaker thanκL ∝ T −1 could not be explained by the conventional PBT. In Cu12 Sb4 S13 , κL is almost tempera-ture independent, which cannot be explained by theconventional PBT. Actually, the HA yields imaginary phonons for Cu12 Sb4 S13 , which hinders a
direct application of Eq. (8).
To understand the origin of the unusual temperature dependence in BGG and tetrahedrite, wehave recently performed SCP calculations of thesecomplex materials. Also, the thermal conductivity
is evaluated based on Eq. (10) for which the SCP
eigenvalues and eigenvectors are used as inputs.
Figures 3 (b) and (e) show the calculated phonon
dispersion curves of BGG and tetrahedrite, respectively. The phonon dispersion curves calculated
within the HA and SCP theory are somewhat similar, particularly for the high-frequency phononmodes. In the low-frequency range, however, we
can see notable differences. In BGG, the frequencyof low-lying optical phonon modes around 4.5 meV
increases with heating. These phonon modes correspond to the collective vibration of Ba guest atoms
inside the large cavities made by the host frame-21

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

work, as indicated by the red arrows of Fig. 3 (a).
These phonon modes are often called “rattling“ motion. In tetrahedrite, the unstable phonon modes exist within the HA. These phonon modes correspondto the collective vibrations of Cu(2) atoms, which
are located on the triangle planes made by sulfur
atoms (see Fig. 3(d)). Since the displacement of theCu(2) atoms occur in the direction perpendicular
to the triangle planes, these characteristic phononmodes were named “planar rattling” [30]. In our
SCP calculations, the planar rattling modes werestabilized by anharmonic renormalization above ∼80 K.
Next, we have calculated κP and κC of Eq. (10)separately; the results are compared with available
experimental data in Figs. 3 (c) and (f). When we
employed the harmonic phonon dispersion in calculating the phonon lifetimes and κP , the thermal
conductivity was significantly underestimated, as
shown by the dashed line of Fig. 3 (c). When the
SCP dispersion was used instead, we could reacha nice agreement with the experimental data, including the temperature dependence weaker thanκL ∝ T −1 . Since the low-lying optical phonons
strongly hybridize with the heat-carrying acoustic phonons, a little temperature-dependence ofthe rattling-mode frequency gives rise to the unusual temperature profile of κL [27]. The coherent
contribution κC calculated from Eq. (11) accounts
for about 10% of total κL of BGG at 300 K. For
tetrahedrite, the temperature dependence of κL is
more exceptional as shown in Fig. 3 (f). Interestingly, the little temperature dependence can be
explained qualitatively by the Peierls term basedon the temperature-dependent SCP wave functions.
Still, the calculated κP value significantly underestimates the experimental thermal conductivity
for the entire temperature range. This discrepancy can be cured by including the coherent termκC ; the sum of two contributions, κP + κC , reproduces the experimental values even quantitatively. Such a large coherent contribution has already been reported for orthorhombic halide perovskite CsPbBr3 [10], where the particle-like thermal transport is hindered by the strong phonon-

phonon interaction. We believe the coherent term
κC should be significant not only in disorderedsolids but also in strongly anharmonic (ordered)solids where the phonon linewidths become broad,and the phonon spectral weights overlap with eachother.

4

Summary

In this short report, we reviewed the recent development of new phonon calculation methods that cantreat the intrinsic effects of lattice anharmonicity
beyond the quasiharmonic level. We showed thatthe self-consistent phonon scheme solves the limitations of the harmonic phonon theory by incorporating the quartic anharmonicity at the mean-fieldlevel. Also, the unified theory of thermal transport,
which can account for the coherent thermal transport in addition to the particle-like contribution,
was shown to better explain the anomalous thermal transport observed in complex thermoelectricmaterials.
All of the methods focused in this report are already available in the alamode software [32]. We
kindly invite interested readers to visit the homepage and try the software.
This work was partially supported by Grant-inAid for Scientific Research (number 16H06345)
from the Ministry of Education, Culture, Sports,
Science and Technology, Japan. A part of the computation in this work was performed using the facilities of the Supercomputer Center, Institute for
Solid State Physics, The University of Tokyo. This
report is based on the recent studies done in collaboration with K. Suekuni, S. Tsuneyuki, J. A.
Flores-Livas, I. Errea, R. Arita, and T. Koretsune.

References[1] O. Hellman and I. A. Abrikosov: Phys. Rev.
B 88 (2013) 144301.
[2] T. Sun, D.-B. Zhang, and R. M. Wentzcovitch:
Phys. Rev. B 89 (2014) 094109.

22

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[3] P. Souvatzis, O. Eriksson, M. I. Katsnelson,
and S. P. Rudin: Comp. Mater. Sci. 44 (2009)888.
[4] B. Monserrat, N. D. Drummond, and R. J.
Needs: Phys. Rev. B 87 (2013) 144302.
[5] I. Errea, M. Calandra, and F. Mauri: Phys.
Rev. B 89 (2014) 064302.
[6] T. Tadano and S. Tsuneyuki: Phys. Rev. B 92
(2015) 054301.
[7] A. van Roekeghem, J. Carrete, and N. Mingo:
Phys. Rev. B 94 (2016) 020303.
[8] R. Bianco, I. Errea, L. Paulatto, M. Calandra,
and F. Mauri: Physical Review B 96 (2017)014111.
[9] T. Tadano and S. Tsuneyuki: Journal of the
Ceramic Society of Japan 127 (2019) 404 .
[10] M. Simoncelli, N. Marzari, and F. Mauri: Nature Physics 15 (2019) 809.
[11] L. Isaeva, G. Barbalinardo, D. Donadio, and
S. Baroni: Nature Communications 10 (2019)3853.
[12] D. J. Hooton: Philos. Mag. 3 (1958) 49.
[13] Y. Oba, T. Tadano, R. Akashi, and
S. Tsuneyuki: Physical Review Materials 3(2019) 033601.
[14] T. Tadano, Y. Gohda, and S. Tsuneyuki: J.
Phys: Condens. Matter 26 (2014) 225402.

[20] F. Zhou, W. Nielson, Y. Xia, and V. Ozolin, š:
Phys. Rev. Lett. 113 (2014) 185501.
[21] T. Tadano and S. Tsuneyuki: Journal of the
Physical Society of Japan 87 (2018) 041015 .
[22] M. Somayazulu, M. Ahart, A. K. Mishra,
Z. M. Geballe, M. Baldini, Y. Meng, V. V.
Struzhkin, and R. J. Hemley: Physical Review Letters 122 (2019) 027001.
[23] A. P. Drozdov, P. P. Kong, V. S. Minkov, S. P.
Besedin, M. A. Kuzovnikov, S. Mozaffari,
L. Balicas, F. F. Balakirev, D. E. Graf, V. B.
Prakapenka, E. Greenberg, D. A. Knyazev,
M. Tkacz, and M. I. Eremets: Nature 569(2019) 528 .
[24] H. Liu, I. I. Naumov, R. Hoffmann,
N. Ashcroft, and R. J. Hemley: Proceedings of the National Academy of Sciences 114(2017) 6990.
[25] F. Peng, Y. Sun, C. J. Pickard, R. J. Needs,
Q. Wu, and Y. Ma: Physical review letters
119 (2017) 107001.
[26] I. Errea, F. Belli, L. Monacelli, A. Sanna,
T. Koretsune, T. Tadano, R. Bianco, M. Calandra, R. Arita, F. Mauri, and J. A. FloresLivas: Nature 578 (2020) 66.
[27] T. Tadano and S. Tsuneyuki: Physical review
letters 120 (2018) 105901.

[15] A. A. Maradudin and A. E. Fein: Phys. Rev.
128 (1962) 2589.

[28] B. C. Sales, B. C. Chakoumakos, R. Jin, J. R.
Thompson, and D. Mandrus: Phys. Rev. B 63
(2001) 245113.

[16] T. Feng, L. Lindsay, and X. Ruan: PhysicalReview B 96 (2017) 72 .

[29] A. F. May, E. S. Toberer, A. Saramat, and G. J.
Snyder: Phys. Rev. B 80 (2009) 125205.

[17] N. K. Ravichandran and D. Broido: PhysicalReview B 98 (2018) 1 .

[30] K. Suekuni, C.-H. Lee, H. I. Tanaka, E. Nishibori, A. Nakamura, H. Kasai, H. Mori,
H. Usui, M. Ochi, T. Hasegawa, M. Nakamura, S. Ohira-Kawamura, T. Kikuchi,
K. Kaneko, H. Nishiate, K. Hashikuni,
Y. Kosaka, K. Kuroki, and T. Takabatake:
Advanced Materials 285 (2018) 1706230 .

[18] P. B. Allen and J. L. Feldman: Physical Review Letters 62 (1989) 645.
[19] P. B. Allen and J. L. Feldman: Physical Review B 48 (1993) 12581 .

23

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[31] T. Tadano, Y. Gohda, and S. Tsuneyuki: Phys.
Rev. Lett. 114 (2015) 095501.
[32] T.
Tadano.
ALAMODE.
https://github.com/ttadano/alamode.

24