A Derivation of Eq. (26)

# Classical heat transport in anharmonic molecular junctions: exact solutions

## Abstract

We study full counting statistics for classical heat transport through anharmonic/nonlinear molecular junctions formed by interacting oscillators. Analytical result of the steady state heat flux for an overdamped anharmonic junction with arbitrary temperature bias is obtained. It is found that the thermal conductance can be expressed in terms of temperature dependent effective force constant. The role of anharmonicity is identified. We also give the general formula for the second cumulant of heat in steady state, as well as the average geometric heat flux when two system parameters are modulated adiabatically. We present an anharmonic example for which all cumulants for heat can be obtained exactly. For a bounded single oscillator model with mass we found that the cumulants are independent of the nonlinear potential.

###### pacs:
05.60.cd, 05.70.Ln, 44.10.+i

## I Introduction

Anharmonicity or nonlinearity plays an important role in many physical processes. For example, it is crucial for umklapp scattering which gives rise to a finite thermal conductivity for bulk systems as pointed out by Peierls Peierls (1955). It is also an essential ingredient required for developing functional phononic devices Li et al. (2012), such as thermal rectifier Terraneo et al. (2002); Chang et al. (2006) and heat pump Segal and Nitzan (2006); Segal (2008); Ren et al. (2010).

Comprehensive understanding of the general behaviors of nonlinearity in transport is a fundamental problem in condensed matter physics and nonequilibrium statistical physics. Although elegant approaches exist to deal with ballistic transport in harmonic systems Dhar and Roy (2006); Wang et al. (2008), solving an anharmonic model analytically is extremely challenging due to the intrinsic non-integrability of equations of motion. As a result, for classical heat transport, the major approaches are limited to molecular dynamics simulations Lepri et al. (2003); Dhar (2008), perturbation theories such as Boltzmann equation Peierls (1955); Pereverzev (2003), mode coupling theory Lepri et al. (1998); Delfini et al. (2007), and/or mean–field theory Li et al. (2006); He et al. (2008). There does not exist analytic solutions of heat transport problem for general anharmonic potentials.

Moreover, in heat transport problems, not only the heat flux, but also the higher order fluctuations (noises) are of significant interest as they contain valuable information about the transport properties and satisfy universal symmetry known as Fluctuation Theorems (FT) Jarzynski (1997); Crooks (1999); Gallavotti and Cohen (1995). Therefore, methods to obtain these higher order moments and to verify this symmetry for arbitrary anharmonic potentials are also highly desirable.

In this article, we present a general approach to nonlinear heat transport by working in an enlarged space, which includes the heat as an extra variable in the Fokker-Planck (FP) equation Ren et al. (2012). The heat current can be expressed as the expectation value of certain operator. If the steady state of the problem is known, then the analytic expression for the current can be obtained. We show this is the case for two different anharmonic oscillator systems, i.e., the overdamped interacting double-oscillator model (Model I) and the bounded single-oscillator model (Model II).

For Model I, the complete cumulants, as well as the geometric contribution to the heat flux, can be obtained if the full spectrum of the original FP equation is solvable. One such example is the V-shaped potential . For Model II, the complete cumulants can always be analytically obtained no matter what the nonlinear interacting potential is. For both models, the Gallavotti-Cohen (GC) symmetry can be verified for arbitrary anharmonic interactions.

In the rest of the article, Sec. II is devoted to Model I and Sec. III is devoted to Model II. After that, we briefly summarize the article in Sec. IV.

## Ii Model I

Model I consists of two coupled Brownian oscillators which are connected to two Langevin heat baths (Fig. 1(a)). The dynamics is described by the Langevin equations: , (), where and are the mass and displacement of oscillator . is the coupling strength between the oscillator and the bath , is the temperature of bath and is a white noise with variance ( is set to 1). depicts the interaction potential, which is assumed to be bounded and depend only on the distance between the two oscillators.

