Nonequilibrium Green's function method for thermal transport in junctions

0.31 MB·English

by  Jian-Sheng Wang

Nonequilibrium Green’s function method for thermal transport in junctions
Jian-Sheng Wang,1,2,3,∗ Nan Zeng,1 Jian Wang,1 and Chee Kwan Gan2
1Center for Computational Science and Engineering,
and Department of Physics, National University of Singapore, Singapore 117542
2Institute of High Performance Computing, 1 Science Park Road, Singapore 117528
7 3Singapore-MIT Alliance, 4 Engineering Drive 3, Singapore 117576
0
(Dated: 9 January 2007)
0
2 We present a detailed treatment of the nonequilibrium Green’s function method for thermal
transportduetoatomicvibrationsinnanostructures. Someofthekeyequations,suchasself-energy
n
a andconductancewithnonlineareffect,arederived. Aself-consistent mean-fieldtheoryisproposed.
J Computational procedures are discussed. The method is applied to a numberof systems including
one-dimensional chains, a benzene ring junction, and carbon nanotubes. Mean-field calculations of
9
theFermi-Pasta-Ulam model arecompared with classical molecular dynamicssimulations. Wefind
] that nonlinearity suppresses thermal transport even at moderately high temperatures.
h
c PACSnumbers: 05.60.Gg, 44.10.+i,63.22.+m,65.80.+n
e
m
- I. INTRODUCTION The nonequilibrium Green’s function formalism de-
t scribed in this article is a serious attempt to be such
a
t a complete theory. This theory is certainly “first princi-
s Fourier’s law describes the transport of heat in a
. macroscopic object. The calculation of the thermal con- ples” giventhe atomic potentials. Severalsimilar formu-
t
a ductivitythatappearsintheFourier’slawisfundamental lations have already been done at the elastic level with-
m out nonlinear interactions [11, 12, 13, 14]. In Ref.[15],
andimportantinunderstandingthepropertiesofmateri-
we presented a new formulation with the nonlinear in-
- als. Sucha calculationmustbe groundeduponquantum
d teraction treated systematically. Mingo also gave a sim-
mechanics with phonons as basic quasi-particles. This is
n
ilar result [16] for the nonlinear interactions. Our ap-
a non-trivial task and was successfully carried out many
o
proach is essentially a generalization of the nonequilib-
c years ago by Peierls [1, 2] using the Boltzmann equation
[ for phonons. Molecular dynamic (MD) simulation is an- rium Green’s function method in electronic transport
other approach. However, MD results are not correct at [17, 18], with fermions replaced by bosons as basic en-
1
tities. Although the techniques have been extensively
v low temperatures as it is purely a classical approach.
described in the literature, most of them are centered
4 The early works are mostly concerned with bulk ma-
around fermions. In this paper, we give a description
6 terials with many degrees of freedom and periodic lat-
1 tices [3]. In recent years, more attentions are paid to of the method, emphasizing on both the theoretical for-
1 mulation and computational implementation. We begin
heat transport in small or low-dimensional systems such
0 withanelementaryintroductionoftheGreen’sfunctions,
as carbon nanotubes and molecular junctions [4, 5]. In
7 including the contour-ordered Green’s functions. We
nano- to meso-scales, new features come in. In most of
0
thendiscussequationsofmotionoftheGreen’sfunctions,
/ these conditions, the discreteness of the atom is impor-
t Feynmandiagrammaticexpansion,andDysonequations.
a tant. The conceptof phonons developed for the periodic
Wediscusstheheatcurrentandderiveaformulaforeffec-
m latticesissomewhatdifficulttoapplyifasystemdoesnot
tive transmission when there are nonlinear interactions.
- possesstranslationalinvariance,suchasanano-junction.
d Insuchsituation,theBoltzmann-Peierls’approachisnot After the introduction of the method, we present results
n ofone-dimensional(1D)chains,benzenering,andcarbon
applicable.
o nanotubes, and discuss some of the interesting features
There are several alternative methods that have over-
c in such systems.
: come the above problems. First, the Landauer formula
v
[6]is asimpleandcleardescriptionofthe heattransport
i
X in purely ballistic regime at low temperatures. There
II. NONEQUILIBRIUM GREEN’S FUNCTION
are also a few other approaches [7, 8] with different ap-
r METHOD
a proximations such as the description based on density
matrices. One of the features of the existing theories is
A. Definition of Green’s functions
that they work in either the ballistic regime or diffusive
regime, but not both. It is rather difficult to have a
completetheorythatcanencompassbothregimes,apart The nonequilibriumGreen’sfunction methods are dis-
from phenomenological treatments [9, 10]. cussedinRefs.[18,19],andequilibriumGreen’sfunctions
are explained in many textbooks such as Refs. [20, 21].
For completeness, we introduce our notations and give
a quick review of the definitions and properties of vari-
∗http://staff.science.nus.edu.sg/∼phywjs/ ous Green’s functions in this subsection. We define the
2
retarded Green’s function as Also, the relations Gr = Gt G< and Ga = G< Gt¯
− −
are useful. Out of the six Green’s functions, only three
i
Gr(t,t′)= θ(t t′) [u(t),u(t′)T] , (1) of them are linearly independent. However, in systems
−¯h − h i with time translational invariance, the functions Gr and
Ga are Hermitian conjugate of one other:
where the column vector u with component u is not
j
the usual displacement operator from equilibrium but a
Ga[ω]=(Gr[ω])†. (12)
renormalized one, i.e., uj = xj√mj, such that the ki-
netic energy is of the form (1/2)u˙Tu˙. The superscript Soingeneralnonequilibriumsteady-statesituations,only
T stands for matrix transpose. The square brackets are two of them are independent. In this paper, we consider
the commutators. Gr is a square matrix with elements them to be Gr and G<. There are other relations in the
Grjk(t,t′) = −(i/¯h)θ(t − t′)h[uj(t),uk(t′)]i. The phys- frequency domain as well:
ical dimension of Gr is time. By definition, Gr(t,t′)
equals zero when t t′. In equilibrium or nonequilib- G<[ω]† = G<[ω], (13)
rium steady states,≤the Green’s function depends only Gr[ ω] = −Gr[ω]∗, (14)
on the difference in time, t t′. The Fourier transform G<[−ω] = G>[ω]T= G<[ω]∗+Gr[ω]T Gr[ω]∗.(15)
of Gr(t t′)=Gr(t,t′) is de−fined as − − −
−
The last two equations show that we only need to com-
+∞
pute the positive frequency part of the functions.
Gr[ω]= Gr(t)eiωtdt. (2)
Z Inthermalequilibrium,thereisanadditionalequation
−∞
relating Gr and G<:
Notethatweusesquarebracketstodelimittheargument
for the Green’s functions in frequency domain. The di- G<[ω]=f(ω) Gr[ω] Ga[ω] , (16)
mension of the retarded Green’s function in frequency (cid:16) − (cid:17)
domain is time squared. The inverse transform is given where
by
1
f(ω)= (17)
1 +∞ exp(β¯hω) 1
Gr(t)= Gr[ω]e−iωtdω. (3) −
2π Z
−∞ istheBose-Einsteindistributionfunctionattemperature
T = 1/(k β). Equation (16) is obtained by writing the
For notationalsimplicity, we’llset h¯ =1. We canalways B
Green’s functions as a sum of energy eigenstates (known
get the final required expressions by a simple dimension
as the Lehmann representation). Thus in equilibrium,
analysis.
there is only one independent Green’s function; we take
The rest of the definitions are the advanced Green’s
it to be Gr.
function
The contour-ordered Green’s function is a convenient
Ga(t,t′)=iθ(t′ t) [u(t),u(t′)T] , (4) book-keeping to treat the different Green’s functions in
− h i a concise notation. We can consider a contour-ordered
the “greater than” Green’s function Green’sfunctionasafunctionG(τ,τ′)withargumentsτ
and τ′ defined on the complex plane. The contour runs
G>(t,t′)= i u(t)u(t′)T , (5) from slightly above the real axis to + and loops
− h i −∞ ∞
back from + slightly below the real axis to . The
∞ −∞
the “less than” Green’s function contour-ordered Green’s function can be mapped onto
four different normal Green’s functions by Gσσ′(t,t′) =
G<(t,t′)= i u(t′)u(t)T T, (6) lim G(t+iǫσ,t′+iǫσ′), where σ = (1), and G++ =
− h i Gtǫi→s0t+he time-ordered Green’s functio±n, G−− = Gt¯ is
the time-ordered Green’s function the anti-time-ordered Green’s function, G+− =G<, and
G−+ =G>.
Gt(t,t′)=θ(t t′)G>(t,t′)+θ(t′ t)G<(t,t′), (7)
− − In a noninteracting harmonic system with a Hamilto-
nian
and anti-time-ordered Green’s function
1 1
Gt¯(t,t′)=θ(t′ t)G>(t,t′)+θ(t t′)G<(t,t′). (8) H0 = u˙Tu˙ + uTKu, (18)
2 2
− −
Thefollowinglinearrelationsholdbothinfrequencyand theretardedGreen’sfunctioninthe frequencydomainis
time domains from the basic definitions: given by
Gr Ga = G> G<, (9) Gr[ω]= (ω+iη)2I K −1, (19)
− − −
Gt+Gt¯ = G>+G<, (10) where I is an identit(cid:0)y matrix, and η(cid:1) 0+. Adding
Gt Gt¯ = Gr +Ga. (11) a small η helps to choose correctly the→inverse Fourier
−
3
transform integration path in the complex ω-plane, so Governing H+H+H +V+H
Hamiltonians L C R n
that the retarded Green’s function has the required
causal property, Gr(t) = 0 for t < 0. Equation (19) can HL+HC+HR +V G
be derived by an equation of motion method. We note HL+HC+HR G Green’s functions
thattheretardedGreen’sfunctionisasymmetricmatrix g 0
since K is symmetric, but this feature is not preserved t = −∞
with nonlinear interactions. Equilibrium at Tα t = T
t= 0
Nonequilibrium
steady state
B. Model and adiabatic switch-on
FIG. 1: A schematic to illustrate the two adiabatic switch-
ons.
Our system is a junction with a central region and
two leads which serve as heat reservoirs. We treat the
leads explicitly as quasi-one-dimensional periodic lat- turned onslowly,and a steady state of the linear system
tices. This setup is experimentally relevant and con- isestablishedatsometimeT 0. TheGreen’sfunctions
≪
ceptually useful for computation. We consider non- of the linear nonequilibrium system will be denoted by
conducting solid and treat only the vibrational degrees G . For this linear problem, the result does not depend
0
offreedomforheattransport. Letthe(mass-normalized) on T . Finally, the nonlinear interaction H is turned
C n
displacement from some equilibrium position for the jth on, and at time t = 0, a nonequilibrium steady state is
degree of freedom in the region α be uα; α = L,C,R, established,see Fig.1 for illustration. The full nonlinear
j
for the left, center, and right regions, respectively. The Green’s functions will be denoted by G.
quantum Hamiltonian is given by Thedensitymatricesattimet= ,T,andt=0are
−∞
related in the following way in the interaction picture:
= H +(uL)TVLCuC+(uC)TVCRuR+H , (20)
α n
H
α=XL,C,R ρ(T) = S0(T, )ρ( )S0( ,T), (23)
−∞ −∞ −∞
where H = 1(u˙α)Tu˙α+ 1(uα)TKαuα, uα is a column S0(t,t′) = e−i tt′V(t′′)dt′′, (24)
α 2 2 T R
vector consisting of all the displacement variables in re- ρ(0) = S(0,T)ρ(T)S(T,0), (25)
tguiomn.α,Kaαndisu˙αthies stphreincgorrceosnpsotnandtinmg acotrnijxugaantde mVoLmCen=- S(t,t′) = e−i tt′Hn(t′′)dt′′, (26)
T R
(VCL)T isthecouplingmatrixoftheleftleadtothecen-
where is the time-order operator.
tralregion;similarlyforVCR. Wenotethatthedynamic T
matrix of the full linear system is
C. Equation of motion method
KL VLC 0
K =VCL KC VCR . (21)
0 VRC KR In this subsection, we derive the equations of motion
  for Green’s functions. The equations of motion can be
The nonlinear part of the interaction will take the form used to perform systematic perturbation expansion for
the Green’s functions, or as a starting point for mean-
1 1
H = T uCuCuC+ T uCuCuCuC. (22) field approximations. There are at least two ways to de-
n 3Xijk ijk i j k 4Xijkl ijkl i j k l rive the equations of motion for nonequilibrium Green’s
functions. The first is to derive the equations of motion
Thequarticinteractionisimportantinstablizingthesys- for the time-ordered Green’s function and then gener-
tem, as a purely cubic nonlinear interaction makes the alize it to the contour-ordered version by evoking the
energy unbounded from below. structure isomorphism of the two sets of Green’s func-
Weneedtoansweranimportantquestion: whatisthe tions [17]. Another possibility is to consider directly the
distribution of the system? The distribution enters the contour-orderedGreen’sfunction[22]. We’lltakethelat-
definition of the Green’s functions as a density matrix ter approach.
ρ(0) for the average . For an equilibrium problem, We define a general n-point contour-ordered Green’s
h···i
thisisjusttheBoltzmannfactor,butinanonequilibrium function as
situation,itisnotknownandmustbecomputedinsome
way. The adiabatic switch-on gives us a clear concep- Gjα11,,jα22,·,····,·j,nαn(τ1,τ2,···,τn)=
tual framework how this problem can be solved, at least i uα1(τ )uα2(τ )...uαn(τ ) . (27)
formally. We imagine that at t = the system has − hTτ j1 1 j2 2 jn n i
−∞
three decoupled regions, each at separate temperatures, The index α = L,C,R labels the region, j labels the
T , T , and T . The nonlinear interactions are turned degrees of freedom in that region, and τ is the contour
L C R
off. The Green’s functions gα are known and take the variable. The function is symmetric with respect to si-
form of Eq. (19). The couplings VLC and VCR are then multaneouspermutationsofthetriplet(α,j,τ). Weview
4
the function as an analytic function (at least along the ∂ gL(τ′,τ′′) = 0, at the upper limit τ′ = iǫ.
∂τ′ −∞ −
path) in variable τ which varies on the Keldysh contour This requirement is equivalent to gt¯( ) = 0, and
[23]. We also introduce a generalized θ(τ,τ′) function g>( ) = 0, or gr( ) = 0 (as w−ell∞as that their
which is defined to be 1 if τ is later than τ′ along the first−d∞erivatives are zer−o∞). The condition gr( ) = 0 is
contour and 0 otherwise. Its derivative is a generalized consistentwiththe causalityrequirementoft−he∞retarded
δ-function,δ(τ τ′)=σδσ,σ′δ(t t′). Theθ functioncan Green’s functions, but the conditions for gt¯or g> arein-
− −
be used to represent the contour-order operator τ. For sufficient, as the information about temperature of the
T
convenience,we allow for n=0, which is just a constant system is unspecified. We fix the ambiguity by requiring
i. Our definition is slightly different from that of the that
−
usualquantumfieldtheorywhereourprefactorisalways
i (instead of the usual ( i)n/2). g>(t=0)= i uL(0)uL(0)T (33)
− − − h i
To start with, we consider one of the important cor-
relation functions for heat transport, GC,L. The first be the equilibrium value at inverse temperature βL. In
other words, gL should exactly be the contour-ordered
derivative with respect to the second argument always
Green’sfunctionofthefreeleadasdefinedinsubsec.IIA.
leads to a direct differentiation inside the contour-order
Such a choice is consistent with the adiabatic switch-on.
operator since the equal-time coordinates u commute,
A full justification of the final result,
∂GC,L(τ,τ′)
j,l∂τ′ =−ihTτuCj (τ)u˙Ll (τ′)i. (28) GC,L(τ,τ′′)=Z GC,C(τ,τ′)VCLgL(τ′,τ′′)dτ′, (34)
ThesecondderivativeissimilarinthiscaseasuC andu˙L can also be given from a perturbation expansion, as has
alsocommuteatequaltime. Thus,substitutingtheequa- been done in Ref. [17].
tion of motion for uL (having identical form in quantum Equation(34)canbe generalizedtothecaseofGC,C,L
and classical cases), we get with a similarresult. We now considerthe Green’s func-
tions of the form GC,C,···,C involving only the central
∂2GC,L(τ,τ′) part. The aim is to derive a recursion relation for the
j,l = GC,L(τ,τ′)KL GC,C(τ,τ′)VCL.
∂τ′2 − j,m ml− j,m ml central part Green’s functions by eliminating references
Xm Xm
tofullGreen’sfunctionsinvolvingtheleads. Thistypeof
(29)
Green’sfunctionsinvolvestwonewfeatures,(1)wemust
In matrix notation, it is
take care of the derivatives of the θ functions which pro-
∂2GC,L(τ,τ′) ducecommutators[uC,u˙C],(2)wegeneratehigherorder
+GC,L(τ,τ′)KL = GC,C(τ,τ′)VCL. Green’s functions due to the nonlinear interactions.
∂τ′2 −
(30) The contour order leads to n! terms, each of which is
To solve this differential equationonthe contour,we de- a permutation of the displacement variable u from the
fine the (contour-ordered) Green’s function of the left original order 1,2, ,n, multiplied by a term of per-
···
lead to satisfy mutation of the arguments τ from the standard order
θ(τ τ )θ(τ τ ) θ(τ τ ). Differentiation of
1 2 2 3 n−1 n
∂2gL(τ′,τ′′) these−θ function−s lead··s·to a sum−of contractions between
+KLgL(τ′,τ′′)= Iδ(τ′ τ′′). (31)
∂τ′2 − − the variable that we are differentiating and all the other
variables. This reduces the orderofthe Green’s function
We can think of gL as the inverse of the operator
by two. We note that
−
∂2/∂τ′2+KL. BymultiplyinggL(τ′,τ′′)toEq.(30)from
the right, GC,L(τ,τ′) to Eq. (31) from the left, and then [uC(τ),u˙C(τ)]=iδ . (35)
j l jl
subtracting the two equations and integrating over the
variable τ′ along the contour, we get (after integration For the moment, we’ll consider only the cubic nonlin-
by part) earity. Four terms are produced when u¨C is substituted
inside the order operator, due to linear coupling KC,
VCL, VCR, and nonlinear T . The first three terms
GC,L(τ,τ′′)= GC,C(τ,τ′)VCLgL(τ′,τ′′)dτ′+ ijk
Z are similar in structure to that in Eq. (29). The cubic
∂GC,L(τ,τ′) ∂gL(τ′,τ′′) τ′=−∞−iǫ nonlinear interaction increases the value of n by 1. We
gL(τ′,τ′′)+GC,L(τ,τ′) (.32) note that the equation for u¨C is evaluated at the same
∂τ′ ∂τ′ (cid:12)τ′=−∞+iǫ time τ. This would be complicated if we have to keep
(cid:12)
The reason that the terms in the second li(cid:12)ne are iden- trackofwhichτ’sareequalinthe Green’sfunctions. We
tically zero involves some considerations. First, since transfer this information to the coupling Tijk to define
the central part and left lead are decoupled at the time the three-body interaction in contour variable form (for
τ = +iǫ, we must have GC,L(τ, +iǫ) = 0 and the force)
∂ G−C,∞L(τ, +iǫ) = 0. This takes−c∞are of the lower
∂τ′ −∞
bound. Next, we require that the free left lead Green’s T (τ,τ′,τ′′)uC(τ′)uC(τ′′)dτ′dτ′′. (36)
function satisfies the boundary condition gL(τ′,τ′′)=0, −Xjk Z Z ijk j k
5
1 n 1 n 1 n 1 n ture evolutionoperator(scatteringmatrix operator)and
is expressed as:
= + i { +
2 2 2 2
G (τ,τ′) = i uH(τ)uH(τ′)
3 3 3 3 jk − hTτ j k i
= i uI(τ)uI(τ′)e−i HnI(τ′′)dτ′′ , (39)
1 n 1 n − hTτ j k R i0
+ + . . . + } where the displacements refer to the central region; we
2 2 have dropped the superscript C for brevity. The opera-
3 3 tors in the top line are in the Heisenberg picture, while
thatinthesecondlineareintheinteractionpicture. The
Green’s function G of the linear system (when H =0)
0 n
FIG. 2: Recursive expansion rule for Green’s functions. A can be computed from that of the free subsystems (an
vertex with n double lines denotes an n-point Green’s func- integral equation form of Eq. (38)):
tion. ThesinglelinedenotesG0. Asingle-linethree-terminal
vertex is associated with T(τ,τ′,τ′′).
G (τ,τ′)=gC(τ,τ′)+ dτ dτ gC(τ,τ )Σ(τ ,τ )G (τ ,τ′),
0 1 2 1 1 2 0 2
Z
(40)
+∞
Since the contour integral is dτ = σ dt, we
σ −∞ where Σ=ΣL+ΣR,
must define R P R
Σ =VCLgLVLC, (41)
L
Tijk(τ,τ′,τ′′)=Tijkσ′δσ,σ′δ(t t′)σ′′δσ,σ′′δ(t t′′). (37)
− −
similarlyforΣ . ThisDysonequationcanbe derivedby
R
Thisisbecauseourδ-functionhasthestandardproperty considering the first step of the adiabatic switching-on
that δ(τ τ′)dτ =1, butwhen writtenin the ordinary process. Since it is a linear system, the self-energy Σ is
−
timeRt and index σ, we must have the extra σ factor to known exactly.
get the required value of +1. By expanding the exponential in Eq. (39), a series in
The two terms involving the leads produce Green’s the nonlinear interaction strength is obtained. The re-
functions of the formGC,C,...,L and GC,C,...,R, which can ductionintermsoftheunperturbedGreen’sfunctionG
0
bereplacedbyGC,C,...,C usingtheresultEq.(34)derived relies on the fact that Wick’s theorem [24] is applicable
earlier. We can now solve for the center-only Green’s here. Although a formal proof of Wick’s theorem is dif-
function using the Green’s function of the linear system: ficult, we can see that the density matrix ρ(0) must be
quadratic in exponential. This is because the ρ( ) is
∂2 −∞
+KC G (τ,τ′)+ dτ′′ VCLgL(τ,τ′′)VLC + quadratic and the system is linear. The fact that the
(cid:18)∂τ2 (cid:19) 0 Z (cid:16) equation of motion method and the perturbation expan-
VCRgR(τ,τ′′)VRC G (τ′′,τ′)= Iδ(τ τ′). (38) siongiveidenticalresultsisaconfirmationofthevalidity
0
(cid:17) − − of Wick’s theorem.
We then obtain the following simple recursive rules for The expansion contains connected as well as discon-
Green’s functions (as shown in Fig. 2): nected Feynman diagrams. The disconnected diagrams
(1) Replace leg 1 by inserting a nonlinear coupling are constant in time, and give rise to a thermal expan-
T(τ,τ′,τ′′) such that the outer leg is immediately con- sion effect. We can show that these diagrams do not
nected with G while the other two terminals increase contribute to the thermal transport, as they are propor-
0
the order of the Green’s function by 1. Quartic interac- tional to δ(ω) in the frequency domain. The thermal
tion is similar, but the process will increase the order by current formula has a factor of ω which makes it zero.
2. Finally,theconnectedpartoftheGreen’sfunctionsat-
(2) Add imaginary unit i times a sum of the n 1 isfies a similar contour-ordered Dyson equation relating
graphs formed by pairing each leg with leg 1 and c−on- Gc to G0 through a nonlinear self-energy Σn:
necting with the propagator G , multiplied by a (n 2)
0
order remaining Green’s function. − Gc(τ,τ′)=G0(τ,τ′)+ dτ1dτ2G0(τ,τ1)Σn(τ1,τ2)Gc(τ2,τ′).
Z
(3)Symmetrizethegraphs,ifdesired,i.e.,dosteps(1)
(42)
and(2) for every leg,1, 2, , n,add them up, and then
··· In ordinary Green’s functions and in frequency domain
divide by n.
(ω argument suppressed), the Dyson equations have so-
The above rules may be conveniently implemented in
lutions [18]:
a symbolic computer language, such as Mathematica.
Gr = Ga† = (ω+iη)2I KC Σr −1, (43)
0 0 − −
D. Feynman diagrams and Dyson equations G<0 = Gr0Σ<G(cid:0)a0, (cid:1) (44)
−1
Gr = Gr−1 Σr , (45)
The contour-ordered Green’s function can also be ob- c (cid:16) 0 − n(cid:17)
tainedbyaperturbationexpansionoftheinteractionpic- G< = Gr Σ<+Σ< Ga. (46)
c c n c
(cid:0) (cid:1)
6
we have
= 2 i + 2 i + (-8) + (-8) + (-8) σ1δσ1,σ2δ(t1−t2) Tj1,j2,j3Tj4,j5,j6
j3,j4X,j5,j6,σ4
Gσ1,σ4 (t t )Gσ4,σ4 (0)σ dt . (49)
Z 0,j3,j4 1− 4 0,j5,j6 4 4
+ (-4) + (-4) + (-2) + 3 i + (-6) + (-9)
The result can now be transformed to frequency domain
if desired.
+ (-6) + (-12) + (-12) + (-6) + (-6)
E. Thermal current and conductance
+ (-12)
+ (-6) We define the current flow from the left lead to the
+ (-6) + (-6) + (-3) + (-6) + (-4) central region as
I = H˙ (t) . (50)
L L
−h i
FIG. 3: The Feynman diagrams for nonlinear self-energy Σ
n
with the prefactors for graphs of cubic and quartic interac- The interpretation of the current is somewhat subtle as
tions to second order in ¯h. the leads are semi-infinite in extent. The definition is
meaningful. By the Heisenberg equation of motion, we
obtain, at t = 0, I = (u˙L)TVLCuC . The expecta-
L
The connected part of the Green’s function Gc will be tion value can be expresshed in terms ofia Green’s func-
used for G in Eq. (52) below, and the extra subscript c tion G< (t,t′) = i uL(t′)uC(t)T T. Using the fact
will be dropped from this point. that opCeLrators u an−dhu˙ are relatediin Fourier space as
In Fig. 3 we give the Feynman diagrams to “second u˙[ω] = ( iω)u[ω], we can eliminate the derivative and
−
order”ofthenonlinearinteractionsforcubicandquartic get,
terms, where the graphs have explicitly h¯ and h¯2 in the
coefficients;theneglectedgraphshavehigherpowersofh¯. I = 1 ∞Tr VLCG< [ω] ωdω. (51)
ThesediagramsarethesameasthoseinRef.[25],butthe L −2π Z CL
−∞ (cid:0) (cid:1)
interpretation is different. It is for the contour-ordered
Using the result derived earlier relating the mixed lead-
Green’sfunctionshereandfortheretardedGreen’sfunc-
center Green’s function to the center-only Green’s func-
tions of the creation/annihilation operators in Ref. [25].
tionandapplyingthe Langreththeorem[18, 26]towrite
To find the actual expressions for each of the diagrams,
contour-ordered Green’s function in ordinary Green’s
wehavethefollowingsimplerules: (1)labeleachinterac-
functions in frequency domain, we have G< [ω] =
tionvetexlinebytheτ variablesandsiteindices,withthe CL
Gr [ω]VCLg<[ω]+G< [ω]VCLga[ω]. The final expres-
interaction function Tijk(τi,τj,τk) (or Tijkl(τi,τj,τk,τl) CC L CC L
sion for the energy current is
forthe quarticinteraction). (2)Eachline connectingthe
labels is associated with a propagator G . (3) Integrate
0 1 +∞
over all the internal variables along the contour; sum I = dωωTr Gr[ω]Σ<[ω]+G<[ω]Σa[ω] , (52)
over the internal site indices. Most of the integrations L −2πZ−∞ (cid:16) L L (cid:17)
are easy, as T ( ) contains two δ-functions in τ. For
ijk ··· wheretheself-energyduetotheinteractionwiththelead
example, the second graph is
isΣ =VCLg VLC. Forsimplicity,wehavedroppedthe
L L
subscriptC ontheGreen’sfunctionsdenotingthecentral
dτ dτ dτ dτ T (τ ,τ ,τ ) region. Thisresultisthe phononanalogofthe electronic
Z 3 4 5 6 j1,j2,j3 1 2 3
j3,jX4,j5,j6 current formula [17, 18].
G (τ ,τ )T (τ ,τ ,τ )G (τ ,τ ),(47) Since the energy current must be conserved, we have
0,j3j4 3 4 j4,j5,j6 4 5 6 0,j5j6 5 6
I = I . In addition, I is real. We can obtain a
L R L
−
where τ and τ are external variables and the rest are symmetrized expression for the current
1 2
internal. After integrations, we get
1 1 ∞
I = (I +I∗ I I∗)= dω
4 L L− R− R 4π Z
0
dτ T δ(τ τ )G (τ ,τ )T G (τ ,τ ).
j3,jX4,j5Z,j6 4 j1,j2,j3 1− 2 0,j3j4 1 4 j4,j5,j6 0,j5j6 4 4 ωTrn(Gr−Ga)(Σ<R−Σ<L)+iG<(ΓR−ΓL)o, (53)
(48)
NdeoxwσwuesicnagnthreewmriateppiinngreGal(τt,iτm′e) t aGndσσ′t(hte,tb′)r,ancdhτin- whWereedΓeαfin=eit(hΣerαt−heΣrmaαa).l conductance as
σ=±1 −∞∞σdt, and δ(τ −τ′)→σ→δσ,σ′δ(t−t′).RNotin→g κ¯ = lim I = δI , (54)
tPhatG0(Rτ4,τ4)=G0(0)by time translationalinvariance, ∆T→0∆T δT
7
showninFig.4,wherethedoublelinerepresentsthenon-
linear full Green’s function G, while a single line is the
usual bare Green’s function G . With such a choice, we
0
= 2 i + 2 i + 3 i + (-6) find that to next order, all the cubic interaction graphs
are reproduced exactly (including the prefactor), except
that two of the diagrams are missing.
FIG. 4: Mean-field approximation to the self-energy. The Withtheself-consistentmean-fieldapproximation,the
double line denotes the full Green’s function G; the single solution has to be found iteratively. The convergence of
line is for G0. the iterative process becomes an issue. With only cubic
interactions,thesystemofequationswillnotconvergeat
high temperatures. We understand that this is simply
where ∆T is the difference of the temperatures in leads, a consequence of the unstable nature when only cubic
such that TL = T +∆T/2 and TR = T ∆T/2. Using terms are used. We thus add the fourth order potential
−
Eq. (53), we can actually take the limit ∆T 0 by in- whichis essentialforstability. However,the convergence
→
troducingvariationalderivativesoftheGreen’sfunctions, is still difficult at high temperatures, mainly due to the
self-energies, or current as, e.g., singular nature of the Bose distribution.
δG G(T ,T ) G(T,T)
L R
= lim − . (55)
δT ∆T→0 TL TR G. Classical limit, molecular dynamics
−
The difference in the numerator is the value of a func-
At high temperatures, the classical molecular dynam-
tionwhenthe leadsareattwodifferenttemperaturesT
L ics (MD) is a correct and numerically exact method for
and T , respectively, minus the value when the system
R the heat conduction problem. However, the standard
is in equilibrium at T = T = T. With this notation,
L R heat baths used in MD such as Nos´e-Hooverdo not cor-
we can express the conductance in a form similar to the
respondto exactly the actionof the leads in our models.
Landauer formula [27] for the ballistic transport as
A correct heat bath should follow a Langevin dynamics
1 ∞ ∂f(ω) with colored noises. Dhar [29] has given a derivation for
κ¯ = dωωT˜[ω] (56) the junction problem. The idea is to solve the lead de-
2π Z ∂T
0 grees of freedom uL and uR formally first, which is easy
since it is a linear system. The resultis then substituted
with an effective transmission coefficient,
into the equation for the center region. This gives the
1 1 following equation:
T˜[ω] = Tr Gr(Γ + Γ S)GaΓ +
L n R
2 2 − t
(cid:8) (cid:9)
1 1 u¨C = KCuC +F (uC) Σr(t,t′)uC(t′)dt′+ξ +ξ ,
2Tr GaΓLGr(ΓR+ 2Γn+S) , (57) − n −ZT L R
(cid:8) (cid:9) (59)
where the nonlinear effectis reflectedin the extra terms, where Fn is the nonlinear force, Σr is the same retarded
Γ =i(Σr Σa), and self-energy of the lead, Σr = Σr + Σr, used in the
n n− n nonequilibrium Green’s function cLalculatRion, but in the
δΣr δΣa δΣ< ∂f −1 time domain. By integration by part, we could also cast
S =i f n n n . (58) Eq. (59) in a standard generalized Langevin form [30]
(cid:20) δT − δT − δT (cid:21)(cid:18)∂T(cid:19)
(cid:0) (cid:1) involving a damping force, but KC has to be shifted to
KC+Σr[ω=0]. An extra contribution from the left lead
Equation (57) is a generalization of the Caroli formula
is
[28] for ballistic transport. We’ll discuss in Sec. IIIB
how the variational derivatives of the self-energy can be ∂gr(t,T)
computed. ξL(t)=VCL(cid:18)gLr(t,T)u˙L(T)− L∂T uL(T)(cid:19),
(60)
wheregr isthe retardedGreen’sfunctionofthe left lead
L
F. Mean-field approximation
in the time domain. The expression for the right lead
ξ is similar. A time T in the lower limit in the inte-
R
We haveseenin Ref.[15]that high-orderFeynmandi- gral and Eq. (60) is some time in the remote past where
agramsareimportantat hightemperatures. However,it the system is decoupled and the leads are in respective
is computationally very expensive to include these high- thermal equilibrium. We’ll take the limit T . We
→−∞
order diagrams,although we can do calculations on very turn Eq. (59) into a stochastic differential equation by
small systems such as a three-site model. The standard requiring that uL(T) and u˙L(T) are random variables
technique is to use a mean-field approximation which in satisfying the Boltzmann-Gibbs distribution
some sense partially takes into account these high-order
graphs. We propose an approximationto the self-energy P(uL,u˙L) e−βL(21(u˙L)Tu˙L+12(uL)TKLuL). (61)
∝
8
Thus, ξ is colored noise with the properties s k
L 00
←
e k
11
←
ξ (t) = 0, (62) α k
L 01
h i 1 β ←αT
βLhξL(t)ξL(t′)Ti = VCL(cid:16)g˙Lr(t,T)KL g˙Lr(t′,T)T g ←←[(ω+iη)2I−e]−1
do
+gLr(t,T)gLr(t′,T)T VLC (63) s′ s+αgβ
∞ (cid:17) if ←s′ s <ǫ exit
= −Zt ΣrL(t′′,t′)dt′′, (t>t′). (64) αe′′|←−eα+g|ααgβ+βgα
←
β′ βgβ
The last line of the above equation is a kind of
s ←s′
“fluctuation-dissipation” theorem. We note that the dot
←
e e′
in g denotes the derivative with respect to T. For a sen-
α← α′
sible heat bath, the noise should not depend on T, but
←
only on the difference t t′. Indeed, this can be done if β β′
− g ←[(ω+iη)2I e]−1
werewritetheresultinvibrationaleigen-modes. We can
← −
end do
write the Fourier transform of the noise correlation as
gr [(ω+iη)2I s]−1.
← −
1
Im2VCL (KL)−1/2 VLC, (65)
(cid:18) (ω2 iη)I KL(cid:19)
− − B. Fast Fourier transform and variational
derivatives
where the squarerootofa matrixis defined by its eigen-
value decomposition taking only the positive eigen fre-
quencies. This is well defined since KL is positive defi- The diagrams in the mean-field approximation can be
nite. The properties of the right lead noise ξ are anal- computed in the time domain easily and efficiently; the
R
ogous. Such a set of stochastic differentialequations can convolutions in frequency domain become simple mul-
be simulated on a computer. tiplications in time domain. Thus we have adapted a
fast Fourier transform method for the calculation. This
reduces the computational complexity, particularly for
the quartic interactionterms. Typically 103 Fourier grid
III. IMPLEMENTATION DETAILS
points are needed for good convergence. Using the Mes-
sage Passing Interface (MPI), a parallel implementation
A. Surface Green’s functions
with data partition is made. Good linear speed-up is
obtained upto 64 processes.
The Green’s functions of the free leads gr arerequired Inthemean-fieldtreatmentoftheheattransport,there
in calculation for the lead self energies Σ. Although the isaneedtofindthevariationalderivativesoftheGreen’s
lead Green’s function is a semi-infinite matrix, due to a functionsandself-energyfortheconductancecalculation.
localized coupling VLC or VCR, only a corner set of ele- This is done by solving iteratively the set of equations
ments is required. These submatrix elements are termed
δGr =GrδΣrGr, (67)
surface Green’s functions. The surface Green’s function n
can be computed rather quickly by recursive iterations δG< =GrδΣrG<+Gr δΣ<+δΣ< Ga+G<δΣaGa(.68)
n n n
[31] or decimation. For completeness, we give an algo- (cid:0) (cid:1)
The expressions for δΣr and δΣ< can be obtained by
rithm in pseudo-code form. The readeris referredto the n n
differentiating the self-energy results directly. For con-
literature for the derivation of such algorithm.
vergencecontrol,linearmixingofthenew andoldvalues
The programtakes three square matrices of equal size
is used.
for k , k , and k . The matrix k is the block imme-
00 11 01 00
diately adjacent to the center, while k , and k = kT
11 01 10
arerepeatedforthe semi-infinitechainofthe lead. They
IV. APPLICATIONS
formthe whole dynamicmatrixofthe leadinthe follow-
ing way:
A. One-dimensional cubic onsite model
k k 0
00 01
··· One important question to ask about the nonequilib-
 k10 k11 k01 0 
KR = 0 k k k . (66) rium formulation is that whether the theory gives bal-
10 11 01
  listic and diffusive heat transport in a single framework.
 0 k10 ...  Clearly,if the nonlinear self-energy Σn can be computed
 ··· 
accurately,thereisnodoubtthattheanswerisyes. How-
The left arrow ‘ ’ below denotes assignment; ǫ is an ever, this is not clear if only the leading order perturba-
←
error tolerance. tiontheoryisused. InRef.[32],Spohngivesaderivation
9
of a Boltzmann equationfor the phonons in a simple 1D Σ = Σ + Σ ; Σ = VCLg VLC, similarly for Σ .
L R L L R
chain with onsite cubic nonlinearity. It is clear that in The matrix elements of Σ< are all zero except that
the Boltzmann equation framework, we obtain Fourier’s Σ = K2g< and Σ = K2g<. Using these results,
11 L NN R
lawanddiffusiveheattransport. Wedemonstrateinthis after some calculation, we can write
subsection that diffusive heat transport is not observed
if only perturbative result is used. This subsection also G <[ω]= fLλj1−l+fRλ1l−jΘ(ω). (76)
serves as an example of illustration of the concepts dis- 0jl (λ λ )K
1 2
−
cussed earlier by a simple and concrete model.
The Θ(ω) function is defined to be 1 inside the phonon
We consider a 1D chainwith inter-particlespring con-
bandandzerooutsidetheband;moreprecisely,Θ(ω)=1
stant K and onsite spring constant K . We take a sim-
ple cubic onsite potential for the inte0raction, 1t u3, if K0 <ω2 <4K+K0, and is 0 otherwise.
3 j j Due to the onsite cubic nonlinearity, the leading or-
for j belonging to the central region. This meanPs that
der nonlinear self-energy calculation is simple. The first
T =tδ δ . Apart from the nonlinear interaction, the
ijk ij ik
graph contribution is
left lead, central region, and the right lead are identical.
The classical equation of motion is 1Σ σσ′(t)=2it2 G σσ′(t) 2. (77)
njl 0jl
u¨ =Ku +( 2K K )u +Ku t u2, (69) (cid:2) (cid:3)
j j−1 − − 0 j j+1− j j The contribution from the second graph is a constant in
ω:
where t = t for 1 j N and t = 0 for j < 1 or
j j
≤ ≤
j >ThNe.retarded Green’s function Gr can be obtained by 2Σnσjlσ′[ω]=2it2σδσ,σ′δjl σ′′G0σjsσ′′[0]G0σss′′σ′′(0).
0 σX′′,s
solving the linear equation[(ω+iη)2−K˜]Gr0 =I, where (78)
(infinite in both directions) matrix K˜ is 2K+K0 on the Both the summation over σ′′ and over index s can be
diagonal and K on the first off-diagonals. The explicit performed analytically. We find
−
solution is
λ|j−l| 2Σnσjlσ′[ω]=2it2σδσ,σ′δjlG¯rj[0]G¯t(0), (79)
G r [ω]= 1 , (70)
0jl (λ λ )K where
1 2
−
1+λ λj λN−j+1
where λ1 and λ2 are the roots of the equation G¯r[ω] = 1− 1− 1 , (80)
j (1 λ )(λ λ )K
1 1 2
Kλ−1+Ω+Kλ=0, Ω=(ω+iη)2 2K K . (71) − −
0 1+(f +f )Θ(ω)
− − G¯t[ω] = L R , (81)
Note that λ1λ2 = 1. We define the two roots such that (λ1−λ2)K
|λ1T|h<e1suarnfadce|λG2r|e>en1’s. functiongrsatisfiesasimilarequa- G¯t(0) = 1 +∞G¯t[ω]dω. (82)
2π Z
tion as Gr except that it is semi-infinite in extent. We −∞
0
consider the left lead (j 0). The result for the right For numerical accuracy, we compute the conductance
≤
lead is identical. Since the matrix VLC is nonzero only directly instead of using finite differences. Let the tem-
for one corner element, we need the g0 =g00 component perature on the left lead be TL = T +∆T/2 and right
of the matrix gr. Consider only the j = 0 column, the lead be TR =T ∆T/2. We define the thermal conduc-
−
equations for gr in component form are tance as that in Eq. (54) and the effective transmission
by Eq. (56). For the expression for effective transmis-
Ωg0+Kg−1 =1, (72) sion T˜[ω], we also need to take the limit ∆T 0. The
→
advantage of doing this is that the leading terms corre-
Kg +Ωg +Kg =0, j = 1, 2, (73) spondingtoequilibriumGreen’sfunctionscancelexactly,
j−1 j j+1
− − ··· since there should be no heat current when T = T .
L R
Substituting the trial solution g =cλ−j, we find that We can also equivalently define the temperature varia-
j 1
tion through the expansion:
λ
1
g0 =−K. (74) Gr =Gr + δGr∆T +O(∆T2), (83)
eq δT
The “lessthan” Green’s function for the left leadis then
andsimilarlyforG<. Wecansimplifytheeffectivetrans-
(c.f. Eq. (16))
mission as
Imλ
gL< =2ifLImg0 =−2ifL K 1, (75) T˜[ω] = 21Trn(Gr−Ga)(Σ<R−Σ<L)+
where fL is the Bose-Einstein distribution function at iG<(Γ Γ ) /(f f )
the temperature of the left lead. The G<0 Green’s func- R− L o L− R
tion can be obtained with the equation G< = GrΣ<Ga, = iK(Imλ )(FL FR ), (84)
0 0 0 1 11− NN
10
0.3 0.25
0 40
0.25 0.2
0.2
2
0.2
K) K)0.15
W/ W/
-90 0.15 0.5 -90 20 7
k (1 k (1 0.1 4
0.1 0.7
1.0 0.05 10
0.05
2.0
0 0
0 500 1000 1500 2000 0 500 1000 1500 2000 2500
T (K) T (K)
FIG. 5: Thermal conductance of the cubic onsite model as a FIG. 6: Thermal conductance as a function of temperature
function of temperature for several values of the cubic cou- for several system lengths N with t=0.5 eV/(˚A3u3/2).
pling strength t = 0, 0.2, 0.5, 0.7, 1.0, and 2.0 eV/(˚A3u3/2).
The length of center chain is N =5.
where the F function is defined as
of the thermal conductance κ¯ of a short nonlinear chain
1 δ(Gr Ga) δG< ∂f −1 of length N = 5. The t = 0 curve is the ballistic Lan-
FL = (Gr Ga)+ f − .
2 − h δT − δT i(cid:18)∂T(cid:19) dauerformularesult,therestofthecurvesshowtheeffect
(85) of nonlinearity. The prominent feature here is that the
FR is the same except that the first term is negative. In conductance decreases quickly as the nonlinear coupling
evaluating the F functions, we should set T = T = tincreasesfrom0to2. Inparticular,thesystembehaves
L R
T as we have already taken the ∆T 0 limit. The likeathermalinsulatoriftheon-sitenonlinearityislarge
variations required in Eq. (85) can be c→omputed by the enough. The thermal conductance is small at low tem-
variation in nonlinear self-energy, which in turn can be peraturesandathightemperatures,aswasfoundexper-
computed by the variationin the linear Green’s function imentally. At the extremely high temperatures, the con-
G . Due to the Dyson equation, the variation of the ductance decreases with temperature according to 1/T
0
retarded Green’s function has a simple result: (data not shown). This agrees with several other theo-
ries of thermal conductance at high temperatures.
δGr δΣr
=Gr nGr. (86)
δT δT Can we see diffusive transport in this model (under
the perturbative approximation)? The answer is no. In
The variation in G< is rather complicated:
Fig.6,weshowthesizedependenceforafixedt. Attem-
δG< peratures below 200K, the results are essentially inde-
δG< = δT ∆T =Gr ΣrnδG<0 + pendent of the lengths, showing ballistic heat transport.
(cid:8) At around 2000K,from sizes 2 to 10, we see decrease of
δΣr G¯<+Gr(Σ<Ga+ΣrG¯<) h.c.
n 0 n n 0 − the thermal conductance roughly as 1/N, but this trend
+δG(cid:2)<+Gr(δΣ<+ΣrδG<Σa)(cid:3)G(cid:9)a, (87) doesnotcontinue,andtheconductancesaturatestosome
0 n n 0 n
value for systems larger than 20. This feature is generic;
whereG¯<0 =G<0(I+ΣanGa),h.c. standsfortheHermitian the FPU model with first-order perturbation treatment
conjugate of the preceding term, and δGr0 =0, alsohasasimilarbehavior. Thecurvesalsodisplaysome
oscillatorystructure. Thenumberofoscillationsappears
δG0<jl = λ1j−l−λ1l−j Θ(ω)∂f, (88) to be roughly equal to half of the chain length N. The
δT 2(λ1 λ2)K ∂T feature is damped out as the chain becomes longer. It
−
4it2 +∞ would be nice if a molecular dynamics simulation could
δΣnrj,l<[ω] = 2π Z dω′G0rj,l<[ω−ω′]δG0<jl[ω′].(89) confirm these features.
−∞
We have made analytic calculations as far as pos- The failure to see diffusive behavior may be under-
sible. The rest of integrations and matrix inversions standable as our result takes into account only the first-
are calculated numerically. We use physical units with order perturbation in the nonlinear interaction. Other
K = 0.625eV/(˚A2u), K = 0.1K, (u for unified atomic features seem qualitatively correct in comparison with
0
massunits). Figure5showsthetemperaturedependence real systems; this is very encouraging.