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

Abstract

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.

1. 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 $\mathbb{T}=\mathrm{Tr}\left[G{\Gamma }_L{G}^{†}{\Gamma }_R\right]$, where G and $\Gamma$ 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 $\omega$ 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 high-temperature 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).

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 lowest-order 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 NEGF-based 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.

2. 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 $\alpha$ is represented by Equation (2):

$$\begin{array}{cc}\hfill H_D& =\sum _{i\in D}\frac{p_i^2}{2m_i}+\sum _{ij\in D}\frac{1}{2!}{\varphi }_{ij}{u}_i{u}_j+\sum _{ijk\in D}\frac{1}{3!}{\psi }_{ijk}{u}_i{u}_j{u}_k+\dots \hfill \end{array}$$

(1)

$$\hfill H_{\alpha ,D}& =\sum _{i\in D}\sum _{l\in \alpha }{W}_{\alpha ,il}{u}_i{u}_l\hfill$$

(2)

The dynamical variable $u_i\left(t\right)$ 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 $\alpha$ are semi-infinite harmonic chains. After the standard change of variables to $x_i\left(t\right)=\sqrt{m_i}{u}_i\left(t\right)$, and to

$${\Phi }_{ij}=\frac{{\varphi }_{ij}}{\sqrt{m_i{m}_j}}=\frac{1}{\sqrt{m_i{m}_j}}\frac{{\partial }^2{H}_D}{\partial u_i\partial u_j}(u=0)$$

$V_{\alpha ,il}=W_{\alpha ,il}/\sqrt{m_i{m}_{l\alpha }}$ and ${\Phi }_{\alpha ,l{l}^{\prime }}={\varphi }_{\alpha ,l{l}^{\prime }}/m_{l\alpha }$ we arrive at the following equation of motion for atoms in the central region:

$$\frac{d^{2}x}{d{t}^2}=-\Phi x-\sum _{\alpha }{V}_{\alpha }{x}_{\alpha }+a$$

(3)

Note we have used capitalized Greek letters ($\Phi ,\Psi ,\dots$) for mass-rescaled force constants, and lower-case Greek letters ($\varphi ,\psi ,\dots$) for the bare potential energy derivatives. The letter $\alpha$ refers to the leads, and the dynamical variable $x=(x_1,\dots ,x_N)$ can be thought of as an array containing the displacements of all the atoms in the central region (also called device), $\Phi$ as the force constant matrix between such atoms, and $V_{\alpha }$ as the force constant matrix connecting atoms of the lead $\alpha$ to atoms in the device. Finally, $a=-1/2\Psi x^2+\dots$ 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 ${\gamma }_{\alpha }$ and noise ${\zeta }_{\alpha }$. The equation of motion for atoms with displacements $x_{\alpha }$ is as follows:

$$\frac{d^2{x}_{\alpha }}{d{t}^2}=-{\Phi }_{\alpha }{x}_{\alpha }-V_{\alpha }^{T}x-{\gamma }_{\alpha }\frac{d{x}_{\alpha }}{dt}+{\zeta }_{\alpha }$$

(4)

where the superscript T stands for transpose. Here again all atomic degrees of freedom are stored in the array $x_{\alpha }=(x_{1\alpha },x_{2\alpha },\dots )$. The force constants $Ph{i}_{\alpha }$ 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_{\alpha }$ 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:

$$X\left(\omega \right)=\int dtx\left(t\right)e^{i\omega t};x\left(t\right)=\int \frac{d\omega }{2\pi }X\left(\omega \right)e^{-i\omega t}$$

(5)

In the frequency domain, the equations of motion (3) and (4) become:

$$\begin{array}{cc}\hfill -{\omega }^{2}X& =-\Phi X-{\mathsf{\Sigma }}_{\alpha }{V}_{\alpha }{X}_{\alpha }+A\hfill \end{array}$$

(6)

