Theory of Non-Equilibrium Heat Transport in Anharmonic Multiprobe Systems at High Temperatures

We consider the problem of heat transport by vibrational modes between Langevin thermostats connected by a central device. The latter is anharmonic and can be subject to large temperature difference and thus be out of equilibrium. We develop a classical formalism based on the equation of motion method, the fluctuation–dissipation theorem and the Novikov theorem to describe heat flow in a multi-terminal geometry. We show that it is imperative to include a quartic term in the potential energy to insure stability and to properly describe thermal expansion. The latter also contributes to leading order in the thermal resistance, while the usually adopted cubic term appears in the second order. This formalism paves the way for accurate modeling of thermal transport across interfaces in highly non-equilibrium situations beyond perturbation theory.


Introduction
Transport theories of non-interacting quantum systems based on the Keldysh formalism, which treats non-equilibrium flow of charge or heat carriers in a one-dimensional (1D) geometry have been developed in the past. To the best of our knowledge, the first such development was performed in 1971 by Caroli et al. [1] where a Green's function formalism was used to describe dynamics of electrons in a 1D crystal. 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, now widely used for both non-interacting electrons and phonons. The equilibrium version of it, namely T = Tr[GΓ L G † Γ R ] , where G and Γ are, respectively, the retarded Green's function and the escape rates to the leads, was established by Meir and Wingreen [2] and in a similar form by Pastawski [3] in 1991. This formula holds for a non-interacting (or harmonic, in the case of phonons) system near equilibrium, meaning the chemical potential or temperature gradients are to be 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 validity of the equilibrium formulas and more accurate description in the case of large interactions and large driving fields.
The general and practical problem of interest is the description of the heat transport across interfaces between different materials; however, to keep it more general, we adopt a multi-probe geometry where the system in which scatterings occur is connected to multiple energy reservoirs, which impose their temperature, and cause the flow of heat carriers (see Figure 1). This model is used for mesoscopic systems where the carrier mean free path could be on the order of the system length, implying Ohm's law of addition of resistances in series does not necessarily hold, and coherence can play an important role. 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 performed and we will be dealing with strictly 1D semi-infinite harmonic leads, although the central device is arbitrary in shape and structure and maybe connected to multiple 1D probes, assumed to be Langevin thermostats and assumed to be isothermal.
For simplicity, we will 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 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 help us solve a differential equation in the frequency domain. In the hightemperature limit both theories lead in principle to the same results. To model interfacial thermal resistance, typical 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 established heat flow as a result of the temperature difference (see Figure 1). Figure 1. (Top) a general "molecular" multiprobe geometry where the device defined by atoms within the ellipse are connected to 4 reservoirs imposing a temperature or measuring a current. (Bottom) the 2-probe geometry. Atoms in each lead are connected to a harmonic thermostat at fixed temperature T L and T R . In the Buttiker probe, also called the self-consistent reservoirs geometry, to measure the "local" temperature, each layer of the central device may also be weakly connected to a fictitious thermostat at a temperature to be (self-consistently) determined so that the net heat flow from that probe to the device is zero.
The non-equilibrium anharmonic phonon problem has been addressed in the past by Mingo [4] and separately by Wang [5] in 2006. They used the many-body perturbation approach of non-equilibrium quantum systems based on the Keldysh formalism, (also called the Non-Equilbrium Green's Function or NEGF method) and derived a lowestorder approximation for the transmission function. Dai and Tian [6] recently implemented this rigorous formulation to calculate the effect of cubic anharmonicity on the phonon transmission function through an ideal interface, applied to Si/Ge and Al/Al with two different masses. A similar model, which explicitly incorporates transverse momentum dependence, was also developed recently by Guo et al. [7], based on previous work by Luisier [8]. Polanco has a recent review of these methods based on NEGF [9]. These NEGFbased models, although fully quantum mechanical, do not include anharmonicity beyond cubic order nor any thermal expansion effects. Other calculations of the transmission in the non-equilibrium regime [10] based on the Green's function method, have been based on the self-consistent reservoirs (also called Buttiker probe method), which was first proposed by Bolsterli et al. in 1970 [11]. In this method, as shown in Figure 1, every layer is connected, with a very weak coupling, to a fictitious probe at a given temperature with which it becomes into equilibrium. The probe temperature, which is also the assigned temperature of that layer, is obtained from the constraint that the net heat current from the system to the added fictitious probe should be zero. Note that the system could be out of equilibrium, so that a temperature is not really well-defined at a given layer. This shortcoming is also present in NEMD simulations where the assigned "temperature" of a layer is just the average kinetic energy of atoms in that layer, 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 [12,13], implying that a definition of local temperature is not really justified near an interface. In another work, based on the MD simulation, which fully includes anharmonicity, Saaskilahti et al. [14] extended the harmonic formulation of the transmission function based on the actual MD trajectories by Chalopin et al. [15,16], to include anharmonic corrections. 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, in order to deduce the interfacial thermal conductance. 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, but otherwise do not suffer from noise and treat ensemble averaging analytically. In this work, we try to overcome these limitations by adopting a non-perturbative many-body approach by fully including in the current the effect of anharmonic terms introduced in the Hamiltonian, without using Buttiker probes. 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, in order to properly describe both the thermal expansion and the dominant temperature dependence effect in the heat current.

Dynamics
We start by defining our model and the assumptions. A multiprobe geometry is assumed as shown in Figure 1 in which a central region, also called the "device" is connected to many semi-infinite one-dimensional (1D) leads, playing the role of thermostats imposing a temperature at the boundaries of the system. We will not be concerned with temperature drops in the thermostats, which are assumed to be harmonic and follow Langevin dynamics. If needed, parts of the leads can be incorporated in the device region to illustrate the temperature drops near the interface. The device (D) Hamiltonian is a fourth-degree polynomial in powers of atomic displacements, and it is harmonically coupled to the lead α is represented by Equation (2): The dynamical variable u i (t) refers to the displacement of atom i about the zero-temperature equilibrium position (the force on atom is zero for u = 0). The leads labeled by α are semiinfinite harmonic chains. After the standard change of variables to x i (t) = √ m i u i (t), and to 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 have used 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 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 is of Langevin type [17], where a set of identical coupled harmonic oscillators are subject to damping γ α and noise ζ α . The equation of motion for atoms with displacements x α is as follows: where the superscript T stands for transpose. Here again all atomic degrees of freedom are stored in the array x α = (x 1α , x 2α , . . .). The force constants Phi α 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. We will proceed by eliminating the lead variables x α in Equation (3) using the Green's function method. To this end, we start by taking the Fourier transform of the above two equations according to: In the frequency domain, the equations of motion (3) and (4) become: Note that the frequency-domain variables are represented with capitalized letters, and the transient part of the solution will be omitted as we are interested in the steady-state solution. 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 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 (7) after the transients have decayed to zero, can be written as: where In Equation (6) we need −V α X α , which we can obtain from Equation (8): 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 shows the effect of this lead on the spectrum of the device, which is given by the poles of G. Its real part provides a correction to the eigenvalue spectrum ω 2 , and its imaginary part, divided by 2ω, gives the inverse lifetime of an excitation of the central region caused by interactions with the lead. It is the rate at which the excitation leaks into the lead. 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 α ζ α .

Physical Observables
The equations of motion we derived are deterministic for every realization of the random forces. To simulate the real thermodynamical behavior of the baths, so that a temperature can be assigned to them, one needs to perform an "ensemble" average, denoted by . . . over all realizations of the forces subject to the constraint imposed by the fluctuation-dissipation (FD) theorem [17]: The noise is white and different sites i, j of leads α, α are uncorrelated with each other. Physical observables are then obtained after an ensemble average is performed over forces. This is where irreversibility is introduced in this deterministic formalism, as a result of which, entropy is generated in the device. The latter can be understood as the log of the distribution function of the device as the random forces are varied within the constraints imposed by the FD theorem.

Entropy Generation Rate
An important quantity of interest is the entropy generation rate in the device, which can be expressed 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 thermostats temperature and the incoming currents which will be discussed next.

Heat Current
The main quantity of interest is the heat current. The heat from lead α can be defined as the net rate at which energy is flowing into the device from that lead. It is the work performed from lead α on the device's degrees of freedom per unit time, which is the product of the velocity degrees of freedom of the device times the force from lead acting on them: where the trace is taken over the device degrees of freedom. Substituting the expressions for x and x α from Equations (10) and (11) we can obtain the final expression of the current in the frequency domain. Since the current depends on the stochastic functions η, we will take its ensemble average to find the response of the system to an applied temperature difference. We will start by taking the Fourier transform of j α (t) and call it J α (Ω). Note due to time integration, we must have J α (Ω = 0) = τ j α , where τ → ∞ is the time integration window. Finally the static DC current is given by: where we used the notation Γ α = −i(σ α − σ † α ) = 2 (σ α ) for twice the imaginary part of the lead α self-energy. The DC response is found by taking the Ω → 0 limit. Note that because we are interested in the DC response, only diagonal terms of correlations (Ω = 0) are needed here. Thus the calculation of the heat current is reduced to the 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 per unit time of the stochastic forces on the device, while the second one X(ω) X † (ω) Γ † α , proportional to a displacement autocorrelation, is the work of the lead dampers trying to reabsorb some of the excess energy injected from the device in order to re-establish thermal equilibrium in the lead.
The average denoted by is an average over the stochastic forces in the thermostats that have white noise characteristics. When it is performed over the device degrees of freedom, the leads being at different temperatures, it becomes a non-equilibrium average and can only be calculated using the equation of motion (11) and the statistical properties of the Langevin thermostats, i.e., the fluctuation-dissipation theorem. For anharmonic interactions involving higher powers of displacements in A, the calculation of current will lead to a hierarchy of equations, each containing higher powers of displacements, and has so far been computed using different approximations [4,5,18].

Constraints
Before proceeding to the calculation of the correlation functions, we recall two constraints that the heat current needs to satisfy. The first is the detailed balance relation, which states that if all leads are at the same temperature T, the net current j α should be identically zero for all leads α. The second constraint is that of current conservation, which in steady state Ω → 0, and under no additional heat generation in the device, reduces to Σ α j α = 0. This is also known as the Kirchhoff's law in the context of electrical circuits. Note that AC components of the current need not satisfy this constraint as they reflect the information on transient currents during the relaxation process, and depend on the heat capacity of the system. Any physically correct description of transport should exactly satisfy these two constraints.

Thermal Expansion
As the temperature of a system is raised, there can be thermal expansion due anharmonicity. The equilibrium position of the atoms is shifted, and this will also cause a change in the force constants as bond lengths have changed. To take these effects into account, while simplifying the notations, we will slightly modify the formalism as follows: the displacement variable X is changed to Y(ω) = X(ω) − X or y(t) = x(t) − x which has zero average by construction. The resulting nonlinear equations satisfied by X are derived by taking the average of the equation of motion (11) or equivalently setting the average force on each atom to zero (see Appendix B for more details).
The resulting equations will depend on correlations such as Y i Y j , Y i Y j Y k and higher powers. Accordingly, the potential energy derivatives will be evaluated at the zero of Y, and will be denoted with a bar sign on top of them (Φ →Φ etc. . . ). While the variable X satisfies the equation of motion: the new variable Y satisfies: Next we will linearize the forces with respect to y: where we have set ( ∂V ∂y ) y=0 = 0 to define the thermal expansion x = x 0 (see Appendix B). As we will show, the effect of temperature will be to renormalize the FCs, not only through thermal expansion but also due the thermal fluctuations as we will show using the nonequilibrium mean-field approximation (NEMF) also detailed in Section 6. Next, we define a renormalized Green's function using renormalized force constantsΦ = ∂ 2 V ∂y 2 y=0 as With this GF, the displacements Y satisfy Note one can add any constant λ to the force constantΦ in the above GF, provided λY is also added to the anharmonic force A. We will make use of this freedom in the next section to further simplify the formalism.

Force Constant Renormalization
Given the form of the above equations, we can add −Y ∂A ∂Y to A and add − ∂A ∂Y to 1/G so that now the renormalized GF becomes: while the renormalized anharmonic force now becomes A = A − Y ∂A ∂Y in the right hand side of Equation (16). This renormalization of harmonic force constants captures a major part of anharmonicity (because ∂A/∂Y = 0), and is in spirit very similar to the lowest-order self-consistent phonon theory, which in the past has been applied to equilibrium systems. The advantage of this renormalization is that, as we will see, the lowest anharmonic correction in Z α disappears by construction since ∂A/∂Y = 0, leading to the smallest variance and higher moments of ∂A/∂Y.
We will refer to this choice of the reference GF as the Non-equilibrium mean-field approximation (NEMF). With this choice, the equation of motion for Y becomes: The Feynman diagram associated with the new Green's function G is shown in Figure 2. The thick line with opposite arrows represents the displacement autocorrelation C(ω 1 ). Internal frequency ω 1 is integrated over.

Displacement-Noise Correlations
One can see from Equation (13) that the calculation of the heat current requires the calculation of the displacement-noise correlation Z α and displacement autocorrelations C. We proceed to the calculation of these quantities first within the harmonic approximation, and then in the presence of anharmonic forces of the form a = −ψ y 2 /2 − χ y 3 /6.
Let us start with the noise autocorrelation, which will appear in the calculation of displacement-noise correlation Z α . In the frequency domain, using the fluctuationdissipation theorem Equation (12), one can derive (see Appendix C.2): where, for brevity, we have replaced the "occupation factor" k B T α /ω by f α , and τ represents the integration time which goes to infinity and cancels the τ in the expression for the current j α = J α /τ. With a Factor ofh in the denominator of f α , the latter would be a real Bose-Einstein occupation factor taken to the classical limit.
In the case of a white noise, we show in Appendix C.2 that the result will not depend on the thermostat damping parameter γ, and thus we adopt this type of noise for the thermostats.
Next, we need to calculate displacement-noise correlations: Z α = Yη † α . Since −iωY is a velocity, ω (Z α ) is the power exerted by the random force η on the device and therefore can be interpreted as the heat injected per unit time and unit frequency (mode) from lead α into the device.
In the harmonic case (A = A = 0; G = G = G), using the equation motion (18), this expression is simplified to: For non-zero anharmonicity, the three GFs are different, and adopting G as the reference GF, we have the NEMF approximation to the displacement-noise correlation function as: To go one step further and include the effect of anharmonicity in Z α , we will use the Novikov-Furutsu-Donsker identity [19] (for a proof also see Appendix D), 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: (18) and the chain rule, we find

From the equation of motion Equation
so that finally, Note that, by construction, the first term involving ∂A/∂Y is identically zero. This justifies the choice of G rather than G as the reference Green's function.
In the language of many-body theory, this is the expanded form of Dyson's equation involving thermal average of powers of the anharmonic force derivatives. This exact result cannot be calculated, because the average of the inverse of powers of Y cannot be calculated exactly, and instead one may add the power series term by term provided the sum is convergent. This is one of the main results of this paper, providing a more accurate and likely convergent expression for the displacement-noise correlations, and leading to a simple calculation of heat currents if cubic terms are to be neglected (using the NEMF reference Green's function G).
The next-order correction consists in adding the effect of anharmonicity to second order, i.e., writing the displacement-noise correlation function as: . 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 has therefore the smallest higher moments. When raised to the 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 can be written as: An explicit form of this equation is provided in the following section and in Appendix C.3 Equation (A17).

