# Mode-coupling theory of lattice dynamics for classical and quantum crystals

Aloïs Castellano  J. P. Alvarinhas Batista  Matthieu J. Verstraete  Nanomat
group, QMAT center, CESAM research unit and European Theoretical Spectroscopy
Facility, Université de Liège, allée du 6 août, 19, B-4000 Liège, Belgium

###### Abstract

The dynamical properties of nuclei, carried by the concept of phonon
quasiparticles (QP), are central to the field of condensed matter. While the
harmonic approximation can reproduce a number of properties observed in real
crystals, the inclusion of anharmonicity in lattice dynamics is essential to
accurately predict properties such as heat transport or thermal expansion. For
highly anharmonic systems, non perturbative approaches are needed, which
result in renormalized theories of lattice dynamics. In this article, we apply
the Mori-Zwanzig projector formalism to derive an exact generalized Langevin
equation describing the quantum dynamics of nuclei in a crystal. By projecting
this equation on quasiparticles in reciprocal space, and with results from
linear response theory, we obtain a formulation of vibrational spectra that
fully accounts for the anharmonicity. Using a mode-coupling approach, we
construct a systematic perturbative expansion in which each new order is built
to minimize the following ones. With a truncation to the lowest order, we show
how to obtain a set of self-consistent equations that can describe the
lineshapes of quasiparticles. The only inputs needed for the resulting set of
equations are the static Kubo correlation functions, which can be computed
using (fully quantum) path-integral molecular dynamics or approximated with
(classical or ab initio) molecular dynamics. We illustrate the theory with an
application on fcc 4He, an archetypal quantum crystal with very strong
anharmonicity.

††preprint: AIP/123-QED

##  I Introduction

In the condensed phases of matter, many physical properties are profoundly
related to the vibrations of nuclei. The understanding of thermodynamical,
dynamical or transport phenomena in crystals requires the inclusion of the
correlated motion of atoms around their equilibrium positions. To describe
such dynamics, the cornerstone of modern theories is the harmonic
approximation [1]. Using a second-order Taylor expansion of the Born-
Oppenheimer surface (BOS), this model describes the lattice vibrations as a
non-interacting quantum gas of quasiparticles called phonons. Despite its
simplicity, the harmonic approximation has been tremendously successful in
explaining a large number of phenomena observed in the solid state, including
excitation spectra, the phase stability of various materials, elastic
properties or even zero point motion [2, 3]. Consequently, it naturally became
one of the most important tools in the field of condensed matter.
Nevertheless, the early truncation of the BOS Taylor expansion can be a severe
limitation to the accuracy of the harmonic approximation [4, 3, 5]. In some
cases these inaccuracies are only quantitative, but many phenomena cannot be
described even qualitatively with a purely harmonic theory of lattice
dynamics. Indeed, the higher order terms of the Taylor expansion, which are
denoted as anharmonic, are responsible for a number of properties that cannot
be predicted at all with the non-interacting phonon picture. Two of the most
striking examples of these limitations are probably the dynamical
stabilisation of some structures [6], and the computation of the thermal
conductivity, which diverges in the harmonic approximation [2]. These
limitations are a good illustration of the intrinsic problem with the non-
interacting phonon picture. Indeed, within the harmonic approach, phonons are
quasiparticles with temperature-independent energies and infinite lifetimes.
Such features are incompatible with the frequency changes associated with a
dynamical stabilisation, or the finite lifetimes associated with a finite
thermal conductivity. In order to reconcile the temperature-dependent
frequencies with a harmonic approximation, the quasi-harmonic approximation is
often used [3, 7], in which the temperature evolution is included indirectly
through a volume dependence. Unfortunately, this approach is biased in the
type of anharmonicity it includes[8], is often insufficient to accurately
describe the frequency evolution, and is still unable, by construction, to
incorporate a linewidth [9].

Nevertheless, a quasiparticle picture still seems to be an adequate
representation of the finite temperature vibrational excitations for most
solids. As a matter of fact, in most vibrational spectroscopy measurements on
pure crystals one can observe reasonably sharp Lorentzian-like resonances
[10]. These resonances can be associated with vibrational quasiparticles, and
it is these excitations that are approximated with harmonic phonons. However,
in contrast with the harmonic picture, these quasiparticles have temperature-
dependent frequencies $\Omega_{\nu}(T)$, as well as a linewidth
$\Gamma_{\nu}(T)$, which can be interpreted as a manifestation of their finite
lifetimes.

In recent years, a number of frameworks have been developed to account for the
anharmonicity of materials. For example, it is nowadays common to use
perturbation theory [11, 12, 13] to describe the effects of finite temperature
on phonons within an ab initio setting [14, 15, 16, 17]. In this approach, the
anharmonic contributions to the BOS are responsible for phonon-phonon
interactions, which leads to temperature-dependent frequencies as well as the
finite lifetimes missing in the harmonic approximation. However, if
perturbation theory can improve the agreement between harmonic phonons and
more realistic vibrational QP, it is based on the assumption that the higher
order anharmonic contributions will be of progressively lower importance, the
(T = 0K) harmonic contribution being dominant. This hypothesis is often too
drastic, and perturbative approaches are consequently known to fail in some
cases [18, 19, 20, 21]. Moreover, a significant drawback of such approaches is
that they can not be applied at all in strongly anharmonic systems with
unstable harmonic-level phonons. These failures suggest that approaches beyond
standard perturbation theory are necessary to obtain a satisfactory
description of vibrational properties of materials.

Starting from this observation, a significant number of works have focused on
developing non-perturbative methods to treat anharmonicity without the
drawbacks of the perturbative expansion. Most of these approaches are
constructed on the idea of renormalized phonons and phonon-phonon
interactions. The foremost non-perturbative theory of anharmonicity is the
Self-Consistent Harmonic Approximation (SCHA) [22, 23, 24, 25, 26, 27, 28, 29,
30, 31], also known as self consistent phonons, first imagined by Born and
since then rederived several times [4]. The idea behind the SCHA is that
interactions between harmonic phonons renormalize their frequencies. Using the
Gibbs-Bogoliubov free energy inequality [32, 23, 30], a Dyson equation for the
Green function [27] or even a path-integral formalism [33], the SCHA can be
derived as a set of self-consistent equations which uses the density matrix of
an effective harmonic Hamiltonian. In essence, the renormalized phonons are
obtained from an average computed on the distribution associated with the
converged Hamiltonian. However, the quadratic form of the effective
Hamiltonian means that the average is always computed on a multivariate
Gaussian distribution, which can differ substantially from the true
distribution. For instance, the renormalized phonons obtained as a result of
the self-consistent equations still need a correction to approximate the QP
frequencies [32, 29].

Fortunately, there exists a clear way to lift the multivariate Gaussian
distribution approximation. Indeed, given a potential $V(\mathbf{R})$, which
can be computed using ab initio methods such as DFT or a classical model, it
is possible to sample the exact canonical distribution in the classical limit
for a given temperature $T$ with methods such as molecular dynamics (MD) or
Markov-chain Monte-Carlo. Consequently, it seems natural to use these methods
to obtain a non-perturbative description of lattice dynamics to go beyond the
effective harmonic distribution of the SCHA. For instance, it has been
proposed to use the displacements covariance tensor computed from MD to map it
on an effective harmonic model through effective interatomic force constants
(IFC) [34, 35, 36]. However, the prevailing MD-based lattice dynamics method
is probably the temperature dependent effective potential (TDEP) [37, 38, 39].
This method uses configurations extracted from MD to fit an effective harmonic
potential as well as third and sometimes fourth order anharmonic terms [39].
An important feature of this potential is that each new order is fitted on the
residual forces from the previous order [40]. The goal of this procedure is to
minimize the importance of high order contributions compared to the effective
harmonic one. By using perturbation theory on top of this effective anharmonic
Hamiltonian, TDEP has proven its ability to describe anharmonic properties
with a great number of applications on the free energy and transport of
materials [37, 38, 39, 41, 42, 43]. Unfortunately, the renormalized anharmonic
Hamiltonian at the heart of TDEP lacks a theoretical justification. This
feature hinders a clear vision of its limits and raises questions on the
justification for the use of perturbation theory with the effective
Hamiltonian.

For instance, an important question raised by this approach is the easily
overlooked problem of going from the classical description of MD to the
quantum description of perturbation theory. Indeed, the QP description of a
many-body system is intrinsically constructed on correlation functions.
However, while the definition of classical correlation functions is
unambiguous, this is not the case for a quantum system, where there exist an
infinite number of ways to define them [44]. Consequently, the classical to
quantum transformation of correlation functions is a non-trivial task that
needs strong justifications [45, 46, 47]. It would seem that to dispose of
this problem, a solution could be to use the true quantum distribution of
displacements using for example path-integral molecular dynamics.
Unfortunately, the freedom of choice in how to define quantum correlation
functions is still an obstacle to a clear vision of how to apply TDEP. Indeed,
from all the different types of quantum correlation functions, which one
should be used as input for the least-squares problem ? In a recent
publication, it has been proposed to use the standard correlation functions
[48]. However, this choice was justified with the goal of having a good
description of the BOS, which unfortunately doesn’t answer the remaining
questions on how to approximate the QP with the effective Hamiltonian.

A more satisfying answer comes from the recent method introduced by Morresi et
al [49, 50] which uses the static Kubo correlation function (KCF) as input to
obtain the QP. With simple manipulations, the authors were able to show that
QPs computed this way match the exact zero frequency of the Matsubara QP
Green’s function. While their QPs feature an infinite lifetime, this result is
a strong suggestion that the KCF might be an answer to the questions raised
previously. From this type of correlation functions, the missing step is how
to obtain the dynamical QP properties, but the use of the KCF points to the
field of linear response theory, in which established theoretical tools exist
to solve such problems.

In this work, we build a foundation for an exact non-perturbative theory of
lattice dynamics that can describe strongly anharmonic quantum systems. With
this objective, we derive a dynamical theory for solids based on the Mori-
Zwanzig projector formalism of linear response theory [51, 52, 53, 54, 55, 56,
57]. This approach projects the dynamics on the slow variables of the system,
resulting in a generalized Langevin equation (GLE) and thus splitting the
equation of motion for the KCF into a conservative and in a memory-dependent
dissipative part. After defining a projection on quasiparticle operators in
reciprocal space, we use linear response theory to obtain an exact formulation
of the spectral function. By using a mode-coupling theory [58, 59, 60, 61] to
approximate the memory kernel inherent to the Mori-Zwanzig formalism, we
derive a theory in which the renormalized interatomic force constants of all
orders appear naturally. This approximation allows us to obtain a set of self-
consistent equations to compute vibrational spectra, for which we also propose
a simplification based on a scattering approximation. Moreover, we show how
our findings can justify a theory of temperature dependent phonons, thus
giving an explanation of the success of renormalized phonon methods. We also
discuss the implications of the method, showing for instance how it can be
used to approximate nuclear quantum effects on classical simulations, while
alleviating ambiguities related to the multiple definitions of quantum
correlation functions. We show that well established methods such as
perturbation theory or the SCHA can be recovered as approximations of our
mode-coupling theory, and that TDEP, used with (PI)MD, is a direct application
of our approach. Another interesting result of the present theory is that it
shows that the often-used Lorentzian approximation of spectral lineshapes can
never be an exact representation.

The paper is organized as follows. In section II, we introduce the notation
used throughout the paper as well as the Kubo correlation functions (KCF),
before deriving the generalized Langevin equation for the displacement KCF,
using the Mori-Zwanzig projector formalism. We then construct the projection
in reciprocal space in section III, before using the GLE to derive exact
equations describing the spectral properties in section IV. Using a mode-
coupling approach, we then derive in section V a set of approximations to the
memory kernel appearing in the GLE. Next, we discuss how to compute from ab
initio calculations the quantities needed as input for the method in section
VI. Section VII is devoted to the discussion of the classical limit as well as
the comparison to other theories of anharmonicity. We illustrate the theory
with an application on 4He in section VIII before concluding in section IX.

##  II Derivation of the Generalized Langevin Equation

###  II.1 Notations

In this work, we consider a three-dimensional crystal in the Born-Oppenheimer
approximation, so that the Hamiltonian governing the dynamics of the system is
taken with the general form

$H=\sum_{i}\frac{p_{i}^{2}}{2M_{i}}+V(\mathbf{R})$ (1)

where $\mathbf{R}$ and $\mathbf{p}$ are respectively the position and momentum
operators, and $V(\mathbf{R})$ is the potential. In this work, we will
consider that there is no atomic diffusion in the system, so that we can use
the displacements $\mathbf{u}_{i}=\mathbf{R}_{i}-\braket{\mathbf{R}_{i}}$ as
dynamical variables. In this definition, $\braket{O(t)}$ is the thermal
average of the operator $O$, meaning that $\braket{\mathbf{R}}$ are the
equilibrium positions. It should be noted that with this definition, the
equilibrium positions, and consequently the definitions of displacements, are
temperature dependent. While it is common to expand the potential in eq.(1) in
a Taylor expansion, here we make no assumptions on its form, apart from its
time independence.

To describe the time evolution of the system, we use the Heisenberg
representation, where the time dependence and derivatives of operators are
obtained with

$\begin{split}A(t)=&e^{\frac{i}{\hbar}Ht}Ae^{-\frac{i}{\hbar}Ht}=e^{i\mathcal{L}t}A\\\
\dot{A}(t)=&\frac{i}{\hbar}\big{[}H,A(t)\big{]}=i\mathcal{L}A(t)\end{split}$
(2)

where $A=A(0)$ and $\mathcal{L}=\frac{1}{\hbar}\big{[}H,\cdot]$ is the quantum
Liouville operator. We will denote $p_{i}(t)=M_{i}\dot{u}_{i}(t)$ the time
derivative of the displacement of atom $i$ with mass $M_{i}$. Furthermore, we
will define mass scaled coordinates $\tilde{u}_{i}=\sqrt{M_{i}}u_{i}$ and
$\tilde{p}_{i}=\frac{1}{\sqrt{M_{i}}}p_{i}$.

###  II.2 Kubo correlation function

