35th Computational Materials Design (CMD) Workshop

Phonon anharmonicity
from first principles
Terumasa TADANONational Institute for Materials Science (NIMS), Tsukuba, Japan
Sep. 6, 2019

Outline
1. Introduction to phonon anharmonicity
2. Perturbative treatment of anharmonicity
3. Non-perturbative treatment of anharmonicity
4. Efficient calculation of anharmonic force constants
5. Applications
6. Summary

2

Phonon (lattice vibration)• Density functional theory (DFT) is theory at absolute zero (T = 0 K)
• Lattice or molecular vibration affects physics at finite temperature.

How to treat lattice vibrationin crystals?

Ab initio phonon calculation
Ab initio molecular dynamics
3Phonon physics
Phonon dispersion relation• Raman, IR, INS, IXS measurements• Dynamical stability

Thermodynamics & mechanics• Vibrational entropy, Helmholtz free energy, …
• Lattice constant, elastic constants at finite temperature• Phase transformation at finite temperature
Transport• Phonon thermal transport• Scattering of electrons by el-ph coupling• Phonon-mediated superconductivity
Dielectrics• Dielectric constant• Band-gap renormalization by phonons
4

Lattice vibration in energy materials
Thermalconductivity

Thermal stability
Structuralphase transition

MOF

Photovoltaic
materials
Thermoelectricmaterials

Electron-phonon
coupling(Piezo, ferro)electricmaterials• Quantitative prediction & robust understanding of thermophysicalproperties is desired.

• Challenging due to anharmonic effects
5

Thermoelectric (TE) materials• Thermoelectric (TE) devices efficiently
convert a temperature difference to anelectrical potential, and vice versa.
• Electrical power generation (Seebeck
effect)• Cooling (Peltier effect)
J.-F. Li et al., NPG Asia Mater 2, 152 (2010).

• About 60% of energy is lost as heatNASA Voyager

BMW 530i

6

Thermoelectric figure of merit
For high energy efficiency, new and abundant materials that possess highfigure-of-merit are necessary
Figure of merit

Electrical conductivity> 103 / (Ω･m)Seebeck coefficient> 100 μV/K

Lattice Thermal Conductivity(LTC) < 1.0 W/mK

ZT > 1 is necessary, ZT > 2 or higher is desired.
Design principle looks simple but realizing higher ZT is challenging.
7

Strong anharmonicity in thermoelectric materials• Chalcolagenide TE materials

Bi2Te3
AgGaTe2• PbTe, Bi2Te3, SnSe, AgGaTe2, Cu2Se, …
• Lone pair

• Host-guest structures
SnSe• Clathrates, Skutterudites, Pyrochrore• “Rattling” guest

Large atomic displacement factor
Strong anharmonicity
Ultralow κL & its unusual T-dependence

Large atomic
displacement

Stronganharmonicity

Clathrate

Skutterudite

Breakdown of
harmonic approx.
8

Temperature Dependence of the Energy Levels of Methylammonium
Lead Iodide Perovskite from First-Principles

Strong anharmonicity in other materials
Wissam A. Saidi,*,† Samuel Poncé,‡ and Bartomeu Monserrat¶,§†

Department of Mechanical Engineering and Materials Science, University of Pittsburgh, Pittsburgh, Pennsylvania, 15261, UnitedStates‡Department of Materials, University of Oxford, Parks Road, Oxford OX1 3PH, United Kingdom¶
Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854-8019, United States§
Purposeof our study
TCM Group, Cavendish Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE, 1.2
UnitedKingdom

High-Tc superconductors*
H3S
S

Supporting Information

Negative
thermal expansion
ScF3

Perovskite solar cell
MAPbI3

1.2 Purpose of our study
(a)

(b)

ABSTRACT: Environmental eﬀects and intrinsic energy-loss processes lead to(b)ﬂuctuations in the operational temperature of solar cells, which can profoundly(a)inﬂuence their power conversion eﬃciency. Here we determine from ﬁrstprinciples the eﬀects of temperature on the band gap and band edges of thehybrid pervoskite CH3NH3PbI3 by accounting for electron−phonon coupling
and thermal expansion. From 290 to 380 K, the computed band gap change of
40 meV coincides with the experimental change of 30−40 meV. The calculationof electron−phonon coupling in CH3NH3PbI3 is particularly intricate as thecommonly used Allen−Heine−Cardona theory overestimates the band gapchange with temperature, and excellent agreement with experiment is only
obtained when including high-order terms in the electron−phonon interaction.
We also ﬁnd that spin−orbit coupling enhances the electron−phonon couplingstrength but that the inclusion of nonlocal correlations using hybrid functionals
Figure1.3: (a):Figure1.3:
Structure(a):of Structure
metal trifluorides.
This structure
is similar
to thetoperovskite
but butthere
has little eﬀect. We reach similar conclusions in the metal−halide perovskiteof metal trifluorides.
This structure
is similar
the perovskiteth
CsPbI3. Our results unambiguously conﬁrm for the ﬁrst time the importance of high-order terms in the iselectron−phononno A-site atoms.
Fluorides
the octahedra
around metals.
(b): Lattice
parameter
andandtheis no A-site
atoms.form
Fluorides
form the octahedra
around metals.
(b): Lattice
parameter
coupling by direct comparison with experiment.
CTETheof ScF
Theisblack
is fitted
to thelattice parameterred dashedline
is linear
the lin
CTE of ScF3 [5].
black
line
fittedline
to thelatticeparameter
and theand
redthedashedline isthe
3 [5].

J. Phys. Chem. Lett. 7, 5247 (2016).

Nature 525, 73 (2015).

CTE.
green
line is thebetween
diﬀerencethe
betweenthe experimentally
determinedparameter
CTE. The green
lineTheis thediﬀerenceexperimentally
determinedlatticelattice
parameter
anda
therefore the performance of the solar cell, as was demonstrated
thesplinefit.
the spline fit.
for MAPbI3-based solar cells, which showed an open-circuit
voltage drop from 1.01 to 0.83 V as the temperature increased
from 300 to 360 K.14 Recently, two experimental studies haveinvestigated the temperature-dependent properties of MAPbI3
and reported changes in the band gap reaching 40 meV in
the Purpose of our study
1.2
1.2and Purpose
even larger of our study
temperature range of solar cell operation15,16
changes exceeding 100 meV in the individual band edges.15From the previous studies, it is known that anharmonicity is important for understanding NTE
While the properties of MAPbI3 are well-characterizedboth studies, it is known that anharmonicity is important for understanding NTE in
From the previous
ScF3 . However,
how each-order anharmonicity plays a role in NTE is unclear and the microsco
theoretically and experimentally, there is an increasing researchwiththeeﬀort to ﬁnd alternatives to this light absorber,
ScF3 . 17−20
However,how each-orderanharmonicity
a role in NTE
is still
unclearand Hence,
the microscopic
comprehension
neededfor designing plays
and controlling
NTE islacking.
in this stu
aim of tuning the physical properties of the solar cells andwe aim to forunderstandthe andrelation
between NTE
anharmonicity
more deeply
by study,
using
comprehension
designing
controlling
NTE andis stilllacking. Hence,in thismoving toward lead-free devices. In this context,
it wouldneededbeeﬃcient first-principles
where
we andanalyzethe eﬀect of cubic
quarticby
anharmonicit
desirable to accurately calculate the eﬀects
of temperature
we aimto understandthe relationmethod,
betweenNTE
anharmonicity
moreanddeeply
using anchanges to aid in the computational design of novel perovskite
separately. This calculation method is independent of the structure of the target materials aeﬃcient first-principles
method, where we analyze the eﬀect of cubic and quartic anharmonicities
solar cell materials under real working conditions.
Unfortufor the discussion of phase stability.
with
nately, the organic−inorganic nature of MAPbI3, togetherapplicablecalculation method is independent of the structure of the target materials andthe presence of heavy Pb atoms, makes it aseparately.
challengingThis
system21−28
to study from ﬁrst-principles.
applicable for the discussion of phase stability.

Metal−halide hybrid perovskites are promising candidates
as light absorbers for the next generation of solid-statehybrid photovoltaic devices.1−7 Their solar power conversion
eﬃciency has reached 22.1%,8 an impressive achievement
considering that the ﬁrst perovskite-based solar cell was only
made 7 years ago.1 These light-absorbing materials belong to
the perovskite family ABX3, where A is a cation located at thecenter of a cube formed from corner-sharing BX6 octahedra.
Methylammonium lead iodide, MAPbI 3 , where MA
represents the organic cation CH3NH+3 , has an optimal band
gap for light absorption of 1.6 eV and exhibits large electron
and hole diﬀusion lengths, which make it a prime candidate forsolar cell applications.3,4 Another factor determining theeﬃciency of a solar cell is the band alignment between thelight absorber and its transport layers.9 As an example, acommon electron transport layer is TiO2, which has aconduction band that is only 0.1 eV lower than that of
MAPbI3.10
Solar cell devices are subject to signiﬁcant temperature
variations under normal operating conditions.11 Ambient
temperatures can oscillate between 250 and 310 K, andenergy-loss mechanisms associated with the conversion of solar
power to electric power can increase the operating temperatureof the cell above 340 K.11−13 These temperature ﬂuctuations
can aﬀect both the absorption onset and band alignment andStructuralphase transition

SrTiO3Anharmonicity affects many
properties

105 K

Tetragonal

Cubic© 2016 American Chemical Society

JACS 132, 15496 (2010).

Thermodynamical stability
Optical properties
1.3 Outline of the thesis
Transport
properties
1.3 Outlineof the thesis
Elastic properties

Received: November 2, 2016
Accepted: December 3, 2016
Published: December 4, 2016
5247

In Chapter 2, we introduce some conventional methods to calculate CTE and the one we used

this study. In Chapter 3 and 4, we show the results of the electronic structure, phonon frequenc

DOI: 10.1021/acs.jpclett.6b02560
and CTEs of ScF3 and detailed discussion on these results. The conclusion of my thesis is giv
J. Phys. Chem. Lett. 2016, 7, 5247−5252

In Chapter 2, we introduce some conventional methods to calculate CTE and the one we used in

in Chapter3 5.
this study. In Chapter
and 4, we show the results of the electronic structure, phonon frequencies

and CTEs of ScF3 and detailed discussion on these results. The conclusion of my thesis is given
in Chapter 5.

9How to include anharmonic effects?
Many-body theory···✓ Quantum statistics✓ Thermodynamic limit is accessible···✗ Anharmonic effects are included only partiallyanharmonic correction✗ Requires anharmonic coupling constants

Molecular Dynamics (MD) methods✓ Full anharmonicity

MI R̈I = FI<latexit sha1_base64="eCctYDp+jL7gpv8EB0PtwKycHO0=">AAAC0HichVE9S+xQED3Gz+fn6mse2CwuPqyWiQiKIIiCaCHo6qrgypJkr2swXyR3FzUs8kr9AxZWChZiYau1jWBt4U8QSwUbizfJBh5PUSckd+bMnMmZO7pnmYEkemhQGpuaW1rbfrR3dHZ196R6+1YCt+IbIm+4luuv6VogLNMReWlKS6x5vtBs3RKr+vZ0lF+tCj8wXWdZ7npiw9bKjrlpGppkqJjKzBfDuVqhVHJlWNDtMFerRUB6Ih1FM3HAVZSl2NIfHTVxMkhswU3doYASXBiowIaAA8m+BQ0BP+tQQfAY20DImM+eGecFamhnboWrBFdojG7zt8zReoI6HEc9g5ht8F8sfn1mpjFI93ROz3RLF/RIb5/2CuMekZZdPvU6V3jFnsNfS6/fsmw+Jbb+sb7ULLGJsVirydq9GImmMOr86t7R89J4bjD8Taf0xPpP6IFueAKn+mKcLYrc8Rd6dJ5b41hwz89vLqqzozwvUn2/to/OynBWpay6OJKZnEpW2oZ+DGCI9zaKScxiAXnuf4BLXOFaySk7yr7yp16qNCScn/jPlMO/FFynTw==</latexit>✓ (Dynamical) disorders✗ Classical statistics✗ Long-time & large-scale simulation is necessary✗ DFT-MD is very costly

Y. He et al., Phys. Chem. Chem. Phys. 14, 16209 (2012).
10

Outline of our approach
Anharmonic
lattice model (ALM)

DFT calculation

Many-body theory······

compressive
sensing

Thermal properties
Molecular dynamics(DFT, Machinelearning FF)

•
•
•
•
•
•

Phonons at finite temperature
Vibrational free-energy (beyond QHA)
Structural phase transition
Thermal conductivity
Thermal expansion…
11

Outline
1. Introduction to phonon anharmonicity
2. Perturbative treatment of anharmonicity
3. Non-perturbative treatment of anharmonicity
4. Efficient calculation of anharmonic force constants
5. Applications
6. Summary

12

Phonons (harmonic)kLinear chain

dispersion relationapotential energy V =

X1
i

2

k(ui

ui 1 ) 2
displacement!q

Eq. of motionq

Generalization to 3D
dispersion relation

potential energy

Harmonic force constants

Eq. of motion
13

Taylor expansion of potential energy surface
U0 = U2 + U3 + U4 + · · ·1 X=µ1 µ2 (`1 1 ; `2 2 ) ⇥ uµ1 (`1 1 )uµ2 (`2 2 )2U{`,,µ}1 X+
3!

{`,,µ}1 X+
4!