This model is a minimal model for heat transport in interacting oscillator systems. At low Reynolds number, we can neglect the inertia and thus obtain the overdamped dynamics Purcell (1977):

 γ1˙x1+V′(x1−x2)=ξ1; (1) γ2˙x2−V′(x1−x2)=ξ2, (2)

in which the prime denotes the derivative with respect to its argument. The total heat transferred from bath 1 to bath 2 up to time , denoted as , satisfy

 ˙Q=V′(x1−x2)˙x1. (3)

By introducing the relative coordinate , we obtain two stochastic differential equations for the random variables and :

 ˙y =−V′(y)γ1−V′(y)γ2+ξ1γ1−ξ2γ2; (4) ˙Q =V′(y)γ1(ξ1−V′(y)). (5)

Their joint probability distribution then evolves according to the FP equation Risken (1998); Kundu et al. (2011)

 ∂ρ∂t=[Tγ∂2∂y2+1γ∂∂yV′+T1γ1∂2∂Q2V′2+∂∂Q(V′2γ1−T1V′′γ1)+2T1γ1∂2∂Q∂yV′]ρ, (6)

where and .

By taking the Fourier transformation on in Eq. (6), we obtain the equation for the characteristic function (CF) , where is short for throughout the paper. Moreover, we introduce and use tilde ( ) to express the -dependence (in such a notation, is omitted in the parameter list). Then the equation reads

 ∂t~z(y,t)=~L~z(y,t) (7)

with twisted FP operator

 ~L=Tγ∂2∂y2+1γ∂∂yV′+λ(T1V′′γ1−V′2γ1−2T1γ1∂∂yV′)+λ2T1V′2γ1≡L(0)+λL(1)+λ2L(2), (8)

For any quantity or operator , We shall also use to be short for evaluated at . For example, and is the probability distribution for . Moreover, for , the equation actually depicts a single Brownian oscillator connected with a single heat bath at an effective temperature and subjected to a potential (see Fig. 1(b)). Therefore, the stationary solution will be in the form of Boltzmann distribution with as the partition function.

We are only interested in the heat . Therefore, by integrating over , we obtain the CF of as , with . The cumulant generating function (CGF) reads , which generates the th-order cumulant as .

Obviously, if we are able to solve Eq. (7) for general potential , then the problem can be completely solved, which is the case for harmonic potential Ren et al. (2012). For anharmonic cases, a general method is not available.

However, using separation of variables , Eq. (7) can still be casted into an eigenvalue problem

 ~F(t)=e−~μt,~L~φ(y)=−~μ~φ(y). (9)

with natural boundary conditions for . The structure of the solutions is clearer if we introduce a transformation with . By direct substitution, it turns out that satisfy the following time independent Schrödinger equation

 ~LS~Ψ≡1γ(−T∂2∂y2+~VS)~Ψ=~μ~Ψ, (10)

which describes a single particle moving in the effective potential (see Fig. 1(c)), given by

 ~VS=V′24T[1+4T1T2γγ1+γ2λ(1T1−1T2−λ)]−12V′′. (11)

Then the well-known results in quantum mechanics for single particle can be directly evoked to discuss our situation. Two of them are that all the eigenvalues are real (regard as a real quantity) which can be ordered as , and upon normalization, all the eigenfunctions form a complete orthonormal basis set, which satisfy and We use to denote the inner product .

The original operator is not Hermitian, i.e., . However, it can be shown that satisfy , namely, is an eigenfunction of with eigenvalue . When normalized properly, we will have

 (~ϕn,~φm)=δnm,∞∑n=0~ϕn(y)~φn(y′)=δ(y−y′). (12)

Then we can write down the formal solution of Eq. (7) as

 ~z(y,t)=∞∑n=0~cne−~μnt~φn(y), (13)

where is time independent and can be determined using the initial condition, .

### ii.1 Dynamic Heat Fluxes

