DFT2024@RIKEN, Feb. 22nd, 2024

Anharmonic phonon theory in condensed matter physics
with applications based on DFT
Terumasa Tadano
National Institute for Materials Science, Tsukuba, Japan

E-mail: TADANO.Terumasa@nims.go.jp

Outline
1. Introduction
2. History of anharmonic phonon theory
3. Perturbative expansion method and applications to thermal transport
4. Self-consistent phonon theory and applications to phase transition
5. Summary & outlook

What is phonon?

Quantum description of lattice vibration

1D chainInteraction
displacement

Eq. of motion

Phonon dispersion (1D)Why/when is phonon important?

Thermal conductivity

Band-gap renormalization

Superconductivity

H3S

LaH10R. Berman. Thermal conduction in Solids.

qqʹ
qʺ

Phonon-phonon catteringkq

PRL 105, 265501 (2010)- thermal expansionkʹ
Phonon-electron cattering

Other phenomena include

Electron-phonon interaction- structural phase transition

Non-interacting picture
Hamiltonian of non-interacting phonon (harmonic approximation)✗ Phonon is non-interacting→ infinite lifetime→ infinite thermal conductivity✗ Phonon is temperature-independent→ no displacive structural phase transition✗ Phonon is volume-independent→ no thermal expansion

Non-interacting picture
Hamiltonian of non-interacting phonon (quasi-harmonic approximation)✗ Phonon is non-interacting→ infinite lifetime→ infinite thermal conductivity✗ Phonon is temperature-independent→ no displacive structural phase transition✗ Phonon is volume-dependent→ finite thermal expansion

Many-body problem of phonons
Phonon Hamiltonian with interaction terms (anharmonicity)✓ Phonon is interacting via anharmonicity→ finite lifetime→ finite thermal conductivity✓ Phonon is volume-dependent→ finite thermal expansion✓ Phonon is temperature-dependent→ explains structural phase transition

How to solve the many-body problem?

Perturbative expansion•

perturbation

Dyson equation

Free energy expansion (cumulant)

Variational approach (self-consistent phonon)•

Find a good one-body approximation
based on the free-energy minimization•

Monte Carlo or Molecular dynamics (if classical)

History of many-body phonon theory (and others)
1910

1920

1930

1927
BornOppenheimer

1912
Debye
Born-vontheory
Karman theory
Birth of “latticedynamics”

1940

1950

1980

1990

Hartree-Fock

1929
PeierlsBoltzmann
theoryRenaissance of anharmonic phonon theory
1950

1960

1955
Self-consistentphonon (SCP)
theory(Hooton)

1957
Kubo
formula

1962
Phononscattering(Maradudin-Fein)

1970

1964
DFT(HohenbergKohn)

CASTEP(first release)

1981
PerdewZunger(LDA)1965
LDA(KohnSham)

Anharmonic phonon calculation based on DFT
2000

1990

1993
VASP

1992
PerdewWang(GGA)1996
PBE(GGA)

2001
Quantum
ESPRESSO(first release)

2010

2008
SCAILDFirst ab initio
2007
DFT-based implementationof SCP
thermal(Souvatzis)
transport(Broido)

2020

2015
2014
FC-SCPSSCHA (Tadano)(Errea)2017
2011
TDEP(Hellman)

DFT-based4-phscattering(Ruan)

2019
2022
Wignerthermal Structure opt.
transport within SCP(Masuki)(Simoncelli)2023
A beyond
2022
SCP
Quasiparticle (Xiao)
SCP(Tadano)Perturbation theory

perturbation 𝐻 !
Dressed phonon Green’s function

Dyson eq.

“tadpole”“loop”“bubble”

3-phononinteraction(cubic anharmonicity)

4-phononinteraction(quartic anharmonicity)Anharmonic self-energies“tadpole”“loop”“bubble”

Complex“4ph”

Complexeigen-energy shift
linewidth (inverse scattering rate)

3-phononinteraction(cubic anharmonicity)

4-phononinteraction(quartic anharmonicity)

Bubble self-energy
Phonon scattering by the three-phonon processd-GeGe46

Ba8Ga16Ge30

Linewidth10-fold
reduction

Lifetime

10x larger phonon-phonon scattering rate of heat-carrying acousticphonons induced by low-energy rattling guest modes in BGG
Maradudin and Fein, Phys. Rev. 128, 2589 (1962).

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

