The dynamical properties of nuclei, carried by the concept of phonon quasiparticles , 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 ^{4}He, an archetypal quantum crystal with very strong anharmonicity.

## 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.^{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 stabilization 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 stabilization 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 Ω_{ν}(*T*), as well as a linewidth Γ_{ν}(*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–13} to describe the effects of finite temperature on phonons within an *ab initio* setting.^{14–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 = 0 K) harmonic contribution being dominant. This hypothesis is often too drastic, and perturbative approaches are consequently known to fail in some cases.^{18–21} Moreover, a significant drawback of such approaches is that they cannot 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–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,^{23,30,32} 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.^{29,32}

Fortunately, there exists a clear way to lift the multivariate Gaussian distribution approximation. Indeed, given a potential *V*(**R**), which can be computed using *ab initio* methods, such as density functional theory (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–36} However, the prevailing MD-based lattice dynamics method is probably the temperature dependent effective potential (TDEP).^{37–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–39,41–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–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 does not 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–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–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 (NQE) 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 [(Path-integral) Molecular Dynamics], 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.

This paper is organized as follows. In Sec. 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 Sec. III, before using the GLE to derive exact equations describing the spectral properties in Sec. IV. Using a mode-coupling approach, we then derive in Sec. 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 Sec. 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 ^{4}He in Sec. VIII before concluding in Sec. IX.

## II. DERIVATION OF THE GENERALIZED LANGEVIN EQUATION

### A. Notations

**R**and

**p**are, respectively, the position and momentum operators, and

*V*(

**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

**u**

_{i}=

**R**

_{i}− ⟨

**R**

_{i}⟩ as dynamical variables. In this definition, ⟨

*O*(

*t*)⟩ is the thermal average of the operator

*O*, meaning that ⟨

**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.

*A*=

*A*(0) and $L=1\u210fH,\u22c5$ is the quantum Liouville operator. We will denote $pi(t)=Miu\u0307i(t)$ the time derivative of the displacement of atom

*i*with mass

*M*

_{i}. Furthermore, we will define mass scaled coordinates $u\u0303i=Miui$ and $p\u0303i=1Mipi$.

### B. Kubo correlation function

*A*and

*B*, is defined as

^{53,62}has the advantage of being easily approximated in the path-integral molecular dynamics formalism.

^{49,63}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}

^{,}

^{47,53}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

*ℏ*→ 0.

*N*

_{at}× 3

*N*

_{at}symmetric matrix defined with

### C. The generalized Langevin equation

^{51}and Zwanzig

^{52}allows to express Eq. (7) in a form amenable to approximations. This formalism is based on operators $P$ and $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

^{55}which results in the properties

*i*,

^{55,59}to expand the $Q$ projection.

**Θ**(note that it is not in mass scaled coordinates)

*i*, as well as

*δ*

**f**(

*t*), which is often called the random force,

**K**(

*t*)

^{58–61}

^{,}

^{55}though this is not the case.

^{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.

**D**(

**q**) as the Fourier transform of the mass weighted

**Θ**,

*i*→ (

*τ*,

*i*) so that $\Theta \delta ,i\tau ,j$ is the frequency matrix relating atom

*i*in the unit cell

*δ*to the atom

*j*in unit cell

*τ*. Similarly to the harmonic phonon case, we can now obtain QP polarization vectors

*ϵ*^{ν}(

**q**) and frequencies Ω

_{q,ν}by diagonalizing

**D**(

**q**),

*ν*is a QP mode. Without loss of generality, we can now express the mass-weighted atomic displacements as

*A*

_{λ}(

**q**,

*t*) are unitless quasiparticle operators with the property $A\lambda (q),A\lambda \u2032(q\u2032)=2kBT\u210f\Omega \lambda (q)\delta \lambda \lambda \u2032\delta q,q\u2032$

^{65}and

*ρ*

_{q,λ}is a prefactor determined below. Using the static solution of Eq. (17) as well as the property $G\u0308ij(0)=\u2212kBT/Mi\delta ij$,

^{66,67}we can express the static mass-weighted displacement covariance tensor as

*i*can be expressed as

## IV. VIBRATIONAL SPECTRA FROM THE KUBO CORRELATION FUNCTION

^{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 Φ

_{AB}, which we define as

*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

_{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 responsefunction

*A*at

*t*= 0.

*A*and

*B*, and $G\u0307\u0303AB(i\omega )$ is the Laplace transform of its time derivative. The dynamical susceptibility can be decomposed into a real and imaginary part $\chi AB(\omega )=\chi AB\u2032(\omega )+i\chi AB\u2032\u2032(\omega )$, where $\chi AB\u2032(\omega )$ and $\chi AB\u2032\u2032(\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 $G\u0303AB(\omega )$. Since the KCF is an even and real function of time, this real part is proportional to the power spectrum

*G*

_{AB}(

*ω*), 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

**χ**″(

**q**,

*ω*) is particularly important since this quantity can be linked to the cross section

*d*

^{2}

*σ*/

*d*Ω

*dE*from spectroscopy measurements, allowing thus a direct comparison between theory and experiment. For the QP KCF, using the initial conditions [Eqs. (30) and (31)], the Laplace transform is expressed as

**G**(

**q**,

*ω*), we start by separating the memory kernel into a real and an imaginary part as $K\u0303(q,i\omega )=\Gamma \u0303(q,\omega )\u2212i\Delta \u0303(q,\omega )$. Since

**K**(

**q**,

*t*) is a real and even function of time, the real part $\Gamma \u0303(q,\omega )$ of its Laplace transform is proportional to its Fourier transform

**Γ**(

**q**,

*ω*), with $\Gamma \u0303(q,\omega )=2\Gamma (q,\omega )$. Moreover, $\Gamma \u0303(q,\omega )$ and $\Delta \u0303(q,\omega )$ are not independent, and one can recover the other using Kramers–Kronig relations.

^{69}For instance, $\Delta \u0303(q,\omega )$ can be computed with

In the end, the dynamical response of the system can be entirely characterized using **Ω**(**q**), **Γ**(**q**, *ω*), and **Δ**(**q**, *ω*). 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

**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 $Pn$,

^{70}which project an operator

*O*(

*t*) on the subspace of

*n*QP operators, as well as the corresponding orthogonal projection $Qn$,

*n*! is here to account for all the permutations of

*A*

_{λ}operators so that $PnPn=Pn$. Going up to order

*N*and inserting several times the identity $Pn+Qn=1$, we can rewrite the random forces

*δ*

_{n}

*A*

_{λ}, and

*δ*

_{N}

*A*

_{λ}is the final residual random force we will neglect. We can now inject this expression into the memory kernel to obtain

*n*,

*n*′) order memory kernel

*n*, the $Qn$ projectors minimize the corresponding contribution of the random forces

*δ*

_{n}

*A*

_{λ}. 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

**K**(

**q**,

*t*).

*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 two point correlations in real time.

^{58,59,61}For instance, a four point correlation function is given by

^{58,59}

*n*≠

*n*′,

*λ*≠

*λ*′). 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 −

*λ*), our approximation for the memory kernel is given by

*δ*

_{G}(

**q**

_{λ}+

**q**

_{μ}+

**q**

_{ν}) = 1 if

**q**

_{λ}+

**q**

_{μ}+

**q**

_{ν}is 0 modulo a reciprocal lattice vector

**G**.

**K**, we now have all the ingredients to compute the power spectrum. Fourier transforming Eq. (49) in frequency domain, and using Eqs. (38) and (39), we obtain the following set of self-consistent equations:

*G*

_{λ,0}(

*ω*), which allows to obtain the linewidth

*s*,

*s*′ = ±1. From this result, a first approximation to the power spectrum can be obtained by using Eqs. (53) and (54) with $\Gamma \lambda s(\omega )$ as an input.

## VI. COMPUTING THE FREQUENCIES AND VERTEX FROM (PI)MD SIMULATIONS

*A*and

*B*can be computed exactly in the path-integral molecular dynamics framework as

^{49,63}

*P*is the number of beads in the classical polymer, and $A\u0304(t)=\u2211pAp(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.

*P*= 1 and can be computed as

*N*of configurations, the frequency matrix can be computed with

**Θ**of the linear problem

**Θu**= −

**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

**Θ**might necessitate long simulations, which could become prohibitive. Fortunately,

**Θ**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,

**Θ**follows the usual symmetries found in the IFC tensor of the harmonic approximation for phonons. Consequently, imposing them on

**Θ**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–74}A similar analysis can be done for the second order vertex in real space, which can be computed as

*δf*

_{i}=

*f*

_{i}+

*∑*

_{j}Θ

_{ij}

*u*

_{j}. A simple variable transformation allows to rewrite this equation as a least-squares solution of the equation

**Θ**

^{T}

**U**=

*δ*

**f**, where

**U**=

**uu**

^{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 $Pn$.

## VII. DISCUSSION

### A. Classical-like correlation functions and classical limit

*j*of energy

*E*

_{j}and

*n*(

*ω*) = 1/(

*e*

^{βℏω}− 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 Eq. (52)–(54). Therefore, the only difference between the quantum and classical description comes from the type of correlation functions used to compute the frequency matrix **Θ** 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.

*n*

_{μ}=

*n*(Ω

_{μ}), and where we used lim

_{ω→0}1/[

*ℏωn*(

*ω*)] =

*β*. Once Ψ 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.

### B. 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.^{49,75} While such an approximation would be accurate at low temperature, it does not hold for high temperature and completely misses the exact character associated with the quantity.^{75}

*ω*→ 0 limit of Eq. (35), the static susceptibility can be written as

*χ*= −

*βG*(0). From Eq. (13) and using the property $ui,fj=\u2212kBT\delta ij$, we obtain a proportionality relation between the frequency matrix and the inverse static susceptibility in real space

*ϵ*_{λ}(

**q**) are also eigenvectors of the Fourier transform of the susceptibility matrix, with eigenvalues $\u210f\Omega \lambda \u22122(q)$. From this observation, we can interpret the eigenvectors

*ϵ*^{λ}(

**q**) of the quasiparticles as displacement patterns associated with a perturbation with a momentum

**q**. The

*A*

_{λ}(

**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.

*G*

_{λ}(

**q**) and $G\u0308\lambda (q)$,

*ω*

_{jk}= (

*E*

_{j}−

*E*

_{k})/

*ℏ*, and where we used $\u3008j|A\u0308\lambda (q)|k\u3009=\u2212\omega jk2\u3008j|A\lambda (q)|k\u3009$, we obtain a formulation of the frequencies in terms of a complete set of eigenstates of the Hamiltonian

_{λ}.

### C. 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 I.

. | Perturbation theory . | SCHA . | Mori–Zwanzig . |
---|---|---|---|

Sampling | Perturbative through finite | Effective harmonic | Exact through MD (classical case) |

of displacements | difference or DFPT | (multivariate Gaussian) | or PIMD (quantum case) |

Static quasiparticles | Harmonic phonons | Effective harmonic phonons | Eigenstates of the inverse |

or eigenstates of the inverse susceptibility | susceptibility | ||

Vertex in the bubble | Perturbative | Gaussian approximation | Exact |

contribution | |||

Implementations | Phono3py^{15} | SSCHA^{30} | MD-TDEP^{37,38} |

Alamode^{27,28} | |||

ShengBTE^{16} | sTDEP^{76} | ||

qSCAILD^{31} |

. | Perturbation theory . | SCHA . | Mori–Zwanzig . |
---|---|---|---|

Sampling | Perturbative through finite | Effective harmonic | Exact through MD (classical case) |

of displacements | difference or DFPT | (multivariate Gaussian) | or PIMD (quantum case) |

Static quasiparticles | Harmonic phonons | Effective harmonic phonons | Eigenstates of the inverse |

or eigenstates of the inverse susceptibility | susceptibility | ||

Vertex in the bubble | Perturbative | Gaussian approximation | Exact |

contribution | |||

Implementations | Phono3py^{15} | SSCHA^{30} | MD-TDEP^{37,38} |

Alamode^{27,28} | |||

ShengBTE^{16} | sTDEP^{76} | ||

qSCAILD^{31} |

#### 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.

#### 2. Perturbation theory

^{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

**Θ**

^{pt}and the memory kernel can be written

*λ*,

*μ*,

*ν*). 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

*κ*∝

*T*

^{−α}, as shown in Refs. 78–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 Sec. II C 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 Ω_{λ}(**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.

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.^{17,21,81} On the contrary, the projection operators of Eqs. (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.

#### 3. Self-consistent harmonic approximation

*O*⟩

_{0}is an average computed within the distribution associated with $V\u0303(R)$. As the name of the method implies, the effective potential has a harmonic form

**Θ**

^{SCHA}and equilibrium positions $R$ are parameters to optimize. The advantage of such a potential is that the distribution of displacements around $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 $F<F\u0303$

^{82}so that optimizing the parameters means minimizing $F\u0303$. For the effective force constants, the minimum is obtained when

^{26–31,76}

**Ξ**of the free energy

^{29,30}

**Ξ**can be understood as the inverse static susceptibility. However, compared to the mode-coupling approach, where this quantity is exact,

*χ*^{−1}is approximated in the SCHA, as only the contribution from the first-order free energy cumulant enters in $F\u0303$.

^{29,30}in which the vertices are computed using the average of the third-order derivative of the potential energy

**Θ**

^{SCHA}. With this assumption, the frequency matrix is actually given by the SCHA force constants

**Θ**

^{SCHA}, and the random forces are given with by $\delta fi=fi+\u2211j\Theta ijSCHAuj$. Applying the mode-coupling decomposition of Sec. V to approximate

**K**

^{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.

#### 4. Average of the Hessian

^{49}

_{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 Π

_{ij}= Θ

_{ij}that only holds in a harmonic system.

^{49}

^{,}

**K**is a positive definite matrix,

^{69}this equation tells us that

**Π**is greater that

**Θ**, meaning that frequencies computed from the diagonalization of the average Hessian will overestimate the Ω

_{λ}. 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.

#### 5. Temperature dependent effective potential

_{ij}are fit using ordinary least-squares. Then, the third order force constants Ψ

_{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} which 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 Sec. VI, the force constants **Θ**^{TDEP} and **Ψ**^{TDEP} computed from the TDEP methods are equal to the frequency matrix **Θ** and the real space vertex **Ψ**. 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.^{40,87–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.

### D. Markovian approximation and Lorentzian lineshape

^{7,15,90}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 Γ(

*ω*) is approximated with a constant given by Γ

_{λ}= Γ

_{λ}(Ω

_{λ}).

^{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

**K**(

*t*) decays very quickly compared to

**G**(

*t*), the memory kernel can be replaced by a constant friction term $\Gamma M(q)=\beta \u222b0\u221edtK(q,t)$ so that the Kubo transformed Green’s function is written as

Nevertheless, for mildly 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.

## VIII. APPLICATIONS

To demonstrate the accuracy of the theory developed here, we apply it to the fcc phase of ^{4}He. 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., $\u27e8u2\u27e9/a$ with ⟨*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 mol^{3}/mol, which correspond to a lattice parameter of $3.915$ Å.

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 Large-scale Atomic/Molecular Massively Parallel Simulator (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 Sec. 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 A. For this step, a 12 × 12 × 12 **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 C 5 with Bose–Einstein occupations in order to capture NQE.

^{4}He, 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

*σ*

^{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

*σ*

^{A}as a measure of the dissipative contribution to the dynamics of the system. For fcc

^{4}He, we find

*σ*

^{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

*σ*

^{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, 0 K 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 Γ − *U* and Γ − *L* directions.

It is interesting to ask why, even with the strong Gaussian approximation, the SCHA still gives an accurate description of 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 second 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 second 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.

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 **Θ**, and a dissipative part, driven by the memory kernel **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 ** χ**″(

*ω*). 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 that 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 ^{4}He 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,^{57,71,90,96,97} which will be the focus of future work.

## ACKNOWLEDGMENTS

The authors acknowledge the Fonds de la Recherche Scientifique (FRS-FNRS Belgium) and Fonds Wetenschappelijk Onderzoek (FWO Belgium) for EOS project CONNECT (Grant No. G.A. 40007563), and Fédération Wallonie Bruxelles and ULiege for funding ARC project DREAMS (Grant No. G.A. 21/25-11).

Simulation time was awarded by the Belgian share of EuroHPC in LUMI hosted by CSC in Finland, PRACE on Discoverer at SofiaTech in Bulgaria (optospin Project No. 2020225411), the CECI (FRS-FNRS Belgium, Grant No. 2.5020.11), and the Zenobe Tier-1 of the Fédération Wallonie-Bruxelles (Walloon Region, Grant Agreement No. 1117545).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Aloïs Castellano**: Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). **J. P. Alvarinhas Batista**: Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). **Matthieu J. Verstraete**: Methodology (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in github at https://github.com/OrbitalC/ModeCouplingTheoryPhononsData.

## REFERENCES

*Dynamical Theory of Crystal Lattices*

*Introduction to Lattice Dynamics*

*Nanoscale Energy Transport*

*Statistical Physics II*

While the choice of prefactor here can seem arbitrary, we use it to simplify comparison with other methods of anharmonic lattice dynamics.

*Nonequilibrium Statistical Mechanics*

*Advances in Chemical Physics*