$$-{\omega }^2{X}_{\alpha }=-{\Phi }_{\alpha }{X}_{\alpha }-V_{\alpha }^{T}X+i\omega {\gamma }_{\alpha }{X}_{\alpha }+{\zeta }_{\alpha }\left(\omega \right)$$

(7)

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_{\alpha }\left(\omega \right)$ be the retarded (causal) Green’s function associated with the lead $\alpha$. The positivity of the damping factor ${\gamma }_{\alpha }$ insures causality. The solution to Equation (7) after the transients have decayed to zero, can be written as:

$$X_{\alpha }=g_{\alpha }\left(\omega \right)({\zeta }_{\alpha }-V_{\alpha }^{T}X)$$

(8)

where

$$g_{\alpha }^{-1}=[-{\omega }^2-i\omega {\gamma }_{\alpha }+{\Phi }_{\alpha }]$$

(9)

In Equation (6) we need $-V_{\alpha }{X}_{\alpha }$, which we can obtain from Equation (8):

$$-V_{\alpha }{X}_{\alpha }={\eta }_{\alpha }+{\sigma }_{\alpha }X,$$

(10)

Here ${\eta }_{\alpha }\left(\omega \right)=-V_{\alpha }{g}_{\alpha }\left(\omega \right){\zeta }_{\alpha }\left(\omega \right)$, and ${\sigma }_{\alpha }=V_{\alpha }{g}_{\alpha }{V}_{\alpha }^T$. Likewise, defining the retarded Green’s function of the central region as $G^{-1}=[-{\omega }^2+\Phi -{\sum }_{\alpha }{\sigma }_{\alpha }\left(\omega \right)]$, we can write the solution to the central region as:

$$X\left(\omega \right)=G\left(\omega \right)\left[\sum _{\alpha }{\eta }_{\alpha }\left(\omega \right)+A\left(\omega \right)\right]$$

(11)

The function ${\sigma }_{\alpha }\left(\omega \right)=V_{\alpha }{g}_{\alpha }\left(\omega \right)V_{\alpha }^T$ is traditionally called the self-energy of lead $\alpha$, 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 ${\omega }^2$, and its imaginary part, divided by $2\omega$, 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 ${\eta }_{\alpha }=-V_{\alpha }{g}_{\alpha }{\zeta }_{\alpha }$.

3. 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 $\langle \cdots \rangle$ over all realizations of the forces subject to the constraint imposed by the fluctuation–dissipation (FD) theorem [17]:

$$\langle {\zeta }_{i,\alpha }\left(t\right)\rangle =0$$

$$\langle {\zeta }_{i,\alpha }\left(t\right){\zeta }_{j,{\alpha }^{\prime }}\left(t^{\prime }\right)\rangle =2{\gamma }_{\alpha }{k}_B{T}_{\alpha }\delta (t-t^{\prime }){\delta }_{\alpha ,{\alpha }^{\prime }}{\delta }_{ij}$$

(12)

The noise is white and different sites $i,j$ of leads $\alpha ,{\alpha }^{\prime }$ 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.

3.1. Entropy Generation Rate

An important quantity of interest is the entropy generation rate in the device, which can be expressed as: $\dot{S}=-{\sum }_{\alpha }\langle j_{\alpha }\rangle /T_{\alpha }$. 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.

3.2. Heat Current

The main quantity of interest is the heat current. The heat from lead $\alpha$ 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 $\alpha$ 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: $j_{\alpha }\left(t\right)=\mathrm{Tr}[\dot{x}{(-V_{\alpha }{x}_{\alpha })}^T]$, where the trace is taken over the device degrees of freedom. Substituting the expressions for x and $x_{\alpha }$ 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 $\eta$, 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_{\alpha }\left(t\right)$ and call it $J_{\alpha }(\Omega )$. Note due to time integration, we must have $\langle J_{\alpha }(\Omega =0)\rangle =\tau \langle j_{\alpha }\rangle$, where $\tau \to \infty$ is the time integration window. Finally the static DC current is given by:

$$\langle j_{\alpha }\rangle =\int \frac{d\omega }{2\pi \tau }\omega \mathrm{Tr}\left[\Im \langle X\left(\omega \right){\eta }_{\alpha }^{†}\left(\omega \right)\rangle -\langle X\left(\omega \right)X^{†}\left(\omega \right)\rangle {\Gamma }_{\alpha }\left(\omega \right)/2\right]$$

(13)

where we used the notation ${\Gamma }_{\alpha }=-i({\sigma }_{\alpha }-{\sigma }_{\alpha }^{†})=2\Im \left({\sigma }_{\alpha }\right)$ for twice the imaginary part of the lead $\alpha$ self-energy. The DC response is found by taking the $\Omega \to 0$ limit. Note that because we are interested in the DC response, only diagonal terms of correlations ($\Omega =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 ${\sigma }_{\alpha }\left(\omega \right)$, followed by a frequency integration. Note that this current from lead $\alpha$ is the sum of two terms: the first one, proportional to $Z_{\alpha }=\langle X\left(\omega \right){\eta }_{\alpha }^{†}\left(\omega \right)\rangle$, is the work per unit time of the stochastic forces on the device, while the second one $\langle X\left(\omega \right)X^{†}\left(\omega \right)\rangle {\Gamma }_{\alpha }^{†}$, 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 $\langle \rangle$ 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].

4. 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 $\langle j_{\alpha }\rangle$ should be identically zero for all leads $\alpha$. The second constraint is that of current conservation, which in steady state $\Omega \to 0$, and under no additional heat generation in the device, reduces to ${\mathsf{\Sigma }}_{\alpha }\langle j_{\alpha }\rangle =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.

5. 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\left(\omega \right)=X\left(\omega \right)-\langle X\rangle$ or $y\left(t\right)=x\left(t\right)-\langle x\rangle$ which has zero average by construction. The resulting nonlinear equations satisfied by $\langle X\rangle$ 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).

$$\langle \partial \mathcal{V}/\partial x\rangle =\Phi \langle x\rangle +\Psi \langle xx\rangle /2+\dots =0$$

(14)

The resulting equations will depend on correlations such as $\langle Y_i{Y}_j\rangle$, $\langle Y_i{Y}_j{Y}_k\rangle$ 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 ($\Phi \to \overline{\Phi }$ etc…). While the variable X satisfies the equation of motion:

$$-{\omega }^{2}X=-\partial \mathcal{V}/\partial X-\sum _{\alpha }{V}_{\alpha }{X}_{\alpha },$$

the new variable Y satisfies:

$$-{\omega }^{2}Y=-\partial \mathcal{V}/\partial X+\langle \partial \mathcal{V}/\partial X\rangle -\sum _{\alpha }{V}_{\alpha }{X}_{\alpha }$$

Next we will linearize the forces with respect to y:

$$\begin{array}{c}\hfill -\frac{\partial \mathcal{V}}{\partial y}=-{\left(\frac{\partial \mathcal{V}}{\partial y}\right)}_{y=0}-{\left(\frac{{\partial }^2\mathcal{V}}{\partial y^2}\right)}_{y=0}y+a\left(y\right)=-\overline{\Phi }y+a\left(y\right)\end{array}$$

where we have set ${\left(\frac{\partial \mathcal{V}}{\partial y}\right)}_{y=0}=0$ to define the thermal expansion $\langle x\rangle =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 non-equilibrium mean-field approximation (NEMF) also detailed in Section 6. Next, we define a renormalized Green’s function using renormalized force constants $\overline{\Phi }={\left(\frac{{\partial }^2\mathcal{V}}{\partial y^2}\right)}_{y=0}$ as

$${\mathcal{G}}^{-1}=[-{\omega }^2+\overline{\Phi }-\sum _{\alpha }{\sigma }_{\alpha }]$$

(15)

With this GF, the displacements Y satisfy

$$Y=\mathcal{G}(\sum _{\alpha }{\eta }_{\alpha }+A)$$