{`,,µ}<latexit sha1_base64="iEPquvQ+RQxcDv2rE7KKQYvu8Zg=">AAAFW3icvVJNa9RQFL2dNFrH1k4VQXDzdOjQ0jq8SQKKUii4cTltTVvoKyGTvmnD5It8DNSQjUt178KVgoj4M9z4B1x04Q8Ql1XcuPDmTWaGTutnxReSc3PuPTfnvpdW4NhRTOnBWEkal8+cnThXPj85dWG6MnNxPfKT0OK65Tt+uNkyI+7YHtdjO3b4ZhBy0205fKPVuZvnN7o8jGzfux/vB3zbNXc9u21bZoyUMVP6opMbRDdSmpHaUh4oGVnIUS1Qy5FZO34cEcbKNbJEWDs0rbSRYS2LEtdIWcq44yyyjhkE5iJzE5ZlhDX3bEy5iYGlApQsm8sLBSFq8+hOj1IGlJLNExbbLo9IMmhwknJ+kB92HmmTG14YGlav/aljAervGi8odUCp/2yWfnroZeQrI6NqfzuqAO0UExeUNqC0/7YJ/fzQ/lEX/T/ZqFRpnYpFjgeNIqhCsZp+5RUw2AEfLEjABQ4exBg7YEKE1xY0gEKA3DakyIUY2SLPIYMyahOs4lhhItvB5y6+bRWsh+95z0ioLfyKg3eISgKz9D19TQ/pO/qGfqTfftgrFT1yL/uIrZ6WB8b0oytrX3+pchFj2Buqfuo5hjbcEl5t9B4IJp/C6um7D54ert1enU1r9AX9hP6f0wP6Fifwup+tlyt89RmU8QAao9t9PFhX6g21rqxo1WWtOIoJuArXYQ73+yYswz1ogg6WxKSH0mPpyfgHWZLL8mSvtDRWaC7BkSVf/g7E9aK+</latexit>µ1 µ2 µ3 (`1 1 ; `2 2 ; `3 3 ) ⇥ uµ1 (`1 1 )uµ2 (`2 2 )uµ3 (`3 3 )µ1 µ2 µ3 µ4 (`1 1 ; `2 2 ; `3 3 ; `4 4 ) ⇥ uµ1 (`1 1 )uµ2 (`2 2 )uµ3 (`3 3 )uµ4 (`4 4 ) + · · ·R0 (`) is atomic displacement from the equilibrium position R0 (`)u(`) = R(`)<latexit sha1_base64="JRsbITa7iZJ19o6uo/0VCFggPEc=">AAAC7XichVFNSxtRFD2OtvWjrVE3gptgsOii4SYNVAQh0I3LGBsVjA0z0xc75M0HMy8BHfIH3EsXLYJCK6U/w41gtxH8CeIyQjdd9M5koEQxvcPMO+/ce+47b67hSStQRNdD2vDIk6fPRsfGJ56/eDmZmpreDNymb4qK6UrX3zb0QEjLERVlKSm2PV/otiHFltF4F+W3WsIPLNd5r/Y9sWvre45Vt0xdMVVLFauGHTbbi1UhZbWhe56+lF5NR2S5n3ydkB9C6kvUUhnKUhzphyCXgAySKLmpS1TxES5MNGFDwIFiLKEj4GcHORA85nYRMuczsuK8QBvjrG1yleAKndkGf/d4t5OwDu+jnkGsNvkUya/PyjQWqEM/qEsX9JNu6M+jvcK4R+Rln1ejpxVebfJwduP3f1U2rwqf/qkGelaoYzn2arF3L2aiW5g9fevgc3djpbwQvqJTumX/J3RN53wDp3VnflsX5S8D/Bjs5fE/FuV9xg3GKqrjQebuj+0h2Mxnc2+y+fVCplhIRjqKOcxjkef2FkWsoYQKn/Mdv9DBleZqR9pX7bhXqg0lmhn0hXb2F73gscs=</latexit><latexit sha1_base64="3Yasg3XIfy1X8OM6b5eaNs3reTg=">AAACxHichVHLSsNQEB3jq9ZX1YWCm2JR6qZMakFxVRDEZVttFayWJF419OZBcluooX6AP+DClYoLUb/CjeDaRT9BXCq4ceEkDYiKdUJyzz0zZ3LuHdXmuisQW11Sd09vX39kIDo4NDwyGhsbL7lWzdFYUbO45Wypisu4brKi0AVnW7bDFEPlbFOtrvj5zTpzXN0yN0TDZjuGcmDq+7qmCKIqscl4WTW8QnPXw2ayzDgvVxXbVuYrsQSmMIj4byCHIAFh5KzYI5RhDyzQoAYGMDBBEOaggEvPNsiAYBO3Ax5xDiE9yDNoQpS0NapiVKEQW6XvAe22Q9akvd/TDdQa/YXT65AyDrP4hNf4ig94g8/48WcvL+jhe2nQqra1zK6Mnkytv/+rMmgVcPil6uhZwD4sBV518m4HjH8Kra2vH52+ri8XZr05vMAX8n+OLbynE5j1N+0qzwpnHfyo5OXvG/PzDuEqYeHX0SDln2P7DUrplLyQSucziWwmHGkEpmEGkjS3RcjCGuSgGPi/hFu4k1YlLrlSrV0qdYWaCfgW0vEnGCahtA==</latexit><latexit sha1_base64="DyBmRYxdA3RUZQKpjSxfJQOWn0g=">AAACZHichVFNSwJBGH7cvswsLQmCICQxOsloQtFJ6NLRj/wAE9ndRltcd5fdVTDpD9S16NCpICL6GV36Ax38A0F0NOjSodd1IUqqd5iZZ555n3eemZEMVbFsxnoeYWx8YnLKO+2b8c/OBYLzCwVLb5kyz8u6qpslSbS4qmg8byu2ykuGycWmpPKi1NgZ7Bfb3LQUXduzOwavNMW6ptQUWbSJymjVYITFmBPhURB3QQRupPXgLfZxAB0yWmiCQ4NNWIUIi1oZcTAYxFXQJc4kpDj7HMfwkbZFWZwyRGIbNNZpVXZZjdaDmpajlukUlbpJyjCi7IndsT57ZPfshX38Wqvr1Bh46dAsDbXcqAZOlnLv/6qaNNs4/FL96dlGDVuOV4W8Gw4zuIU81LePLvq57Wy0u8au2Sv5v2I99kA30Npv8k2GZy/how+I/3zuUVBIxOIbsUQmGUkl3a/wYhmrWKf33kQKu0gjT+dynOIM555nwS+EhMVhquBxNSF8C2HlE99kieQ=</latexit>n th-order interatomic force constant (IFC):<latexit sha1_base64="W1RiuunzGr0mjbyfX+5ylcnJDo0=">AAADLnichVFLaxRBEK4ZX3F9ZKMXwcvgEkkua08MRAyBgBePm8RJAulk6On0bprteTDTs5BM+g/4Bzx4UhARL/4HL/oDFHIQr4rHCF48WDM7bIjro4bp/uqr+qqqu4NEyUwTcmTZZ86eO39h4mLj0uUrVyebU9fWszhPufB4rOJ0M2CZUDISnpZaic0kFSwMlNgI+g/K+MZApJmMo0d6PxHbIetFsis500j5zfe0syf9goa5X7iG7sY6q3BkzAwVSlVsnyUJK9FilbA4DESjQGRmncaSQ7sp4wVNWKolUztIO54Z+U4+avOn0rOUl7WdsfSTSU41NDSQvd4hJhU5NUvE+M0WaZPKnHHg1qAFtXXi5kugsAsxcMghBAERaMQKGGT4bYELBBLktqFALkUkq7gAAw3U5pglMIMh28e1h95WzUbolzWzSs2xi8I/RaUD0+QDeUWOyTvymnwlP/9aq6hqlLPs4x4MtSLxJx/fWPvxX1WIu4a9E9U/Z9bQhXvVrBJnTyqmPAUf6gcHT47X7q9OF7fJc/IN539GjshbPEE0+M5frIjVp9DAB3B/v+5xsD7Xdu+251bmW8vz9VNMwE24BTN43wuwDA+hAx5w647lWTuWb7+xP9qf7M/DVNuqNdfhlNlffgFDg9sB</latexit>@nUµ1 ...µn (`1 1 ; . . . ; `n n ) =@uµ1 (`1 1 ) · · · @uµn (`n n ) {u}=014

Limitation of the harmonic approx.
Ĥ0 =<latexit sha1_base64="QArOEQDrw0J49OmlCf5dhn6fOJ0=">AAADCXichVFLTxRBEC7GF4LKghcTLxM2GIzJpmYxgZCYELlw5OECCYOT7qF3tsO86OndBDvzB/QHeOCkiQfjyZvxKBcTzxz4BxpOiokXDtTMTuKDgDWZ6aqv6qv5qounocw04uGAdenylavXBq8PDd+4eWukNjq2miVd5YuWn4SJWucsE6GMRUtLHYr1VAkW8VCs8e35Ir/WEyqTSfxE76ZiM2JBLNvSZ5ogr+a6HabNQu4ZzO1Htpt1I8/s5G6HM+UmkQhYGXIZTJaVPC/ip8bdYkEgVP4H+MBtK+YbJzfN3C4Y971aHRtYmn3WcSqnDpUtJrUv4MIWJOBDFyIQEIMmPwQGGT0b4ABCStgmGMIUebLMC8hhiLhdqhJUwQjdpm9A0UaFxhQXPbOS7dNfQnoVMW2YwAN8i8f4Gd/hNzw5t5cpexRadunkfa5IvZHnd1Z+/ZcV0amh85t1oWYNbZgptUrSnpZIMYXf5/eevTxemV2eMPfwNR6R/ld4iPs0Qdz76b9ZEst7F+jhNDejWFDP82+uqIuKPC3S+XdtZ53VZsOZajSXHtbnHlcrHYS7MA6TtLdpmIMFWIQW9f8EX+E7/LBeWO+tD9bHfqk1UHFuw19m7Z8CjlHAFw==</latexit>

X
q

1
†
~!q b̂q b̂q +2✗ Phonon is non-interacting✗ Phonon is volume-independent→ no thermal expansion✗ Phonon is temperature-independent→ no structural phase transition<latexit sha1_base64="uQeaQbAgjaNj0YTMpMo3e/RtVgM=">AAACtXichVG7SgNRED1ZXzE+ErURbIIhYhVuohCxCthY5mEeqCHsrjdxyb7c3QRiyA9YCxaioGAhfoaNYG2RTxDLCDYWzm6CokGdZffOnJkze+aOZKqK7TDW9Qkjo2PjE/7JwNT0zGwwNDdfsI2GJfO8bKiGVZJEm6uKzvOO4qi8ZFpc1CSVF6X6lpsvNrllK4a+47RMXtbEmq5UFVl0CNrdNzReEyvto04lFGEx5ln404n/dCIYWNoIPWIfBzAgowENHDoc8lWIsOnZQxwMJmFltAmzyFO8PEcHAeI2qIpThUhonb41ivYGqE6x29P22DL9RaXXImYYUfbEblmPPbA79szef+3V9nq4Wlp0Sn0uNyvBk8Xc278sjU4Hh1+sPzU7qGLD06qQdtND3CnkPr95fNbLbWaj7RV2zV5I/xXrsnuaQG++yjcZnj3/Q49Ec4sUc+r5+825dZqbp0UOrW3YKSRi8bVYIrMeSaUGK/VjCctYpb0lkcI20shTfx2nuMClkBTKwoFQ7ZcKvgFnAd9MMD4AppqczQ==</latexit><latexit sha1_base64="q101utA0uruwIUW76iQ25hO/kEU=">AAACxXichVHLSsNQED2Nr/qOuim4EYtSEcptFRRXBQVd+mhVsKUk8VqDeZncFmqp7v0BF64UXYjoT7gRXLvoJ4hLBTcunKQVFalOSO7MmTmTM3dUx9A9wVgtJLW0trV3hDu7unt6+/rlgcF1zy66Gs9otmG7m6ricUO3eEbowuCbjssVUzX4hro37+c3Stz1dNtKi7LDc6ZSsPQdXVMEQXk5krVNXlDylf3qZHaBGyJwY+mJvBxNxFlgI+yX85mKomHLtvyALLZhQ0MRJjgsCPINKPDo2UICDA5hOVQIc8nTgzxHFV3ELVIVpwqF0D36FijaaqAWxX5PL2Br9BeDXpeYIxhjj+yKvbB7ds2e2HvTXpWgh6+lTKda53In338cWXv7l2XSKbD7xfpTs8AOZgOtOml3AsSfQqvzSwcnL2tzq2OVcXbOnkn/GauxO5rAKr1qlyt89fQPPSrNrVDMqWfzm/PrTD//fZHNnfVkPDEVT65MR1OpxkrDGMYoYrS3GaSwhGVkqP8hLnCDW2lRMiUhleqlUqjBGcIPk44+AL2+omg=</latexit>!q +q (T )

Phonon spectrum→ infinite lifetime→ infinite thermal conductivity!q

2 q (T )<latexit sha1_base64="2BdExlXttN360VT0HDjixZ/NH3A=">AAACuXichVHLSsNQED2N7/po1Y3gpliUuim3VVB0U3ChS+1DCz7KTbxqNC+TtKDBH/AHFFy14EL8DDeCaxd+grhUcOPCSRpRkdYJyZ05M2dy5o5saarjMvYUkTo6u7p7evui/QODQ7H48Mi6Y1ZtRZQUUzPtsswdoamGKLmqq4myZQuuy5rYkI+W/PxGTdiOahpF98QS2zrfN9Q9VeEuQTvZrWWu67ziHZ+litOVeDKTZoEl2B/nK5VEaKtm/AFb2IUJBVXoEDDgkq+Bw6FnExkwWIRtwyPMJk8N8gJniBK3SlWCKjihR/Tdp2gzRA2K/Z5OwFboLxq9NjETmGSP7Ia9snt2y57ZR8teXtDD13JCp9zkCqsSOx8rvP/L0ul0cfDNaqvZxR7mA60qabcCxJ9CafJrpxevhYX8pDfFGuyF9NfZE7ujCYzam3K9JvJXbfTINDenWFDP1jfn1+l+/uciWzvr2XRmJp1dm03mcuFKezGOCaRob3PIYQWrKFF/G5eooyEtSlw6kA6bpVIk5Izil0nOJ9+DnaY=</latexit>

Energy15Anharmonicity (phonon-phonon interaction)Ĥ0 =<latexit sha1_base64="QArOEQDrw0J49OmlCf5dhn6fOJ0=">AAADCXichVFLTxRBEC7GF4LKghcTLxM2GIzJpmYxgZCYELlw5OECCYOT7qF3tsO86OndBDvzB/QHeOCkiQfjyZvxKBcTzxz4BxpOiokXDtTMTuKDgDWZ6aqv6qv5qounocw04uGAdenylavXBq8PDd+4eWukNjq2miVd5YuWn4SJWucsE6GMRUtLHYr1VAkW8VCs8e35Ir/WEyqTSfxE76ZiM2JBLNvSZ5ogr+a6HabNQu4ZzO1Htpt1I8/s5G6HM+UmkQhYGXIZTJaVPC/ip8bdYkEgVP4H+MBtK+YbJzfN3C4Y971aHRtYmn3WcSqnDpUtJrUv4MIWJOBDFyIQEIMmPwQGGT0b4ABCStgmGMIUebLMC8hhiLhdqhJUwQjdpm9A0UaFxhQXPbOS7dNfQnoVMW2YwAN8i8f4Gd/hNzw5t5cpexRadunkfa5IvZHnd1Z+/ZcV0amh85t1oWYNbZgptUrSnpZIMYXf5/eevTxemV2eMPfwNR6R/ld4iPs0Qdz76b9ZEst7F+jhNDejWFDP82+uqIuKPC3S+XdtZ53VZsOZajSXHtbnHlcrHYS7MA6TtLdpmIMFWIQW9f8EX+E7/LBeWO+tD9bHfqk1UHFuw19m7Z8CjlHAFw==</latexit>

X
q

1
†
~!q b̂q b̂q +2+Û3 + Û4 + · · ·<latexit sha1_base64="yR5FL83aMztyzTUSO1f6BlKRFsE=">AAACiXichVHLSsNAFD2Nr1pfVTeCm2KtCEKZaMHiSnDjUquthbaUJJ1qaJqEZFrQ0B/wB1y4UhARN25168YfcOEniMsKblx4mwZ8FPWGzJx75p47Z2ZU29BdwdhTSOrrHxgcCg9HRkbHxieik1M512o4Gs9qlmE5eVVxuaGbPCt0YfC87XClrhp8T61tdNb3mtxxdcvcFYc2L9WVfVOv6poiiCpH55dixQNFeNlW2Vtpxb5kKT/TKpZwy9E4SzI/Yr1ADkAcQWxZ0UsUUYEFDQ3UwWFCEDagwKWvABkMNnEleMQ5hHR/naOFCGkbVMWpQiG2RuM+ZYWANSnv9HR9tUa7GPQ7pIwhwR7ZFWuzB3bNntn7r708v0fHyyHNalfL7fLE8czO27+qOs0CB5+qPz0LVJH2verk3faZzim0rr55dNLeWcskvAV2zl7I/xl7Yvd0ArP5ql1s88wpIvQA8s/r7gW55aS8klzeTsXXU8FThDGLOSzSfa9iHZvYQpb2PcYNbnEnjUiylJbWuqVSKNBM41tIGx9UjpaB</latexit>✓ Phonon is interacting via anharmonicity✓ Phonon is volume-dependent→ finite thermal expansion✓ Phonon is temperature-dependent→ explains structural phase transition<latexit sha1_base64="uQeaQbAgjaNj0YTMpMo3e/RtVgM=">AAACtXichVG7SgNRED1ZXzE+ErURbIIhYhVuohCxCthY5mEeqCHsrjdxyb7c3QRiyA9YCxaioGAhfoaNYG2RTxDLCDYWzm6CokGdZffOnJkze+aOZKqK7TDW9Qkjo2PjE/7JwNT0zGwwNDdfsI2GJfO8bKiGVZJEm6uKzvOO4qi8ZFpc1CSVF6X6lpsvNrllK4a+47RMXtbEmq5UFVl0CNrdNzReEyvto04lFGEx5ln404n/dCIYWNoIPWIfBzAgowENHDoc8lWIsOnZQxwMJmFltAmzyFO8PEcHAeI2qIpThUhonb41ivYGqE6x29P22DL9RaXXImYYUfbEblmPPbA79szef+3V9nq4Wlp0Sn0uNyvBk8Xc278sjU4Hh1+sPzU7qGLD06qQdtND3CnkPr95fNbLbWaj7RV2zV5I/xXrsnuaQG++yjcZnj3/Q49Ec4sUc+r5+825dZqbp0UOrW3YKSRi8bVYIrMeSaUGK/VjCctYpb0lkcI20shTfx2nuMClkBTKwoFQ7ZcKvgFnAd9MMD4AppqczQ==</latexit><latexit sha1_base64="q101utA0uruwIUW76iQ25hO/kEU=">AAACxXichVHLSsNQED2Nr/qOuim4EYtSEcptFRRXBQVd+mhVsKUk8VqDeZncFmqp7v0BF64UXYjoT7gRXLvoJ4hLBTcunKQVFalOSO7MmTmTM3dUx9A9wVgtJLW0trV3hDu7unt6+/rlgcF1zy66Gs9otmG7m6ricUO3eEbowuCbjssVUzX4hro37+c3Stz1dNtKi7LDc6ZSsPQdXVMEQXk5krVNXlDylf3qZHaBGyJwY+mJvBxNxFlgI+yX85mKomHLtvyALLZhQ0MRJjgsCPINKPDo2UICDA5hOVQIc8nTgzxHFV3ELVIVpwqF0D36FijaaqAWxX5PL2Br9BeDXpeYIxhjj+yKvbB7ds2e2HvTXpWgh6+lTKda53In338cWXv7l2XSKbD7xfpTs8AOZgOtOml3AsSfQqvzSwcnL2tzq2OVcXbOnkn/GauxO5rAKr1qlyt89fQPPSrNrVDMqWfzm/PrTD//fZHNnfVkPDEVT65MR1OpxkrDGMYoYrS3GaSwhGVkqP8hLnCDW2lRMiUhleqlUqjBGcIPk44+AL2+omg=</latexit>!q +q (T )

Phonon spectrum→ finite lifetime→ finite thermal conductivity!q

2 q (T )<latexit sha1_base64="2BdExlXttN360VT0HDjixZ/NH3A=">AAACuXichVHLSsNQED2N7/po1Y3gpliUuim3VVB0U3ChS+1DCz7KTbxqNC+TtKDBH/AHFFy14EL8DDeCaxd+grhUcOPCSRpRkdYJyZ05M2dy5o5saarjMvYUkTo6u7p7evui/QODQ7H48Mi6Y1ZtRZQUUzPtsswdoamGKLmqq4myZQuuy5rYkI+W/PxGTdiOahpF98QS2zrfN9Q9VeEuQTvZrWWu67ziHZ+litOVeDKTZoEl2B/nK5VEaKtm/AFb2IUJBVXoEDDgkq+Bw6FnExkwWIRtwyPMJk8N8gJniBK3SlWCKjihR/Tdp2gzRA2K/Z5OwFboLxq9NjETmGSP7Ia9snt2y57ZR8teXtDD13JCp9zkCqsSOx8rvP/L0ul0cfDNaqvZxR7mA60qabcCxJ9CafJrpxevhYX8pDfFGuyF9NfZE7ujCYzam3K9JvJXbfTINDenWFDP1jfn1+l+/uciWzvr2XRmJp1dm03mcuFKezGOCaRob3PIYQWrKFF/G5eooyEtSlw6kA6bpVIk5Izil0nOJ9+DnaY=</latexit>

Energy16Vibrational many-body theory
Ĥ0 =

X

1
†
~!q b̂q b̂q +2q<latexit sha1_base64="QArOEQDrw0J49OmlCf5dhn6fOJ0=">AAADCXichVFLTxRBEC7GF4LKghcTLxM2GIzJpmYxgZCYELlw5OECCYOT7qF3tsO86OndBDvzB/QHeOCkiQfjyZvxKBcTzxz4BxpOiokXDtTMTuKDgDWZ6aqv6qv5qounocw04uGAdenylavXBq8PDd+4eWukNjq2miVd5YuWn4SJWucsE6GMRUtLHYr1VAkW8VCs8e35Ir/WEyqTSfxE76ZiM2JBLNvSZ5ogr+a6HabNQu4ZzO1Htpt1I8/s5G6HM+UmkQhYGXIZTJaVPC/ip8bdYkEgVP4H+MBtK+YbJzfN3C4Y971aHRtYmn3WcSqnDpUtJrUv4MIWJOBDFyIQEIMmPwQGGT0b4ABCStgmGMIUebLMC8hhiLhdqhJUwQjdpm9A0UaFxhQXPbOS7dNfQnoVMW2YwAN8i8f4Gd/hNzw5t5cpexRadunkfa5IvZHnd1Z+/ZcV0amh85t1oWYNbZgptUrSnpZIMYXf5/eevTxemV2eMPfwNR6R/ld4iPs0Qdz76b9ZEst7F+jhNDejWFDP82+uqIuKPC3S+XdtZ53VZsOZajSXHtbnHlcrHYS7MA6TtLdpmIMFWIQW9f8EX+E7/LBeWO+tD9bHfqk1UHFuw19m7Z8CjlHAFw==</latexit>+Û3 + Û4 + · · ·<latexit sha1_base64="yR5FL83aMztyzTUSO1f6BlKRFsE=">AAACiXichVHLSsNAFD2Nr1pfVTeCm2KtCEKZaMHiSnDjUquthbaUJJ1qaJqEZFrQ0B/wB1y4UhARN25168YfcOEniMsKblx4mwZ8FPWGzJx75p47Z2ZU29BdwdhTSOrrHxgcCg9HRkbHxieik1M512o4Gs9qlmE5eVVxuaGbPCt0YfC87XClrhp8T61tdNb3mtxxdcvcFYc2L9WVfVOv6poiiCpH55dixQNFeNlW2Vtpxb5kKT/TKpZwy9E4SzI/Yr1ADkAcQWxZ0UsUUYEFDQ3UwWFCEDagwKWvABkMNnEleMQ5hHR/naOFCGkbVMWpQiG2RuM+ZYWANSnv9HR9tUa7GPQ7pIwhwR7ZFWuzB3bNntn7r708v0fHyyHNalfL7fLE8czO27+qOs0CB5+qPz0LVJH2verk3faZzim0rr55dNLeWcskvAV2zl7I/xl7Yvd0ArP5ql1s88wpIvQA8s/r7gW55aS8klzeTsXXU8FThDGLOSzSfa9iHZvYQpb2PcYNbnEnjUiylJbWuqVSKNBM41tIGx9UjpaB</latexit>!q<latexit sha1_base64="uQeaQbAgjaNj0YTMpMo3e/RtVgM=">AAACtXichVG7SgNRED1ZXzE+ErURbIIhYhVuohCxCthY5mEeqCHsrjdxyb7c3QRiyA9YCxaioGAhfoaNYG2RTxDLCDYWzm6CokGdZffOnJkze+aOZKqK7TDW9Qkjo2PjE/7JwNT0zGwwNDdfsI2GJfO8bKiGVZJEm6uKzvOO4qi8ZFpc1CSVF6X6lpsvNrllK4a+47RMXtbEmq5UFVl0CNrdNzReEyvto04lFGEx5ln404n/dCIYWNoIPWIfBzAgowENHDoc8lWIsOnZQxwMJmFltAmzyFO8PEcHAeI2qIpThUhonb41ivYGqE6x29P22DL9RaXXImYYUfbEblmPPbA79szef+3V9nq4Wlp0Sn0uNyvBk8Xc278sjU4Hh1+sPzU7qGLD06qQdtND3CnkPr95fNbLbWaj7RV2zV5I/xXrsnuaQG++yjcZnj3/Q49Ec4sUc+r5+825dZqbp0UOrW3YKSRi8bVYIrMeSaUGK/VjCctYpb0lkcI20shTfx2nuMClkBTKwoFQ7ZcKvgFnAd9MMD4AppqczQ==</latexit><latexit sha1_base64="q101utA0uruwIUW76iQ25hO/kEU=">AAACxXichVHLSsNQED2Nr/qOuim4EYtSEcptFRRXBQVd+mhVsKUk8VqDeZncFmqp7v0BF64UXYjoT7gRXLvoJ4hLBTcunKQVFalOSO7MmTmTM3dUx9A9wVgtJLW0trV3hDu7unt6+/rlgcF1zy66Gs9otmG7m6ricUO3eEbowuCbjssVUzX4hro37+c3Stz1dNtKi7LDc6ZSsPQdXVMEQXk5krVNXlDylf3qZHaBGyJwY+mJvBxNxFlgI+yX85mKomHLtvyALLZhQ0MRJjgsCPINKPDo2UICDA5hOVQIc8nTgzxHFV3ELVIVpwqF0D36FijaaqAWxX5PL2Br9BeDXpeYIxhjj+yKvbB7ds2e2HvTXpWgh6+lTKda53In338cWXv7l2XSKbD7xfpTs8AOZgOtOml3AsSfQqvzSwcnL2tzq2OVcXbOnkn/GauxO5rAKr1qlyt89fQPPSrNrVDMqWfzm/PrTD//fZHNnfVkPDEVT65MR1OpxkrDGMYoYrS3GaSwhGVkqP8hLnCDW2lRMiUhleqlUqjBGcIPk44+AL2+omg=</latexit>!q +q (T )Bai

Dqla)

Phonon spectrum

Dyson equationem

2 q (T )<latexit sha1_base64="2BdExlXttN360VT0HDjixZ/NH3A=">AAACuXichVHLSsNQED2N7/po1Y3gpliUuim3VVB0U3ChS+1DCz7KTbxqNC+TtKDBH/AHFFy14EL8DDeCaxd+grhUcOPCSRpRkdYJyZ05M2dy5o5saarjMvYUkTo6u7p7evui/QODQ7H48Mi6Y1ZtRZQUUzPtsswdoamGKLmqq4myZQuuy5rYkI+W/PxGTdiOahpF98QS2zrfN9Q9VeEuQTvZrWWu67ziHZ+litOVeDKTZoEl2B/nK5VEaKtm/AFb2IUJBVXoEDDgkq+Bw6FnExkwWIRtwyPMJk8N8gJniBK3SlWCKjihR/Tdp2gzRA2K/Z5OwFboLxq9NjETmGSP7Ia9snt2y57ZR8teXtDD13JCp9zkCqsSOx8rvP/L0ul0cfDNaqvZxR7mA60qabcCxJ9CafJrpxevhYX8pDfFGuyF9NfZE7ujCYzam3K9JvJXbfTINDenWFDP1jfn1+l+/uciWzvr2XRmJp1dm03mcuFKezGOCaRob3PIYQWrKFF/G5eooyEtSlw6kA6bpVIk5Izil0nOJ9+DnaY=</latexit>

.

Dqlw)

EnergySpectral function

q (!) =1
Aq (!) = ImDq (!)
T⇡

Re⌃q (!)

q (!) = Im⌃q (!)
T<latexit sha1_base64="jFNPJflIRo02pN/64SJhHsb/7VI=">AAAC43ichVFNT9RQFD1URT4EBtmQuGmYYHAzeQUDxoQE0IXs+HCAhJLJa3kzvNDXlvbNJND0D+iOmLgQFpK4MP4L3JjA1gU/gbDExI0LbztNFAl4m/bde+49t+e+64SejDVjZx3Gnbv3Ou93dff0PujrHygNPlyJg2bkiqobeEG05vBYeNIXVS21J9bCSHDleGLV2X6R5VdbIopl4L/Wu6HYULzhy7p0uSaoVpqcrSU76ZgdKNHgT8xp065H3E2sNLFDmdqK661IJfMqNV/+XVgrlVmF5WZed6zCKaOwhaB0AhubCOCiCQUBH5p8DxwxPeuwwBAStoGEsIg8mecFUvQQt0lVgio4odv0bVC0XqA+xVnPOGe79BeP3oiYJkbZd/aZXbJv7As7Z79u7JXkPTItu3Q6ba4IawNvhpd//pel6NTY+sO6VbNGHc9yrZK0hzmSTeG2+a2995fLz5dGk8fsiF2Q/o/sjH2lCfzWD/fTolj6cIseh+bmFAvqefPNZXUqy9MirX/Xdt1ZGa9YE5XxxaflmblipV14hBGM0d6mMINXWECV+h/gGCc4NYTx1tg33rVLjY6CM4QrZhz+BtAerok=</latexit>

÷
•

I
.

T

.

•

T
•

T

.

p

17

of
*
It
of
*
It
of
*It
T

÷
•

Anharmonic self-energy“tadpole”

T

÷
•

.

T

q

÷
q

.

•

.

“loop”

T

I

•

T

T

T

.

T

.

“bubble”

TT

q

T

T
•

T

4
.

*

•

*

•

•

4

•

p

I

a.

.

•

''

*

T

T

t

•

t

-

t

a.

T
-

-

-

.

a.

•

tT

T

T

T

T

HH

X X

3-phononinteraction(cubic anharmonicity).

.

4-phononinteraction(quartic anharmonicity)
18

at
*It

÷÷
Bubble self-energy
T

I

•

.

I

.

T
T
.

q

.

q

T
T

•

•

•

.

T
T

p
p
.

''

I

''

Phonon scatteringby thethree-phonon
process
t
X.
t
*
t
*
t
•

I

•

*

T
T

a.

T
T
b
a.

Linewidth
p
'

q

-

-

•

-

•

-

-

-

-

.

-

.

Maradudin and Fein, Phys. Rev. 128,
2589 (1962).

•

:q

X.
H
X.
H

"

q

I

⇐

(in high-T region

÷

÷
.
.

got
got ..
qq

qq

)
of
of

X
#

a
a

'
'

"
"

Lifetime
19

Phonon lifetimes in Ge-based materials
d-Ge

Ge46

Ba8Ga16Ge30

10-fold
reduction

10x larger phonon-phonon scattering rate of heat-carrying acoustic
phonons induced by low-energy rattling guest modes in BGG
TT, Y. Gohda, and S. Tsuneyuki, PRL 114, 095501 (2015).

20

Lattice thermal conductivity
Fourier’s law j =<latexit sha1_base64="fn3aeVuEYPFZfpGA4iSalacNdFw=">AAACxnichVHLSsNQED3GV3226kZxEywVN5bbKiiCILpx6atVsKXcxGuNzYskrWgRXfsDLlwpuBBBP8KN4NpFP0FcKrhx4SQNiJbqhOTOnJkzOXNHsXXN9RirtUitbe0dnZGu7p7evv5obGAw61plRxUZ1dItZ0vhrtA1U2Q8zdPFlu0Ibii62FRKS35+syIcV7PMDe/QFnmDF01tV1O5R1AhNpJTjOr+sTwvT+ZK3La5nDO5onN5oxCLsyQLTG50UqETR2grVuwJOezAgooyDAiY8MjXweHSs40UGGzC8qgS5pCnBXmBY3QTt0xVgio4oSX6FinaDlGTYr+nG7BV+otOr0NMGQn2zG7YG3tkt+yFfTbtVQ16+FoO6VTqXGEXomfD6x//sgw6Pex9s/7U7GEXs4FWjbTbAeJPodb5laPzt/W5tUR1nF2xV9J/yWrsgSYwK+/q9apYu/hDj0Jzc4oF9Wx+c36d4edpkanfa2t0sulkaiqZXp2OLyyGK41gFGOYoL3NYAHLWEGG+p/gGne4l5YlUypLB/VSqSXkDOGHSadfwP2h9g==</latexit>

rT

Heat flux in “phonon gas” modelj ⇡ jph =<latexit sha1_base64="6Mdbbjkoo5uO5r6ldTUAEuD1j68=">AAADJHichVHPa9RAFP4aq/2htateCr0El4qnZbaKLYJQ8CL00h9uW2hKmMTZ3dhMMp3MLtaQf8B/wIMnCwrFk3+DFw+9SZFeSnsUjxV66cGXbERs6TphZt773ve9eS/PU2GQGMYOBqwrg1evDQ2PjF6/MXZzvHLr9koSd7QvGn4cxnrN44kIg0g0TGBCsaa04NILxaq3+TSPr3aFToI4em62ldiQvBUFzcDnhiC3Mu94Mn2ZOVwpHb+ye56bOpKbtpapameZ/cR2ko50063MaXtcO7EULV64xO5muRXlh1upshorln3RqJdGFeVaiCvHcPACMXx0ICEQwZAdgiOhbx11MCjCNpASpskKirhAhlHSdogliMEJ3aSzRd56iUbk5zmTQu3TKyFtTUobU2yf7bIT9pV9Yj/Y2aW50iJHXss23V5PK5Q7/mZi+fS/Kkm3Qfuvqm/NBk3MFrUGVLsqkLwLv6fvvn57svx4aSq9x3bYT6r/PTtgX6iDqPvL/7Aolt71qcejvjn5gnJe/udynuwb/9ORpmiTNnFp6PXzI75orEzX6g9q04sPq3OPyvEPYxJ3cZ9mPIM5PMMCGvTWZ3zHEY6tj9ae9c3a71GtgVJzB/8s6/A3is3EEA==</latexit>

Boltzmann transport equation (BTE)X~!q vq nq

q

✓
◆
✓
◆
◆
8 ✓
X
@nq X
@fk @nq>>vq · rT= = Pkk0 +Pkk0 q + · · ·>>@T@t scatt@t scatt
<
k0
k0 q
✓
◆
X
X
X
>
@n
>
k
>
=
P̃qq0 +
P̃qkk0 +P̃qq0 q00 + · · ·>:@t scatt
00 00
0<latexit sha1_base64="B5QL1JhiT9rxT143MnhGAPNwAQs=">AAADSnicjVHNahRBEK6ZjSaOP1n1EvAyuESSQ5bejagIQiAXj/nZTQLpMPT09m6a9PzY3buwDvMCvoAH8aDgQXwML4JePeQRxJsRRPRgzexI0JBoDTNd9VV9NV91hamSxhJy6Li1qXPnp2cueBcvXb4yW796bcskQ81Flycq0TshM0LJWHSttErspFqwKFRiOzxYLfLbI6GNTOKOHadiL2KDWPYlZxahoP6c9qRJFRsbO1bCX6JhlI3yIHuUU95LLI1ZqJjf8WkoB4MFn/Y14xlNmbaSKT8uCvPjuJNPChe9B/7Sf3Lsb06Q0YjZfR1lBtXZPPe8oN4gTVKaf9JpVU4DKltL6u+BQg8S4DCECATEYNFXwMDgswstIJAitgcZYho9WeYF5OAhd4hVAisYogf4HWC0W6ExxkVPU7I5/kXhq5Hpwzz5SF6TI/KOvCGfyM9Te2Vlj0LLGM9wwhVpMPtkbvPbP1kRnhb2j1lnarbQh3ulVona0xIppuAT/ujx06PN+xvz2S3yknxG/S/IIXmLE8Sjr/zVuth4doaeEOdmGAvsefrNFXVRkcdFtv5e20lnq91sLTfb67cbK3eqlc7ADbgJC7i3u7ACD2ENusCdmrPotJ1l94P7xf3u/piUuk7FuQ5/WG3qF1311Gc=</latexit><latexit sha1_base64="spXx6ZmzsThdl6pAac6si5bkIsI=">AAAEEniclVHNaxNBFH/b9aOuH031InhZDDUVMUyiqAhCwYvHtDVtoVPC7OwkHbLfMwnEZf8Bz4IHTwoexP/Am/Qi1KuH/gniSSoI6sG3m8SahlacZfe993vv995v9jmRJ5UmZM+YMU+cPHV69ox19tz5C3Ol+YtrKuwlXDR56IXJhsOU8GQgmlpqT2xEiWC+44l1p/swz6/3RaJkGDzWg0hs+awTyLbkTCPUmjdc6oiODFKOTVRmUVeqyGMDpQeesGlfRYyLtO77mXXTpo7sdBZt2k4YT2nEEi2ZZ7dbaTfLDmKdDQuvt1LqM72d+KnCcTrL7Ac2VT0f6yuZ3UCT2xt/sHgMxgXK3VArm9JJTZZ9lI7gP3QMR8Y4nmrpuSJtZBjGf8vpHkpOiI0rcWWKXEBj4RYVgTv+q61SmVRJcexppzZyyjA6jbC0CxRcCIFDD3wQEIBG3wMGCp9NqAGBCLEtSBFL0JNFXkAGFnJ7WCWwgiHaxW8Ho80RGmCc91QFm+MUD98EmTYskE/kDdknH8hb8pn8OrJXWvTItQzQOkOuiFpzTy+vfv8ny0erYfuAdaxmDW24V2iVqD0qkPwWfMjvP3m+v3p/ZSG9Rl6RL6j/JdkjO3iDoP+Nv14WKy+O0ePgvRnGAnse/efyOj/P4yJrh9c27azVq7Vb1fry7fLSndFKZ+EKXIVF3NtdWIJH0IAmcGPX+Gr8MH6az8x35ntzZ1g6Y4w4l2DimB9/A87GH+A=</latexit>qkk: n ⇡ n0 + n
q
q
q
1
0
nq = ~!qe
1<latexit sha1_base64="dmc6UdojLJwbraIkSay15k46K0I=">AAADQnichVFNb9NAEB2br9Z8NMClEhdDVISoFG0CAsSpUjlw7AdpK9XFWjuTsMraXrybiGD5D3CvOHACiQPiZ8ABiTOH/gTghILEgR4YO0ZVqFrWsnfmzbzntzuBkkIbxvYs+8TJU6fPzMw6Z8+dvzBXu3hpQyeDNMR2mMgk3Qq4RilibBthJG6pFHkUSNwM+stFfXOIqRZJ/MiMFO5EvBeLrgi5Iciv7XoB9kSchaShc8frCK0kH2kzkuh6Q614iFkrinKn62f93PW4UmnyzC2zxxnL3UXXe4DScLdq8KZFHDf2s6cHxDKbJpaQ42Hc+WvDr9VZg5XLPRw0q6AO1VpJap/Bgw4kEMIAIkCIwVAsgYOmZxuawEARtgMZYSlFoqwj5OAQd0BdSB2c0D59e5RtV2hMeaGpS3ZIf5H0psR0YYF9Ye/YmH1i79lXtn+kVlZqFF5GtAcTLip/7sX8+q//siLaDTw5YB3r2UAX7pVeBXlXJVKcIpzwh89fjtfvry1k19kb9p38v2Z77AOdIB7+DN+u4tqrY/wEdG5OOZLm0TdX9EVFnQbZ/Hdsh4ONVqN5q9FavV1fulONdAauwDW4QXO7C0vwEFagTfq/ravWTWvR/mh/s3/Y40mrbVWcyzC17P0/v8vSqQ==</latexit><latexit sha1_base64="+gDcJlU8tjqyV3B19UQhJp1vlKU=">AAAC3XichVFNT9RgEH6oiIgCq15MvGzcYLiwmQJRY2JC4sUjXwskFJq3dXZp6Jftu5tg06MXYzzKwZMmHgh/gZuRmHj2wE8wHCHx4sFptwlRAk7TvjPPzDN95h0n9r1UEx0NGFcGrw5dG74+cuPm6Nh47dbtlTTqJi633MiPkjVHpex7Ibe0p31eixNWgePzqrP9rMiv9jhJvShc1jsxbwSqE3ptz1VaILtGoZ29zDczyutP61Y7UW5m5hlvZpbDWllbjkqsKOCOKsryKTO3aw1qUmn1845ZOQ1UNh/VvsPCC0Rw0UUARggtvg+FVJ51mCDEgm0gEywRzyvzjBwjwu1KFUuFEnRbvh2J1is0lLjomZZsV/7iy5sIs44J+kF7dELfaJ9+0u8Le2Vlj0LLjpxOn8uxPf7m7tKv/7ICOTW2zliXatZo43Gp1RPtcYkUU7h9fu/V7snSk8WJ7AF9omPR/5GO6ItMEPZO3c8LvPjhEj2OzK0kZul58c0VdUGRl0Wa/67tvLMy3TRnmtMLs425h9VKh3EP9zEpe3uEOTzHPFrS/z0O8BWHhm28Nt4a7/qlxkDFuYO/zNj9AzUgrMs=</latexit>

q q

k
q

8
<linearization
fk ⇡ fk0 + fk

kʹ

q

qʹ

Relaxation-time approx.

qʺ<latexit sha1_base64="DtKOQ231aHMMY6CmcBA0Vt1eSfE=">AAADPXichVFLaxRBEP4yvuL6yKoXwcvgEokHl94oKoIQSA455uEmgUwYeia9myY9j/T0LsRh/oB/IIecEsgh5Ed48CLePeTowYMoIkQISg6pmZ0YNGTtoaervvq+6qouL1YyMYztD1gXLl66fGXwauXa9Rs3h6q3bs8lUUf7oulHKtILHk+EkqFoGmmUWIi14IGnxLy3Op7H57tCJzIKX5n1WCwFvB3KlvS5IcitLtuPbMeT7faI7bQ091Mn5tpIruzQTdey7NQ3WY/40E2dgJsVHaQJZTFZZtsvT8QTQhn+R2p4x13LKhW3WmN1Viz7rNEojRrKNRVVP8HBMiL46CCAQAhDtgJHQt8iGmCICVtCSpgmSxZxgQwV0naIJYjBCV2lf5u8xRINyc9zJoXap1sUbU1KG8PsI9tlB+w922Nf2NG5udIiR17LOp1eTytid+jN3dnD/6oCOg1WTlV9azZo4XlRq6Ta4wLJu/B7+u7rjYPZFzPD6QO2zb5S/Vtsn72jDsLuT39nWsxs9qnHo745+YJynv9yOS/oGz/pSFO0RZu4NPTGvyM+a8yN1huP66PTT2pjT8vxD+Ie7mOEZvwMY5jEFJp01wf8wC/8tt5an61v1vce1RooNXfw17IOjwGh68tl</latexit>

✓

@nq
@t

◆

=
scatt

nq
⌧q

BTE-RTA21

Lattice thermal conductivity from Peierls-Boltzmann theory
Boltzmann transport eq. (BTE)✓ Correct treatment of zero-point vibration✓ Reasonably accurate & efficient✗ Neglect 4th-order anharmonicity✗ Breakdown when an imaginary mode exists

22

LTC reduction by rattler

empty clathrate

20-fold reduction

filled clathrate

TT, Y. Gohda, and S. Tsuneyuki, PRL 114, 095501 (2015).

23

FIG. 1. (Color online) Lattice thermal conductivities calculated
at 300 K with respect to volume per formula unit.

First-principles calculation ofpart of the self-energy. For the Hamiltonian H the expansion
can be written as(1)

H = !0 + T + H 2 + H 3 + · · · ,• VASP code

with

BTE-RTA!
1PBEsol2•T =mκ [u̇xc-functionalα (lκ)] ,

(2)

2 lκα

where
! is theatomicannihilawave ve
vector W
a dynam• ALAMODE package

1 !!
H2 =
!αβ (lκ,l ′ κ ′ )uα (lκ)uβ (l ′ κ ′ ),2 lκα l ′ κ ′ βLines : calculation
Points: experiment

PbTe

Ba8Ga16Ge30

Bi2Te3

κ calculation (W/m-K)

Mg2Si

Si

103

with

BN

Zincblende-type
Wurtzite-type

BP
GaN

Dαβ (

SiC
BeO
AlN

GaP
AlAs
AlP InP
ZnO
ZnS GaSb
GaAs
ZnTe
InAs
CdS ZnSe
InSb

102

Using E
creationpart as a

AlSb

101
CuI

100

(3)

CdTe

and theA. Togo, L. Chaput, and I. Tanaka,collision
Phys. Rev. B 91, 094306 (2015).

CuBr
AgI
CuCl

100

phono3py software
101

102

κ experiment (W/m-K)

103

FIG. 2. (Color online) Comparison in lattice thermal conductivities at room temperature between experiments and calculations with
RTA and without isotope scattering.

H3 =Here thated by24
of intera

Softwares
Code
Developer
ShengBTE(almaBTE)

phonopy,phono3py

Language,
LicenseFeatures•N. Mingo(CEA, France)

Fortran (C++),
GPL

•
•

•

A. Togo(Kyoto)Python & C,
BSD 3-clause

•
•
•

•

PhonTSA.
Chernatynskiy(U. Florida)

Fortran, GPL

•
•
•

•

ALAMODE

T. Tadano(NIMS)•

C++, MIT•Interfaces

Iterative solution to BTE
Thermal transport innanostructure
MPI parallelization

VASP, QE, GULP

Direct solution to BTE
QHA
API for python and C
Easy to use

VASP, ABINIT, QE,
SIESTA, Elk, FHI-aims,
WIEN2k, CRYSTAL

Iterative and variational
solution to BTE
QHA
Includes empirical force fields
MPI parallelization

VASP, QE, LAMMPS

RTA
Self-consistent phonon
calculation
MPI + OpenMP parallelization

VASP, QE, xTAPP,
OpenMX, LAMMPS,
GULP

Outline
1. Introduction to phonon anharmonicity
2. Perturbative treatment of anharmonicity
3. Non-perturbative treatment of anharmonicity
4. Efficient calculation of anharmonic force constants
5. Applications
6. Summary

26

Limitation of the perturbative approach
The perturbation approachtends to become inaccurate when the anharmonicity is strong• breaks down when the harmonic phonon is unstable ( !q2 < 0 )•<latexit sha1_base64="r9TvkVoZo7yCU8+ovOsmtbhFG0s=">AAAC63ichVHLLgRRED3a+z3YSGwmJsRqUjMEEQuJjaXXIPGYdLc7ozP90t0zQWd+wNZCgkhIBPEZNn7AwtJSLElsLFT3tAhiVOf2rXuqTt1TtxRb11yP6KFGqq2rb2hsam5pbWvv6Ix1dS+5VtFRRUa1dMtZUWRX6JopMp7m6WLFdoRsKLpYVgrTQXy5JBxXs8xFb9cW64acN7WcpsoeQ9lY95pliLyc9bfLG366HJ+MUzaWoCSFFv/tpCIngchmrdgj1rAJCyqKMCBgwmNfhwyXv1WkQLAZW4fPmMOeFsYFymhhbpGzBGfIjBb4n+fTaoSafA5quiFb5Vt0Xg4z4xige7qmF7qjG3qi9z9r+WGNQMsu70qFK+xs537vwtu/LIN3D1tfrKqaPeQwHmrVWLsdIkEXaoVf2jt8WZiYH/AH6ZyeWf8ZPdAtd2CWXtWLOTF/XEWPwn3LfBZc8++XC/KMqvHPjhyO5nhxLg899XPEv52ldDI1nEzPjSSmRqPxN6EP/RjiGY9hCjOYRYbv2sEpLnElGdKBdCSdVFKlmojTg28mXXwAmDSpSA==</latexit>

Si

PbTe

HarmonicĤ0 =<latexit sha1_base64="QArOEQDrw0J49OmlCf5dhn6fOJ0=">AAADCXichVFLTxRBEC7GF4LKghcTLxM2GIzJpmYxgZCYELlw5OECCYOT7qF3tsO86OndBDvzB/QHeOCkiQfjyZvxKBcTzxz4BxpOiokXDtTMTuKDgDWZ6aqv6qv5qounocw04uGAdenylavXBq8PDd+4eWukNjq2miVd5YuWn4SJWucsE6GMRUtLHYr1VAkW8VCs8e35Ir/WEyqTSfxE76ZiM2JBLNvSZ5ogr+a6HabNQu4ZzO1Htpt1I8/s5G6HM+UmkQhYGXIZTJaVPC/ip8bdYkEgVP4H+MBtK+YbJzfN3C4Y971aHRtYmn3WcSqnDpUtJrUv4MIWJOBDFyIQEIMmPwQGGT0b4ABCStgmGMIUebLMC8hhiLhdqhJUwQjdpm9A0UaFxhQXPbOS7dNfQnoVMW2YwAN8i8f4Gd/hNzw5t5cpexRadunkfa5IvZHnd1Z+/ZcV0amh85t1oWYNbZgptUrSnpZIMYXf5/eevTxemV2eMPfwNR6R/ld4iPs0Qdz76b9ZEst7F+jhNDejWFDP82+uqIuKPC3S+XdtZ53VZsOZajSXHtbnHlcrHYS7MA6TtLdpmIMFWIQW9f8EX+E7/LBeWO+tD9bHfqk1UHFuw19m7Z8CjlHAFw==</latexit>

X

Cubic SrTiO3

Anharmonic1†~!q b̂q b̂q +q2“imaginary phonon”
27Temperature dependence of phonon frequency
Temperature dependence of phonon frequency 𝜔q and eigenvector eq is significant inenergy harvesting materials!

G. Shirane and Y.Yamada,
Phys. Rev.177, 858 (1969).

O. Delaire et al, Nature Mater. 10, 614 (2011)
K. Suekuni et al., PRB 81, 205207 (2010)28

Frequency renormalization
X1(a) ˜ ij (T ) = ij +ul)iijkl(huk3̄m
正方晶 (I4/mcm)
立方晶
P
m
4

(b)

thermal fluctuation<latexit sha1_base64="5PDZgMpVqnJZUeRozFdYrcgA9vE=">AAACtXichVG7SgNRED1Z3/GRqI1gI4aIVbhJBEUsBBvLmBgVNYS7641esi92NwFd/AFrwUIUFCzEz7ARrC38BLFUsLFwdrMgGqKz7N6ZM3Nmz9xRbV26HmPPMaWru6e3r38gPjg0PJJIjo5tuFbD0URZs3TL2VK5K3RpirInPV1s2Y7ghqqLTbW+EuQ3m8JxpWWue4e2qBh835Q1qXGPoO2Csatyx88fG9VkimVYaFPtTjZyUoisYCUfsYs9WNDQgAEBEx75OjhcenaQBYNNWAU+YQ55MswLHCNO3AZVCarghNbpu0/RToSaFAc93ZCt0V90eh1iTiHNntgte2MP7I69sM+OvfywR6DlkE61xRV2NXEyUfr4l2XQ6eHgm/WnZg81LIRaJWm3QySYQmvxm0dnb6XFYtqfYdfslfRfsWd2TxOYzXftZk0Uz//Qo9LcnGJBPTvfXFBnBHlaZPb32tqdjVwmm8/k1uZSy0vRSvsxiWnM0t7msYxVFFCm/iZOcYFLZV6pKHtKrVWqxCLOOH6YYn0BAC2cgQ==</latexit><latexit sha1_base64="2Nfz1nj920h18AOVRq8nY/MfYnA=">AAACsXichVFNLwNRFD3GV9VXsZHYNBpiVa8lIWIhsWGHapFqmHlea3S+MvPapBp/wMpOsCKxED/DRmJt0Z8gliQ2Fu5MJxEEdzLz7j33njvnvqs5hu5JxhotSmtbe0dnpCva3dPb1x8bGMx5dsXlIsttw3Y3NdUThm6JrNSlITYdV6imZogNrbzo5zeqwvV021qXNUcUTLVk6UWdq5Kg3PL0pMnNnViCJVlg8Z9OKnQSCG3Fjj1gG3uwwVGBCQELknwDKjx68kiBwSGsgDphLnl6kBc4QpS4FaoSVKESWqZviaJ8iFoU+z29gM3pLwa9LjHjGGOP7Ia9sHt2y57Y+6+96kEPX0uNTq3JFc5O//Fw5u1flkmnxP4n60/NEkXMBlp10u4EiD8Fb/Krh6cvmbm1sfo4u2LPpP+SNdgdTWBVX/n1qli7+EOPRnOrFAvq+fvN+XWmn6dFpr6v7aeTSydTU8n06nRiYT5caQQjGMUE7W0GC1jCCrLU/wAnOMO5MqVsKbuK1ixVWkLOEL6YUv4A0tWaXA==</latexit>

O<latexit sha1_base64="qhTsIL2L2FejlDegTZUJpGzqhWg=">AAADa3ichZHPaxNBFMffZv1R649GexH1sBgjFSXM1oIiCAEv4iltkraQDcvudJKOmf3B7mygDvMP+A948KTgQfwfvHgRPHvonyDiqYIXD76dJNhaUt+yM2/ee5+339kXpoLnkpB9q2KfOn3m7MK5xfMXLl5aql6+spknRUZZlyYiybbDIGeCx6wruRRsO81YEIWCbYWjJ2V+a8yynCdxR+6lrB8Fw5gPOA0khvzqR+WZJr1sGPaV2yDG7pGGeY1pT3Kxw5TX2uXaV/y5Xunc0c5j5whKDqETx9UlYQB91/EGWUCVq9Wa9vIi8tVIaGdWgP48HYdkhFkwYlIVyGpchNbar9Zmeee4M2tTg6m1kuoX8GAHEqBQQAQMYpDoCwggx6cHLhBIMdYHhbEMPW7yDDQsIltgFcOKAKMjXId46k2jMZ7LnrmhKX5F4Jsh6UCdfCXvyQH5TD6Qb+T33F7K9Ci17OEeTliW+ksvr7Z//ZeKcJew+5c6UbOEATw0WjlqT02kvAWd8OMXrw7ajzbq6jZ5S76j/jdkn3zCG8Tjn/TdOtt4fYKeEO8d4Jlhz/l/rqyLyjwO0v13bMedzdWGe7+xur5WazanI12A63ATVnBuD6AJT6EFXaDWLeuZ1bY6lR/2sn3NvjEprVhTZhmOmF3/A9g41zw=</latexit>

Sr

kl

Ti

2

⌦ /T
<latexit sha1_base64="4ekRC13T5Ly+Ik6FlfHJVa5WbyM=">AAAC0XichVE9SyNRFD3O+rV+ZtdmwWYwRmwML3FhF6uAzXabmEQFE8Ob8RkH54uZl0AcAmKZ/QFbbLWChdjYamsjWFvkJyxbKthYeGcyICrqHWbevefec+fcdzXXNHzJWLdP+dA/MDg0/HFkdGx8YjLx6fOq7zQ8XZR1x3S8dY37wjRsUZaGNMW66wluaaZY03aXw/xaU3i+4dgl2XJF1eJ129g2dC4JqiVmKz8tUeebQbatVlzPcaWjlhbUUi2oWFzueFagt9u1RJKlWWTqSycTO0nElncSV6hgCw50NGBBwIYk3wSHT88GMmBwCasiIMwjz4jyAm2MELdBVYIqOKG79K1TtBGjNsVhTz9i6/QXk16PmCpS7Jodsxt2yU7YP3b/aq8g6hFqadGp9bjCrU12vhTv3mVZdErsPLLe1Cyxje+RVoO0uxESTqH3+M293zfFpZVUMMcO2X/S/5d12QVNYDdv9aOCWPnzhh6N5uYUC+r5+s2FdVaYp0Vmnq/tpbOaTWcW09nC12QuF690GNOYwTzt7Rty+IE8ytS/g1Oc4VwpKi1lXznolSp9MWcKT0z59QBhr6dt</latexit>

hQR25 i =6 0<latexit sha1_base64="fz2G7azXf9awvbJHreZ2c0/A8SQ=">AAACz3ichVG/T9tQEP5w+d1C0nZB6hIRUXWKLoGKqmJA6tKRBAKRCLKezQOs+Fftl1RguerY9h9gYAKpA0JipXsXJOYO+ROqjlTq0qFnxxIqUeAs+919d9/5u3eGb1uhIuoOaQ+GR0bHxicmHz6ams7lHz9ZD712YMq66dle0DBEKG3LlXVlKVs2/EAKx7DlhtF6k+Q3OjIILc9dU/u+3HLErmvtWKZQDOn52aYRiJZUUVWPmo5Qe4ET1WI9qryM47jpyncF0vNFKlFqhX6nnDlFZLbi5a/QxDY8mGjDgYQLxb4NgZCfTZRB8BnbQsRYwJ6V5iViTDK3zVWSKwSjLf7ucrSZoS7HSc8wZZv8F5vfgJkFzNEPOqVruqQz+kl/B/aK0h6Jln0+jR5X+nruy8zqn3tZDp8KezesOzUr7OBVqtVi7X6KJFOYPX7n4PB69XVtLnpOJ/SL9R9Tl77zBG7nt/m1KmtHd+gxeG7BseSeg28uqXOSPC+yfHtt/c56pVSeL1WqC8XlpWyl43iGWbzgvS1iGW+xgjr3/4RzXOCbVtXeax+0j71SbSjjPMV/pn3+B2+0pw4=</latexit><latexit sha1_base64="IWcLZ2gJYTRolKYQLJluzxHzguw=">AAACx3ichVFNLwNRFD0d39/FRtiIhlg1ryUhVsKGFS3VJkrzZrzWxHxlZtqoSRO2/oCFFQmJWPAfbCTWFn6CWFZiY+HOtIJIuZN5795z77lz7lzZ0lTHZewpJDU1t7S2tXd0dnX39PaF+wc2HLNoKyKlmJppZ2TuCE01RMpVXU1kLFtwXdZEWt5b9PPpkrAd1TTW3bIltnReMNS8qnCXoFx4OJu3ueLFKl68kl3RRYFvk5fwj1w4EouywEbZL+czFUHdVs3wA7LYgQkFRegQMOCSr4HDoWcTMTBYhG3BI8wmTw3yAhV0ErdIVYIqOKF7dBYo2qyjBsV+TydgK/QVjV6bmKMYZ4/silXZPbtmz+y9YS8v6OFrKdMt17jCyvUdD629/cvS6Xax+8X6U7OLPGYDrSpptwLEn0Kp8UsHJ9W1ueS4N8HO2QvpP2NP7I4mMEqvykVCJE//0CPT3JxiQT0b/zm/Tvfz3xfZ2NmIR2NT0XhiOjK/UF9pO0Ywhkna2wzmsYRVpKj/IS5xg1tpWTKlkrRfK5VCdc4gfph09AEn5qNk</latexit>

Tc

hQR25 i = 0<latexit sha1_base64="WtlUfBgPNvyKnurly8ItwMn8r7s=">AAACx3ichVFNLwNRFD3GV9VXsRE2oiFWzWtJiJWwYaeqSNpq3ozXmpivzEwbTJqw9QcsrEhIxIL/YCOxtuhPEEsSGwt3pg1BWncy79177j13zp0rW5rquIxVW6TWtvaOzlBXuLunt68/MjC44ZglWxFpxdRMe0vmjtBUQ6Rd1dXElmULrsua2JT3lvz8ZlnYjmoa6+6BJXI6LxpqQVW4S1A+MpIt2Fzx4hUvUcmauijybfKS/pGPRFmMBTb25cR/O1HUbdWMPCKLHZhQUIIOAQMu+Ro4HHoyiIPBIiwHjzCbPDXIC1QQJm6JqgRVcEL36CxSlKmjBsV+TydgK/QVjV6bmGOYYE/smr2yB3bDntlHw15e0MPXckC3XOMKK99/Mpx6/5el0+1i95vVVLOLAuYCrSpptwLEn0Kp8cuHp6+p+bUJb5JdsBfSf86q7J4mMMpvymVSrJ010SPT3JxiQT0b/zm/TvfztMg/a/vrbCRi8elYIjkTXVisrzSEUYxjivY2iwUsYxVp6n+EK9ziTlqRTKks7ddKpZY6Zwg/TDr+BHJ0o4Q=</latexit>

1 2 2
⌦ Q
2

1 2 2
! Q
2

<latexit sha1_base64="uTCHZNgtuEhsVTeBqcvQpvbpqUw=">AAACzXichVFNS9xgEH6M9bvVbXspeAkuFk/L7KpYSgXBS2+6q6uCK+FNfNWw+SJ5d8HGeKzgH/DgScGDtNf2D/RS6LkHf0LxaKEXD06yAbGiTkjemWfmmTzzjhk4dqSILrq07mc9vX39A4NDz18MjxRevlqJ/FZoybrlO364ZopIOrYn68pWjlwLQilc05GrZnM+za+2ZRjZvresdgO54Yptz96yLaEYMgp6wwxFU6q4asQNV6id0I1riRFXppMk0Wd1MgpFKlFm+n2nnDtF5LboF36hgU34sNCCCwkPin0HAhE/6yiDEDC2gZixkD07y0skGGRui6skVwhGm/zd5mg9Rz2O055Rxrb4Lw6/ITN1jNNvOqcr+klf6A9dP9grznqkWnb5NDtcGRgjh2+W/j3JcvlU2LllPapZYQvvMq02aw8yJJ3C6vDbn46ult7XxuO3dEqXrP+ELugHT+C1/1pnVVk7fkSPyXMLjiX3fPjm0jo3zfMiy/+v7b6zUimVJ0uV6lRx7kO+0n6MYgwTvLcZzOEjFlHn/p/xFd/wXVvQWtqett8p1bpyzmvcMe3gBickpbc=</latexit>

QR25
Q<latexit sha1_base64="CVGGTcjKZO30lVuns1t2JD5Suvg=">AAAC2HichVE7SwNBEP5yvuI7aiPYBEPEKmxUVKwCNpZ5GA3GEO7OjR65F3eXQAyCnaidKBZWChbiz7DxD1hYWoplBBsL5y4noiFxjr2d/Wa+2W92JFNVbIex54DQ1d3T2xfsHxgcGh4ZDY2Nb9hGxZJ5VjZUw8pJos1VRedZR3FUnjMtLmqSyjel8qob36xyy1YMfd2pmbygibu6UlJk0SEolSqGIizGPAu3OnHficC3pBF6wTZ2YEBGBRo4dDjkqxBh05dHHAwmYQXUCbPIU7w4xwEGiFuhLE4ZIqFl+u/SKe+jOp3dmrbHlukWlZZFzDCi7IndsQZ7ZPfslX22rVX3arhaarRLTS43i6Mnk5mPf1ka7Q72flgdNTsoYdnTqpB200PcLuQmv7p/2cispKP1GXbD3kj/NXtmD9SBXn2Xb1M8fdVBj0R9i3TmVLP9y7l5Wsf4d0cWRUu0KJeGHv874lZnYy4Wn4/NpRYiiUV//EFMYRqzNOMlJLCGJLJ0F8cpznEhbAmHwpFw3EwVAj5nAr9MOPsCYcaiXQ==</latexit><latexit sha1_base64="WrG2IyClBwk6HXEMru3AXafIiMU=">AAACwHichVHLSsNQED2Nr/pstRvBTbEorsptVRRXpW5c2taqoBKSeLWheZncFtrQH/AHXLgQBUHxM9wIrl34CeJSwY0LJ2lAVNQJyZ05M2dy5o7qGLonGHuMST29ff0D8cGh4ZHRsURyfGLTsxuuxquabdjutqp43NAtXhW6MPi243LFVA2+pdZXg/xWk7ueblsbouXwPVM5tPQDXVMEQXIyVZL9XVMRNdf0yx3Zzy92OnIyw7IstPRPJxc5GUS2bifvsYt92NDQgAkOC4J8Awo8enaQA4ND2B58wlzy9DDP0cEQcRtUxalCIbRO30OKdiLUojjo6YVsjf5i0OsSM40Z9sCu2Qu7Yzfsib3/2ssPewRaWnSqXS535MTxZOXtX5ZJp0Dtk/WnZoEDLIdaddLuhEgwhdblN9snL5WV8ow/yy7YM+k/Z4/sliawmq/aZYmXT//Qo9LcCsWcev5+c0GdGeRpkbnva/vpbOazuflsvrSQKRSjlcYxhWnM0d6WUMAa1lGl/i2c4QrXUlGqSbZ01C2VYhEnhS8mtT8AgPCgwQ==</latexit>

Self-consistent phonon theory:
Stochastic implementation : Errea et al., PRB 2014
Deterministic implementation: TT and S. Tsuneyuki, PRB 2015
Effective force constant from MD: Hellman et al., PRB 201329

Self-consistent phonon theory
Effective one-body Hamiltonian : H̃0 =<latexit sha1_base64="oGBddNXx0XKboUKeYRQ2L1DySpI=">AAAC+HichVFLa9VQEJ7GVx9qr7oRugleKhXhMikFiyAU3BRc2Ie3LTRtOMmdmx6alyfnXmgP+QP+ARdurCBSutK/4EZwLdKdLovLCm5cOMkNFC2tE5KZ+Wa+yXfO+Fkkc414OGRduHjp8pXhkdGxq9eujzdu3FzJ054KqB2kUarWfJFTJBNqa6kjWssUidiPaNXfflzWV/ukcpkmz/RORhuxCBPZlYHQDHmNJ66WUYfMfOEZLOxHtpv3Ys88L9wtXyj3aUyhKFN7qnKbxu2IMCRVDND7bleJwDiFmS7ueY0mtrAy+3Tg1EETaltIG5/BhQ6kEEAPYiBIQHMcgYCcn3VwACFjbAMMY4ojWdUJChhlbo+7iDsEo9v8DTlbr9GE83JmXrED/kvEr2KmDZP4BffxGD/hAR7h7zNnmWpGqWWHvT/gUuaNv7i9/Ou/rJi9hq0T1rmaNXRhttIqWXtWIeUpggG/v/vyePnh0qS5i2/wB+vfw0P8yCdI+j+Dt4u09OocPT6fW3BOPPPsmyv74rLOi3T+XdvpYGW65WDLWZxpzs3WKx2GCbgDU7y3BzAH87AAbZ7/Hr7CN/hu7VqvrXfW/qDVGqo5t+Avsz78AcwHt/Q=</latexit>

N. R. Werthamer, Phys. Rev. B 1, 572 (1970)X

1†~⌦q (aq aq + )2

qLet us introduce an operator ρ defined as ⇢ =<latexit sha1_base64="NHK/aShZwifOwtROxBgNvZNb5RE=">AAAC+HichVHPS9xQEB5jtVZb3dpLoZfQxeKly0QEpVAQehF68deqYOzyEmfdh/nFy9sFfeQf6D/QQy9toYh4qv+Cl0LPpXhrj6VHBS8eOskGRGV1QvJmvplv8s0bLwlkqhFP+qz+ewOD94ceDI88fDQ6Vnk8vprGbeVT3Y+DWK17IqVARlTXUge0nigSoRfQmrfzJs+vdUilMo5W9G5Cm6HYjmRT+kIz1Ki8dVUrtl/bblMJ39A789L1SAvb1TLYIjOfNQxmWWbcUOiWCs2KynoWNSpVrGFh9k3HKZ0qlLYQV36AC1sQgw9tCIEgAs1+AAJSfjbAAYSEsU0wjCn2ZJEnyGCYuW2uIq4QjO7wd5ujjRKNOM57pgXb578E/Cpm2jCBP/EAT/E7HuJfvOjZyxQ9ci27fHpdLiWNsfdPl8/vZIV8amhdsm7VrKEJs4VWydqTAsmn8Lv8zt6H0+VXSxPmBX7Bf6z/M57gMU8Qdc78r4u09PEWPR7PLTgm7tn75vK6MM/zIp3ra7vprE7VHKw5i9PVudlypUPwDJ7DJO9tBuZgHhagzv2/wS/4DX+sPeuTtW8ddEutvpLzBK6YdfQfvIa4ZQ==</latexit>

e
Tre

H̃0
H̃0

Then, the following Gibbs inequality holds:
F [⇢] = F̃0 + hH
<latexit sha1_base64="qnZl+5r40PntJ69k+UMCVAg3fYw=">AAAC6nichVFNa9RgEH4av2r96KqXgpfgUhHEZSKCRSgUhGWP/XDbwmZZknS6fdk3HybvLtSwf8CjF5Ge6sdB/BleFG+Kh/4E8VjBiwcn2YBoaZ2QvDPPM89k5h0/0SozRAdT1qnTZ86emz4/c+HipcuztStX17N4mAbcDmIdp5u+l7FWEbeNMpo3k5S90Ne84Q8eFvzGiNNMxdEjs5twN/T6kdpWgWcE6tUWmx033Ym79qLtGqW3OG+OezmN7du266fegE3eulMxrZKZ0G6fH9vNXq1ODSrNPuo4lVNHZctx7RNcbCFGgCFCMCIY8TU8ZPJ04ICQCNZFLlgqnip5xhgzoh1KFkuGJ+hAvn2JOhUaSVzUzEp1IH/R8qaitDFPX+ktHdIHekff6NextfKyRtHLrpz+RMtJb/bp3NrP/6pCOQ12/qhO7NlgGwtlr0p6T0qkmCKY6EdPnh+uPVidz2/SS/ou/e/TAb2XCaLRj+DNCq/undCPL3N7ErPUPP7miryw4GWRzr9rO+qs32041HBW7tWXFqqVTuM6buCW7O0+ltDCMtpS/xU+4jO+WNp6Zr2w9iap1lSluYa/zHr9G5PesJ0=</latexit>

H̃0 i0F

The self-consistent phonon (SCP) frequency is obtained via<latexit sha1_base64="9qCDAdrfimgU3xkg0OUFAQJiZNE=">AAAC03ichVHPSxtREP5cbf3RH0ntRfCyNCg9hVkRFKEgCNKjJk0UkhDeri9xyf5i9yVgl1yKJ9FzDz210EPx5tXeeil47iF/QulRwUsPzm5WSiuxs+y+mW/mm/3mjRk4dqSIBmPa+MSDh5NT0zOPHj95mss/m61Gfje0ZMXyHT/cNUUkHduTFWUrR+4GoRSu6cgds7OR5Hd6Moxs33ujDgLZcEXbs1u2JRRDzfxivRUKK67vSUcJfbNWD/f9Rv82TqK+/kqnZr5ARUpNv+sYmVNAZlt+/gJ17MGHhS5cSHhQ7DsQiPipwQAhYKyBmLGQPTvNS/Qxw9wuV0muEIx2+NvmqJahHsdJzyhlW/wXh9+QmToW6Ad9oUv6Tqf0k36P7BWnPRItB3yaQ64MmrmjufL1f1kunwr7f1j3alZoYTXVarP2IEWSKawhv/f2/WV5rbQQL9In+sX6P9KAvvEEXu/K+rwtSx/u0WPy3IJjyT1H31xS5yZ5XqTx79ruOtWlokFFY3u5sL6arXQK83iBl7y3FazjNbZQ4f7HOMM5vmoVLdbeaYfDUm0s4zzHX6ad3ACJk6fF</latexit>

F [⇢]
=0
⇢

“Hatree-Fock” theory for phonons. (mean-field approximation)

30

Self-consistent phonon theory
How to update ρ (effective Hamiltonian)?
Stochastic self-consistent harmonic approximation (SSCHA)
I. Errea et al., Phys. Rev. B 89, 064302 (2014).

Use conjugate gradient method to update ρ and F[ρ] byrepeatedly calculating atomic forces in supercells.
• The gradient of F[ρ] is calculated stochastically.
• Full anharmonicity included•

Self-consistent phonon approach based on force constants (SCP)
TT and S. Tsuneyuki, Phys. Rev. B 92, 054301 (2015).
•

Compute δF[ρ]/δρ = 0 analytically and derive the self-consistentequation by truncating the higher-order anharmonic terms.
hH<latexit sha1_base64="NnASG8sEhQKPuDjgUjnJ0Zfa1ko=">AAADNnichVHNS9xQEJ/EWj9bV0uhUNDgogjSZaKCIhSEXvboR1cFI+El+7SPzRfJ20Ub9g+w/0APnlroofQf6L2XgmcPXryJFA8KCl48OHkbFCvqhLyZ+c38fpmXcSJPJBLxQNPbnrU/7+js6u7pffGyr9A/sJyE9djlFTf0wnjVYQn3RMArUkiPr0YxZ77j8RWn9iGrrzR4nIgw+Ci3I77us81AbAiXSYLswo7lxKzGZVp+Z0nhVXlabtopNtVhvDduyiq/22OMGxU7ncz9VOataigTQ1UtFkVxuPWEgiKq2C4UsYTKjPuBmQdFyG0+LOyBBVUIwYU6+MAhAEmxBwwSetbABISIsHVICYspEqrOoQndxK1TF6cORmiNzk3K1nI0oDzTTBTbpa949MbENGAE9/EnnuNf/IX/8OpBrVRpZLNsk3daXB7ZfV/eLF0+yfLJS/h0y3p0ZgkbMKNmFTR7pJDsFm6L3/j89XxpdnEkHcXveELzf8MD/EM3CBoX7o8Fvrj7yDwO3ZtRzknz4T+X9flZnRZp/r+2+8HyRMnEkrkwVZybyVfaCW9hGMZob9MwB2WYhwrpn2qvtUFtSP+tH+pH+nGrVddyziu4Y/rZNZS0y/k=</latexit>

H̃0 i0 = hH0

H̃0 + U3 + U4 + . . .i0 ⇡ hH0

H̃0 + U4 i0

Efficient calculation of ρ and F[ρ] at various temperature• Higher-order terms (sixth, eight, …) are neglected.
•

31Self-energy in SCP theory
We only consider the frequency-independent first-order self-energies (loopand tadpole):
Self-energy in perturbation theory

**

T

T

55
L•
L•

,

,

T

T

I
I

•

•

T

•

T

•

•

•

T

TT

T

T

T

tadpoleloop

Self-energy in self-consistent phonon theory
-

-

÷÷
¥
i*¥
÷÷
i.*.ff

**

.

T

=

=

T
µ
.

.

T

t

.

.

t

.

µ

T
.

T

T
.

.

T
.

If

.

.

+

T

ItIt →
y→
y
.

If

T
.

t

.

.

.

T
+

Phonon version of HartreeFock approximation.
cf. Fetter & Walecka chap. 4

. .
.

t

"

.

"

;I÷::
t

::

;I÷
t

+

•

.

•

+

T

;+

.in

T

.

. •

;

•

,

t

..

.
. •

•

,

.

t

..

.

.

+

.in

.

32

Temperature dependent phonon in cubic SrTiO3

Green line:
SCP result@300K
Red circle:
INS result@300 K

HSE hybrid functional
Anharmonic stabilization of soft modes
TT and S. Tsuneyuki, J. Phys. Soc. Jpn. 87, 041015 (2018)33

DicusWhat about the bubble self-energy?
Dglw )The SCP propagator does not include the bubble self-energy, which gives thelowest-order contribution to the phonon line width. How can we include it?
"

i÷n;÷=÷+÷:+:÷
:.int.im:7.

Bf

Dglw )

Dicus

900
800

Frequency (cm−1 )700"i÷n;÷=÷+÷:+

600

:

500

1
Aq (!) = ImDq (!)⇡

400<latexit sha1_base64="jFNPJflIRo02pN/64SJhHsb/7VI=">AAAC43ichVFNT9RQFD1URT4EBtmQuGmYYHAzeQUDxoQE0IXs+HCAhJLJa3kzvNDXlvbNJND0D+iOmLgQFpK4MP4L3JjA1gU/gbDExI0LbztNFAl4m/bde+49t+e+64SejDVjZx3Gnbv3Ou93dff0PujrHygNPlyJg2bkiqobeEG05vBYeNIXVS21J9bCSHDleGLV2X6R5VdbIopl4L/Wu6HYULzhy7p0uSaoVpqcrSU76ZgdKNHgT8xp065H3E2sNLFDmdqK661IJfMqNV/+XVgrlVmF5WZed6zCKaOwhaB0AhubCOCiCQUBH5p8DxwxPeuwwBAStoGEsIg8mecFUvQQt0lVgio4odv0bVC0XqA+xVnPOGe79BeP3oiYJkbZd/aZXbJv7As7Z79u7JXkPTItu3Q6ba4IawNvhpd//pel6NTY+sO6VbNGHc9yrZK0hzmSTeG2+a2995fLz5dGk8fsiF2Q/o/sjH2lCfzWD/fTolj6cIseh+bmFAvqefPNZXUqy9MirX/Xdt1ZGa9YE5XxxaflmblipV14hBGM0d6mMINXWECV+h/gGCc4NYTx1tg33rVLjY6CM4QrZhz+BtAerok=</latexit>

300
200

Green line:
SCP result@300K

100
0
Γ

X

M

Γ

R

M

34

Improved theory for lattice thermal conductivity

SCP + BTE

µ⌫
ph =

Crystalstructure

1 X
cq (T )vqµ (T )vq⌫ (T )˜⌧q (T )Nq V qSelf-consistentphonon calc.

loop (quartic)
direct method

2nd4th
order
2nd,
order
force constant

3rd order
force constant•• phonon
T-dep. phonon•• group
T-dep.vel.
group vel.

Lattice thermal
conductivity

lifetime

bubble (cubic)
35

Phonon transport in cubic SrTiO3
SCP + BTE

µ⌫
ph =

1 X
cq (T )vqµ (T )vq⌫ (T )˜⌧q (T )Nq V q

New theory(cubic & quartic terms)
Theory(only cubic term)

TT and S. Tsuneyuki, J. Phys. Soc. Jpn. 87, 041015 (2018)36

Outline
1. Introduction to phonon anharmonicity
2. Perturbative treatment of anharmonicity
3. Non-perturbative treatment of anharmonicity
4. Efficient calculation of anharmonic force constants
5. Applications
6. Summary37

Anharmonic lattice model (ALM)
Taylor expansion of potential energy:
VALM = V2 + V3 + V4 + · · ·1X
1 X=
ij ui uj +
2 i,j
3!

i,j,k

<latexit sha1_base64="4dD5J1t6Vw3JP8LzcAdmmuthi2k=">AAAE3HicnVI7b9NQFD41Bkp4NIUFicUQtUJqVF0/kAEJqYiFAaS0JWmlOFi2e5O4fsqPSMW6GxOIFQYmkBgQP4CJiYU/wNCZCQFTkVgYOHYcoKloaa7lcz6fe75zvutzzdC144SQ7SnuCH/02PHpE5WTp06fmanOnm3FQRpZtGkFbhCtm0ZMXdunzcROXLoeRtTwTJeumc6tfH9tQKPYDvx7yVZIO57R8+2ubRkJhvRZrtXSM80zkn7kZTfv3GWsMn9DwJjEhIXcy6VXcq9ZG0ESC5qWJ2Va0b4d9cxORhZJseq/gVgCpnUjw8pEhiWZFqeentn1TZaN0dUrknhNRrqoSqqq5ECSRVlkWqNvI2WTMfa/HVPMzw1yKguH1ylf/Eto3Rlve6BSZxKlaJzJ5Cq75dY16rqH11ywJtWNZsivjO6IXq2NWMJeMKpTg3I1gupb0GADArAgBQ8o+JAgdsGAGJ82iEAgxFgHMoxFiOxinwKDCnJTzKKYYWDUQdvDr3YZ9fE7rxkXbAu7uPhGyBRgjnwkr8kO+UDekM/k5z9rZUWNXMsWenPIpaE+8+j86o8DWR76BPp/WPtqTqALVwutNmoPi0h+CmvIHzx4trN6fWUumycvyRfU/4Jsk/d4An/w3Xq1TFee76PHRC34x3BA4vg49oKWtCjiXVlWaktKOappuACX4DLOQ4UluA0NaILFveM+cV+5b/x9/iH/mH8yTOWmSs452LX4p78Aa5VNWw==</latexit>

parameters=b·Parameters:
maximum order, cutoff radii

1 X
ijk ui uj uk +4!

i,j,k,`

ijk` ui uj uk u` + · · ·basis<latexit sha1_base64="dazZZfy5YGSQDcP29I4L/zpjVfE=">AAAC2XichVE9axRRFD0Zv2L8yKpNwGZxiVgsy53dldGAEEhjufnYJLCzLPMmL5tH5ouZtwtxmCKdpLWwsFKwEH+AYGvjH7BIZ5ukjGBj4Z3ZATHBeIeZd96599w5710ReSrRRIdTxqXLV65em74+c+PmrduzlTt315NwFLuy64ZeGG8KJ5GeCmRXK+3JzSiWji88uSF2l/L8xljGiQqDNb0Xyb7vDAO1rVxHMzWorD2rpnbRphcPRT+lOtXNzBZ+KrLMdrdCfSbfsB43zaetOjVMq2lZ7Rw0W2ZrIrI7OyrLBpUaNaiI6nlglqCGMjph5RNsbCGEixF8SATQjD04SPjpwQQhYq6PlLmYkSryEhlmWDviKskVDrO7/B3yrleyAe/znkmhdvkvHr8xK6uYp2/0gU7pK32kI/r1z15p0SP3ssermGhlNJg9mFv9+V+Vz6vGzh/VhZ41tvGk8KrYe1Qw+SnciX784vXp6sLKfPqQ3tEJ+39Lh/SFTxCMf7jvl+XKmwv8CPbCN8YDMs+O4zxYbzZMnu5yu7bYLkc1jft4gEc8DwuLeI4Outz/M77jCMdGz9g3XhoHk1JjqtTcw19hvPoNfjWuOQ==</latexit>

M: The number of independent parameters
FALM =<latexit sha1_base64="l4Ru8317G+XTYpZ2PLtddpvGfmM=">AAADE3ichVFNTxRBEK0ZUWH9YNGLkcvEFePFTY0hwZiQQEyMB0mWj11IGNj0DL1sh/lKT+8mOOk/wIEb8cAJEw+GH8CBm3rwwpEDP8F4hMSLB2tmJxIXxJ7M1OtX9Wped7mxLxKFeGKY1wau37g5OFS6dfvO3eHyyL1GEnWkx+te5EdyyWUJ90XI60oony/FkrPA9fmiu/Eqyy92uUxEFC6ozZivBGw9FC3hMUVUs7ztuEH6WjdTJ2CqLYN0+u2M1tak9cxpSealTsykEsy3Gn0l+jyVtehoXSKR1afKUq5eTRcuqc+CU2sLPTn9BzbLFaxivqyLwC5ABYpVi8oH4MAaROBBBwLgEIIi7AODhJ5lsAEhJm4FUuIkIZHnOWgokbZDVZwqGLEb9F2n3XLBhrTPeia52qO/+PRKUlowhsf4CU/xG+7jd/z1z15p3iPzsknR7Wl53BzeejD/87+qgKKC9rnqSs8KWvAi9yrIe5wz2Sm8nr777v3p/Mu5sfQJfsAf5H8PT/ALnSDsnnkfZ/nc7hV+XPJCN0YDsvvHcRE0nldtrNqz45Wp8WJUgzAKj+ApzWMCpuAN1KBO/c+Mh0bFeGzumIfmZ/Nrr9Q0Cs19+GuZR78BR9DI6Q==</latexit>

@VALM
=
@u
2

<latexit sha1_base64="+7U/RwvLDBl7CZoBAa0rkl7EGOQ=">AAAECHichVE9T9tAGH5j0xbSjwS6ULG4jYjoQHR2kYoqIYG6MGTgowEkTCP7cgnX+Ev2JSq1vDDyBzp0aiUGxNINdYWlf6ADPwF1pFKXDn3jGJKYht7Jvufee573fe5e07N4IAg5z0jyyJ2790bHsvcfPHyUy49PbARuy6esQl3L9bdMI2AWd1hFcGGxLc9nhm1abNNsvu6cb7aZH3DXeSP2PLZjGw2H1zk1BIaq45n9JWVB0U3W4E5o2obw+fsoO6u0qiGP3oYz6vNIKSqzet03aKhGoRb1nSB89y/Si6dDWAibPYFOa64IFF3vq6cNraf1MqVJffUGWVf1tFS9a1zEOWigPNRAuZc6TeozMMi6MlBOGdCZU7t+7mq+QEokHspNoCagAMlYcfMnoEMNXKDQAhsYOCAQW2BAgHMbVCDgYWwHQoz5iHh8ziCCLGpbyGLIMDDaxH8Dd9tJ1MF9J2cQqylWsfDzUanANPlBjsgl+U6OyQX5MzRXGOfoeNnD1exqmVfNHUyu//6vysZVwG5PdatnAXWYj71y9O7Fkc4taFff/vDxcv3V2nRYJF/IT/T/mZyTM7yB0/5FD1fZ2qdb/JjoBV8MG6Sm23ETbGgllZTU1bnC4lzSqlGYgmcwg/14CYuwDCtQAZq5kHLSpPREPpC/yifyty5VyiSaxzAw5NO/KRkWEw==</latexit>

6
6
A=6
4

@bT
@u

=A

(1)
ui
(2)
ui

1 (1) (1)
2 ui uj
1 (2) (2)2 ui uj

(L)

ui

···
1 (L) (L)
2 ui uj

1 (1) (1) (1)3! ui uj uk
1 (2) (2) (2)3! ui uj uk
1 (L) (L) (L)3! ui uj uk

3

···
7
· · ·77
5···38

Parameter estimation by least squares
The parameter vector is determined so that the ALM force best reproduces the DFT
force (a.k.a. force matching method).

All parameters can be determined uniquely if matrix A is overdetermined.
Necessary displacement patterns can be generated systematically.
@2V
=
ij =
@ui @uj@Fj⇡@ui[Fj (ui = h)Fj (ui =
2h

h)]

@3V
@Fk
=
=
ijk
@ui @uj @uk@ui @uj1⇡[Fk (ui = h, uj = h) Fk (ui = h, uj = h)24h
Fk (ui = h, uj = h) + Fk (ui = h, uj = h)]2 patterns

4 patterns

39

Rapid increase of displacement patterns
Silicon
FC6

104

FC4

FC5

# static calculations# static calculations

104

NaCl

103

FC3

102

101

FC6

FC4

FC5

103

FC3
102

101

FC2

FC2
10

0

100
1

2

3

4

# shell’s cuto↵

5

6

7

1

2

3

4

5

6

7

# shell’s cuto↵

More than 1,000 DFT calculations will be necessary to compute higher-orderanharmonic terms (>3) if the finite-displacement approach is employed.

40Compressive sensing
LASSO

˜ = arg min kA

PRL 113, 185501 (2014)

F DFT k22 + k k1

PHYSICAL REVIEW LETTERS

week ending
31 OCTOBER 2014

Lattice Anharmonicity and Thermal Conductivity from Compressive Sensingof First-Principles Calculations
Fei Zhou (周非)Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, California 94550, USA

Weston Nielson, Yi Xia, and Vidvuds OzoliņšDepartment of Materials Science and Engineering, University of California, Los Angeles, California 90095-1595, USA(Received 22 April 2014; published 27 October 2014)
First-principles prediction of lattice thermal conductivity κ L of strongly anharmonic crystals is a longstanding challenge in solid-state physics. Making use of recent advances in information science, we
propose a systematic and rigorous approach to this problem, compressive sensing lattice dynamics.
Compressive sensing is used to select the physically important terms in the lattice dynamics model anddetermine their values in one shot. Nonintuitively, high accuracy is achieved when the model is trainedon first-principles forces in quasirandom atomic configurations. The method is demonstrated for Si, NaCl,
and Cu12 Sb4 S13 , an earth-abundant thermoelectric with strong phonon-phonon interactions that limit theroom-temperature κL to values near the amorphous limit.

•

Avoid over-fitting•

Improve the generalization ability•

DOI: 10.1103/PhysRevLett.113.185501

To a large extent, thermal properties of crystalline solids
are determined by the vibrations of their constituent atoms.
Hence, an accurate description of lattice dynamics isessential for fundamental understanding of the structure,
thermodynamics, phase stability, and thermal transport
properties of solids. The seminal work of Born and
Huang [1] forms the theoretical basis of our understandingof harmonic vibrations and their relation to elastic properties. With the advent of efficient density-functional-theory(DFT-)based methods for solving Schrödinger’s equation,
several ab initio methods for studying harmonic phonon
properties of solids have been proposed, such as thefrozen phonon approach [2,3], supercell small displacement method [4,5], and the density-functional perturbation
theory (DFPT) [6]. Because of these developments,
ab initio calculations of the harmonic phonon dispersion
curves and phonon mode Grüneisen parameters havebecome routine.
A systematic approach to anharmonicity has been more

G. James, D. Witten, T. Hastie, and R. Tibshirani,“An Introduction to Statistical Learning with
Applications in R”

PACS numbers: 63.20.Ry, 63.20.dk, 66.70.−faccuracy of these methods are well established [11–14].
Unfortunately, PT tends to be computationally expensive
for solids with large, complex unit cells, and its ability tohandle strong anharmonicity is insufficient, especiallywhen the harmonic phonon dispersion contains imaginary
frequencies or when phonon scattering becomes so intensethat κ L saturates at its theoretical minimum [15]. Manyinteresting and technologically relevant materials belong
to this class, e.g., ferroelectrics and thermoelectrics withultralow κ L [16]. In these cases, a general and efficient
nonperturbative approach that can accurately describe
four-phonon and higher-order interactions is needed. The“2n þ 1” theorem of DFPT [17] can be used to calculate
the fourth- and higher-order terms, but the computations
are cumbersome and require specialized codes which, tothe best of our knowledge, are not available for n > 1.
In this Letter, we introduce an approach to building
lattice dynamical models that can treat compounds withlarge, complex unit cells and strong anharmonicity, includ-

Improve the computational efficiency

Each row of the sensing matrix A should be uncorrelated.
Use DFT-MD + random displacement.
41

Workflow of compressive sensing
Silicon (64 atoms), VASP (PBEsol)
Add 5 new configurations per each iteration
5-fold CV
2747 total features••

CV score (mean, dev)

6
4

Converge in 4-5th
iterations(20-25 DFT
calculations).

2

1500

1000

L @ 300 K (W/mK)# non-zero features

Error in forces (%)

•

130

120

110

1

2

3

4

Iteration

5

6

42

Cubic SrTiO3
2x2x2 supercell (40 atoms)• 3,018 anharmonic force constants (3rd, 4th, 5th, 6th)•

finite displacement methodwith central difference

more than
5,000 DFT calc.
100x more efficient!

ab initio MD + LASSO

40 DFT calc.

Accurate DFT calculations based on hybrid functions arealso possible
43Model selection (SrTiO3)
˜ = arg min kA

F DFT k22 + k k1

Relative error

0.5
Training
Validation

0.4
0.3

SrTiO3, energy

0.2
0.1
0.0

2500

∥Φ∥0

2000
1500
1000
500
0 −4
10

10−3
opt.

10−2

10−1

λ

• 40 training structures• 2024 coefficients out of 3018 are nonzero (~67%)

SrTiO3, force

44

Cs
Pb

o-

r3

3

3

Br

bB

t-C
sP

Br

3

3

iO

rT

Cs
Pb

c-

t-S

O

Ti

Sr

c-

Se

0

e3

Sn

6G

a1

8G

Ba

Te

Pb

3

F

Sc

105

Si

No. of static DFT calcs.

Efficiency improvement
100x speedup is possible

Conventional
LASSO

104

103

102

45

Outline
1. Introduction to phonon anharmonicity
2. Perturbative treatment of anharmonicity
3. Non-perturbative treatment of anharmonicity
4. Efficient calculation of anharmonic force constants
5. Applications
6. Summary

46

Intermetallic clathrate• A8X16Y30 (A= Ba, Sr, Eu; X, Y = Si, Ga, Ge, Sn)• sp3 host framework & guest atom
• “Rattling motion” of guest atoms
Large amplitude (~0.4 Å) thermal motion
Strong anharmonicity
Type-I clathrate• High thermoelectric figure-of-merit ZT (~1.3)• Low lattice thermal conductivity (< 1 W/mK)

47

b 60

d 400

Intensity (cts min–1)

Failure of conventional theory at high-T–1

q = 0.25 Å : (6, 0.31, 0.31)

300

40

2

3

5
4
E (meV)

E1

200

Peierls-Boltzmann
theory200

q = 0.87 Å–1: (6, 1.06, 1.06)
E2

E3

E4

100
6

7

cubic anharmonicity for•0 0 Consider
4
8
12
16
calculating
phonon lifetimes
E (meV). 6 Transverse acoustic (TA) phonons proﬁles measured at 10 and 300 K in Ba7.81Ge40.67Au5.33. a–c Comparison between constant q-scans of the TA100
011
Neglect
higher-orderanharmonicity
anch measured at 10 K (full blue circles) and 300 K (red squares) on the triple axis spectrometer (called 2T@LLB) with an effective resolution of 200 μeV.
imilar comparison between scans with a larger resolution of 1.2 meV measured at the wave vector (6, 1.06, 1.06) close to the centre of the BZ where the
10 1 modes, labeled E1,2,3,4 in the text, are clearly visible (a, b). In a–d, the intensities are corrected for the Bose population factor. Vertical error
ur optic phonons are s.d. from counting statistics
Expt. (Sales, PRB2001)
Ba7.81Ge40.67Au5.33
Lattice thermal conductivity (W/mK)•

Kph (Wm–1K–1)

e long wavelength limit roughly scales with ω2 (see Fig. 5 and10
Expt.
pplementary Note 5).
Experiment(May, PRB2009)
Simulation
The origin of the discrepancy on the temperature dependence
Simulation
10 0
d the absolutevalues of phonon lifetimes and thermalnductivity needs thus to be sought elsewhere, and requires
rther measurements of both temperature and energy depen1
Ab
initionce of phonons lifetime, coupled to atomistic or molecularsimulation
namics simulations taking into account all the present results,
1/T
well as a possibly
force ﬁeld.
Batemperature-dependent
8Ga16Ge30
Indeed, the10 −1 1/T decay in simulations is related to the
10 1
10 2
10 3
se–Einstein phonon
population factor,
with the assumption
0.1
Temperature (K)10100
at the force ﬁeld does not change with temperature (see ref. 9 ).
T (K)this assumption is released, the temperature dependence of themulated thermal
conductivity
could
be different.
Ab initio
data: TT,
Y. Gohda,
and S. Tsuneyuki,
F. Loryconductivity
et al., Natincommun.
8, Au
491
(2017)
Fig. 7 LatticeP.thermal
Ba7.81Ge40.67
5.33. Comparison
Concerning theabsolute
value of(2015).
phonon lifetimes and thermal between the experiment (full circles) and the anharmonic three phonon
PRL.
114, 095501
nductivity, the current phonon Umklapp simulations might simulation (red dotted line) carried out using the experimental lattice
derestimate them. Indeed is worth noticing that a slight change parameter. The blue solid line includes boundary effect with a size of 500
the lattice parameter used in simulations can strongly affect them, nm. The dashed black line indicates the 1/T decay. Other scattering48
shown by Tadano et al.10 Going from the relaxed (1.10 nm) to theF0Quartic
anharmonicity in clathrates
X

Γ

X

Γ

X

X

(b)
1.4

Harmonic +T2g mode
Quartic +SexticHarmonic +++2%Harmonic
Quartic
Quartic +Sextic

1.0

XX

(c)

Quartic +SexticHarmonic +Harmonic
Quartic

1.2X(c) 8 XX16 Γ 30

-2%Harmonic +opt.

Harmonic +Quartic

Γ

Ba Ga Ge (BGG)

Potential energy surface
XX Γ

20
X
Γ

E/N (eV)

Γ

0.8
20

0.6

0

10

20

10 +4% 20

01

η (u 2 Å)

0.4

e dynamics of type-I clathrate Ba8 Ga16 Ge30 with various lattice constants. (a) Phonon dispersion curves in thehe Γ–X line. The dotted lines show the harmonic phonon frequencies and the solid lines are the result of the SCP
0.2 “×” indicate the frequencies of the+6%
ross symbols
lowest T 2g modes. (b) Potential energy surface of the lowest
epresent the results of DFT calculations and the lines are calculated using the harmonic and anharmonic force
Cage size
re of Ba0.0
pattern of the T 2g mode.
8 Ga16 Ge30 . The red arrows indicate the displacement-2%
+6%

20

00

10
10

20

1
we also calculated the potentialen2
η (umodeÅ)byRaman active T 2g guest= M −1/2 e η with η being the am-

0

10

20

0

10

1
tic anharmonicity generally
increase the frequencies of low2
⌘ (u Intuitively,Å)lying rattling modes.
this can be attributed to thedominant and positive contribution from the diagonal term of

20
49

Γ

Finite-temperature phonon dispersion
Calculated with the optimized lattice constant (VASP, PBE)

Gray lines: Harmonic phonon dispersion
Red lines : Anharmonic phonon dispersion

2

Phonon frequency shift in Ba8Ga16Ge30(a)80Frequency (cm−1)

60
40
20
0
X

−20

Γ

Γ

X
X

Γ

X

Γ

X

(b) (c)X

ΓX

Phonon frequency (c)(cm-1)

1.4

Harmonic

1.2

Harmonic +Quartic

1.0

E/N (eV)Γ

LatticeHarmonic +Quartic + Harmonic
constants

-2%

36.2

41.5

opt.

28.6

35.5

+2%

19.7

29.5

+4%

7.0

23.9

20
15.3i

18.8

0.8
0.6
0.4

SexticSCP (300K)

0.2
0.0
0

T10
2g mode
20

0

10

+6%0

20
1
2

η (u Å)

1051

Change of high-temperature behavior of κL
SCP + BTE
µ⌫
ph =

1
Nq V

• Consider cubic anharmonicity forXcalculating phonon lifetimes

cq (T )vqµ (T )vq⌫ (T )˜⌧q (T )

• Additionally consider quarticqanharmonicity for frequency
renormalization

Lattice thermal conductivity (W/mK)

10 1

Thermal conductivity (300 K)
100

SCP+BTE

Conventional
BTE

Conventional
BTE

SCP+BTE

Expt.

0.58

0.97
?

1.06

Ba8Ga16Ge30
10 −1

10 1

10 2

10 3

Temperature (K)TT and S. Tsuneyuki, PRL 120, 105901 (2018)52

ational modes(d)at the (f) M and (g) R points involve coordinated staggered rotations of octahedra in the patterns shown. The rhombohedral
se can be described as a static tilt according to the R4+ pattern.

Description of negative thermal expansion

1.2
gely preserve the internal dimensions of the stiff metal-anion
lecular subunits. In order to retain local constraints of bondances and angles, staggered rotations induce a shrinking
he cell dimensions [Figs. 1(c) and 1(d) show this for the+
distortion], establishing a long-range coupling between
al transverse linkage
3 fluctuations and lattice volume thatlong been ascribed as the cause of structural NTE [16–19].
R is an important mode branch in perovskites because
mperature-driven condensation of modes on this branch
describe many structural transitions which accompany

Purposeof ourofstudy
ordering
transitions
other types [1–3]. These soft modes
circumscribe the BZ edges [Fig. 1(e)] within a very narrowwindow of energy set by the M-R branch dispersion, with the
R point ≃200 µeV lower at all temperatures measured, consis(a) symmetry-lowered rhombohedral
tent with an approach to theR3c structure, as observed in other 3d BF3 systems.
The rhombohedral R3c structure is related to the hightemperature cubic P m3m structure via a staggered octahedral
tilt around the ⟨111⟩ axis [Fig. 1(a) inset], corresponding toa frozen R4+ lattice distortion [Fig. 1(g)]. The staggered 111

ScF shows large negative thermal expansion,
which can be attributed to the rigid unit modes(RUMs) of fluorine octahedra.

(a)(g)(b)(b)

Figure1.3: (a): Structure of metal trifluorides. This structure is similar to the perovskite but there

PRB 92, 134101 (2015)

is no A-site atoms. Fluorides form the octahedra around metals. (b): Lattice parameter and the

The quasi-harmonic approx. (QHA) breaks
down due to strong anharmonic effects of
1.2 Purpose of our study
the RUMs.

CTE (ppm/K)

5
undary vs the B-site mean radius in 3d
CTE of ScF3 [5]. The black line is fitted to the lattice parameter and the red dashed line is the linears. [5,8–11]. (b) Simple cubic structure
CTE.determined
The greenline isthethe
diﬀerencethe experimentally determined lattice parameter and
FIG. 2. (Color online) (a) Lattice volume as a function of temperaturefrom tracking
Braggpeak (3,0,0)between
in single-crystalline
terofregularcorner-linked
octahedral
0
3 . Inset: The lattice volume over a larger thermal window determined on powder samples [7]. (b) Overview of the lattice mode dispersion
fit. by arrows. The longitudinal (LA) and transverse (TA)placementof Strongmodesbranch
attributed
electmomenta patternand energies.
softening along thethe
M-Rspline
cut is shown
usticsimple
branchescubic
are shownnear ". Blacksolid lines
are a guide to the eye. Red circles in the inset mark the momenta R + δ and M + ϵ.
hestructure,
showingthe5π,π ) points. Soft collective rigid unit
134101-2
as indicated by shadowing. The lowest
10
Present study (LDA)
he patterns shown. The rhombohedral
Present study (PBE)
Present study (PBEsol)

