Equilibrium and Non-Equilibrium Lattice Dynamics of Anharmonic Systems

In this review, motivated by the recent interest in high-temperature materials, we review our recent progress in theories of lattice dynamics in and out of equilibrium. To investigate thermodynamic properties of anharmonic crystals, the self-consistent phonon theory was developed, mainly in the 1960s, for rare gas atoms and quantum crystals. We have extended this theory to investigate the properties of the equilibrium state of a crystal, including its unit cell shape and size, atomic positions and lattice dynamical properties. Using the equation-of-motion method combined with the fluctuation–dissipation theorem and the Donsker–Furutsu–Novikov (DFN) theorem, this approach was also extended to investigate the non-equilibrium case where there is heat flow across a junction or an interface. The formalism is a classical one and therefore valid at high temperatures.


Introduction
There has recently been intense research on developing high-temperature materials that can operate at temperatures of around 1500 K or above. These materials can potentially be used in jet engines, petrochemical and material processing, and power plants based on concentrated solar power. Another application would be in hypersonic vehicles, for which the surface temperature can exceed 2000 • C. These materials are typically borides or carbides made of refractory elements along with their high-entropy alloy combinations [1,2]. At such high temperatures, the atomic vibrational amplitude, being proportional to (kT/mω 2 ) 1/2 , grows beyond the region where harmonic approximation is valid, and therefore, more sophisticated models are required to describe such dynamics. A successful yet relatively simple model that can handle anharmonic systems is the self-consistent phonon (SCP) theory originally proposed by Born and Hooton [3]. This theory was then extended and further developed by Boccara and Sarma [4] using a variational approach based on free-energy minimization in order to study displacive phase transitions. Later, Choquard used the cluster expansion theory [5] and Ranninger used many-body diagrammatic theory [6] to derive similar equations. The first complete SCP theory was, however, initially worked out as a mean-field theory by Gillis et al. in 1968 [7]. It was then completed to second-order by Werthamer in 1970 [8].
In the first section of this review, we go over this theory, its variational formulation and its modern implementations based on first-principles Density Functional Theory (DFT) calculations. The second part of this review deals with the non-equilibrium case where an anharmonic device is connected to multiple probes or thermostats with different temperatures. We then derive, based on the equation of motion of atoms in the device, an expression for the thermal current and entropy generation in the device.

The Self-Consistent Phonon Theory
Following the variational formulation first proposed by Boccara and Sarma [4], we ask ourselves: For a system at a given temperature, what is the "best" harmonic Hamiltonian that describes the thermodynamics of this system (the "best" being defined as the one that minimizes a variational free energy)? According to Bogoliubov's inequality [9], if the Hamiltonian is a sum of a solvable Hamiltonian H 0 = T + V 0 and a remaining potential energy term V − V 0 , then its exact free energy satisfies the following inequality: The average in this inequality is with respect to the density matrix associated with the solvable Hamiltonian H 0 [9]: A 0 = Trρ 0 A; ρ 0 = e −H 0 /kT /Z; Z = Tr e −H 0 /kT ; F 0 = −kTlnZ (2) In the solvable Hamiltonian H 0 , which, in our case, is a harmonic Hamiltonian, the force constants between pairs of atoms are free parameters that can be determined from the above minimization once the potential energy V is known and its average can be calculated.
In what follows, we assume that the potential energy V is a polynomial in powers of atomic displacements about their equilibrium positions, and we truncate this expansion at the fourth-order: where i, j, k, l refer to atoms, and u i , p i , and m i refer to the displacement, momentum, and mass, respectively, of atom i. The Greek letters φ, ψ, χ are, respectively, the second, third and fourth derivatives of the potential energy evaluated at the equilibrium positions. They are also called the force constants, and are harmonic in the case of φ and anharmonic for ψ and χ.
We have omitted the Cartesian indices for brevity. In case there are multiple atoms in a unit cell, the label i can more generally be extended to include the labels τ of the atoms in the cell labeled by R, as well as the Cartesian component of u. In other words, i = (Rτα). In order to include structural phase transitions and thermal expansion, we introduce two additional sets of variational parameters: the strain tensor η and u 0 i , which is the internal relaxation of atoms beyond the strain effect (captured by ηR i ) due to changes in the unit cell shape or volume. Thus, the general displacement of atom i from its reference position as a result of cell deformation can be written as: u i (t) = ηR i + u 0 i + (1 + η)y i (t) = s i + (1 + η)y i (t), where y i (t) is the dynamical displacement of atom i about its equilibrium, so that by definition, y i 0 = 0. The vector R i refers to the lattice translation vector of the cell containing i. Further, s i = u i = ηR i + u 0 i is the static displacement of atom i due to strain and internal relaxations, and is introduced for brevity of notations. The coefficients of this expansion can be obtained from a zero-temperature DFT calculation, for instance, as in [10,11].
The trial harmonic potential containing the variational parameters K ij is defined as: For the trial harmonic Hamiltonian, the free energy can be calculated analytically, as it is the sum of free energies associate with each harmonic mode: where β = 1/k B T and ω kλ are, respectively, the harmonic frequency of mode λ and the wave-vector k obtained from diagonalizing the dynamical matrix associated with the trial force constants: where k refers to a k vector in the first Brillouin zone, and the indices τ refer to atoms in the primitive cell. Note that the existence of the term τ − τ in the exponent does not change the eigenvalues and only causes a phase shift in the eigenfunctions. The matrix D, being Hermitian, has real eigenvalues denoted by ω 2 kλ and eigenvectors e τα,λ (k), where α is the Cartesian coordinate and λ refers to the vibrational mode (τα can be understood as the line index and λ as the column index of the unitary eigenvector matrix e for each k). We also need the thermal averages of the trial and anharmonic potentials, V 0 and V, respectively. In terms of the eigenvalues and eigenvectors of the above dynamical matrix and the phonon creation and annihilation operators a and a † , respectively, the displacements y can be written as [12] so that for the thermal average of displacement autocorrelations, we have where a † kλ a k λ = n kλ δ λ,λ δ k,k ; a k λ a † kλ = (n kλ + 1)δ λ,λ δ k,k were used. Using the identity 2n kλ + 1 = coth βhω kλ 2 , the equilibrium average of the harmonic potential can be written as: We see, therefore, that the trial free energy and harmonic potential are both only functions of the harmonic frequencies ω kλ , and thus K ij , while the average of the anharmonic potential V depends also on η and u 0 i . As for the model anharmonic potential expressed in Equation (3), since the variable y has a Gaussian distribution, the average of its odd powers is zero, and thus, the thermal average of the actual potential V can be presented in Equation (9) as: where for the sake of simplicity of notation, all the Cartesian indices have been omitted, and we used C ij = u i u j = s i s j + (1 + η) 2 y i y j . We also used the invariance of the force constants under permutation of indices. Further, note that in C, each term (1 + η) is always contracted with one y i .