(16)

Note one can add any constant $\lambda$ to the force constant $\overline{\Phi }$ in the above GF, provided $\lambda 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.

6. Force Constant Renormalization

Given the form of the above equations, we can add $-Y\langle \frac{\partial A}{\partial Y}\rangle$ to A and add $-\langle \frac{\partial A}{\partial Y}\rangle$ to $1/\mathcal{G}$ so that now the renormalized GF becomes:

$${\mathbb{G}}^{-1}={\mathcal{G}}^{-1}-\langle \frac{\partial A}{\partial Y}\rangle =[-{\omega }^2+\overline{\Phi }-\langle \frac{\partial A}{\partial Y}\rangle +\sum _{\alpha }{\sigma }_{\alpha }]$$

(17)

while the renormalized anharmonic force now becomes $\mathbb{A}=A-Y\langle \frac{\partial A}{\partial Y}\rangle$ in the right hand side of Equation (16). This renormalization of harmonic force constants captures a major part of anharmonicity (because $\langle \partial \mathbb{A}/\partial Y\rangle =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_{\alpha }$ disappears by construction since $\langle \partial \mathbb{A}/\partial Y\rangle =0$, leading to the smallest variance and higher moments of $\partial \mathbb{A}/\partial 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:

$$\begin{array}{|c|}\hline Y=\mathbb{G}({\displaystyle \sum _{\alpha }}{\eta }_{\alpha }+\mathbb{A})\\ \hline\end{array}$$

(18)

The Feynman diagram associated with the new Green’s function $\mathbb{G}$ is shown in Figure 2.

7. 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_{\alpha }$ 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=-\psi y^2/2-\chi y^3/6$.

Let us start with the noise autocorrelation, which will appear in the calculation of displacement-noise correlation $Z_{\alpha }$. In the frequency domain, using the fluctuation–dissipation theorem Equation (12), one can derive (see Appendix C.2):

$$\begin{array}{c}\hfill \begin{array}{|c|}\hline \langle {\eta }_{\alpha }\left(\omega \right){\eta }_{\alpha }^{†}\left(\omega \right)\rangle ={\Gamma }_{\alpha }\left(\omega \right)\frac{k_B{T}_{\alpha }}{\omega }\tau ={\Gamma }_{\alpha }\left(\omega \right)f_{\alpha }\tau \\ \hline\end{array}\end{array}$$

(19)

where, for brevity, we have replaced the “occupation factor” $k_B{T}_{\alpha }/\omega$ by $f_{\alpha }$, and $\tau$ represents the integration time which goes to infinity and cancels the $\tau$ in the expression for the current $j_{\alpha }=J_{\alpha }/\tau$. With a Factor of ℏ in the denominator of $f_{\alpha }$, 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 $\gamma$, and thus we adopt this type of noise for the thermostats.

Next, we need to calculate displacement-noise correlations: $Z_{\alpha }=\langle Y{\eta }_{\alpha }^{†}\rangle$. Since $-i\omega Y$ is a velocity, $\omega \Im \left(Z_{\alpha }\right)$ is the power exerted by the random force $\eta$ 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=\mathbb{A}=0;G=\mathcal{G}=\mathbb{G}$), using the equation motion (18), this expression is simplified to:

$$Z_{\alpha }^H\left(\omega \right)=\mathcal{G}\langle {\eta }_{\alpha }{\eta }_{\alpha }^{†}\rangle =\mathcal{G}{\Gamma }_{\alpha }{f}_{\alpha }\tau$$

(20)

For non-zero anharmonicity, the three GFs are different, and adopting $\mathbb{G}$ as the reference GF, we have the NEMF approximation to the displacement-noise correlation function as:

$$Z_{\alpha }^{NEMF}\left(\omega \right)=\mathbb{G}\langle {\eta }_{\alpha }{\eta }_{\alpha }^{†}\rangle =\mathbb{G}{\Gamma }_{\alpha }{f}_{\alpha }\tau$$