15

Li et al. (LDA)er types [1–3]. These soft modes
Liu et al. (PBE)
Expt.
[Fig. 1(e)] within a very narrow
From the previous studies, it is known that anharmonicity is important for understanding
NTE in20
e M-R branch dispersion, with the
0
200
400
600
800
1000
1200
ScF3 . However, how each-order anharmonicity plays a role in NTE
is unclear
and the microscopic
all temperatures measured, consisTemperature(K)symmetry-lowered rhombohedral
comprehension needed for designing and controlling NTE is still lacking. Hence, in this study,
in other 3d BF3 Y.
systems.
Ooba, TT, R. Akashi and S. Tsuneyuki, Phys. Rev. Materials 3, 033601 (2019). 53
structure is related to the high-Failure of QHA: Pressure induced softening
Pressure induced softening (calculated with PBEsol)

Frequency@ R (1/2,1/2,1/2)

800

6

700

4

Frequency (meV)

Phonon freq. (cm 1 )

600

500

400

300

LDA
PBE
PBEsol

2
0
2
4
6
8

200

10
100

62

64

66

68

3

Volume (Å )
0

100

X

M

R

X|M

R

Y. Ooba, TT, R. Akashi and S. Tsuneyuki, Phys. Rev. Materials 3, 033601 (2019).

