#  Ab initio Investigation of Thermal Transport in Insulators: Unveiling the
Roles of Phonon Renormalization and Higher-Order Anharmonicity

Soham Mandal  Manish Jain  Centre for Condensed Matter Theory, Department of
Physics, Indian Institute of Science, Bangalore-560012, India  Prabal K. Maiti
<EMAIL_ADDRESS>Centre for Condensed Matter Theory, Department of Physics,
Indian Institute of Science, Bangalore-560012, India

###### Abstract

The occurrence of thermal transport phenomena is widespread, exerting a
pivotal influence on the functionality of diverse electronic and thermo-
electric energy-conversion devices. The traditional first-principles theory
governing the thermal and thermodynamic characteristics of insulators relies
on the perturbative treatment of interatomic potential and ad-hoc displacement
of atoms within supercells. However, the limitations of these approaches for
highly anharmonic and weakly bonded materials, along with discrepancies
arising from not considering explicit finite temperature effects, highlight
the necessity for a well-defined quasiparticle approach to the lattice
vibrations. To address these limitations, we present a comprehensive numerical
framework in this study, designed to compute the thermal and thermodynamic
characteristics of crystalline semiconductors and insulators. The self-
consistent phonon renormalization method we have devised reveals phonons as
quasiparticles, diverging from their conventional characterization as bare
normal modes of lattice vibration. The extension of the renormalization impact
to interatomic force constants (IFCs) of third and fourth orders is also
integrated and demonstrated. For the comprehensive physical insights, we
employed an iterative solution of the Peierls-Boltzmann transport equation
(PBTE) to determine thermal conductivity and carry out Helmholtz free energy
calculations, encompassing anharmonicity effects up to the fourth order. In
this study, we utilize our numerical framework to showcase its applicability
through an examination of phonon dispersion, phonon linewidth, anharmonic
phonon scattering, and temperature-dependent lattice thermal conductivity in
both highly anharmonic materials (NaCl and AgI) and weakly anharmonic
materials (cBN and 3C-SiC). Our study finds that ignoring higher-order
(4-phonon) phonon scattering processes is not viable for materials with strong
anharmonicity in their interatomic potential. Meanwhile, renormalization
demonstrates negligible impact on materials characterized by weak
anharmonicity. The investigation also explores the effect of fourth-order
Helmholtz free energy and pressure-dependent thermal conductivity of NaCl
crystal. The theoretical framework developed in this work will help us to
understand the physical insight of phonon and phonon-driven thermal and
thermodynamic properties of materials and offer valuable guidance for the
strategic development of efficient thermal management techniques.

##  I  INTRODUCTION

Lattice vibrations are the primary heat carriers in crystalline insulators and
semiconductors. The quasi-particles of quantized lattice vibration is known as
a phonon. Understanding phonon and phonon-driven properties in materials are
of immense importance in the domain of material science, solid-state physics,
and chemistry. This is due to the evergrowing technological application of
such materials. Phonons play a role in several physical phenomena, including
thermal transport, thermal expansion of solids, the stability of crystalline
solids, and structural phase transitions, either directly or indirectly.
Phonons and several phonon-assisted physical phenomena have been theoretically
analyzed and well-established for more than a century with their experimental
validations. Rapid increase of high performance computational power has
enhanced our ability to use first-principle computation to predict phonon-
driven physical phenomena through numerical techniques. These traditional ab
initio schemes range from simple models taking into account only the harmonic
approximation [1, 2] to predict phonon dispersions, the density of states, and
heat capacity to more sophisticated approaches to include anharmonicity in the
interatomic potential to predict thermal expansion and thermal conductivity of
materials[3, 4]. These harmonic and quasi-harmonic approximations are
excellent in capturing the aforementioned phonon-driven properties of weakly
anharmonic strong covalently bonded materials, like silicon and diamond[5].

In the case of materials with a significant degree of anharmonicity in their
interatomic potentials, such as alkali halides and heavy-metal chalcogenides,
simple models of calculating the normal modes of lattice vibration are
inadequate. Materials with extreme anharmonicities, such as solid helium,
typically have bare phonon frequencies having imaginary values[6].
Renormalization makes these squared frequencies positive as anharmonicity is
sufficiently large. Furthermore, at a higher temperature (depending on the
Debye temperature of a specific material), conventional ab initio approaches
fail to comply with experimental measurements. Therefore, a well-defined and
robust quasi-particle methodology is crucial to cover weakly bonded, highly
anharmonic materials and to fit measured data over the entire experimental
temperature ranges.

The interatomic force constants (IFCs) calculation in traditional ab initio
approaches involves the density-functional perturbation (DFPT)[7, 8] method or
the finite-difference supercell[9, 10, 11, 3] method. The determination of
IFCs in DFPT relies on analyzing the analytical derivatives of the potential
energy. It can efficiently determine harmonic IFCs and related properties of
materials. While the utilization of DFPT has been reported for calculating
third-order anharmonic IFCs[12], the computational cost of determining higher-
order anharmonic IFCs is generally prohibitive. Furthermore, DFPT does not
explicitly incorporate finite temperature effects. The finite-difference
supercell method approximates the derivative of the potential energy by
employing finite differences, and by calculating the forces on atoms within
the perturbed supercell at 0K, the IFCs can be determined. The perturbation is
accomplished by ad hoc displacement of atoms within the supercell. To compute
higher-order anharmonic IFCs (beyond 3rd order), force calculation for a
considerable number of perturbed snapshots is necessary for the finite-
difference supercell method[13]. Recent studies[14, 15, 16, 17] have explored
the incorporation of fourth-order IFCs and four-phonon scattering processes
across various materials to reduce discrepancies between theoretical
predictions and experimental observations of the thermal transport properties.
Classical force field (FF) based molecular dynamics-based approaches for
thermal conductivity calculation, using Green-Kubo response function-based
formalism[18, 19] or the non-equilibrium heating-cooling method[20] can take
into account finite temperature effects and the inherent anharmonicity. In our
current work, we will employ the temperature-dependent effective potential
(TDEP)[21, 22] method, which is a computationally efficient framework capable
of addressing the aforementioned challenges. TDEP has been used in recent
studies to determine the phonons and associated physical properties of
materials[23, 24, 25].

In this work, we have developed a unified numerical scheme to calculate IFCs
up to the fourth order utilizing the temperature-dependent force-displacement
data from ab initio molecular dynamics (AIMD). We have also developed a
Thermal Stochastic Snapshot (TSS) technique to avoid the computational
overhead of AIMD. Our current study also involves a numerical approach for the
self-consistent phonon renormalization and calculation of thermal transport
properties of the insulating crystalline solids. In this renormalization
method, bare second-order IFCs are renormalized only up to first-order using
the statistical operator renormalization technique. This renormalization
process can also be realized through the one-phonon propagator approach in the
many-body theory[26]. In the results section, we will demonstrate the phonon-
driven physical properties of materials with pronounced anharmonicity (NaCl
and AgI) and those with weak anharmonicity (cBN and 3C-SiC) in their
interatomic potentials. This work thoroughly investigates the impact of self-
consistent renormalization of harmonic and anharmonic IFCs and higher-order
phonon scattering processes on thermal transport properties. Finally, our
investigation will examine the influence of fourth-order Helmholtz free energy
on the lattice constant and pressure-dependent thermal conductivity of
crystalline NaCl.

##  II  METHODS

The potential energy of a crystal due to the interaction among the constituent
atoms, $U$, can be expanded in the displacements from the initial
configuration, assuming $U$ is an analytic function of the positions of the
atoms, as:

$\displaystyle
U=U_{0}+\sum\limits_{{N_{1}\mu}}\sum\limits_{\alpha}\Phi_{1}^{\alpha}({N_{1}\mu})u_{\alpha}(N_{1}\mu)+\sum\limits_{{N_{1}\mu},{N_{2}\nu}}\sum\limits_{\alpha,\beta}\frac{1}{2!}\Phi_{2}^{\alpha\beta}({N_{1}\mu},{N_{2}\nu})u_{\alpha}(N_{1}\mu)u_{\beta}(N_{2}\nu)+$
(1)
$\displaystyle\sum\limits_{{N_{1}\mu},{N_{2}\nu},{N_{3}\pi}}\sum\limits_{\alpha,\beta,\gamma}{\frac{1}{3!}}\Phi_{3}^{\alpha\beta\gamma}({N_{1}\mu},{N_{2}\nu},{N_{3}\pi})u_{\alpha}(N_{1}\mu)u_{\beta}(N_{2}\nu)u_{\gamma}(N_{3}\pi)+$
$\displaystyle\sum\limits_{{N_{1}\mu},{N_{2}\nu},{N_{3}\pi},{N_{4}\rho}}\sum\limits_{\alpha,\beta,\gamma,\delta}\frac{1}{4!}\Phi_{4}^{\alpha\beta\gamma\delta}({N_{1}\mu},{N_{2}\nu},{N_{3}\pi},{N_{4}\rho})u_{\alpha}(N_{1}\mu)u_{\beta}(N_{2}\nu)u_{\gamma}(N_{3}\pi)u_{\delta}(N_{4}\rho)+\ldots$

Here $U_{0}$ is the potential energy at the equilibrium. $\Phi_{1}$,
$\Phi_{2}$, $\Phi_{3}$, and $\Phi_{4}$ are first-, second-, third-, and
fourth-order interatomic force constants (IFCs), respectively. $N_{1}$,
$N_{2}$, $N_{3}$, and $N_{4}$ denote the lattice points of the crystal. $\mu$,
$\nu$, $\pi$, and $\rho$ run over the basis atoms; $\alpha$, $\beta$,
$\gamma$, and $\delta$ are cartesian indices. The first-order derivative of
crystal potential $U$ gives the force on atoms:

$\displaystyle F^{\alpha}({N_{1}\mu})=-\frac{\partial U}{\partial
u_{\alpha}(N_{1}\mu)}$ (2)

Material | Order of IFC | Nearest neighbour | Cutoff distance($\AA$) | Total IFCs | Independent IFCs
---|---|---|---|---|---
| 2nd | 11 | 9.88 | 2790 | 55
NaCl | 3rd | 5 | 6.55 | 175446 | 1240
| 4th | 2 | 4.40 | 1111158 | 2627
| 2nd | 11 | 10.75/6.10/7.43 | 2700 | 74
AgI/cBN/3C-SiC | 3rd | 6 | 8.05/4.56/5.56 | 272214 | 3860
| 4th | 2 | 4.90/2.79/3.39 | 795906 | 3758
Table 1: Number of independent IFCs up to 4th order after applying space group
symmetry, permutation symmetry, and acoustic sum rule of the IFCs.