Determination of the State of Thermal Equilibrium at a Given T
Having the equilibrium average of the anharmonic potential, we can proceed to minimize the variational or trial free energy F trial = F 0 (K) + V(C, s) − V 0 (K) in an iterative manner by setting its gradients with respect to different variational parameters to zero and solving the resulting set of non-linear equations in (u 0 , η, K). If we think of K ij and C ij = u i u j = s i s j + (1 + η) 2 y i y j as independent variables in order to avoid the use of the chain rule to express this dependence, minimization of the trial free energy F trial (u 0 , η, K, C) leads to an additional, fourth non-linear equation. Thus, we have the following gradient formulas resulting from the minimization condition: The first two equations express the equilibrium conditions: no force on atoms and no stress in the unit cell. In the third Equation (12), since the anharmonic part V does not depend on K, we recover a general thermal property of harmonic Hamiltonians that leads, using the chain rule, to the expression of the displacement correlations in terms of the eigenvalues and eigenvectors of the dynamical matrix, as we have expressed in Equation (7). Finally, the fourth equation defines what the effective or optimal force constants K need to be in terms of the parameters of the potential energy V. Boccara and Sarma [4] showed that effective force constants are given by the thermal average of the second derivative of the potential energy: K ij = ∂ 2 V ∂u i ∂u j 0 , which gives them a clear physical interpretation: The effective harmonic force constants are the thermal average of the second derivative of the potential, i.e., instead of evaluating the second derivative at u i = 0, it needs to be averaged over all values of the dynamical variables u i (t) weighted by their Boltzmann distribution probability.
Equation (13) provides a relationship between C ij and K ij , allowing elimination of K as a variational parameter. From Equation (9), we have: so that Equation (13) reduces to the thermal-averaged second derivative: Therefore, the set of non-linear equations to solve for (u 0 , η, K, C) are Equations (7), (10), (11) and (15). Note, this is a set of self-consistent equations, as the knowledge of K ij depends, through Equation (15), on C ij or y i y j , which, in turn, depends on the eigenvalues and eigenvectors of K through Equation (7).
If one wants to impose an external pressure or a general stress tensor and find the equilibrium state under applied external stress σ ext , the following thermodynamic potential F(η) − d 3 r σ ext αβ η αβ needs to be minimized with respect to η. Since we have considered C and s to be independent, we have: For convenience, we also reproduce the derivatives of the quantities s and C with respect to u 0 , η, in full notation: (17) ∂s Rτµ In conclusion, for the adopted model potential in Equation (3), the four non-linear equations to be solved self-consistently are Equations (7), (15), (21) and (22). The latter two are reproduced below, including all Cartesian indices for the sake of clarity and completeness: The self-consistent procedure may be implemented as follows: (0) Start from initial guesses for (u 0 , η, K ij ); (1) Diagonalize the dynamical matrix obtained from K, and calculate y i y j from Equation (7); (2) Use y i y j and η in Equation (21) to solve for u 0 contained in s and C; (3) Use Equation (22) to solve for η from the knowledge of u 0 , K and yy ; (4) Use (s, C) to calculate the new K from Equation (15); (5) Eventually mix the result with old Ks to avoid numerical instabilities; (6) Go to 1) until the process converges. Note the non-linear Equations above, (21) and (22), are only quadratic in u 0 and η. This procedure is not always straightforward, and care needs to be taken when solving numerically: for instance, the dynamical matrix obtained, K, in Step 1 may not be positive definite. Another subtlety is in solving Equation (16), which is a non-linear equation in S and may have multiple solutions. In this review, we do not discuss numerical issues that will be elaborated on in a future publication. Before ending this section, one last note is on the actual implementation of the algorithm to achieve self-consistency. The algorithm outlined above is described just for the conceptual clarity of the approach, but there are more efficient ways to reach self-consistency. Namely, we are using Broyden's method [13,14] to solve this non-linear set of equations. The equations to solve are basically the free-energy gradients with respect to (u 0 , η, K, C) being equal to zero. Broyden's method precisely solves a non-linear system of the form G( X) = 0, with the variables (u 0 , η, K, C) being stored in an array X, and the array G containing the free energy gradients with respect to X. After every iteration of Broyden, where all (u 0 , η, K, C) are updated at the same time, the displacement autocorrelation matrix y i y j is calculated from the eigenvalues and eigenvectors of K ij using Equation (7) and is used to calculate the gradients G, which are fed back into the Broyden algorithm for the next iteration. Non-positivity of the dynamical matrix can be dealt with by taking atomic displacements along the softest modes and shear deformations along − ∂F trial ∂η αβ to lower the trial free energy until a minimum is reached along that direction. The flowchart of this approach is displayed in Figure 1. Using this scheme, not only is it possible to obtain the phonon spectrum at finite temperatures, but one may also obtain the equilibrium shape of the unit cell and atomic positions defined by η, u 0 . Thus, this method can predict solid-solid phase transitions as a function of temperature (and also imposed pressure). Two simple implementations of this method have been included in our previous paper [15,16]. One is a single particle in an anharmonic quartic potential in which the change of the equilibrium position and the vibrational frequency versus temperature is illustrated. The other is the treatment of the phase change in a one-dimensional chain with up to sixth-order onsite potential energy and up to quartic nearest-neighbor pair interactions. At low temperatures, the equilibrium position is away from the lattice site, while above the critical temperature, it has a discontinuous jump to the lattice site.