(21)

To go one step further and include the effect of anharmonicity in $Z_{\alpha }$, 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\left[\eta \right]$, we have

$$\begin{array}{|c|}\hline \langle f\left[\eta \right]{\eta }_{\alpha }^{†}\left(\omega \right)\rangle =\langle \frac{\delta f\left[\eta \right]}{\delta {\eta }_{\alpha }}\rangle \langle {\eta }_{\alpha }{\eta }_{\alpha }^{†}\rangle \\ \hline\end{array}$$

(22)

It has the advantage of lowering the powers of $\eta$ in f. Using this theorem, we have: $\langle Y{\eta }_{\alpha }^{†}\rangle =\langle \frac{\partial Y}{\partial {\eta }_{\alpha }}\rangle \langle {\eta }_{\alpha }{\eta }_{\alpha }^{†}\rangle$.

From the equation of motion Equation (18) and the chain rule, we find

$$\begin{array}{c}\hfill \frac{\partial Y}{\partial {\eta }_{\alpha }}={(1-\mathbb{G}\frac{\partial \mathbb{A}}{\partial Y})}^{-1}\mathbb{G}=\mathbb{G}+\mathbb{G}\frac{\partial \mathbb{A}}{\partial Y}\mathbb{G}+\mathbb{G}\frac{\partial \mathbb{A}}{\partial Y}\mathbb{G}\frac{\partial \mathbb{A}}{\partial Y}\mathbb{G}+\dots \end{array}$$

so that finally,

$$\begin{array}{ccc}\hfill Z_{\alpha }& =& \langle Y{\eta }_{\alpha }^{†}\rangle =\langle \frac{\partial Y}{\partial {\eta }_{\alpha }}\rangle \langle {\eta }_{\alpha }{\eta }_{\alpha }^{†}\rangle =\langle {\left(1-\mathbb{G}\frac{\partial \mathbb{A}}{\partial Y}\right)}^{-1}\rangle Z_{\alpha }^{NEMF}\hfill \end{array}$$

(23)

$$\begin{array}{ccc}& =& \left(1+\mathbb{G}\langle \frac{\partial \mathbb{A}}{\partial Y}\mathbb{G}\frac{\partial \mathbb{A}}{\partial Y}\rangle +\mathbb{G}\langle \frac{\partial \mathbb{A}}{\partial Y}\mathbb{G}\frac{\partial \mathbb{A}}{\partial Y}\mathbb{G}\frac{\partial \mathbb{A}}{\partial Y}\rangle +\dots \right)Z_{\alpha }^{NEMF}\hfill \end{array}$$

(24)

Note that, by construction, the first term involving $\langle \partial \mathbb{A}/\partial Y\rangle$ is identically zero. This justifies the choice of $\mathbb{G}$ rather than $\mathcal{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 $\mathbb{G}$).

The next-order correction consists in adding the effect of anharmonicity to second order, i.e., writing the displacement-noise correlation function as:

$$Z_{\alpha }\approx \left(1+\mathbb{G}\langle \frac{\partial \mathbb{A}}{\partial Y}\mathbb{G}\frac{\partial \mathbb{A}}{\partial Y}\rangle \right)Z_{\alpha }^{NEMF}$$

(25)

with $\mathbb{A}=A-Y\langle \frac{\partial A}{\partial Y}\rangle =-\frac{1}{2}\overline{\Psi }YY-\frac{1}{6}\chi (YYY-3Y\langle YY\rangle )$ and $\frac{\partial \mathbb{A}}{\partial Y}=\frac{\partial A}{\partial Y}-\langle \frac{\partial A}{\partial Y}\rangle =-\overline{\Psi }Y-\frac{1}{2}\chi (YY-\langle YY\rangle )$. We can note that the major part of the quartic anharmonicity has been removed, and the expansion is in powers of $\partial \mathbb{A}/\partial 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_{\alpha }$, 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:

$$Z_{\alpha }\approx \left(1+\langle \mathbb{G}\overline{\Psi }Y\mathbb{G}\overline{\Psi }Y\rangle \right)Z_{\alpha }^{NEMF}$$

(26)

An explicit form of this equation is provided in the following section and in Appendix C.3 Equation (A17).

8. Displacement Autocorrelations

Finally, the last correlation function needed is $C\left(\omega \right)=\langle Y{Y}^{†}\rangle$. Note the autocorrelation matrix C is Hermitian. The function $\omega C\left(\omega \right)$ can be interpreted as the (non-equilibrium) number of excitations of frequency $\omega$ 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\left({\mathsf{\Sigma }}_{\alpha }{\Gamma }_{\alpha }\right)G^{†}=2\Im \left(G\right)$, which is the total (equilibrium) number of excitations in the device. In the heat current, this term appears as $-\omega C{\Gamma }_{\alpha }/2$ and can be interpreted as the heat current going from (because of the negative sign) the device into this lead as ${\Gamma }_{\alpha }/2\omega$ is the escape rate into the lead $\alpha$.

Similar to the treatment of $Z_{\alpha }$, we will start with the equation of motion, Equation (18), and the FD theorem, Equation (19), to write C as:

$$\begin{array}{cc}& C\left(\omega \right)=\langle Y\left(\omega \right)Y^{†}\left(\omega \right)\rangle =\hfill \\ & \mathbb{G}\left(\omega \right)\left[{\mathsf{\Sigma }}_{\alpha }\left(\langle {\eta }_{\alpha }{\eta }_{\alpha }^{†}\rangle +\langle \mathbb{A}{\eta }_{\alpha }^{†}\rangle +\langle {\eta }_{\alpha }{\mathbb{A}}^{†}\rangle \right)+\langle \mathbb{A}{\mathbb{A}}^{†}\rangle \right]{\mathbb{G}}^{†}\left(\omega \right)\hfill \end{array}$$

(27)

The first term is the NEMF contribution and is reduced to a known result, with now the renormalized Green’s function $\mathbb{G}$ being used instead of the standard harmonic one (G or $\mathcal{G}$), in order to include thermal effects to some extent:

$$C^{NEMF}\left(\omega \right)=\mathbb{G}{\mathsf{\Sigma }}_{\alpha }\langle {\eta }_{\alpha }{\eta }_{\alpha }^{†}\rangle {\mathbb{G}}^{†}={\mathsf{\Sigma }}_{\alpha }\left(\mathbb{G}{\Gamma }_{\alpha }{\mathbb{G}}^{†}\right)f_{\alpha }\tau$$

(28)

The NEMF approximation, which consists in using $\mathbb{G}$ and neglecting the contributions of anharmonicity included in $\mathbb{A}$, is very similar to the harmonic approximation. Within this approximation, where $C\approx C^{NEMF}$ and $Z_{\alpha }\approx Z_{\alpha }^{NEMF}$, the transmission becomes $\mathrm{Tr}\left[\mathbb{G}{\Gamma }_L{\mathbb{G}}^{†}{\Gamma }_R\right]$, very similar to the harmonic result (see Appendix E), but with the GF substituted by the renormalized $\mathbb{G}$, and includes the effect of (quartic) anharmonicity to the lowest order.

Linear terms in $\mathbb{A}$ lead to terms similar to $Z_{\alpha }$, which has already been discussed. The only remaining difficulty is with the $\langle \mathbb{A}{\mathbb{A}}^{†}\rangle$ 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_{\alpha }$ in Equation (26), we have to use:

$${\langle \mathbb{A}{\mathbb{A}}^{†}\rangle }^{\left(3\right)}=C^{\left(3\right)}=\frac{1}{4}{\overline{\Psi }}^2\sum _{{\omega }_1}C(\omega -{\omega }_1)C\left({\omega }_1\right)$$