Our calculation of the interatomic force constants (up to the fourth order)
involved conducting ab initio molecular dynamics (AIMD) simulations at a
designated temperature on a supercell of the crystal under consideration.
Unless explicitly stated otherwise, all MD simulations conducted in this study
utilize a $5\times 5\times 5$ supercell configuration. MD trajectories provide
us with the displacements ($u_{\alpha}(N_{1}\mu)$) and forces
($F^{\alpha}({N_{1}\mu})$) on each atom in Eq. 2 at every time step. Eq. 2 can
be rewritten as a matrix equation, $\mathcal{D}x=f$, with these displacement
and force datasets. Here, the matrix elements of $\mathcal{D}$ are constructed
from the atomic displacements; $f$ is a row of atomic forces; $x$ is the row
of IFCs to calculate. The matrix equation is usually an overdetermined system.
A least-squares technique is used to calculate the IFCs by inverting the
displacement matrix $\mathcal{D}$ through the Moore-Penrose pseudoinverse[27].

###  II.1 Symmetry

In general, the number of IFCs for the $d$-th order is $(3N)^{d}$, where $N$
is the total number of atoms in the crystal. As such, solving the matrix
equations for determining higher order IFCs (beyond fourth order) is
computationally prohibitive. To date even calculations including fourth-order
IFCs have only been attempted for simplest materials like Silicon and NaCl. We
make several approximations to make these calculations tractable. While in
principle, one should include all higher-order terms in Eq. 1, we include the
terms till fourth order. To further reduce the number of IFCs, we use cutoff
distances, as the third-, fourth- and non-polar contribution of the second-
order IFCs fall rapidly with distance. Ewald-Summation technique is used to
separate the long-range part of second-order IFCs caused by dipole-dipole
interactions in polar materials, as explained in a later section. This allows
us to identify a suitable set of IFCs for a specific temperature that can
accurately match the forces obtained from molecular dynamics simulations. We
reduce the computational cost further by exploiting permutation symmetry,
point-group symmetry, and acoustic sum rule (ASR) of the IFCs. The specifics
of these symmetries are discussed in Ref[9]. These symmetries provide
constrained equations among the IFCs, from which we extract a set of
independent IFCs. After determining the independent IFCs, the rest of the IFCs
can be reconstructed as their linear combinations. We reconstruct the whole
set of IFCs for the calculation of material properties. For determining point-
group symmetries, we employ the SPGLIB[28] crystal-symmetry library to find
the transformation matrices of the point group operations. The number of
independent IFCs drastically reduces when we combine these symmetries, as
illustrated in Table 1, making it more computationally feasible to calculate
up to the fourth order.

###  II.2 Long-range correction

The vibration-induced dipole-dipole interaction contribution to the second-
order IFCs diminishes at a rate of $\sim 1/r^{3}$ with distance. It poses a
problem in the case of polar materials, as we cannot set a cutoff for second-
order IFCs and, in principle, must extend them to infinity. Most of the
available packages address this issue by incorporating a non-analytic
component into the Dynamical matrix[2, 4]. In this mixed-space approach[29],
the short-range second-order IFCs are calculated in real space separately. To
incorporate long-range correction, a direction-dependent non-analytical term
(in reciprocal space) is introduced to the Dynamical matrix. This correction
is applicable only for the long wavelength $(q\to 0)$ limit. With a reasonable
level of accuracy, this method illustrates the LO-TO splitting of polar
materials along the high symmetry path. In our case, we used the Ewald
Summation technique, as outlined in Ref.[30, 31], to compute the harmonic IFCs
in reciprocal space from the Born effective charge and dielectric tensor of
the material. We then use an inverse Fourier transform to convert the
reciprocal space IFCs resulting from the dipole-dipole interaction into real-
space IFCs. The long-range contribution to the force is then calculated by
multiplying these newly created real-space IFCs with the displacement dataset
already obtained from MD. We eliminated this long-range force component from
the MD force to fit only the short-range component of the IFCs. Subsequently,
the harmonic IFCs (in reciprocal space) obtained through the Ewald Summation
technique, which includes the ionic contribution, are incorporated into the
Dynamical matrix to accurately determine the phonon dispersion, group
velocity, and density of states. Therefore, the complete dynamical matrix is
defined as:

$\displaystyle
D^{\alpha\beta}(\mu\nu,q)=\frac{1}{\sqrt{M_{\mu}M_{\nu}}}\sum_{N}\Phi_{2(short-
range)}^{\alpha\beta}(0\mu,N\nu)e^{i\mathbf{q}\cdot R(N)}$
$\displaystyle+\>Ewald\>term.$

where $M_{\mu}$ is the mass of the basis atom type $\mu$. The phonon frequency
($\omega_{qs}$) and eigenvector ($\epsilon(qs)$) are obtained by performing
diagonalization of the dynamical matrix.

###  II.3 Thermal Stochastic Snapshots

The method we discussed is based on the Temperature-dependent Effective
Potential (TDEP) method described by Hellman et al. [21] TDEP requires
expensive ab initio molecular dynamics (AIMD) on a supercell at a finite
temperature to map it into the model potential energy given by Eq. 1. AIMD is
cumbersome computationally as it requires a sufficiently long run to achieve
equilibrium at a specific temperature, followed by the selection of well-
spaced snapshots for proper sampling of the Born-Oppenheimer (BO) surface of
the supercell. To circumvent the computation-intensive nature of this method
while still retaining all finite temperature effects, we opt for the thermal
stochastic snapshot (TSS) technique, as described in Ref[32, 33, 34]. In the
TSS scheme, the atoms $(Nv)$ in the supercell are displaced following the
equation:

$\displaystyle
u_{\alpha}(N\nu)=\sqrt{\frac{\hbar}{M_{\nu}N_{0}}}\sum_{qs}\sqrt{\frac{2n_{qs}+1}{\omega_{qs}}}\cos(2\pi\xi_{(1,qs)})$
$\displaystyle\times\sqrt{-\ln(1-\xi_{(2,qs)})}\epsilon_{\alpha}(qs,\nu)e^{iq\cdot
R(N)}$ (4)

Where $\xi_{(1,qs)}$ and $\xi_{(2,qs)}$ are two mode-dependent random numbers
generated from a uniform distribution between 0 and 1, with $(qs)$ referring
to phonon modes. To ensure that the displacement $u_{\alpha}(N\nu)$ is a real
value, the random numbers are restricted such that
$\xi_{(1,qs)}=\xi_{(1,-qs)}$ and $\xi_{(2,qs)}=\xi_{(2,-qs)}$. $\omega_{qs}$
and $\epsilon_{\alpha}(qs,\nu)$ are the frequency and eigenvector of the
phonon mode $(qs)$, respectively. $N_{0}$ is the number of q-points in the
Brillouin zone (BZ) sampled. $n_{qs}$ is the Bose-Einstein distribution
corresponding to the temperature $T$. The displacement is consistent with the
canonical ensemble, which takes into account complete quantum statistics,
thereby including all finite temperature effects. Additionally, it also
accounts for the zero-point motions of atoms, which is not achievable in AIMD.
These supercells of displaced atoms are input into an ab initio force
calculator to calculate the force on each atom. We employed the SIESTA[35]
package to carry out the force computation at the Gamma point for each
stochastically perturbed supercell. A limitation of utilizing Eq. II.3 is that
it necessitates prior knowledge of phonon frequencies ($\omega_{qs}$)and
eigenvectors ($\epsilon_{\alpha}(qs,\nu)$), which are the calculations we aim
to perform. As such, a self-consistent approach is utilized. The procedure is
initiated using the IFCs obtained through TDEP by conducting AIMD at a very
low temperature (80K for NaCl, AgI, cBN and 3C-SiC studied here). We initiate
TSS calculations from the TDEP-generated IFCs to obtain an initial set of
IFCs. By using these IFCs as feedback, we generate a new set of a force-
displacement dataset through TSS and obtain a new set of IFCs. This procedure
is repeated until the desired accuracy is achieved self consistently. To
determine the IFCs of the next higher temperature, we utilized the IFCs from
the previous temperature and subsequently executed the self-consistent cycle.
By doing so, we only had to run AIMD once to initiate the process.

###  II.4 Renormalization

In the self-consistent phonon renormalization approach, an effective harmonic
IFC is obtained by combining unrenormalized (bare) second-order IFCs and
fourth-order anharmonic IFCs. The statistical operator renormalization
technique produces the equation for the renormalized harmonic IFC, which is
stated as follows[6, 34, 36]:

$\displaystyle\Psi_{2}^{\alpha\beta}(0\mu,N\nu)=\Phi_{2}^{\alpha\beta}(0\mu,N\nu)+$
$\displaystyle\frac{\hbar}{4N_{0}}\sum_{N_{3}\pi}\sum_{N_{4}\rho}\sum_{\gamma\delta}\sum_{qs}\Phi_{4}^{\alpha\beta\gamma\delta}(0\mu,N\nu,N_{3}\pi,N_{4}\rho)$
$\displaystyle\times
e^{i\mathbf{q}\cdot(R(N_{3})-R(N_{4}))}E_{\gamma}(qs,\rho)E_{\delta}(qs,\pi)\frac{2n_{qs}+1}{\Omega_{qs}\sqrt{M_{\rho}M_{\pi}}}$
(5)

Where $\Phi_{2}^{\alpha\beta}(0\mu,N\nu)$ and
$\Phi_{4}^{\alpha\beta\gamma\delta}(0\mu,N\nu,N_{3}\rho,N_{4}\pi)$ stands for
bare second- and fourth-order IFCs, and $\Psi_{2}^{\alpha\beta}(0\mu,N\nu)$ is
the renormalized second-order IFCs. The renormalization process can also be
realized by utilizing the one-phonon propagator approach[26]. The expression
mentioned above is obtained by truncating an infinite series that includes
contributions from all even-order IFCs. As the statistical average is
calculated over the effective harmonic phonon states, $\Omega_{qs}$ and
$E_{\alpha}(qs,\nu)$ in Eq. II.4, represent the eigenvalue and eigenvector of
renormalized phonon, which are unknown initially. As such, the problem is
solved self-consistently. We start the self-consistent loop from the
eigenvalue and eigenvectors of bare second-order IFCs
($\Phi_{2}^{\alpha\beta}(0\mu,N\nu)$). During the following iterations, the
second term in Eq. II.4 is revised using the eigenvalues and eigenvectors of
the renormalized harmonic IFCs ($\Psi_{2}^{\alpha\beta}(0\mu,N\nu)$) from the
previous iteration. The loop is repeated until the desired level of precision
is achieved. It is important to note that, thus far, no interaction among the
phonon states has been considered. The present study does not delve into the
impact of the interaction among self-consistent phonon states in the energy
levels, mediated by the odd terms in the crystal potential, as it is outside
the scope of this work. Following the procedure described in Ref.[34], we
extend the renormalization effect to anharmonic IFCs of third and fourth
order. The initial step to achieve this involves the multiplication of
recently constructed renormalized second-order IFCs
($\Psi_{2}^{\alpha\beta}(0\mu,N\nu)$) with the available displacement dataset
($u_{\alpha}(N\nu)$) from TSS (or AIMD), which enables the calculation of the
force ($f_{ren.}^{2nd}$) attributed to the renormalized harmonic IFCs.
Afterward, we subtract the calculated force from the actual ab initio force
($f_{\textit{abinitio}}$), which has already been obtained from TSS (or AIMD).
The rest, which remains in Eq. 2, can be written as:

$\displaystyle
f_{\textit{abinitio}}-f_{ren.}^{2nd}=-\sum\limits_{{N_{2}\nu},{N_{3}\pi}}\sum\limits_{\beta,\gamma}\frac{1}{2!}\Psi_{3}^{\alpha\beta\gamma}({0\mu},{N_{2}\nu},{N_{3}\pi})u_{\beta}(N_{2}\nu)u_{\gamma}(N_{3}\pi)-\sum\limits_{{N_{2}\nu},{N_{3}\pi},{N_{4}\rho}}\sum\limits_{\beta,\gamma,\delta}\frac{1}{3!}\Psi_{4}^{\alpha\beta\gamma\delta}({0\mu},{N_{2}\nu},{N_{3}\pi},{N_{4}\rho})$
(6) $\displaystyle\times
u_{\beta}(N_{2}\nu)u_{\gamma}(N_{3}\pi)u_{\delta}(N_{4}\rho)$

By fitting Eq. 6 with the displacement dataset already obtained from TSS (or
AIMD), we can determine the renormalized third-
($\Psi_{3}^{\alpha\beta\gamma}({0\mu},{N_{2}\nu},{N_{3}\pi})$) and fourth-
order
($\Psi_{4}^{\alpha\beta\gamma\delta}({0\mu},{N_{2}\nu},{N_{3}\pi},{N_{4}\rho})$)
interatomic force constants.

###  II.5 Free Energy

The expression of Helmholtz free energy up to fourth order anharmonicity is
given by[6]:

$\displaystyle F=U_{0}+F_{H}+F_{A},$ (7)

where harmonic free energy,

$\displaystyle
F_{H}=\sum\limits_{q,s}\left[\frac{1}{2}\hbar\omega_{qs}+k_{B}Tln\left(1-e^{-\frac{\hbar\omega_{qs}}{k_{B}T}}\right)\right],$
(8)

and the anharmonic contribution to free energy:

$\displaystyle
F_{A}=12\sum\limits_{qs,q^{{}^{\prime}}s^{{}^{\prime}}}\phi(qs,-qs,q^{{}^{\prime}}s^{{}^{\prime}},-q^{{}^{\prime}}s^{{}^{\prime}})\left(n_{qs}+\frac{1}{2}\right)\left(n_{q^{{}^{\prime}}s^{{}^{\prime}}}+\frac{1}{2}\right)-$
(9)
$\displaystyle\frac{18}{\hbar}\sum\limits_{qs,q^{{}^{\prime}}s^{{}^{\prime}},q^{{}^{\prime\prime}}s^{{}^{\prime\prime}}}\Bigg{[}|\phi(qs,q^{{}^{\prime}}s^{{}^{\prime}},q^{{}^{\prime\prime}}s^{{}^{\prime\prime}})|^{2}\left\\{\frac{n_{qs}n_{q^{{}^{\prime}}s^{{}^{\prime}}}+n_{qs}+\frac{1}{3}}{(\omega_{qs}+\omega_{q^{{}^{\prime}}s^{{}^{\prime}}}+\omega_{q^{{}^{\prime\prime}}s^{{}^{\prime\prime}}})_{p}}+\frac{2n_{qs}n_{q^{{}^{\prime\prime}}s^{{}^{\prime\prime}}}-n_{qs}n_{q^{{}^{\prime}}s^{{}^{\prime}}}+n_{q^{{}^{\prime\prime}}s^{{}^{\prime\prime}}}}{(\omega_{qs}+\omega_{q^{{}^{\prime}}s^{{}^{\prime}}}-\omega_{q^{{}^{\prime\prime}}s^{{}^{\prime\prime}}})_{p}}\right\\}+$
$\displaystyle
2\phi(qs,-qs,q^{{}^{\prime\prime}}s^{{}^{\prime\prime}})\phi(q^{{}^{\prime}}s^{{}^{\prime}},-q^{{}^{\prime}}s^{{}^{\prime}},-q^{{}^{\prime\prime}}s^{{}^{\prime\prime}})\frac{n_{qs}n_{q^{{}^{\prime}}s^{{}^{\prime}}}+n_{qs}+\frac{1}{4}}{(\omega_{q^{{}^{\prime\prime}}s^{{}^{\prime\prime}}})_{p}}\Bigg{]}.$

Here
$\phi(qs,q^{{}^{\prime}}s^{{}^{\prime}},q^{{}^{\prime\prime}}s^{{}^{\prime\prime}})$
and
$\phi(qs,q^{{}^{\prime}}s^{{}^{\prime}},q^{{}^{\prime\prime}}s^{{}^{\prime\prime}},q^{{}^{\prime\prime\prime}}s^{{}^{\prime\prime\prime}})$
are scattering matrix elements related to third and fourth-order IFCs:

$\displaystyle\phi(qs,q^{{}^{\prime}}s^{{}^{\prime}},q^{{}^{\prime\prime}}s^{{}^{\prime\prime}})=\frac{1}{3!}\left(\frac{\hbar}{2}\right)^{3/2}\frac{1}{\sqrt{N_{0}}}\frac{1}{\sqrt{\omega_{qs}\omega_{q^{{}^{\prime}}s^{{}^{\prime}}}\omega_{q^{{}^{\prime\prime}}s^{{}^{\prime\prime}}}}}\sum\limits_{\mu,{N_{2}\nu},{N_{3}\pi}}\sum\limits_{\alpha,\beta,\gamma}\Phi_{3}^{\alpha\beta\gamma}({0\mu},{N_{2}\nu},{N_{3}\pi})e^{i\vec{q^{{}^{\prime}}}\cdot\vec{R(N_{2})}}e^{i\vec{q^{{}^{\prime\prime}}}\cdot\vec{R(N_{3})}}\times$
(10)
$\displaystyle\frac{\epsilon_{\alpha}(qs,\mu)\epsilon_{\beta}(q^{{}^{\prime}}s^{{}^{\prime}},\nu)\epsilon_{\gamma}(q^{{}^{\prime\prime}}s^{{}^{\prime\prime}},\pi)}{\sqrt{M_{\mu}M_{\nu}M_{\pi}}},$

and,

$\displaystyle\phi(qs,q^{{}^{\prime}}s^{{}^{\prime}},q^{{}^{\prime\prime}}s^{{}^{\prime\prime}},q^{{}^{\prime\prime\prime}}s^{{}^{\prime\prime\prime}})=\frac{1}{4!}\left(\frac{\hbar}{2}\right)^{2}\frac{1}{N_{0}}\frac{1}{\sqrt{\omega_{qs}\omega_{q^{{}^{\prime}}s^{{}^{\prime}}}\omega_{q^{{}^{\prime\prime}}s^{{}^{\prime\prime}}}\omega_{q^{{}^{\prime\prime\prime}}s^{{}^{\prime\prime\prime}}}}}\sum\limits_{\mu,{N_{2}\nu},{N_{3}\pi},{N_{4}\rho}}\sum\limits_{\alpha,\beta,\gamma,\delta}\Phi_{4}^{\alpha\beta\gamma\delta}({0\mu},{N_{2}\nu},{N_{3}\pi},{N_{4}\rho})\times$
(11) $\displaystyle
e^{i\vec{q^{{}^{\prime}}}\cdot\vec{R(N_{2})}}e^{i\vec{q^{{}^{\prime\prime}}}\cdot\vec{R(N_{3})}}e^{i\vec{q^{{}^{\prime\prime\prime}}}\cdot\vec{R(N_{4})}}\times\frac{\epsilon_{\alpha}(qs,\mu)\epsilon_{\beta}(q^{{}^{\prime}}s^{{}^{\prime}},\nu)\epsilon_{\gamma}(q^{{}^{\prime\prime}}s^{{}^{\prime\prime}},\pi)\epsilon_{\delta}(q^{{}^{\prime\prime\prime}}s^{{}^{\prime\prime\prime}},\rho)}{\sqrt{M_{\mu}M_{\nu}M_{\pi}M_{\rho}}}.$

The given expressions for scattering matrix elements incorporate quasi-
momentum conservation, where the relations
$q+q^{{}^{\prime}}+q^{{}^{\prime\prime}}=G$ and
$q+q^{{}^{\prime}}+q^{{}^{\prime\prime}}+q^{{}^{\prime\prime\prime}}=G$ hold,
with $G$ representing the reciprocal lattice vector of the crystal.

###  II.6 Thermal Conductivity

The non-equilibrium distribution function obtained from the solution of the
linearized Peierls-Boltzmann transport equation (PBTE) is used to calculate
thermal conductivity ($\kappa$). The sources of the scattering term in PBTE
are currently limited to three-phonon (3-ph), four-phonon (4-ph), and isotopic
impurity[37] scattering in our implementation. The three-phonon scattering
rate, as given by Fermi golden rule:

$\displaystyle\frac{1}{\tau_{3-ph}(qs)}=\frac{2\pi}{\hbar^{2}}\sum\limits_{q_{1}s_{1},q_{2}s_{2}}\Bigg{[}\lvert\phi(qs,q_{1}s_{1},-q_{2}s_{2})\rvert^{2}(n_{q_{1}s_{1}}-n_{q_{2}s_{2}})\delta(\omega_{qs}+\omega_{q_{1}s_{1}}-\omega_{q_{2}s_{2}})+\frac{1}{2}\lvert\phi(qs,-q_{1}s_{1},-q_{2}s_{2})\rvert^{2}\times$
(12)
$\displaystyle(1+n_{q_{1}s_{1}}+n_{q_{2}s_{2}})\delta(\omega_{qs}-\omega_{q_{1}s_{1}}-\omega_{q_{2}s_{2}})\Bigg{]},$

the four-phonon scattering rate:

$\displaystyle\frac{1}{\tau_{4-ph}(qs)}=\frac{2\pi}{\hbar^{2}}\sum\limits_{q_{1}s_{1},q_{2}s_{2},q_{3}s_{3}}\Bigg{[}\frac{1}{6}\lvert\phi(qs,-q_{1}s_{1},-q_{2}s_{2},-q_{3}s_{3})\rvert^{2}\frac{n_{q_{1}s_{1}}n_{q_{2}s_{2}}n_{q_{3}s_{3}}}{n_{qs}}\delta(\omega_{qs}-\omega_{q_{1}s_{1}}-\omega_{q_{2}s_{2}}-\omega_{q_{3}s_{3}})+$
(13)
$\displaystyle\frac{1}{2}\lvert\phi(qs,q_{1}s_{1},-q_{2}s_{2},-q_{3}s_{3})\rvert^{2}\frac{(1+n_{q_{1}s_{1}})n_{q_{2}s_{2}}n_{q_{3}s_{3}}}{n_{qs}}\delta(\omega_{qs}+\omega_{q_{1}s_{1}}-\omega_{q_{2}s_{2}}-\omega_{q_{3}s_{3}})+$
$\displaystyle\frac{1}{2}\lvert\phi(qs,q_{1}s_{1},q_{2}s_{2},-q_{3}s_{3})\rvert^{2}\frac{(1+n_{q_{1}s_{1}})(1+n_{q_{2}s_{2}})n_{q_{3}s_{3}}}{n_{qs}}\delta(\omega_{qs}+\omega_{q_{1}s_{1}}+\omega_{q_{2}s_{2}}-\omega_{q_{3}s_{3}})\Bigg{]},$

and the phonon-isotopic impurity scattering rate:

$\displaystyle\frac{1}{\tau_{ph-
iso}(qs)}=\frac{\pi\omega_{qs}^{2}}{2N_{0}}\sum\limits_{q_{1}s_{1}}\sum\limits_{\mu}g_{2}(\mu)\>\lvert\epsilon(qs,\mu)^{*}\cdot\epsilon(q_{1}s_{1},\mu)\rvert^{2}\delta(\omega_{qs}-\omega_{q_{1}s_{1}}),$
(14)

where $g_{2}(\mu)=\sum\limits_{i}f_{\mu i}(1-M_{\mu i}/\bar{M_{\mu}})^{2}$,
with $M_{\mu i}$ and $f_{\mu i}$ being the mass and fractional concentration
of the ith isotope of $\mu$th basis atom and
$\bar{M_{\mu}}=\sum\limits_{i}f_{\mu i}M_{\mu i}$ is the average mass. The
total scattering rate, as given by Matthiessen’s rule:

$\displaystyle\frac{1}{\tau(qs)}=\frac{1}{\tau_{3-ph}(qs)}+\frac{1}{\tau_{4-ph}(qs)}+\frac{1}{\tau_{ph-
iso}(qs)}.$ (15)

The expression for thermal conductivity tensor ($\bm{\kappa}$) in relaxation
time approximation (RTA) is given by:

$\displaystyle\bm{\kappa^{RTA}}=\frac{1}{N_{0}V_{0}}\sum\limits_{q,s}\frac{\partial
n_{qs}}{\partial
T}\>\hbar\omega_{qs}\>\tau(qs)\>\vec{v}_{qs}\otimes\vec{v}_{qs}.$ (16)

The thermal conductivity tensor from the iterative solution of PBTE is
represented by the following expression:

$\displaystyle\bm{\kappa}=\frac{k_{B}T^{2}}{N_{0}V_{0}}\sum\limits_{q,s}\frac{\partial
n_{qs}}{\partial T}\vec{v}_{qs}\otimes\vec{\Delta}_{qs}$ (17)

Here $N_{0}$ is the number of q-points in the BZ sampled, $V_{0}$ is the unit
cell volume, and
$\vec{v}_{qs}=\frac{1}{2\omega_{qs}}\langle\epsilon(qs)|\frac{\partial
D}{\partial\vec{q}}|\epsilon(qs)\rangle$ is phonon group velocity with $D$
being dynamical matrix (Eq. II.2). The self-consistent equation governing the
vector function ($\vec{\Delta}_{qs}$) is interconnected with the phonon
scattering rate and 3-ph, 4-ph, and phonon-isotope scattering probabilities.
For further details regarding these terms, readers can consult the following
references[38, 39, 40, 34]. The novel aspects of our theoretical framework are
as follows: (1) The iterative solution of PBTE[41], which differentiates the
momentum-conserving normal processes as non-resistive. (2) The calculation of
phonon group velocity involves the exact analytical differentiation of the
long-range Ewald term[31] for polar materials. This distinguishes our
implementation from the others, which typically differentiate direction-
dependent non-analytic terms for determining phonon group velocity in polar
materials. (3) In the treatment of energy conservation within the Brillouin
zone (BZ) sum, the tetrahedron method[42] is employed. The tetrahedron method
offers several advantages over the Gaussian or Lorentzian approximation of the
Dirac-delta function. One notable advantage is that it does not rely on any
adjustable parameters and exhibits slightly improved convergence with the
number of q-points. In this study, the convergence of BZ integration for
thermal conductivity is achieved using a q-mesh grid of $17\times 17\times 17$
for the 3-phonon scattering processes. The phonon linewidth ($\Gamma(qs)$)
discussed in this paper is defined as:

$\displaystyle\Gamma(qs)=\frac{1}{2\>\tau(qs)}.$ (18)

For all the ab initio calculations, we utilized SIESTA[35], a code based on
atomic-orbital basis. To expand the wave functions, we employ the polarized
double-zeta basis set. Norm-conserving pseudopotentials are utilized, and
exchange-correlation treatment follows the GGA-PBE[43] approach.

Figure 1: (a) Phonon dispersion of cBN along the high symmetry path calculated
from QHA and renormalized IFCs at 300K. Observe that there is no significant
difference between QHA and renormalization. The experimental data are
extracted from Exp.1[44], Exp.2[45], and Exp.3[46]. (b) - (e) Variation in
phonon line width along the $\Gamma\to X$ high-symmetry path for one acoustic
branch (TA) and three optic branches (TO, ZO, and LO) in cBN across
temperatures ranging from 100K to 500K. (f) Thermal conductivity of cBN
calculated across the temperature range of 100K to 500K. Renormalization
results in approximately a 2% to 3% difference in the $\kappa$ value compared
to the QHA. Experimental data are obtained from Exp. (2020)[47] and Exp.
(2018)[48]. Inset plot: Contribution from the three acoustic phonon branches
to the thermal conductivity at different temperatures. (g) Phonon scattering
rates, considering 3-phonon and phonon-isotope scattering, calculated for cBN
at 300K and 500K using renormalized IFCs. The elevated scattering rates at
higher temperatures contribute to a reduction in thermal conductivity. (h)
Spectrally resolved thermal conductivity of cBN at various temperatures,
providing insights into the scattering rate plot. Notably, the renormalized
spectral $\kappa$ exhibits a slight increase across all temperatures.

##  III  RESULTS

###  III.1  High $\kappa$ Material: cBN and 3C-SiC

We selected cubic-Boron Nitride (cBN) and cubic-Silicon Carbide (3C-SiC) as
materials characterized by low anharmonicity in their interatomic potential.
High thermal conductivity in these materials is attributed to the combination
of weak anharmonicity and the light atomic mass of their constituent atoms.
Figure 1(a) and Figure 2(a) display phonon dispersion for cBN and 3C-SiC,
respectively, presenting results from both the QHA and renormalized IFCs at
300K. Note that in QHA, the temperature effect on the lattice constant is
taken into account, and the displacement-force dataset (Eq. 2) is fitted up to
the fourth order. There is little to no noticeable difference between the
phonon frequencies obtained through the QHA and the renormalized approach.
This signifies the existence of weak higher-order anharmonic terms in these
materials, which appear in the renormalization process. The figures also
include the experimental phonon frequency for cBN[44, 45, 46] and 3C-SiC[49],
demonstrating a close agreement with our calculations. Figure 1(b-e) and
Figure 2(b-e) illustrate the temperature-dependent phonon linewidth for one
acoustic branch and optic branches along the high-symmetry $\Gamma\to X$ path
of cBN and 3C-SiC, respectively. As depicted in the figures, the phonon
linewidth exhibits a monotonic increase with temperature, resulting in reduced
thermal transport as the temperature rises. Observe that the phonon linewidth
for the acoustic branch is significantly smaller compared to that of the optic
branches, suggesting lower scattering for the acoustic phonons.

Figure 2: (a) Phonon dispersion of 3C-SiC along the high symmetry path
calculated from QHA and renormalized IFCs at 300K. Observe that there is no
significant difference between QHA and renormalization. The experimental data
are extracted from Ref.[49] (b) - (e) Variation in phonon linewidth along the
$\Gamma\to X$ high-symmetry path for one acoustic branch (LA) and three optic
branches (TO, ZO, and LO) in 3C-SiC across temperatures ranging from 100K to
500K. (f) Thermal conductivity of 3C-SiC calculated across the temperature
range of 100K to 500K. Renormalization results in a marginally higher $\kappa$
value compared to the QHA. Experimental data are obtained from Exp.1[50] and
Exp.2[51]. Note that Exp.2[51] corresponds to a polycrystalline sample, which
tends to match better with our calculation at higher temperatures. Inset plot:
Contribution from the three acoustic phonon branches to the thermal
conductivity at different temperatures. (g) Phonon scattering rates,
considering 3-phonon + phonon-isotope scattering, calculated for 3C-SiC at
300K and 500K using renormalized IFCs. The elevated scattering rates at higher
temperatures contribute to a reduction in thermal conductivity. (h) Spectrally
resolved thermal conductivity of 3C-SiC at various temperatures. Notably, the
renormalized spectral $\kappa$ exhibits a slight increase across all
temperatures.