54

Failure of QHA: Equation of states(QH)
Fvib (V, T ) =F (V, T ) = E0 (V ) + Fvib (V, T )<latexit sha1_base64="Pup2avc8LDeg3Wk6NcOT4Z2Ea8s=">AAAC1nichVG/SxxREP5cY2LUxDM2Qpolh+EO5ZiVgEcgIIiS0l93Hqgsu+tTH7e/2H13oMulE8HGdClSJZAipLY1RZqAtYV/gqQ0YGOR2b2FkIg6y+6b+Wa+2W/e2KErY0V03qP1Puh7+Kj/8cDg0JOnw4WRZ/U4aEWOqDmBG0QN24qFK31RU1K5ohFGwvJsV6zazdk0v9oWUSwDf0XthmLDs7Z9uSUdSzFkFsrzpfrkSll/o8+ZCXVK9bI+oc+bybpnqZ3IS9rS7nSyErNQpAplpt90jNwpIreFoHCKdWwigIMWPAj4UOy7sBDzswYDhJCxDSSMRezJLC/QwQBzW1wluMJitMnfbY7WctTnOO0ZZ2yH/+LyGzFTxzid0Ve6pJ/0jS7o+tZeSdYj1bLLp93litAcPhxbvrqX5fGpsPOXdadmhS1UM62StYcZkk7hdPntvQ+Xy6+XxpOX9Jl+sf5PdE4/eAK//dv5siiWPt6hx+a5LY4F97z95tI6L83zIo3/13bTqU9VDKoYi6+KM9V8pf14jhco8d6mMYO3WECN+x/hGCf4rjW0d9q+dtAt1Xpyzij+Me39H/wHpqQ=</latexit>