Given the many-body nature of atomic vibrations in a crystal, a natural tool
to describe the dynamics are correlation functions. In this work, we will be
using the Kubo correlation function (KCF), which, for two operators $A$ and
$B$, is defined as

$\begin{split}G_{AB}(t)&=\mathrm{k}_{\mathrm{B}}\mathrm{T}\int_{0}^{\beta}d\lambda\braket{A^{\dagger}(i\lambda\hbar)B(t)}\\\
&=\frac{\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\mathcal{Z}}\int_{0}^{\beta}d\lambda\mathrm{Tr}\bigg{[}e^{-(\beta-\lambda)H}A^{\dagger}e^{-\lambda
H}B(t)\bigg{]}\\\ &=\big{(}A,B(t)\big{)}\end{split}$ (3)

where $\beta=\frac{1}{\mathrm{k}_{\mathrm{B}}\mathrm{T}}$ is the inverse
temperature and $\mathcal{Z}$ is the partition function. This correlation
function, which can be derived from linear response theory [53, 62], has the
advantage of being easily approximated in the path-integral molecular dynamics
formalism [63, 49]. Moreover, the KCF has a number of properties which
simplify its manipulation. First, similarly to classical correlation
functions, it is a real and even function of time and it obeys the relations
$G_{AB}(t)=G_{AB}(-t)$ as well as $G_{AB}(t)=G_{BA}(t)$. Using the Lehmann
representation, it can also be shown that derivatives of the static Kubo
correlation function are related [64]

$\displaystyle\frac{d^{k}}{dt^{k}}\big{(}A,B(t)\big{)}|_{t=0}=$
$\displaystyle(-1)^{k}\frac{d^{k}}{dt^{k}}\big{(}B,A(t)\big{)}|_{t=0}$ (4)
$\displaystyle\frac{d^{k}}{dt^{k}}\big{(}A,A(t)\big{)}|_{t=0}=$ $\displaystyle
0\leavevmode\nobreak\ \leavevmode\nobreak\ \forall k=1,3,5,\ldots$ (5)

Naturally, all the different quantum correlation functions are related by
simple relations [53, 47], so that knowledge of any type of correlation
function allows to obtain any other one. Finally, we will note that the KCF,
as the standard or any other quantum correlation function, reduces to the
classical correlation function in the limit $\hbar\rightarrow 0$.

In this work, the KCF of interest is the mass-weighted displacement KCF, which
is a $3N_{\mathrm{at}}\times 3N_{\mathrm{at}}$ symmetric matrix defined with

$\begin{split}G_{ij}(t)=&\sqrt{M_{i}M_{j}}\frac{\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\mathcal{Z}}\int_{0}^{\beta}d\lambda\mathrm{Tr}\bigg{[}e^{-(\beta-\lambda)H}u_{i}^{\dagger}e^{-\lambda
H}u_{j}(t)\bigg{]}\\\ =&\sqrt{M_{i}M_{j}}\big{(}u_{i},u_{j}(t)\big{)}\\\
=&\big{(}\tilde{u}_{i},\tilde{u}_{j}(t)\big{)}\end{split}$ (6)

with an equation of motion given by

$\begin{split}\ddot{G}_{ij}(t)=&\sqrt{M_{i}M_{j}}\big{(}u_{i},\ddot{u}_{j}(t)\big{)}\\\
=&\big{(}\tilde{u}_{i},\dot{\tilde{p}}_{j}(t)\big{)}\end{split}$ (7)

Our goal in this work is to derive a solution for this equation of motion.

###  II.3 The Generalized Langevin Equation

Presented in the form of eq.(7), solving for the dynamics of the many-body
system is an intractable problem. Fortunately, the projection operator
formalism derived by Mori [51] and Zwanzig [52] allows to express eq.(7) in a
form amenable to approximations. This formalism is based on operators
$\mathcal{P}$ and $\mathcal{Q}$ which project an observable $O(t)$ onto the
subspace of “slow" dynamical variables of interest. In this work, we will use
the mass scaled coordinates as “slow" variables, so that the projectors are
defined as

$\displaystyle\mathcal{P}O(t)=$
$\displaystyle\sum_{ij}\frac{\big{(}\tilde{u}_{j},O(t)\big{)}}{\big{(}\tilde{u}_{i},\tilde{u}_{j}\big{)}}\tilde{u}_{i}+\sum_{ij}\frac{\big{(}\tilde{p}_{j},O(t)\big{)}}{\big{(}\tilde{p}_{i},\tilde{p}_{j}\big{)}}\tilde{p}_{i}$
(8) $\displaystyle\mathcal{Q}O(t)=$ $\displaystyle(1-\mathcal{P})O(t)$ (9)

By construction, it should be noted that the projectors $\mathcal{P}$ and
$\mathcal{Q}$ are orthogonal, in the sense that they follow
$\mathcal{P}+\mathcal{Q}=1$ and
$\mathcal{P}\mathcal{Q}=\mathcal{Q}\mathcal{P}=0$. Being projection operators,
it is natural that $\mathcal{P}$ and $\mathcal{Q}$ are idempotent, so that
$\mathcal{P}\mathcal{P}=\mathcal{P}$ and $\mathcal{Q}\mathcal{Q}=\mathcal{Q}$.
Moreover, it can be shown that these projectors, as the Liouville operators,
are hermitian with respect to the static Kubo correlation function [55], which
results in the properties

$\displaystyle\big{(}X,\mathcal{P}Y\big{)}=$
$\displaystyle\big{(}\mathcal{P}X,Y\big{)}$ (10)
$\displaystyle\big{(}X,\mathcal{Q}Y\big{)}=$
$\displaystyle\big{(}\mathcal{Q}X,Y\big{)}$ (11)

Figure 1: Schematic illustration of the Mori-Zwanzig projection operator
formalism. Within this formalism, the dynamics of a dynamical variable,
represented here with a black arrow, happen in the full Hilbert space. The
operators $\mathcal{P}$, represented as the red arrows, projects this dynamic
onto a “slow" subspace depicted here as a the blue rectangle. By construction,
the starting point of the dynamics, visualized with the green dot, is located
on the projected subspace.

Equipped with these operators, we can now express the equation of motion of
the mass scaled momentum of atom $i$

$\begin{split}\dot{\tilde{p}}_{i}(t)=&e^{i\mathcal{L}t}i\mathcal{L}\tilde{p}_{i}\\\
=&e^{i\mathcal{L}t}(\mathcal{P}+\mathcal{Q})i\mathcal{L}\tilde{p}_{i}\\\
=&e^{i\mathcal{L}t}\mathcal{P}i\mathcal{L}\tilde{p}_{i}+e^{i\mathcal{Q}\mathcal{L}t}\mathcal{Q}i\mathcal{L}\tilde{p}_{i}\\\
&+\int_{0}^{t}dse^{i\mathcal{L}(t-s)}\mathcal{P}i\mathcal{L}e^{i\mathcal{Q}\mathcal{L}s}\mathcal{Q}i\mathcal{L}\tilde{p}_{i}\end{split}$
(12)

where we used the Dyson identity [59, 55] to expand the $\mathcal{Q}$
projection.

To simplify the notation of eq.(12), we introduce the frequency matrix
$\mathbf{\Theta}$ (note that it is not in mass scaled coordinates)

$\Theta_{ij}=-\sum_{k}\frac{\big{(}u_{k},f_{i}\big{)}}{\big{(}u_{j},u_{k}\big{)}}$
(13)

where $f_{i}=M_{i}\mathcal{L}^{2}u_{i}$ is the force acting on atom $i$, as
well as $\delta\mathbf{f}(t)$, which is often called the random force,

$\delta
f_{i}(t)=e^{i\mathcal{Q}\mathcal{L}t}\big{(}f_{i}+\sum_{j}\Theta_{ij}u_{j}\big{)}$
(14)

and the memory kernel $\mathbf{K}(t)$

$K_{ij}(t)=\frac{\big{(}\delta f_{i},\delta
f_{j}(t)\big{)}}{\sqrt{M_{i}M_{j}}}$ (15)

so that the equations of motion for the momentum can be written as [58, 59,
60, 61]

$\begin{split}\dot{\tilde{p}}_{i}(t)=&-\sum_{j}\frac{\Theta_{ij}}{\sqrt{M_{i}M_{j}}}\tilde{u}_{j}(t)-\frac{\delta
f_{i}(t)}{\sqrt{M_{i}}}\\\
&-\beta\sum_{j}\int_{0}^{t}dsK_{ij}(s)\tilde{p}_{j}(t-s)\end{split}$ (16)

Formally, the Mori-Zwanzig projection splits the evolution of the dynamical
variables into a “slow" contribution, driven by the frequency matrix and a
“fast" one, given by the random force and the memory kernel. As depicted in
Fig.1, the “slow" part of the dynamics represents a subspace of the full
Hilbert space of the dynamical variables. Thus, the goal of the operator
$\mathcal{P}$ is to project any dynamical variable into this subspace. The
random forces, which account for the effects of other degrees of freedom,
evolve in the “fast" orthogonal subspace, as indicated by the evolution
operator $e^{i\mathcal{Q}\mathcal{L}t}$, and its orthogonal nature is
recovered thanks to the property $\big{(}\tilde{u}_{i},\delta
f_{i}(t)\big{)}=0$. It should be noted that the randomness of this term is
only apparent. By contracting information through the projector $\mathcal{Q}$,
the time evolution of this term is more difficult to apprehend and its
dynamics appears “random" [55], though this is not the case.

We can now inject $\dot{\tilde{p}}_{i}(t)$ into the equations of motion of
eq.(7), to obtain the generalized Langevin equation of motion of the mass
weighted displacement KCF

$\ddot{\mathbf{G}}(t)=-\frac{\mathbf{\Theta}}{\sqrt{\mathbf{M}^{T}\mathbf{M}}}\mathbf{G}(t)-\beta\int_{0}^{t}ds\mathbf{K}(s)\dot{\mathbf{G}}(t-s)$
(17)

An interesting result of this formulation is the separation of the dynamics
into two parts. The first one, driven by the static frequency matrix is
conservative, while the second part, driven by the memory kernel is
dissipative [61]. The formulation of the memory kernel as the autocorrelation
of the random forces is a result of the equilibrium state of the system and is
known as the second fluctuation-dissipation theorem [53]. It is important to
note that there is no approximation leading to eq.(17), so that it is an exact
description of the dynamics of the system. However, the main difficulty in
solving this equation comes from the memory kernel of eq.(15), for which the
dynamics is not generated by the Liouville operator, but by its orthogonal
projection.

##  III The GLE in QP space

Up to this point, we have been working in real space. To formulate the problem
in reciprocal space, there is a freedom of choice for the basis set into which
we will project the dynamical variables. A natural projection is one where the
conservative part of the GLE is block diagonal in momentum and completely
diagonal in vibrational modes, all while taking advantage of the periodicity
of the crystal. Such a transformation can be achieved using concepts from the
toolkit of harmonic phonon theory.

With this goal, we define a generalized dynamical matrix
$\mathbf{D}(\mathbf{q})$ as the Fourier transform of the mass weighted
$\mathbf{\Theta}$

$\mathbf{D}_{ij}(\mathbf{q})=\sum_{\tau}\frac{\Theta_{0,i}^{\tau,j}}{\sqrt{M_{i}M_{j}}}e^{i\mathbf{q}(\braket{\mathbf{R}_{\tau,j}}-\braket{\mathbf{R}_{0,i}})}$
(18)

where we introduced a unit cell based notation $i\rightarrow(\tau,i)$, so that
$\mathbf{\Theta}_{\delta,i}^{\tau,j}$ is the frequency matrix relating atom
$i$ in the unit cell $\delta$ to the atom $j$ in unit cell $\tau$. Similarly
to the harmonic phonon case, we can now obtain QP polarization vectors
$\pmb{\epsilon}^{\nu}(\mathbf{q})$ and frequencies $\Omega_{\mathbf{q},\nu}$
by diagonalizing $\mathbf{D}(\mathbf{q})$

$\mathbf{D}(\mathbf{q})\pmb{\epsilon}^{\nu}(\mathbf{q})=\Omega_{\nu}^{2}(\mathbf{q})\pmb{\epsilon}^{\nu}(\mathbf{q})$
(19)

where $\nu$ is a QP mode. Without loss of generality, we can now express the
mass-weighted atomic displacements as

$\tilde{u}_{i}(t)=\sum_{\mathbf{q},\lambda}\rho_{\mathbf{q},\lambda}\epsilon_{i}^{\lambda}(\mathbf{q})A_{\lambda}(\mathbf{q},t)$
(20)

where $A_{\lambda}(\mathbf{q},t)$ are unitless quasiparticle operators with
the property
$\big{(}A_{\lambda}(\mathbf{q}),A_{\lambda^{\prime}}(\mathbf{q}^{\prime})\big{)}=\frac{2\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\hbar\Omega_{\lambda}(\mathbf{q})}\delta_{\lambda\lambda^{\prime}}\delta_{\mathbf{q},\mathbf{q}^{\prime}}$111While
the choice of prefactor here can seem arbitrary, we use it to simplify
comparison with other methods of anharmonic lattice dynamics and
$\rho_{\mathbf{q},\lambda}$ is a prefactor determined below. Using the static
solution of eq.(17) as well as the property
$\ddot{G}_{ij}(0)=-\mathrm{k}_{\mathrm{B}}\mathrm{T}/M_{i}\delta_{ij}$ [66,
67], we can express the static mass-weighted displacement covariance tensor as

$\big{(}\tilde{u}_{i},\tilde{u}_{j}\big{)}=\mathrm{k}_{\mathrm{B}}\mathrm{T}\sqrt{M_{i}M_{j}}\Theta^{-1}_{ij}$
(21)

Injecting eq.(20) into this equation and with some reorganisation, it can be
shown that
$\rho_{\mathbf{q},\lambda}=\sqrt{\hbar/\Omega_{\lambda}(\mathbf{q})}$, so that
the displacements can be expressed as