It can be shown that this approximation, taken to a self-consistent level satisfies current conservation, meaning ${\sum }_{\alpha }{j}_{\alpha }=0$, however the spectral components of the current: $\Im Z_{\alpha }\left(\omega \right)-C\left(\omega \right){\Gamma }_{\alpha }\left(\omega \right)/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 $\mathbb{G}$ are:

$$\begin{array}{cc}& Z_{\alpha }=Z_{\alpha }^{NEMF}+\mathbb{G}({\mathsf{\Sigma }}^{\left(3\right)}+{\mathsf{\Sigma }}^{\left(4\right)})Z_{\alpha }\hfill \\ & {\mathsf{\Sigma }}^{\left(3\right)}\left(\omega \right)=2{\overline{\Psi }}^2\sum _{{\omega }_1}\mathbb{G}(\omega -{\omega }_1)C\left({\omega }_1\right)\hfill \\ & {\mathsf{\Sigma }}^{\left(4\right)}\left(\omega \right)=\frac{3}{2}{\overline{\chi }}^2\sum _{{\omega }_1,{\omega }_2}\mathbb{G}(\omega -{\omega }_1-{\omega }_2)C\left({\omega }_1\right)C\left({\omega }_2\right)\hfill \\ & C=\sum _{\alpha }(\delta Z_{\alpha }{\mathbb{G}}^{†}+\mathbb{G}\delta Z_{\alpha }^{†})+C^{NEMF}+\mathbb{G}P{\mathbb{G}}^{†}\hfill \\ & \delta Z_{\alpha }=Z_{\alpha }-Z_{\alpha }^{NEMF}\hfill \\ & P=C^{\left(3\right)}+C^{\left(4\right)}\hfill \\ & C^{\left(3\right)}=\frac{1}{2}{\overline{\Psi }}^2\sum _{{\omega }_1}C(\omega -{\omega }_1)C\left({\omega }_1\right)\hfill \\ & C^{\left(4\right)}=\frac{1}{6}{\overline{\chi }}^2\sum _{{\omega }_1,{\omega }_2}C(\omega -{\omega }_1-{\omega }_2)C\left({\omega }_1\right)C\left({\omega }_2\right)\hfill \end{array}$$

(29)

The equations for $Z_{\alpha }$ and C can also be represented using Feynman diagrams as shown in Figure A1 and Figure A2 in Appendix C.3 and Appendix 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_{\alpha }$ 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_{\alpha }$ requires $\partial \mathbb{A}/\partial Y$, the terms C require $\mathbb{A}$ itself, but these equations contain only second powers of $\mathbb{A}$ and $\partial \mathbb{A}/\partial Y$. Once iterations converge, the obtained C and $Z_{\alpha }$ 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_{\alpha }=k_{B}T/\omega$ 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 $\Delta T$ is held small). This is in agreement with previous MD simulations [13,14].

9. 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 of ℏ in the formalism and $\omega$ 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 $\mathbb{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 $\mathbb{A}$ up to quartic terms and only included up to second powers of $\mathbb{A}$ and $\partial \mathbb{A}/\partial 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, $\mathbb{G}$, has two corrections: one due to thermal expansion implying changes in bond length and strength ($\overline{\Phi }$ instead of $\Phi$), and the other due to thermal fluctuations about the average position ($\langle \partial \mathbb{A}/\partial Y\rangle$, 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 fluctuation–dissipation 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_i^{NEMD}\approx F_i^{Harmonic}=-{\overline{\Phi }}_{ij}{y}_j$ in order to extract the effective (non-equilibrium) harmonic force constants $\overline{\Phi }$. The remainder can then be defined as the anharmonic force: $F_i^{NEMD}=-{\overline{\Phi }}_{ij}{y}_j+a_i\left(y\right)$, 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 $\mathbb{G}\langle \partial \mathbb{A}/\partial Y\rangle \approx 1$, signaling resonances, in which case the Taylor expansion in Equation (24) is not appropriate.

Applications to nanoscale systems will appear in future publications.