This difference contributes to the high thermal conductivity of acoustic
phonons, given the inverse relationship between thermal conductivity and
linewidth. More than 90% of the total $\kappa$ value for both materials is
contributed by three acoustic modes, as depicted in the inset plot of Figure
1(f) and Figure 2(f). Also, observe that certain peak positions occur within
the linewidth curve along the high-symmetry path. These peaks emerge as both
momentum conservation (Eq. 10) and energy conservation (Eq. 12) are met to a
greater extent, leading to heightened scattering for these phonon modes
($qs$). Consequently, the linewidth of these phonon modes becomes more
pronounced, manifesting as peaks within the linewidth curve. The temperature-
dependent thermal conductivity of cBN and 3C-SiC is presented in Figure 1(f)
and Figure 2(f) across the 100K to 500K range. Note that for these materials
characterized by high thermal conductivity, renormalization of the IFCs has
minimal to no impact on the $\kappa$ values. For instance, cBN exhibits a
thermal conductivity of 648.31 W/m.K at 300K, according to the QHA. Meanwhile,
the renormalization of IFCs predicts $\kappa$ at 666.97 W/m.K, representing a
roughly 3% increase compared to the QHA-derived value. The experimental
determination of the $\kappa$ of cBN is 648$\pm$52 W/m.K as extracted from
Ref.[48] and 850 $\pm$ 90 W/m.K as extracted from Ref.[47] Note that in our
calculation of phonon-isotope scattering, we considered the natural abundance
of the isotopes of the constituent atoms, and the experimental values
extracted also correspond to the natural isotopic concentration in the sample.
For 3C-SiC, our calculation using the QHA yields a thermal conductivity value
of 364.25 W/m.K, while the renormalization predicts a slightly elevated
$\kappa$ of 378.68 W/m.K at 300K. The experimental $\kappa$ of single-
crystalline 3C-SiC as extracted from the Ref.[50] is 437 W/m.K. Figure 2(f)
also includes the thermal conductivity of polycrystalline 3C-SiC, as extracted
from the Ref.[51], which tends to match our calculation at higher
temperatures. The bar plots inset in Figure 1(f) and Figure 2(f) demonstrate
the contribution of the three acoustic phonon branches to thermal conductivity
at various temperatures. With increasing temperature, higher-frequency phonon
branches and optic branches tend to contribute more to thermal conduction.
Figure 1(g) and Figure 2(g) illustrate the three-phonon + phonon-isotope
scattering rates at temperatures of 300K and 500K for cBN and 3C-SiC,
respectively. The elevation in scattering rate as temperature increases
accounts for the reduction in thermal conductivity. To translate the
scattering rate plot into thermal conductivity, we computed spectral thermal
conductivity, as shown in Figure 1(h) and Figure 2(h). Observe that the
spectral thermal conductivity tends to decrease as temperature rises,
resulting in a reduction in the overall thermal conductivity. Additionally,
note that the spectral thermal conductivity for renormalized IFCs is slightly
higher than that of the QHA for all the temperatures. This accounts for the
marginal increase in thermal conductivity observed with the renormalized
approach in our calculations.

###  III.2  Low $\kappa$ Material: NaCl and AgI

To demonstrate the practicality of our numerical framework for materials with
low thermal conductivity and significant anharmonicity, we selected sodium
chloride (NaCl) and silver iodide ($\gamma$-AgI, space group - $F\bar{4}3m$)
as representative materials. Being a member of the alkali halides family, NaCl
exhibits high anharmonicity in its inter-atomic potential due to the weakly
bonded, highly ionic nature of its constituent atoms. Furthermore, it displays
remarkably low thermal conductivity, making it an ideal candidate for studying
the anharmonicity within our current framework. Figure 3(a) presents the
phonon dispersion of NaCl crystal at 300K, calculated from QHA and through the
renormalization of the phonon quasiparticle as described in the method
section. Our findings indicate that the renormalization process significantly
increases the phonon frequency of the optical modes of NaCl. Such behavior is
not observed in cBN and 3C-SiC. The acoustic branches remain largely
unaffected after renormalization. The experimentally measured data[52] through
the inelastic neutron scattering demonstrates a close agreement with the
renormalized phonon frequency. This emphasizes the significance of treating
phonons as quasiparticles through the renormalization process, rather than
considering them solely as normal modes of lattice vibrations, especially for
highly anharmonic materials. Figure 3(b) to Figure 3(d) show the dependence of
phonon dispersion with temperature from 100K to 400K, both for bare phonon
within QHA and renormalized phonon. The outcome reveals that, in comparison to
the renormalized phonon, the bare phonon experiences a higher degree of
softening in the optical phonon branches as the temperature increases. In
Figure 3(d), the plot illustrates the variation of phonon frequency at the
$\Gamma$-point of the Brillouin zone (BZ) with temperature for both the
transverse optical (TO) and longitudinal optical (LO) branches. In contrast to
QHA, the renormalized phonon demonstrates a minimal temperature dependence
while remaining consistent with the experimental data. Figure 3(e) illustrates
the distribution of phonon linewidth or inverse phonon lifetime along the
$\Gamma$-X high-symmetry path in the BZ due to the three-phonon scattering
process at 300K. The quantitative value of the linewidth can be obtained from
the vertical extent of the highlighted area around the phonon branch, with the
unit being the same as the phonon frequency (THz). In Figure 3(f), it can be
observed that for the temperature range of 100K-400K, the linewidth of the LO
branch at the $\Gamma$-point is greater in the case of the renormalized
phonon. Conversely, the renormalized linewidth of the TO branch consistently
remains lower compared to that of the bare phonon. It should be noted that, in
our approach, anharmonic third-order IFCs are also subject to renormalization,
according to Eq. 5. The increased linewidth of the LO branch indicates a
higher rate of 3-phonon scattering, resulting in a reduced contribution to
thermal conduction.

Figure 3: (a) Phonon dispersion of NaCl along the high-symmetry path at 300K.
The renormalization leads to an upward shift in the frequency of the optical
branches, aligning with the experimental data.[52] (b) Temperature dependence
of the bare phonon dispersion within the 100K-400K range. Observe the
pronounced softening of optical branches with increasing temperature. (c) The
renormalization leads to weaker temperature dependence phonon dispersion
consistent with the experimental data. (d) Temperature dependence of the
phonon frequency at the $\Gamma$ point for LO and TO branches. Notice the
reduced extent of softening with increasing temperature due to
renormalization, closely matching the experimental data.[53] (e) Distribution
of phonon linewidth along the $\Gamma$-X path at 300K. The vertical spread of
the shaded region surrounding each phonon branch determines the linewidth. (f)
Temperature dependence of the phonon linewidth at the $\Gamma$ point for LO
and TO branches. Figure 4: (a) Thermal conductivity of NaCl calculated using
the QHA and renormalization. Notably, renormalization consistently
overestimates the thermal conductivity across all temperatures. Experimental
values are extracted from Ref.[54, 55, 56] (b) Three-phonon scattering rate
calculated from bare and renormalized IFCs at 300K. (c) Spectral distribution
of thermal conductivity at 100K and 300K calculated through QHA and
renormalization.

Figure 4(a) presents the thermal conductivity ($\kappa$) in the 100K to 400K
temperature range calculated from bare and renormalized IFCs. The thermal
conductivity calculated from the three-phonon scattering process of the
renormalized phonon consistently overestimates the values obtained from QHA
and experimental data. For example, the renormalized $\kappa$ at 100K is 28.78
W/mK, exhibiting an approximately 17% increase compared to the unrenormalized
$\kappa$ of 24.62 W/mK. At 300K, the renormalization results in an
overprediction of about 25%.

Figure 5: (a) Phonon dispersion of AgI along the high symmetry path calculated
from QHA and renormalized IFCs at 300K. (b) - (e) Variation in phonon
linewidth for one acoustic branch (ZA) and three optic branches (TO, ZO, and
LO) in AgI across temperatures ranging from 100K to 450K. (f) Thermal
conductivity of AgI calculated over the temperature range of 100K to 450K
using 3-phonon scattering. Results from four-phonon scattering processes are
presented for the $11^{3}$ q-grid for 200 K, 300 K, and 400 K temperatures.
Experimental data are sourced from Exp. (1982)[57] and Exp. (2023)[58]. It is
noteworthy that both the QHA and renormalized values tend to over-predict the
measured $\kappa$ values in three-phonon, whereas the four-phonon prediction
closely aligns with the experimental data. Inset plot: Convergence of $\kappa$
with q-grid in BZ for three-phonon and four-phonon processes. (g) Phonon
scattering rates, taking into account 3-phonon + phonon-isotope scattering at
$11^{3}$ and $17^{3}$ q-grid, and 3-phonon + 4-phonon + phonon-isotope
scattering processes at $11^{3}$ q-grid. The highlighted regions indicate a
significant increase in the scattering rate, approximately an order of
magnitude higher for four-phonon processes compared to three-phonon processes.
Notably, the two dips in the scattering rate in the highlighted region
correspond to the two peaks in the spectral $\kappa$. (h) Spectrally resolved
thermal conductivity of AgI for three-phonon and four-phonon scattering,
illustrating the notable decrease in $\kappa$ when fourth-order processes are
considered.

Importantly, our implementation of the tetrahedron method for energy
conservation delta function does not involve any adjustable parameters for the
thermal conductivity calculation, distinguishing it from the Gaussian or
Lorentzian approximations. On the other hand, the thermal conductivity
obtained from the bare IFCs exhibits a close agreement with the experimental
data, even in the scenario involving minimal three-phonon scattering
processes. At a temperature of 300K, the $\kappa$ derived from the QHA is in
close proximity to the experimental results, with a deviation of less than 7%
from the findings in Ref[55] and less than 16% from Ref[54]. To interpret this
anomalous behavior of renormalized phonon, we calculated the three-phonon
scattering rate from renormalized and bare IFCs, as presented in Figure 4(b).
As illustrated in the figure, the scattering rate derived from the bare IFCs
surpasses the scattering rate of the renormalized IFCs. Furthermore, the bare
phonon exhibits a greater extent of temperature-induced softening, as depicted
in Figure 3(a-d). By combining these two factors, the thermal conductivity
predicted from the unrenormalized IFCs consistently remains lower than that
from the renormalized IFCs, since $\kappa$ is directly proportional to the
phonon frequency and inversely proportional to the scattering rate. This can
be observed in Figure 4(c), which illustrates the spectral distribution of
$\kappa$ at 100K and 300 K. At both temperatures, the renormalized spectral
$\kappa$ consistently surpasses the values obtained from the unrenormalized
spectral $\kappa$ across all phonon frequency ranges, reinforcing the higher
thermal conductivity exhibited by the renormalized phonons. The peculiar
characteristics of the renormalized phonon can be addressed by incorporating
higher-order phonon scattering processes, as elaborated in the Ref.[34],
especially for highly anharmonic materials like NaCl. The inclusion of four-
phonon scattering processes provides a means to address this issue.