$u_{i}(t)=\sqrt{\frac{\hbar}{2M_{i}}}\sum_{\mathbf{q},\lambda}\frac{\epsilon_{i}^{\lambda}(\mathbf{q})}{\sqrt{\Omega_{\lambda}(\mathbf{q})}}A_{\lambda}(\mathbf{q},t)$
(22)

We can now define a Kubo transformed QP Green’s function

$\mathbf{G}_{\lambda,\mu}(\mathbf{q},t)=\big{(}A_{\lambda}(\mathbf{q}),A_{\mu}(\mathbf{q},t)\big{)}$
(23)

To obtain the expression of the memory kernel projected on the QP, we note
that from eq.(22), the force acting on nucleus $i$ can be expressed as

$f_{i}=\sqrt{\frac{\hbar
M_{i}}{2}}\sum_{\mathbf{q}\lambda}\frac{\epsilon_{i}^{\lambda}(\mathbf{q})}{\sqrt{\Omega_{\lambda}(\mathbf{q})}}\ddot{A}_{\lambda}(\mathbf{q})$
(24)

The part of the random forces involving the frequency matrix has a structure
resembling mass-scaled harmonic forces, which can be written in terms of
quasiparticle operators as

$\sum_{j\beta}\frac{\Theta_{ij}^{\alpha\beta}}{\sqrt{M_{i}M_{j}}}u_{j}^{\beta}=\sqrt{\frac{\hbar
M_{i}}{2}}\sum_{\mathbf{q}\lambda}\frac{\epsilon_{i}^{\lambda}(\mathbf{q})}{\sqrt{\Omega_{\lambda}(\mathbf{q})}}\Omega_{\lambda}^{2}(\mathbf{q})A_{\lambda}(\mathbf{q})$
(25)

Consequently, we can define the random forces in terms of QP operators as

$\begin{split}\delta f_{i}=&\sqrt{\frac{\hbar
M_{i}}{2}}\sum_{\mathbf{q},\lambda}\frac{\epsilon_{i}^{\lambda}(\mathbf{q})}{\sqrt{\Omega_{\lambda}(\mathbf{q})}}\bigg{(}\ddot{A}_{\lambda}(\mathbf{q})-\Omega_{\lambda}^{2}(\mathbf{q})A_{\lambda}(\mathbf{q})\bigg{)}\\\
=&\sqrt{\frac{\hbar
M_{i}}{2}}\sum_{\mathbf{q},\lambda}\frac{\epsilon_{i}^{\lambda}(\mathbf{q})}{\sqrt{\Omega_{\lambda}(\mathbf{q})}}\delta
A_{\lambda}(\mathbf{q})\end{split}$ (26)

where we defined the QP random forces operator $\delta
A_{\lambda}(\mathbf{q})=\ddot{A}_{\lambda}(\mathbf{q})-\Omega_{\lambda}^{2}(\mathbf{q})A_{\lambda}(\mathbf{q})$,
which follows orthogonal time dynamics

$\delta A_{\lambda}(\mathbf{q},t)=e^{i\mathcal{Q}\mathcal{L}t}\delta
A_{\lambda}(\mathbf{q})$ (27)

We are now able to define the memory kernel in reciprocal space

$\mathbf{K}(\mathbf{q},t)=\frac{\hbar}{\mathrm{k}_{\mathrm{B}}\mathrm{T}}\mathbf{\Omega}(\mathbf{q})^{-1}\big{(}\delta\mathbf{A}(\mathbf{q}),\delta\mathbf{A}(\mathbf{q},t)\big{)}$
(28)

which allows us to write the GLE for the Kubo transformed Green’s function as

$\ddot{\mathbf{G}}(\mathbf{q},t)=-\mathbf{\Omega}^{2}(\mathbf{q})\mathbf{G}(\mathbf{q},t)-\int_{0}^{t}ds\mathbf{K}(\mathbf{q},s)\dot{\mathbf{G}}(\mathbf{q},t-s)$
(29)

The dynamical properties of the many-body system can now be obtained by
solving eq.(29) with the initial conditions

$\displaystyle G_{\lambda\mu}(\mathbf{q},0)=$
$\displaystyle\frac{2\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\hbar\Omega_{\lambda}(\mathbf{q})}\delta_{\lambda,\mu}$
(30) $\displaystyle\dot{G}_{\lambda\mu}(\mathbf{q},0)=$ $\displaystyle 0$ (31)

##  IV Vibrational spectra from the Kubo correlation function

From linear response theory, it is known that the correlations of an
unperturbed system are related to its response to a perturbation [53, 62, 68].
For a perturbation applied to an observable $A$, the response of the system
through the observable $B$ can be computed using the response function
$\Phi_{AB}$, that we define as

$\Phi_{AB}(t)=-\frac{i}{\hbar}\braket{[A,B(t)]}$ (32)

where the average of the commutator of $A$ and $B(t)$ is called the
symmetrized correlation function. The dynamical susceptibility of the system
is then simply given by the Laplace transform of the response function

$\raisebox{2.0pt}{\mbox{\large$\chi$}}_{AB}(\omega)=\int_{0}^{\infty}dt\Phi_{AB}(t)e^{-i\omega
t}$ (33)

It should be noted that the use of a Laplace transform with the definition
eq.(32) allows to include causality in the response of the system. This
formulation is equivalent to more common definitions of the susceptibility
where causality is introduced with the use of the Heaviside function. To
connect the preceding derivation to the response function, we need to express
$\Phi_{AB}(t)$ as a function of the KCF. This can be achieved using the cyclic
property of the trace, which gives a classical-like formulation of the
response function

$\begin{split}\Phi_{AB}(t)=&-\frac{i}{\hbar}\frac{1}{\mathcal{Z}}\bigg{(}\mathrm{Tr}\big{[}e^{-\beta
H}AB(t)]\big{]}-\mathrm{Tr}\big{[}e^{-\beta H}B(t)A]\big{]}\bigg{)}\\\
=&\frac{i}{\hbar}\frac{1}{\mathcal{Z}}\bigg{(}\mathrm{Tr}\big{[}e^{-\beta
H}e^{\beta H}Ae^{-\beta H}B(t)]\big{]}\\\ &-\mathrm{Tr}\big{[}e^{-\beta
H}AB(t)]\big{]}\bigg{)}\\\
=&\int_{0}^{\beta}d\lambda\frac{1}{\mathcal{Z}}\mathrm{Tr}\big{[}e^{-\beta
H}\dot{A}(i\hbar\lambda)B(t)]\\\ =&\beta\dot{G}_{AB}(t)\end{split}$ (34)

where $\dot{A}=i[H,A]/\hbar$ is the time derivative of $A$ at $t=0$.

From this response function, the dynamical susceptibility of the system is
simply obtained by applying a Laplace transform

$\begin{split}\raisebox{2.0pt}{\mbox{\large$\chi$}}_{AB}(\omega)=&\beta\int_{0}^{\infty}dt\big{(}\dot{A},B(t)\big{)}e^{-i\omega
t}\\\ =&\beta\widetilde{\dot{G}}_{AB}(i\omega)\\\
=&i\beta\omega\widetilde{G}_{AB}(i\omega)-\beta G_{AB}(0)\end{split}$ (35)

with $\widetilde{G}_{AB}(i\omega)$ the Laplace transform of the KCF of $A$ and
$B$, and $\widetilde{\dot{G}}_{AB}(i\omega)$ is the Laplace transform of its
time derivative. The dynamical susceptibility can be decomposed into a real
and imaginary part
$\raisebox{2.0pt}{\mbox{\large$\chi$}}_{AB}(\omega)=\raisebox{2.0pt}{\mbox{\large$\chi$}}^{\prime}_{AB}(\omega)+i\raisebox{2.0pt}{\mbox{\large$\chi$}}^{\prime\prime}_{AB}(\omega)$,
where $\raisebox{2.0pt}{\mbox{\large$\chi$}}_{AB}^{\prime}(\omega)$ and
$\raisebox{2.0pt}{\mbox{\large$\chi$}}_{AB}^{\prime\prime}(\omega)$ are
related by a Kramers-Kronig transform, meaning that the knowledge of one is
sufficient to recover the dynamical response of the system. From eq. (35), we
see that the imaginary part of the dynamical susceptibility, which corresponds
to the dissipative part of the response, is simply proportional to the real
part of $\widetilde{G}_{AB}(\omega)$. Since the KCF is an even and real
function of time, this real part is proportional to the power spectrum
$G_{AB}(\omega)$, which corresponds to the Fourier transform of $G_{AB}(t)$,
so that the dissipative part of the response, also called the spectral
function, can be written as

$\raisebox{2.0pt}{\mbox{\large$\chi$}}_{AB}^{\prime\prime}(\omega)=\frac{\omega}{2\mathrm{k}_{\mathrm{B}}\mathrm{T}}G_{AB}(\omega)$
(36)

This equation is simply the fluctuation-dissipation theorem, written using the
Kubo transformed correlation function. The trace of the dissipation
$\mathbf{\raisebox{2.0pt}{\mbox{\large$\chi$}}}^{\prime\prime}(\mathbf{q},\omega)$
is particularly important since this quantity can be linked to the cross
section $d^{2}\sigma/d\Omega dE$ from spectroscopy measurements, allowing thus
a direct comparison between theory and experiment. For the QP KCF, using the
initial conditions eq. (30) and (31), the Laplace transform is expressed as

$\widetilde{\mathbf{G}}(\mathbf{q},i\omega)=\frac{2\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\hbar\mathbf{\Omega}(\mathbf{q})}\frac{-i\omega-\tilde{\mathbf{K}}(\mathbf{q},i\omega)}{\omega^{2}-\mathbf{\Omega}^{2}(\mathbf{q})-i\omega\tilde{\mathbf{K}}(\mathbf{q},i\omega)}$
(37)

To obtain a formula for the Fourier transform of
$\mathbf{G}(\mathbf{q},\omega)$, we start by separating the memory kernel into
a real and an imaginary part as
$\tilde{\mathbf{K}}(\mathbf{q},i\omega)=\tilde{\mathbf{\Gamma}}(\mathbf{q},\omega)-i\tilde{\mathbf{\Delta}}(\mathbf{q},\omega)$.
Since $\mathbf{K}(\mathbf{q},t)$ is a real and even function of time, the real
part $\tilde{\mathbf{\Gamma}}(\mathbf{q},\omega)$ of its Laplace transform is
proportional to its Fourier transform $\mathbf{\Gamma}(\mathbf{q},\omega)$,
with
$\tilde{\mathbf{\Gamma}}(\mathbf{q},\omega)=2\mathbf{\Gamma}(\mathbf{q},\omega)$.
Moreover, $\tilde{\mathbf{\Gamma}}(\mathbf{q},\omega)$ and
$\tilde{\mathbf{\Delta}}(\mathbf{q},\omega)$ are not independent, and one can
recover the other using Kramers-Kronig relations [69]. For instance,
$\tilde{\mathbf{\Delta}}(\mathbf{q},\omega)$ can be computed with

$\tilde{\mathbf{\Delta}}(\mathbf{q},\omega)=\frac{1}{\pi}\int
d\omega^{\prime}\frac{\tilde{\mathbf{\Gamma}}(\mathbf{q},\omega^{\prime})}{\omega^{\prime}-\omega}$
(38)

where the absence of the usual minus sign comes from the previous definition
of the memory kernel. This allows us to take the real part of eq.(37) and
obtain the power spectrum as

$\mathbf{G}(\mathbf{q},\omega)=\frac{\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\pi\hbar}\frac{8\mathbf{\Omega}(\mathbf{q})\mathbf{\Gamma}(\mathbf{q},\omega)}{(\omega^{2}-\mathbf{\Omega}^{2}(\mathbf{q})-2\omega\mathbf{\Delta}(\mathbf{q},\omega))^{2}+4\omega^{2}\mathbf{\Gamma}^{2}(\mathbf{q},\omega)}$
(39)

Finally, the dissipative part of the response function can be computed using
the fluctuation-dissipation theorem of eq.(36), giving

$\mathbf{\raisebox{2.0pt}{\mbox{\large$\chi$}}}^{\prime\prime}(\mathbf{q},\omega)=\frac{1}{\pi\hbar}\frac{4\omega\mathbf{\Omega}(\mathbf{q})\mathbf{\Gamma}(\mathbf{q},\omega)}{(\omega^{2}-\mathbf{\Omega}^{2}(\mathbf{q})-2\omega\mathbf{\Delta}(\mathbf{q},\omega))^{2}+4\omega^{2}\mathbf{\Gamma}^{2}(\mathbf{q},\omega)}$
(40)

In the end, the dynamical response of the system can be entirely characterised
using $\mathbf{\Omega}(\mathbf{q})$, $\mathbf{\Gamma}(\mathbf{q},\omega)$ and
$\mathbf{\Delta}(\mathbf{q},\omega)$. Consequently, at this point the only
missing ingredient needed to compute the dynamical properties of the system is
the real part of the memory kernel.

##  V Mode-coupling approximation for the memory kernel

For the moment, the formulation is exact and consists only in a rewriting of
the dynamical problem, by framing the unsolvable part into a memory kernel.
The exact computation of $\mathbf{K}(t)$ is a formidable task and we have to
resort to approximations to be able to apply the theory to real simulations.
To obtain a tractable formulation of the memory kernel, we use a cluster
expansion of the random forces. With this goal, we define the $n$-th order
projection operators $\mathcal{P}_{n}$ [70], which project an operator $O(t)$
on the subspace of $n$ QP operators, as well as the corresponding orthogonal
projection $\mathcal{Q}_{n}$

$\displaystyle\mathcal{P}_{n}O(t)=$
$\displaystyle\frac{1}{n!}\sum_{\begin{subarray}{c}\lambda_{1}\ldots\lambda_{n}\\\
\lambda^{\prime}_{1}\ldots\lambda^{\prime}_{n}\end{subarray}}\frac{\big{(}A_{\lambda^{\prime}_{1}}\ldots
A_{\lambda^{\prime}_{n}},O(t)\big{)}}{\big{(}A_{\lambda_{1}}\ldots
A_{\lambda_{n}},A_{\lambda^{\prime}_{1}}\ldots
A_{\lambda^{\prime}_{n}}\big{)}}A_{\lambda_{1}}\ldots A_{\lambda_{n}}$ (41)
$\displaystyle\mathcal{Q}_{n}=$ $\displaystyle 1-\mathcal{P}_{n}$ (42)

