††thanks<EMAIL_ADDRESS>
# Theory of Non-equilibrium Heat transport in anharmonic multiprobe systems at
high temperatures

Keivan Esfarjani  Department of Mechanical and Aerospace Engineering,
University of Virginia, Charlottesville, Virginia 22904  Department of
Materials Science and Engineering, University of Virginia, Charlottesville,
Virginia 22904  Department of Physics, University of Virginia,
Charlottesville, Virginia 22904

###### Abstract

We consider the problem of heat transport by vibrational modes (conduction)
between Langevin thermostats connected by a central device. The latter is
anharmonic and can be subject to large temperature differences and thus be out
of equilibrium. We develop a self-consistent Green’s function formalism to
describe high-temperature and non-equilibrium transport, and derive a formula
for the heat current for up to quartic anharmonicity (4th-order in the
potential energy). We show the importance of including quartic terms in the
anharmonic potential in order to properly describe thermal expansion and
temperature dependence to leading order in anharmonicity. This formalism paves
the way for accurate and efficient modeling of thermal transport in highly
non-equilibrium situations beyond perturbation theory.

##  I  Introduction

Transport theories of non-interacting quantum systems based on the Keldysh
formalismKeldysh (1965), which treats nonequilibrium 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 done in 1971 by
Caroli et al.Caroli  _et al._ (1971) 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={\rm Tr}[G\Gamma_{L}G^{\dagger}\Gamma_{R}]$ , where $G$ and $\Gamma$
are respectively the retarded Green’s function and the escape rates to the
leads, was established by Meir and WingreenMeir and Wingreen (1992) and in a
similar form by PastawskiPastawski (1991) 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 strong interactions and large driving fields.

The basic geometry of our problem is a multi-probe one where the system in
which scatterings occur is connected to multiple reservoirs or contacts, which
impose their chemical potential and/or temperature, and cause flow of charge
or heat carriers (see Fig. 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 done and we will be dealing with strictly 1D semi-infinite
leads, although the central device is arbitrary in shape and structure and
maybe connected to multiple 1D probes. In this paper, we will be interested in
transport of anharmonic phonons, or more generally, vibrational modes, in
mesoscopic systems under large temperature differences. 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. 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
interfacial thermal conductance (see Fig. 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 MingoMingo (2006) and separately by WangWang  _et al._ (2006) 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 TianDai and Tian (2020) 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. Guo  _et al._ (2020), based on previous work by LuisierLuisier
(2012). Polanco has a recent review of these methods based on NEGFPolanco
(2021). 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
regimeSääskilahti  _et al._ (2013) 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 1970Bolsterli  _et
al._ (1970). In this method, as shown in Fig. 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 interfaceLe  _et al._ (2017a);
Feng  _et al._ (2019), 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.Sääskilahti
_et al._ (2014) extended the harmonic formulation of the transmission function
based on the actual MD trajectories by Chalopin et al.Chalopin  _et al._
(2012); Chalopin and Volz (2013), 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.

##  II Dynamics

We start by defining our model and the assumptions. A multiprobe geometry is
assumed as shown in Fig. 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
profile. The Hamiltonian of the device (D) and its coupling to the lead
$\alpha$ are:

$\displaystyle H_{D}$ $\displaystyle=\sum_{i\in
D}\frac{p_{i}^{2}}{2m_{i}}+\sum_{ij\in
D}\,\frac{1}{2!}\phi_{ij}\,u_{i}\,u_{j}$ (1) $\displaystyle+\sum_{ijk\in
D}\,\frac{1}{3!}\psi_{ijk}\,u_{i}\,u_{j}\,u_{k}+...$ (2) $\displaystyle
H_{\alpha,D}$ $\displaystyle=\sum_{i\in
D}\sum_{l\in\alpha}\,W_{\alpha,il}\,u_{i}\,u_{l}$ (3)

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 $\alpha$ are semi-infinite chains following
Langevin dynamics. After the standard change of variables to
$x_{i}(t)=\sqrt{m_{i}}u_{i}(t)$, and to

$\Phi_{ij}=\frac{\phi_{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_{\alpha}}$ and
$\Phi_{\alpha,ll^{\prime}}=\phi_{\alpha,ll^{\prime}}/m_{\alpha}$ we arrive at
the following equation of motion for atoms in the central region:

$\frac{d^{2}x}{dt^{2}}=-\Phi\,x-\sum_{\alpha}V_{\alpha}\,x_{\alpha}+a$ (4)

Note we have used capitalized Greek letters ($\Phi,\Psi,...$) for mass-
rescaled force constants, and lower-case Greek letters ($\phi,\psi,...$) for
the bare potential energy derivatives. The letter $\alpha$ 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), $\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}+...$ 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, where a set of identical
coupled harmonic oscillators are subject to damping $\gamma_{\alpha}$ and
noise $\zeta_{\alpha}$ as follows:

$\frac{d^{2}x_{\alpha}}{dt^{2}}=-\Phi_{\alpha}x_{\alpha}-V_{\alpha}^{T}x-\gamma_{\alpha}\frac{dx_{\alpha}}{dt}+\zeta_{\alpha}$
(5)

where the superscript $T$ stands for transpose. The force constants
$Phi_{\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 Eq. 4 using the Green’s
function method. To this end, we start by taking the Fourier transform of the
above two equations according to:

$X(\omega)=\int dt\,x(t)\,e^{i\omega
t};\,\,x(t)=\int\,\frac{d\omega}{2\pi}\,X(\omega)\,e^{-i\omega t}$ (6)

The equations of motion become:

$\displaystyle-\omega^{2}X$ $\displaystyle=-\Phi
X-\Sigma_{\alpha}V_{\alpha}\,X_{\alpha}+A$ (7)
$\displaystyle-\omega^{2}X_{\alpha}$
$\displaystyle=-\Phi_{\alpha}X_{\alpha}-V_{\alpha}^{T}X+i\omega\gamma_{\alpha}X_{\alpha}+\zeta_{\alpha}(\omega)$
(8)

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_{\alpha}(\omega)$ 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 Eq. 8 after the transients have decayed to
zero, can be written as:

$X_{\alpha}=g_{\alpha}(\omega)\,(\zeta_{\alpha}-V_{\alpha}^{T}X)$ (9)

where

$g_{\alpha}^{-1}=[-\omega^{2}-i\omega\gamma_{\alpha}+\Phi_{\alpha}]$ (10)

In Eq. 7 we need $-V_{\alpha}X_{\alpha}=\eta_{\alpha}+\sigma_{\alpha}X$, which
we obtain from Eq. 9. Here
$\eta_{\alpha}(\omega)=-V_{\alpha}g_{\alpha}(\omega)\zeta_{\alpha}(\omega)$,
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}(\omega)]$, we can write
the solution to the central region as:

$X(\omega)=G(\omega)\,\Big{[}\sum_{\alpha}\eta_{\alpha}(\omega)\,+A(\omega)\Big{]}$
(11)

The function
$\sigma_{\alpha}(\omega)=V_{\alpha}\,g_{\alpha}(\omega)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}$.

##  III 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, we need to perform an
“ensemble” average, denoted by $\langle\ldots\rangle$ over all realizations of
the forces subject to the constraint imposed by the fluctuation-dissipation
(FD) theorem:

$\langle\zeta_{i,\alpha}(t)\rangle=0$
$\langle\zeta_{i,\alpha}(t)\zeta_{j,\alpha^{\prime}}(t^{\prime})\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 with the
constraints imposed by the FD theorem..

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 done 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}(t)={\rm
Tr}[\dot{x}\,(-V_{\alpha}x_{\alpha})^{T}]$, where the trace is taken over the
device degrees of freedom after using Eq. 9. Since the current depends on the
stochastic functions $\zeta$, 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}(t)$ 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:

$\langle J_{\alpha}(\Omega)\rangle=\int\frac{d\omega}{2\pi}\,\omega\,{\rm
Tr}\,\Big{[}\Im\langle
X(\omega)\eta_{\alpha}^{\dagger}(\omega-\Omega)\rangle-\langle
X(\omega)X^{\dagger}(\omega-\Omega)\rangle\,\Gamma_{\alpha}(\omega-\Omega)/2\Big{]}$
(13)

where we used the notation
$\Gamma_{\alpha}=-i(\sigma_{\alpha}-\sigma_{\alpha}^{\dagger})=2\Im(\sigma_{\alpha})$
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}(\omega)$, 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(\omega)\,\eta_{\alpha}^{\dagger}(\omega)\rangle$, is the work per unit time
of the stochastic forces on the device, while the second one $\langle
X(\omega)\,X^{\dagger}(\omega)\rangle\Gamma_{\alpha}^{\dagger}$, 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 which 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
approximationsMingo (2006); Wang  _et al._ (2006); Tian  _et al._ (2014).
Another 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}$

##  IV 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 $\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.

##  V 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(\omega)=X(\omega)-\langle X\rangle$ or $y(t)=x(t)-\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+...=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\bar{\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$:

$\displaystyle-\frac{\partial\mathcal{V}}{\partial
y}=-\big{(}\frac{\partial\mathcal{V}}{\partial
y}\big{)}_{y=0}-\big{(}\frac{\partial^{2}\mathcal{V}}{\partial
y^{2}}\big{)}_{y=0}\,y+a(y)=-\bar{\Phi}y+a(y)$

where we have set $(\frac{\partial\mathcal{V}}{\partial y})_{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 VI. Next, we define a renormalized Green’s function using
renormalized force constants
$\bar{\Phi}=\big{(}\frac{\partial^{2}\mathcal{V}}{\partial
y^{2}}\big{)}_{y=0}$ as

$\mathcal{G}^{-1}=[-\omega^{2}+\bar{\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 $\bar{\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.

##  VI 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}+\bar{\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 Eq. 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 $\langle\partial\mathbb{A}/\partial Y\rangle$.

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:

$\boxed{Y=\mathbb{G}(\sum_{\alpha}\eta_{\alpha}+\mathbb{A})}$ (18)

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

Figure 2: Feynman diagram associated with $\mathbb{G}$ represented by thick
dashed lines. The phonon Green’s function $\mathcal{G}$ is shown with thin
solid lines, and the quartic vertex $\chi$ with the solid square. The thick
line with opposite arrows represents the displacement autocorrelation
$C(\omega_{1})$. Internal frequency $\omega_{1}$ is integrated over.

##  VII Displacement-noise correlations

One can see from Eq. 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 Eq. 12, one can derive (see
appendix C.2):

$\displaystyle\boxed{\langle\eta_{\alpha}(\omega)\eta_{\alpha}^{\dagger}(\omega)\rangle=\Gamma_{\alpha}(\omega)\frac{k_{B}T_{\alpha}}{\omega}\tau=\Gamma_{\alpha}(\omega)f_{\alpha}\tau}$
(19)

where, for brevity, we have replaced the “occupation factor”hba
$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$. 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}^{\dagger}\rangle$. Since $-i\omega Y$ is a
velocity, $\omega\Im(Z_{\alpha})$ 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 $\alpha$ 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}(\omega)=\mathcal{G}\langle\eta_{\alpha}\,\eta_{\alpha}^{\dagger}\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}(\omega)=\mathbb{G}\langle\eta_{\alpha}\,\eta_{\alpha}^{\dagger}\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 identityRahimi Tabar
(2019) (for a proof also see appendixD), which states that for any functional
of the white noise $f[\eta]$, we have

$\boxed{\langle
f[\eta]\,\eta_{\alpha}^{\dagger}(\omega)\rangle=\langle\frac{\delta
f[\eta]}{\delta\eta_{\alpha}}\rangle\langle\eta_{\alpha}\,\eta_{\alpha}^{\dagger}\rangle}$
(22)

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

From the equation of motion Eq. 18 and the chain rule, we find

$\displaystyle\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}+...$

so that finally,

$Z_{\alpha}=\langle Y\eta_{\alpha}^{\dagger}\rangle=\langle\frac{\partial
Y}{\partial\eta_{\alpha}}\rangle\langle\eta_{\alpha}\,\eta_{\alpha}^{\dagger}\rangle=\langle\big{(}1-\mathbb{G}\frac{\partial\mathbb{A}}{\partial
Y}\big{)}^{-1}\rangle\,Z_{\alpha}^{NEMF}=\big{(}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+...\big{)}\,Z_{\alpha}^{NEMF}$ (23)

Note that, by construction, the first term involving
$\langle\partial\mathbb{A}/\partial Y\rangle$ is identically zero.

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\big{(}1+\mathbb{G}\langle\frac{\partial\mathbb{A}}{\partial
Y}\mathbb{G}\frac{\partial\mathbb{A}}{\partial
Y}\rangle\big{)}\,Z_{\alpha}^{NEMF}$ (24)