1X

ln 2 sinhq<latexit sha1_base64="XEMqtAfnQc/d2qsd1P3u3jUdCL4=">AAADNnichVHNa9RAFH+JX3X96KoIgqDBpZIFWSaLYBGEgiA9dtvutrCzhpk4mwzNl5PZhRryB+g/4MGTggfxH/DuRfDsoRdvIuJBoYIXD75kQ0VL6wvJe+/33u/lN/N4GspME7JjmEeOHjt+Yu5k49TpM2fnm+fOD7JkojzR95IwUZucZSKUsehrqUOxmSrBIh6KDb51t6xvTIXKZBKv6+1UjCLmx3IsPaYRcpuPrXtuTiOmAxXlU8mL4v5eaveW20VhD26st607Fh0r5uVOkVMuNCtoNonc/GFBw7hh0VCM9dDq0kzGAeXS9+299m5REWjAmaJJJHxW0uxB26oa0SnpB3rkNlukQyqz9gdOHbSgtpWk+R4oPIAEPJhABAJi0BiHwCDDZwgOEEgRG0GOmMJIVnUBBTSQO8EugR0M0S38+pgNazTGvJyZVWwP/xLiq5BpwQL5QF6RXfKOvCZfyK8DZ+XVjFLLNno+44rUnX9yae3nf1kReg3BH9ahmjWMYbHSKlF7WiHlKbwZf/ro6e7a7dWF/Dp5Qb6i/udkh7zFE8TTH97Lnlh9dogejudmmAucefDNlX1RWcdFOv+ubX8w6HYc0nF6N1tLi/VK5+AyXAMb93YLlmAZVqCP878ZF40rxlXzjfnR/GR+nrWaRs25AH+Z+f03P73NZw==</latexit>✓1~!q (V )2◆