where the prefactor $1/n!$ is here to account for all the permutations of
$A_{\lambda}$ operators, so that
$\mathcal{P}_{n}\mathcal{P}_{n}=\mathcal{P}_{n}$. Going up to order $N$ and
inserting several times the identity $\mathcal{P}_{n}+\mathcal{Q}_{n}=1$, we
can rewrite the random forces

$\displaystyle\begin{split}\delta
A_{\lambda}=&(\mathcal{P}_{2}+\mathcal{Q}_{2})\delta A_{\lambda}\\\
=&\mathcal{P}_{2}\delta
A_{\lambda}+(\mathcal{P}_{3}+\mathcal{Q}_{3})\mathcal{Q}_{2}\delta
A_{\lambda}\\\
=&\sum_{n}^{N}\mathcal{P}_{n}\delta_{n-1}A_{\lambda}+\delta_{N}A_{\lambda}\end{split}$
(43) $\displaystyle\delta_{N}A_{\lambda}=$
$\displaystyle\mathcal{Q}_{N}\mathcal{Q}_{N-1}\ldots\mathcal{Q}_{2}\delta
A_{\lambda}$ (44)

where each successive application of $\mathcal{Q}_{n}$ will produce an ever
smaller $\delta_{n}A_{\lambda}$, and $\delta_{N}A_{\lambda}$ is the final
residual random force we will neglect. We can now inject this expression into
the memory kernel to obtain

$\begin{split}K^{\lambda,\lambda^{\prime}}(t)=&\big{(}\delta
A_{\lambda},\delta A_{\lambda^{\prime}}(t)\big{)}\\\
=&\sum_{n=2}^{N}\sum_{n^{\prime}=2}^{N}\big{(}\mathcal{P}_{n}\delta_{n-1}A_{\lambda},\mathcal{P}_{n^{\prime}}\delta_{n^{\prime}-1}A_{\lambda^{\prime}}(t)\big{)}\\\
&+\big{(}\delta_{N}A_{\lambda^{\prime}},\delta_{N}A_{\lambda^{\prime}}(t)\big{)}\\\
=&\sum_{n=2}^{N}\sum_{n^{\prime}=2}^{N}K_{nn^{\prime}}^{\lambda,\lambda^{\prime}}(t)+\big{(}\delta_{N}A_{\lambda},\delta_{N}A_{\lambda^{\prime}}(t)\big{)}\end{split}$
(45)

where we define the $(n,n^{\prime})$ order memory kernel

$\begin{split}K_{n,n^{\prime}}^{\lambda,\lambda^{\prime}}(t)=&\big{(}\mathcal{P}_{n}\delta_{n-1}A_{\lambda},\mathcal{P}_{n^{\prime}}\delta_{n^{\prime}-1}A_{\lambda^{\prime}}(t)\big{)}\\\
=&\sum_{\begin{subarray}{c}\lambda_{1},\ldots,\lambda_{n}\\\
\lambda^{\prime}_{1},\ldots,\lambda^{\prime}_{n^{\prime}}\end{subarray}}\Psi_{n}(\lambda,\lambda_{1},\ldots,\lambda_{n})\Psi_{n^{\prime}}(\lambda^{\prime},\lambda^{\prime}_{1},\ldots,\lambda^{\prime}_{n^{\prime}})\\\
&\times\big{(}A_{\lambda_{1}}\ldots
A_{\lambda_{n}},e^{i\mathcal{Q}\mathcal{L}t}A_{\lambda^{\prime}_{1}}\ldots
A_{\lambda^{\prime}_{n^{\prime}}}\big{)}\end{split}$ (46)

as well as the vertex

$\Psi_{n}(\lambda_{1},\ldots,\lambda_{n})=\frac{1}{n!}\sum_{\mu_{1},\ldots,\mu_{n}}\frac{\big{(}A_{\mu_{1}}\ldots
A_{\mu_{n}},\delta_{n-1}A_{\lambda}\big{)}}{\big{(}A_{\lambda_{1}}\ldots
A_{\lambda_{n}},A_{\mu_{1}}\ldots A_{\mu_{n}}\big{)}}$ (47)

The preceding manipulations serve to shift the orthogonal space time
dependence from the random forces into multiple displacement Kubo correlation
functions. By subtracting the contributions up to order $n$, the
$\mathcal{Q}_{n}$ projectors minimize the corresponding contribution of the
random forces $\delta_{n}A_{\lambda}$. Consequently, this “clustering"
procedure ensures that the contribution of each new order in the memory kernel
is less important than the previous one. While this is still impossible to
solve in practice, this rewriting opens a way to approximate
$\mathbf{K}(\mathbf{q},t)$.

As a first approximation, we truncate the expansion of the memory kernel at
some order $N$ and we neglect the $N$-th order correlation of the random
forces. The second approximation consists in decoupling the $n$ point
correlations in orthogonal time into $n/2$ products of $2$ point correlations
in real time [58, 59, 61]. For instance, a four point correlation function is
given by

$\begin{split}\big{(}A_{\lambda_{1}}A_{\lambda_{2}},e^{i\mathcal{Q}\mathcal{L}t}A_{\lambda_{3}}A_{\lambda_{4}}\big{)}\approx&\big{(}A_{\lambda_{1}},A_{\lambda_{3}}(t)\big{)}\big{(}A_{\lambda_{2}},A_{\lambda_{4}}(t)\big{)}\\\
&+\big{(}A_{\lambda_{1}},A_{\lambda_{4}}(t)\big{)}\big{(}A_{\lambda_{2}},A_{\lambda_{3}}(t)\big{)}\end{split}$
(48)

These two approximations are at the foundation of the mode-coupling theory of
the liquid-glass transition [58, 59].

In this work, we will limit ourselves to the second order in the memory kernel
expansion. To simplify further the expression of the memory kernel, we will
neglect its off-diagonal terms ($n\neq n^{\prime}$,
$\lambda\neq\lambda^{\prime}$). We note that this approximation does not
remove mode coupling (the sum in eq.(46) runs over all the other modes) and
would not prevent the inclusion of off-diagonal velocity terms in the thermal
conductivity as in Ref. [71]. Finally, taking into account the conservation of
crystal momentum (with the $-\lambda$), our approximation for the memory
kernel is given by

$K_{\lambda}(t)\approx\frac{\hbar}{2\mathrm{k}_{\mathrm{B}}\mathrm{T}\Omega_{\lambda}}\sum_{\mu,\nu}|\Psi(-\lambda,\mu,\nu)|^{2}G_{\mu}(t)G_{\nu}(t)$
(49)

The only inputs are contained in the vertex $\Psi$. We define the third-order
tensor $\Theta$ as

$\Theta_{ijk}=-\frac{1}{2}\frac{\sum_{l,m}\big{(}u_{l}u_{m},\delta
f_{i}\big{)}}{\sum_{l,m}\big{(}u_{j}u_{k},u_{l}u_{m}\big{)}}$ (50)

Using the projection of displacements and random forces on QP, the vertex can
then be expressed as

$\begin{split}\Psi(\lambda,\mu,\nu)=&\sqrt{\frac{\hbar}{2}\frac{\Omega_{\lambda}}{\Omega_{\mu}\Omega_{\nu}}}\sum_{i,j,k}\frac{\pmb{\epsilon}_{i}^{\lambda}\pmb{\epsilon}_{j}^{\mu}\pmb{\epsilon}_{k}^{\nu}}{\sqrt{M_{i}M_{j}M_{k}}}\Theta_{ijk}\\\
&\times
e^{i(\mathbf{R}_{i}\mathbf{q}_{\lambda}+\mathbf{R}_{j}\mathbf{q}_{\mu}+\mathbf{R}_{k}\mathbf{q}_{\nu})}\delta_{\mathbf{G}}(\mathbf{q}_{\lambda}+\mathbf{q}_{\mu}+\mathbf{q}_{\nu})\\\
\end{split}$ (51)

where
$\delta_{\mathbf{G}}(\mathbf{q}_{\lambda}+\mathbf{q}_{\mu}+\mathbf{q}_{\nu})=1$
if $\mathbf{q}_{\lambda}+\mathbf{q}_{\mu}+\mathbf{q}_{\nu}$ is $0$ modulo a
reciprocal lattice vector $\mathbf{G}$.

With the preceding approximation of $\mathbf{K}$, we now have all the
ingredients to compute the power spectrum. Fourier transforming eq.(49) in
frequency domain, and using eq.(38) and eq.(39), we obtain the following set
of self-consistent equations

$\displaystyle\Gamma_{\lambda}(\omega)=$
$\displaystyle\frac{\hbar}{2\mathrm{k}_{\mathrm{B}}\mathrm{T}\Omega_{\lambda}}\sum_{\mu,\nu}|\Psi(-\lambda,\mu,\nu)|^{2}(G_{\mu}\ast
G_{\nu})(\omega)$ (52) $\displaystyle\Delta_{\lambda}(\omega)=$
$\displaystyle\frac{1}{\pi}\int
d\omega^{\prime}\frac{\Gamma_{\lambda}(\omega^{\prime})}{\omega^{\prime}-\omega}$
(53) $\displaystyle G_{\lambda}(\omega)=$
$\displaystyle\frac{\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\pi\hbar}\frac{8\Omega_{\lambda}\Gamma_{\lambda}(\omega)}{(\omega^{2}-\Omega_{\lambda}^{2}-2\omega\Delta_{\lambda}(\omega))^{2}+4\omega^{2}\Gamma_{\lambda}^{2}(\omega)}$
(54)

In order to obtain a starting point to the self-consistent cycle, we introduce
an approximation inspired by scattering theory. We begin by defining a bare
Green’s function as the solution of a memory-less GLE

$\ddot{G}_{\lambda,0}(t)=-\Omega^{2}_{\lambda}G_{\lambda,0}(t)$ (55)

which is simply given by
$G_{\lambda,0}(t)=\frac{\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\hbar\Omega_{\lambda}}\cos(\Omega_{\lambda}t)$.
The scattering theory approximation is obtained by replacing the Green’s
function in eq.(52) by $G_{\lambda,0}(\omega)$, which allows to obtain the
linewidth

$\Gamma_{\lambda}^{s}(\omega)=\frac{\pi\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\hbar\Omega_{\lambda}\Omega_{\mu}\Omega_{\nu}}\sum_{\mu,\nu,s,s^{\prime}}|\Psi(-\lambda,\mu,\nu)|^{2}\delta(\omega+s\Omega_{\mu}+s^{\prime}\Omega_{\nu})$
(56)

with $s,s^{\prime}=\pm 1$. From this result, a first approximation to the
power spectrum can be obtained by using eq.(53) and (54) with
$\Gamma_{\lambda}^{s}(\omega)$ as an input.

##  VI Computing the frequencies and vertex from (PI)MD simulations

One advantage of the formalism derived previously is that only static averages
are needed as input. This is particularly interesting for applications given
that the static KCF of two variables $A$ and $B$ can be computed exactly in
the path-integral molecular dynamics framework as [63, 49]:

$\begin{split}G_{AB}(0)=&\lim_{P\rightarrow\infty}\lim_{t\rightarrow\infty}\int_{0}^{t}\frac{e^{-\beta
H(t)}}{\mathcal{Z}}\big{(}\sum_{p}^{P}A_{p}(t)\big{)}\big{(}\sum_{p}^{P}B_{p}(t)\big{)}dt\\\
\approx&\frac{1}{N}\sum_{n}^{N}\bar{A}(t_{n})\bar{B}(t_{n})\end{split}$ (57)

where $P$ is the number of beads in the classical polymer,
$\bar{A}(t)=\sum_{p}A_{p}(t)$ is the mean of the property $A$ over the
classical polymer. The second line of eq.(57) corresponds to the approximation
with a finite number of configurations and beads from path integral molecular
dynamics.

The static correlation function of classical MD simulations corresponds to the
limit $P=1$ and can be computed as

$\begin{split}G_{AB}(0)\approx&\frac{1}{N}\sum_{n}^{N}A(t_{n})B(t_{n})\end{split}$
(58)

We note that this can be used as a way to approximate the static KCF, as will
be discussed on the next section.

In the limit of a finite number $N$ of configurations, the frequency matrix
can be computed with

$\Theta_{ij}=-\frac{\sum_{n}^{N}\sum_{k}\bar{u}_{k}(t_{n})\bar{f}_{i}(t_{n})}{\sum_{n}^{N}\sum_{k}\bar{u}_{k}(t_{n})\bar{u}_{j}(t_{n})}$
(59)

where we used the property
$\big{(}u_{i},f_{j}\big{)}=-\mathrm{k}_{\mathrm{B}}\mathrm{T}\delta_{ij}$ to
decouple the sums in the numerator and denominator of eq.(13). In this
equation, one can recognize the least squares solution $\mathbf{\Theta}$ of
the linear problem $\mathbf{\Theta}\mathbf{u}=-\mathbf{f}$. This implies that
the computation of the frequency matrix can be mapped to the problem of
fitting an effective harmonic Hamiltonian using the distribution of forces and
displacements at the studied temperature. This corresponds to the TDEP method
at second order [37], for which our framework offers an exact formal
justification and a path for generalization. In practice, the average needed
in the computation of $\mathbf{\Theta}$ might necessitate long simulations
which could become prohibitive. Fortunately, $\mathbf{\Theta}$ is a sparse
matrix, with many coefficients related by symmetry, as it must be invariant
under the symmetry group of the crystal. Moreover, since the displacements and
forces are invariant with respect to global translation and rotation of the
system, their static KCF should also respect these invariances. This means
that, by definition, $\mathbf{\Theta}$ follows the usual symmetries found in
the IFC tensor of the harmonic approximation for phonons. Consequently,
imposing them on $\mathbf{\Theta}$ before applying eq.(59) will considerably
decrease the length of the simulations needed, similarly to what is done in
lattice dynamic fitting methods [38, 72, 73, 74]. A similar analysis can be
done for the second order vertex in real space, which can be computed as