with $\mathbb{A}=A-Y\langle\frac{\partial A}{\partial
Y}\rangle=-\frac{1}{2}\bar{\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=-\bar{\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 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\big{(}1+\langle\mathbb{G}\bar{\Psi}Y\mathbb{G}\bar{\Psi}Y\rangle\big{)}\,Z_{\alpha}^{NEMF}$
(25)

An explicit form of this equation is provided in the later section and in
appendix Eq. 45.

##  VIII Displacement autocorrelations

Finally, the last correlation function needed is $C(\omega)=\langle
YY^{\dagger}\rangle$. Note the autocorrelation matrix $C$ is Hermitian. The
function $\omega C(\omega)$ 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(\Sigma_{\alpha}\Gamma_{\alpha})G^{\dagger}=2\Im(G)$, 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, Eq. 18, and the FD theorem, Eq. 19, to write $C$ as:

$\displaystyle C(\omega)=\langle Y(\omega)Y^{\dagger}(\omega)\rangle=$
$\displaystyle\mathbb{G}(\omega)\Big{[}\Sigma_{\alpha}\big{(}\langle\eta_{\alpha}\eta^{\dagger}_{\alpha}\rangle+\langle\mathbb{A}\eta^{\dagger}_{\alpha}\rangle+\langle\eta_{\alpha}\mathbb{A}^{\dagger}\rangle\big{)}+\langle\mathbb{A}\mathbb{A}^{\dagger}\rangle\Big{]}\mathbb{G}^{\dagger}(\omega)$
(26)

The first term is the NEMF contribution and is reduced to a known result, with
now the renormalized GF $\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}(\omega)=\mathbb{G}\,\Sigma_{\alpha}\langle\eta_{\alpha}\eta^{\dagger}_{\alpha}\rangle\,\mathbb{G}^{\dagger}=\Sigma_{\alpha}\,(\mathbb{G}\Gamma_{\alpha}\mathbb{G}^{\dagger})\,f_{\alpha}\tau$
(27)

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
${\rm Tr}\,[\mathbb{G}\Gamma_{L}\mathbb{G}^{\dagger}\Gamma_{R}]$, 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 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}^{\dagger}\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 Eq. 25, we have to
use:

$\langle\mathbb{A}\mathbb{A}^{\dagger}\rangle^{(3)}=C^{(3)}=\frac{1}{4}\bar{\Psi}^{2}\sum_{\omega_{1}}C(\omega-\omega_{1})C(\omega_{1})$

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}}(\omega)-C(\omega)\Gamma_{\alpha}(\omega)/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:

$\displaystyle
Z_{\alpha}=Z_{\alpha}^{NEMF}+\mathbb{G}(\Sigma^{(3)}+\Sigma^{(4)})Z_{\alpha}$
$\displaystyle\Sigma^{(3)}(\omega)=2\bar{\Psi}^{2}\sum_{\omega_{1}}\mathbb{G}(\omega-\omega_{1})C(\omega_{1})$
$\displaystyle\Sigma^{(4)}(\omega)=\frac{3}{2}\bar{\chi}^{2}\sum_{\omega_{1},\omega_{2}}\mathbb{G}(\omega-\omega_{1}-\omega_{2})C(\omega_{1})C(\omega_{2})$
$\displaystyle C=\sum_{\alpha}(\delta
Z_{\alpha}\,\mathbb{G}^{\dagger}+\mathbb{G}\,\delta
Z_{\alpha}^{\dagger})-C^{NEMF}+\mathbb{G}P\mathbb{G}^{\dagger}$
$\displaystyle\delta Z_{\alpha}=Z_{\alpha}-Z_{\alpha}^{NEMF}$ $\displaystyle
P=C^{(3)}+C^{(4)}$ $\displaystyle
C^{(3)}=\frac{1}{2}\bar{\Psi}^{2}\sum_{\omega_{1}}C(\omega-\omega_{1})C(\omega_{1})$
$\displaystyle
C^{(4)}=\frac{1}{6}\bar{\chi}^{2}\sum_{\omega_{1},\omega_{2}}C(\omega-\omega_{1}-\omega_{2})C(\omega_{1})C(\omega_{2})$
(28)

The equations for $Z_{\alpha}$ and $C$ can also be represented using Feynman
diagrams as shown in Figs. 3 and 4 in appendices C.3 and C.4. More explicit
forms of these equations are also reproduced in this appendix as Eqs. 45 and
46.

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 Eq 30 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 of the overall conductance
with the temperature (if $\Delta T$ is held small). This is in agreement with
previous MD simulations Sääskilahti  _et al._ (2014); Le  _et al._ (2017b).

##  IX Conclusion