28.08028.235

Imaginarymode exists

28.24028.085

Free energy (eV)28.24528.090

No imaginarymodes within
HA

28.25028.09528.25528.10028.26028.10528.26528.11028.2700 KK
TT == 300

28.11528.275
63.0

63.5

64.0

64.5

65.0

65.5

66.0

66.5

Volume (Å3 )Y. Ooba, TT, R. Akashi and S. Tsuneyuki, Phys. Rev. Materials 3, 033601 (2019).

55

SCP free energy(SCP)
Fvib (V, T ) =

1X
q



ln 2 sinh

1 Xh 2⌦q (V, T )4 q<latexit sha1_base64="CV2PJVjUmkRG1aSSQNIC4it4Ct0=">AAADz3ichVHNa9RQEJ80ftT1o1u9CF6eLpVdsMvLUrAIQnVBehDc7X60sG+7vGTfJmnzZZJdqSHiUf0HPHhS8CCCV717ETx76J8gHit48eAkGxu1tL6QZOY385v5zRvVs8wgpHRXmpGPHT9xcvZU4fSZs+fmivPnu4E79jXR0VzL9TdUHgjLdEQnNENLbHi+4LZqiXV1u57E1yfCD0zXaYc7nujbXHfMkanxEKHBvHT5ziBiNg8N344mphrHm1H5t9+qN+JKXO5ea1fITcJGPtciJY6YKkIeExaM7UF0PyYFZjmEWWIU9kgNYdMxInbb1Ms5pRanJGao3Gf3bKHzhImVSbuSpFawnG/qRtgnjBUWc+JSvN8m65DTN7EsmapbJOU6zqHaCCPOhlzXhR+zu3gTQ/5HZH+21Vs4Ws6pDKKtrZjkKrjlGWmX6qBYolWaHnLQUDKjBNlpuMXPwGAILmgwBhsEOBCibQGHAJ8eKEDBQ6wPEWI+WmYaFxBDAbljzBKYwRHdxq+OXi9DHfSTmkHK1rCLha+PTAIL9At9Q/foJ/qWfqU/D60VpTUSLTv4V6dc4Q3mnl1s/fgvy8Z/CEbOOlJzCCNYTrWaqN1LkWQKbcqfPHy+17qxthBdpa/oN9T/ku7SjziBM/muvW6KtRdH6FFxbo6+wJqH31ySZydxXKTy79oOGt1aVaFVpblUWlnOVjoLl+AKlHFv12EFVqEBHdCkJ9I76b30QW7KD+RH8uNp6oyUcS7AX0d++gtd+gNm</latexit>⇣12~⌦q (V, T )i⌘(Cq† ⇤(HA)Cq )jj ↵q .
q28.235