$\Theta_{ijk}=-\frac{1}{2}\sum_{n}^{N}\frac{\sum_{l,m}\bar{u}_{l}(t_{n})\bar{u}_{m}(t_{n})\delta\bar{f}_{i}(t_{n})}{\sum_{l,m}\bar{u}_{l}(t_{n})\bar{u}_{m}(t_{n})\bar{u}_{j}(t_{n})\bar{u}_{k}(t_{n})}$
(60)

with $\delta f_{i}=f_{i}+\sum_{j}\Theta_{ij}u_{j}$. A simple variable
transformation allows to rewrite this equation as a least-squares solution of
the equation $\mathbf{\Theta}^{T}\mathbf{U}=\delta\mathbf{f}$, where
$\mathbf{U}=\mathbf{u}\mathbf{u}^{T}$ [39, 74]. As in the frequency matrix
case, this means that this interaction can be constructed exactly as the
third-order IFC in the TDEP method. It can be shown that this TDEP-like
construction arises naturally for every order, by definition of the $n$-th
order projectors $\mathcal{P}_{n}$.

##  VII Discussion

###  VII.1 Classical-like correlation functions and classical limit

Looking at the main results in the preceding derivation, it would be easy to
miss the quantum character of the dynamics described, since most of the
equations seem to lack the usual quantum-related prefactors. The quantum
behavior of the system is hidden in the Kubo correlation function, as can be
inferred from its Lehmann representation given by

$G_{\lambda}(\omega)=\mathrm{k}_{\mathrm{B}}\mathrm{T}\sum_{jk}p_{j}\frac{|\bra{j}A_{\lambda}\ket{k}|^{2}}{\hbar\omega(n(\omega)+1)}\delta\bigg{(}\omega-\frac{E_{j}-E_{k}}{\hbar}\bigg{)}$
(61)

where $p_{j}=e^{-\beta E_{j}}/\mathcal{Z}$ is the Boltzmann weight associated
with state $j$ of energy $E_{j}$ and $n(\omega)=1/(e^{\beta\hbar\omega}-1)$ is
the Bose-Einstein distribution. Consequently, our approach allows to alleviate
any ambiguity related to the type of quantum correlation functions one should
use in PIMD simulations, all while giving a consistent formulation of the
dynamics of the system.

Interestingly, the classical limit of this mode-coupling approach is
completely analogous to its quantum counterpart [60]. Indeed, the classical
version of the Mori-Zwanzig formalism is obtained by simply replacing the KCF
with classical correlation functions. Using the similarity between the
properties of Kubo and classical correlation functions, the same derivation
can be used to obtain the classical version of the set of self-consistent
equations (52), (53) and (54). Therefore, the only difference between the
quantum and classical description comes from the type of correlation functions
used to compute the frequency matrix $\mathbf{\Theta}$ and the vertex. This
classical-like behavior is particularly visible in the KCF formulation of the
quantum fluctuation-dissipation theorem given by eq.(36), which is exactly the
same as its classical limit [62]. This observation brings a new justification
to the interpretation of classical correlation functions as approximated KCF
[47]. From a practical point of view, and as discussed in the previous
section, this means that the same set of equations can be used in the two
limits, the only difference stemming from the type of simulation used to
compute the correlation functions.

Further, our approach allows for a simple inclusion of quantum behavior when
doing classical simulations. For this, we introduce the lesser Green’s
function $G^{<}_{\lambda}(t)=\braket{A^{\dagger}_{\lambda}(t)A_{\lambda}}$,
given in the Lehmann representation by

$G^{<}_{\lambda}(\omega)=\sum_{jk}p_{j}e^{-\beta\hbar\omega}|\bra{j}A_{\lambda}\ket{k}|^{2}\delta\bigg{(}\omega-\frac{E_{j}-E_{k}}{\hbar}\bigg{)}$
(62)

From this representation, it is easy to see that the Kubo and lesser Green’s
functions are related by
$G_{\lambda}^{<}(\omega)=\frac{\hbar\omega}{\mathrm{k}_{\mathrm{B}}\mathrm{T}}n(\omega)G_{\lambda}(\omega)$.
Replacing the Kubo Green’s functions in eq.(52) with the lesser Green’s
functions makes the Bose-Einstein distributions visible in the equation

$\begin{split}\Gamma_{\lambda}(\omega)=&\frac{\hbar}{2\mathrm{k}_{\mathrm{B}}\mathrm{T}\Omega_{\lambda}}\sum_{\mu,\nu}|\Psi(-\lambda,\mu,\nu)|^{2}\\\
&\times\int_{-\infty}^{\infty}d\omega^{\prime}\frac{\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\hbar(\omega-\omega^{\prime})n(\omega-\omega^{\prime})}G_{\lambda}^{<}(\omega-\omega^{\prime})\\\
&\times\frac{\mathrm{k}_{\mathrm{B}}\mathrm{T}}{\hbar\omega^{\prime}n(\omega^{\prime})}G_{\lambda}^{<}(\omega^{\prime})\end{split}$
(63)

In the scattering theory approximation,
$G_{\lambda}^{<}(t)=n(\Omega_{\lambda})\cos(\Omega_{\lambda}t)$, so that
$\Gamma_{\lambda}^{s}(\omega)$ reduces to

$\displaystyle\Gamma_{\lambda}^{s}(\omega)=$
$\displaystyle\frac{\hbar\pi}{4\mathrm{k}_{\mathrm{B}}\mathrm{T}\Omega_{\lambda}}\sum_{\mu,\nu}|\Psi(-\lambda,\mu,\nu)|^{2}F(\omega,\mu,\nu)$
(64) $\displaystyle\begin{split}F(\omega,\mu,\nu)=&\sum_{s=\pm
1}\bigg{[}\big{(}n_{\mu}+n_{\nu}+1\big{)}\delta(\omega+s\Omega_{\mu}+s\Omega_{\nu})\\\
&+\big{(}n_{\mu}-n_{\nu}\big{)}\delta(\omega+s\Omega_{\mu}-s\Omega_{\nu})\bigg{]}\end{split}$
(65)

with $n_{\mu}=n(\Omega_{\mu})$, and where we used $\lim_{\omega\rightarrow
0}1/[\hbar\omega n(\omega)]=\beta$. Once $\Psi$ has been calculated, using
eq.(64) to process inputs from MD simulations is a simple zero-cost and
legitimate way to approximate quantum effects without needing to do expensive
PIMD simulations. Such an approximation should be valid for systems with low
anharmonicity or at high temperature.

###  VII.2 Interpretation of the static quasiparticles

With the previous discussions, we focused on dynamical properties. It is
interesting to also consider static properties associated with the formalism,
in particular the quasiparticle frequencies. Interestingly, these temperature-
dependent frequencies have often been compared with great success to
experiments [35, 37, 43], but these comparisons lacked a rigorous theoretical
justification. It has been proposed to interpret these frequencies as an
approximation of the first excitation of the system [75, 49]. While such an
approximation would be accurate at low temperature, it doesn’t hold for high
temperature and completely misses the exact character associated with the
quantity [75].

Taking the $\omega\rightarrow 0$ limit of eq.(35), the static susceptibility
can be written as $\raisebox{2.0pt}{\mbox{\large$\chi$}}=-\beta G(0)$. From
eq.(13) and using the property
$\big{(}u_{i},f_{j}\big{)}=-\mathrm{k}_{\mathrm{B}}\mathrm{T}\delta_{ij}$, we
obtain a proportionality relation between the frequency matrix and the inverse
static susceptibility in real space

$\frac{\mathbf{\Theta}}{\sqrt{\mathbf{M}^{T}\mathbf{M}}}=\mathbf{\raisebox{2.0pt}{\mbox{\large$\chi$}}}^{-1}$
(66)

meaning that the frequency matrix is directly related to the instantaneous
response of the system to an external perturbation. By doing the
transformation in section III, and due to the Hermitian character of the
matrices involved, we can conclude that the quasiparticle polarizations
$\pmb{\epsilon}_{\lambda}(\mathbf{q})$ are also eigenvectors of the Fourier
transform of the susceptibility matrix, with eigenvalues
$\hbar\Omega_{\lambda}^{-2}(\mathbf{q})$. From this observation, we can
interpret the eigenvectors $\pmb{\epsilon}^{\lambda}(\mathbf{q})$ of the
quasiparticles as displacement patterns associated to a perturbation with a
momentum $\mathbf{q}$. The $A_{\lambda}(\mathbf{q})$ are thus the amplitude of
these displacement patterns, which can be seen as a temperature-dependent
generalization of the normal modes of the harmonic theory.

To have a better understanding of the meaning of these frequencies, it is
useful to focus on the second order time derivative of the KCF. From the
relation $\ddot{G}_{ij}=\mathrm{k}_{\mathrm{B}}\mathrm{T}$, and using eqs.
(22) and (24), we have
$\ddot{G}_{\lambda}(\mathbf{q})/G_{\lambda}(\mathbf{q})=\Omega_{\lambda}^{2}(\mathbf{q})$.
Using then the Lehmann representation of the KCF $G_{\lambda}(\mathbf{q})$ and
$\ddot{G}_{\lambda}(\mathbf{q})$

$\displaystyle G_{\lambda}(\mathbf{q})$
$\displaystyle=\mathrm{k}_{\mathrm{B}}\mathrm{T}\sum_{jk}\frac{|\bra{j}A_{\lambda}(\mathbf{q})\ket{k}|^{2}(p_{k}-p_{j})}{\omega_{jk}}$
(67)
$\displaystyle\begin{split}\ddot{G}_{\lambda}(\mathbf{q})&=\mathrm{k}_{\mathrm{B}}\mathrm{T}\sum_{jk}\frac{\bra{j}A_{\lambda}(\mathbf{q})\ket{k}\bra{k}\ddot{A}_{\lambda}(\mathbf{q})\ket{j}(p_{k}-p_{j})}{\omega_{jk}}\\\
&=\mathrm{k}_{\mathrm{B}}\mathrm{T}\sum_{jk}\frac{\omega_{jk}^{2}|\bra{j}A_{\lambda}(\mathbf{q})\ket{k}|^{2}(p_{k}-p_{j})}{\omega_{jk}}\\\
\end{split}$ (68)

with $\omega_{jk}=(E_{j}-E_{k})/\hbar$, and where we used
$\bra{j}\ddot{A}_{\lambda}(\mathbf{q})\ket{k}=-\omega_{jk}^{2}\bra{j}A_{\lambda}(\mathbf{q})\ket{k}$,
we obtain a formulation of the frequencies in terms of a complete set of
eigenstates of the Hamiltonian

$\displaystyle c_{jk}(\lambda,\mathbf{q})$
$\displaystyle=\frac{|\bra{j}A_{\lambda}(\mathbf{q})\ket{k}|^{2}\big{(}p_{k}-p_{j})}{\omega_{jk}}$
(69) $\displaystyle\Omega_{\lambda}(\mathbf{q})$
$\displaystyle=\sqrt{\frac{\sum_{jk}\omega_{jk}^{2}c_{jk}(\lambda,\mathbf{q})}{\sum_{jk}c_{jk}(\lambda,\mathbf{q})}}$
(70)

This result shows that the quasiparticle frequencies correspond to a
thermodynamic average of the transitions between states of the system
projected on the displacement patterns.

It should be noted that these frequencies do not necessarily align with peaks
in the spectral functions. Indeed, the relation between the static frequencies
and the spectral function is given by the sum rule

$\int_{0}^{\infty}d\omega\raisebox{2.0pt}{\mbox{\large$\chi$}}^{\prime\prime}_{\lambda}(\mathbf{q},\omega)\omega=\Omega_{\lambda}(\mathbf{q})$
(71)

Consequently, the peaks and the static frequencies should be aligned only in
the case of a symmetric lineshape, or at least in cases where the spectral
weight is evenly distributed around $\Omega_{\lambda}$.

###  VII.3 Comparison with other methods

To help understanding the implication of the mode-coupling theory of lattice
dynamics, it is informative to compare it to well established methods. In this
section, we show that the perturbation theory, the SCHA and the TDEP method
can be understood as different approximations of the mode-coupling approach,
that we summarize in table 1.

Table 1: Summary of different approach of anharmonic lattice dynamics.  | Perturbation theory | SCHA | Mori-Zwanzig
---|---|---|---
Sampling of displacements |  Perturbative through finite difference or DFPT |  Effective harmonic (multivariate Gaussian) |  Exact through MD (classical case) or PIMD (quantum case)
Static quasiparticles |  Harmonic phonons |  Effective harmonic phonons or Eigenstates of the inverse susceptibility |  Eigenstates of the inverse susceptibility
Vertex in the bubble contribution | perturbative | Gaussian approximation | exact
Implementations |  Phono3py [15] ShengBTE [16] |  SSCHA [30] Alamode [27, 28] sTDEP [76] qSCAILD [31] |  MD-TDEP [37, 38]

####  VII.3.1 The harmonic limit

A good theory of lattice dynamics should reduce to the exactly known harmonic
limit. This can be trivially shown to be true for the mode-coupling theory.
Indeed, for a harmonic system, the frequency matrix presented in eq.(13) is
equal to the harmonic force constants, so that the random forces cancel out.
Consequently, the GLE reduces to the expected memory-less harmonic equations
of motion.

####  VII.3.2 Perturbation theory

More interestingly, it can be shown that the perturbative limit of the mode-
coupling approach allows us to retrieve the usual linewidth formulation
obtained from a diagrammatic expansion [11, 13]. To show this, we suppose that
the potential energy of the system can be approximated by an anharmonic
potential truncated to third order