Figure 5 presents the thermal transport-related physical properties of AgI.
Research interest has recently grown in the thermal transport of materials
exhibiting ultralow thermal conductivity[59, 60], primarily due to their
potential applications in thermoelectric devices. Among these materials, AgI
stands out as a promising candidate. Understanding the ultralow thermal
transport in AgI has been enhanced by recent studies[58, 61] highlighting the
significance of higher-order phonon scattering processes. Figure 5(a) presents
the phonon dispersion of AgI along the high symmetry path of BZ at 300K
calculated from both the unrenormalized and renormalized IFCs. Note that the
phonon frequencies experience a slight stiffening across both the acoustic and
optical branches when considering renormalized IFCs. The increase in phonon
frequency in AgI resulting from renormalization is of a smaller magnitude than
that observed in NaCl. Figure 5(b-e) depicts temperature-dependent phonon
linewidth for one acoustic branch and three optic branches along the high
symmetry $\Gamma\to X$ path for the 100K to 450K temperature range. The
monotonic increase of the linewidth with temperature implies higher phonon
scattering as temperature rises. Figure 5(f) presents the thermal conductivity
of AgI for the temperature range 100K to 450K, calculated using both QHA and
renormalized IFCs. The thermal conductivity calculations, relying solely on
anharmonic 3-phonon and phonon-isotope scattering as the primary sources of
scattering in PBTE, are represented by solid red and blue lines for QHA and
renormalized IFCs, respectively. The experimental values, as extracted from
the Ref.[57, 58], are also presented in the figure. Observe that our
calculated $\kappa$ values from three-phonon scattering processes consistently
exceed the experimental values across all temperature ranges, and the process
of renormalization leads to an additional overprediction compared to the QHA
values. For example, our calculated $\kappa$ of AgI at 300K is 1.03 W/m.K,
whereas the experimental $\kappa$ is 0.36 W/m.K, as reported in Ref.[58] and
0.41 W/m.K as extracted from the Ref[57]. The renormalization predicts
$\kappa$ as 1.17 W/m.K, which is about 13% higher than the QHA value. The
discrepancy between the calculated thermal conductivity and experimental
observations highlights the significance of higher-order phonon scattering
processes (specifically, 4-phonon interactions) for highly anharmonic
materials such as AgI. To enhance the robustness of our numerical framework in
handling materials with extreme anharmonicity, we incorporated the computation
of four-phonon scattering processes for AgI. Figure 5(f) displays the computed
thermal conductivity, taking into account four-phonon scattering, for
temperatures of 200 K, 300 K, and 400 K. Take note of the impressive alignment
between the computed $\kappa$ and experimental values achieved simply by
integrating the fourth-order scattering process. Our calculation predicts the
$\kappa$ of AgI to be 0.41 W/m.K at 300 K, closely resembling the experimental
values of 0.41 W/m.K[57] and 0.36 W/m.K[58]. The computational expense
increases significantly by several orders of magnitude when incorporating
four-phonon scattering. To address this computational bottleneck, we mitigated
it by leveraging symmetries in the calculation of fourth-order scattering
matrix elements (Eq. 11). In Eq. 11, the first index ($qs$) is constrained to
the irreducible Brillouin Zone (BZ), and the scattering matrix elements for
permutations among the remaining three indices
($q^{{}^{\prime}}s^{{}^{\prime}},q^{{}^{\prime\prime}}s^{{}^{\prime\prime}},q^{{}^{\prime\prime\prime}}s^{{}^{\prime\prime\prime}}$)
are computed only once, as they yield the same values. We also utilized
specific computer hardware acceleration techniques to enhance the speed of the
computation. To provide a comparison of computational costs, our three-phonon
calculation for a $17\times 17\times 17$ ($17^{3}$) q-grid in the BZ required
approximately 6 minutes, utilizing 240 CPU cores distributed across 5 compute
nodes. In contrast, our four-phonon calculation for an $11\times 11\times 11$
($11^{3}$) q-grid took 15 hours, employing 960 CPU cores distributed across 20
compute nodes. The convergence of $\kappa$ with q-grid is depicted in the
inset plot in Figure 5(f). While it requires a $17^{3}$ q-grid for three-
phonon convergence, no significant difference in $\kappa$ values for four-
phonon has been observed between $7^{3}$ and $11^{3}$ q-grid. The convergence
of $\kappa$ at a lower q-grid for four-phonon processes has been well-
established in previous studies[13, 16]. In our calculation, we address energy
conservation using the parameter-free tetrahedron method[42], which achieves
convergence with a smaller q-grid compared to Gaussian or Lorentzian
approximations. In Figure 5(g), the scattering rate is depicted for both
three-phonon+phonon-isotope scattering processes and the additional
contribution from fourth-order scattering mechanisms calculated through the
Fermi golden rule. Note that the scattering rates, encompassing four-phonon
processes, are computed with an $11^{3}$ q-grid. To facilitate a direct
comparison with three-phonon scattering, we have also included the three-
phonon scattering rate for both $11^{3}$ q-grid and $17^{3}$ q-grid. Note that
the inclusion of four-phonon scattering results in an increased scattering
rate compared to the conventional three-phonon scattering. In the highlighted
region of Figure 5(g), the four-phonon scattering rate is approximately an
order of magnitude greater than that of the three-phonon. The rise in phonon
scattering rates results in a decrease in the thermal transport coefficient of
AgI during the four-phonon process, as depicted in Figure 5(f). For a more in-
depth exploration of the impact of higher-order four-phonon scattering on
thermal conductivity, we computed the spectral distribution of $\kappa$, as
depicted in Figure 5(h). Observe that the two peaks in the spectral $\kappa$
align with the two dips in the highlighted region of the phonon scattering
rate in Figure 5(g). The spectral $\kappa$ exhibits a notable reduction in the
four-phonon process across all phonon frequency ranges, leading to an overall
decline in the total thermal conductivity compared to the three-phonon
scenario.

Figure 6: (a) Variation of harmonic and fourth-order free energy with unit
cell volume at 300K. Observe the overlap between $F_{2nd}$ and $F_{4th}$ for
QHA, resulting from the third-order and fourth-order term cancellation. The
inset plot illustrates pressure calculated using the volumetric derivative of
free energy. (b) Thermal conductivity as a function of pressure at 300K,
computed using QHA and renormalized IFCs. Experimental data extracted from
Ref[54, 62]. (c) Phonon dispersion of NaCl at 0 GPa and 3 GPa pressure,
computed at 300K using renormalized Interatomic Force Constants (IFCs).
Observe the significant phonon stiffening, particularly in the optical
branches, under increased pressure. (d) Three-phonon scattering rate computed
at 0 GPa and 3 GPa pressure. The reduction in anharmonic scattering rate with
increased pressure results in higher thermal conductivity at elevated
pressures. (e) Spectral decomposition of thermal conductivity at 0 GPa and 3
GPa pressure. Notice that the spectral $\kappa$ rises due to the lower
scattering rates of phonon modes as pressure increases. (f) and (g) Three-
phonon scattering phase space (AAA, AA0, and AOO) of NaCl at two different
pressures. Note that there is no significant change in AAA and AOO phase space
as pressure increases. The AAO phase space is enhanced for higher-frequency
phonon modes due to the widening of the acoustic-optic gap in the phonon
dispersion, enabling more higher-frequency modes to participate in scattering
events while energy conservation is maintained.

###  III.3  Free energy and pressure dependent thermal conductivity

Figure 6(a) depicts the per unit cell free energy variation as a function of
unit cell volume at 300K. The figure illustrates that the free energy
calculated up to the second and fourth order in QHA is nearly identical across
all volumes. In our observations using QHA, we find that the contribution from
the third-order and fourth-order anharmonic terms to the free energy mutually
nullifies. Conversely, in the renormalized case, such cancellation does not
occur. The points of lowest free energy on the curve correspond to the
equilibrium lattice constant of the crystal at zero pressure and the specific
temperature. By minimizing the fourth-order free energy using the renormalized
IFCs, the lattice constant of NaCl at 300K is determined to be 5.637 $\AA$.
The lattice constant obtained from the renormalized harmonic free energy is at
5.624 $\AA$. In contrast, the QHA predicts the lattice constant to be 5.632
$\AA$. The experimentally measured lattice constant at NaCl at 300K is 5.645
$\AA$[63]. While both the QHA and renormalization techniques underpredict the
lattice constant, the renormalized fourth-order approach shows better
agreement with the measured data. In Figure 6(a), the inset plot illustrates
the variation of unit cell volume with pressure for NaCl at 300K. The pressure
($P$) is determined using the thermodynamic relation involving the free
energy: $P=-\left(\frac{\partial F}{\partial V}\right)_{T}$. The pressure-
dependent thermal conductivity of NaCl at 300K is presented in Figure 6(b). As
the pressure rises from 0 to 3.4 GPa, the thermal conductivity obtained from
QHA experiences an increase, going from 7.53 W/mK to 11.12 W/mK. As explained
in the earlier section, the thermal conductivity obtained from renormalized
IFCs tends to overpredict, primarily because we considered only up to three-
phonon scattering processes. The experimental data[54, 62] also exhibit a
similar increasing trend for thermal conductivity, as seen in the QHA values.
The phonon modes stiffen with the increase in pressure, as depicted in Figure
6(c). It is important to note that the rate of phonon stiffening is more
significant in the optic branches compared to the acoustic branches. In order
to comprehend the rise in thermal conductivity with pressure, we calculated
three-phonon scattering rates at 0 GPa and 3GPa pressure, as depicted in
Figure 6(d). As the pressure increases, the total three-phonon scattering rate
decreases, leading to higher thermal conductivity at elevated pressure levels.
The spectral thermal conductivity, as shown in Figure 6(e), also shifts to
higher values as pressure increases due to reduced scattering rates of phonon
modes. To investigate the anharmonic three-phonon scattering processes
further, we calculated the scattering phase space involving different
polarizations of phonon modes. The scattering phase space considered here is
defined as[64]:

$\displaystyle\mathcal{P}(\omega_{qs})=\frac{2}{3N_{p}^{2}N_{0}}\sum\limits_{q^{{}^{\prime}},s^{{}^{\prime}}}\sum\limits_{q^{{}^{\prime\prime}},s^{{}^{\prime\prime}}}\Big{[}\delta(\omega_{qs}+\omega_{q^{{}^{\prime}}s^{{}^{\prime}}}-\omega_{q^{{}^{\prime\prime}}s^{{}^{\prime\prime}}})+$
$\displaystyle\frac{1}{2}\delta(\omega_{qs}-\omega_{q^{{}^{\prime}}s^{{}^{\prime}}}-\omega_{q^{{}^{\prime\prime}}s^{{}^{\prime\prime}}})\Big{]},$

where $N_{p}$ is the number of polarizations. The above expression implicitly
incorporates quasi-momentum conservation: $q\pm
q^{{}^{\prime}}+q^{{}^{\prime\prime}}=G$, which leads to considering only the
summation over $q^{{}^{\prime}}$. Multiple scattering processes are feasible,
involving various acoustic (A) and optic (O) phonon branch combinations.
Figures 3(f) and 3(g) illustrate the three-phonon scattering phase space for
AAA, AAO, and AOO scattering processes. Due to energy conservation
constraints, scattering processes involving three optical branches (OOO) are
not permitted. As depicted in Figure 6(f), the AAA scattering phase space
remains relatively constant as pressure increases. AAA scattering phase space
at 0 GPa also reveals more feasible scattering events conserving momentum and
energy than at 3 GPa. As depicted in Figure 6(g), the AOO scattering phase
space demonstrates minimal to no variation with the increase in pressure.
Conversely, the AAO scattering phase space displays an increment with
increasing pressure. This is attributed to the stiffening of phonon
dispersion, which results in an increased gap between acoustic and optic
phonons. According to energy conservation principles, only acoustic phonons
with frequencies higher than this are allowed to take part in AAO scattering
processes[65, 66]. So, a greater number of higher frequency phonon modes
participate in the scattering process, as evident from Figure 6(g).