Imaginarymode exists

28.240Free energy (eV)28.245

No imaginarymodes within
HA

28.25028.25528.26028.26528.270T = 300 K

28.275
63.0

63.5

64.0

64.5

65.0

65.5

66.0

66.5

Volume (Å3 )Y. Ooba, TT, R. Akashi and S. Tsuneyuki, Phys. Rev. Materials 3, 033601 (2019).

56

Eq. (7) asij

A. Quasiharmonic theory"%free energy#$
Anharmoniccorrection
to
!
1
1
(PT)
!!
1
"
# 1∆q nq +$ % .
! −
1
2
2
(QH)
Fvib
(V, T) =ln 2 sinhq β!ωq (V) .
(2)
β

q

2

(a)

10

(b) 100(8)
53This result is the same as the one reported by Allen [29]. In thisHere,βalways=term
1/kTusetheBoltzmann
constantk, generalized
andEq.
ωq(6)(V)as andstudy,
weEq.(5)(7)
becauseitis morethird
inwithEq.
canbe removedby theoryusing
Self-consistentphonon(SCP)
the phonon tofrequency
volumeunstable
V obtained(imaginary)
within the harevenisapplicablesystemsat where
phonon
"
%
$ is
monic
approximation
eﬀect
1 !While the #anharmonic
1
modes
exist
within the (HA).
HA.
(SCP)
(V, T) =2 sinhdependence
β!Ωq (V,ofT)ωq (V),partiallyFvib
incorporated
via the ln
volume
(b)
(b)
β 100
2
C. Improved self-consistentphonontheory100

