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 4He, an archetypal quantum crystal with very strong anharmonicity.

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 theory11–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 theory58–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 4He in Sec. VIII before concluding in Sec. IX.

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=ipi22Mi+V(R),
(1)
where 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 ui = Ri − ⟨Ri⟩ 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.
To describe the time evolution of the system, we use the Heisenberg representation, where the time dependence and derivatives of operators are obtained with
A(t)=eiHtAeiHt=eiLtA,Ȧ(t)=iH,A(t)=iLA(t),
(2)
where A = A(0) and L=1H, is the quantum Liouville operator. We will denote pi(t)=Miu̇i(t) the time derivative of the displacement of atom i with mass Mi. Furthermore, we will define mass scaled coordinates ũi=Miui and p̃i=1Mipi.
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
GAB(t)=kBT0βdλA(iλ)B(t)=kBTZ0βdλTre(βλ)HAeλHB(t)=A,B(t),
(3)
where β=1kBT is the inverse temperature and 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.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 GAB(t) = GAB(−t) as well as GAB(t) = GBA(t). Using the Lehmann representation, it can also be shown that derivatives of the static Kubo correlation function are related64,
dkdtkA,B(t)|t=0=(1)kdkdtkB,A(t)|t=0,
(4)
dkdtkA,A(t)|t=0=0k=1,3,5,.
(5)
Naturally, all the different quantum correlation functions are related by simple relations47,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.
In this work, the KCF of interest is the mass-weighted displacement KCF, which is a 3Nat × 3Nat symmetric matrix defined with
Gij(t)=MiMjkBTZ0βdλTre(βλ)HuieλHuj(t)=MiMjui,uj(t)=ũi,ũj(t)
(6)
with an equation of motion given by
G̈ij(t)=MiMjui,üj(t)=ũi,p̃̇j(t).
(7)
Our goal in this work is to derive a solution for this equation of motion.
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 Mori51 and Zwanzig52 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
PO(t)=ijũj,O(t)ũi,ũjũi+ijp̃j,O(t)p̃i,p̃jp̃i,
(8)
QO(t)=(1P)O(t).
(9)
By construction, it should be noted that the projectors P and Q are orthogonal, in the sense that they follow P+Q=1 and PQ=QP=0. Being projection operators, it is natural that P and Q are idempotent so that PP=P and QQ=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
X,PY=PX,Y,
(10)
X,QY=QX,Y.
(11)
Equipped with these operators, we can now express the equation of motion of the mass scaled momentum of atom i,
p̃̇i(t)=eiLtiLp̃i=eiLt(P+Q)iLp̃i=eiLtPiLp̃i+eiQLtQiLp̃i+0tdseiL(ts)PiLeiQLsQiLp̃i,
(12)
where we used the Dyson identity55,59 to expand the Q projection.
To simplify the notation of Eq. (12), we introduce the frequency matrix Θ (note that it is not in mass scaled coordinates)
Θij=kuk,fiuj,uk,
(13)
where fi=MiL2ui is the force acting on atom i, as well as δf(t), which is often called the random force,
δfi(t)=eiQLtfi+jΘijuj
(14)
and the memory kernel K(t)
Kij(t)=δfi,δfj(t)MiMj
(15)
so that the equations of motion for the momentum can be written as58–61,
p̃̇i(t)=jΘijMiMjũj(t)δfi(t)Miβj0tdsKij(s)p̃j(ts).
(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 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 eiQLt, and its orthogonal nature is recovered thanks to the property ũi,δfi(t)=0. It should be noted that the randomness of this term is only apparent. By contracting information through the projector Q, the time evolution of this term is more difficult to apprehend and its dynamics appears “random,”55 though this is not the case.
FIG. 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 P, represented as the red arrows, projects this dynamic onto a “slow” subspace depicted here as a blue rectangle. By construction, the starting point of the dynamics, visualized with the green dot, is located on the projected subspace.

FIG. 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 P, represented as the red arrows, projects this dynamic onto a “slow” subspace depicted here as a blue rectangle. By construction, the starting point of the dynamics, visualized with the green dot, is located on the projected subspace.

Close modal
We can now inject 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
G̈(t)=ΘMTMG(t)β0tdsK(s)Ġ(ts).
(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.

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 D(q) as the Fourier transform of the mass weighted Θ,
Dij(q)=τΘ0,iτ,jMiMjeiq(Rτ,jR0,i),
(18)
where we introduced a unit cell based notation i → (τ, i) so that Θδ,iτ,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),
D(q)ϵν(q)=Ων2(q)ϵν(q),
(19)
where ν is a QP mode. Without loss of generality, we can now express the mass-weighted atomic displacements as
ũi(t)=q,λρq,λϵiλ(q)Aλ(q,t),
(20)
where Aλ(q, t) are unitless quasiparticle operators with the property Aλ(q),Aλ(q)=2kBTΩλ(q)δλλδq,q65 and ρq,λ is a prefactor determined below. Using the static solution of Eq. (17) as well as the property G̈ij(0)=kBT/Miδij,66,67 we can express the static mass-weighted displacement covariance tensor as
ũi,ũj=kBTMiMjΘij1.
(21)
Injecting Eq. (20) into this equation and with some reorganization, it can be shown that ρq,λ=/Ωλ(q) so that the displacements can be expressed as
ui(t)=2Miq,λϵiλ(q)Ωλ(q)Aλ(q,t).
(22)
We can now define a Kubo transformed QP Green’s function
Gλ,μ(q,t)=Aλ(q),Aμ(q,t).
(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
fi=Mi2qλϵiλ(q)Ωλ(q)Äλ(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
jβΘijαβMiMjujβ=Mi2qλϵiλ(q)Ωλ(q)Ωλ2(q)Aλ(q).
(25)
Consequently, we can define the random forces in terms of QP operators as
δfi=Mi2q,λϵiλ(q)Ωλ(q)Äλ(q)Ωλ2(q)Aλ(q)=Mi2q,λϵiλ(q)Ωλ(q)δAλ(q),
(26)
where we defined the QP random forces operator δAλ(q)=Äλ(q)Ωλ2(q)Aλ(q), which follows orthogonal time dynamics
δAλ(q,t)=eiQLtδAλ(q).
(27)
We are now able to define the memory kernel in reciprocal space
K(q,t)=kBTΩ(q)1δA(q),δA(q,t),
(28)
which allows us to write the GLE for the Kubo transformed Green’s function as
G̈(q,t)=Ω2(q)G(q,t)0tdsK(q,s)Ġ(q,ts).
(29)
The dynamical properties of the many-body system can now be obtained by solving Eq. (29) with the initial conditions
Gλμ(q,0)=2kBTΩλ(q)δλ,μ,
(30)
Ġλμ(q,0)=0.
(31)
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 ΦAB, which we define as
ΦAB(t)=i[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
χAB(ω)=0dtΦAB(t)eiω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 Φ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
ΦAB(t)=i1ZTreβHAB(t)TreβHB(t)A=i1ZTreβHeβHAeβHB(t)TreβHAB(t)=0βdλ1ZTreβHȦ(iλ)B(t)=βĠAB(t),
(34)
where Ȧ=i[H,A]/ 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
χAB(ω)=β0dtȦ,B(t)eiωt=βĠ̃AB(iω)=iβωG̃AB(iω)βGAB(0),
(35)
with G̃AB(iω) the Laplace transform of the KCF of A and B, and Ġ̃AB(iω) is the Laplace transform of its time derivative. The dynamical susceptibility can be decomposed into a real and imaginary part χAB(ω)=χAB(ω)+iχAB(ω), where χAB(ω) and χAB(ω) 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̃AB(ω). Since the KCF is an even and real function of time, this real part is proportional to the power spectrum GAB(ω), which corresponds to the Fourier transform of GAB(t), so that the dissipative part of the response, also called the spectral function, can be written as
χAB(ω)=ω2kBTGAB(ω).
(36)
This equation is simply the fluctuation–dissipation theorem, written using the Kubo transformed correlation function. The trace of the dissipation χ″(q, ω) is particularly important since this quantity can be linked to the cross section d2σ/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,iω)=2kBTΩ(q)iωK̃(q,iω)ω2Ω2(q)iωK̃(q,iω).
(37)
To obtain a formula for the Fourier transform of G(q, ω), we start by separating the memory kernel into a real and an imaginary part as K̃(q,iω)=Γ̃(q,ω)iΔ̃(q,ω). Since K(q, t) is a real and even function of time, the real part Γ̃(q,ω) of its Laplace transform is proportional to its Fourier transform Γ(q, ω), with Γ̃(q,ω)=2Γ(q,ω). Moreover, Γ̃(q,ω) and Δ̃(q,ω) are not independent, and one can recover the other using Kramers–Kronig relations.69 For instance, Δ̃(q,ω) can be computed with
Δ̃(q,ω)=1πdωΓ̃(q,ω)ωω,
(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
G(q,ω)=kBTπ8Ω(q)Γ(q,ω)(ω2Ω2(q)2ωΔ(q,ω))2+4ω2Γ2(q,ω).
(39)
Finally, the dissipative part of the response function can be computed using the fluctuation–dissipation theorem of Eq. (36), giving
χ(q,ω)=1π4ωΩ(q)Γ(q,ω)(ω2Ω2(q)2ωΔ(q,ω))2+4ω2Γ2(q,ω).
(40)

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.

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 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 nth 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,
PnO(t)=1n!λ1λnλ1λnAλ1Aλn,O(t)Aλ1Aλn,Aλ1AλnAλ1Aλn,
(41)
Qn=1Pn,
(42)
where the prefactor 1/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
δAλ=(P2+Q2)δAλ=P2δAλ+(P3+Q3)Q2δAλ=nNPnδn1Aλ+δNAλ,
(43)
δNAλ=QNQN1Q2δAλ,
(44)
where each successive application of Qn will produce an ever smaller δnAλ, and δNAλ is the final residual random force we will neglect. We can now inject this expression into the memory kernel to obtain
Kλ,λ(t)=δAλ,δAλ(t)=n=2Nn=2NPnδn1Aλ,Pnδn1Aλ(t)+δNAλ,δNAλ(t)=n=2Nn=2NKnnλ,λ(t)+δNAλ,δNAλ(t),
(45)
where we define the (n, n′) order memory kernel
Kn,nλ,λ(t)=Pnδn1Aλ,Pnδn1Aλ(t)=λ1,,λnλ1,,λnΨn(λ,λ1,,λn)Ψn(λ,λ1,,λn)×Aλ1Aλn,eiQLtAλ1Aλn
(46)
as well as the vertex
Ψn(λ1,,λn)=1n!μ1,,μnAμ1Aμn,δn1AλAλ1Aλn,Aμ1Aμn.
(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 Qn projectors minimize the corresponding contribution of the random forces δnAλ. 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).
As a first approximation, we truncate the expansion of the memory kernel at some order N and we neglect the Nth 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
Aλ1Aλ2,eiQLtAλ3Aλ4Aλ1,Aλ3(t)Aλ2,Aλ4(t)+Aλ1,Aλ4(t)Aλ2,Aλ3(t).
(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 (nn′, λλ′). 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
Kλ(t)2kBTΩλμ,ν|Ψ(λ,μ,ν)|2Gμ(t)Gν(t).
(49)
The only inputs are contained in the vertex Ψ. We define the third-order tensor Θ as
Θijk=12l,mulum,δfil,mujuk,ulum.
(50)
Using the projection of displacements and random forces on QP, the vertex can then be expressed as
Ψ(λ,μ,ν)=2ΩλΩμΩνi,j,kϵiλϵjμϵkνMiMjMkΘijk×ei(Riqλ+Rjqμ+Rkqν)δG(qλ+qμ+qν),
(51)
where δG(qλ + qμ + qν) = 1 if qλ + qμ + qν is 0 modulo a reciprocal lattice vector G.
With the preceding approximation of 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:
Γλ(ω)=2kBTΩλμ,ν|Ψ(λ,μ,ν)|2(Gμ*Gν)(ω),
(52)
Δλ(ω)=1πdωΓλ(ω)ωω,
(53)
Gλ(ω)=kBTπ8ΩλΓλ(ω)(ω2Ωλ22ωΔλ(ω))2+4ω2Γλ2(ω).
(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
G̈λ,0(t)=Ωλ2Gλ,0(t),
(55)
which is simply given by Gλ,0(t)=kBTΩλcos(Ωλt). The scattering theory approximation is obtained by replacing the Green’s function in Eq. (52) by Gλ,0(ω), which allows to obtain the linewidth
Γλs(ω)=πkBTΩλΩμΩνμ,ν,s,s|Ψ(λ,μ,ν)|2δ(ω+sΩμ+sΩν)
(56)
with s, s′ = ±1. From this result, a first approximation to the power spectrum can be obtained by using Eqs. (53) and (54) with Γλs(ω) as an input.
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 as49,63
GAB(0)=limPlimt0teβH(t)ZpPAp(t)pPBp(t)dt1NnNĀ(tn)B̄(tn),
(57)
where P is the number of beads in the classical polymer, and Ā(t)=pAp(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
GAB(0)1NnNA(tn)B(tn).
(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
Θij=nNkūk(tn)f̄i(tn)nNkūk(tn)ūj(tn),
(59)
where we used the property ui,fj=kBTδij to decouple the sums in the numerator and denominator of Eq. (13). In this equation, one can recognize the least squares solution Θ 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
Θijk=12nNl,mūl(tn)ūm(tn)δf̄i(tn)l,mūl(tn)ūm(tn)ūj(tn)ūk(tn)
(60)
with δfi = fi + jΘijuj. A simple variable transformation allows to rewrite this equation as a least-squares solution of the equation ΘTU = δf, where U = uuT.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 nth order projectors Pn.
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λ(ω)=kBTjkpj|j|Aλ|k|2ω(n(ω)+1)δωEjEk,
(61)
where pj=eβEj/Z is the Boltzmann weight associated with state j of energy Ej 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.

Furthermore, our approach allows for a simple inclusion of quantum behavior when doing classical simulations. For this, we introduce the lesser Green’s function Gλ<(t)=Aλ(t)Aλ, given in the Lehmann representation by
Gλ<(ω)=jkpjeβω|j|Aλ|k|2δωEjEk.
(62)
From this representation, it is easy to see that the Kubo and lesser Green’s functions are related by Gλ<(ω)=ωkBTn(ω)Gλ(ω). Replacing the Kubo Green’s functions in Eq. (52) with the lesser Green’s functions makes the Bose–Einstein distributions visible in the following equation:
Γλ(ω)=2kBTΩλμ,ν|Ψ(λ,μ,ν)|2×dωkBT(ωω)n(ωω)Gλ<(ωω)×kBTωn(ω)Gλ<(ω).
(63)
In the scattering theory approximation, Gλ<(t)=n(Ωλ)cos(Ωλt) so that Γλs(ω) reduces to
Γλs(ω)=π4kBTΩλμ,ν|Ψ(λ,μ,ν)|2F(ω,μ,ν),
(64)
F(ω,μ,ν)=s=±1nμ+nν+1δ(ω+sΩμ+sΩν)+nμnνδ(ω+sΩμsΩν)
(65)
with nμ = nμ), and where we used limω→01/[ℏω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.

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 

Taking the ω → 0 limit of Eq. (35), the static susceptibility can be written as χ = −βG(0). From Eq. (13) and using the property ui,fj=kBTδij, we obtain a proportionality relation between the frequency matrix and the inverse static susceptibility in real space
ΘMTM=χ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 Sec. III, and due to the Hermitian character of the matrices involved, we can conclude that the quasiparticle polarizations ϵλ(q) are also eigenvectors of the Fourier transform of the susceptibility matrix, with eigenvalues Ωλ2(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.
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 G̈ij=kBT, and using Eqs. (22) and (24), we have G̈λ(q)/Gλ(q)=Ωλ2(q). Using then the Lehmann representation of the KCF Gλ(q) and G̈λ(q),
Gλ(q)=kBTjk|j|Aλ(q)|k|2(pkpj)ωjk,
(67)
G̈λ(q)=kBTjkj|Aλ(q)|kk|Äλ(q)|j(pkpj)ωjk=kBTjkωjk2|j|Aλ(q)|k|2(pkpj)ωjk
(68)
with ωjk = (EjEk)/, and where we used j|Äλ(q)|k=ωjk2j|Aλ(q)|k, we obtain a formulation of the frequencies in terms of a complete set of eigenstates of the Hamiltonian
cjk(λ,q)=|j|Aλ(q)|k|2pkpjωjk,
(69)
Ωλ(q)=jkωjk2cjk(λ,q)jkcjk(λ,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
0dωχλ(q,ω)ω=Ωλ(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 Ωλ.

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.

TABLE I.

Summary of different approach of anharmonic lattice dynamics.

Perturbation theorySCHAMori–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 Phono3py15  SSCHA30  MD-TDEP37,38 
Alamode27,28 
ShengBTE16  sTDEP76  
qSCAILD31  
Perturbation theorySCHAMori–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 Phono3py15  SSCHA30  MD-TDEP37,38 
Alamode27,28 
ShengBTE16  sTDEP76  
qSCAILD31  

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

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
V(R)12ijΘijptuiuj+13!ijkΨijkptuiujuk,
(72)
Θijpt=2V(R)RiRj,
(73)
Ψijkpt=3V(R)RiRjRk,
(74)
where we assume that contributions from Ψijkpt are perturbations with respect to Θijpt 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 Θpt and the memory kernel can be written
Kijpt(t)=klmnΨiklptΨjmnptukul,eiQLtumun,
(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 Ψ(λμ, ν). 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.

FIG. 2.

Feynman diagrams appearing in the self-energy in perturbation theory to lowest order. Each line represents 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.

FIG. 2.

Feynman diagrams appearing in the self-energy in perturbation theory to lowest order. Each line represents 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.

Close modal

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

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 F̃ defined as
F̃=F̃0+V(R)Ṽ(R)0,
(76)
where Ṽ(R) is an effective potential, F̃0 is its associated free energy, and ⟨O0 is an average computed within the distribution associated with Ṽ(R). As the name of the method implies, the effective potential has a harmonic form
Ṽ(R)=(RR)TΘSCHA(RR),
(77)
where the effective force constants Θ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̃82 so that optimizing the parameters means minimizing F̃. For the effective force constants, the minimum is obtained when
ΘijSCHA=2V(R)uiuj0,
(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–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 Ξ of the free energy29,30
Ξij=2F̃RiRj.
(79)
By seeing the Gibbs–Bogoliubov free energy as a Landau expansion with an order parameter R, the free energy Hessian Ξ 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̃.
On top of the quasiparticles computed from the free energy Hessian, the dynamical susceptibility can be obtained with a bubble-like contribution29,30 in which the vertices are computed using the average of the third-order derivative of the potential energy
ΨijkSCHA=1MiMjMk3V(R)RiRjRk0=1MiMjMklmulumfk0uiul0ujum0,
(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 ΘSCHA. With this assumption, the frequency matrix is actually given by the SCHA force constants ΘSCHA, and the random forces are given with by δfi=fi+jΘijSCHAuj. Applying the mode-coupling decomposition of Sec. V to approximate KSCHA(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

In the literature, it has been proposed to define a generalization of the harmonic force constants through an average of the potential energy Hessian49 
Πij=2V(R)RiRj,
(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 Π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.
To highlight the differences between this formulation and the mode-coupling theory, the potential energy Hessian average can be expanded in terms of force–force KCF,49,
Πij=fi,fjkBT.
(82)
Using this equation, the static part of memory kernel can be formulated as
Kij(0)=ΠijΘij.
(83)
Since 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

TDEP is a method based on the construction of an effective anharmonic potential of the form (here presented up to third order)
Ṽ(R)=12ijΘijTDEPuiuj+13!ijkΨijkTDEPuiujuk.
(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 Θ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 GeTe84 or the phase diagram of uranium85,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.

In the limit of small anharmonicity, the spectral lineshapes are often approximated with Lorentzians.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 ΓM(q)=β0dtK(q,t) so that the Kubo transformed Green’s function is written as
G̈(q,t)=Ω2(q)G(q,t)ΓM(q)Ġ(q,t)
(85)
with solutions
G(q,t)=cosΩ(q)teΓM(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 dG(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 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.

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., u2/a with ⟨u2⟩ 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 Å.

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.

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 anharmonicity5 to our Mori–Zwanzig formalism as
σA=iδfi2ifi2.
(87)
We note, however, that used in this way σ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 4He, 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.

FIG. 3.

Spectral function χ″(ω) 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 

FIG. 3.

Spectral function χ″(ω) 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 

Close modal

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.

FIG. 4.

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

FIG. 4.

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

Close modal
FIG. 5.

Coefficients entering the (a) frequency matrix and (b) the second order vertex in the first shell of interactions computed with PIMD, MD, SCHA, and perturbation theory.

FIG. 5.

Coefficients entering the (a) frequency matrix and (b) the second order vertex in the first shell of interactions computed with PIMD, MD, SCHA, and perturbation theory.

Close modal

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.

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 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,57,71,90,96,97 which will be the focus of future work.

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

The authors have no conflicts to disclose.

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

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

1.
M.
Born
and
K.
Huang
,
Dynamical Theory of Crystal Lattices
(
Oxford University Press
,
1954
).
2.
M. T.
Dove
,
Introduction to Lattice Dynamics
(
Cambridge University Press
,
1993
).
4.
M. L.
Klein
and
G. K.
Horton
,
J. Low Temp. Phys.
9
,
151
(
1972
).
5.
F.
Knoop
,
T. A. R.
Purcell
,
M.
Scheffler
, and
C.
Carbogno
,
Phys. Rev. Mater.
4
,
083809
(
2020
).
6.
G.
Grimvall
,
B.
Magyari-Köpe
,
V.
Ozoliņš
, and
K. A.
Persson
,
Rev. Mod. Phys.
84
,
945
(
2012
).
7.
P. B.
Allen
,
Mod. Phys. Lett. B
34
,
2050025
(
2019
).
8.
A.
Glensk
,
B.
Grabowski
,
T.
Hickel
, and
J.
Neugebauer
,
Phys. Rev. Lett.
114
,
195901
(
2015
).
9.
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
).
10.
11.
A. A.
Maradudin
and
A. E.
Fein
,
Phys. Rev.
128
,
2589
(
1962
).
13.
14.
D. A.
Broido
,
A.
Ward
, and
N.
Mingo
,
Phys. Rev. B
72
,
014308
(
2005
).
15.
A.
Togo
,
L.
Chaput
, and
I.
Tanaka
,
Phys. Rev. B
91
,
094306
(
2015
).
16.
W.
Li
,
J.
Carrete
,
N.
A Katcho
, and
N.
Mingo
,
Comput. Phys. Commun.
185
,
1747
(
2014
).
17.
N. K.
Ravichandran
and
D.
Broido
,
Phys. Rev. X
10
,
021063
(
2020
).
18.
T.
Sun
,
X.
Shen
, and
P. B.
Allen
,
Phys. Rev. B
82
,
224304
(
2010
).
19.
T.
Sun
and
P. B.
Allen
,
Phys. Rev. B
82
,
224305
(
2010
).
20.
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
(
2020
).
21.
X.
Yang
,
T.
Feng
,
J.
Li
, and
X.
Ruan
,
Phys. Rev. B
105
,
115205
(
2022
).
22.
23.
24.
N. K.
Ravichandran
and
D.
Broido
,
Phys. Rev. B
98
,
085205
(
2018
).
25.
K.
Esfarjani
and
Y.
Liang
,
Nanoscale Energy Transport
(
IOP Publishing
,
2020
), pp.
7-1
7-35
.
26.
P.
Souvatzis
,
O.
Eriksson
,
M.
Katsnelson
, and
S.
Rudin
,
Comput. Mater. Sci.
44
,
888
(
2009
).
27.
T.
Tadano
and
S.
Tsuneyuki
,
Phys. Rev. B
92
,
054301
(
2015
).
28.
T.
Tadano
and
S.
Tsuneyuki
,
J. Phys. Soc. Jpn.
87
,
041015
(
2018
).
29.
R.
Bianco
,
I.
Errea
,
L.
Paulatto
,
M.
Calandra
, and
F.
Mauri
,
Phys. Rev. B
96
,
014111
(
2017
).
30.
L.
Monacelli
,
R.
Bianco
,
M.
Cherubini
,
M.
Calandra
,
I.
Errea
, and
F.
Mauri
,
J. Phys.: Condens. Matter
33
,
363001
(
2021
).
31.
A.
van Roekeghem
,
J.
Carrete
, and
N.
Mingo
,
Comput. Phys. Commun.
263
,
107945
(
2021
).
32.
P.
Choquard
,
The Anharmonic Crystal
(
WA Benjamin
,
1967
).
33.
V.
Samathiyakanit
and
H. R.
Glyde
,
J. Phys. C: Solid State Phys.
6
,
1166
(
1973
).
34.
R. M.
Levy
,
A. R.
Srinivasan
,
W. K.
Olson
, and
J. A.
McCammon
,
Biopolymers
23
,
1099
(
1984
).
35.
L. T.
Kong
,
G.
Bartels
,
C.
Campañá
,
C.
Denniston
, and
M. H.
Müser
,
Comput. Phys. Commun.
180
,
1004
(
2009
).
36.
L. T.
Kong
,
Comput. Phys. Commun.
182
,
2201
(
2011
).
37.
O.
Hellman
,
I. A.
Abrikosov
, and
S. I.
Simak
,
Phys. Rev. B
84
,
180301(R)
(
2011
).
38.
O.
Hellman
,
P.
Steneteg
,
I. A.
Abrikosov
, and
S. I.
Simak
,
Phys. Rev. B
87
,
104111
(
2013
).
39.
O.
Hellman
and
I. A.
Abrikosov
,
Phys. Rev. B
88
,
144301
(
2013
).
40.
J.
Klarbring
,
O.
Hellman
,
I. A.
Abrikosov
, and
S. I.
Simak
,
Phys. Rev. Lett.
125
,
045701
(
2020
).
41.
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
).
42.
A. H.
Romero
,
E. K. U.
Gross
,
M. J.
Verstraete
, and
O.
Hellman
,
Phys. Rev. B
91
,
214310
(
2015
).
43.
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. Mater.
5
,
035004
(
2021
).
45.
46.
S. A.
Egorov
,
K. F.
Everitt
, and
J. L.
Skinner
,
J. Phys. Chem. A
103
,
9494
(
1999
).
47.
R.
Ramírez
,
T.
López-Ciudad
,
P.
Kumar P
, and
D.
Marx
,
J. Chem. Phys.
121
,
3973
(
2004
).
48.
H. Y.
Geng
,
J. Phys. Chem. C
126
,
19355
(
2022
).
49.
T.
Morresi
,
L.
Paulatto
,
R.
Vuilleumier
, and
M.
Casula
,
J. Chem. Phys.
154
,
224108
(
2021
).
50.
T.
Morresi
,
R.
Vuilleumier
, and
M.
Casula
,
Phys. Rev. B
106
,
054109
(
2022
).
51.
H.
Mori
,
Prog. Theor. Phys.
33
,
423
(
1965
).
54.
C.
Hijón
,
P.
Español
,
E.
Vanden-Eijnden
, and
R.
Delgado-Buscalioni
,
Faraday Discuss.
144
,
301
(
2010
).
55.
A.
Carof
,
R.
Vuilleumier
, and
B.
Rotenberg
,
J. Chem. Phys.
140
,
124103
(
2014
).
56.
Z.
Li
,
X.
Bian
,
X.
Li
, and
G. E.
Karniadakis
,
J. Chem. Phys.
143
,
243128
(
2015
).
57.
A.
Fiorentino
and
S.
Baroni
, “
From Green-Kubo to the full Boltzmann kinetic approach to heat transport in crystals and glasses
,”
Phys. Rev. B
107
,
054311
(
2022
).
58.
J.
Bosse
,
W.
Götze
, and
M.
Lücke
,
Phys. Rev. A
17
,
434
(
1978
).
59.
D. R.
Reichman
and
P.
Charbonneau
,
J. Stat. Mech.: Theory Exp.
2005
,
P05013
.
60.
T. E.
Markland
,
J. A.
Morrone
,
K.
Miyazaki
,
B. J.
Berne
,
D. R.
Reichman
, and
E.
Rabani
,
J. Chem. Phys.
136
,
074511
(
2012
).
61.
62.
R.
Kubo
,
M.
Toda
, and
N.
Hashitsume
,
Statistical Physics II
(
Springer
,
Berlin, Heidelberg
,
1991
).
63.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
121
,
3368
(
2004
).
65.

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

66.
B. J.
Braams
and
D. E.
Manolopoulos
,
J. Chem. Phys.
125
,
124105
(
2006
).
67.
M.
Rossi
,
M.
Ceriotti
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
140
,
234116
(
2014
).
68.
R.
Zwanzig
,
Nonequilibrium Statistical Mechanics
(
Oxford University Press
,
2001
).
69.
B. J.
Berne
and
G. D.
Harp
,
Advances in Chemical Physics
(
John Wiley and Sons, Inc.
,
1970
), pp.
63
227
.
70.
G.
Szamel
,
J. Chem. Phys.
127
,
084515
(
2007
).
71.
M.
Simoncelli
,
N.
Marzari
, and
F.
Mauri
,
Nat. Phys.
15
,
809
(
2019
).
72.
F.
Eriksson
,
E.
Fransson
, and
P.
Erhart
,
Adv. Theory Simul.
2
,
1800184
(
2019
).
73.
F.
Zhou
,
W.
Nielson
,
Y.
Xia
, and
V.
Ozoliņš
,
Phys. Rev. B
100
,
184308
(
2019
).
74.
F.
Bottin
,
J.
Bieder
, and
J.
Bouchet
,
Comput. Phys. Commun.
254
,
107301
(
2020
).
75.
R.
Ramírez
and
T.
López-Ciudad
,
J. Chem. Phys.
111
,
3339
(
1999
).
76.
N.
Shulumba
,
O.
Hellman
, and
A. J.
Minnich
,
Phys. Rev. B
95
,
014302
(
2017
).
77.
Y.
Xia
,
K.
Pal
,
J.
He
,
V.
Ozoliņš
, and
C.
Wolverton
,
Phys. Rev. Lett.
124
,
065901
(
2020
).
78.
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
,
Research
2020
,
4589786
.
79.
Z.
Zeng
,
C.
Zhang
,
Y.
Xia
,
Z.
Fan
,
C.
Wolverton
, and
Y.
Chen
,
Phys. Rev. B
103
,
224307
(
2021
).
80.
X.
Yang
,
J.
Tiwari
, and
T.
Feng
,
Mater. Today Phys.
24
,
100689
(
2022
).
81.
T.
Feng
,
L.
Lindsay
, and
X.
Ruan
,
Phys. Rev. B
96
,
161201(R)
(
2017
).
82.
A.
Decoster
,
J. Phys. A: Math. Gen.
37
,
9051
(
2004
).
83.
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
(
2014
).
84.
D.
Dangić
,
O.
Hellman
,
S.
Fahy
, and
I.
Savić
,
npj Comput. Mater.
7
,
57
(
2021
).
85.
J.
Bouchet
and
F.
Bottin
,
Phys. Rev. B
95
,
054113
(
2017
).
86.
V.
Ladygin
,
P.
Korotaev
,
A.
Yanilkin
, and
A.
Shapeev
,
Comput. Mater. Sci.
172
,
109333
(
2020
).
87.
J.
Ding
,
T.
Lanigan-Atkins
,
M.
Calderón-Cueva
,
A.
Banerjee
,
D. L.
Abernathy
,
A.
Said
,
A.
Zevalkink
, and
O.
Delaire
,
Sci. Adv.
7
,
eabg1449
(
2021
).
88.
L.
Xie
,
J. H.
Feng
,
R.
Li
, and
J. Q.
He
,
Phys. Rev. Lett.
125
,
245901
(
2020
).
89.
X.
He
,
D.
Bansal
,
B.
Winn
,
S.
Chi
,
L.
Boatner
, and
O.
Delaire
,
Phys. Rev. Lett.
124
,
145901
(
2020
).
90.
L.
Isaeva
,
G.
Barbalinardo
,
D.
Donadio
, and
S.
Baroni
,
Nat. Commun.
10
,
3853
(
2019
).
91.
J.
Eckert
,
W.
Thomlinson
, and
G.
Shirane
,
Phys. Rev. B
16
,
1057
(
1977
).
92.
C.
Cazorla
and
J.
Boronat
,
Rev. Mod. Phys.
89
,
035003
(
2017
).
93.
V.
Goldman
,
G.
Horton
, and
M.
Klein
,
Phys. Rev. Lett.
21
,
1527
(
1968
).
94.
R. A.
Aziz
,
V. P. S.
Nain
,
J. S.
Carley
,
W. L.
Taylor
, and
G. T.
McConville
,
J. Chem. Phys.
70
,
4330
(
1979
).
95.
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
).
96.
M.
Simoncelli
,
N.
Marzari
, and
F.
Mauri
,
Phys. Rev. X
12
,
041011
(
2022
).
97.
G.
Caldarelli
,
M.
Simoncelli
,
N.
Marzari
,
F.
Mauri
, and
L.
Benfatto
,
Phys. Rev. B
106
,
024312
(
2022
).