##  IV  CONCLUSIONS

Material | Temperature (K) |  $\kappa$ (QHA) ($Wm^{-1}K^{-1}$) |  $\kappa$ ($Wm^{-1}K^{-1}$) (Renormalized) |  $\kappa$ (Exp.) ($Wm^{-1}K^{-1}$)
---|---|---|---|---
| 300K | 648.31 | 666.97 |  $\sim$ 648$\pm$52[48], 850$\pm$90[47]
cBN | 400K | 512.95 | 527.42 |  $\sim$460$\pm$40[48], $\sim$603$\pm$70[47]
| 500K | 423.64 | 437.94 |  $\sim$360$\pm$42[48], $\sim$483$\pm$62[47]
| 300K | 364.25 | 378.68 |  $\sim$437[50]
3C-SiC | 200K | 677.71 | 692.97 |  $\sim$768[50]
| 250K | 484.02 | 499.28 |  $\sim$552[50]
| 300K | 7.53 | 9.95 | 6.02[54], $\sim$5.8[56], $\sim$8.1[55]
NaCl | 200K | 10.90 | 13.75 |  $\sim$10.27[54], $\sim$11.35[55]
| 100K | 24.62 | 28.78 |  $\sim$22.47[54], $\sim$22.79[55]
| 300K | 1.03 | 1.17, 0.41(4-ph) | 0.36[58], $\sim$0.41[57]
AgI | 200K | 1.54 | 1.74, 0.77(4-ph) |  $\sim$0.61[58], $\sim$0.72[57]
| 400K | 0.78 | 0.88, 0.25(4-ph) |  $\sim$0.21[57]
Table 2:  Comparison of thermal conductivity ($\kappa$) at various
temperatures for the materials examined in this study, calculated via the QHA
and renormalized IFCs, alongside experimental ($\kappa$) values for reference.
The symbol $\sim$ denotes data extracted from corresponding references’
figures, which may contain inherent errors.

In summary, we have developed a comprehensive numerical framework for
understanding phonon and phonon-driven physical properties of crystalline
semiconductors and insulators. The phonon quasiparticle approach forms the
foundation for our self-consistent phonon renormalization technique, which
effectively characterizes both the vibrational attributes of strongly
anharmonic solids and those of harmonic crystals. Our renormalization
technique also extends to the anharmonic third- and fourth-order IFCs. Both
the TDEP and TSS utilized in this study include explicit variations in
temperature when computing IFCs. The investigation of the NaCl crystal in our
study demonstrates a notable alignment between the phonon dispersion obtained
from the renormalized IFCs and the experimental data. On the other hand, the
thermal conductivity values of NaCl exhibit a close correspondence with the
experimental measurements when considering the IFCs obtained through the
quasiharmonic approximation (QHA). Our investigation reveals that both the QHA
and the renormalization approach consistently overestimate the experimentally
measured thermal conductivity of AgI when considering three-phonon scattering
processes. This implies that higher-order phonon scattering processes are
crucial considerations in materials with ultralow thermal conductivity, like
AgI. Consequently, we undertook the computationally demanding task of
calculating fourth-order scattering processes to address these discrepancies.
We observe that materials exhibiting high $\kappa$ experience minimal impact
on thermal conductivity and phonon-driven physical properties through the
renormalization of the IFCs. Both cBN and 3C-SiC, investigated in this study
as materials with high thermal conductivity, demonstrate a 2-3% difference in
$\kappa$ values between QHA and renormalization. These results closely match
experimental observations, even when considering only up to 3-phonon
scattering processes. In our current investigation, the Helmholtz free energy
takes into account anharmonicity up to the fourth order and demonstrates a
close agreement with the experimental lattice constant of the NaCl crystal.
Additionally, we have examined how the thermal conductivity of NaCl changes
under varying pressure conditions. Anharmonic three-phonon scattering rates
and scattering phase space are also analyzed to elucidate the behavior of the
thermal transport coefficient with temperature and pressure.

While the phonon dispersion and its associated properties exhibit favorable
alignment with experimental data through the renormalization approach,
anharmonic characteristics like lattice thermal conductivity consistently
exhibit an overprediction compared to experimental measurements when IFCs are
subjected to renormalization, both for NaCl and AgI. The primary reason behind
these discrepancies is that we have considered only up to the three-phonon
scattering processes. The significance of higher-order anharmonicity cannot be
overlooked in the case of strongly anharmonic materials such as NaCl and AgI.
While our current exploration focuses solely on the four-phonon scattering
processes for AgI, we have omitted the consideration of higher-order
scattering processes for other materials in this study due to computational
limitations. The findings of our present investigation are summarized in Table
2. The effects of electron-phonon scattering processes[67] are not accounted
for in our current study. For certain semiconductors with higher carrier
concentrations at standard room temperature, ignoring the influence of
electron-phonon scattering processes is not viable.

In conclusion, we offer an extensive numerical framework that enables the
exploration of mechanical, thermodynamic, and thermal properties driven by
crystalline vibrations in insulators and semiconductors. This framework is
built upon a sturdy foundation of interatomic force constant (IFC) calculation
and phonon renormalization techniques. The present development enables the
efficient and higher-accuracy numerical and theoretical study for
understanding phonon and phonon-driven physical properties for better and more
efficient application of materials.

##  V  ACKNOWLEDGMENTS

We thank TUE-CMS, IISc, funded by DST, India, as well as Param-Pravega at
SERC, IISc, for providing the computational facilities. The authors thank
Navaneetha Krishnan Ravichandran for valuable discussions. S.M. acknowledges
the generous fellowship provided by MHRD, India. PKM acknowledges funding
through SERB, IRHPA (No. IPA/2020/000034).