Displacement Autocorrelations
Finally, the last correlation function needed is C(ω) = YY † . Note the autocorrelation matrix C is Hermitian. The function ωC(ω) can be interpreted as the (non-equilibrium) number of excitations of frequency ω present in the device due to its contact with the leads. If all lead temperatures are equal, we recover the equilibrium occupation times the total DOS, G(Σ α Γ α )G † = 2 (G), which is the total (equilibrium) number of excitations in the device. In the heat current, this term appears as −ωCΓ α /2 and can be interpreted as the heat current going from (because of the negative sign) the device into this lead as Γ α /2ω is the escape rate into the lead α.
Similar to the treatment of Z α , we will start with the equation of motion, Equation (18), and the FD theorem, Equation (19), to write C as: The first term is the NEMF contribution and is reduced to a known result, with now the renormalized Green's function G being used instead of the standard harmonic one (G or G), in order to include thermal effects to some extent: The NEMF approximation, which consists in using G and neglecting the contributions of anharmonicity included in A, is very similar to the harmonic approximation. Within this approximation, where C ≈ C NEMF and Z α ≈ Z NEMF α , the transmission becomes Tr [GΓ L G † Γ R ], very similar to the harmonic result (see Appendix E), but with the GF substituted by the renormalized G, and includes the effect of (quartic) anharmonicity to the lowest order.
Linear terms in A lead to terms similar to Z α , which has already been discussed. The only remaining difficulty is with the AA † terms, which do not explicitly contain any noise term, but have higher powers of displacements. Such terms can only be calculated approximately as the use of the equations of motion will involve higher powers of displacements.
To have a second-order approximation consistent with that used for Z α in Equation (26), we have to use: It can be shown that this approximation, taken to a self-consistent level satisfies current conservation, meaning ∑ α j α = 0, however the spectral components of the current: Z α (ω) − C(ω)Γ α (ω)/2 do not necessarily lead to zero when summed over all leads. Including for completeness both the cubic and quartic components of the anharmonic forces to second order, the self-consistent set of equations to be solved with the reference Green's function G are: The equations for Z α and C can also be represented using Feynman diagrams as shown in Figures A1 and A2 in Appendices C.3 and C.4. More explicit forms of these equations are also reproduced in this appendix as Equations (A17) and (A18).
Starting inputs for C and Z α could be their NEMF values in the right-hand sides of the above equations, and the latter can be solved iteratively until convergent. Note that while the term Z α requires ∂A/∂Y, the terms C require A itself, but these equations contain only second powers of A and ∂A/∂Y. Once iterations converge, the obtained C and Z α functions can then be inserted in Equation (13) to compute the heat currents from each lead.
These equations would be the same as the ones obtained from the many-body non-equilibrium Keldysh formalism with the difference that the "occupation factors" f α = k B T/ω are classical ones, instead of Bose-Einstein functions. In this sense, they can directly be compared to results from classical non-equilibrium MD simulations, which are exact in anharmonicity but have inherent statistical noise in them.
Another interesting feature to note is the increase in the overall conductance with the temperature (if ∆T is held small). This is in agreement with previous MD simulations [13,14].