$\displaystyle V(\mathbf{R})\approx$
$\displaystyle\frac{1}{2}\sum_{ij}\Theta_{ij}^{\mathrm{pt}}u_{i}u_{j}+\frac{1}{3!}\sum_{ijk}\Psi_{ijk}^{\mathrm{pt}}u_{i}u_{j}u_{k}$
(72) $\displaystyle\Theta_{ij}^{\mathrm{pt}}=$
$\displaystyle\frac{\partial^{2}V(\mathbf{R})}{\partial R_{i}\partial R_{j}}$
(73) $\displaystyle\Psi_{ijk}^{\mathrm{pt}}=$
$\displaystyle\frac{\partial^{3}V(\mathbf{R})}{\partial R_{i}\partial
R_{j}\partial R_{k}}$ (74)

where we assume that contributions from $\Psi_{ijk}^{\mathrm{pt}}$ are
perturbations with respect to $\Theta_{ij}^{\mathrm{pt}}$, so that averages
can be taken with respect to the harmonic contribution. In this approximation,
the frequency matrix is equal to the second order force constants
$\mathbf{\Theta}^{\mathrm{pt}}$ and the memory kernel can be written

$K_{ij}^{\mathrm{pt}}(t)=\sum_{klmn}\Psi_{ikl}^{\mathrm{pt}}\Psi_{jmn}^{\mathrm{pt}}\big{(}u_{k}u_{l},e^{-i\mathcal{Q}\mathcal{L}t}u_{m}u_{n}\big{)}$
(75)

which corresponds to the $(2,2)$-order memory kernel defined in real space,
but with third-order force constants instead of the full vertex. From this
observation, the derivation of the phonon quasiparticles follows the same
route as the more general Mori-Zwanzig theory. In the end, the main difference
between the perturbation and mode-coupling theories lies in the definition of
the quasiparticles and in the strength of the vertex $\Psi(\lambda,\mu,\nu)$.
In order to quantitatively reconcile experimental observations and theoretical
predictions for properties such as thermal conductivity, these
renormalizations are important. Indeed, while the temperature dependence of QP
frequencies is known to have a large impact on this kind of quantity [77], the
deviations brought by scattering renormalization can be of similar magnitude
and have the potential to remove the expected high temperature trend
$\kappa\propto T^{-\alpha}$, as shown in Refs. [78, 79, 80].

From a diagrammatic derivation to lowest-order, this dynamical contribution is
given by the bubble diagram represented in Fig.2 a). Doing a naive comparison,
it would seem that the other lowest-order diagrams, known as the tadpole
(Fig.2 b)) and the loop (Fig.2 c)), are missing from the mode-coupling theory.
However, the absence of these diagrams is correct, and a consequence of some
exact properties of this approach. By construction, the Mori-Zwanzig
projection formalism used in section II.3 results in a conservative
contribution that gives the _exact_ static part of the Kubo Green’s function.
The tadpole and loop diagrams are purely real, so that their contributions to
the Green’s function only result in a shift of the static terms [11, 28]. From
this observation, one can infer that the frequencies
$\Omega_{\lambda}(\mathbf{q})$ already include all static diagrams to all
orders. Adding any static diagram, such as the loop or the tadpole, would
result in an over-correction, thus explaining their absence from the lowest-
order mode-coupling approximation.

Figure 2: Feynman diagrams appearing in the self-energy in perturbation theory
to lowest order. Each line represent a propagator while dots and squares are
respectively third and fourth order vertices. a) The bubble b) the tadpole and
c) the loop. In the mode-coupling approach presented here, we observe only
bubble-like contributions in the memory kernel.

Even if the approaches are quite similar, we should expect the mode-coupling
theory to have a wider range of validity than bare perturbation theory.
Indeed, being constructed on a truncated expansion of the Hamiltonian with
respect to displacements, perturbation theory is known to fail for systems
with strong anharmonicity. In these systems, the high-order contributions to
the Hamiltonian can be of similar - or even higher - magnitude than the
harmonic contribution. This is because the perturbative nature of the high
order contributions is only a supposition, and there is no guarantee that this
holds in general [81, 17, 21]. On the contrary, the projection operators of
eq.(8) and (41) are constructed so that effects of higher order correlations
are minimized. Such a construction allows the mode-coupling theory to be valid
even in strongly anharmonic systems. One should note that the use of linear
response theory for this kind of system is not a problem for the Mori-Zwanzig
approach, as the linearity is assumed in the response of the system to a weak
probe, not in the behavior of the system itself.

####  VII.3.3 Self-consistent harmonic approximation

The self consistent harmonic approximation is a workhorse for the study of
anharmonic crystals. The method is based on the minimization of the Gibbs-
Bogoliubov free energy $\widetilde{\mathcal{F}}$ defined as

$\widetilde{\mathcal{F}}=\widetilde{\mathcal{F}}_{0}+\braket{V(\mathbf{R})-\widetilde{V}(\mathbf{R})}_{0}$
(76)

where $\widetilde{V}(\mathbf{R})$ is an effective potential,
$\widetilde{\mathcal{F}}_{0}$ is its associated free energy and
$\braket{O}_{0}$ is an average computed within the distribution associated
with $\widetilde{V}(\mathbf{R})$. As the name of the method implies, the
effective potential has a harmonic form

$\widetilde{V}(\mathbf{R})=(\mathbf{R}-\pmb{\mathcal{R}})^{\mathrm{T}}\mathbf{\Theta}^{\mathrm{SCHA}}(\mathbf{R}-\pmb{\mathcal{R}})$
(77)

where the effective force constants $\mathbf{\Theta}^{\mathrm{SCHA}}$ and
equilibrium positions $\pmb{\mathcal{R}}$ are parameters to optimize. The
advantage of such a potential is that the distribution of displacements around
$\mathbf{\mathcal{R}}$ is a multivariate Gaussian. Consequently, the sampling
of atomic positions to compute averages is greatly simplified compared to the
true potential. Using Jensen’s inequality, it can be shown that the Gibbs-
Bogoliubov is an upper bound to the real free energy
$\mathcal{F}<\widetilde{\mathcal{F}}$ [82], so that optimizing the parameters
means minimizing $\widetilde{\mathcal{F}}$. For the effective force constants,
the minimum is obtained when

$\Theta_{ij}^{\mathrm{SCHA}}=\Braket{\frac{\partial^{2}V(\mathbf{R})}{\partial
u_{i}\partial u_{j}}}_{0}$ (78)

while the effective equilibrium positions are optimized when the average
difference of true and harmonic forces is zero. In practice, the minimization
is done self-consistently, and several approaches and implementations can be
found in the literature [26, 29, 30, 27, 28, 31, 76].

After the minimization, while the effective force constants can be used to
approximate quasiparticles, a more consistent result is obtained by using the
Hessian $\mathbf{\Xi}$ of the free energy [29, 30]

$\Xi_{ij}=\frac{\partial^{2}\widetilde{\mathcal{F}}}{\partial\mathcal{R}_{i}\partial\mathcal{R}_{j}}$
(79)

By seeing the Gibbs-Bogoliubov free energy as a Landau expansion with an order
parameter $\mathcal{R}$, the free energy Hessian $\mathbf{\Xi}$ can be
understood as the inverse static susceptibility. However, compared to the
mode-coupling approach, where this quantity is exact,
$\mathbf{\raisebox{2.0pt}{\mbox{\large$\chi$}}}^{-1}$ is approximated in the
SCHA, as only the contribution from the first-order free energy cumulant
enters in $\widetilde{\mathcal{F}}$.

On top of the quasiparticles computed from the free energy Hessian, the
dynamical susceptibility can be obtained with a bubble-like contribution [29,
30] in which the vertices are computed using the average of the third-order
derivative of the potential energy

$\begin{split}\Psi_{ijk}^{\mathrm{SCHA}}=&\frac{1}{\sqrt{M_{i}M_{j}M_{k}}}\Braket{\frac{\partial^{3}V(\mathbf{R})}{\partial
R_{i}\partial R_{j}\partial R_{k}}}_{0}\\\
=&-\frac{1}{\sqrt{M_{i}M_{j}M_{k}}}\sum_{lm}\frac{\braket{u_{l}u_{m}f_{k}}_{0}}{\braket{u_{i}u_{l}}_{0}\braket{u_{j}u_{m}}_{0}}\end{split}$
(80)

where the second line can be derived with an integration by parts. This
procedure can also be compared with our mode-coupling theory. To show this, we
start by assuming that we can approximate the real distribution with the
Gaussian associated with $\mathbf{\Theta}^{\mathrm{SCHA}}$. With this
assumption, the frequency matrix is actually given by the SCHA force constants
$\mathbf{\Theta}^{\mathrm{SCHA}}$, and the random forces are given with by
$\delta f_{i}=f_{i}+\sum_{j}\Theta_{ij}^{\mathrm{SCHA}}u_{j}$. Applying the
mode-coupling decomposition of section V to approximate
$\mathbf{K}^{\mathrm{SCHA}}(t)$, it is then simple to show that the vertex
appearing in the $(2,2)$-order memory kernel in real space is given by
eq.(80), as the average of an odd number of displacements cancels out with the
Gaussian average. From this observation, we can conclude that the bubble-like
dynamical contribution in the SCHA is an approximation of the mode-coupling
theory.

To summarize, two approximations enter in the SCHA compared to the mode-
coupling theory. The first one concerns the static frequencies, which are
approximated by the effective harmonic frequencies corrected by a first-order
cumulant of the free energy. The second one is in the vertices in the bubble
contribution, which are computed on the basis of a Gaussian approximation.
From a computational point-of-view, the SCHA bears the important advantage of
a lower cost due to its simplified sampling compared to (PI)MD. This is
particularly true for the inclusion of nuclear quantum effects which can be
added in a simple way. However higher-order corrections and applications to
strongly anharmonic (and hence non-Gaussian) systems are more delicate to
justify.

####  VII.3.4 Average of the Hessian

In the literature, it has been proposed to define a generalization of the
harmonic force constants through an average of the potential energy Hessian
[49]

$\begin{split}\Pi_{ij}=&\bigg{(}\frac{\partial^{2}V(\mathbf{R})}{\partial
R_{i}\partial R_{j}}\bigg{)}\end{split}$ (81)

where the only difference with the SCHA force constants of eq.(78) lies in the
Hamiltonian used to compute the average. From the harmonic limit, where
$\Pi_{ij}$ reduces to the harmonic force constants, and the similarity with
eq.(78), this generalization seems like a sensible choice. However, it is
constructed on the equality $\Pi_{ij}=\Theta_{ij}$ that only holds in a
harmonic system.

To highlight the differences between this formulation and the mode-coupling
theory, the potential energy Hessian average can be expanded in term of force-
force KCF [49]

$\Pi_{ij}=\frac{\big{(}f_{i},f_{j}\big{)}}{\mathrm{k}_{\mathrm{B}}\mathrm{T}}$
(82)

Using this equation, the static part of memory kernel can be formulated as

$K_{ij}(0)=\Pi_{ij}-\Theta_{ij}$ (83)

Since $\mathbf{K}$ is a positive definite matrix [69], this equation tells us
that $\mathbf{\Pi}$ is greater that $\mathbf{\Theta}$, meaning that
frequencies computed from the diagonalization of the average Hessian will
overestimate the $\Omega_{\lambda}$. This result shows that while formulations
based on the SCHA and on the true distribution by means of (PI)MD are similar,
they are not interchangeable.

####  VII.3.5 Temperature dependent effective potential

TDEP is a method based on the construction of an effective anharmonic
potential of the form (here presented up to third order)

$\widetilde{V}(\mathbf{R})=\frac{1}{2}\sum_{ij}\Theta_{ij}^{\mathrm{TDEP}}u_{i}u_{j}+\frac{1}{3!}\sum_{ijk}\Psi_{ijk}^{\mathrm{TDEP}}u_{i}u_{j}u_{k}$
(84)

The temperature dependence in this potential is introduced through a
successive fitting of each order. First, using a set of positions and forces
distributed according to a given temperature, the second order interatomic
force constants $\Theta_{ij}$ are fit using ordinary least-squares. Then the
third order force constants $\Psi_{ijk}$ are fit on the residual of the
forces. In principle, this successive fitting can be continued to any order,
though usual applications of the method truncate the potential at the third or
fourth order.

In applications of TDEP, two approaches have been used to generate the set of
positions and forces. In the original works [37, 38], that we will call MD-
TDEP, MD simulations are used. In the second one, sometimes called sTDEP (for
stochastic TDEP), the displacements are computed from an effective harmonic
model. In the latter case, the set of forces and displacements have to be
constructed self-consistently, by refitting the second order potential from
the set of the previous iteration.

It is important to differentiate the approaches. On one hand, sTDEP is simply
an implementation of the SCHA, since the self-consistent least-squares fit
allows to minimize the Gibbs-Bogoliubov free energy in eq.(77). On the other
hand, by using the real distribution, MD-TDEP is an implementation of the
mode-coupling theory. Indeed, as has already been discussed in section VI, the
force constants $\mathbf{\Theta}^{\mathrm{TDEP}}$ and
$\mathbf{\Psi}^{\mathrm{TDEP}}$ computed from the TDEP methods are equal to
the frequency matrix $\mathbf{\Theta}$ and the real space vertex
$\mathbf{\Psi}$. Consequently, the perturbation theory usually applied to
compute the spectral function is equivalent to the scattering approximation of
mode-coupling theory eq.(64) if only the bubble term is used. As this is the
case for most applications of TDEP in the literature, we have a large number
of examples showing the accuracy of the mode-coupling approach. For instance,
the method has been used successfully to explain the anomalous neutron
scattering of SnTe and PbTe [83], to study the phase transition of GeTe [84]
or the phase diagram of uranium [85, 86] as well as many other phenomena for
highly anharmonic systems [87, 40, 88, 89]. Our work provides a theoretical
justification for the results of these studies, as well as a generalization,
allowing to include quantum effects through the use of PIMD and to include
higher order corrections.

###  VII.4 Markovian approximation and Lorentzian lineshape