Prediction of thermal conductivity

Lines : calculation
Points: experiment

Boltzmann transport theory

SiMg2Si
PbTe
Bi2Te3

Lifetime

SrTiO3

Ba8Ga16Ge30

Good agreement between theory and experiments

useful for computational discovery of novel
thermal management materials

4-phonon scatteringFirst theoretical prediction

Experimental validation

Phys. Rev. Lett. 111, 025901 (2013)

Science 361,575-578 (2018)

Boltzmann theory
only with 3ph process

Boltzmann theory
with 3ph & 4phprocesses

Inclusion of 4ph scattering lead to better agreement with theexperimental thermal conductivity
However, the same approach does not necessarily improve theprediction by Boltzmann+3ph in general cases.

Limitation of the perturbation method̂•"Series expansion may not converge if the perturbation term 𝐻 is reasonably large.

This can happen more likely for strongly anharmonic materialŝ

•The non-interacting Hamiltonian 𝐻# is ill-defined in many cases (high-symmetry phases of solids).

Phonon(Harmonic)
Frequency is pure imaginary(“imaginary phonon”)

Need to construct a better one-body starting point

Anharmonic renormalization solves imaginary phonon problem
Effective harmonic force constant at finite temperature

SrTiO3 tetragonal (T < 105K)

thermal/quantum fluctuation
SrTiO3 cubic (T > 105 K)

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 2011

Self-consistent phonon (SC1) theory
N. R. Werthamer, Phys. Rev. B 1, 572 (1970)

̂𝐻 : Exact Hamiltonian of the interacting phonon system with full anharmonicity.
Density operator:
Exact free energy:

infeasible to calculate: Effective one-body Hamiltonian of the interacting phonon system
Density operator:
Approximated free energy:
Feynmann–Gibbs–Bogoliubov inequality:

computation isfeasibleSC1 free energy and SC1 frequency
Feynmann–Gibbs–Bogoliubov inequality:
SC1 vibrational free energy
SC1 frequency
anharmonic
force constants

even-term anharmonicity

mean-squaredisplacement

A diagrammatic view of SC1 theory
The self-consistent phonon (SC1) equation is obtained via“Hatree-Fock” theory for phonons
A diagrammatic representation of the SC1 theory:
Self-energy in self-consistent phonon theory

Self-energy in perturbation theory

Tadpole (T)Loop (L)Force constant calculation: A brute-force approach

2 patterns

4 patterns

SiliconFC6FC5NaCl
FC4FC3FC2FC6FC5FC4FC3FC2

Force constants from sparse modelingparametersbasis

M: The number of independent parameters

Shrinkage

p=1: LASSO

p=2: ridge regression

Benchmark: Efficiency of sparse modelingResult for cubic CsPbBr3✓ Adaptive LASSO gives more sparse solution✓ 50 structure snapshots is enough to reach convergenceof Fvib within 1 meV/atom error

3294 coefs. in total

Brute force
LASSO

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)

First-principles implementation of SCP (SC1)
How to update Φ (effective 2nd-order force constants )?
Stochastic methods-SSCHA

http://sscha.eu-QSCAILD-

HiPhivehttps://github.com/vanroeke/qscaild
https://hiphive.materialsmodeling.org-Phonopy

https://phonopy.github.io/

Force constant based approach (FC-SCP)-ALAMODE

https://alamode.readthedocs.io/From, derive the self-consistentequation analytically. The Taylor series is truncated.

Stochastically displace atoms in a supercell and update
Φ so as to minimize

Pros. Full anharmonicity (at the mean-field level)• No explicit computation of anharmonic
force constants

Pros.

•Cons. • More expensive computationallyMore efficient particularly for scanningtemperature• Faster convergence•

Cons.
Omit higher-order anharmonicity• Requires anharmonic force constants as input•

Consistency between two approaches

Stochastic method (SSCHA)

Force-constant based approach (FC-SCP)

I. Errea et al., Phys. Rev. B 89, 064302 (2014).

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

H3S

H3 S
Phys. Rev. Lett. 114, 157004 (2015)

H3 S
Phys. Rev. B. 93, 094525 (2016).

Another comparison on anharmonic solidSoft phonon frequency

High-energy phonon frequency

CsPbBr3

MD EHM = TDEP法•

Does not satisfy the variational principles. Looks heuristic.