Conclusions
To summarize, we developed a self-consistent approximation for transport in anharmonic systems out of equilibrium in the high-temperature (classical) regime. There is therefore no factors ofh in the formalism and ω is to be interpreted as frequency only, not energy. Although the set of derived equations for the current Equation (13) the equation of motion Equation (18), and the Equations (24) and (27) defining the correlation functions, were formally exact, one has to develop approximations to solve the Dyson's Equation (24) and the Equation (27) defining C. One, because the anharmonic force A is an infinite Taylor expansion and is usually truncated, and two, because its derivative appears in the denominator of Equation (24), 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 (24) and (27).
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 (Φ instead of Φ), and the other due to thermal fluctuations about the average position ( ∂A/∂Y , which usually involves the quartic term and the autocorrelation C), similar in spirit to the self-consistent phonon theory, except that one is not at thermal equilibrium. This is the leading-order anharmonic correction, and cubic anharmonicity contributes to second-order correction terms as shown in the self-consistent equations (29). Non-equilibrium averages were possible to calculate with the use of the fluctuationdissipation theorem (Equation (19)), the equations of motion (Equation (18)) and the NFD theorem (Equation (22)).
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 and the position averages-this, however, has inherent noise in it.
One way to extract the effective force constants is to fit from the knowledge of the forces on atoms and their positions in each MD snapshot, the forces to a linear model F NEMD i ≈ F Harmonic i = −Φ ij y j in order to extract the effective (non-equilibrium) harmonic force constantsΦ. The remainder can then be defined as the anharmonic force: F NEMD i = −Φ ij y j + a i (y), and the present results may be used. Despite the approximations used in this work, the advantage of this formalism over MD simulations, which 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 make a comparison of the results with 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 (24) if G ∂A/∂Y ≈ 1, signaling resonances, in which case the Taylor expansion in Equation (24) is not appropriate.
Applications to nanoscale systems will appear in future publications.
Funding: Internal support at UVa from the Hobby Fund is greatly acknowledged.