In the limit of small anharmonicity, the spectral lineshapes are often
approximated with Lorentzians [15, 90, 7]. An interesting conclusion to be
taken from the presented theory is that such a solution cannot be an exact
representation of the dynamics of the nuclei. In the Lorentzian approximation,
the frequency dependency of the lifetime is neglected and $\Gamma(\omega)$ is
approximated with a constant given by
$\Gamma_{\lambda}=\Gamma_{\lambda}(\Omega_{\lambda})$ [15]. In the Mori-
Zwanzig formalism presented here, such Lorentzian spectra would appear in the
context of a Markovian approximation of the memory kernel [54]. In this case,
where it is assumed that $\mathbf{K}(t)$ decays very quickly compared to
$\mathbf{G}(t)$, the memory kernel can be replaced by a constant friction term
$\mathbf{\Gamma}^{M}(\mathbf{q})=\beta\int_{0}^{\infty}dt\mathbf{K}(\mathbf{q},t)$,
so that the Kubo transformed Green’s function is written as

$\ddot{\mathbf{G}}(\mathbf{q},t)=-\mathbf{\Omega}^{2}(\mathbf{q})\mathbf{G}(\mathbf{q},t)-\mathbf{\Gamma}^{M}(\mathbf{q})\dot{\mathbf{G}}(\mathbf{q},t)$
(85)

with solutions

$\mathbf{G}(\mathbf{q},t)=\cos\left(\mathbf{\Omega}(\mathbf{q})t\right)e^{-\mathbf{\Gamma}^{M}(\mathbf{q})t}$
(86)

However, it can be observed that this solution is not an even function of
time, as requested for the KCF. For instance, in a Markovian approximation,
the property $\frac{d\mathbf{G}(t)}{dt}|_{t=0}=0$ would not be respected. This
means that a Markovian approximation, and consequently a Lorentzian
approximation, will never fully represent the dynamics of the system.

Nevertheless, for midly anharmonic systems, the Markovian limit of the GLE
should still be a reliable approximation, and Lorentzian-like lineshapes are
frequently observed in spectroscopy measurements. Consequently, a derivation
of the inputs needed for such an approximation is of interest, in order to
provide an efficient way to compute the dynamical properties of a system of
nuclei from (PI)MD.

Figure 3: Spectral function
$\raisebox{2.0pt}{\mbox{\large$\chi$}}^{\prime\prime}(\omega)$ of fcc 4He at
$9.03$ cm3/mol and $38$ K computed with the mode-coupling approach in the
scattering-like approximation. The frequency matrix and second order vertex
were computed using PIMD for a), classical MD for b), the SCHA for c) and
through perturbation theory for d). The green dots are from inelastic neutrons
scattering experiment[91].

##  VIII Applications

To demonstrate the accuracy of the theory developed here, we apply it to the
fcc phase of 4He. Due to its very small atomic mass, helium is known to show
strong nuclear quantum effects (NQE) making it an archetypal example of a
quantum crystal [92]. Indeed, the zero-point motion in solid helium is large
enough to make the system explore the BOS far beyond the range of
applicability of the harmonic approximation. This unusual behavior manifests
itself in the largest Lindemann ratio (i.e. $\sqrt{\braket{u^{2}}}/a$ with
$\braket{u^{2}}$ the mean squared atomic displacement and $a$ the lattice
parameter) of any material observed at cryogenic temperature [92]. The high
anharmonicity in helium, which can be found in other rare gas solids, is what
started the field of anharmonic lattice dynamics and spawned the development
of theories such as the SCHA [4, 22, 23, 93].

Given these properties, this system represents a perfect case to illustrate
that both NQE and strong anharmonicity can be captured with our mode-coupling
approach. To emphasize this point, we performed numerical simulations in a
fully quantum setting using PIMD as well as classical MD simulations. For the
sake of comparison, we also performed simulations using both the SCHA and
perturbation theory. To allow for a comparison with inelastic neutron
scattering results [91], all simulations were run at a temperature of $38$ K
and with a molar volume of $9.03$ mol3/mol, which correspond to a lattice
parameter of $3.915\leavevmode\nobreak\ \mathring{a}$.

To describe the BOS of helium, we use the Aziz pair potential [94], as it
offers an accurate representation of the solid phases of helium at low
pressure [92]. PIMD and MD simulations were done using the implementation of
the Langevin thermostat in LAMMPS [95], for a duration of $125$ ps with a
timestep of $0.5$ fs in both cases and 64 beads for the PIMD. To compute the
frequency matrix and the second order vertex in real space, we used the
recursive least-square fit presented in section VI. After projecting these
tensors on QP following subsection III, we solved the mode-coupling theory
equations in the scattering approximation using the greater and lesser Green’s
function formalism of subsection VII.1. For this step, a $12\times 12\times
12$ $\mathbf{q}$-point grid was sufficient to obtain converged results. We
note that for both PIMD and MD results, the memory kernel is obtained within a
quantum formalism, though NQE are included in the frequency matrix and second
order vertex only in the PIMD case. For the SCHA, we used the stochastic TDEP
approach presented in subsection VII.3.5 with Bose-Einstein occupations in
order to capture NQE.

Before presenting the spectral properties of 4He, we first assess the
anharmonicity of the system. To this end, we adopt a recently proposed measure
of anharmonicity [5] to our Mori-Zwanzig formalism as

$\sigma^{\mathrm{A}}=\sqrt{\frac{\sum_{i}\braket{\delta
f_{i}^{2}}}{\sum_{i}\braket{f_{i}^{2}}}}$ (87)

We note, however, that used in this way $\sigma^{\mathrm{A}}$ is not strictly
a measure of anharmonicity, since the frequency matrix already goes beyond the
harmonic contribution of the BOS. It is more accurate to see
$\sigma^{\mathrm{A}}$ as a measure of the dissipative contribution to the
dynamics of the system. For fcc 4He, we find $\sigma^{\mathrm{A}}=0.62$ when
computed with PIMD, which would indicate a very dissipative system. Similarly,
using a formulation based on the frequency matrix [74], we find a Debye
temperature of $174$ K, close to the indirect experimental estimation of $154$
K [91]. At $38$ K, this means that zero-point motion is expected to be
important. Coupled with the high $\sigma^{\mathrm{A}}$, this confirms that the
He system is a stringent test for the ability of the mode-coupling theory to
capture both NQE and anharmonicity.

For each of the simulations, the theoretical spectra are presented in Fig.3
where they are compared with inelastic neutron scattering results [91]. Among
all the results, PIMD mode coupling presents the best agreement with
experiment, with peaks of the theoretical spectral function fully aligning
with the experimental phonon frequencies. On the contrary, the classical MD
case severely underestimates phonon frequencies for all reciprocal space
directions. This underestimation is even more important when using (harmonic,
0K based) perturbation theory, showing the need to go to non-perturbative
methods for an accurate description of this system. The spectra computed with
the SCHA presents a good agreement with experimental data, even though a
slight overestimation is observed compared to the PIMD mode-coupling results,
especially for the transverse modes in the $\Gamma-U$ and $\Gamma-L$
directions.

Figure 4: Distribution of first neighbor distances in fcc 4He at 38K. The
vertical grey bar correspond to the equilibrium distance.

It is interesting to ask why, even with the strong Gaussian approximation, the
SCHA still gives an accurate description of the the quasiparticle frequencies,
while the classical MD completely disagrees with experimental results. To
answer this question, we plot in Fig. 4 the density of first neighbor
distances computed with MD, PIMD and the SCHA. This figure allows us to
highlight the effects of NQE on the system. Indeed, comparing MD and PIMD
results, the former presents a high skewness directed to large distances,
while the inclusion of zero-point motion in the latter results in a partial
re-symmetrization of the density. Because of this, the tails of the PIMD and
SCHA densities are similar, even though their areas differ, indicating a close
to symmetric yet non-Gaussian distribution for the displacements in PIMD. In
the end, the differences in exploration of the BOS strongly impact the
generalized IFC. This is illustrated in Fig. 5 which compares the symmetry-
irreducible coefficients entering the frequency matrix and 2nd order vertex in
the first shell of interactions. If the SCHA and PIMD frequency matrix
coefficients are quite close, their magnitudes are almost twice those computed
with MD, which explains the softer frequencies observed within the classical
mode-coupling theory. The difference of ensemble sampling used is even more
marked for the 2nd order vertex. Comparing the vertex coefficients for the
PIMD and SCHA cases, we observe that the harmonic perturbation theory is
coincidentally close to the full PIMD results, while the SCHA has an almost
doubled magnitude compared to PIMD. This implies that even if the SCHA has
extended validity for effective QP frequencies thanks to NQE, higher order
quantities will still be strongly biased by the Gaussian starting point.

Figure 5: Coefficients entering the a) frequency matrix and b) the 2nd order
vertex in the 1st shell of interactions computed with PIMD, MD, SCHA and
perturbation theory.

From this application, it is clear that the method used to compute the
frequency matrix and the interaction vertex significantly impact the resulting
generalized susceptibility. In future works, it will be interesting to
meticulously compare the effects of the sampling strategies in different
regimes of anharmonicity. Moreover, while we limited ourselves to the
scattering-like approximation, it remains to test the effects of the self-
consistent solution of the mode-coupling equations, especially on systems
where the QP picture is broken and the spectral function is strongly non-
Lorentzian in shape.

##  IX Conclusion

In summary, we have developed a non-perturbative theory of lattice dynamics
based on the Mori-Zwanzig projection operator formalism. We derived a
Generalized Langevin Equation for the displacement Kubo correlation functions,
resulting in an exact description of the dynamics of the system in terms of a
conservative part, driven by the frequency matrix $\mathbf{\Theta}$, and a
dissipative part, driven by the memory kernel $\mathbf{K}(t)$. From the real
space formulation of the dynamics, we were able to project this GLE into a
quasiparticle basis, thus giving an equation of motion in reciprocal space.
Using linear response theory, we then derived an exact formulation for the
vibrational spectrum, in terms of the spectral function
$\mathbf{\raisebox{2.0pt}{\mbox{\large$\chi$}}}^{\prime\prime}(\omega)$. In
order to approximate the memory kernel, we developed a systematic expansion
inspired by mode-coupling theory, in which every new order is constructed to
minimize higher-order contributions. By truncating this expansion to lowest
order, we obtain a set of self-consistent equations for the memory kernel.
Within this mode coupling approach, the only inputs needed are the _static_
KCFs, and we show how to compute them from path-integral or classical
molecular dynamics simulations.

Our Mori-Zwanzig approach to lattice dynamics, being based on the KCF, takes
into account the quantum nature of nuclei. Consequently, this derivation
demonstrates how to obtain vibrational properties of quantum systems from
PIMD. Moreover, the classical-like nature of the equations derived allows us
to justify the use of classical MD to perform vibrational spectroscopy when
quantum effects are not relevant, for instance for systems with heavy nuclei
or at high temperature. We also compare our mode-coupling scheme to more
standard perturbation theory, which allows us to derive a scattering-like
approximation to the memory kernel. Our approach gives physical meaning to the
temperature-dependent phonons which appear in numerous recent works. Moreover,
we are able to relate the mode-coupling theory to other approaches for
anharmonic lattice dynamics. Notably, we show that in the perturbative limit,
the presented theory is equivalent to the usual diagrammatic formulation to
lowest order, and that the SCHA can be seen as a Gaussian averaged
approximation of it. We also showed that the TDEP method used in conjunction
with (PI)MD is a direct application of the mode-coupling theory, provided that
static corrections to the frequencies are _not_ applied, as they would double-
count interactions. These comparisons treat all of the main theories for
anharmonic lattice dynamics in a unified formalism, and give a better vision
of their inter-relationships, strengths and limits. The Mori-Zwanzig formalism
also allows us to show that while the commonly used Lorentzian approximation
for lineshapes can be a useful approximation, it can never be an exact
description of the lattice dynamics. Finally, we apply the mode-coupling
theory on fcc 4He to demonstrate its ability to describe both strong
anharmonicity and NQE in quantum crystals, in comparison with other theories
of anharmonic lattice dynamics.

The mode-coupling theory of lattice dynamics should apply to a more general
class of strongly anharmonic and/or quantum systems, where standard methods
are known to fail. For instance, while we limited ourselves here to
crystalline materials, the formalism only requires well defined equilibrium
positions so that it would be straightforward to adapt the theory to
disordered solids such as alloys or glasses. This framework can serve as a
starting point for new formulations of related properties such as the thermal
conductivity [71, 90, 96, 97, 57], which will be the focus of future work.

###### Acknowledgements.

The authors acknowledge the Fonds de la Recherche Scientifique (FRS-FNRS
Belgium) and Fonds Wetenschappelijk Onderzoek (FWO Belgium) for EOS project
CONNECT (G.A. 40007563), and Fédération Wallonie Bruxelles and ULiege for
funding ARC project DREAMS (G.A. 21/25-11). Simulation time was awarded by the
Belgian share of EuroHPC in LUMI hosted by CSC in Finland, by PRACE on
Discoverer at SofiaTech in Bulgaria (optospin project id. 2020225411), by the
CECI (FRS-FNRS Belgium Grant No. 2.5020.11), and by the Zenobe Tier-1 of the
Fédération Wallonie-Bruxelles (Walloon Region grant agreement No. 1117545).