CTE (ppm/K)

The
theoryis the standard
Fvibquasiharmonic(V, T) = (QH) ln
2 sinhβ!Ωqapproxima(V, T)
β q the vibrational2free energy is givention of Fvib (V, T), in whichas Quasiharmonic (QH) theory
)
*

0

90

Bulk modulus (GPa)

q,i j

80

7770

−5

60

Bulk
Bulk modulus
modulus (GPa)(GPa)

q the intrinsic anharmonic eﬀects
the QHA completely neglects
−10
50
&
'
!
which are responsible formakingthetemperaturedependence
1
2† (HA)
Expt. (Greve et al.)
−
Ω
(V,
T)
−
(C
Λ
C
)
α
.
(7)
q
j
j
q
q
of
phononfrequencies.
Nevertheless,the
QH
theoryturns
qq90
The SCP theory accounts
4 q90 for the quartic anharmonicity
−15
40
0
250
500
750
1000
0
out to be a good approximation at temperatures far below thenonperturbatively but neglects the eﬀect of the cubic anharTemperature (K)melting point and has been80employed to predict the thermal
80(SCP)
monicity.
The
ISC
theoryaccounts
theadditional
expansivity
and
phase
boundaryof various
materials
basedon (V, T)
Hence,
evaluating
the [27]
anharmonic
freeforenergyFvib
SCP (ISC)
three-phonon
term
perturbatively
astheFIG. 6. (a) CTESCP
and (LDA)(b)
isothermal bulkISC
modulus
DFT.Improved
the temperaturereaches
melting
pointtoorthe
theinterpo(LDA)
inWhenthe thermodynamic
limittheoryisfeasible
owing
SCP
(LDA)
ISC
(LDA)calculated w
70 firstthe
et toal [6] and Morelock
et al. [49].
70anharmonic,
QHA
reliSCP(5)structure
ISC
lationis strongly
technique.
The
term
inbecomes
Eq. (7) lesscorresponds
SCP (PBE)(PBE)
ISC(PBE)
(PBE)
(ISC)
(SCP)
(B)
able.theMoreover,theT)
QHA
not validforphononSCF
ISC
contribution[Eq.
(2)]T)
but
the where
harmonic
frequency
SCF (PBEsol)(PBEsol)
FQH(V,= Fisvib(V,+casesF(V, T),(9)
ISC(PBEsol)(PBEsol)
vib
vib
modes
withinHA as will The
be discussed
60 theis becomereplacedunstablewith theSCP
frequency.
second interm is a60
D. Temperature dependence of soft mode frequenc
Sec. correction
IVnecessary to satisfy the correct thermodynamic re(B)B.

where Flationship(V, T) istheHelmholtz
freeenergyfrom thebubblevib
S=−dF/dT
[29],
where
the
entropyisgiven asboundary
Improved
prediction
of
phase
50
50
(
More
diagram,S associatedwithln(ntheq +1)−ncubicqanharmonicity,
given
as
=
k
[(n
+1)
ln
n
]
with
n
=
n(Ω
).
When
Finally, we discuss the influence of quartic and cubqqqqΩaccurate)]
q al.)
Expt. (Greve(Greve etetal.)q Expt.
Expt.
(Morelocketal.)Expt. (Morelock et al.)the frequency of the R4+ soft mode.
to wecomputefree-energy
the anharmonic eﬀectAble
is small,can assumeΩ ≈ ω +harmonicities∆of high-Tonphases

q
q
q
40
tion.
40
B.≈ Self-consistentphonon00
750and C
1000
0 we
250 theory
500perturbation
750
1000
δ
.
Then,
obtainthe500limitureof 7 shows the temperature dependence of the sq
00
750
1000
0
250
750
1000
q
,i
j
i
j
rmal
perature(K)
Temperature (K)phonon frequency of the R4+ soft mode calculated wi
Eq.
(7)
as
erature(K)
Temperature
(K)
, see
Y. theory
Ooba,isTT,
Akashiand S. Tsuneyuki,
Phys.
Rev.functional.
MaterialsThe
3, 033601(2019).
PBEsolresults shownin Fig. 57
7(a) wa
The SCPone R.
of themost
successful
approaches
"
%
#
$
!

Anharmonic effects in HxS
H3S (Im-3m) @ 250 GPa

H2S (P-1) @ 140 GPa

300

250

200

200

Frequency (meV)

Frequency (meV)

250

150

100

SSCHA
Phys. Rev. Lett. 114,
157004 (2015)

100

50

50

0
Γ

150

Harmonic
SCPH

H

N

P

Γ

N

0
Γ

Harmonic
SCPH

X

U

R

Γ

Z

SCP (our method)
Tc based on the Migdal-Eliashberg theory*

Tc = 202 K (Harmonic)

Tc = 44 K(Harmonic)

Tc = 181 K (Anharmonic)

Tc = 34 K(Anharmonic)

Expt.

Tc ⇠ 160 K*W. Sano, T. Koretsune, TT, R. Akashi, and R. Arita, PRB 93, 094525 (2016).

58

Anharmonic effects in LaH10

2Superconductivity at ~250 K is observed under high pressure (~140—220 GPa)
A.P. Drozdov et al., Nature 569, 528 (2019).

For LaH10, harmonic approximation breaks down.
The SCP theory gives stable phonons at T = 0 K.

Fm-3m

R-3m

250

300
Harmonic (2x2x2)
250
Anharmonic (2x2x2, T = 0 K)

200

200

Frequency (meV)

Frequency (meV)

300

150
100
50
0
50

300

Anharm

150
100
50
0

X U|K

L

W

X

50

T

Harmonic (2x2x2)
250

HarmonAnharmonic (2x2x2, T = 0 K)

H2|H0

L

S2|S0 F59

Nuclear quantum effects in LaH10
2

I. Errea, TT et al., arXiv:1907.11916

Comparison of energies (eV / LaH10)

Figure 1. Quantum eﬀects stabilize the symmetric F m-3m p
pressure for diﬀerent structures of LaH10 calculated neglecting the zer
from V (R), neglecting quantum eﬀects on it. The crystal structure of th
of a Born-Oppenheimer energy surface V (R) exemplifying the presence
represents the positions of atoms treated classically as simple points. B
surface including quantum eﬀects. R represents the quantum centroid
functions, i.e., the average atomic positions. All phases collapse to a si

Structure
E0 + PV
Ezp-vib
Total
by an independent group that measured a T of 250 K
from 137 to 218 GPa in a structure with fcc arrangement-6350.471
2.266-6348.205
1) Fm-3mof the La atoms and suggested a LaH stoichiometry [3].
Even if it is tempting to assign the record supercon-6350.517
2.423-6348.094
2) R-3m
gure 1. Quantum eﬀects stabilize the symmetric F m-3m phase of LaH . Top panel: Enthalpy
as function
ductivity
to theof F m-3m phase predicted previously [2, 3],ssure for diﬀerent structures of LaH calculated neglecting the zero-point energy. The pressure in thethere
figure isis acalculated
clear problem: the F m-3m structure is predicted-0.046
0.157
Δ=X
2-X
1 on it. The crystal structurem V (R), neglecting
quantum
eﬀectsof the diﬀerent
phases found are shown. 0.111
Bottom left: Sketch
to be dynamically unstable in the whole pressure range

tions
F m-3
for th
10
We
struc
10
10
of-th
The
Figure 1. theQuantumstabilize
thesymmetric
F m-3mRphase of LaH10 . Top panel: Enthalpy as functiona Born-Oppenheimer energy surface V (R) exemplifying
presence of eﬀects
many local
minima for
many
distorted structures.
pressurefor diﬀerent
structures
LaHofcalculated
neglecting
theenergy. The
pressure
in thatthe figure
where
a 250
Kzero-point
Tc was observed.
This
implies
this is calcula
not i
10the
resents the positions of atoms treated classicallyas simple
points. Bottom
right: ofsketchconfigurationalE(R)energyface including quantum eﬀects. R represents
the V
quantum
centroid positions,determinethe centerof structure
theisionic
from
(R),
neglecting
quantumwhicheﬀectson it. The
crystalthe diﬀerent
phasesfound are shown. energy
Bottom left:
Ske
phase
not waveaofminimumof theBorn-Oppenheimer
Born
Fm-3m
structure
gives
the
60
nctions, i.e., the average atomic positions. All
collapse to a singleenergy
phase, surface
the highly
symmetric
F m-3m. the presence of many local minima for many distorted structuresof phases
a Born-Oppenheimer
V (R)exemplifying
surface. Consequently, no Tc has been estimated for thispositc

ALAMODE (Anharmonic Lattice Model)
Software for anharmonicity and thermal transport
Open source (MIT license)
Can be combined with any DFT codes and classical MD codes
Inputs data: crystal structure and training data set (displacement & force)

Main features(1) Extraction of force constants by linear
regression (fourth- and higher-order terms
are also supported)(2) Phonon calculation (harmonic & selfconsistent phonon)(3) Thermal conductivity calculation
http://alamode.readthedocs.io/

Please visit the document page for details
61Summary
Perturbative treatment of lattice anharmonicity
Non-perturbative treatment of fourth-order anharmonicity based on
the self-consistent phonon theory
Able to compute phonon lifetimes and thermodynamic functions ofseverely anharmonic materials and high-temperature phases.
Compressive sensing approach to estimate force constants
100x speed up is expected!
Software development
Features presented today are available in ALAMODE software.
Review paper:
TT and S. Tsuneyuki, J. Phys. Soc. Jpn. 87, 041015 (2018)
ALAMODE hands-on materials:
https://github.com/ttadano/alamode/tree/gh-pages/files

62Anharmonic self-energy
Second-order diagrams

First-order diagram“loop”
Real“tadpole”
Real“bubble”
Complex

Higher-order diagrams
Phonon scattering by the three-phonon process
Linewidth…
Lifetime(in high-T region)64Self-consistent phonon theory from Dyson eq.
The self-consistent phonon (SCP) theory is a non-perturbative approach tocompute anharmonic phonon frequency.
In this study, we estimate

from the pole of the Green’s function.

65

Iterative update of phonon frequency
SCP eq.

In this work, we only consider the loop diagram, which is
initial value-independent.

Efficient implementation

TT and S. Tsuneyuki,
Phys. Rev. B 92, 054301 (2015).
interpolation

diagonalization

converged

SCP solution

66

Loop diagram✓ Correct treatment of zero-point vibration✓ Easy to check the finite-size effect by changing the q1-grid density✓ Renormalized frequencies at various T can be calculated efficiently✓ (Probably) systematically improvable✗ Need to calculate fourth-order force constants

67Temperature-dependent effective potential (TDEP)
A MD-based approach to compute effective force constants.
O. Hellman, I. A. Abrikosov, and S. I. Simak, Phys. Rev. B 84, 180301 (2011).

Fi =

X

⇥ij (T )uj (T )j<latexit sha1_base64="SRc8ixeodrbejbTUtiekke7vzd0=">AAAC1nichVG/SxxREP5ck2jMD0/TBGyWHAYtcsyGgCIIQkAs/XGnB54su5unPt1f7L490OXsRLDRzsJKwUKs02phI6RO4Z8QUhpIkyKzewtBRTPL7pv5Zr7Zb97YoStjRXTdoXU+efqsq/t5z4uXr173lvr65+MgiRxRcwI3iOq2FQtX+qKmpHJFPYyE5dmuWLDXP2f5haaIYhn4VbURiiXPWvHlsnQsxZBZGp40U9nSx/UPjTjxzHSt1aiuCmUxutbSh6rDepKB7JilMlUoN/2+YxROGYVNB6VvaOALAjhI4EHAh2LfhYWYn0UYIISMLSFlLGJP5nmBFnqYm3CV4AqL0XX+rnC0WKA+x1nPOGc7/BeX34iZOgbpO53SDV3RGf2gPw/2SvMemZYNPu02V4Rm7+7bud//ZXl8Kqz+Yz2qWWEZo7lWydrDHMmmcNr85ubBzdzY7GD6no7pJ+s/omu65An85i/nZEbMHj6ix+a5LY4F93z45rI6L8vzIo27a7vvzH+sGFQxZj6VJ0aLlXZjAO8wxHsbwQSmMI0a99/DV5zjQqtrW9q2ttMu1ToKzhvcMm3/Lx0SqFg=</latexit>⇥ = arg min kF TDEP⇥

<latexit sha1_base64="LsjvdLMPE89LdEfTrUXC2v+rzPA=">AAADDHichVHPS9xQEB5jW6394aoXwUvoonjpMpFCRRAsatGDsOquCkbDS3zuPsxLQvJ2wab5B3rotQdPCh7Ea2+lp9JS6LkH/4TiUYsXD06yAbWLdkLyZr5vvsnMGztwRaQQTzq0zgcPH3V1P+558vTZ895CX/9K5DdCh1cd3/XDNZtF3BUeryqhXL4WhJxJ2+Wr9s50yq82eRgJ36uo3YBvSFbzxLZwmCLIKlimLWOzUueKJfqkbrKwpptSeFZ8g0h0872exm+TzdiUTNVDGVdmZsvEvGwj3swvzGSSzXgssQpFLGFmertj5E4Rciv7hV9gwhb44EADJHDwQJHvAoOInnUwACEgbANiwkLyRMZzSKCHtA3K4pTBCN2hb42i9Rz1KE5rRpnaob+49Iak1GEYf+MRnuFPPMY/eHlnrTirkfayS6fd0vLA6v0wuHzxX5WkU0H9WnVvzwq2YTzrVVDvQYakUzgtffPdp7PliaXheAQP8JT638cT/EYTeM1z53CRL+3d049NczOKOdW8++bSPJnytEjj37W1OytjJQNLxuKr4tR4vtJuGIIXMEp7ew1TMAdlqFL973AK5/BX+6h91r5oX1upWkeuGYBbpv24AgQev68=</latexit>

F AIMD k2

Fit displacement-force data sets sampled by AIMD at finite
temperature with effective harmonic force• Full anharmonicity included• Neglect the effect of zero-point motion•

The TDEP method can be extended to include third-order terms:
Fi =X

˜ ij (T )uj (T )⇥j<latexit sha1_base64="ykYFt2g0uApIa8pqRfbPlJc8GLM=">AAADHXichVFLS91QEJ4b67Nar3YjdBO8KLrwMhFBEQpCobj0ca8KRsJJPFePNy+Skws25A/0D3ThQip0UfwZUhAs3bXFnyBdKrhx4ZzcSKmvTkjmm2/mm8w5Y4euiCXieUnreNHZ1d3T2/eyf+DVYHloeC0OksjhdSdwg2jDZjF3hc/rUkiXb4QRZ57t8nW7+U7l11s8ikXg1+R+yLc8tuOLhnCYJMoqh++tVGT6W33KjBPPSvcyUwp3m6dmbZdLllF2L9MnapN6opIKTOlmI2JOamTpdFaoWPMRXVOV38nI57FVrmAVc9MfAqMAFShsKSifgQnbEIADCXjAwQdJ2AUGMT2bYABCSNwWpMRFhESe55BBH2kTquJUwYht0neHos2C9SlWPeNc7dBfXHojUuowhj/xK17iKR7jBd482SvNe6hZ9snbbS0PrcGPI6vX/1V55CXs/lU9O7OEBszlswqaPcwZdQqnrW99+HS5Or8ylo7jEf6h+T/jOZ7QCfzWlfNlma8cPDOPTedmFHPq+fTNqTpP5WmRxv21PQRr01UDq8byTGVhrlhpD7yBUZigvc3CAizCEtSp/y+4KXWWurRD7Zv2XfvRLtVKheY1/GPa71vf/sQl</latexit>

1X˜⇥ijk (T )uj (T )uk (T )
2
jak

˜ ij
⇥ij 6= ⇥<latexit sha1_base64="74AlIz+ag5B/fNdNK+vy0W7pvqs=">AAAC0XichVE9SxxhEH5cP6OJnqYJpDlyUVIdsxJQrASblOp5Knhy7K6jvnG/svvewbkciKX5ASmsDKQIadJqayNYW/gTJKWCjUVm9xaCimaW3XfmmXlmn3nHDl0Va6LLLqO7p7evf+DF4NDLV8MjhdGx5ThoRA5XncANolXbitlVPle10i6vhhFbnu3yir0zl+ZXmhzFKvCXdCvkdc/a8tWmciwtUL3wvra0zdqqJ+pzu1jz+UuxppW7wUkHb2eJeqFEZcqs+Ngxc6eE3OaDwjlq2EAABw14YPjQ4ruwEMuzBhOEULB1JIJF4qksz2hjULgNqWKpsATdke+WRGs56kuc9owztiN/ceWNhFnEOF3QT7qmM/pFV3T3ZK8k65Fqaclpd7gc1kcO3lRu/8vy5NTY/sd6VrPGJqYzrUq0hxmSTuF0+M3db9eVmcXxZIK+0x/Rf0SXdCoT+M0b58cCLx4+o8eWuS2JWXo+fXNpnZfmZZHmw7U9dpYnyyaVzYWPpdnpfKUDeIt3+CB7m8IsPmEeVel/gN84xolRMVrGnrHfKTW6cs5r3DPj61/jm6gG</latexit>

68