Informed Consent Statement: Not applicable.
Acknowledgments: I would like to thank M. R. RahimiTabar for useful discussions on the handling of correlations, the NFD theorem, and a review of the manuscript, and J. Shiomi and H. Cheraghchi for discussions at the early stages of this work. I specifically thank V. Chiloyan for introducing me to the Langevin thermostat method. This paper is dedicated to the memory of Rouzbeh Rastgarkafshgarkolaei with whom I had several related discussions.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

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 However, 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πδ(ω)/τ = 1 or 0. As an example, to calculate the average (DC) current one simply needs to take the diagonal terms in frequencies, omitting factors of τ appearing in the correlation functions:

Appendix B. Change of Position Variables Due to Thermal Expansion
As the temperature of a system is raised, there can the thermal expansion due to anharmonicity. The equilibrium position of the atoms is shifted, and this can also cause a change in the force constants as bond lengths have changed. To take these effects into account, we will slightly modify the formalism as follows: Let us Taylor-expand the interaction potential in the device in powers of rescaled atomic displacements x = √ m u, about the zero temperature equilibrium positions: where, for brevity, we omitted the summation sign over repeated indices, and where the residual force on atom i is zero because one starts from a fully relaxed configuration. The coefficients of this expansion can be obtained from a zero-temperature DFT calculation for instance [20,21].
One way, and actually the correct way to compute the non-equilibrium average of the force constants, is to perform a non-equilibrium molecular dynamics (NEMD) simulation in which the system is subject to two or more thermostats at different temperatures and zero pressure. Average positions x maybe computed from the runs, and then one may fit the forces with an effective harmonic model as we will describe below. In a thermally expanded system where x is non-zero, we Taylor expand the potential about the "non-equilibrium" values x instead of zero. This leads to a change in the force constants. The new dynamical variables are chosen so that their average is zero, i.e., they oscillate about the new thermal non-equilibrium positions: The average position is defined to be the solution of the average force being zero. In terms of these new variables, and up to the second order in powers of displacements, the force on atom i can be written as: The first term, the residual force, will be used to define the thermal expansion x (by setting the residual or average force to zero: , which is the Equation (14) in the main text. The second term is the effective, or temperature-dependent, harmonic force, and the last term is the effective anharmonic force. In the case where the potential energy is Taylor-expanded up to quartic order such as in Equation (A3), one can explicitly work out the equation satisfied by x , and derive explicit expressions for the effective harmonic and anharmonic FCs. The expression for the force in this case is as follows: where use was made of the symmetry of Ψ and χ under permutation of their j, k, l indices.

Appendix B.1. Thermal Expansion
To find the average positions x , we set the average force on each atom i, F i = − ∂V /∂x i x= x to zero, and x α = 0 (leads are harmonic and will not thermal-expand), and solve for x . Thus, the coupling term with the leads, which is linear in lead degrees of freedom and has zero average, will vanish, and we have the following set of non-linear equations in x in terms for force constants and averages of yy and yyy : χ ijkl y j y k y l (A5)

Appendix B.2. The New Force Constants
The new force constants are simply defined to the be potential derivatives evaluated at the new average positions. Up to quartic order, they are defined as: Finally, the average positions can be re-expressed in terms of the new force constants as: These coupled set of non-linear equations (for all i) in ( x ,Φ,Ψ,χ must be solved in terms of yy and yyy covariance matrices, to find the average positions and effective force constants. It requires a self-consistent iterative solution as the yy correlations will in turn depend on x .

Appendix B.3. New Equations of Motion
Finally, changing to new dynamical variables y(t) and new force constants, the equation of motion for y becomes: Adopting this change of variable, the GF keeps the same form with Φ substituted byΦ, and anharmonic forces in the right-hand side of the equation of motion become by definition: 6χ ijkl y j y k y l (A8)

Appendix C. Explicit Form of the Correlation Functions including the Atomic and Cartesian Indices
In this section, we will give explicit formulas for the correlation functions needed in the average heat current expression in Equation (A2). First let us label by roman letters i, j, k the dynamical degrees of freedom in the device. If there are N atoms in the device in three dimensions, each of these labels refers to an atom and one of the three Cartesian components of its displacements. It therefore varies from 1 to 3N. As a result the matrices C, σ α and Z α are 3N × 3N matrices. Taking their trace in Formula (A2) gives the heat current, which is a scalar. Note that since transport is along the length of the leads, which are one-dimensional, in principle the components of the power perpendicular to the lead direction should yield zero. More specifically, if the current in lead α is given by J α = ẋF x +ẏF y +żF z , and the lead is infinite along the z direction for instance, then we should have J α = żF z and ẋF x = ẏF y = 0. So in each lead, one may take three partial traces along and perpendicular to the lead direction and confirm these relations. The total trace should still yield the correct result. One final note is the matrices ofΦ ij andΨ ijk are invariant under permutations of their indices.
Appendix C.1. Lead Self-Energies σ α and Escape Rates Γ α The surface Green functions of the leads were defined by an inverse, as stated in Equation (9). The self-energy for lead α is accordingly defined as σ α,ij (ω) = V α,ik γ α,kk (ω) V α,jk where the two indices (k, k ) refer to the lead degrees of freedom. Note the matrices V α connecting the device to lead α need not be square.
The escape rates were defined by: For Langevin thermostats with white noise, these results do not depend on the thermostat damping factors γ.

Appendix C.2. Noise Autocorrelation Functions
Let us start with the noise autocorrelation that will appear in the calculation of displacement-noise correlation Z α . This autocorrelation can be obtained from the fluctuationdissipation theorem, a relation that 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 [22] 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 Novikov-Furutsu-Donsker (NFD) identity (see Appendix D) 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: 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 A, Equation (A2). 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.
Using the equation of motion (18), the expression defining this correlation function to second order in ∂A/∂Y can be derived to satisfy Equation (25). After substitution, the explicit relations become: where, as usual, the summation over repeated indices is implied. Figure A1. Feynman diagrams associated with Z α up to second order in anharmonic forces. Dashed lines represent the phonon Green's function G, the circle represents Γ α , the triangle represents the third-order vertexΨ and the square the fourth-order vertex χ. The thick line with opposite arrows represents the displacement autocorrelation C(ω). Frequency must be conserved at each vertex.

Appendix C.4. Displacement Autocorrelations C(ω)
These functions, whose trace represents the number of excitations in the device, are defined via Equation (27). The static distortion in x x does not contribute to the current and will thus be omitted. To the second order in A, the expression for C is given by: kk (ω) = 1 6χ klmnχk l m n 1,2 C ll (ω − ω 1 − ω 2 ) C mm (ω 1 ) C nn (ω 2 ) (A18) Figure A2. Feynman diagrams associated with C(ω) up to the second order in anharmonic forces. Conventions are the same as in Figure A1.

Appendix D. Statement and Proof of the Novikov-Furutsu-Donsker (NFD) Relation
The Novikov-Furutsu-Donsker relation [19,[23][24][25] relates the correlation function of any functional A[η] of the noise η to the noise-noise correlation times the expectation value of the derivative of that functional: implying that in the frequency domain, we have This can easily be extended to higher dimensions where η is an array:

Appendix E. Heat Current within the Harmonic Approximation
For a harmonic system (G = G = G), we have derived the following expressions: The frequency in the denominator will cancel the frequency in the current coming from the time derivative of positions, so that the harmonic current becomes: Below we will show that this current satisfies both detailed balance and current conservation.
Detailed balance: If all leads are at the same temperature, using Equation (A13), we have Tr Making this substitution for (G) in the heat current formula, we end up with: This is equation manifestly shows detailed balance as it is linear in temperature differences. It is also easy to see that this harmonic part of the heat current satisfies current conservation, due to the antisymmetric (under exchange of α and β) form of the following sum: In a two-terminal device geometry (β = R; α = L), the expression for the current reduces to the well-known formula: In general, Tr [Γ α G Γ β G † ] maybe interpreted as the harmonic transmission from lead α to lead β. We see that even in the non-equilibrium regime (large ∆T), the harmonic approximation leads to the same transmission function regardless of how large the temperature difference is. A simple extension to quantum case: Given the correspondence between the quantum and classical versions of noise autocorrelation, to recover the quantum limit, one may replace f = k B T/ω byh(1 + f BE (hω/k B T)). The constant term 1 is irrelevant and disappears due to the principle of detailed balance, and the temperature difference is replaced by the difference in the distribution functions times the phonon energy: k B (T α − T β ) <=> hω( f α − f β ).