Previous Works/Implementations by Other Groups
We need to mention other similar works in which self-consistent phonon approximation has been implemented from a set of first-principles DFT calculations on forces as a function of displacements, and we point out the small variations compared to the present work. To the best of our knowledge, these include, in chronological order:

•
The SCAILD (Self Consistent Ab Initio Lattice Dynamics) method by Souvatzis et al. in 2008 [17,18], in which the effective dynamical matrix was extracted from the Fourier transform of the forces and displacements in the reciprocal space. • TDEP (Temperature-Dependent Effective Potential) method by Hellman et al. in 2011 [19][20][21], in which both harmonic and cubic force constants are extracted from a molecular dynamics simulation performed at several temperatures in a supercell. While they all enforce self-consistency between the atomic displacements and the effective harmonic force constants, they differ in the way the phase space is sampled and force constants are calculated. The numerical challenges we outlined, especially near a phase transition, still exist in all these methods. In addition, many of these methods, with the exception of SSCHA, do not implement thermal expansion, meaning the effective temperature-dependent strain tensor η and residual atomic relaxations u 0 i are often not calculated self-consistently unless supercell DFT energy minimization is performed between updates.

Summary
To conclude, we have developed a formalism based on the self-consistent phonon theory in which thermodynamic properties of a crystalline material are computed from the knowledge of the higher-order anharmonic force constants (up to quartic has been illustrated in this paper). It is therefore possible to predict the phonon dispersion, thermal expansion and equilibrium shape of the unit cell and the atoms therein (including phase change) as a function of temperature. The novelty of the approach is that all thermodynamic calculations are done from knowledge of only the force constants of the high-symmetry phase. The limits of the validity of this approach will be tested in a future publication. It is worth mentioning that the implementation as outlined in this paper is much faster than DFT calculations in a supercell, as all averages are Gaussian and can be calculated analytically. The disadvantage is the accuracy of the Taylor expansion force field, which is truncated to fourth-order and up to some nearest-neighbors interactions. In this case, either the order of expansion has to be increased to sixth or higher, or new force constants φ, ψ, χ need to be calculated for the new phase in order to calculate its thermodynamic properties (and not its phase change).

