Effects of electron-phonon interaction on thermal and electrical transport through molecular nano-conductors

3.3 MB·English

by  Jing-Tao Lü

Effects of electron-phonon interaction on thermal and electrical transport through
molecular nano-conductors
Jing-Tao Lu¨,1,∗ Hangbo Zhou,2,3 Jin-Wu Jiang,4 and Jian-Sheng Wang2
1School of Physics, Huazhong University of Science and Technology, 430074 Wuhan, People’s Republic of China
2Department of Physics and Center for Computational Science and Engineering,
National University of Singapore, 117551 Singapore, Republic of Singapore
3NUS Graduate School for Integrative Sciences and Engineering,
National University of Singapore, 117456 Singapore, Republic of Singapore
4Shanghai Institute of Applied Mathematics and Mechanics,
Shanghai Key Laboratory of Mechanics in Energy Engineering,
Shanghai University, 200072 Shanghai, People’s Republic of China
(Dated: January 26, 2015)
The topic of this review is the effects of electron-phonon interaction (EPI) on the transport
5 propertiesofmolecularnano-conductors. Anano-conductorconnectstotwoelectronleadsandtwo
1 phonon leads, possibly at different temperatures or chemical potentials. The EPI appears only
0 in the nano-conductor. We focus on its effects on charge and energy transport. We introduce
2 three approaches. For weak EPI, we use the nonequilibrium Green’s function method to treat it
perturbatively. We derive the expressions for the charge and heat currents. For weak system-lead
r
a couplings, we use the quantum master equation approach. In both cases, we use a simple single
M level model to study the effects of EPI on the system’s thermoelectric transport properties. It is
also interesting to look at the effect of currents on the dynamics of the phonon system. For this,