E. Fransson, P. Rosander, F. Eriksson, J. M. Rahm, TT, and P. Erhart, Commun. Phys. 6, 1 (2023).

Application to el-ph problem: band-gap renormalization
HarmoniccubicSrTiO3Anharmonic

Allen–Heine–Cardona (AHC) theory

Finite difference method using non-diagonal supercell
Monserrat, J Phys Condens Matter 30, 083001 (2018).
Y. Wu, TT et al., J. Phys. Chem. Lett. 11, 2518 (2020).

History of superconducting Tc record
Somayazulu et al.,PRL 122, 027001 (2019).

p =180 GPa!

https://kittaka.r.chuo-u.ac.jp/contents/others/tc-history/index.html

Distortion and phonon instability in 𝐹𝑚3𝑚-LaH10

𝐹𝑚3𝑚-LaH10 distorts into alower-symmetry structure
below ~250 GPa when atoms
are treated as classicalparticles.

Stabilization of 𝐹𝑚3𝑚-LaH10 by anharmonic renormalization
Classical particle(at absolute rest)
No anharmonic renormalization

Quantum fluctuation(> 10% of interatomic distance)

Large anharmonic renormalization

Nature 578, 66-69 (2020)

Tc calculated with renormalized phonon

Marques et al., PRB 72, 024545 (2005).
Sanna et al., JPSJ 87, 041012 (2018).
Sanna et al., PRL 125, 057001 (2020).

Migdal–Eliashberg Theory

Superconducting Tc (K)DFT for superconductors (SCDFT)

Pressure (GPa)Sano et al., PRB 93, 094525 (2016).

Nature 578, 66-69 (2020)

Phonon linewidth from first principles
SC1 is still one-body (infinite lifetime)harmonic

SC1 + Bubble gives finite line width

SC1“bubble” self-energy

4th-order anharmonicity

3rd-order anharmonicity

SrTiO3: phonon linewidth & thermal conductivity
Phonon spectral function

Green line: SCP result@300K

Boltzmann transport theory

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

Quasiparticle approximation

Quasiparticle (QP) approx.
SrTiO3•Full-frequency dependence•Better one-body approximation•(often) inconvenient for el-ph calculation•

Dropped frequency dependenceTT and W. Saidi, Phys. Rev. Lett. 129, 185901 (2022)

Quasiparticle approximation

Quasiparticle (QP) approx.
SrTiO33rd-order
anharmonicity
4th-order
anharmonicity•

Full-frequency dependence•

Better one-body approximation•(often) inconvenient for el-ph calculation•

Dropped frequency dependence

TT and W. Saidi, Phys. Rev. Lett. 129, 185901 (2022)

GWA and SCP: A comparison

starting(unperturbed)
G0Interaction

GWASCP

Kohn–Sham (DFT)SC1(Hartree–Fock for phonons)

Screened Coulomb W

Renormalized
anharmonicity Φ"QSGW, scGW

QSGW-like treatment,
another approach (e.g., PRB 107, 094303 (2023))

Self-energy

QP approx.
(“one-shot”)
Beyondone-shotBut, accurate prediction of Tc is challenging

Quasiparticle (QP) approx.

500 K
CsPbBr3

SC1

QP approx.

SCP theory: Tc = 198 K
QP approx.: Tc =423 K
Expt.: Tc = 403 K

Influence of lattice parameter

Expt. Tc
(403 K)

AFD mode

TT and W. Saidi, Phys. Rev. Lett. 129, 185901 (2022)

Structure optimization at finite temperatures
Self-consistent phonon (SC1) theory
Able to compute renormalized phonons at finite temperatures
Able to compute crystal structures at finite temperatures
R. Masuki, TT et al., Phys. Rev. B 106, 224104 (2022)
R. Masuki, TT et al., Phys. Rev. B 107, 134119 (2023)

Ryota Masuki(Univ. Tokyo)

SC1 free energy and SC1 frequency
Feynmann–Gibbs–Bogoliubov inequality:
SC1 vibrational free energy
SC1 frequency
anharmonicforce constants

even-term anharmonicity

mean-squaredisplacement

Atomic force at finite temperature
SC1 free energy

Free energy surface

Atomic force
DFT force

SC1 correction

displacement

odd-term anharmonicity
See also stochastic implementation (SSCHA): Monacelli et al., Phys Rev B 98, 024106 (2018).

DFT and SC1: A comparison
DFT (0 K)
What tominimize