Non-Equilibrium Heat Transport in Anharmonic Systems
The previous section dealt with equilibrium properties at finite temperatures. In this section, however, we are interested in heat current flow under an applied temperature difference across a device. Thus, this will be a non-equilibrium situation, but the formalism will have some similarities to the equilibrium case in terms of thermal expansion and effective force constants. Thermal averages need to be computed, but they will be nonequilibrium thermal averages, for which we use a different approach based on the equationof-motion method, the fluctuation-dissipation theorem and the Donsker-Furutsu-Novikov theorem [32][33][34], which we describe in detail in the following.
Transport theories of non-interacting quantum systems based on the Keldysh formalism, which treats non-equilibrium flow of charge or heat, have been developed in the past for a one-dimensional (1D) geometry. In this review, we propose a simplified version in the classical or high-temperature limit using a lighter formalism based on the equation-of-motion method.
To the best of our knowledge, the first quantum development of an electron transport theory applied to nanostructures across a junction was done in 1971 by Caroli et al. [35] where a Green's function formalism was used to describe transmission of electrons in a 1D system. Following the seminal work of Caroli et al, many other groups worked on similar formalisms and proved a formula for the transmission through the system, which is now widely used for both non-interacting electrons and phonons. The equilibrium version of it, namely T = Tr[GΓ L G † Γ R ] , where T, G and Γ are, respectively, the transmission, the retarded Green's function and the escape rates to the leads, was established by Meir and Wingreen [36], and in a similar form, by Pastawski [37] in 1991. This formula holds for a non-interacting (in the case of electrons) or harmonic (in the case of phonons) system near equilibrium, meaning the chemical potential or temperature gradients are also infinitesimally small. These assumptions might not always be realistic, especially in small (mesoscopic) systems subject to temperature differences over fractions of a micrometer, and a formulation for non-equilibrium situations and interacting systems is preferable for the sake of testing the domain of the validity of the equilibrium formulas and for more accurate description in the case of large interactions and large driving fields.
In this review, we show a semi-classical atomistic derivation of the transport theory of heat carriers in a solid across a junction connected to thermal reservoirs based on Green's function formalism. The basic geometry of our problem is multi-probe where the system in which scatterings occur is connected to multiple reservoirs, contacts or heat sinks that impose their temperature and cause flow of heat carriers (see Figure 2). This model is used for mesoscopic systems where the carrier mean free path could be on the order of the system length, implying that Ohm's law of addition of resistances in series does not necessarily hold, and coherence can play an important role. It is also worth noting that since we are in the non-equilibrium regime in general, temperature or chemical potential are not necessarily well-defined concepts, and thus will be avoided. The geometry of the reservoirs is fundamentally one-dimensional (1D), and if there is translational symmetry perpendicular to the current flow, one can use Bloch's theorem to decouple the 3D system into many non-interacting 1D systems, each labeled by a quantum number, which is the transverse momentum. So in what follows, we assume such decoupling has been done, and that we will be dealing with strictly 1D semi-infinite leads, although the central device is arbitrary in shape and structure and may be connected to multiple 1D probes. In this paper, we are interested in transport of anharmonic phonons, or, more generally, vibrational modes, in mesoscopic systems under large temperature differences. For simplicity, we use a classical description. A generalization to the quantum case will be inferred at the end. So in our classical treatment, the frequency ω is just the frequency of vibrational modes and not the energy of phonons. This classical formalism avoids fancier mathematics involving commutation relations and concepts such as time-ordered or contour-ordered Green's functions. It will only involve "retarded" or causal Green's functions, which helps us solve a differential equation in the frequency domain. Typically considered geometries will be identical to a non-equilibrium molecular dynamics (NEMD) setup, where the two ends of the system are attached to two thermostats at different temperatures, and one is interested in measuring the interfacial thermal conductance (see Figure 2). The ND atom anharmonic structure is connected to two semi-infinite reservoirs in this example (called L and R). After elimination of the reservoirs' degrees of freedom, the infinite system is replaced by an isolated "cluster" subject to thermostating forces η L and η R applied on the boundary, reflecting the equation of motion in Equation (33). As a result of this boundary condition, force constants φ, ψ, χ are replaced with effective force constants Φ, Ψ, X and a harmonic self energy σ.
The non-equilibrium anharmonic phonon problem has been addressed in the past by Mingo [38] and separately by Wang [39] in 2006. They used the many-body perturbation approach of non-equilibrium quantum systems based on the Keldysh formalism, (also called the Non-Equilibrium Green's Function or NEGF method) and derived a lowest-order approximation for the transmission function. Dai and Tian [40] recently implemented this rigorous formulation to calculate the effect of cubic anharmonicity on the phonon transmission function through an ideal interface by applying it to Si/Ge and Al/Al with two different masses. A similar model that explicitly incorporates transverse momentum dependence was also developed recently by Guo et. al. [41] based on previous work by Luisier [42]. Polanco has a recent review of these methods based on NEGF [43]. These NEGF-based models, although fully quantum mechanical, do not include anharmonicity beyond cubic order nor any thermal expansion effects. In this work, we intend to go beyond cubic and include thermal expansion effects self-consistently.
Other calculations of transmission in the non-equilibrium regime [44] based on the Green's function method have been based on self-consistent reservoirs (also called the Buttiker probe method,) which was first proposed by Bolsterli et al. in 1970 [45]. In this method, to obtain the local temperature, an atom is connected with a very weak coupling to a fictitious probe at a given temperature with which it enters equilibrium. The coupling to the probe is weak so that it does not disturb the heat flow within the system. The probe temperature is changed until the net heat current from the system to the added fictitious probe is zero. Therefore, because there is no flow, the probe is in equilibrium with that atom, and the probe temperature defines the temperature of the atom. Note that the system could be out of equilibrium, so the temperature is not really well-defined on a given atom. This shortcoming is also present in NEMD simulations where the assigned "temperature" is just the average kinetic energy of the atom although there is no evidence of local thermal equilibrium. So even though non-equilibrium effects are included through this Buttiker probe method, it does not include anharmonicity. We should mention at this point that there is numerical evidence of absence of equipartition in the vibrational modes near the interface [46,47], implying that a definition of local temperature is not really justified near an interface where the drop could be large and therefore ∇T is not a small perturbation.
In another work, to include anharmonic corrections in the transmission, obtained from MD trajectories, Saaskilahti et al. [48] extended the harmonic formulation of the transmission function by Chalopin et al. [49,50]. In the actual calculation, arguing that the anharmonic part of the current is usually small, they only used its harmonic formula but with velocities and positions coming from the full anharmonic atomic trajectories. The advantage of this approach over NEMD is that the heat current and thermal conductance can be decomposed in the frequency domain. Approaches based on MD trajectories, while including full anharmonicity, suffer from noise and would require a large number of simulations in order to perform proper ensemble averaging, whereas many-body approaches might be inaccurate if a perturbative expansion in powers of anharmonicity is used; otherwise, they do not suffer from noise and treat ensemble averaging analytically. In this review, we show how some of these limitations can be overcome by adopting a non-perturbative self-consistent many-body approach by fully including in the current the effect of anharmonic terms introduced in the Hamiltonian. Furthermore, using a classical method to derive an expression for the heat current, we argue that to leading order, it is necessary to include quartic terms in the Hamiltonian to ensure stability and properly describe both the thermal expansion and the dominant temperature dependence effect on the heat current.

Equations of Motion and Langevin Thermostats
We start by defining our model and the assumptions. A multiprobe geometry is assumed. Figure 2 shows a two-probe example, in which a central region, also called the "device", is connected to two semi-infinite one-dimensional (1D) leads that play the role of thermostats imposing a temperature at the boundaries of the system. For a 3D geometry where periodic boundary conditions are used in the transverse direction, i.e., perpendicular to the heat flow direction, as shown in the figure, a Fourier transformation will reduce it to m 2 one-dimensional problems, where m is the number of transverse kpoints. We are not concerned with temperature drops in the thermostats, which are assumed to be harmonic heat sinks and follow Langevin dynamics.
The anharmonic Hamiltonian of the device (D) and its coupling to the lead or probe α are: The dynamical variable u i (t) refers to the displacement of atom i about the zero-temperature equilibrium position (the force on the atom is zero for u = 0). The leads labeled by α are semi-infinite chains following Langevin dynamics. After the standard change of variables to x i (t) = √ m i u i (t), and with V α,il = W α,il / √ m i m lα and Φ α,ll = φ α,ll /m lα , we arrive at the following equation of motion for atoms in the central region: Note we use capitalized Greek letters (Φ, Ψ, ...) for mass-rescaled force constants and lower-case Greek letters (φ, ψ, ...) for the bare potential energy derivatives. The letter α refers to the leads, and the dynamical variable x = (x 1 , ..., x N ) can be thought of as an array containing the displacements of all the atoms in the central region (also called the device), Φ as the force constant matrix between such atoms, and V α as the force constant matrix connecting atoms of the lead α to atoms in the device. Finally, a = −1/2Ψx 2 + ... is the anharmonic part of the force, which, for now, we keep as a for brevity. The dynamics of atoms in the lead are of Langevin type, where a set of identical coupled harmonic oscillators are subject to damping γ α and noise ζ α as follows: where the superscript T stands for transpose. The force constants Φ α can be thought of as effective FCs at the temperature of interest, so that we do not need to introduce anharmonicity in the leads, which merely play the role of absorbing phonons from the device and reinjecting thermalized phonons into the device. To insure each thermostat has temperature T α , the damping coefficient γ and the random forces ζ must satisfy the following fluctuation-dissipation relation in Equation (27): We proceed by eliminating the lead variables x α in Equation (25) using the Green's function method. To this end, we start by taking the Fourier transform of the above two equations according to: The equations of motion satisfied by the Fourier transform of the displacements become: Note that the frequency-domain variables are represented with capitalized letters. Now that the differential equations are transformed to algebraic ones, one can easily proceed to eliminate the lead degrees of freedom in the main equation of motion by using the Green's functions. Let g α (ω) be the retarded (causal) Green's function associated with the lead α. The positivity of the damping factor γ α insures causality. The solution to Equation (30) after the transients have decayed to zero can be written as: where In Equation (29), we need −V α X α = η α + σ α X, which we obtain from Equation (31). Here, η α (ω) = −V α g α (ω)ζ α (ω), and σ α = V α g α V T α . Likewise, defining the retarded Green's function of the central region as G −1 = [−ω 2 + Φ − ∑ α σ α (ω)], we can write the solution to the central region as: The function σ α (ω) = V α g α (ω)V T α is traditionally called the self-energy of lead α, and it seems like an added correction to the harmonic FC's Φ. Its effect is to shift the vibrational frequencies of the device and give them a a finite lifetime, as these modes can leak into the leads. Omitting the transient contribution of the initial conditions, this is the solution to the equations of motion, which depend on the stochastic functions η α = −V α g α ζ α .
To summarize, as is also illustrated in Figure 2, the thermostats have been eliminated and replaced by Langevin forces η α applied directly on the device atoms connected to a thermostat, and the harmonic FCs have been renormalized to become Φ − ∑ α σ α .

Entropy Generation Rate
An important quantity of interest is the entropy generation rate in the device, which in steady state can be expressed in terms of heat fluxes j α as:Ṡ = − ∑ α j α /T α . This quantity results from the interaction between the device degrees of freedom and the random noise due to Langevin thermostats. It involves the thermostat's temperature and the incoming currents, which will be discussed next.

Heat Current
The heat flow from leads to the device is microscopically defined as the average rate lead atoms do work on the device atoms. The work rate is the product of the velocity degrees of freedom of the device multiplied by the force from the lead acting on them. With this link being harmonic, the heat current expression is simply: where the trace is taken over the device degrees of freedom, and averaging is over the thermostat degrees of freedom. Since the current depends on the stochastic functions ζ, we need to take the ensemble average over all realizations of the random forces subject to the fluctuation-dissipation theorem [51]. To get the steady-state current, we also take the time average. In terms of its Fourier transform, the time-integrated current is J α (Ω = 0) = ∞ −∞ j α (t) dt. The current may be calculated using Equations (31) and (33), and the result can be written as: The DC response is found by taking the Ω → 0 limit and dividing by the timeintegration parameter τ → ∞. As shown later, the results for j α do not depend on τ. Note that because we are interested in the DC response, only diagonal terms of correlations (Ω = 0) are needed here. Thus, calculation of the heat current is reduced to calculation of the two correlation functions and the so-called lead self-energy σ α (ω), followed by a frequency integration. Note that this current from lead α is the sum of two terms: the first one, proportional to Z α = X(ω) η † α (ω) , is the work of the stochastic forces on the device; while the second one, proportional to a displacement autocorrelation X(ω) X † (ω) , is the work of the lead dampers trying to reabsorb some of the excess energy injected from the device in order to reestablish thermal equilibrium in the lead.
The average denoted by is an average of the stochastic forces in the thermostats, which have white noise characteristics. With the leads being at different temperatures, the average becomes a non-equilibrium average and can only be calculated using the equation of motion (33) and the statistical properties of the Langevin thermostats following the fluctuation-dissipation theorem. For anharmonic interactions involving higher powers of displacements in A, the calculation of current leads to a hierarchy of equations, each containing higher powers of displacements, and has, so far, been computed using different approximations [38,39,52].
One should note that at equal lead temperatures, one is in equilibrium and no net current enters the device: j α = 0. If the lead temperatures are different and one is in the non-equilibrium regime, in the absence of heat generation in the device, we should still have Σ α j α = 0. This is also known as Kirchhoff's law of current conservation in the context of electrical circuits.

Thermal Expansion and Renormalization of the Force Constants
At finite temperatures, equilibrium or not, the average position of the atoms is shifted from its zero-temperature reference value because of anharmonicity. Treatment of the thermal expansion is similar to that of the equilibrium case, where solving Equation (16) for given y i y j leads to the shifts s, with the exception that the thermal averages are now non-equilibrium averages. So first, we proceed to calculate the thermal expansion and the resulting renormalization of the force constants.
The set of equations expressing that the force on each atom is, on average (over time and over random Langevin forces), zero ∂H D ∂x i = 0, leads to the average atomic positions x i = s i satisfying the following cubic equation in the parameters s: Due to these position shifts, as well as thermal fluctuations of atomic positions, the force constants effectively change, similar to the SCP case. One can show that the effective harmonic FCs will still be given by Equation (15), again, with the exception that the thermal average is a non-equilibrium one. In this case, we change displacement variables from x to y such that y(t) = x(t) − x = x(t) − s, which has zero average by construction. The resulting equations of motion satisfied by y can be shown to be: where the renormalized anharmonic force is and the renormalized cubic force constant isΨ = Ψ + χs. (The quartic term is affected if there are fifth or higher derivatives of the potential. In this approximation where we limit ourselves to the quartic term in the expansion of the potential energy, the latter is not renormalized.) The renormalized Green's function corresponding to this equation of motion is therefore defined as leading to the general solution: One could call this the non-equilibrium mean-field, or the non-equilibrium self-consistent phonon (NESCP) approximation. In our original paper on this topic [53], we explain that this is the "best" mean-field approximation, as it includes the mean-field corrections − ∂A ∂Y = 1 2 χ YY in the harmonic force constants and leads to leading-order cancellation of anharmonic terms in the expansion of the correlation functions Yη † and YY † in powers of ∂A ∂Y and A, respectively.

Correlation Functions Yη † and YY †
Now, we proceed to the calculation of the two needed correlation functions Yη † and YY † in order to obtain the heat current j α , given as below.
is the power exerted by the random force η on the device and can therefore be interpreted as the heat injected per unit time and unit frequency (mode) from lead α into the device. We can calculate this term using the solution to the equation of motion, Equation (38).
In the NESCP approximation, A is neglected, and we have Z NESCP α = G η α η † α /τ. The noise autocorrelation is obtained from the fluctuation-dissipation theorem, Equation (27), and it leads to (see Appendix B.2, Equation (A8)): where, for brevity, we replaced the "occupation factor" k B T α /ω with f α , (With a factor of h in the denominator f α would be a real Bose-Einstein occupation factor taken to the classical limit) and τ represents the integration time, which goes to infinity and is canceled in the final expressions. For white noise, one can show that the result does not depend on the thermostat damping parameter γ, and thus, we adopt this type of noise for the thermostats. Finally, within NESCP, the displacement noise correlation function is given by:

Beyond NESCP
To go one step further beyond SCP and include the effect of anharmonicity in Z α , we use the Donsker-Furutsu-Novikov (DFN) identity [54] (for a proof, see Appendix C), which states that for any functional of the white noise f [η], we have It has the advantage of lowering the powers of η in f . Using this theorem, we have: From the equation of motion for Y and the chain rule, we find so that finally, This is an exact relation satisfied by Z α . Note that, by construction, the first term involving G ∂A/∂Y G is identically zero because ΨY = 0 and χYY /2 is already included in K, which appears in the denominator of G. More details about the proof can be found in [53]. In the language of many-body theory, this is the expanded form of Dyson's equation involving the thermal average of powers of the anharmonic force derivatives. This result, even though exact, cannot be calculated because the average of the inverse of powers of Y cannot be calculated exactly. Instead, in the spirit of perturbation theory, one may add the power-series term-by-term provided the sum is convergent. Another non-perturbative approach, which we adopt here, is to sum a series of terms up to infinite order.
and ∂A ∂Y = ∂A ∂Y − ∂A ∂Y = −ΨY − 1 2 χ(YY − YY ), which has a zero average by construction. We can note that the major part of the quartic anharmonicity has been removed, and the expansion is in powers of ∂A/∂Y, which is centered at zero and therefore has the smallest higher moments. When raised to second power in the above formula for Z α , cubic and quartic terms become decoupled if we neglect averages of terms of odd power in Y, which are expected to be small if non-zero. At low temperatures or weak anharmonicity, the dominant contribution to the second-order terms comes from cubic terms and can be written as:

Displacement Autocorrelations YY †
Finally, the last correlation function needed is C(ω) = Y(ω)Y † (ω) /τ. Note the autocorrelation matrix C is Hermitian. As before, within the SCP, where the effective force constant K appears in the denominator of G, one arrives at a simple expression for C: Let us first give a physical interpretation of this function as it is related to the heat current. As shown in Equation (A6) in the Appendix B.2 , if (Σ α Γ α ) = Γ is twice the imaginary part of the denominator of the Green's function G, then GΓG † = 2 (G), so that at equilibrium, where all leads are at the same temperature, C NESCP = 2 G f eq . On the other hand, G is related to the vibrational density of states (DOS) through: Therefore, if all lead temperatures are equal, the function DOS(ω) f eq (ω)/h = −ωTr C NESCP (ω)/πh can be interpreted as the equilibrium number of vibrational excitations of frequency ω present in the device due to its contact with the leads (Since f α = kT α /ω, then f α /h = kT α /hω is the classical (high-temperature) limit of the equilibrium Bose-Einstein occupation factor, and multiplied by the DOS, it can be interpreted as the number of vibrational excitations of frequency ω.). On the other hand, as the imaginary part of the device's Green's function Γ α /2ω can be understood as the rate of absorption of these vibrational excitations by the leads, we can say that −TrC NESCP Γ α /2πh is the rate of flow of phonons of frequency ω going from the device into the lead α. The heat current can be obtained by multiplying this particle current by the energyhω of each carrier, leading to ωTrC NESCP Γ α /2π being the heat flowing from device to lead α. This is indeed the second term in the expression of the heat current in Equation (34), the first term being the heat flow from the lead into the device.

Beyond NESCP
Now, we can go one step beyond NESCP and, using the solution to the equation of motion (Equation (38)), express the exact formula for the displacement autocorrelation as: The first term is the NESCP contribution discussed above. The second and third terms involve terms linear in A that also appeared in the calculation of Z α , and finally, the last term involves the second power of the anharmonic force, which includes powers of displacement equal to or higher than four (see Equation (45)), and as such, cannot be exactly calculated unless the non-equilibrium distribution function of displacements Y is known.
As for the displacement autocorrelations, they have the terms Aη † and AA † sandwiched between G and G † . The former can be shown to exactly satisfy: The last term, to the lowest-order in cubic anharmonicity, is: AA † =ΨΨ YYYY /4. Not having the distribution function of Y, this expression cannot be calculated; but again, to the leading order in anharmonicity, we use a decoupling approximation to write it as: This approximation becomes exact in the limit where all leads have the same temperature, so it should be reliable for small temperature differences. It is also good in the limit where the four displacement operators are not the same. Finally, collecting all quadratic terms inΨ, we have: The symbolic diagrams summarizing this approximation are shown in Figure 3. 3.6. The Particular Case of Non-Equilibrium Self-Consistent Phonon Approximation: NESCP To summarize, within the NESCP, the renormalized anharmonic force A is neglected, and we end up with simpler expressions for the correlation functions, similar to the harmonic case: The frequency in the denominator cancels the frequency in the current coming from the time derivative of positions, so the NESCP current becomes: It can be shown that this current satisfies both detailed balance and current conservation [53].
As a result, this expression can be simplified to the familiar form: This equation is linear in temperature differences and clearly satisfies Kirkhoff's law of current conservation, as it is antisymmetric with respect to swapping of lead indices.
In general, Tr [Γ α G Γ β G † ] may be interpreted as the transmission from lead α to lead β. While in a harmonic theory, the transmission is temperature-independent, within NESCP, it actually depends on the leads' temperatures through G, which contains K, which itself contains the temperature-dependent term YY or C (see Equation (15)).
To obtain the heat current within the NESCP model, one needs to self-consistently solve Equations (15), (35), (37) and (47) defining, respectively, C NESCP , thermal expansion s, effective force constant K and the effective Green's function G. Once converged, these quantities can be used in Equation (54) to calculate the current.

Conclusions
To summarize, in this review, we presented a formalism to describe thermodynamic and thermal transport properties of solids at high temperatures within a lattice dynamics approach. The first part of this review of our recent work concerned the thermal equilibrium properties of crystalline solids. Starting from calculations performed in a supercell, one can extract the force constants [10] and thus define the anharmonic model Hamiltonian. The variational formulation then allowed definition of the thermal equilibrium state: the free energy needs to be minimized with respect to the parameters defining the equilibrium state, namely the primitive cell translation vectors, atomic positions and effective harmonic force constants. Thermodynamic properties, including the phonon frequencies, can then be obtained from the latter.
In the second part, we developed a self-consistent current-conserving approximation for anharmonic systems out of equilibrium in the high-temperature (classical) regime. Although the set of derived Equations for the current (Equation (34)), the Equation of motion (Equation (38)), and for defining the correlation functions (Equations (43) and (48)) were formally exact, one has to develop approximations to solve Dyson's Equation (Equation (43)) and the Equation defining displacement autocorrelations C (Equation (48)). One reason is because the anharmonic force A is an infinite Taylor expansion and is usually truncated, and another reason is because its derivative appears in the denominator of Equation (43), which cannot be exactly inverted. In this work, we truncated the Taylor expansion of A up to quartic terms and only included up to second powers of A and ∂A/∂Y in Equations (43) and (48), leading to the approximation summarized in Equation (51).
Furthermore, we showed that thermal expansion needs to be included using both cubic and quartic terms (to avoid any divergence), and it has the effect of renormalizing FCs as T is increased. The reference GF to work with, G, has two corrections: one due to thermal expansion implying changes in bond length and strength, and the other due to thermal fluctuations about the average position ( 1 2 χ YY , similar in spirit to the equilibrium self-consistent phonon theory Equation (15), except that one is not at thermal equilibrium, and averages depend on the temperature of all attached probes). Non-equilibrium averages were possible to calculate with the use of the fluctuationdissipation theorem (Equation (40)), the equations of motion (Equation (38)) and the DFN theorem (Equation (42)). Within the NESCP, where A is dropped from the equation of motion, we end up with a simple description that includes only the average of quartic anharmonicity in the effective force constants. We argue that besides the thermal expansion, this is the leading anharmonic correction to the thermal current.
An alternative approach to investigate non-equilibrium effects at and near interfaces would be to perform a non-equilibrium molecular dynamics simulation (NEMD) of the system attached to thermostats at different temperatures and sample the atomic trajectories in the phase space to find the distribution functions, position averages and force constant averages by fitting the forces to a linear model F NEMD i ≈ F Harmonic i = −K ij y j based on the knowledge of the forces on atoms and their positions in each MD snapshot in order to extract the effective (non-equilibrium) harmonic force constants K. The remainder can then be defined as the anharmonic force F NEMD i = −K ij y j + a i (y), and the present results maybe used. This approach, however, has inherent noise in it. If error in the averages is small, one still can use Equation (54) to compute current if there is no access to the heat current expression in terms of interactions. Another possibility is to use the approximate harmonic formula of the heat current in order to evaluate it [48].
Despite the approximations used in this work, the advantage of this formalism over MD simulations that includes anharmonicity to all orders is that it is analytical and therefore fast and free of simulation noise, although reaching self-consistency can be challenging for some model systems. It would be desirable to compare the results with those of NEMD to validate these approximations for a given system. The accuracy also relies on the force field and strength of higher-order terms: sources of divergence would be in the denominator of Equation (43)  Acknowledgments: K.E. would like to thank M. R. Rahimi Tabar for useful discussions and for introducing him to the DFN theorem, and J. Shiomi, H. Cheraghchi, N. Ravichandran and T. Hamakawa for discussions at the early stages of this work. He also specifically thanks V. Chiloyan for introducing him to the Langevin thermostat method.

Conflicts of Interest:
The authors declare no conflict of interest. The sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Calculation of Time Averages
When calculating time average of a product such as a(t)b(t) in terms of their Fourier transform, some care needs to be taken: where one set of brackets is for time averaging and the second set is a thermodynamic average over different initial conditions. To simplify the notations, we have however used only one set of brackets.
In terms of their Fourier transform, we have On the other hand, taking τ → ∞ in the boundaries of integral of e −iωt leads to 2πδ(ω) which should cancel the τ in the denominator in order to give 1. This means, when taking time averages, one simply needs to take the convolution of the Fourier transforms at Ω = 0: When calculating diagonal terms of autocorrelations, such as in ζ(ω)ζ † (ω) = 2γk B Tτ, the result is proportional to the integration time τ. The latter cancels the τ in the denominator coming from time averaging. Loosely speaking, 2πδ(ω)/τ δ ω,0 .
For Langevin thermostats with white noise, these results do not depend on the thermostat damping factors γ.

Appendix B.2. Noise Autocorrelation Functions
Let us start with the noise autocorrelation which will appear in the calculation of displacement-noise correlation Z α . This autocorrelation can be obtained from the fluctuation-dissipation theorem [51], a relation which has to hold if the noise is such that the leads are to behave as a Langevin thermostat at temperature T: The noise is white and different sites are uncorrelated with each other. This relation implies that in the frequency domain the autocorrelation of ζ and that of η satisfy the following relations: In the above, g α is a diagonal matrix of size equal to the number of degrees of freedom in the lead α (which is infinity!) At this point, we will use a convenient identity satisfied by any Green's function G. If G −1 = a + ib, with (a, b) real matrices, then Since (g −1 α ) = −ωγ α we can write the η-autocorrelation as: where we used the notation Γ α = −i(σ α − σ † α ) = 2 (σ α ) for twice the imaginary part of the lead α self-energy. Alternatively, the diagonal elements in frequency can be written as: where τ is the integration time which goes to infinity. Dhar and Roy [55] have shown that in the quantum limit, when taking a semi-infinite harmonic lead and averaging over all possible initial conditions sampled from a canonical ensemble, the quantum noise term has an autocorrelation given by where f (x) = [e x − 1] −1 is the equilibrium Bose-Einstein distribution function, which, in the classical (high-temperature) limit, reduces to k B T/hω. Our classical result based on properties of Langevin thermostats is thus consistent with the quantum one. We can note that the explicit dependence on the damping factor γ α has been replaced by twice the imaginary part of the lead self-energy Γ α , which one may interpret as 2ω times a "rate". For the adopted white noise, the dependence on its damping factor has gone away! Furthermore, in the calculation of Z α , the notation ηη † implies the diagonal terms ω = −ω must be taken. Using the Donsker-Furutsu-Novikov (DFN) identity (see Appendix C) we can see that Aη † ∝ ηη † and therefore only diagonal terms in frequency will appear in the frequency integral of the heat current, and assuming thermostat temperatures are steady, the heat current can only have a DC (Ω = 0) component. Note that this component is infinite the way we defined it: J α (Ω = 0) = ∞ −∞ j α (t) dt. To find the average heat current, we have to divide this by the integration time τ, and then take the limit τ → ∞. This division cancels the τ appearing in the numerator of noise autocorrelations. This issue is further discussed in the appendix, Equation (39). Final results do not depend on τ.
The factor 2πδ(ω + ω )/τ can now be excluded from the current as it is essentially 1 for the diagonal terms, and zero otherwise.