Based on the formalism above, we come to discuss the heat fluxes in nonequilibrium steady state. It should be pointed out that the discussions here are also applicable to the dynamic heat fluxes Ren et al. (2010, 2012) when system parameters are modulated adiabatically.

If the long-time behavior is of main interest, then the component with smallest eigenvalue () will dominate, so that and . Consequently, the CGF and the th cumulant of the heat fluctuations are given by

 Missing or unrecognized delimiter for \right (14)

Of course, one would not be able to get all the cumulants unless can be obtained as a function of . However, if one is only interested in the first cumulant, i.e., the average heat flux , then here is an approach.

We differentiate Eq. (9) with respect to and obtain

 (∂λ~L)~φ0+~L∂λ~φ0=−(∂λ~μ0)~φ0−~μ0∂λ~φ0. (15)

Then we perform on both sides. Noticing that , we obtain

 ¯J=−∂λ~μ0|λ=0=(~ϕ0,(∂λ~L)~φ0)∣∣λ=0=(ϕ0,L(1)φ0). (16)

This equation is the analogue of the Hellmann-Feynman theorem in quantum mechanics. Note that and in Eq. (16) are evaluated at , therefore, are associated with the FP operator . It can be shown that the smallest eigenvalue for is and the corresponding eigenfunctions are and respectively Risken (1998).

After normalization, Eq. (16) becomes . In this sense, serves as the heat flux operator, which gives the average heat flux when acted on the stationary state. Using , we finally obtain

 Missing or unrecognized delimiter for \left (17)

where is the canonical ensemble average of at temperature and we have used . Remember that .

Equation (17) is exact for a general potential . It is also applicable to the cases with arbitrarily large temperature differences, i.e., the cases in a far-from-equilibrium state. As far as we know, this is the first exact result for anharmonic heat transport with general interaction, although it is only for a small system with overdamped dynamics.

For harmonic potential , , which is just the force constant and is independent of the temperature. For anharmonic potentials, will generally be temperature dependent. In analogy with the harmonic case, we define as the effective force constant.

We notice that the temperature dependence of the effective force constant plays a key role in many of the phenomena for heat control and management. For example, a thermal rectifier is a device that directs the flow of heat Terraneo et al. (2002); Chang et al. (2006). To use this oscillator system as a thermal rectifier, it is required that

 ∣∣∣¯J(T1=TH,T2=TL)¯J(T1=TL,T2=TH)∣∣∣=⟨V′′⟩T=γ2TH+γ1TLγ1+γ2⟨V′′⟩T=γ2TL+γ1THγ1+γ2≠1. (18)

Obviously, this condition is satisfied only when is temperature dependent and the system is asymmetric, namely, .

As another example, in Li et al. (2009), it was demonstrated that heat can be shuttled across a symmetric system when the temperature of one heat bath is modulated adiabatically but with zero average temperature bias, e.g., with . If our system is modulated in such a way with and being a periodic function of , one can obtain an average heat flux in one period. Again, is nonzero only if depends on . Actually, when , we can expand into Taylor series and obtain the leading term as , where denotes derivative and denotes time average over one period. This pump effect is more significant if is larger; and is vanishing when .

Since the effective force constant is independent of for harmonic potential, the phenomena mentioned above cannot be observed using harmonic potential. On the other hand, is harmonic potential the only potential that has such a property? The answer is no. An example is the Toda potential for which . Interestingly, both the harmonic lattice and Toda lattice are models displaying ballistic heat transport. These facts imply a clue that temperature independent is related to ballistic transport.

It is worth mentioning that there are other approximated mean–field theories trying to convert an anharmonic potential in lattices to effective harmonic potential with temperature dependent force constant. For example, for the FPU- model with , the effective phonon theory (EPT) estimates it as Li et al. (2006, 2010) while the self-consistent phonon theory (SCPT) estimates it as He et al. (2008). In Fig. 2, we compare these results with the present exact result for the two-oscillator model, i.e., . It is shown that both the EPT and SCPT underestimate the nonlinear effect, although at low temperature the deviations are tiny. At the high temperature regime, all these theories result in the dependence. For , the coefficient is for the exact value. For the EPT and for the SCPT , which are 27% and 14% less than the exact value, respectively.