To summarize, we developed a self-consistent current-conserving approximation
for anharmonic systems out of equilibrium in the high-temperature (classical)
regime. There is therefore no factors of $\hbar$ in the formalism and $\omega$
is to be interpreted as frequency only, not energy. Although the set of
derived equations for the current Eq. 13 the equation of motion Eq. 18, and
the Eq. 23 and Eq. 26 defining the correlation functions, were formally exact,
one has to develop approximations to solve the Dyson’s equation 23 and the
equation 26 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 Eq. 23 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 Eqs. 23 and 26.

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 ($\bar{\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 28.

Non-equilibrium averages were possible to calculate with the use of the
fluctuation-dissipation theorem (Eq. 19), the equations of motion (Eq. 18) and
the NFD theorem (Eq. 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}=-\bar{\Phi}_{ij}y_{j}$ in
order to extract the effective (non-equilibrium) harmonic force constants
$\bar{\Phi}$. The remainder can then be defined as the anharmonic force:
$F_{i}^{NEMD}=-\bar{\Phi}_{ij}y_{j}+a_{i}(y)$, and the present results maybe
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 Eq. 23 if
$\mathbb{G}\langle\partial\mathbb{A}/\partial Y\rangle\approx 1$, signaling
resonances, in which case the Taylor expansion in Eq. 23 is not appropriate.

Applications to nanoscale systems will appear in future publications.

##  X Acknowledgments

We would like to thank Prof. M. R. RahimiTabar for useful discussions on the
handling of correlations, the NFD theorem and a review of the manuscript, and
Profs. J. Shiomi and H. Cheraghchi for discussions at the early stages of this
work. I specifically thank Dr. V. Chiloyan for introducing me to the Langevin
thermostat method. Internal support at UVa from the Hobby Fund is also greatly
acknowledged.

This paper is dedicated to the memory of Rouzbeh Rastgarkafshgarkolaei with
whom I had several related discussions.

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

$\langle\langle
a(t)b(t)\rangle\rangle=\frac{1}{\tau}\int_{-\tau/2}^{\tau/2}\langle
a(t)b(t)\rangle\,dt;(\tau\to\infty)$

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

$\langle\langle
a(t)b(t)\rangle\rangle=\int_{-\tau/2}^{\tau/2}\int\frac{d\omega}{2\pi}\frac{d\omega^{\prime}}{2\pi}\,e^{-i(\omega+\omega^{\prime})t}\,\frac{dt}{\tau}\,\langle
A(\omega)B(\omega^{\prime})\rangle$

But

$\langle e^{-i\omega t}\rangle=\int_{-\tau/2}^{\tau/2}\,e^{-i\omega
t}\,\frac{dt}{\tau}=\frac{{\rm
sin}\,\omega\tau/2}{\omega\tau/2}\underset{\tau\to\infty}{\to}\delta_{\omega,0}(=1\,{\rm
or}\,\,0)$

On the other hand, taking $\tau\to\infty$ in the boundaries of integral of
$e^{-i\omega t}$ leads to $2\pi\delta(\omega)$ which should cancel the $\tau$
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
$\Omega=0$:

$\langle a(t)b(t)\rangle=\frac{1}{\tau}\int\frac{d\omega}{2\pi}\,\langle
A(\omega)B(-\omega)\rangle$ (29)

When calculating diagonal terms of autocorrelations, such as in
$\langle\zeta(\omega)\zeta^{\dagger}(\omega)\rangle=2\gamma k_{B}T\tau$, the
result is proportional to the integration time $\tau$. The latter cancels the
$\tau$ in the denominator coming from time averaging. Loosely speaking,
$2\pi\delta(\omega)/\tau=1\,{\rm or}\,0$. As an example, to calculate the
average (DC) current one simply needs to take the diagonal terms in
frequencies, omitting factors of $\tau$ appearing in the correlation
functions:

$\langle j_{\alpha}\rangle={\rm
lim}_{\stackrel{{\scriptstyle\tau\to\infty}}{{\Omega\to 0}}}\frac{\langle
J_{\alpha}(\Omega)\rangle}{\tau}=\Re\int\frac{d\omega}{2\pi}\frac{(-i\omega)}{\tau}\Big{[}\langle
Y(\omega)\eta_{\alpha}^{\dagger}(\omega)\rangle+\langle
Y(\omega)Y^{\dagger}(\omega)\rangle\,\sigma_{\alpha}^{\dagger}(\omega)\Big{]}$
(30)

##  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={\sqrt{m}}\,u$, about the zero temperature
equilibrium positions:

$\mathcal{V}(x_{1},...,x_{n})=\mathcal{V}(0,...0)+\Pi_{i}x_{i}+\frac{1}{2!}\Phi_{ij}\,x_{i}\,x_{j}+\frac{1}{3!}\,\Psi_{ijk}\,x_{i}x_{j}x_{k}+\frac{1}{4!}\,\chi_{ijkl}\,x_{i}x_{j}x_{k}x_{l}+...$
(31)

where, for brevity, we omitted the summation sign over repeated indices, and
where $\Pi_{i}$ is the residual force on atom $i$ and is zero if one starts
from a fully relaxed configuration. The coefficients of this expansion can be
obtained from a zero-temperature DFT calculation for instanceEsfarjani and
Stokes (2008); Tadano  _et al._ (2014).

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 $\langle x\rangle$
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 $\langle x\rangle$ is non-zero, we Taylor expand the potential
about the ”non-equilibrium” values $\langle x\rangle$ 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:

$y(t)=x(t)-\langle x\rangle=>Y(\omega)=X(\omega)-2\pi\delta(\omega)\langle
x\rangle;\langle y\rangle=0$

The average position is defined to be the solution of the average force being
zero. In terms of these new variables, and up to second-order in powers of
displacements, the force on atom $i$ can be written as:

$-\frac{\partial\mathcal{V}}{\partial x_{i}}=F_{i}=F_{i}(x=\langle
x\rangle)+\Big{(}\frac{\partial F_{i}}{\partial y_{j}}\Big{)}_{x=\langle
x\rangle}\,y_{j}+\frac{1}{2}\Big{(}\frac{\partial^{2}F_{i}}{\partial
y_{j}y_{k}}\Big{)}_{x=\langle x\rangle}\,y_{j}y_{k}+...$ (32)

The first term, the residual force, will be used to define the thermal
expansion $\langle x\rangle$ (by setting the residual or average force to
zero: $F_{i}(x=\langle x\rangle)=-(\frac{\partial\mathcal{V}}{\partial
y_{i}})_{x=\langle x\rangle}=0$), which is the Eq. 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 Eq. 31, one
can explicitly work out the equation satisfied by $\langle x\rangle$, and
derive explicit expressions for the effective harmonic and anharmonic FCs. The
expression for the force in this case is as follows:

$\displaystyle-\frac{\partial\mathcal{V}}{\partial x_{i}}$
$\displaystyle=F_{i}=-\Phi_{ij}x_{j}-\frac{1}{2}\Psi_{ijk}x_{j}x_{k}-\frac{1}{6}\chi_{ijkl}\,x_{j}x_{k}x_{l}-\Sigma_{\alpha}V_{\alpha}\,x_{\alpha}$
$\displaystyle=-\Phi_{ij}\langle x_{j}\rangle-\frac{1}{2}\Psi_{ijk}\langle
x_{j}\rangle\langle x_{k}\rangle-\frac{1}{6}\chi_{ijkl}\langle
x_{j}\rangle\langle x_{k}\rangle\langle
x_{l}\rangle-\Phi_{ij}y_{j}-\frac{1}{2}\Psi_{ijk}(y_{j}y_{k}+y_{j}\langle
x_{k}\rangle+\langle x_{j}\rangle y_{k})$
$\displaystyle-\frac{1}{6}\chi_{ijkl}\,(y_{j}y_{k}y_{l}+3y_{j}y_{k}\,\langle
x_{l}\rangle+3y_{j}\langle x_{k}\rangle\langle
x_{l}\rangle)-\Sigma_{\alpha}V_{\alpha}\,x_{\alpha}$

where use was made of the symmetry of $\Psi$ and $\chi$ under permutation of
their $j,k,l$ indices.

###  B.1 Thermal expansion

To find the average positions $\langle x\rangle$, we set the average force on
each atom $i$, $\langle
F_{i}\rangle=-\langle\big{(}\partial\mathcal{V}/\partial
x_{i}\big{)}_{x=\langle x\rangle}\rangle$ to zero, and $\langle
x_{\alpha}\rangle=0$ (leads are harmonic and will not thermal-expand), and
solve for $\langle x\rangle$. 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 $\langle x\rangle$ in terms
for force constants and averages of $\langle yy\rangle$ and $\langle
yyy\rangle$:

$\Big{[}\Phi_{ij}\,+\frac{1}{2}\chi_{ijkl}\,\langle
y_{k}y_{l}\rangle\Big{]}\langle x_{j}\rangle+\frac{1}{2}\Psi_{ijk}\,\langle
x_{j}\rangle\langle x_{k}\rangle+\frac{1}{6}\chi_{ijkl}\,\langle
x_{j}\rangle\langle x_{k}\rangle\langle
x_{l}\rangle=-\frac{1}{2}\Psi_{ijk}\,\langle
y_{j}y_{k}\rangle-\frac{1}{6}\chi_{ijkl}\,\langle y_{j}y_{k}y_{l}\rangle$ (33)

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

$\displaystyle\begin{aligned}
\bar{\Phi}_{ij}&=\Big{(}\frac{\partial^{2}\mathcal{V}}{\partial x_{i}\partial
x_{j}}\Big{)}_{\langle x\rangle}=\Phi_{ij}+\Psi_{ijk}\,\langle
x_{k}\rangle+\frac{1}{2}\chi_{ijkl}\,\langle x_{k}\rangle\langle
x_{l}\rangle\\\
\bar{\Psi}_{ijk}&=\Big{(}\frac{\partial^{3}\mathcal{V}}{\partial x_{i}\partial
x_{j}\partial x_{k}}\Big{)}_{\langle x\rangle}=\Psi_{ijk}+\chi_{ijkl}\,\langle
x_{l}\rangle\\\ \bar{\chi}_{ijkl}&=\chi_{ijkl}\end{aligned}$

Finally, the average positions can be re-expressed in terms of the new force
constants as:

$\boxed{\bar{\Phi}_{ij}\,\langle
x_{j}\rangle-\frac{1}{2}\bar{\Psi}_{ijk}\,\langle x_{j}\rangle\langle
x_{k}\rangle+\frac{1}{6}\bar{\chi}_{ijkl}\,\langle x_{j}\rangle\langle
x_{k}\rangle\langle x_{l}\rangle=-\frac{1}{2}\bar{\Psi}_{ijk}\,\langle
y_{j}y_{k}\rangle-\frac{1}{6}\bar{\chi}_{ijkl}\,\langle
y_{j}y_{k}y_{l}\rangle}$ (34)

These coupled set of non-linear equations (for all $i$) in (${\langle
x\rangle},\bar{\Phi},\bar{\Psi},\bar{\chi}$ must be solved in terms of
$\langle yy\rangle$ and $\langle yyy\rangle$ covariance matrices, to find the
average positions and effective force constants. It requires a self-consistent
iterative solution as the $\langle yy\rangle$ correlations will in turn depend
on $\langle x\rangle$.

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

$\displaystyle\frac{d^{2}y_{i}}{dt^{2}}$
$\displaystyle=-\Big{(}\frac{\partial^{2}\mathcal{V}}{\partial y_{i}\partial
y_{j}}\Big{)}_{y=0}\,y_{j}-\frac{1}{2}\Big{(}\frac{\partial^{3}\mathcal{V}}{\partial
y_{i}\partial y_{j}\partial
y_{k}}\Big{)}_{y=0}\,y_{j}y_{k}-\frac{1}{6}\Big{(}\frac{\partial^{4}\mathcal{V}}{\partial
y_{i}\partial y_{j}\partial y_{k}\partial
y_{l}}\Big{)}_{y=0}\,y_{j}y_{k}y_{l}-...-\Sigma_{\alpha}V_{\alpha}\,x_{\alpha}$
$\displaystyle=-\bar{\Phi}_{ij}\,y_{j}-\frac{1}{2}\bar{\Psi}_{ijk}\,y_{j}y_{k}-\frac{1}{6}\bar{\chi}_{ijkl}\,y_{j}y_{k}y_{l}-\Sigma_{\alpha}V_{\alpha}\,x_{\alpha}$
(35)

Adopting this change of variable, the GF keeps the same form with $\Phi$
substituted by $\bar{\Phi}$, and anharmonic forces in the right-hand side of
the equation of motion become by definition:

$a\stackrel{{\scriptstyle\rm
def}}{{=}}-\frac{1}{2}\bar{\Psi}_{ijk}\,y_{j}y_{k}-\frac{1}{6}\bar{\chi}_{ijkl}\,y_{j}y_{k}y_{l}$
(36)

##  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 Eq. 30. 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 3 dimensions, each of these labels refers to an
atom and one of the 3 cartesian components of its displacements. It therefore
varies from 1 to 3N. As a result the matrices $C$, $\sigma_{\alpha}$ and
$Z_{\alpha}$ are $3N\times 3N$ matrices. Taking their trace in formula 30
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 $\alpha$ is given by
$J_{\alpha}=\langle\dot{x}F_{x}+\dot{y}F_{y}+\dot{z}F_{z}\rangle$, and the
lead is infinite along the $z$ direction for instance, then we should have
$J_{\alpha}=\langle\dot{z}F_{z}\rangle$ and
$\langle\dot{x}F_{x}\rangle=\langle\dot{y}F_{y}\rangle=0$. So in each lead,
one may take 3 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 $\bar{\Phi}_{ij}$ and
$\bar{\Psi}_{ijk}$ are invariant under permutations of their indices.

###  C.1 Lead Self-energies $\sigma_{\alpha}$ and escape rates
$\Gamma_{\alpha}$

The surface green functions of the leads were defined by an inverse, as stated
in Eq. 10. The self-energy for lead $\alpha$ is accordingly defined as
$\sigma_{\alpha,ij}(\omega)=V_{\alpha,ik}\,\gamma_{\alpha,kk^{\prime}}(\omega)\,V_{\alpha,jk^{\prime}}$
where the two indices $(k,k^{\prime})$ refer to the lead degrees of freedom.
Note the matrices $V_{\alpha}$ connecting the device to lead $\alpha$ need not
be square.

The escape rates were defined by:

$\Gamma_{\alpha,ij}=-i(\sigma_{\alpha,ij}-\sigma_{\alpha,ij}^{\dagger})=2\Im(\sigma_{\alpha,ij})$
(37)

For Langevin thermostats with white noise, these results do not depend on the
thermostat damping factors $\gamma$.

###  C.2 Noise autocorrelation functions

Let us start with the noise autocorrelation which will appear in the
calculation of displacement-noise correlation $Z_{\alpha}$. This
autocorrelation can be obtained from the fluctuation-dissipation theorem, a
relation which has to hold if the noise is such that the leads are to behave
as a Langevin thermostat at temperature $T$:

$\langle\zeta_{\alpha}(t)\zeta_{\alpha^{\prime}}^{T}(t^{\prime})\rangle\stackrel{{\scriptstyle
FD}}{{=}}2\gamma_{\alpha}k_{B}T_{\alpha}\,\delta(t-t^{\prime})\delta_{\alpha,\alpha^{\prime}}$

The noise is white and different sites are uncorrelated with each other. This
relation implies that in the frequency domain the autocorrelation of $\zeta$
and that of $\eta$ satisfy the following relations:

$\langle\zeta_{\alpha}(\omega)\zeta_{\alpha^{\prime}}^{T}(\omega^{\prime})\rangle=2\pi\delta(\omega+\omega^{\prime})\,\delta_{\alpha,\alpha^{\prime}}\,2\gamma_{\alpha}k_{B}T_{\alpha}$
(38)

$\displaystyle\langle\eta_{\alpha}(\omega)\,\eta_{\alpha^{\prime}}^{T}(\omega^{\prime})\rangle=V_{\alpha}g_{\alpha}\
,\langle\zeta_{\alpha}\zeta_{\alpha^{\prime}}^{T}\rangle\,g_{\alpha^{\prime}}^{T}V_{\alpha^{\prime}}^{T}$
(39)
$\displaystyle=2\pi\delta(\omega+\omega^{\prime})\delta_{\alpha,\alpha^{\prime}}\,k_{B}T_{\alpha}\,[V_{\alpha}g_{\alpha}\,(2\gamma_{\alpha})\,g_{\alpha}^{T}V_{\alpha}^{T}]$
(40)

In the above, $g_{\alpha}$ is a diagonal matrix of size equal to the number of
degrees of freedom in the lead $\alpha$ (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

$i(G-G^{\dagger})=iG(\frac{1}{G^{\dagger}}-\frac{1}{G})G^{\dagger}=G\,(2b)\,G^{\dagger}$
(41)

Since $\Im(g_{\alpha}^{-1})=-\omega\gamma_{\alpha}$ we can write the
$\eta$-autocorrelation as:

$\displaystyle\langle$
$\displaystyle\eta_{\alpha}(\omega)\,\eta_{\alpha^{\prime}}^{T}(\omega^{\prime})\rangle=-iV_{\alpha}(g_{\alpha}-g_{\alpha}^{\dagger})V_{\alpha}^{T}\times
2\pi\delta(\omega+\omega^{\prime})\delta_{\alpha,\alpha^{\prime}}\,\frac{k_{B}T_{\alpha}}{\omega}$
$\displaystyle=-i(\sigma_{\alpha}-\sigma_{\alpha}^{\dagger})\times
2\pi\delta(\omega+\omega^{\prime})\delta_{\alpha,\alpha^{\prime}}\,\frac{k_{B}T_{\alpha}}{\omega}=\Gamma_{\alpha}\times
2\pi\delta(\omega+\omega^{\prime})\delta_{\alpha,\alpha^{\prime}}\,\frac{k_{B}T_{\alpha}}{\omega}$
(42)

where we used the notation
$\Gamma_{\alpha}=-i(\sigma_{\alpha}-\sigma_{\alpha}^{\dagger})=2\Im(\sigma_{\alpha})$
for twice the imaginary part of the lead $\alpha$ self-energy. Alternatively,
the diagonal elements in frequency can be written as:

$\displaystyle\langle\eta_{\alpha}(\omega)\eta_{\alpha^{\prime}}^{T}(-\omega)\rangle=\Gamma_{\alpha}\times\delta_{\alpha,\alpha^{\prime}}\,\frac{k_{B}T_{\alpha}}{\omega}\tau$
(43)

where $\tau$ is the integration time which goes to infinity. Dhar and RoyDhar
and Roy (2006) 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

$\langle\eta_{\alpha}(\omega)\eta_{\alpha^{\prime}}^{T}(\omega^{\prime})\rangle=h\Gamma_{\alpha}\delta(\omega+\omega^{\prime})\delta_{\alpha,\alpha^{\prime}}\Big{(}1+f(\hbar\omega/k_{B}T)\Big{)}$

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/\hbar\omega$. 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
$\gamma_{\alpha}$ has been replaced by twice the imaginary part of the lead
self-energy $\Gamma_{\alpha}$, which one may interpret as $2\omega$ times a
“rate”. For the adopted white noise, the dependence on its damping factor has
gone away! Furthermore, in the calculation of $Z_{\alpha}$, the notation
$\langle\eta\eta^{\dagger}\rangle$ implies the diagonal terms
$\omega^{\prime}=-\omega$ must be taken. Using the Novikov-Furutsu-Donsker
(NFD) identity (see appendix D) we can see that $\langle
A\eta^{\dagger}\rangle\propto\langle\eta\eta^{\dagger}\rangle$ 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 ($\Omega=0$) component. Note that this component is
infinite the way we defined it:
$J_{\alpha}(\Omega=0)=\int_{-\infty}^{\infty}\langle
j_{\alpha}(t)\rangle\,dt$. To find the average heat current, we have to divide
this by the integration time $\tau$, and then take the limit $\tau\to\infty$.
This division cancels the $\tau$ appearing in the numerator of noise
autocorrelations. This issue is further discussed in the appendix, Eq. 30.
Final results do not depend on $\tau$.

The factor $2\pi\delta(\omega+\omega^{\prime})/\tau$ can now be excluded from
the current as it is essentially 1 for the diagonal terms, and zero otherwise.

$\boxed{\langle\eta_{\alpha}(\omega)\eta_{\alpha}^{\dagger}(\omega)\rangle=\Gamma_{\alpha}(\omega)\frac{k_{B}T_{\alpha}}{\omega}\tau=2\Im(\sigma_{\alpha}(\omega))\frac{k_{B}T_{\alpha}}{\omega}\tau}$
(44)

###  C.3 Noise-displacement correlations $Z_{\alpha}$

Using the equation of motion 18, the expression defining this correlation
function to second order in $\partial\mathbb{A}/\partial Y$ can be derived to
satisfy Eq. 24. After substitution, the explicit relations become:

$\displaystyle
Z_{\alpha,ij}(\omega)=Z^{NEMF}_{\alpha,ij}(\omega)+\mathbb{G}_{ik}(\omega)\,\Big{(}\Sigma^{(3)}_{kk^{\prime}}+\Sigma^{(4)}_{kk^{\prime}}\Big{)}\,Z_{\alpha,k^{\prime}j}(\omega)$
$\displaystyle
Z^{NEMF}_{\alpha,ij}(\omega)=\mathbb{G}_{ik}(\omega)\,\Gamma_{\alpha,kj}(\omega)\frac{k_{B}T_{\alpha}}{\omega}$
$\displaystyle\Sigma^{(3)}_{kk^{\prime}}(\omega)=2\bar{\Psi}_{klm}\,\bar{\Psi}_{k^{\prime}l^{\prime}m^{\prime}}\Big{(}\int_{1}\mathbb{G}_{ll^{\prime}}(\omega-\omega_{1})\,C_{mm^{\prime}}(\omega_{1})\Big{)}$
$\displaystyle\Sigma^{(4)}_{kk^{\prime}}(\omega)=\frac{3}{2}\bar{\chi}_{klmn}\,\bar{\chi}_{k^{\prime}l^{\prime}m^{\prime}n^{\prime}}\Big{(}\int_{1,2}\mathbb{G}_{ll^{\prime}}(\omega-\omega_{1}-\omega_{2})\,C_{mm^{\prime}}(\omega_{1})\,C_{nn^{\prime}}(\omega_{2})\Big{)}$
(45)

where, as usual, the summation over repeated indices is implied.

Figure 3: Feynman diagrams associated with $Z_{\alpha}$ up to second order in
anharmonic forces. Dashed lines represent the phonon Green’s function
$\mathbb{G}$, the circle represents $\Gamma_{\alpha}$, the triangle represents
the third-order vertex $\bar{\Psi}$ and the square the fourth-order vertex
$\chi$. The thick line with opposite arrows represents the displacement
autocorrelation $C(\omega)$. Frequency must be conserved at each vertex.

###  C.4 Displacement autocorrelations $C(\omega)$

These functions, whose trace represents the number of excitations in the
device, are defined via Eq 26. The static distortion in $\langle
x\rangle\langle x\rangle$ does not contribute to the current and will thus be
omitted. To second-order in $\mathbb{A}$, the expression for $C$ is given by:

$\displaystyle
C_{ij}(\omega)=\sum_{\alpha}(Z_{\alpha,ik}-Z^{NEMF}_{\alpha,ik})\mathbb{G}_{kj}^{\dagger}+\mathbb{G}_{ik}(Z_{\alpha,kj}-Z^{NEMF}_{\alpha,kj})^{\dagger}-C_{ij}^{NEMF}(\omega)+\mathbb{G}_{ik}\,\Big{(}C^{(3)}_{kk^{\prime}}+C^{(4)}_{kk^{\prime}}\Big{)}\,\mathbb{G}_{k^{\prime}j}^{\dagger}$
$\displaystyle
C_{ij}^{NEMF}(\omega)=\Sigma_{\alpha}\,(\mathbb{G}_{ik}\,\Gamma_{\alpha,kk^{\prime}}\,\mathbb{G}^{\dagger}_{k^{\prime}j})\,\frac{k_{B}T_{\alpha}}{\omega}$
$\displaystyle
C^{(3)}_{kk^{\prime}}(\omega)=\frac{1}{2}\bar{\Psi}_{klm}\,\bar{\Psi}_{k^{\prime}l^{\prime}m^{\prime}}\Big{(}\int_{1}C_{ll^{\prime}}(\omega-\omega_{1})\,C_{mm^{\prime}}(\omega_{1})\Big{)}$
$\displaystyle
C^{(4)}_{kk^{\prime}}(\omega)=\frac{1}{6}\bar{\chi}_{klmn}\,\bar{\chi}_{k^{\prime}l^{\prime}m^{\prime}n^{\prime}}\Big{(}\int_{1,2}C_{ll^{\prime}}(\omega-\omega_{1}-\omega_{2})\,C_{mm^{\prime}}(\omega_{1})\,C_{nn^{\prime}}(\omega_{2})\Big{)}$
(46)

Figure 4: Feynman diagrams associated with $C(\omega)$ up to second order in
anharmonic forces. Conventions are the same as in Fig. 3

##  Appendix D Statement and proof of the Novikov-Furutsu-Donsker (NFD)
relation

The Novikov-Furutsu-Donsker relationDonsker  _et al._ (1962); Furutsu (1964);
Novikov (1965); Rahimi Tabar (2019) relates the correlation function of any
functional $A[\eta]$ with the noise $\eta$ to the noise-noise correlation
times the expectation value of the derivative of that functional:

$\langle A[\eta]\,\eta_{\alpha}^{\dagger}(t)\rangle=\sum_{\beta}\int
dt^{\prime}\langle\eta_{\beta}(t^{\prime})\eta_{\alpha}^{\dagger}(t)\rangle\langle\frac{\delta
A[\eta]}{\delta\eta_{\beta}(t^{\prime})}\rangle$ (47)

implying that in the frequency domain, we have

$\langle
A_{i}[\eta]\,\eta_{j\alpha}^{\dagger}(\omega)\rangle=\sum_{l,\beta}\langle\frac{\delta
A_{i}[\eta]}{\delta\eta_{l,\beta}}(\omega)\rangle\langle\eta_{l,\beta}(\omega)\eta_{i,\alpha}^{\dagger}(\omega)\rangle$
(48)

Proof of the NFD relation: Let $\eta$ be a random variable with Gaussian
distribution of mean 0 and variance $\sigma^{2}$. Consider the average
$S=\langle A[\eta]\eta\rangle$ where $A$ is a functional of $\eta$. By
definition,
$S=\int_{-\infty}^{+\infty}d\eta\frac{e^{-\eta^{2}/2\sigma^{2}}}{\sqrt{2\pi\sigma^{2}}}\eta
A(\eta)$. After performing an integration by parts, we find:

$S=\int_{-\infty}^{+\infty}-\sigma^{2}\frac{1}{\sqrt{2\pi\sigma^{2}}}d(e^{-\eta^{2}/2\sigma^{2}})A(\eta)=\sigma^{2}\int_{-\infty}^{+\infty}\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\eta^{2}/2\sigma^{2}}\frac{dA(\eta)}{d\eta}d\eta=\sigma^{2}\langle\frac{dA(\eta)}{d\eta}\rangle=\langle\eta\eta\rangle\langle\frac{dA(\eta)}{d\eta}\rangle$

This can easily be extended to higher dimensions where $\eta$ is an array:
$\langle A[\eta]\eta_{i}\rangle=\sum_{j}\langle\frac{\partial
A}{\partial\eta_{j}}\rangle\langle\eta_{j}\eta_{i}\rangle$.

##  Appendix E Heat current within the Harmonic approximation

For a harmonic system ($G=\mathcal{G}=\mathbb{G}$), we have derived the
following expressions :

$\displaystyle Z_{\alpha}^{H}(\omega)$ $\displaystyle=\langle
Y\eta^{\dagger}_{\alpha}\rangle=\mathcal{G}\Gamma_{\alpha}\frac{k_{B}T_{\alpha}}{\omega}\tau;$
(49) $\displaystyle C^{H}(\omega)$ $\displaystyle=\langle
YY^{\dagger}\rangle=\Sigma_{\beta}(\mathcal{G}\Gamma_{\beta}\mathcal{G}^{\dagger})\frac{k_{B}T_{\beta}}{\omega}\tau$
(50)

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:

$\langle j_{\alpha}^{H}\rangle=k_{B}\int\frac{d\omega}{2\pi}\,{\rm
Tr}\,\Big{(}\Im(\mathcal{G})\,\Gamma_{\alpha}T_{\alpha}-\Sigma_{\beta}(\mathcal{G}\,\Gamma_{\beta}T_{\beta}\,\mathcal{G}^{\dagger})\frac{\Gamma_{\alpha}}{2}\,\Big{)}$

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 Eq. 41, we
have

${\rm Tr}\,\Gamma_{\alpha}\,\Im(\mathcal{G})=\frac{1}{2}{\rm
Tr}\,\Gamma_{\alpha}\Sigma_{\beta}(\mathcal{G}\Gamma_{\beta}\mathcal{G}^{\dagger})$
(51)

Making this substitution for $\Im(\mathcal{G})$ in the heat current formula,
we end up with:

$\langle j_{\alpha}^{H}\rangle=\frac{k_{B}}{2}\int\frac{d\omega}{2\pi}\,{\rm
Tr}\,\Gamma_{\alpha}\,\mathcal{G}\,\Big{(}\Sigma_{\beta}\Gamma_{\beta}(T_{\alpha}-T_{\beta})\Big{)}\mathcal{G}^{\dagger}$
(52)

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 $\alpha$ and $\beta$) form of the following sum:

$\Sigma_{\alpha}\langle
j_{\alpha}^{H}\rangle=\frac{k_{B}}{2}\int\frac{d\omega}{2\pi}\,\Sigma_{\alpha\beta}\,(T_{\alpha}-T_{\beta})\,{\rm
Tr}\,\Big{[}\Gamma_{\alpha}\,\mathcal{G}\,\Gamma_{\beta}\mathcal{G}^{\dagger}\Big{]}=0$
(53)

In a two-terminal device geometry ($\beta=R;\alpha=L$), the expression for the
current reduces to the well-known formula:

$\langle
j_{L}^{H}\rangle=\frac{k_{B}}{2}(T_{L}-T_{R})\int\frac{d\omega}{2\pi}\,{\rm
Tr}\,[\Gamma_{L}\,\mathcal{G}\,\Gamma_{R}\,\mathcal{G}^{\dagger}]=-\langle
j_{R}^{H}\rangle$ (54)

In general, ${\rm
Tr}\,[\Gamma_{\alpha}\,\mathcal{G}\,\Gamma_{\beta}\,\mathcal{G}^{\dagger}]$
maybe interpreted as the harmonic transmission from lead $\alpha$ to lead
$\beta$. We see that even in the non-equilibrium regime (large $\Delta 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/\omega$ by $\hbar(1+f_{\rm BE}(\hbar\omega/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_{\alpha}-T_{\beta})<=>\hbar\omega(f_{\alpha}-f_{\beta})$.

## References

  * Keldysh (1965) L. V. Keldysh, Soviet Physics JETP 20, 1018 (1965).
  * Caroli  _et al._ (1971) C. Caroli, R. Combescot, D. Lederer, P. Nozieres, and D. Saint-James, Journal of Physics C: Solid State Physics 4, 2598 (1971).
  * Meir and Wingreen (1992) Y. Meir and N. S. Wingreen, Phys. Rev. Lett. 68, 2512 (1992).
  * Pastawski (1991) H. M. Pastawski, Phys. Rev. B 44, 6329 (1991).
  * Mingo (2006) N. Mingo, Physical Review B 74, 125402 (2006).
  * Wang  _et al._ (2006) J. Wang, J. Wang, and N. Zeng, Phys. Rev. B 74, 033408 (2006).
  * Dai and Tian (2020) J. Dai and Z. Tian, Phys. Rev. B 101, 041301 (2020).
  * Guo  _et al._ (2020) Y. Guo, M. Bescond, Z. Zhang, M. Luisier, M. Nomura, and S. Volz, Physical Review B 102, 195412 (2020), 2007.14046 .
  * Luisier (2012) M. Luisier, Phys. Rev. B 86, 245407 (2012).
  * Polanco (2021) C. A. Polanco, Nanoscale and Microscale Thermophysical Engineering 25, 1 (2021), https://doi.org/10.1080/15567265.2021.1881193 .
  * Sääskilahti  _et al._ (2013) K. Sääskilahti, J. Oksanen, and J. Tulkki, Phys. Rev. E 88, 012128 (2013).
  * Bolsterli  _et al._ (1970) M. Bolsterli, M. Rich, and W. M. Visscher, Phys. Rev. A 1, 1086 (1970).
  * Le  _et al._ (2017a) N. Q. Le, C. A. Polanco, R. Rastgarkafshgarkolaei, J. Zhang, A. W. Ghosh, and P. M. Norris, Phys. Rev. B 95, 245417 (2017a).
  * Feng  _et al._ (2019) T. Feng, Y. Zhong, J. Shi, and X. Ruan, Phys. Rev. B 99, 045301 (2019).
  * Sääskilahti  _et al._ (2014) K. Sääskilahti, J. Oksanen, J. Tulkki, and S. Volz, Phys. Rev. B 90, 134312 (2014).
  * Chalopin  _et al._ (2012) Y. Chalopin, K. Esfarjani, A. Henry, S. Volz, and G. Chen, Phys. Rev. B 85, 195302 (2012).
  * Chalopin and Volz (2013) Y. Chalopin and S. Volz, Applied Physics Letters 103, 051602 (2013), https://doi.org/10.1063/1.4816738 .
  * Tian  _et al._ (2014) Z. Tian, K. Esfarjani, and G. Chen, Physical Review B - Condensed Matter and Materials Physics 89, 235307 (2014).
  * (19) With a factor of $\hbar$ in the denominator $f_{\alpha}$ would be a real Bose-Einstein occupation factor taken to the classical limit .
  * Rahimi Tabar (2019) M. Rahimi Tabar,  _Analysis and data-based reconstruction of complex nonlinear dynamical systems_ (Springer, 2019).
  * Le  _et al._ (2017b) N. Q. Le, C. A. Polanco, R. Rastgarkafshgarkolaei, J. Zhang, A. W. Ghosh, and P. M. Norris, Phys. Rev. B 95, 245417 (2017b).
  * Esfarjani and Stokes (2008) K. Esfarjani and H. T. Stokes, Phys. Rev. B 77, 144112 (2008).
  * Tadano  _et al._ (2014) T. Tadano, Y. Gohda, and S. Tsuneyuki, Journal of Physics Condensed Matter 26, 225402 (2014).
  * Dhar and Roy (2006) A. Dhar and D. Roy, Journal of Statistical Physics 125, 805 (2006), arXiv:0606465 [cond-mat] .
  * Donsker  _et al._ (1962) M. Donsker, J. Lions, _et al._ , Acta Mathematica 108, 147 (1962).
  * Furutsu (1964) K. Furutsu,  _On the Statistical Theory of Electromagnetic Waves in a Fluctuating Medium (II) Mathematical Basis of the Analogies to Quantum Field Theory_ (Government Press, 1964).
  * Novikov (1965) E. A. Novikov, Sov. Phys. JETP 20, 1290 (1965).