ForcesStress tensor

Enthalpy

SC1 (finite temperature)
Free energyComputational workflow

basis

parameters
for each temperature

Generate displacement-force dataset using DFT calculation (~ 100 training snapshots)
Estimate force constants using compressive sensing (linear regression)
adaptive LASSO
F. Zhou, W. Nielson, Y. Xia, V. Ozoliņš, PRL 113, 185501 (2014).
TT and Tsuneyuki, PRB 92, 054301 (2015)

[001]

[110]

[111]
initial displacement
works as a small
perturbation thatbreaks the symmetry

𝑃𝑚3𝑚

𝑃4𝑚𝑚

𝐴𝑚𝑚2

Figures from J. Fujioka. et al, Sci Rep 5, 13207 (2015).

𝑅3𝑚

“Renormalization” of force constants
nth-order IFC
in a displaced structurenth-, (n+1)th-, … order IFC
in the reference structure
Normal mode
displacement (variable)

nth-order IFC
in a distorted lattice

Strain tensor (variable)

Cartesian coordinate

R. Masuki, TT et al., Phys. Rev. B 106, 224104 (2022)

How good is “renormalization”?
Harmonic phonon of BaTiO3

Thermal expansion of Si

purple: compute harmonic and anharmonic
IFCs for 14 different lattice constants✓

Good agreement between the renormalization
technique and full DFT calculation.

green: compute harmonic and anharmonic IFCs
for a0 only. The IFCs at different lattice
constants are estimated by “renormalization”.
R. Masuki, TT et al., Phys. Rev. B 106, 224104 (2022)

Application to consecutive transition in BaTiO3•••••

VASP code for DFT (PAW)
PBEsol xc functional
2x2x2 supercell for force constant calculation
Truncate Taylor series at the fourth-order
80 training structures (80 static DFT calculations)
https://alamode.readthedocs.io/

SC1 frequency for cubic BaTiO3Ferroelectric soft mode

Cubic-to-Tetragonal transition•

Start the calculation with small displacements along [001]✓Succeeded in predicting structural phase transition
with a temperature hysteresis✓

Predicted 𝑇! ~ 600 K✓

Overestimate Tc due to slightly larger (+0.01 Å) lattice constantby PBEsol

Expt.: 𝑇! ~ 393 K

Previous theoretical studies:•MC, effective Hamiltonian (no ZPE, LDA): 297 K [1]•MC, effective Hamiltonian (ZPE, LDA) : 265 K [2]•MD, DFT-fitted model potentials (PBEsol): 160–170 K [3][1] Zhong, Vanderbilt, Rabe, PRL1994.
[2] Zhong, Vanderbilt, PRB1996.
[3] Qi, Liu, Grinberg, Rappe, PRB2016.

R. Masuki, TT et al., Phys. Rev. B 106, 224104 (2022)

Successive phase transition•••

Tetra.: small displacements along [001]
Ortho.: small displacements along [110]Rhombo.: small displacements along [111]✓ Successive phase transition can be reproduced with the correct sequence
R. Masuki, TT et al., Phys. Rev. B 106, 224104 (2022)

Lattice parameters and spontaneous polarization✓ Able to compute finite-temperature structural parameters and spontaneous polarization
R. Masuki, TT et al., Phys. Rev. B 106, 224104 (2022)

Finite temperature + finite pressure

Our methodExpt.
Phys. Rev. Lett. 78, 2397 (1997).

Comparison between SC1 and MDMolecular dynamics (MD)SC1•

Lower computational costs•

Finite-temperature structures canbe obtained deterministicallyPros.
•Phonons can be obtainedsimultaneously•

Include zero-point motion•

Approximate the atomic
distribution by Gaussian
assuming displacive transition

Cons.
•Anharmonicity is included at
the mean-field level•

Full anharmonicity•

Able to simulate orderdisorder transition•

High computational cost•(usually) ignores the zeropoint motion

Summary & perspective
Perturbative expansion has been successful in predicting phononphonon scattering rates and thermal conductivity of various solids.
However, it may breakdown for strongly anharmonic solids.
Variational approach (Self-consistent phonon) is a powerful tool forobtaining a better one-body description of phonon. It even enables usto calculate structures at finite temperatures.

More application studies
Transport properties near the phase transition
Higher-order anharmonicity (5th, …)Comparison to MD using machine-learning FF