Upto now we only discussed the first cumulant, i.e., the average heat flux. However, using the standard perturbation theory, all higher order cumulants can be derived. For example, the second cumulant reads (see for example Szabo and Ostlund (1996))

 Missing or unrecognized delimiter for \right (19)

where for short. Note that all the terms in the formula are again evaluated at .

The excited states are involved in the expression for higher order cumulants. Therefore, to analytically obtain the higher order cumulants, we should be able to completely solve all the eigenstates of the Schrödinger equation Eq. (10) with , which describes a single particle with mass moving in the potential (see Fig. 1(c)). This is of course an impossible mission for most of the anharmonic potentials. However, numerically solving it by diagonalizing the operator is straightforward. In this sense, our method provides a new numerical approach other than molecular dynamics to investigate the nonlinear heat transport problems.

It can also be noticed that, even we are not able to obtain the CGF as a function of for general , we can still verify the GC symmetry

 ~G(λ)=~G(Δβ−λ), (20)

where . This result can be seen from Eqs. (10) and (11). Because possess this symmetry, all its eigenvalues and eigenfunctions should possess this symmetry. However, and do not follow it because the transformation contains which is not invariant and hence GC symmetry is not satisfied in the transient stage. It is only valid in the long-time limit when the largest eigenvalue dominates.

### ii.2 Geometric Heat Fluxes

It has been discovered that in classic heat transport, the geometric curvature might contribute to heat transport when two or more system parameters are modulated adiabatically Ren et al. (2012).

When two parameters and are modulated, the average geometric heat flux reads

 ¯Jgeom=−1Tp∬u1u2∂~Fu1u2∂λ∣∣ ∣∣λ=0du1du2, (21)

where is the geometric curvature in the space

 ~Fu1u2=(∂~ϕ0∂u1,∂~φ0∂u2)−(∂~ϕ0∂u2,∂~φ0∂u1). (22)

In the expression, and should satisfy the normalization condition Eq. (12). However, we are free of using different conventions to normalize these eigenfunctions. It can be shown that the curvature shown above will not be affected by how we normalize the eigenfunctions, which is the analogue to “gauge invariance” in quantum mechanics when dealing with Berry curvature.

Typically, in quantum mechanics, to normalize an unnormalized eigenfunction , we need to choose . For the Fokker-Planck equation discussed here, we can still use this method and then transform to . But this will make the later derivation difficult. A more convenient way in our situation is to choose

 φ0(y)=PstT=1ZTe−V(y)Tandϕ0(y)=1. (23)

In other words, is the Boltzmann distribution and is a constant. We call this the Boltzmann convention. In this convention, it can be derived that

 ~F∘u1u2=(∂~ϕ∘0∂u1,∂φ0∂u2)−(∂~ϕ∘0∂u2,∂φ0∂u1), (24)

in which for any quantity , we use the notation

 ~A∘≡∂~A∂λ∣∣∣λ=0. (25)

#### Modulating T1 and T2

It has been discovered that when the two bath temperatures are modulated, there is no geometric contribution to the average heat flux for the harmonic potential. However, for the anharmonic FPU- model, this effect has been observed numerically Ren et al. (2012).

To speculate the role of anharmonic potential on geometric heat flux when and are modulated, we derive the following equation for the curvature (see Appendix A)

 ~FT1T2∘=∂~FT1T2∂λ∣∣λ=0=1T2(γ1+γ2)∞∑m=1V′′0mVm0μm−μ0, (26)

where we still use for any . However, since both and are functions of instead of operators, here we can also evaluate them as .