## References

  * Born and Huang [1954] M. Born and K. Huang,  _Dynamical theory of crystal lattices_ (Oxford University Press, 1954).
  * Dove [1993] M. T. Dove,  _Introduction to Lattice Dynamics_ (Cambridge University Press, 1993).
  * Fultz [2010] B. Fultz, Prog. in Mater. Sci. 55, 247 (2010).
  * Klein and Horton [1972] M. L. Klein and G. K. Horton, J. Low Temp. Phys. 9, 151 (1972).
  * Knoop  _et al._ [2020] F. Knoop, T. A. R. Purcell, M. Scheffler, and C. Carbogno, Phys. Rev. Materials 4, 083809 (2020).
  * Grimvall  _et al._ [2012] G. Grimvall, B. Magyari-Köpe, V. Ozoliņš, and K. A. Persson, Rev. Mod. Phys. 84, 945 (2012).
  * Allen [2019] P. B. Allen, Modern Phys. Lett. B 34, 2050025 (2019).
  * Glensk  _et al._ [2015] A. Glensk, B. Grabowski, T. Hickel, and J. Neugebauer, Phys. Rev. Lett. 114, 195901 (2015).
  * Kim  _et al._ [2020] D. S. Kim, O. Hellman, N. Shulumba, C. N. Saunders, J. Y. Y. Lin, H. L. Smith, J. E. Herriman, J. L. Niedziela, D. L. Abernathy, C. W. Li, and B. Fultz, Phys. Rev. B 102, 174311 (2020).
  * Allen [2015] P. B. Allen, Phys. Rev. B 92, 064106 (2015).
  * Maradudin and Fein [1962] A. A. Maradudin and A. E. Fein, Phys. Rev. 128, 2589 (1962).
  * Pathak [1965] K. N. Pathak, Phys. Rev. 139, A1569 (1965).
  * Cowley [1968] R. A. Cowley, Rep. Prog. Phys. 31, 123 (1968).
  * Broido  _et al._ [2005] D. A. Broido, A. Ward, and N. Mingo, Phys. Rev. B 72, 014308 (2005).
  * Togo  _et al._ [2015] A. Togo, L. Chaput, and I. Tanaka, Phys. Rev. B 91, 094306 (2015).
  * Li  _et al._ [2014a] W. Li, J. Carrete, N. A. Katcho, and N. Mingo, Comput. Phys. Commun. 185, 1747 (2014a).
  * Ravichandran and Broido [2020] N. K. Ravichandran and D. Broido, Phys. Rev. X 10, 021063 (2020).
  * Sun  _et al._ [2010] T. Sun, X. Shen, and P. B. Allen, Phys. Rev. B 82, 224304 (2010).
  * Sun and Allen [2010] T. Sun and P. B. Allen, Phys. Rev. B 82, 224305 (2010).
  * Xia  _et al._ [2020a] Y. Xia, V. I. Hegde, K. Pal, X. Hua, D. Gaines, S. Patel, J. He, M. Aykol, and C. Wolverton, Phys. Rev. X 10, 041029 (2020a).
  * Yang  _et al._ [2022a] X. Yang, T. Feng, J. Li, and X. Ruan, Phys. Rev. B 105, 115205 (2022a).
  * Koehler [1966] T. R. Koehler, Phys. Rev. Lett. 17, 89 (1966).
  * Werthamer [1970] N. R. Werthamer, Phys. Rev. B 1, 572 (1970).
  * Ravichandran and Broido [2018] N. K. Ravichandran and D. Broido, Phys. Rev. B 98, 085205 (2018).
  * Esfarjani and Liang [2020] K. Esfarjani and Y. Liang, in  _Nanoscale Energy Transport_ (IOP Publishing, 2020) pp. 7–1–7–35.
  * Souvatzis  _et al._ [2009] P. Souvatzis, O. Eriksson, M. Katsnelson, and S. Rudin, Comput. Mater. Sci. 44, 888 (2009).
  * Tadano and Tsuneyuki [2015] T. Tadano and S. Tsuneyuki, Phys. Rev. B 92, 054301 (2015).
  * Tadano and Tsuneyuki [2018] T. Tadano and S. Tsuneyuki, J. Phys. Soc. Jpn 87, 041015 (2018).
  * Bianco  _et al._ [2017] R. Bianco, I. Errea, L. Paulatto, M. Calandra, and F. Mauri, Phys. Rev. B 96, 014111 (2017).
  * Monacelli  _et al._ [2021] L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea, and F. Mauri, J. Phys.: Condens. Matter 33, 363001 (2021).
  * van Roekeghem  _et al._ [2021] A. van Roekeghem, J. Carrete, and N. Mingo, Comput. Phys. Commun. 263, 107945 (2021).
  * Choquard [1967] P. Choquard,  _The anharmonic crystal_ (WA Benjamin, 1967).
  * Samathiyakanit and Glyde [1973] V. Samathiyakanit and H. R. Glyde, J. Phys. C: Solid State Phys. 6, 1166 (1973).
  * Levy  _et al._ [1984] R. M. Levy, A. R. Srinivasan, W. K. Olson, and J. A. McCammon, Biopolymers 23, 1099 (1984).
  * Kong  _et al._ [2009] L. T. Kong, G. Bartels, C. Campañá, C. Denniston, and M. H. Müser, Comput. Phys. Commun. 180, 1004 (2009).
  * Kong [2011] L. T. Kong, Comput. Phys. Commun. 182, 2201 (2011).
  * Hellman  _et al._ [2011] O. Hellman, I. A. Abrikosov, and S. I. Simak, Phys. Rev. B 84, 180301(R) (2011).
  * Hellman  _et al._ [2013] O. Hellman, P. Steneteg, I. A. Abrikosov, and S. I. Simak, Phys. Rev. B 87, 104111 (2013).
  * Hellman and Abrikosov [2013] O. Hellman and I. A. Abrikosov, Phys. Rev. B 88, 144301 (2013).
  * Klarbring  _et al._ [2020] J. Klarbring, O. Hellman, I. A. Abrikosov, and S. I. Simak, Phys. Rev. Lett. 125, 045701 (2020).
  * Mei  _et al._ [2015] A. B. Mei, O. Hellman, N. Wireklint, C. M. Schlepütz, D. G. Sangiovanni, B. Alling, A. Rockett, L. Hultman, I. Petrov, and J. E. Greene, Phys. Rev. B 91, 054101 (2015).
  * Romero  _et al._ [2015] A. H. Romero, E. K. U. Gross, M. J. Verstraete, and O. Hellman, Phys. Rev. B 91, 214310 (2015).
  * Chaney  _et al._ [2021] D. Chaney, A. Castellano, A. Bosak, J. Bouchet, F. Bottin, B. Dorado, L. Paolasini, S. Rennie, C. Bell, R. Springell, and G. H. Lander, Phys. Rev. Materials 5, 035004 (2021).
  * Hele [2017] T. J. H. Hele, Mol. Phys. 115, 1435 (2017).
  * Schofield [1960] P. Schofield, Phys. Rev. Lett. 4, 239 (1960).
  * Egorov  _et al._ [1999] S. A. Egorov, K. F. Everitt, and J. L. Skinner, J. Phys. Chem. A 103, 9494 (1999).
  * Ramírez  _et al._ [2004] R. Ramírez, T. López-Ciudad, P. Kumar P, and D. Marx, J. Chem. Phys. 121, 3973 (2004).
  * Geng [2022] H. Y. Geng, J. Phys. Chem. C 126, 19355 (2022).
  * Morresi  _et al._ [2021] T. Morresi, L. Paulatto, R. Vuilleumier, and M. Casula, J. Chem. Phys. 154, 224108 (2021).
  * Morresi  _et al._ [2022] T. Morresi, R. Vuilleumier, and M. Casula, Phys. Rev. B 106, 054109 (2022).
  * Mori [1965] H. Mori, Prog. Theor. Phys. 33, 423 (1965).
  * Zwanzig [1961] R. Zwanzig, Phys. Rev. 124, 983 (1961).
  * Kubo [1966] R. Kubo, Rep. Prog. Phys. 29, 255 (1966).
  * Hijón  _et al._ [2010] C. Hijón, P. Español, E. Vanden-Eijnden, and R. Delgado-Buscalioni, Faraday Discuss. 144, 301 (2010).
  * Carof  _et al._ [2014] A. Carof, R. Vuilleumier, and B. Rotenberg, J. Chem. Phys. 140, 124103 (2014).
  * Li  _et al._ [2015] Z. Li, X. Bian, X. Li, and G. E. Karniadakis, J. Chem. Phys. 143, 243128 (2015).
  * Fiorentino and Baroni [2022] A. Fiorentino and S. Baroni, From green-kubo to the full boltzmann kinetic approach to heat transport in crystals and glasses (2022).
  * Bosse  _et al._ [1978] J. Bosse, W. Götze, and M. Lücke, Phys. Rev. A 17, 434 (1978).
  * Reichman and Charbonneau [2005] D. R. Reichman and P. Charbonneau, J. Stat. Mech.: Theory Exp. 2005 (05), P05013.
  * Markland  _et al._ [2012] T. E. Markland, J. A. Morrone, K. Miyazaki, B. J. Berne, D. R. Reichman, and E. Rabani, J. Chem. Phys. 136, 074511 (2012).
  * Janssen [2018] L. M. C. Janssen, Front. Phys. 6, 97 (2018).
  * Kubo  _et al._ [1991] R. Kubo, M. Toda, and N. Hashitsume,  _Statistical Physics II_ (Springer Berlin Heidelberg, 1991).
  * Craig and Manolopoulos [2004] I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 121, 3368 (2004).
  * Brown [1994] E. B. Brown, Phys. Rev. B 49, 4305 (1994).
  * Note [1] While the choice of prefactor here can seem arbitrary, we use it to simplify comparison with other methods of anharmonic lattice dynamics.
  * Braams and Manolopoulos [2006] B. J. Braams and D. E. Manolopoulos, J. Chem. Phys. 125, 124105 (2006).
  * Rossi  _et al._ [2014] M. Rossi, M. Ceriotti, and D. E. Manolopoulos, J. Chem. Phys. 140, 234116 (2014).
  * Zwanzig [2001] R. Zwanzig,  _Nonequilibrium statistical mechanics_ (Oxford university press, 2001).
  * Berne and Harp [1970] B. J. Berne and G. D. Harp, in  _Advances in Chemical Physics_ (John Wiley & Sons, Inc., 1970) pp. 63–227.
  * Szamel [2007] G. Szamel, J. Chem. Phys. 127, 084515 (2007).
  * Simoncelli  _et al._ [2019] M. Simoncelli, N. Marzari, and F. Mauri, Nat. Phys. 15, 809 (2019).
  * Eriksson  _et al._ [2019] F. Eriksson, E. Fransson, and P. Erhart, Adv. Theory Simul. 2, 1800184 (2019).
  * Zhou  _et al._ [2019] F. Zhou, W. Nielson, Y. Xia, and V. Ozoliņš, Phys. Rev. B 100, 184308 (2019).
  * Bottin  _et al._ [2020] F. Bottin, J. Bieder, and J. Bouchet, Comput. Phys. Commun. 254, 107301 (2020).
  * Ramírez and López-Ciudad [1999] R. Ramírez and T. López-Ciudad, J. Chem. Phys. 111, 3339 (1999).
  * Shulumba  _et al._ [2017] N. Shulumba, O. Hellman, and A. J. Minnich, Phys. Rev. B 95, 014302 (2017).
  * Xia  _et al._ [2020b] Y. Xia, K. Pal, J. He, V. Ozoliņš, and C. Wolverton, Phys. Rev. Lett. 124, 065901 (2020b).
  * Zhu  _et al._ [2020] Y. Zhu, Y. Xia, Y. Wang, Y. Sheng, J. Yang, C. Fu, A. Li, T. Zhu, J. Luo, C. Wolverton, G. J. Snyder, J. Liu, and W. Zhang, Res. 2020, 1 (2020).
  * Zeng  _et al._ [2021] Z. Zeng, C. Zhang, Y. Xia, Z. Fan, C. Wolverton, and Y. Chen, Phys. Rev. B 103, 224307 (2021).
  * Yang  _et al._ [2022b] X. Yang, J. Tiwari, and T. Feng, Mater. Today Phys. 24, 100689 (2022b).
  * Feng  _et al._ [2017] T. Feng, L. Lindsay, and X. Ruan, Phys. Rev. B 96, 161201(R) (2017).
  * Decoster [2004] A. Decoster, J Phys. A Math. 37, 9051 (2004).
  * Li  _et al._ [2014b] C. W. Li, O. Hellman, J. Ma, A. F. May, H. B. Cao, X. Chen, A. D. Christianson, G. Ehlers, D. J. Singh, B. C. Sales, and O. Delaire, Phys. Rev. Lett. 112, 175501 (2014b).
  * Dangić  _et al._ [2021] Đ. Dangić, O. Hellman, S. Fahy, and I. Savić, npj Comput. Mater. 7, 10.1038/s41524-021-00523-7 (2021).
  * Bouchet and Bottin [2017] J. Bouchet and F. Bottin, Phys. Rev. B 95, 054113 (2017).
  * Ladygin  _et al._ [2020] V. Ladygin, P. Korotaev, A. Yanilkin, and A. Shapeev, Comput. Mater. Sci. 172, 109333 (2020).
  * Ding  _et al._ [2021] J. Ding, T. Lanigan-Atkins, M. Calderón-Cueva, A. Banerjee, D. L. Abernathy, A. Said, A. Zevalkink, and O. Delaire, Sci. Adv. 7, 10.1126/sciadv.abg1449 (2021).
  * Xie  _et al._ [2020] L. Xie, J. H. Feng, R. Li, and J. Q. He, Phys. Rev. Lett. 125, 245901 (2020).
  * He  _et al._ [2020] X. He, D. Bansal, B. Winn, S. Chi, L. Boatner, and O. Delaire, Phys. Rev. Lett. 124, 145901 (2020).
  * Isaeva  _et al._ [2019] L. Isaeva, G. Barbalinardo, D. Donadio, and S. Baroni, Nat. Commun. 10, 3853 (2019).
  * Eckert  _et al._ [1977] J. Eckert, W. Thomlinson, and G. Shirane, Physical Review B 16, 1057 (1977).
  * Cazorla and Boronat [2017] C. Cazorla and J. Boronat, Rev. Mod. Phys. 89, 10.1103/revmodphys.89.035003 (2017).
  * Goldman  _et al._ [1968] V. Goldman, G. Horton, and M. Klein, Phys. Rev. Lett. 21, 1527 (1968).
  * Aziz  _et al._ [1979] R. A. Aziz, V. P. S. Nain, J. S. Carley, W. L. Taylor, and G. T. McConville, J. Chem. Phys. 70, 4330 (1979).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Comput. Phys. Commun. 271, 108171 (2022).
  * Simoncelli  _et al._ [2022] M. Simoncelli, N. Marzari, and F. Mauri, Phys. Rev. X 12, 041011 (2022).
  * Caldarelli  _et al._ [2022] G. Caldarelli, M. Simoncelli, N. Marzari, F. Mauri, and L. Benfatto, Phys. Rev. B 106, 024312 (2022).