## References

  * Alfè [2009] D. Alfè, Phon: A program to calculate phonons using the small displacement method, Computer Physics Communications 180, 2622 (2009).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scripta Materialia 108, 1 (2015).
  * Togo  _et al._ [2015] A. Togo, L. Chaput, and I. Tanaka, Distributions of phonon lifetimes in brillouin zones, Phys. Rev. B 91, 094306 (2015).
  * Li  _et al._ [2014] W. Li, J. Carrete, N. A. Katcho, and N. Mingo, ShengBTE: a solver of the Boltzmann transport equation for phonons, Comp. Phys. Commun. 185, 1747–1758 (2014).
  * Mounet and Marzari [2005] N. Mounet and N. Marzari, First-principles determination of the structural, vibrational and thermodynamic properties of diamond, graphite, and derivatives, Phys. Rev. B 71, 205214 (2005).
  * Wallace [1972] D. C. Wallace,  _Thermodynamics of Crystals_ (John Wiley & Sons, Inc., New York, 1972).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. 73, 515 (2001).
  * Romero  _et al._ [2020] A. H. Romero, D. C. Allan, B. Amadon, G. Antonius, T. Applencourt, L. Baguet, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, F. Bruneval, G. Brunin, D. Caliste, M. Côté, J. Denier, C. Dreyer, P. Ghosez, M. Giantomassi, Y. Gillet, O. Gingras, D. R. Hamann, G. Hautier, F. Jollet, G. Jomard, A. Martin, H. P. C. Miranda, F. Naccarato, G. Petretto, N. A. Pike, V. Planes, S. Prokhorenko, T. Rangel, F. Ricci, G.-M. Rignanese, M. Royo, M. Stengel, M. Torrent, M. J. van Setten, B. Van Troeye, M. J. Verstraete, J. Wiktor, J. W. Zwanziger, and X. Gonze, ABINIT: Overview and focus on selected capabilities, The Journal of Chemical Physics 152, 124102 (2020), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5144261/16711978/124102_1_online.pdf .
  * Li  _et al._ [2012] W. Li, L. Lindsay, D. A. Broido, D. A. Stewart, and N. Mingo, Thermal conductivity of bulk and nanowire mg2sixsn1-x alloys from first principles, Phys. Rev. B 86, 174307 (2012).
  * Togo  _et al._ [2023] A. Togo, L. Chaput, T. Tadano, and I. Tanaka, Implementation strategies in phonopy and phono3py, J. Phys. Condens. Matter 35, 353001 (2023).
  * Togo [2023] A. Togo, First-principles phonon calculations with phonopy and phono3py, J. Phys. Soc. Jpn. 92, 012001 (2023).
  * Paulatto  _et al._ [2013] L. Paulatto, F. Mauri, and M. Lazzeri, Anharmonic properties from a generalized third-order ab initio approach: Theory and applications to graphite and graphene, Phys. Rev. B 87, 214303 (2013).
  * Han  _et al._ [2022] Z. Han, X. Yang, W. Li, T. Feng, and X. Ruan, Fourphonon: An extension module to shengbte for computing four-phonon scattering rates and thermal conductivity, Computer Physics Communications 270, 108179 (2022).
  * Feng and Ruan [2016] T. Feng and X. Ruan, Quantum mechanical prediction of four-phonon scattering rates and reduced thermal conductivity of solids, Phys. Rev. B 93, 045202 (2016).
  * Feng  _et al._ [2017] T. Feng, L. Lindsay, and X. Ruan, Four-phonon scattering significantly reduces intrinsic thermal conductivity of solids, Phys. Rev. B 96, 161201 (2017).
  * Han and Ruan [2023] Z. Han and X. Ruan, Thermal conductivity of monolayer graphene: Convergent and lower than diamond, Phys. Rev. B 108, L121412 (2023).
  * Xia  _et al._ [2020] Y. Xia, V. I. Hegde, K. Pal, X. Hua, D. Gaines, S. Patel, J. He, M. Aykol, and C. Wolverton, High-throughput study of lattice thermal conductivity in binary rocksalt and zinc blende compounds including higher-order anharmonicity, Phys. Rev. X 10, 041029 (2020).
  * Knoop  _et al._ [2023a] F. Knoop, T. A. R. Purcell, M. Scheffler, and C. Carbogno, Anharmonicity in thermal insulators: An analysis from first principles, Phys. Rev. Lett. 130, 236301 (2023a).
  * Knoop  _et al._ [2023b] F. Knoop, M. Scheffler, and C. Carbogno, Ab initio green-kubo simulations of heat transport in solids: Method and implementation, Phys. Rev. B 107, 224304 (2023b).
  * Mandal  _et al._ [2022] S. Mandal, I. Maity, A. Das, M. Jain, and P. K. Maiti, Tunable lattice thermal conductivity of twisted bilayer mos2, Phys. Chem. Chem. Phys. 24, 13860 (2022).
  * Hellman  _et al._ [2011] O. Hellman, I. A. Abrikosov, and S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B 84, 180301 (2011).
  * Hellman and Abrikosov [2013] O. Hellman and I. A. Abrikosov, Temperature-dependent effective third-order interatomic force constants from first principles, Phys. Rev. B 88, 144301 (2013).
  * Bottin  _et al._ [2020] F. Bottin, J. Bieder, and J. Bouchet, a-tdep: Temperature dependent effective potential for abinit – lattice dynamic properties including anharmonicity, Computer Physics Communications 254, 107301 (2020).
  * Jain [2020] A. Jain, Multichannel thermal transport in crystalline ${\mathrm{tl}}_{3}{\mathrm{vse}}_{4}$, Phys. Rev. B 102, 201201 (2020).
  * Tadano  _et al._ [2014] T. Tadano, Y. Gohda, and S. Tsuneyuki, Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations, Journal of Physics: Condensed Matter 26, 225402 (2014).
  * Mahan [2000] G. D. Mahan,  _Many-Particle Physics_ (Springer New York, NY, 2000).
  * Anderson  _et al._ [1999] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK Users’ Guide_ , 3rd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999).
  * Togo and Tanaka [2018] A. Togo and I. Tanaka, Spglib: a software library for crystal symmetry search 10.48550/arXiv.1808.01590 (2018), arXiv:1808.01590 [cond-mat.mtrl-sci] .
  * Wang  _et al._ [2010] Y. Wang, J. J. Wang, W. Y. Wang, Z. G. Mei, S. L. Shang, L. Q. Chen, and Z. K. Liu, A mixed-space approach to first-principles calculations of phonon frequencies for polar materials, Journal of Physics: Condensed Matter 22, 202201 (2010).
  * Gonze  _et al._ [1994] X. Gonze, J.-C. Charlier, D. Allan, and M. Teter, Interatomic force constants from first principles: The case of $\alpha$-quartz, Phys. Rev. B 50, 13035 (1994).
  * Gonze and Lee [1997] X. Gonze and C. Lee, Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Phys. Rev. B 55, 10355 (1997).
  * West and Estreicher [2006] D. West and S. K. Estreicher, First-principles calculations of vibrational lifetimes and decay channels: Hydrogen-related modes in si, Phys. Rev. Lett. 96, 115504 (2006).
  * Shulumba  _et al._ [2017] N. Shulumba, O. Hellman, and A. J. Minnich, Intrinsic localized mode and low thermal conductivity of pbse, Phys. Rev. B 95, 014302 (2017).
  * Ravichandran and Broido [2018] N. K. Ravichandran and D. Broido, Unified first-principles theory of thermal properties of insulators, Phys. Rev. B 98, 085205 (2018).
  * Soler  _et al._ [2002] J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, The siesta method for ab initio order-n materials simulation, Journal of Physics: Condensed Matter 14, 2745 (2002).
  * Pal  _et al._ [2021] K. Pal, Y. Xia, and C. Wolverton, Microscopic mechanism of unusual lattice thermal transport in tlinte2, npj Computational Materials 7, 5 (2021).
  * Tamura [1983] S.-i. Tamura, Isotope scattering of dispersive phonons in ge, Phys. Rev. B 27, 858 (1983).
  * Ziman [2001] J. Ziman,  _Electrons and Phonons: The Theory of Transport Phenomena in Solids_ (Oxford University Press, 2001).
  * Broido  _et al._ [2005] D. A. Broido, A. Ward, and N. Mingo, Lattice thermal conductivity of silicon from empirical interatomic potentials, Phys. Rev. B 72, 014308 (2005).
  * Lindsay and Broido [2008] L. Lindsay and D. A. Broido, Three-phonon phase space and lattice thermal conductivity in semiconductors, Journal of Physics: Condensed Matter 20, 165209 (2008).
  * Omini and Sparavigna [1995] M. Omini and A. Sparavigna, An iterative approach to the phonon boltzmann equation in the theory of thermal conductivity, Physica B: Condensed Matter 212, 101 (1995).
  * MacDonald  _et al._ [1979] A. H. MacDonald, S. H. Vosko, and P. T. Coleridge, Extensions of the tetrahedron method for evaluating spectral properties of solids, Journal of Physics C: Solid State Physics 12, 2991 (1979).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Sanjurjo  _et al._ [1983] J. A. Sanjurjo, E. López-Cruz, P. Vogl, and M. Cardona, Dependence on volume of the phonon frequencies and the ir effective charges of several iii-v semiconductors, Phys. Rev. B 28, 4579 (1983).
  * Reich  _et al._ [2005] S. Reich, A. C. Ferrari, R. Arenal, A. Loiseau, I. Bello, and J. Robertson, Resonant raman scattering in cubic and hexagonal boron nitride, Phys. Rev. B 71, 205201 (2005).
  * Kikkawa  _et al._ [2021] J. Kikkawa, T. Taniguchi, and K. Kimoto, Nanometric phonon spectroscopy for diamond and cubic boron nitride, Phys. Rev. B 104, L201402 (2021).
  * Chen  _et al._ [2020] K. Chen, B. Song, N. K. Ravichandran, Q. Zheng, X. Chen, H. Lee, H. Sun, S. Li, G. A. G. U. Gamage, F. Tian, Z. Ding, Q. Song, A. Rai, H. Wu, P. Koirala, A. J. Schmidt, K. Watanabe, B. Lv, Z. Ren, L. Shi, D. G. Cahill, T. Taniguchi, D. Broido, and G. Chen, Ultrahigh thermal conductivity in isotope-enriched cubic boron nitride, Science 367, 555 (2020), https://www.science.org/doi/pdf/10.1126/science.aaz6149 .
  * Kang  _et al._ [2018] J. S. Kang, M. Li, H. Wu, H. Nguyen, and Y. Hu, Experimental observation of high thermal conductivity in boron arsenide, Science 361, 575 (2018), https://www.science.org/doi/pdf/10.1126/science.aat5522 .
  * Serrano  _et al._ [2002] J. Serrano, J. Strempfer, M. Cardona, M. Schwoerer-Böhning, H. Requardt, M. Lorenzen, B. Stojetz, P. Pavone, and W. J. Choyke, Determination of the phonon dispersion of zinc blende (3C) silicon carbide by inelastic x-ray scattering, Applied Physics Letters 80, 4360 (2002), https://pubs.aip.org/aip/apl/article-pdf/80/23/4360/10192323/4360_1_online.pdf .
  * Morelli  _et al._ [2002] D. T. Morelli, J. P. Heremans, and G. A. Slack, Estimation of the isotope effect on the lattice thermal conductivity of group iv and group iii-v semiconductors, Phys. Rev. B 66, 195304 (2002).
  * Ivanova  _et al._ [2006] L. M. Ivanova, P. A. Aleksandrov, and K. D. Demakov, Thermoelectric properties of vapor-grown polycrystalline cubic sic, Inorganic Materials 42, 1205 (2006).
  * Raunio  _et al._ [1969] G. Raunio, L. Almqvist, and R. Stedman, Phonon dispersion relations in nacl, Phys. Rev. 178, 1496 (1969).
  * Cowley  _et al._ [1983] E. R. Cowley, S. Satija, and R. Youngblood, Line shapes of longitudinal-optic phonons in sodium chloride at 300 and 600 k, Phys. Rev. B 28, 993 (1983).
  * Håkansson and Andersson [1986] B. Håkansson and P. Andersson, Thermal conductivity and heat capacity of solid nacl and nai under pressure, Journal of Physics and Chemistry of Solids 47, 355 (1986).
  * Klein and Caldwell [2004] M. V. Klein and R. F. Caldwell, Low Temperature System for Thermal Conductivity Measurements, Review of Scientific Instruments 37, 1291 (2004), https://pubs.aip.org/aip/rsi/article-pdf/37/10/1291/8352456/1291_1_online.pdf .
  * McCarthy and Ballard [2004] K. A. McCarthy and S. S. Ballard, Thermal Conductivity of Eight Halide Crystals in the Temperature Range 220°K to 390°K, Journal of Applied Physics 31, 1410 (2004), https://pubs.aip.org/aip/jap/article-pdf/31/8/1410/7928241/1410_1_online.pdf .
  * Goetz and Cowen [1982] M. Goetz and J. Cowen, The thermal conductivity of silver iodide, Solid State Communications 41, 293 (1982).
  * Wang  _et al._ [2023] Y. Wang, Q. Gan, M. Hu, J. Li, L. Xie, and J. He, Anharmonic lattice dynamics and the origin of intrinsic ultralow thermal conductivity in agi materials, Phys. Rev. B 107, 064308 (2023).
  * Acharyya  _et al._ [2022] P. Acharyya, T. Ghosh, K. Pal, K. S. Rana, M. Dutta, D. Swain, M. Etter, A. Soni, U. V. Waghmare, and K. Biswas, Glassy thermal conductivity in cs3bi2i6cl3 single crystal, Nature Communications 13, 5053 (2022).
  * Acharyya  _et al._ [2023] P. Acharyya, K. Pal, A. Ahad, D. Sarkar, K. S. Rana, M. Dutta, A. Soni, U. V. Waghmare, and K. Biswas, Extended antibonding states and phonon localization induce ultralow thermal conductivity in low dimensional metal halide, Advanced Functional Materials 33, 2304607 (2023), https://onlinelibrary.wiley.com/doi/pdf/10.1002/adfm.202304607 .
  * Ouyang  _et al._ [2023] N. Ouyang, Z. Zeng, C. Wang, Q. Wang, and Y. Chen, Role of high-order lattice anharmonicity in the phonon thermal transport of silver halide $\mathrm{Ag}x(x=\mathrm{Cl},\mathrm{Br},\mathrm{I})$, Phys. Rev. B 108, 174302 (2023).
  * Hsieh [2021] W.-P. Hsieh, High-pressure thermal conductivity and compressional velocity of nacl in b1 and b2 phase, Scientific Reports 11, 21321 (2021).
  * Srivastava and Merchant [1973] K. Srivastava and H. Merchant, Thermal expansion of alkali halides above 300°k, Journal of Physics and Chemistry of Solids 34, 2069 (1973).
  * Ravichandran and Broido [2019] N. K. Ravichandran and D. Broido, Non-monotonic pressure dependence of the thermal conductivity of boron arsenide, Nature Communications 10, 827 (2019).
  * Lindsay  _et al._ [2013] L. Lindsay, D. A. Broido, and T. L. Reinecke, First-principles determination of ultrahigh thermal conductivity of boron arsenide: A competitor for diamond?, Phys. Rev. Lett. 111, 025901 (2013).
  * Ravichandran and Broido [2020] N. K. Ravichandran and D. Broido, Phonon-phonon interactions in strongly bonded solids: Selection rules and higher-order processes, Phys. Rev. X 10, 021063 (2020).
  * Protik  _et al._ [2022] N. H. Protik, C. Li, M. Pruneda, D. Broido, and P. Ordejón, The elphbolt ab initio solver for the coupled electron-phonon boltzmann transport equations, npj Computational Materials 8, 28 (2022).