Now we are able to distinguish the effect of anharmonicity to the geometric heat flux. For harmonic potential , , which is obviously zero for . Therefore, at any point in the parameter space, and finally . For anharmonic potentials, this term would generally be nonzero, which results in a non-vanishing curvature and geometric heat flux.

Similar to Eq. (19), it can be seen that only when we can obtain all the excited states, are we able to get the final result of Eq. (26). However, for weak nonlinear cases, Eq. (19) and (26) can be used to obtain the approximated values. Consider the case , with a small . By regarding as the perturbation, we can evaluate Eq. (26) to the first order of as

 ~FT1T2∘≈α(Ψ(0)0,U′′Ψ(0)2)(Ψ(0)2,12ky2Ψ(0)0)μ(0)2(γ1+γ2)T2, (27)

where in the superscript indicates the unperturbed correspondences, which are for the harmonic potential. For the FPU- model with , Eq. (27) becomes , which is independent of the temperature. Therefore, the heat transfered in one period should be proportional to the area enclosed by the driving path.

#### Modulating γ1 and γ2

If we fix and module , then we cannot obtain a geometric heat flux, for any . This is because in this case , which is independent of and . As a consequence,

 φ0=e−V(y)/T∫e−V(y)/Tdy (28)

is also independent of and . Inserting this result into Eq. (24) immediately shows and therefore .

The same result has also been found in Ren et al. (2010) for quantum heat transport and in Ren et al. (2012) for classic case with harmonic potential.

### ii.3 An Example

Next, we discuss an anharmonic model that is analytically solvable, that is Model I with a V-shaped potential . Correspondingly, , which is an infinitely deep well potential plus a constant shift. This shift will only lift the eigenvalues and will not affect the eigenstates. Therefore, the full spectrum of the Schrödinger equation with , as well as all the cumulants, should be solvable.

Actually, for this potential, we can directly solve Eq. (10). It is easy to show that only the ground state is bounded and the corresponding eigenvalue is

 ~μ0=−k2T1T2T(γ1+γ2)[λ−(1T1−1T2)]λ, (29)

which is a quadratic form. Therefore, in the long-time limit, will be a Gaussian random variable, which only has the lowest two cumulants: and

To obtain the geometric heat flux, the ground state eigenfunctions are needed, which read

 ~ϕ0 =e−k|y|λT1/γ1T1/γ1+T2/γ2; (30) ~φ0 =k2Te−k|y|1/γ1+1/γ2−λT1/γ1T1/γ1+T2/γ2. (31)

It is then easy to obtain Interestingly, this curvature is proportional to . Therefore, the first cumulant is nonzero and all higher order cumulants are zero.

## Iii Model II

The idea presented above is not restricted to the overdamped two-oscillator model, but generally applicable to all heat transport problems. If the steady state can be solved, then at least the average heat flux can be obtained. This is the situation also for Model II, a bounded single-oscillator model with inertia, following the Langevin dynamics:

 ˙y =u; (32) m˙u =−V′(y)−γ1u−γ2u+ξ1+ξ2. (33)

Meanwhile, the heat satisfies

 ˙Q=(−γ1˙y+ξ1)˙y (34)

and their joint probability satisfies the Fokker-Planck equation

 ∂∂t~P(y,u)=~L~P(y,u)=(L(0)+λL(1)+λ2L(2))~P(y,u)=[−∂∂yu+∂∂u(γmu+V′(y)m)+∂2∂u2γTm2−λ(γ1u2−γ1T1m+2∂∂uγ1T1mu)+λ2γ1T1u2]~P(y,u), (35)

where and .

For , Eq. (35) is a single particle Kramers equation. It has a unique steady state (ground state), which is just the Boltzmann distribution Risken (1998)

 PstT(y,u)=1Ze−H/TwithH=12mu2+V(y), (36)

and being the partition function.

Therefore, we can directly use our formalism to obtain the dynamic heat flux as

 ¯J=∬L(1)PstTdydu=γ1γ2m(γ1+γ2)(T1−T2), (37)