4 we derive a semi-classical generalized Langevin equation to describe the nano-conductor’s atomic
2 dynamics, taking the nonequilibrium electron system, as well as the rest of the atomic degrees of
freedom as effective baths. We show simple applications of this approach to the problem of energy

] transfer between electrons and phonons.
l
l
a
PACSnumbers: 85.35.Gv,85.85.+j,85.65.+h,05.60.Gg,73.63.-b
h
-
s
e I. INTRODUCTION motions15. Itisfundamentallydifferentfromthestochas-
m
tic Joule heating. These advances have motivated the
. development of methods treating current-induced forces
t Electron-phonon interaction (EPI) is one of the most
a and Joule heating on the same footing17–19,34.
m important many-body interactions in condensed-matter
and molecular systems1, responsible for a variety of phe- Equally significantly, there has been an increasing in-
-
d nomena, from electrical, thermal conduction, supercon- terestinthethermoelectricpropertiesoflowdimensional
n ductivitytoRamanscattering,polaronformation,justto systems50–54. A starting point of the theoretical treat-
o listafew2–8. Itseffectsontheelectrical,thermal,andop- ment is to ignore the effect of EPI, and study the trans-
c tical properties of bulk semiconductors and metals have port of electrons and phonons separately. But how im-
[
been intensively studied along with the development of portant the effect of EPI is is a pertinent question, on

2 many-body theories and experimental techniques. Re- which much of recent work is devoted to55–59. Here, we
v cent advances in experimental fabrication of meso- and willlookatthisproblemusingthevariousapproacheswe
3 nano-scopicstructureshavegeneratedtremendousefforts have developed.
4
in understanding the effects of EPI on transport proper-
3 EPI is a genuine many-body interaction, the exact
ties of reduced-dimensional systems9–11.
6 treatment of which is challenging, if possible at all. One
0 Ofspecialinterestarecurrent-inducedforcesandJoule natural approach is to perform perturbation calculation
.
1 heating in low-dimensional systems, especially in molec- over a certain small parameter. In the most common
0 ular nano-conductors12–49. On the one hand, the electri- multi-probe transport setup (see Fig. 1 and Sec. II), this
5 cal transport signature of EPI is an invaluable spectro- smallparametercanbechosenaccordingtothestrength
1
scopictooltostudythestructuralinformationofmolecu- ofEPI.Thisstrengthcanberoughlycharacterizedbythe
:
v lar nano-conductors22,24. On the other hand, these pro- ratio between two time scales: the first one corresponds
i cesses are crucial in maintaining the stability of these to the phonon period, and the second one corresponds
X
conductors25, relevant to the continuous scaling down to the electron dwell time60 in the nano-conductor. If
r
of modern electronic devices. Different theoretical ap- the time electrons spend in the nano-conductor is much
a
proaches have been developed to study these problems, shorterthanthephononperiod,thesystemisintheweak
in many cases separately. Recently, it was realized that EPI regime. The small parameter is the EPI matrix. In
non-conservative nature of current-induced forces pro- the other limit, the coupling of the nano-conductor to
videsanalternative,deterministicwayofenergytransfer electrodes is the small parameter, over which one can
betweenelectronsandphonons,ormoregenerallyatomic perform the perturbation expansion.
2
In this review, we summarize our own effort in devel- where x represents the set of all coordinates of the elec-
opingand/orutilizingdifferenttheoreticalapproachesto trons, R the positions of all the ions, and α is the elec-
study the aforementioned problems in different parame- tronic state quantum number. The eigen-functions and
ter regimes. We discuss some relevant results when pos- the eigenvalues depend on R parametrically.
sible, but we make no effort on reviewing all of them We assume an orthonormal set {φ } that satisfies
α
considering the huge amount of literature. The paper Eq. (2) has been obtained. To take into account the
is organized as follows: In Sec. II, we give a brief in- effects of the ions, we consider a trial full wavefunction
troduction of the EPI problem starting from the Born- in a factored form
Oppenheimerapproximation. Wethenintroducethesys-
tem setup and Hamiltonian we use in this paper. In Ψ(x,R)=φ (x;R)χ (R)=|αβ(cid:105), (3)
α β;α
Sec. III, we briefly summarize our use of the nonequi-
librium Green’s function (NEGF) method to study elec- and consider the variational solution4 of the full Hamil-
tron, phonon transport and their interaction perturba- tonian, min (cid:104)Ψ|H|Ψ(cid:105), subject to the normalization
χ
tively. We consider several applications of the method. (cid:104)χ|χ(cid:105) = 1. This variational approach is equivalent to
The first one is the effects of EPI on the thermoelectric omitting the off-diagonal elements (which is the Born-
transportcoefficientsinasinglelevelmodel. Thesecond Oppenheimerapproximation,seeRef.2,App.VIII),giv-
one is the heat transport between electrons and phonons ing an equation for the ions
due to EPI. The use of simple models enables us to ap-
(cid:16)
proachtheproblemssemi-analytically. Thelastexample Pi+Ee(R)+(cid:104)φ |Pi|φ (cid:105)
α α α
is a numerical study of the Joule heating and phonon-
(cid:126)2 (cid:17)
drag effect in carbon nanotubes. In Sec. IV, we consider − (cid:104)φ |∇ |φ (cid:105)·∇ χ=Eχ, (4)
α R α R
the case of strong EPI using the quantum master equa- mp
tion (QME) approach. After reviewing the earlier work,
where (cid:104)···(cid:105) means the x-dependence is integrated out
the same thermoelectric transport model is re-visited fo-
but still R-dependent; ∇ is a multi-dimensional gradi-
cusing on how the strength of EPI affects the results. In R
ent operator with respect to R. Since the left-hand side
Sec.V,wefocusonthecurrent-induceddynamics. Based
depends on the electronic quantum number α, the full
on the Feynman-Vernon influence functional approach,
eigen-energyE andfunctionsalsodependonαparamet-
we derive the semi-classical Langevin equation, taking
rically, e.g., we may write E .
intoaccounttheequilibriumphononandnonequilibrium β;α
Ifweassumethattheelectronsareinitsinstantaneous
electron baths. The final section is our conclusion and
ground state, the ions move in a potential surface gen-
remarks.
erated by the electrons. There are no explicit electron-
phononinteraction(EPI)terms. ToaccountfortheEPI,
weneedtogobacktothebasis,Eq.(3),andconsiderthe
II. BORN-OPPENHEIMER APPROXIMATION
matrix elements
AND ELECTRON-PHONON INTERACTIONS
(cid:104)αβ|H|α(cid:48)β(cid:48)(cid:105). (5)
To discuss the meaning and formulation of the elec-
tron and phonon systems and their mutual interac- The off-diagonal terms are interpreted as the EPI5,6,
tion, we need to start from the Born-Oppenheimer
whicharesmall. Iftheoff-diagonalsareomitted,theelec-
approximation2,3. Consider an electron-ion system with
trons stay in a given quantum state α. The off-diagonal
a total Hamiltonian H = Pi +He, where Pi is the ki- termsdescribethescatteringoftheelectronstodifferent
netic energy operator for the ions, and He = Pe + U state α(cid:48). If ion displacements are small, the most impor-
is electron Hamiltonian with kinetic energy of the elec-
tant contribution is from the linear term in the displace-
trons, Pe, and potential energy U = Uee +Uei +Uii, ment
whichincludestheCoulombinteractionsamongtheelec-
trons and ions. Since the ions are much heavier than the (cid:126)2
electrons, one can treat the ion kinetic energy term as a − (cid:104)φαχβ;α|∇R|φα(cid:48)(cid:105)·∇R|χβ(cid:48);α(cid:48)(cid:105), (α,β)(cid:54)=(α(cid:48),β(cid:48)).
m
p
small perturbation with the expansion parameter2 (6)
These off-diagonal matrix elements can be used, e.g., in
(cid:18) (cid:19)1/4
m a Fermi-Golden rule calculation of scattering processes.
e , (1)
m However, the identity (in the sense of effective Hamil-
p
tonians) of the electrons and phonons and their mutual
where m is the mass of an electron and m mass of an interaction are not at all clear. Although EPI plays ma-
e p
ion (assuming all have the same mass). If the ions are jor role in many physical processes7, such as electronic
considered infinitely heavy, the ions will not move and transport and superconductivity, its conceptual founda-
the electron wavefunctions satisfy tionisstillnotverysolid. WithintheBorn-Oppenheimer
scheme, it is not clear at all how to transform the origi-
Heφ (x;R)=Ee(R)φ (x;R), (2) nalHamiltonianH intoaformofanelectronsystemand
α α α
3
independent phonon system and their interaction unam- asimilarpartitioningforK usingthenotationofRef.68.
biguously. The problem is related to the fact that in de- The EPI takes the form
riving the phonon Hamiltonian (the potential surfaces),
(cid:88) (cid:88)
the effect of electrons is already used. Thus, putting the H (d,uC)= MkuCd†d = uCd†Mkd. (10)
ep ij k i j k
electrons back amounts to double counting, see Refs. 4 ijk k
and 8 for some of the modern treatments.
Insteadofpursuingaself-consistenttheoryofEPIfrom We assume that the EPI appears only in the central re-
the Born-Oppenheimer approximation, here in this re- gion. A schematic representation of the system setup is
view,andalsoinmanyofthepracticalapplications61–65, shown in Fig. 1.
we adopt a phenomenological point of view, and use the Theseparationoftheelectronandphononleadsmakes
model Hamiltonians as given below in Eqs. (8) and (10). the theoretical development easier. In reality, they could
Focusing only the term linear in the displacements away either be physically separated, or built into one. For ex-
from the equilibrium positions of the ions, we can think ample, one electrode could serve both as an electron and
of the single electron Hamiltonian H below having a R- a phonon lead, but we assume that we have independent
dependence. Taylor expanding it, Re=R0+u/√mp, we control over their temperatures Teα and Tpα, α=L,R.
obtain
1 ∂H (R)
Mk = √ (cid:104)i| e |j(cid:105), (7)
ij m ∂R
p k
where |j(cid:105) is the single particle state when ionic sys-
tem is in equilibrium position R0. The extra factor of
square root of ion mass m is because of our conven-
p
tion of displacement variable u. This form of interaction
is intuitively understandable and originally proposed by
Bloch66. In Chap. 4 of Ref. 6, a derivative from Eq. (6)
to (7) is given, but the reasoning does not seem to be
rigorous.
Thus, our starting point of a derivation is a tight-
FIG. 1. Model system considered in this review. The center
binding Hamiltonian for the electrons, harmonic cou-
device, including both electrons, phonons, and their interac-
plings for the phonons and a standard EPI term. They tions, is coupled with two electron and two phonon leads.
are taken as given and exact. The charge redistributions Each electron lead is characterized by its chemical potential
and self-consistency for the electrons are not part of the µ and temperature Tα, and each phonon lead by tempera-
α e
discussion. Symbolically, the total many-body Hamilto- ture Tα.
p
nian is given as
H =H0+H +H , (8)
tot e p ep
III. WEAK EPI REGIME: NONEQUILIBRIUM
where the electron part is H0 = c†Hc, the phonon GREEN’S FUNCTION METHOD
e
part H = 1(pTp + uTKu) + V (uC). The variable
p 2 √ n
u is mass normalized, uj = mj(Rj − Rj0). Because A. Theory
of this, the conjugate momentum is p = u˙. V is
n
the nonlinear force contribution. c is a column vec-
We first consider the case where EPI is weak, so that
tor of the electron annihilation operators, which we can
we can perform a perturbation expansion over the inter-
separate into three regions, the left, center, and right,
action matrix M. In order to do so, we use the NEGF
c = (cL,d,cR)T, T stands for matrix transpose. Sim-
method. Detailed introduction is given in our previous
ilarly u = (uL,uC,uR)T. Accordingly, the matrices H work32,68–70. This section can be considered as an appli-
and K are partitioned into nine regions (submatrices),
cation of the general approach developed in Refs. 68-69
e.g.,
to the EPI problem. We use similar notations therein,
  and only give a brief outline of the approach here.
HL HLC 0
We denote the electron device Green’s function with-
H =HCL HC HCR , (9) out and with EPI by G and G, the corresponding
0 HRC HR 0
phonon Green’s functions by D and D, and the lead
0
Green’s functions without coupling to the center as g
α
such that He0 = HeL+HeR+HeC +Ve, with Ve = VeL+ and dα, respectively. The couplings of the device with
VR, VL = cL†HLCd+H.c.. Note that we assume no the leads and that between the electrons and phonons
e e
interaction between the left and right leads (See Ref. 67 aredescribedbyself-energies,withΣandΠrepresenting
for transport when there is a lead-lead coupling). We do that of electron and phonon, respectively. For example,
4
we define the time-ordered electron Green’s function in- theinteractionHamiltonian H , using Feynmandiagra-
ep
cluding EPI on the Keldysh contour [Fig. 14 (b)] matics. The interacting Green’s functions are expressed
using similar Dyson equations as Eqs. (13-14),
i
G (τ,τ(cid:48))=− (cid:104)T c (τ)c†(τ(cid:48))(cid:105). (11)
ij (cid:126) C i j G(1,2)=G (1,2)+G (1,3)Σ (3,4)G(4,2), (15)
0 0 ep
Here, τ/τ(cid:48) istimeonthecontour, andi/j isindexofthe D(1,2)=D0(1,2)+D0(1,3)Πep(3,4)D(4,2). (16)
electronic states. The contour time order operator T
C Here, Σ and Π are electron and phonon self-energies
puts the operators later in the contour to the left. The ep ep
due to EPI. Using Eq. (12), at steady state, we can get
average (cid:104)·(cid:105) is with respect to the density matrix of the
thefollowingusefulrelationsinenergy/frequencydomain
full Hamiltonian.
The contour ordered Green’s function can be divided
into different groups according to the spatial position of Gr(ε)=(cid:2)(ε+i0+)I−HC −Σrtot(ε)(cid:3)−1, (17)
i/j, similar to the Hamiltonian. The most interesting Dr(ω)=(cid:2)(ω+i0+)2I−KC −Πr (ω)(cid:3)−1, (18)
one is GC, where i and j are both at the center device tot
Σr (ε)=Σr(ε)+Σr (ε), (19)
region. At the same time, it can be written as a 2×2 tot b ep
matrix in time space Πr (ω)=Πr(ω)+Πr (ω). (20)
tot b ep
(cid:18) Gt(t ,t ) G<(t ,t )(cid:19) Weuseεfortheenergyofelectronandω fortheangular
G(τi,τj)= G>(ti,tj) Gt¯(ti,tj) , (12) frequency of phonon, respectively, and I is the identity
i j i j
matrix. Togetanexpressionforthecurrent,wealsoneed
with Gt, Gt¯, G>, G< the time-ordered, anti-time- the greater and lesser version of the Green’s functions71
ordered, greater and lesser Green’s functions. The re-
tardedandadvancedGreen’sfunctionsareobtainedfrom G>,<(ε)=Gr(ε)Σ>to,t<(ε)Ga(ε), (21)
them, i.e., Gr = Gt−G<, and Ga = G<−Gt¯. For the D>,<(ω)=Dr(ω)Π>,<(ω)Da(ω). (22)
tot
definition and relations among these Green’s functions,
we refer to the book by Haug and Jauho71, and our pre- Theelectricalcurrent(Ie)isexpressedasthechangerate
vious publications32,68,69. oftheelectronnumberinoneoftheleads(Nα)timesthe
TocalculatetheGreen’sfunctions, weuseaprocessof charge of electron (−e). For example,
two-step adiabatic switch on. We start from the decou-
pled system and leads. Each of the electron and phonon dNL(t)
I =−e(cid:104) (cid:105)
e
leads is at its own equilibrium state, characterized by dt
the temperature Tα and/or chemical potential µα. The =−2eImTr(cid:104)VLC(cid:104)cL†(t)d(t)(cid:105)(cid:105)
corresponding equilibrium Green’s functions can thus be (cid:126) e
defined according to the equilibrium canonical distribu- =2eReTr(cid:2)VLCG<,CL(t=0)(cid:3). (23)
e
tion. The initial state of the system is arbitrary and not
important in most cases (e.g., for steady state). Itcanbe expressedbythe Green’sfunctionofthecenter
At the first step, we switch on the interaction of the region and the lead self-energies71–73,
center Hamiltonian with the electron and phonon leads.
Wewaituntiltheelectronandphononsubsystemreaches I = e (cid:90) +∞ dεTr(cid:2)G>Σ<−G<Σ>(cid:3). (24)
theirownnonequilibriumsteadystate, sincethetemper- e (cid:126) 2π L L
−∞
ature and/or chemical potential of each lead can be dif-
ferent. The two subsystems are quadratic and exactly Similarly for the heat current carried by electrons (I )
h
solvable, and we get the non-interacting center Green’s and phonons (I )
p
functions G and D from the Dyson equation (we omit
0 0
the superscript C) I = 1 (cid:90) +∞ dε(ε−µ )Tr(cid:2)G>Σ<−G<Σ>(cid:3), (25)
h (cid:126) 2π L L L
G (1,2)=g (1,2)+g (1,3)Σ (3,4)G (4,2), (13) −∞
D00(1,2)=dCC(1,2)+dCC(1,3)Πbb(3,4)D00(4,2). (14) Ip =−(cid:90) +∞ d4ωπ(cid:126)ωTr(cid:2)D>Π<L −D<Π>L(cid:3). (26)
−∞
Here, we have used a single number to represent the ma-
trix indices and contour time arguments, i.e., G (1,2)= We have defined the positive current direction as elec-
0
G (τ ,τ ). Summation or integration over repeated trons going from the lead to the center. We dropped
0j1j2 1 2
indices is assumed. g (d ) is the center electron the argument of the Green’s functions for simplicity. We
C C
(phonon)Green’sfunctionwithoutcouplingtotheLand ignore the spin degrees of freedom, since it is not rele-
R leads. The self-energy Σ =Σ +Σ includes contri- vanthere. Currentsoutoftherightleadareobtainedby
b L R
butionsfromLandR,withΣ (1,2)=HCαg (1,2)HαC; replacingindexLbyR. Onecansymmetrizetheexpres-
α α
similarly for Π . sions based on energy and charge conservation.
b
Atthesecondstep,weadiabaticallyswitchontheEPI The set of coupled equations Eqs. (17-22) is difficult
inthecenter. Weperformaperturbationexpansionover to solve, due to the many-body EPI. Since the EPI is
5
weak, we consider only the lowest order Feynman dia- constant level-width broadening Γ with energy cutoff
α
grams shown in Fig. 2. The expressions for the self- ε (see Eq. (40) for the general definition). It interacts
D
energies are as follows. The electron Fock self-energies with an isolated phonon mode with frequency ω , and
0
from phonons are H = m d†du. In the linear regime, we introduce an
ep 0
infinitesimal change of the chemical potential or temper-
(cid:90)
dω
ΣFm,n<,>(ε)=i(cid:126) MmkiG<0,,i>j(ε−)D0<,k,>l(ω)Mjln2π, σatu=reeaotrlepa,dµLa,ned.gT., µaLre=thµe+coδrµre,sTpσoLn=dinTg+eqδuTilσib,rwiuitmh
(27) values. We look at the response of the charge and heat
currentduetothissmallperturbation. Theresult, upto
(cid:90) the 2nd order in M, is summarized as follows,
ΣF,r(ε)=i(cid:126) Mk (cid:0)Gr (ε )D< (ω)
mn mi 0,ij − 0,kl (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
+G>0,ij(ε−)D0r,kl(ω)(cid:1)Mjlnd2ωπ. (28) IIehe = LL01 LL12 δδTTµe . (32)
The linear conductance and the Seebeck coefficient are
The Hartree self-energy does not depend on energy
(cid:90) G =e2L , (33)
dε e 0
Σrmn =−iMmi nDirj(ω =0) MkjlG<lk(ε)2π. (29) S ≡−δV =− L1 . (34)
δT eL T
0
The phonon self-energies from electrons are
The coefficients L are
n
(cid:90)
Π<,>(ω)=−i dεTr(cid:2)MmG<,>(ε)MnG>,<(ε )(cid:3), 3
mn 2π 0 0 − L =(cid:88)L(i), (35)
(30) n n
i=1
(cid:90)
Πr,a(ω)=−i dεTr(cid:2)MmGr,a(ε)MnG<(ε ) with
mn 2π 0 0 − 1 (cid:90) dε
+MmG<0(ε)MnGa0,r(ε−)(cid:3), (31) L(n1) = (cid:126) 2π(ε−µ)n(Aα¯Γα)f(cid:48), (36)
withε− =ε−(cid:126)ω. Summationoverrepeatedindicesisas- L(2) = 1 (cid:90) dε(ε−µ)n(∆A(cid:48) Γ )f(cid:48), (37)
sumed here. Different charge and energy conserving ap- n (cid:126) 2π α¯ α
proximations have been developed in the literature. We L(3) = 1 (cid:90) dε(ε−µ)n(cid:0)∆A(cid:48)(cid:48)Γ +2GrImΣr GaΓ (cid:1)f(cid:48).
will use two of them. In Subsec. IIIB, we perform an n (cid:126) 2π α¯ α 0 ep 0 α
expansion of the current up to the second order in M,
(38)
followingtheideaofRef.29. Inthenumericalmodelcal-
culation in Subsec. IIIC, we use the self-consistent Born We have defined
approximation (SCBA), which means we replace G , D
0 0 ∂f
by G, D respectively in the above equations. f(cid:48) =− , (39)
∂ε
Γ =i(Σr −Σa), (40)
α α α
A =GrΓeGa, (41)
α 0 α 0
∆A(cid:48) =GrReΣr A +A ReΣa Ga, (42)
α¯ 0 ep α¯ α¯ ep 0
∆A(cid:48)(cid:48) =iGrImΣr A +iA ImΣa Ga, (43)
α¯ 0 ep α¯ α¯ ep 0
and α¯ means the lead different from α. L(1) is the single
n
electron Landauer result. L(2) is the quasi-elastic term.
n
L(3) is the inelastic term. f is the Fermi-Dirac distribu-
tion function
FIG. 2. Feynman diagrams due to electron-phonon interac-
tion. ThefirsttwoareHartreeandFockdiagramforelectrons, (cid:20) (cid:18)ε−µ (cid:19) (cid:21)−1
and the last one is the polarization bubble for phonons. The fα(ε)= exp k Tαα +1 . (44)
expressions of these diagrams can be found in Eqs. (27-31). B
Since we are looking at the linear response regime, f =
L
f , we dropped the subscript in Eq. (39). We will also
R
use the Bose-Einstein distribution later
B. Thermoelectric transport through a single (cid:20) (cid:18) (cid:126)ω (cid:19) (cid:21)−1
electronic level n (ω,T)= exp −1 . (45)
B
k T
B
We consider a single electronic level HC =ε d†d, cou- When there is no ambiguity, we will also drop the argu-
e 0
pledtotheleftandrightelectrodes,characterizedbythe ment T.
6
In the following, we set the position of the electronic general, it depends on the system parameters. At non-
level to ε = 0, and look at the dependence of the con- zerotemperature, thesharpthresholdbroadensout, and
0
ductance on the chemical potential µ. We firstly write the linear conductance is affected: its correction to the
downtheexpressionsfortheself-energies,andmakesome conductance ∆G is negative for µ ∼ 0, and positive
eFI
observations based on their functional forms. for µ(cid:29)Γ [Fig. 3 (c)].
TheHartreeself-energyisrealanddoesnotdependon Physically, ImΣr gives rise to phonon scattering pro-
energy32 cesses. ItseffectcaFnbeunderstoodasfollows: AtT =0,
Σr =−(cid:88)m20Γα (cid:90) εD fα(ε) dε. (46) fbolresdmuealtlobPiaasu(lieVblo<ck(cid:126)inωg0,),wphhiloenpohnoenmonissaidosnorisptnioont pisonssoit-
H α 2πω02 −εD ε2+Γ2/4 possible due to zero phonon population. So, ImΣrF does
not affect the linear conductance. At high enough tem-
At T =0, we get
perature,bothphononemissionandadsorptionarepossi-
Σr =−(cid:88)m20Γα tan−12ε(cid:12)(cid:12)(cid:12)µα , (47) bdliestervibeuntaiotns,maanldlbfiinaist,edpuoeptuolatthieonbroofapdheonninognomfotdheesF.eTrmhei
H πω2Γ Γ(cid:12)
α 0 −εD phonon scattering process decreases the conductance on
resonance, but increases it far off resonance. As a result,
with Γ=Γ +Γ . For large enough ε , the lower limit
L R D
term turns to −m2/(2ω2), which is the polaron energy the Seebeck coefficient becomes smaller.
0 0
shift. Note here that the 1/2 is due to the fact we use TherealpartReΣr(ε)isobtainedbytheHilberttrans-
F
uC in our definition of H (Eq. 10). This is different form of the imaginary part. At zero temperature, it di-
frkom the common definitioepn that uses the creation and verges logarithmically at ε−µ = ±(cid:126)ω077. Its effects on
annihilation operatore a† + a (Eq. 74). We have sub- the conductance and Seebeck coefficient are difficult to
tracted the polaron shift term in the following calcula- analyze. We rely on the numerical result [Fig. 3 (b)].
tion, since it is a constant. After this subtraction, Σr Figure3showsthecorrectiontothelinearconductance
H
is odd in µ with a negative slope near µ = 0. We fo- of different self-energy terms as a function of µ. These
cus on the µ > 0 regime. It saturates to −m2/(2ω2) for numerical results confirm our qualitative analysis. By
0 0
large µ, e.g., µ(cid:29)Γ. At non-zero temperature, the slope comparing the total conductance at low [Fig. 3 (e)] and
and the saturation value change, but the shape of the hightemperature[Fig.3(f)],weseethat,(1)theHartree
curve is similar to the T =0 case. This means that, the term dominates at low temperature, and the G -µ peak
e
Hartree term shifts the electronic level, and reduces the becomesnarrower. (2)theFocktermbecomesimportant
conductance. On the other hand, when µ(cid:29)Γ, the con- at high temperature, and the G -µ peak broadens out.
e
ductance tends to zero, whether we include the EPI or Their effects on the Seebeck coefficient (S) are shown in
not. Thus, the correction to the conductance due to the Fig. 4. At low temperature, when the EPI is included,
Hartree term ∆GeH ≤ 0. It starts from zero at µ = 0, the magnitude of S gets larger for µ ∼ 0, and smaller
goes back to zero at µ (cid:29) Γ, and it reaches a maximum for |µ|(cid:29)Γ. At high temperature, the effect of the Fock
magnitude at some point in the middle. The described term results in drop of S. In any case, the correction to
behaviourisschematicallyshowninFig.3(a). Thiseffec- S is small for weak EPI. But for the case of strong EPI,
tively reduces the broadening of the single level spectral the correction could be large (see Subsec. IVC).
function, also the conductance peak. Since the Seebeck
coefficient is related to the logarithmic derivative of the
conductance, we expect it to increase the magnitude of
the Seebeck coefficient near resonance, and to reduce it C. Heat transport between electrons and phonons
off resonance [Fig. 4].
The imaginary part of the retarded Fock self-energy
Let us look back at the setup in Fig. 1. We want to
is32
study the heat transport between electrons and phonons
m2 (cid:88) at finite temperature bias, but zero voltage bias. The
ImΣrF(ε)=−4ω0 sAα(ε−s(cid:126)ω0) simplest setup is that the system couples to one electron
0
α,s=± and one phonon lead, each at its own temperature, see
×(cid:2)1+n (sω )−f (ε−s(cid:126)ω )(cid:3). (48) Fig. 5. The expression for the energy current from elec-
B 0 α 0
trons to phonons can be obtained from Eq. (26) and the
It is negative and even in ε. Its role on the differential expressions of the self-energies Eqs. (27-31)32
conductance at the phonon threshold (eV = (cid:126)ω ) has
0
been discussed extensively74–77. The main conclusions (cid:90) (cid:90)
are: it reduces the differential conductance at eV =(cid:126)ω Q=i(cid:126) dε dωω(cid:2)G> (ε)Mk D<(ω)G<(ε )Ml (cid:3),
0 2π 2π nm mi kl ij − jn
for resonant case (µ ∼ 0), where the bare transmission
without EPI (T ∼ 1), while it does the opposite for far (49)
0
offresonancecase(µ(cid:29)Γ),whereT →0. Thetransition
0
point between the two opposite behaviors is T = 1/2 if where, again, summation over repeated indices is as-
0
the electronic density of state (DOS) is flat. But, in sumed. Fortheeaseofanalysis,weperformanexpansion
7
0.00 0.06
(cid:3)(cid:2)(cid:3)(cid:3)GGeH0(cid:2)(cid:2)(cid:2)(cid:2)(cid:2)(cid:2)00000000........0101000120864220 (a) (cid:2)(cid:3)(cid:3)GGGeFR0(cid:2)0000....00000242(b) εε μL,tTL εε t εε HKtep εε K TTKR
(cid:2)4 (cid:2)2 0 2 4 (cid:2)0.04 (cid:2)4 (cid:2)2 0 2 4
ΜΜ(cid:2)(cid:2)(cid:2)(cid:2)ΩΩ00(cid:3)(cid:3) Μ(cid:2)(cid:2)(cid:2)(cid:2)Ω00(cid:3)(cid:3) FIG. 5. The model system we consider to study the energy
0.00 0.00 transport between electrons and phonons.
(cid:2)(cid:3)(cid:3)GGFI0eFI(cid:2)(cid:2)(cid:2)(cid:2)(cid:2)00000.....0000014322 (c) (cid:2)(cid:3)(cid:3)(cid:2)GGGe0(cid:2)(cid:2)(cid:2)(cid:2)00000.....0000042864 (d) where (cid:90) dε (cid:2) (cid:3)
(cid:2)44 (cid:2)22 00 22 44 (cid:2)44 (cid:2)22 00 22 44 Λ(ω,Te)= 2πTr MA(ε)MA(ε−)
1.0 Μ(cid:2)(cid:2)Ω0(cid:3) 1.0 Μ(cid:2)(cid:2)Ω0(cid:3) ×(cid:2)f(ε,Te)−f(ε−,Te)(cid:3), (51)
00..88
00..88
and
(cid:2)(cid:3)GGe00000....2462 ((ee)) (cid:2)(cid:3)GGe00000....2462 ((ff)) A(ω)=i(cid:2)D0r(ε)−D0a(ε)(cid:3). (52)
0.0 (cid:2)4 (cid:2)2 0 2 4 0.0 (cid:2)4 (cid:2)2 0 2 4 Now the question we ask is whether there is a diode
Μ(cid:2)(cid:2)Ω(cid:3) behaviour for the heat transport between electrons and
0 Μ(cid:2)(cid:2)Ω(cid:3)
0 phonons52, e.g., Q(∆T)(cid:54)=Q(−∆T), with∆T =T −T .
e p
This is relevant because in some special situation, e.g.,
atmetal-insulatorinterface,orinsulatingmolecularjunc-
FIG. 3. Chemical potential dependence of the conductance
correctionduetotheΣr (a),ReΣr(b),ImΣr(c),andthesum tions, EPI becomes the bottleneck of heat transfer78–82.
H F F
ofthem(d). Theblue,red,purple,andblacklinescorrespond We can define the rectification ratio as
to k T = 0.05,0.15,0.25,0.35(cid:126)ω , respectively. (e)-(f) The
B 0 Q(∆T)+Q(−∆T)
electricalconductanceasafunctionofchemicalpotential(µ) R= . (53)
at k T = 0.05(cid:126)ω (e) and k T = 0.35(cid:126)ω (f). Solid lines Q(∆T)−Q(−∆T)
B 0 B 0
includeEPI,whilethedashedlinesdonot. Otherparameters
used: Γ =Γ =(cid:126)ω =1, m =1(cid:126)ω /(˚A√u). If we assume a constant electron DOS, Λ(ω) does not
L R 0 0 0
depend on T, we get Q(−∆T) = −Q(∆T), and R = 0.
The physical reason is that in this case, it is possible
to map the electron-hole pair excitation into harmonic
1.0
oscillators79,83,84. Then, it is equivalent to heat trans-
port within a two-terminal harmonic system. We do not
0.5
expect any rectification effect. To make R(cid:54)=0, the elec-
L
e tronicDOShastobeenergy-dependentwithinthebroad-
(cid:144)
B 0.0 ening of the Fermi-Dirac distribution given by k ∆T.
k B
H This effectively introduces anharmonicity into the sys-
S -0.5 tem, consistent with previous studies79–81.
We can go one step further, by making a Taylor ex-
-1.0 pansion of the spectral function A(ε) about the Fermi
-4 -2 0 2 4 energy, we find that the sign of R is determined by the
sign of ∂2A (The 1st order term is zero).
m (cid:209)w0 To ch∂eεc2k this argument, we calculated the heat cur-
rent across one-dimensional (1D) metal-insulator junc-
FIG. 4. Change of the Seebeck coefficient due to EPI at dif- tion. Themetalsideisrepresentedbya1Dtight-binding
ferent temperatures. The parameters and meaning of colors chain, with hopping element t = −0.1 eV, and onsite
are the same as Fig. 3. energy ε = 0. The insulator side is represented by a
0
H L
1D harmonic chain with the spring constant K = 0.1
eV/(˚A2u). Theinsulatorandmetalcouplethroughtheir
last two degrees of freedom. Their interaction matrices
oftheaboveexpressionto2ndorderinM,anditbecomes are
(cid:18) (cid:19)
0 1
Mk =(−1)km , k =1,2. (54)
Q(2) =(cid:126)(cid:90)0+∞ d2ωπωTr(cid:2)Λ(ω,Te)A(ω)(cid:3) Here, m0 = 0.05 eV/(˚A√0 u)1, k0represents the phononic
(cid:2) (cid:3)
× n (ω,T )−n (ω,T ) , (50) degrees of freedom. This means that the system couples
B e B p
8
mechanicalexpansiontechnique,whichrecordstheJoule
((a)) ((bb)) heat induced temperature rise85. The Joule heat can in-
crease the temperature of the molecular junction from
room temperature to 463 K, which has been examined
through the inelastic electron tunneling spectroscopy86.
Grosse et al. investigated the nanoscale Joule heating
in phase change memory devices85. The Joule heating
p e
leads to the temperature rise in the phase change mem-
e p
ory device, which results in an obvious volume expan-
sion. In another experiment, Joule heating is found to
ΔΔTT((KK)) μμ((eV)) be responsible for the correlated breakdown of nanotube
(d) forests39,41.
(c)
For a system without localized phonon modes, all
phonon modes have important contribution to the Joule
heat. The Joule heat contributed by these propagating
phonon modes have important effects on the electric de-
vices. For example, in graphene transistors, the output
electriccurrentwillsaturatewithincreasingsource-drain
voltage87,88. This saturated current density can be re-
duced by 16.5% due to the Joule heating88.
μ(eV) Energy(eV)
Thelocalizedphononmodesexistaroundsomedefects
or nonuniform configurations, such as the free edge, the
FIG. 6. (a) The energy current with µ = 0 (red, solid) and isotropic doping, interface, etc. This particular type of
µ=0.1eV(blue,dotted). WehavesetT =T+∆T/2,T =
e p phonon modes has no direct contribution to the thermal
T −∆T/2, with T = 300 K. (b) Fermi level (µ) dependence
conduction, but localized phonon modes play a particu-
of the energy current at fixed temperature difference ∆T =
larlyimportantroleintheJouleheatphenomenon. They
300 K. (c) Rectification ratio R as a function of µ. (d) The
are characteristic for their exponentially decaying vibra-
electronic spectral function A(ε) as a function of energy ε.
tion displacement; i.e., only a small portion of atoms
are involved in the localized vibration. For instance,
therearesomelocalizededgephononmodesatgraphene
to only one electron and one phonon lead [Fig. 5]. Fig-
nanoribbon’s free edge. In these modes, edge atoms vi-
ure 6 summarizes our result. In Fig. 6 (a), we show the
brate with large amplitude, but the vibrational displace-
energycurrentasafunctionoftemperaturedifferencebe-
mentdecaysexponentiallyfromtheedgeintotheinterior
tween the metal and insulator (∆T), at two Fermi levels
region.
(µ = 0,0.1 eV). The energy current is asymmetric with
The localized-phonon-mode-induced Joule heat was
respect to the sign change of ∆T. But the rectification
observed in graphene nanoribbons in experiment, and
ratio R has opposite sign. This is further highlighted in
explained theoretically. Jia, et al. utilized Joule heat-
theplotoftheenergycurrentandtherectificationratioR
ingtotriggertheedgereconstructionatthefreeedgesin
asafunctionofµforfixedtemperaturebias∆T =±300
thegraphenenanoribbons89. Engelund,etal. attributed
K [Fig. 6 (b) and (c), respectively]. The sign of R is well
correlated with the sign of ∂2A [Fig. 6 (d)]. thisphenomenontotheJouleheatingoftheedgephonon
∂ε2 modes90. There are two conditions for the important
Joule heating of the edge phonon modes. First, these lo-
calized edge phonon modes can spatially confine the en-
D. Effect of EPI on thermal transport in ergyattheedges. Second,theelectronsinteractstrongly
single-walled carbon nanotubes withthelocalizededgephononmodes. Themeansteady-
state occupation of the edge phonon mode can be calcu-
Phonon modes can be excited by the mobile electrons latedfromtheratioofthecurrent-inducedphononemis-
due to the EPI effect; i.e., a high bias over the system sion rate and damping rate. The effective temperature
leadstoself-heating(Jouleheating). Innanoscaleelectric forthefreeedgecanbeextractedbyassumingthisoccu-
devices, the electric current density can be much larger pationtobeBosedistributed. Theeffectivetemperature
than that in the macroscopic system. The high current wasfoundtobeashighas2500Kforbiasaround0.55V.
density will generate strong Joule heating, which may This high effective temperature was proposed to be the
eventually break the device. In this sense, Joule heating origin for the edge reconstruction.
becomes a bottleneck for further increase of the electric Although Joule heating might be used for selectively
currentdensity. Hence,lotsoftheoreticalandexperimen- bond-breaking91,92,itsmostcommonoutcomeisadisas-
tal efforts have been devoted to understanding the Joule ter of device breakdown. The effective temperature is a
heatingphenomenoninthenanoscaleelectricdevices. In suitable quantity to describe the Joule heating. In 1998,
experiment,Jouleheatcanbemeasuredviathethermal- Todorov studied Joule heating problem in a molecular
9
K) 1
W/ 0 T=150K
n (a)
(h−1 ballistic
σp−2 epi
3 0.2 0.8 1.4 2
K)
W/ (b)
(nh 2 T=300K
p
σ
1
K) 6.5
W/
n (c)
(h 6
p T=1000K
σ
5.5
0.2 0.8 1.4 2
µ (eV)
FIG.8. (Coloronline)Thephononthermalconductancever-
sus chemical potential for semiconductor SWCNT (10, 0) at
(a)150K,(b)300K,and(c)1000K.Solidlineisforballistic
phonon thermal conductance without EPI effect. Reprinted
from J. Appl. Phys., 110, 124319 (2011).
FIG.7. (Coloronline)Thephononthermalconductancever-
effect on the thermal conductance in single-walled car-
sus chemical potential for metallic SWCNT (10, 10) at (a)
bon nanotubes (SWCNTs)94. For them, we apply the
150 K and (b) 300 K. Solid line is for ballistic phonon ther-
malconductancewithoutEPIeffect. ReprintedfromJ.Appl. Born approximation to consider the EPI effect using the
Phys., 110, 124319 (2011). NEGF approach, as the SCBA is computationally more
expensive. Thephononthermalcurrentcanbecalculated
by considering the three EPI contributions shown in the
Feynman diagrams in Fig. 2. The phonon thermal cur-
junction12. In his work, the Einstein model is applied
rent flowing from the left lead into the center is given by
to represent the phonon modes in the system, and the
Eq. (26). The expression for the right lead is analogous.
electron-electroninteractionisignoredasthesystemsize
The Joule heat is generated in the system and flows into
is much smaller than the electron mean free path. Part
the leads, so the total Joule heat is the sum of heat cur-
of the EPI-induced Joule heat will be delivered out of
rentsintobothleads,Q=−(IL+IR). Thethermalcur-
the system by the phonon heat conduction, while the re- p p
rentfromEq.(26)alsoincludesthatinducedbythetem-
maining Joule heat gives a high effective temperature.
perature gradient, which satisfies IL = −IR. Hence, Q
For low ambient temperature, the effective temperature p p
givessolelytheJouleheat.FormetallicSWCNT(10,10),
scales with voltage V as T4 ≈γ4V2, with γ as an EPI-
eff both electrons and phonons are important heat carriers.
dependentconstant. Itwasshownthattheeffectivetem-
The EPI only slightly reduces the electron thermal con-
peraturecanbeabove200Kforaverylowambienttem-
ductance, but it has a strong effect on the phonon ther-
peraturearound4K93. Butatveryhighbias,thescaling
mal conductance. More specifically, Fig. 7 (a) shows an
law could differ from this31.
‘electron-drag’effectonthephononthermalconductance
There are several experimental approaches to inves- at150K.Thephononthermalconductancebecomesneg-
tigate the effective temperature of the electric device ativeforhighchemicalpotentialvalueµ>2.0eV,which
induced by Joule heating. The effective temperature indicates that electrons can help to drag phonons from
can be extracted by measuring some quantities that are cold temperature region to the hot temperature region.
temperature-dependent. Forexample,thebreakingforce The‘electron-drag’phenomenonhappensatlowtemper-
ofthesinglemolecularjunctionisrelatedtothetempera- ature and high chemical potential, and it does not hap-
ture. This force-temperature relationship can be used to penatahighertemperature300KasshowninFig.7(b).
estimate the effective temperature30. The Raman spec- For semiconductor SWCNT (10, 0), the electronic ther-
troscopyalsodependsonthetemperature. Hence, itcan mal conductance contributes less than 10% of the to-
be used to deduce the effective temperature of Raman- tal thermal conductance at low bias (e.g. µ = 0.3 eV),
active phonon modes35,37,38. while phonons make most significant contribution to the
It has also been shown that EPI has an important total thermal conductance. Similar ‘electron-drag’ phe-
10
nomenon also exists in the semiconductor SWCNT (10,
0) at low temperature as shown in Fig. 8 (a).
∂ρ˜I(t) = −1(cid:88)(cid:90) tdt(cid:48)[Sα(t),Sβ(t(cid:48))ρ˜ (t )]Cαβ(t−t(cid:48))+H.c.
∂t (cid:126)2 I I I 0
α,β t0
IV. STRONG EPI REGIME: QUANTUM (59)
MASTER EQUATION APPROACH where Cαβ(t−t(cid:48)) = Tr[ρL⊗ρRBα(t)Bβ(t(cid:48))] is the cor-
e e I I
relation function of the leads. Here we have used the
A. Quantum Master equation formulism condition that the expectation value of a single lead op-
erator Bα is zero. We can now transform back to the
In this section, we introduce the QME approach to Schr¨odinger pictureandextendtheinitialtimet0 to−∞
consider the case of strong EPI. Before doing that, we to get the QME of Redfield type101–103
should mention that the NEGF method has also been
∂ρ˜ i 1 (cid:88)
usedtotreatthestrongEPI95–100. Sincetheideabehind ∂t =−(cid:126)[HS,ρ˜]− (cid:126)2 (60)
itisverysimilartothatofthemasterequationapproach, α,β
we choose not to introduce it here. (cid:90) t
To simplify the formula, we ignore the coupling of × dt(cid:48)(cid:8)[Sα,Sβ(t(cid:48)−t)ρ˜]Cαβ(t−t(cid:48))+H.c.(cid:9).
molecularphononmodestothephononleads. Themodel −∞
Hamiltonian simplifies to Here we have replaced ρ˜(t ) by ρ˜, which is essential
0
and correct only when one intends to get the 0th-
H =H +HL+HR+V (55) order reduced density matrix ρ˜0 by solving the above
tot S e e e
QME69,104,105. In the application to the EPI problem,
where H =HC+HC+H denotes the system Hamil- byexactdiagonalizingthesystemHamiltonian,thisRed-
S p e ep
tonian, and V = VL+VR is the system-lead coupling. field QME can take into account the coherence between
e e e
In the QME formalism we assume the system-lead cou- electronsandphonons,incontrasttotheusualrateequa-
pling V is weak so we can do perturbation on it. We tion approach106,107.
e
work in the interaction picture with H = H − V as We write the above equation in the eigenbasis of the
0 e
non-interacting part and Ve as the interaction. For sim- system Hamiltonian HS to obtain102,108
plicity, in this section we use V to represent V since we
mdoant’rtixhafvoelloVwp.sTthheeveoqnuaNtieounmoafnmnoetqiounatfoiornthefuelldensity dρ˜dntm =−(cid:126)i∆nmρ˜nm+(cid:88)Rnijmρ˜ij, (61)
ij
∂ρ (t) where the relaxation tensor reads109
i(cid:126) I =[V (t),ρ (t)]. (56)
I I
∂t 1 (cid:88)(cid:110)
Rij = SαSβ Wαβ
nm (cid:126)2 ni jm ni
Here, the subscript I denotes operator in the interaction
α,β
picture. The time argument in the parentheses means (cid:88) (cid:111)
non-interacting evolution O(t) = eiH0t/(cid:126)Oe−iH0t/(cid:126). The −δjm SnαlSlβiWlαiβ +H.c. (62)
above equation can be written in an integral form as l
The transition coefficients are given by
−i(cid:90) t
ρI(t)= (cid:126) t0 dt(cid:48)[VI(t(cid:48)),ρI(t(cid:48))]+ρI(t0). (57) Wkαjβ =(cid:90) t dt(cid:48)ei∆kj(t(cid:48)−t)/(cid:126)Cαβ(t−t(cid:48)), (63)
−∞
Onecanrecursivelyapplytheaboveequationtogetase-
where ∆ = E − E are the energy spacings of the
ries expansion of the full density matrix in power of V . kj k j
I
system Hamiltonian.
Wetruncatetheseriestothesecondorderanddifferenti-
Since we are only interested in the steady state, we
ateitwithrespecttotimetatbothsidesoftheequation
impose the condition dρ˜/dt = 0 at t = 0 and solve the
to get the following integro-differential equation
above equation order by order with respect to V. One
can find that all the off-diagonal elements of the 0th-
∂ρ (t) −i 1 (cid:90) t order reduced density matrix vanish in steady state and
∂It ≈ (cid:126) [VI(t),ρI(t0)]−(cid:126)2 dt(cid:48)[VI(t),[VI(t(cid:48)),ρI(t0)]]. the diagonal elements can be evaluated via the matrix
t0 (58) equation109(cid:80)iRniinρ˜(ii0) =0,togetherwiththeconstraint
We prepare the initial state as a product state of the of Tr[ρ˜(0)]=1.
system and each lead, ρ(t ) = ρ˜(t )⊗ρL⊗ρR. For the For the calculation of currents, we go through a sim-
0 0 e e
system-lead coupling V we assume it can be written as ilar derivation as the QME. The electronic current op-
a product of system operator S and lead operator B as erator Je and heat current operator Jh can be writ-
V = (cid:80)αSα ⊗Bα. In such cases we can trace over the ten in the form Je(h) = (cid:80)αSα ⊗ Beα(h). The expec-
lead degrees of freedom to get tation value of currents can be calculated according to