which is independent of because is independent of . The immediate question that can be raised is whether the higher cumulants are also independent of .

In fact, for this model, the ground state of Eq. (35), i.e.,

 ~L~φ0(y,u)=(L(0)+λL(1)+λ2L(2))~φ0(y,u)=−~μ0~φ0(y,u), (38)

as well as the left eigen problem

 ~L†~ϕ0(y,u)=−~μ0~ϕ0(y,u) (39)

can be solved using trial solutions

 ~ϕ0(y,u)∝e−H/~TLand~φ0(y,u)∝e−H/~TR, (40)

in which is independent of and .

For the right eigenfunction and eigenvalue, we get

 ~L~φ0=[g(~TR)+h(~TR)u2]~φ0, (41)

where

 g(~TR) =γ(~TR−T)−γ1λT1~TRm~TR; (42) h(~TR) =γ1λ(λT1−1)~T2R+(2γ1λT1−γ)~TR+γT~T2R. (43)

To make it an eigenvalue problem, we demand . Solving it for gives

 ~TR=γ−2γ1λT1±√(γ−2γ1λT1)2−4γγ1λT(λT1−1)2γ1λ(λT1−1). (44)

There are two possible solutions. However, we know that the ground state should converge to the Boltzmann distribution Eq. (36). Therefore, it is required that

 limλ→0~TR=T=γ1T1+γ2T2γ1+γ2. (45)

So we need choose the solution with minus “-” sign. Consequently, we obtain the corresponding eigenvalue

 Missing or unrecognized delimiter for \left (46)

which is indeed independent of . Therefore, all the higher order cumulants are also independent of at long time. Previously, this formula has been obtained for the harmonic case Kundu et al. (2011); Fogedby and Imparato (2011) and numerically tested for the FPU- potential Fogedby and Imparato (2011).

Similarly, the left ground state eigenfunction can be obtained, which is in the form of Eq. (40) with

 ~TL=2γ1λT1−γ−√(γ−2γ1λT1)2−4γγ1λT(λT1−1)2γ1λ(λT1−1), (47)

since the left eigenfunction for should be 1 ().

With the ground state eigenvalue and eigenfunctions known, we can easily calculate other quantities as well, such as the geometric curvature, which is left to the readers, if interested.

## Iv Summary

Using the Fokker-Planck equation for the full counting statistics for heat, we have developed a formalism to obtain the dynamic and geometric heat current for arbitrary anharmonic potential. The importance of anharmonicity in manipulating heat is investigated and it is found that the temperature dependent force constant plays a key role. The higher cumulants can also be obtained if the eigenspectrum of the original FP operator is known. As an example, we choose shaped potential and obtained all the cumulants explicitly. We also extend this method by including the inertia term for a bounded single particle and found that the cumulants are independent of the nonlinear onsite potential. For arbitrary nonlinear potential, the Gallavotti-Cohen symmetry is verified for both models.

In common sense, the nonlinear transport problems are generally thought to be not exactly solvable. We hope the existance of exact solutions, presented here, for nonlinear transport problems will motivate further investigations on solving more general anharmonic systems, such as heat transport in nonlinear chain models which is one of the ultimate goals in the community of heat transport.

###### Acknowledgements.
This work is supported in part by the MOE T2 Grant R-144-000-300-112 and an NUS Grant R-144-000-285-646.

## Appendix A Derivation of Eq. (26)

In this appendix, we derive Eq. (26) from Eq. (24) with .

We first note that can be obtained using first order perturbation theory, which reads

 ~ϕ∘0=∞∑m=1(ϕ0,L(1)φm)μm−μ0ϕm+a00ϕ0, (48)

where depends on the normalization convention used. However, from Eq. (22), it can be easily shown that the projection of on will not contribute to for any . Then we can just omit this term and regard .

Then it can be noticed that all the eigenfunctions , , as well as the eigenvalues , are explicit functions of . Their -dependence comes from . Therefore, for or ,