PHYSICAL REVIEW B 86, 125443 (2012)

Electron-phonon mediated heat flow in disordered graphene
Wei Chen and Aashish A. Clerk
Department of Physics, McGill University, Montreal, Canada H3A 2T8(Received 11 July 2012; published 26 September 2012)
We calculate the heat ﬂux and electron-phonon thermal conductance in a disordered graphene sheet, going
beyond a Fermi’s golden rule approach to fully account for the modiﬁcation of the electron-phonon interaction
by disorder. Using the Keldysh technique combined with standard impurity averaging methods in the regime
kF l  1 (where kF is the Fermi wave vector and l is the mean free path), we consider both scalar potential(i.e., deformation potential) and vector-potential couplings between electrons and phonons. We also consider theeffects of electronic screening at the Thomas-Fermi level. We ﬁnd that the temperature dependence of the heatﬂux and thermal conductance is sensitive to the presence of disorder and screening, and reﬂects the underlying
chiral nature of electrons in graphene and the corresponding modiﬁcation of their diffusive behavior. In the case of
weak screening, disorder enhances the low-temperature heat ﬂux over the clean system (changing the associatedpower law from T 4 to T 3 ), and the deformation potential dominates. For strong screening, both the deformation
potential and vector-potential couplings make comparable contributions, and the low-temperature heat ﬂux obeys
a T 5 power law.
DOI: 10.1103/PhysRevB.86.125443

PACS number(s): 65.80.Ck, 72.10.Di, 44.10.+i

I. INTRODUCTION

The potential to exploit the exceptional thermal propertiesof graphene in applications has recently generated considerable activity:1–4 possible applications include sensitive
bolometry and calorimetry for detecting infrared and terahertz
radiation. Such detectors would ultimately be based on the simple heating of electrons in a graphene sheet by the absorptionof incident photons. Ideally, the detector electrons would bethermally decoupled from their surroundings, thus allowing
any heating produced by the incident radiation to be long
lived. In this respect, graphene provides a potential advantage:
Its low electron density and relatively weak electron-phonon
coupling implies that the expected low-temperature thermal
decoupling between electrons and the lattice5,6 could occur
over a much wider temperature range than in a conventional
metal.7,8 Further, one can effectively suppress the thermal link
between graphene electrons and electrons in the contacts byemploying superconducting leads.9
Given the above, it is crucial to develop a rigorous andquantitative understanding of the thermal link between electrons and phonons in graphene at low temperatures. Several
recent theoretical works have addressed this problem in thecase of clean graphene (i.e., no electronic disorder).7,8,10,11 At
low temperatures, the heat ﬂux (per volume) between electrons
and longitudinal acoustic phonons takes the general formδ,
P (Te ,Tph ) = F (Te ) − F (Tph ) =  Teδ − Tph(1)where F (T ) is called the energy control function,  is a
coupling constant, and Te and Tph are the temperatures ofthe electrons and lattice (i.e., phonons), respectively. Previous
theoretical works7,8,10,11 ﬁnd that δ = 4 in the low-temperature
limit (assuming an unscreened deformation potential electronphonon coupling); this has been conﬁrmed in recent experiments on graphene in the clean limit.2–4 This power law is also
identical to what would be expected for a clean conventional
two-dimensional (2D) metal.11
1098-0121/2012/86(12)/125443(14)In this work, we now ask how the above result is modiﬁed in
the presence of electronic disorder. While great experimental
progress has been made in reducing disorder effects ingraphene,12 the devices studied for bolometric applications in
Refs. 2–4 are sitting on a silicon substrate and have mean freepaths that are 100 nm or less. In conventional metals, electronic
disorder can strongly affect the electron-phonon coupling attemperatures low enough that the wavelength of a thermal
phonon is comparable to (or longer than) the electronic mean
free path.13–16 This deﬁnes a characteristic temperature scale
Tdis below which disorder effects are important,kB Tdis ≡ hs/ l,(2)where s is the speed of sound and l is the electronic mean freepath. As discussed in Refs. 13–16 (and below), the effects ofdisorder are subtle: Depending on the nature of the disorderand the electronic system, the power law δ in Eq. (1) can either
be enhanced by disorder or be suppressed.13–16
Here, we study how the additional complexity arisingfrom the unique electronic properties of graphene modify
the interplay of disorder and the electron-phonon interaction.
The principle new ingredients arise from the effective chiralnature of carriers in graphene, which both modiﬁes electronic
diffusion, and allows for a new kind of effective vectorpotential electron-phonon coupling.17–19 For simplicity, we
will focus on disorder originating from charges in the substrate
below the graphene ﬂake, and thus take the disorder potential
to preserve the symmetries (valley and sublattice) of thelow-energy graphene Hamiltonian.20,21 The impurity potential
is thus also taken to be static, i.e., it does not move withthe graphene sheet. We also consider the case where thegraphene ﬂake has been doped sufﬁciently that kF l  1 (kF is
the Fermi wave vector), meaning that the standard impurityaveraged perturbation theory is appropriate.22,23 Combining
this approach with the Keldysh technique then allows us
to rigorously address the electron-phonon interaction in thepresence of disorder, in a manner analogous to the classic
works looking at this physics in a conventional disordered

125443-1©2012 American Physical Society

WEI CHEN AND AASHISH A. CLERKPHYSICAL REVIEW B 86, 125443 (2012)

metal.13–16 We stress that properly addressing disorder effects
involves going beyond the sort of golden rule calculation used
to address the clean case.11
Note that since we work in the regime kF l  1, the temperature scale Tdis below which disorder effects emerge will
necessarily be well below the Bloch-Grüneisen temperature
kB TBG = 2h̄skF ; we will thus explicitly focus on temperatures
T < TBG . For a typical doped graphene electron density
1012 /cm2 , TBG ∼ 70 K. In contrast, for a typical mean free pathof 100 nm and a graphene acoustic phonon velocity of ∼2 ×
104 m/s, Tdis ∼ 10 K. We note that a recent study examined
disorder effects on electron-phonon interactions above the BG
temperature.1 Somewhat surprisingly, the expression derived
in that work for impurity-assisted electron-phonon coolingfor T  TBG is exactly half of our expression for e-phonon
cooling based on the deformation potential at T  Tdis < TBG[cf. Table I and Eq. (38)].
Our main results for the energy control function F (T ) [cf.
Eq. (1)]) for graphene at low temperatures are summarized inTable I. We consider the contribution to the electron-phonon
heat ﬂux arising from both the standard deformation-potential(DP) coupling, as well as from the effective vector-potential(VP) coupling . In the absence of disorder (or at temperatures
well above Tdis but below TBG ), and in the absence of electronic
screening, one ﬁnds that both these mechanisms contributeindependently and in a similar manner: The respective heat
ﬂuxes are each described by Eq. (1) with δ = 4 (in agreementwith Refs. 10 and 11).
In contrast, for T < Tdis , the two coupling mechanisms areaffected oppositely by the electronic disorder. We ﬁnd thatthe heat ﬂux associated with the VP coupling is suppressedby disorder: Still neglecting screening, it is described now by
Eq. (1) with an enhanced power law of δ = 5. Heuristically,
this is attributed to the disorder broadening of the graphene
energy levels. The effect of disorder on the heat ﬂux associatedwith the DP coupling for T < Tdis is the opposite from theabove: it is enhanced. It is described by a reduced power law
δ = 3 (again, no screening), and will thus dominate the VP at
low temperatures. On a heuristic level, this enhancement is dueto the diffusive charge dynamics, which effectively increases
the time an electron interacts with a given phonon (i.e., thisbecomes the time to diffuse across a phonon wavelength,
as opposed to the time needed to ballistically traverse thisdistance). The absence of any diffusive enhancement of the
VP coupling is (as we will show) a direct consequence of thenonconservation of pseudospin.
We also consider how including screening changes theabove results; the importance of screening the e-ph interaction
has been the subject of several recent studies.24–26 As discussed
extensively in Refs. 24 and 26, the DP coupling is expected
to be screened, whereas the VP coupling is expected to be
unscreened, as it induces no net electronic charge (i.e., theeffective vector potential generated by a phonon ﬁeld has
opposite sign in the two graphene valleys). As a result (seeTable I), even without disorder (i.e., Tdis < T < TBG ), the DP
and VP heat ﬂuxes are not equivalent in the limit of strong
screening: The VP power law remains δ = 4, where the DP
power law is increased to δ = 6. Similarly, at low temperatures
where disorder effects matter, the VP power law is unchanged,
but the DP power law becomes T 5 .

Note that our results suggest that even though the bare
VP coupling strength g2 is believed to be about an order ofmagnitude smaller than the bare DP coupling g1 ,18 if screening
is strong, its contribution to the heat ﬂux could be comparable
to or even larger than that from the DP coupling. This is
despite the relative enhancement of the DP coupling over the
VP coupling by disorder. Further, our results suggest that theelectron-phonon heat ﬂux could be a means for empiricallydetermining if screening is important. In particular, the only
way to obtain a T 3 power law is via an unscreened DP couplingin the diffusive limit. We note that measurements of the phonon
contribution to the resitivity in a clean graphene sheet (asrecently measured27 ) cannot directly resolve this issue, as bothDP and VP couplings contribute to a T 4 dependence.26
This paper is organized as follows. In Sec. II, we present
our model and an outline of the calculational method. This
includes a brief derivation of the kinetic equation of electrons
in graphene in the Keldysh formalism (Sec. II B), as well
as a derivation and discussion of the diffusion propagator ingraphene and the resulting diffusive renormalization of theelectron-phonon vertex (Sec. II C). We present the main resultsin Sec. III, i.e., the heat ﬂux due to electron-phonon interaction
in both weak and strong screening cases. Finally, we brieﬂysummarize the paper in Sec. IV.
II. MODEL AND CALCULATION
A. Model Hamiltonians of impurity andelectron-phonon scattering
1. Electrons in disordered graphene

The low-energy electronic degrees of freedom are describedby a massless Dirac Hamiltonian. Focusing on a single valley(the K+ valley), one has12,19,28,29
H = d 2 r  † (r)[−ivF σ̂j ∂j + U (r)1̂](r),(3)) is awhere vF = 106 m/s is the Fermi velocity, (r) = ( ψψBA (r)(r)
spinor ﬁeld operator describing the amplitude of electrons on
the two sublattices, σ̂j (j = x,y) are Pauli matrices, and U (r)is the disorder potential; we also set h̄ = 1 throughout unlessotherwise indicated. We do not include an index for spin orvalley, as for the physics we consider, each spin and valley
contributes in an equal and independent fashion.
As mentioned, we focus on a smooth disorder potential
originating with impurities in the substrate. We thus treat the
impurity potential U (r) as a scalar potential with respect toboth the valley degree of freedom and the sublattice degree of
freedom (i.e., pseudospin).20,21 In the standard way, U will betreated as δ-correlated Gaussian disorder, with zero mean andcorrelator:U (r)U (r ) = wδ(r − r ).

(4)

The corresponding scattering rate is 1/τ ≡ vF / l = π νw,
where ν = kF /2π vF is the density of states at the Fermi
energy per spin per valley. We note here that for a given
disorder strength w the electronic mean free path is inversely
proportional to the Fermi wave vector in graphene. As a result,
one could tune the disorder temperature scale Tdis by tuning
the electron density. Tdis increases with increasing density.

125443-2

ELECTRON-PHONON MEDIATED HEAT FLOW IN . . .

PHYSICAL REVIEW B 86, 125443 (2012)

TABLE I. Energy control function [cf. Eq. (1)] of deformation potential and vector potential e-phonon couplings in graphene below the
Bloch-Grüneisen temperature TBG = 2h̄skF /kB . g1 and g2 are the deformation potential and vector-potential coupling constant, respectively,[cf. Eqs. (6)]. EF is the Fermi energy with respect to the Dirac point, vF is the Fermi velocity, l is the electronic mean free path, and ρM is themass density of graphene per area. qTF is the Thomas-Fermi wave vector (as deﬁned in Ref. 24), and ζ (n) is the ζ function.
Deformation potential
T < Tdis
Weak screening
F (T )

vF3 ρMEF2 3
2ζ (3) g1 kB
π 2 h̄4 ls 2

T3

Strong screening
2 5
24ζ (5) g1 kB
2 l
π 2 h̄6 s 4 qTF

Vector potential

Tdis < T < TBG
Weak screening
2 4
π 2 g1 kB
15 h̄5 s 3

T5

2. Electron-phonon interaction

The electron-phonon interaction in graphene has beenstudied extensively in Res. 24 and 26. For suspended graphene,
there are both in-plane phonon modes and ﬂexural phononmodes (out of plane). In this work, we consider graphene
on a substrate (as in recent experiments probing thermal
properties2–4 ), such that ﬂexural motion is suppressed; we thus
only focus on in-plane motion. Further, at low to moderatetemperatures, the optical modes are barely excited and thedominant modes participating in cooling of hot electrons are
acoustic modes; we thus focus exclusively on the coupling tothese modes.
Due to the Dirac Hamiltonian of electrons in graphene, there
are two distinct electron-phonon coupling mechanisms.17–19
The ﬁrst is a standard deformation-potential coupling, whichcorresponds to a local dilation of the lattice. In the Dirac
theory, it appears as a scalar potential (with respect topseudospin). The second mechanism is an effective gauge-ﬁeld
coupling or vector-potential coupling. This corresponds tothe change in hopping matrix elements accompanying a pure
shear deformation, and enters the Dirac theory the same way
as an external gauge ﬁeld (the only proviso being that thisphonon-induced vector potential is valley odd, and hence does
not break time-reversal symmetry).
Letting k denote a momentum-space electronic ﬁeld
operator and bη,q a phonon annihilation operator, the totalinteraction between electrons and acoustic phonons can be
written in the general form26
 †
†
Hep =
k+q M̂ η (q) k (bη,q + bη,−q ).
(5)
η=l,t k,qHere η = l,t denote longitudinal (LA) and transverse (TA)acoustic modes, respectively. The 2 × 2 coupling matricesM̂ η (q) take the form
g1
−ig2 e2iφq
M̂ l (q) = iqξql,(6a)
ig2 e−2iφqg1
0
g2 e2iφqtt,(6b)M̂ (q) = iqξq
g2 e−2iφq0
whereξqη = (h̄/2ρM ωqη )1/2 .

(7)ηHere ρM is the mass density of the graphene sheet, ωq is
the phonon frequency for mode η, and φq is the angle ofthe phonon wave vector q with respect to the x axis (which
is taken to be along the armchair direction of the graphenelattice). For simplicity, we take the speed of sound to be theT4Strong screening
8π 4
63

6
g12 kB2h̄7 s 5 qTF

T6

T < Tdis
2 530ζ (5) g2 lkB
π 2 h̄6 s 4

T5

Tdis < T < TBG
2 4

π 2 g2 kB 4
T
15 h̄5 s 3same for the LA and TA phonons and drop the superscriptη in the frequency from now on. g1 (g2 ) is the deformationpotential (vector-potential) coupling constant. Previous works
have estimated g1 ∼ 20–30 eV and g2 ∼ 1.5 eV,18 though we
note that even the value of the deformation-potential couplingis subject to some debate.25,30,31 Our theory is thus not tiedto speciﬁc values of these parameters, and we keep both the
DP and VP couplings in our discussion. Note that transverse
phonons induce only a pure shear deformation, and hence
couple only through the vector potential.
B. Keldysh formalism for the kinetic equationof electrons in graphene

Having established the basic electronic and electronphonon Hamiltonians [cf. Eqs. (3) and (5)], we now turn to
our main goal of calculating the heat ﬂux between electronsand phonons. We consider the standard situation where eachsubsytem is independently in thermal equilibrium at its
own temperature (electrons at Te , phonons at Tph ). In thedisorder-free case, this heat ﬂux can be conveniently calculated
by using Fermi’s golden rule to calculate electron-phonon
scattering rates.7,8,10,11 Including disorder, we need a moregeneral formalism, one that is capable of capturing theinterference between electron-phonon and electron impurity
scatterings [i.e., the vertex correction of the electron-phonon
vertices M̂ η (q) by disorder]. To that end, we make use of
the Keldysh technique,32,33 coupled with standard disorderaveraged perturbation theory. Such an approach was used by
Kechedzhi et al. to study conductance ﬂuctuations (in theabsence of any electron-phonon coupling).21
We start by noting that the heat ﬂux of interest (i.e., energy
lost/gained by the electrons) can be directly related to thecollision integral I0 (ε,p) appearing in a standard Boltzmann
equation describing the dynamics of the electronic phase-space
distribution function n[ε,p; t]:
dn[ε,p].
(8)
I0 (ε,p) ≡dt
e-ph scatt
The collision integral I0 (ε,p) tells us the rate of changeof n[ε,p; t] due to the emission and absorption of acoustic
phonons.
As electronic momentum relaxation is much faster thanenergy relaxation, to describe the latter process we canfocus on times longer than the electron momentum relaxation
time; in addition, the distribution function will be sharply
peaked on shell (magnitude of momentum set by energy).
The relevant kinetics can thus be described by an electron

125443-3

WEI CHEN AND AASHISH A. CLERKPHYSICAL REVIEW B 86, 125443 (2012)

distribution function n(ε; t) that depends only on energy, not on
momentum. The kinetic equation for the electronic distributionfunction n(ε; t) takes the form∂n(ε; t)1
dp
= I¯0 (ε) =I0 (ε,p)A(ε,p),
(9)
∂t
πν
(2π )2where(10)

is the electron spectral function.
The heat ﬂux between electrons and lattice (for one valley
and spin projection) is given by
P (Te ,Tph ) = ν dε εI¯0 (ε).
(11)
The collision integral I0 (ε,p) is obtained in the standardˆmanner by calculating the electronic Keldysh self-energies 
arising from the electron-phonon interaction, to ﬁrst order. Oneﬁnds the general relation21 (see Appendix A)iˆA−ˆ R )}. (12)
ˆ K + [1 − 2n(ε,p; t)](I0 (ε,p) = − Tr{4
ˆ j (j = K, R, and A) are 2 × 2 matrices
The self-energies (in pseudospin space), and are functions of both energy ε andmomentum p; they include the effects of disorder averaging.
As we are considering a quasiequilibrium situation where bothphonons and electrons are individually in thermal equilibrium,n(ε,p; t) in Eq. (12) can be replaced by a Fermi distributionfunction at temperature Te , and the self-energies can becalculated assuming phonons are in thermal equilibrium attemperature Tph .
In the regime of interest (kF l  1), the dominant selfenergy diagram describing the leading-order electron-phonon
contribution to the kinetic equation (in the presence of disorder) is shown in Fig. 1. In this diagram, the wavy line representsa phonon propagator β (ω,q; η), where β = R,A,K denotes
phonon
(a)

, p
(b)

ω, q

e

e

 + ω, p + q
, p

, p
=

R/A (ω,q) =2ωq,2
ω − ωq2 ± iδ

, p
, p

+

 + ω, p + q  + ω, p + q  + ω, p + q
FIG. 1. (Color online) (a) Electron self-energy diagram. The
wavy line represents the phonon propagator, the solid line represents
the electron propagator, and the square block represents the dressed
electron phonon vertex by diffusion. (b) Vertex correction of theelectron-phonon interaction by impurity scattering. The dashed line
represents an impurity average. The square block represents dressedvertex and the dot represents bare vertex.

(13)

with ωq = sq, while the Keldysh propagator isK (ω,q) = [1 + 2N (ω,Tph )][R (ω,q) − A (ω,q)],
2 ε + 2τiA(ε,p) = −Im 2
ε + 2τi − vF2 p2eretarded, advanced, and Keldysh propagators. The retarded
and advanced phonon propagators appearing here are(14)

where N (ω,Tph ) is the Bose-Einstein distribution evaluated at
T = Tph .
The solid line in Fig. 1 represents an impurity-averaged
electronic Green function Ĝβ (ε,p); note that these are 2 ×
2 matrices in pseudospin space. The retarded and advanced
components are given by
ε ± 2τi + vF σ · p
ĜR/A (ε,p) = ,2
ε ± 2τi − vF2 p2(15)

while the Keldysh electron Green function isĜK (ε,q) = [1 − 2n(ε,Te )][ĜR (ε,q) − ĜA (ε,q)],(16)

where n(ε,T ) is now the Fermi-Dirac distribution function
evaluated at T = Te and chemical potential EF , where EF is
the Fermi energy (measured from the Dirac point).
Note that the electron Green function in Eq. (15) has anextremely simple form: It is just a free propagator with thesubstitution ε → ε ± i/2τ , corresponding to disorder-induced
broadening of energy levels. This broadening represents theﬁrst mechanism by which the electron-phonon heat ﬂux will
be modiﬁed due to disorder; this broadening generally causes a
suppression of the heat ﬂux. The second key effect of disorderis via the vertex correction of the electron-phonon vertex
appearing in the self-energy in Fig. 1. Each vertex describes theemission or absorption of a phonon; heuristically, the disorderinduced vertex correction corresponds to the modiﬁcation of
the amplitude of such a process due to the diffusive motionof electrons. The full details of the renormalization of thee-phonon vertices in a conventional (non-Dirac) metal bydisorder using the Keldysh formalism are presented in Ref. 14.
The Keldysh structure and renormalization of the e-phonon
vertices in graphene can be treated in a similar fashion. The
only key difference comes from the 2 × 2 matrix structureassociated with pseudospin; as we will see, this leads tointeresting new physical consequences.
The electron-phonon vertices M̂ (the η superscript is
dropped here) in the Keldysh technique are represented byγ
the form M̂μ,ν , where the upper index is for phonons andthe lower for electrons. The indices γ ,μ,ν each have two
possible values, cl and q, due to the two-component structureof the electron and phonon ﬁelds in the Keldysh formalism.32
The vertices with different indices are renormalized differently
by disorder in the Keldysh technique as shown in Ref. 14. Here,
we only present the simplest case, the renormalization of theqvertex M̂cl,cl . We stress here that our full calculation considers
all relevant Keldysh-space vertices and their disorder-induced
renormalization. Upon summation of all self-energy diagramsin the Keldysh formalism, one ﬁnds that the form of the renorq
malized vertex M̂cl,cl appears directly in the ﬁnal expression
for the collision integral, Eq. (21).13,14 The renormalization ofqvertex M̂cl,cl by disorder is depicted in Fig. 1(b). We focus on

125443-4

ELECTRON-PHONON MEDIATED HEAT FLOW IN . . .

PHYSICAL REVIEW B 86, 125443 (2012)

this speciﬁc vertex in the remainder of this section and drop
the upper and lower indices from now on.
In general, we may describe the renormalization of anelectron-phonon interaction vertex byM̂diff (q,ω) = Ď(q,ω) ◦ M̂0 (q),M̂0 = m0 · (1̂,σ̂x ,σ̂y ,σ̂z ),(18a)M̂diff = mdiff · (1̂,σ̂x ,σ̂y ,σ̂z ).

(18b)

Equation (17) takes the form
mdiff (q,ω) = D(q,ω) · m0 (q),(19)

where D is a 4 × 4 matrix. Using this representation, thesummation of ladder diagrams depicted in Fig. 1 results inthe form
D(q,ω) = [1 − P(q,ω)]−1 ,1d 2k[ĜR (k + q,ω)]ki
σ̂ijα[P(q,ω)]αβ =2π ντ(2π )2β−iω + Dq 2⎜
1 −1
0
D (q,ω) = ⎜
⎝
0
τ
0

0

2
−
iω+Dq
2 τ
0
0

1 1

1 
Tr
2 η=l,t



 η
dq dω
†
R(ε,ω)M̂0 (q) R (−ω; η)(2π )3η
× ĜR (ε + ω,p + q)M̂diff (−ω,−q) + H.c. ,(21)
where R(ε,ω) is the expected combination of Bose-Einstein
and Fermi-Dirac functions appropriate for phonon emission
and absorption processes:
R(ε,ω) = N (ω,Tph )n(ε,Te )[1 − n(ε + ω,Te )]− [1 + N (ω,Tph )][1 − n(ε,Te )]n(ε + ω,Te )= [n(ε,Te ) − n(ε + ω,Te )][N (ω,Tph ) − N (ω,Te )].
(22)
The ﬁrst term of R(ε,ω) describes absorption of a phonon ωfrom energy state ε to ε + ω and the second term describes
emission of a phonon ω from energy state ε + ω to state ε.
As expected, R(ε,ω) vanishes if Te = Tph . We note that apartfrom the matrix structure of the electron-phonon vertices andelectron propagators, the expression for the collision integral
has the same form as that found for conventional diffusive
metals.13,14,16 Nonetheless, we will see that the added matrix
structure (which encodes the chiral nature of the graphene
electronic excitations) gives rise to qualitatively new effects.
C. Diffusion propagator and renormalizationof the e-phonon vertex(20)Here, the indices α,β run from 0 to 3, and repeated indices
are to be summed over; we also use σ̂ 0 to denote the 2 × 2
unit matrix and σ̂ i (i = 1,2,3) are the Pauli matrices. P
describes a single “rung” in a standard diffusion ladder. Its
4 × 4 matrix structure is now directly related to the factthat each propagator carries an initial and ﬁnal charge or
pseudospin index. A similar structure is encountered when
considering diffusive dynamics in a system with strong spinorbit coupling or in 2D helical metals, as studied by Burkov
et al.34,35
We will discuss the properties and physics encoded in
the matrix diffusion propagator D in more detail in the next
section. For now, we only show how it enters in the ﬁnal
expression for the collision integral (and hence the heat ﬂux).
One ﬁnds that due to the causality structure of Keldysh
Green functions, the two vertices in Fig. 1(a) cannot bothbe simultaneously dressed by impurity scattering. Summingup all the self-energy diagrams in the Keldysh formalism, one⎛

I (ε,p) = −(17)

where M̂0 (M̂diff ) is the bare (renormalized) vertex, andĎ(q,ω) is a linear operator acting in the space of 2 × 2
matrices. It represents the diffusion propagator for electrons
in graphene, with the nontrivial matrix structure reﬂecting thefact that charge and pseudospin diffusion are linked together.
It is convenient to write this expression using a basis ofPauli matrices. Deﬁning the four-component vectors m,mdiff
via× [ĜA (k,0)]j l σ̂lk .

ﬁnally obtains

1. Diffusion propagator

It follows from Eqs. (11) and (21) that a key part of thedisorder-induced modiﬁcation of the electron-phonon heat ﬂux
is due to the modiﬁcation of the effective electron-phonon
interaction vertex. This modiﬁcation is in turn directly related
to the chiral diffusive dynamics of electrons in graphene, as
described by Eq. (17). In this section, we discuss the diffusion
propagator in more detail, as well as the forms of the dressed
electron-phonon vertices. Note that we restrict our discussion
here (as we do throughout the paper) on electrons in the K+
valley. While the sign of the chirality will be different for holes,
or for the K− valley, this sign has no impact on the quantity ofinterest, the e-phonon heat ﬂux in the presence of disorder.
We focus in this section on the most interesting diffusive
regime, where ql  1, ωτ  1 (i.e., we are interested inlength scales longer than l and time scales longer than τ ). In
this limit we can work to lowest nonvanishing order in ql andωτ . The inverse diffusion propagator D−1 (q,ω) = 1 − P(q,ω)
simpliﬁes to

0
0


1 1
2
−
iω+Dq
2 τ
0

⎞
⎛
0
0
⎜2ivF qx10⎟⎟+ ⎜0 ⎠ 4 ⎝2ivF qy
10τ

2ivF qy2ivF qx 2Dq
D qx2 − qy2
 x qy 
2Dqx qyD qy2 − qx2
0
0

⎞
0
0⎟
⎟,
0⎠
0
(23)

in the diffusive limit, where D = vF l/2 is the usual diffusion constant in two dimensions.
125443-5

WEI CHEN AND AASHISH A. CLERKPHYSICAL REVIEW B 86, 125443 (2012)or the diffusion equations of Cooperons in graphene.36 The
interpretation here is similar to Ref. 35: The coupling between
charge and pseudospin dynamics in the diffusive limit is a
result of the effective helicity of the electronic eigenstates.
By helicity, we mean that in the ε > 0 eigenstate of Eq. (3),pseudospin will be aligned with momentum. Thus, a positive
gradient in, say, Sx in the x direction implies a corresponding
positive gradient in the density of electronic x momentum. This
will then naturally cause the charge density N to decrease intime: This is the third term in Eq. (24a). Alternatively, writing
Eq. (24a) in the form of a continuity equation,

To gain intuition, it is useful to follow Ref. 34 and consider
the real-space representation of the matrix diffusion propagator, which describes the coarse-grained evolution of charge
and pseudospin densities [N (r,t) and Sj (r,t), respectively,
j = x,y,z]. The ﬁrst term in Eq. (23) would simply lead to
uncoupled equations for each of these quantities: N would
be described by a standard diffusion equation, while Sxand Sy would be described by diffusion equations with anadditional decay term (rate 1/τ ), corresponding to the fact thatpseudospin is not a conserved quantity. Sz has no dynamics inthe limit we consider, as it precesses with frequency EF andaverages away on the time scale 1/τ .
The second term in Eq. (23) complicates the above picture,
as it now links the dynamics of charge and pseudospin
densities. We thus obtain a set of coupled diffusion equations,describing the dynamics of these quantities (note that we havetaken into account the fact that the Pauli matrices are twice thepseudospin matrices):


∂Sy
∂Sx
∂N
2
= D∇ N − vF
+
,
(24a)
∂t
∂x
∂y

∂N
= −∇ · J,∂t
one sees that the charge current density J has the form
J = −D ∇N + vF (Sx x̂ + Sy ŷ).

These equations are analogous (but not identical) to the
diffusive dynamics for charge and spin in a 2D helical metal351(−iω+2Dq 2 )τ
⎜
iqx l
⎜
(iω−2Dq 2 )τD(q,ω) = ⎜
⎜
iqy l
⎝
(iω−2Dq 2 )τiqx l(iω−2Dq 2 )τiωτ +(qx2 −qy2 )l 2
3
+ 2(iω−2Dq2 )τ
2
qx qy l 2(iω−2Dq 2 )τiqy l(iω−2Dq 2 )τqx qy l 2(iω−2Dq 2 )τiωτ −(qx2 −qy2 )l 2
3
+ 2(iω−2Dq
2 )τ
2

0

0

0

Note that the effective diffusion constant (i.e., the coefﬁcient of
q 2 in the diffusion poles appearing above) is twice the value ofthe standardly deﬁned D appearing in Eqs. (24): Deff = 2D =vF l. This effective doubling of the diffusion constant is a directconsequence of the chiral nature of electrons in graphene, andis consistent with the results of previous transport studies.20,37
Also note that as expected from our discussion following
Eq. (23), only the charge-charge and charge-spin componentsof D [i.e., the top row and column] diverge in the limit of smallω and q. The lack of any corresponding large enhancementof the spin components of D is directly tied to the fact thatpseudospin is not a conserved quantity.
2. Disorder vertex correction of deformation-potential
e-phonon vertex

Having discussed the basic form of the matrix diffusion
propagator, we now turn to the renormalization of the e-phonon(26)

The ﬁrst term is the usual diffusive current, while the second
term corresponds to a “drift” current driven by the pseudospin
density following from the Hamiltonian in Eq. (3) (i.e., thecurrent operator is the pseudospin operator).
Turning to the dynamics of pseudospin densities, Eqs. (24b)and (24c) again reﬂect the fact that pseudospin is not conserved, and effectively decays on a time scale τ due to elastic
impurity scattering. In addition, we see that the diffusion ofthese densities is anisotropic: This is also a simple consequenceof helicity, as a net pseudospin density in a speciﬁc direction
also implies a net momentum density in this direction which
reinforces the diffusion in this direction.
Finally, inverting Eq. (23) [using as always the Pauli matrix
representation deﬁned in Eq. (19)], one ﬁnds the diffusion
propagator in Eq. (20) in the diffusive limit to be∂Sx
3D ∂ 2
D ∂2
∂2
Sx
=
Sy −
S
+
S
+
D
x
x
2
2
∂t
2 ∂x
2 ∂y
∂x ∂y
τ
vF ∂N
−
,
(24b)
2 ∂x
Sy
∂Sy
3D ∂ 2
∂2
D ∂2
S
+
S
+
D
=
Sx −
x
x
∂t
2 ∂x 22 ∂y 2
∂x ∂y
τ
vF ∂N
.
(24c)
−
2 ∂y

⎛

(25)

0

⎞

⎟
0⎟
⎟.
⎟
0⎠
1

(27)

interaction vertex due to diffusion, as given by Eqs. (19)and (27). Consider ﬁrst the DP contribution to the vertex. The
bare DP vertex is just a diagonal matrix in the sublattice basis[cf. Eqs. (6)] proportional to the coupling constant g1 . The
impurity-dressed vertex in the diffusive limit takes the form
iqξq g1
M̂DP,diff (q,ω) =τ

1(−iω+2Dq 2 )iqleiφq− (−iω+2Dq
2)

−iφq

iqle
− (−iω+2Dq
2)
1(−iω+2Dq 2 )


.
(28)

The diagonal parts of the vertex acquire a diffusion pole
protected by charge conservation, analogous to the case of
a normal metal. We will be interested in Eq. (28) with ω,qcorresponding to a thermal phonon, h̄ω = h̄sq kB T . Ass  vF in graphene, we thus have that over a wide range

125443-6

ELECTRON-PHONON MEDIATED HEAT FLOW IN . . .

PHYSICAL REVIEW B 86, 125443 (2012)of temperatures
Tdis < T < Tdis ,
vF

(29)

one has ω  Dq , and thus the diagonal parts of the vertex in
Eq. (28) will be enhanced by a factor ∼1/(q 2 l 2 ) compared tothe clean case. This corresponds to the well-known diffusive
enhancement of the electron-phonon interaction: The diffusive
motion effectively enhances the interaction time between an
electron and a long-wavelength phonon.
More surprisingly, Eq. (28) implies that the diffusive
renormalization of the deformation potential induces a vector
potential which is along the direction of wave vector q. This isa direct consequence of the chirality of the graphene electrons,
which links the dynamics of charge and pseudospin, and thus
allows a scalar potential to generate a vector potential (i.e., apotential which couples to pseudospin). In the next section, we
will show that this induced vector potential only gives a small
contribution to the heat ﬂux compared to the renormalized
deformation-potential coupling.

vertex which is a 2 × 2 matrix in pseudospin space. The barevertices areη

2

3. Disorder vertex correction of vector-potential
e-phonon vertex

The bare VP e-phonon vertex is given in Eqs. (6); it ispurely off diagonal in pseudospin space, and implies thata phonon of wave vector q = q q̂ = q(cos φq , sin φq ) andpolarization η generates an effective vector potential Al0 =g2 (sin 2φq , cos 2φq ) (η = l) or At0 = g2 (cos 2φq , sin 2φq )(η = t). It is useful to decompose these vectors into theirlongitudinal and transverse partsηηη

A0 = A,0 q̂ + A⊥,0 q̂⊥ ,(30)

where q̂⊥ is a unit vector perpendicular to q. One ﬁnds
Al,0 = g2 sin 3φq , Al⊥,0 = g2 cos 3φq ,(31a)At,0 = g2 cos φq , At⊥,0 = −g2 sin φq .

(31b)

Given the linearity of the vertex correction described by
Eq. (17), we can separately analyze how disorder changes
the interaction with the transverse and longitudinal phononinduced vector potentials. Each of these will yield an e-phononηM̂,0 = A,0 (cos φq σ̂x + sin φq σ̂y ),

(32a)

η
η
M̂⊥,0 = A⊥,0 (− sin φq σ̂x + cos φq σ̂y ).

(32b)

The corresponding renormalized vertices in the diffusive limit
are given by Eqs. (19) and (27), yieldingηM̂,diff =η2(iω − Dq 2 ) η
iql/τ
η
M̂ +
A σˆ0 , (33a)(iω − 2Dq 2 ) ,0 iω − 2Dq 2 ,0ηM̂⊥,diff = 2M̂⊥,0 .

(33b)

We see that the diffusive renormalization simply doubles thee-phonon vertex associated with the transverse part of thephonon-induced vector potential; there is no diffusion polehere, as there is no charge associated with a transverse vector
potential. In contrast, the longitudinal part acquires a diffusion
pole. The vector-potential part of the renormalized vertex is2)τ
simply the bare vertex multiplied by a factor of 2(iω−Dq.
(iω−2Dq 2 )τThis factor tends to 2 in the dc limit q → 0, ω → 0 and 1
in the regime we are most interested in in this work, i.e.,
ql > s/vF . The renormalization factor of 2 in the dc limit isconsistent with the renormalization of the current vertex ingraphene in the dc limit.21,37 Finally, we see that the vertex
associated with the longitudinal vector potential also acquiresa scalar potential [second term in Eq. (33a)]: This is an induced
deformation-potential coupling, again arising from the helicity
of electrons in graphene.
We stress that in contrast to the renormalized DP [cf.
Eq. (28)], Eqs. (33) explicitly show that there is no large
enhancement of the VP vertex in the ω = 0, q → 0 limit ofinterest. As discussed, the lack of a diffusive enhancement isa direct consequence of the nonconservation of pseudospin
and the consequent lack of a protected diffusion pole. The net
result is that the diffusive vertex correction discussed here does
not signiﬁcantly enhance the heat ﬂux associated with the VP
coupling at low temperatures.
Finally, for completeness, we give the full form of the
dressed VP e-phonon vertex in the diffusive limit. For theinteraction with LA phonons, combining Eqs. (30)–(33) yields

⎛

⎞
3
iql
iωτ/2  2iφqq 2 l 2 /2
−4iφq
−
sin
3φ
−i
+
−
i
e
e
q
2
2
2
2
2
2
2
−iωτ +q liωτ −q l−iωτ +q lη=l⎠.
M̂VP,diff (q,ω) = iqξql g2 ⎝ iωτ/2  −2iφqq 2 l 2 /2
iql
3
4iφq
i 2 + iωτ −q 2 l 2 e+ i −iωτ +q 2 l 2 e− −iωτ +q 2 l 2 sin 3φq

The full interaction vertex for TA phonons can be obtained ina similar fashion. One ﬁnds that vector potentials arising fromTA and LA phonons make identical contributions to the heat
ﬂux.
III. RESULTS AND DISCUSSIONS
A. Heat flux without screening(34)

we have all the necessary ingredients to evaluate Eq. (21) forthe electronic collision integral. From this, Eq. (11) directly
yields the desired electron-phonon heat ﬂux. As with thedisorder-free case, we again ﬁnd that the DP and VP couplings
contribute independently; we can thus meaningfully discuss
the ﬂux associated with each coupling. It is useful to express
each of these heat ﬂuxes in terms of an energy control function
Fα (T ) (α = DP,VP), deﬁned via

Having now determined both the renormalized electronphonon vertices [cf. Eqs. (28) and (34)] as well as thedisorder-averaged electronic Green functions [cf. Eq. (15)],125443-7Pα (Te ,Tph ) = νdε εI¯α,0 (ε) ≡ Fα (Te ) − Fα (Tph ). (35)

WEI CHEN AND AASHISH A. CLERKPHYSICAL REVIEW B 86, 125443 (2012)

Here, I¯α,0 (ε) is the momentum average of the collision integral
corresponding to the coupling mechanism α.
We discuss each mechanism in turn, focusing as always on
the regime vsF Tdis < T  TBG ; as discussed, the lower limit
here allows us to ignore the frequency dependence of therenormalized electron-phonon vertices. We also ﬁrst discuss
our results in the absence of any electronic screening.

G A W K m2

104

The energy control function determining the DP heat ﬂux
through the whole temperature regime vsF Tdis < T  TBG isobtained from Eqs. (C2) and (C3) to be
 

 ∞
1
s
ν
1
dq q 3 ql 
+ 22
FDP (T ) = 4g12vF 2πρM 0
1 + q 2l2 q l


1
1
−
N (ωq ,T ).
1− 
(36)
ql
1 + q 2l2
Here, ωq = sq, and we have included an overall factor of
4 reﬂecting the identical contribution from both valleys andboth physical spin projections. The integrand reﬂects thecontribution from phonons having momentum q to the heat
ﬂux. The ﬁrst term corresponds to the contribution stemming
only from the diagonal parts of the renormalized DP electronphonon vertex, whereas the second term corresponds to offdiagonal terms (i.e., the effective vector potential generated by
the chiral diffusion). The clean limit can easily be obtained bytaking l → ∞, yielding(37)

where we explicitly include factors of h̄ in the expression. This
result is the same as the heat ﬂux for deformation potential inclean graphene obtained in previous theoretical work.7,10,11
For ﬁnite l, both terms in Eq. (36) contribute. For
long-wavelength phonons satisyﬁng ql  1, the ﬁrst term
dominates the second (effective vector potential) term by alarge factor 1/(ql)2 . In contrast, for phonons having ql ∼ 1,
both terms make comparable contributions. In the diffusive
limit T < Tdis , the energy control function reduces to
FDP (T ) =EF2ζ (3) 2(kB T )3 .
g
π 2 1 h̄4 ρM s 2 vF3 l

1
0.01

10 6
0.01

T2

1

10

100

FIG. 2. (Color online) Thermal conductance per unit area G/A
associated with the deformation-potential coupling versus temperature T , including the effects of disorder, but without electronic
screening. We have taken a bare coupling constant g1 = 20 eV, carrier
density n = 1012 /cm2 , and mean free path l = 20 nm. The black solidline is the full result of our theory. The green-dashed line shows theasymptotic T 2 dependence in the low-temperature T  Tdis limit,
whereas the red-dotted line shows the asymptotic T 3 behavior in thehigh-temperature (clean) limit.

but rather a piezoelectric coupling. For these reasons, theseprevious works obtain different power laws for the e-phonon
heat ﬂux than those reported in this work, i.e., a powerlaw of T 5 in the clean limit and T 4 in the disordered limit
for deformation-potential coupling in the weak screening
case.15,16 The same power law of T 5 in the clean case andT 4 in the disordered case were also obtained for piezoelectric
coupling in the strong screening case.38–41
As noted in the Introduction, the expression in Eq. (38) isexactly twice of that of the e-phonon heat ﬂux associatedwith impurity-assisted “supercollisions,” as described bySong et al. in Ref. 1, for temperatures T  TBG . Unlike ourstudy, Ref. 1 considers e-phonon scattering dressed by onlya single impurity scattering event, and simply sums the effectof each impurity. In our approach, the q dependence of the
105

(38)

The power law becomes T 3 instead of T 4 ; this is the resultof the diffusive enhancement of the diagonal parts of therenormalized DP vertex, in complete analogy to what happens
in a conventional disordered metal.13–16,38–40 The corresponding thermal conductance G = (d/d)P (T + ,T ) associatedwith the DP coupling is shown in Figs. 2 and 3.
We note that the T 3 power law in Eq. (38), valid in
the disorder-dominated regime T < Tdis , is identical to what
would be obtained for conventional 2D electron gas coupling
via deformation potential to 2D phonons. Previous works
studying related physics15,16,38–41 consider a conventional 2D
electron gas embedded in a 3D structure, so that electronsinteract with 3D phonons. For the common case of a 2D
electron gas in a GaAs/Alx Ga1−x As heterostructure (asstudied in Refs. 38–41), the dominant e-phonon interaction
at low temperature is not a deformation-potential coupling,

0.1

TK

1000
G A W K m2

l→∞

EF
π2 2
(kB T )4 ,g
15 1 h̄5 ρM vF3 s 3

100

10 4

1. Deformation-potential heat flux

FDP,clean ≡ lim FDP (T ) =T3

l = 20 nm

10

l = 20 nm
l = 100 nm
l = 10 Μm

0.1
0.001
10 5
10 7
0.01

0.1

1

10

100

TK

FIG. 3. (Color online) Thermal conductance per unit area G/A
associated with the deformation-potential coupling versus temperature T , including the effects of disorder, showing the effects ofvarying the mean free path l as indicated; screening is neglected.
Remaining parameters are the same as Fig. 2. Both the enhancementof the low-temperature thermal conductance and shift of the crossover
temperature with increasing disorder are clearly evident.

125443-8

ELECTRON-PHONON MEDIATED HEAT FLOW IN . . .

PHYSICAL REVIEW B 86, 125443 (2012)

integrand in Eq. (36) in the dirty limit l → 0 can ultimately
be traced to the diffusive enhancement of the DP vertex [cf.
Eq. (28)]. In contrast, the corresponding q dependence in Ref. 1
can be traced to the energy of a virtual electronic state in asecond-order process involving both e-phonon and e-impurity
scattering events.

G A W K m2

10

2. Vector-potential coupling heat flux

The energy control function FVP (T ) for vector-potential
coupling in the temperature regime vsF Tdis < T  TBG isobtained from Eqs. (C4) and (C5) to besν
FVP (T ) = 4g22
dq q 3vF 2πρM 
1
1
× ql 
+2(1 + q 2 l 2 )1 + q 2l21
1
N (ωq ,T ).
−
1− 
(39)
2ql
1 + q 2l2
We have again included an overall factor of 4 in front to take
into account the valley and spin degeneracy. The ﬁrst term inthe square bracket is the contribution from the vector potential
after renormalization and the second term is the contribution
from the induced deformation-potential coupling from vector
potential. At ql  1, the second term gives a contribution
about 1/6 of the ﬁrst term, while at ql  1, the second term
is smaller by a factor of 1/2ql. In the clean limit T  Tdis ,the energy control function for vector potential has the samefunctional form as that for the DP coupling, i.e.,EF
π2 2(kB T )4 .
g
15 2 h̄5 ρM vF3 s 3(40)

In the opposite diffusive limit T  Tdis , the energy control
function for the VP coupling reduces to30ζ (5) 2 EF l(kB T )5 .
FVP (T ) =g2 6
π2h̄ ρM s 4 vF3

0.01
10 5
10 8

T4

10 11
0.01

0.1

B. Heat flux with electronic screening

We now consider how the above results are altered if one includes the screening of the e-phonon interaction. As discussed

1

10

100

FIG. 4. (Color online) Thermal conductance per unit area G/A
associated with the vector-potential coupling versus temperature T ,including the effects of disorder, but without electronic screening.
We have taken a bare coupling constant g2 = 1.5 eV, carrier density
n = 1012 /cm2 , and mean free path l = 20 nm. The black solid lineis the full result of our theory. The green-dashed line shows theasymptotic T 4 dependence in the low-temperature T  Tdis limit,
whereas the red-dotted line shows the asymptotic T 3 behavior in thehigh-temperature (clean) limit.

extensively by von Oppen et al.,24 the deformation-potential
coupling will be subject to screening at long wavelengths in
the usual manner, whereas the vector potential will not bescreened, as it does not induce any net charge density.
While the full dynamical screening of the deformation
potential can be important at low temperatures, it has beenshown that a static-screening approximation is sufﬁcient forphonon wave vectors q satisfying ql > s/2vF .15,42,43 Thus, in
our regime of interest 2vsF Tdis < T < Tdis , dynamic screening
will not be important, and a Thomas-Fermi approach is
sufﬁcient. Screening thus implies that the DP coupling constant
g1 in Eq. (6) now becomes q dependent:
g1,sc (q) = g1

(41)

Unlike the DP coupling, we see disorder increases the powerof temperature of the low-temperature heat ﬂux, indicating a
suppression of heat ﬂux from VP coupling. As discussed, thelack of a diffusive enhancement of the VP coupling is directlytied to the nonconservation of pseudospin; the main remaining
effect of disorder is a simple broadening of the electronic
eigenstates. The result is that for T  Tdis , disorder suppresses
VP-mediated heat transport. Note that as the bare VP coupling
constant g2 has been estimated to be more than an orderof magnitude smaller than the corresponding DP coupling
constant g1 ,18 it follows that in the absence of screening,
the heat ﬂux associated with the VP coupling is expected tobe negligible in comparison to that associated with the DP
coupling. The thermal conductance associated with the VP
coupling is shown in Figs. 4 and 5.

T3

TK

10
G A W K m2

FVP (T ) =l = 20 nm

0.01

q
,
q + qTF(42)

l = 20 nm
l = 100 nm
l = 10 Μm

10 5
10 8
10 11
0.01

0.1

1

10

100

TK

FIG. 5. (Color online) Thermal conductance per unit area G/A
associated with the vector-potential coupling versus temperature T ,including the effects of disorder, showing the effects of varying themean free path l as indicated; screening is neglected. Remaining
parameters are the same as Fig. 4. Both the suppression of the
low-temperature thermal conductance and shift of the crossover
temperature with increasing disorder are clearly evident.

125443-9

WEI CHEN AND AASHISH A. CLERKPHYSICAL REVIEW B 86, 125443 (2012)

where the Thomas-Fermi wave vector qTF is given by12,25
e
kF
κh̄vF

and κ is an effective dielectric constant. Using the value of κ
appropriate to graphene on a SiO2 substrate, one has qTF
3.2kF .12 Note that as we focus on the regime kF l  1, theeffects of screening will generally set in at a much higher
temperature Tsc = sqTF /kB than the temperature Tdis below
which disorder effects become important.
One can now easily include the effects of screening
into our theory by making the substitution g1 → g1,sc (q) in
Eq. (36) for the energy control function FDP (T ) determining
the deformation-potential mediated heat ﬂux. One ﬁnds
⎧ 8π 4 g2 E k6
6
1 F B
⎨
if Tdis  T  Tsc ,2 T ,63ρM h̄7 vF3 s 5 qTF
FDP,sc (T ) =⎩ 24 g12 EF ζ (5)T 5 , if T  T ,T .
dis scπ 2 ρ h̄6 v 3 s 4 q 2 l
M

F

TF

(44)
The suppression of the DP coupling by screening at low
temperatures T  Tdis implies that its associated heat ﬂux
can now become comparable or even smaller in magnitude
to that associated with the VP coupling [cf. Eq. (41)]. In the
low-temperature limit, both mechanisms yield energy control
functions F (T ) ∝ T 5 , with
FDP,sc (T )g2 1
∼ 12 2 2 .
FVP (T )g2 qTF l(45)

We see that the largeness of qTF l can compensate for the
relative smallness of g2 with respect to g1 , leading bothmechanisms to make comparable contributions. This behavior
is demonstrated in Fig. 6, where the thermal conductance

G A W K m2

104
10

DP unscreened
DP screened
VP

0.01
10 5
10 8
10 11
0.01

0.1

1

100

(43)

10

100

TK

FIG. 6. (Color online) Thermal conductance per unit area G/Aversus temperature T , showing the effects of electronic screening.
g1 = 20 eV, g2 = 1.5 eV, a carrier density n = 1012 /cm2 , and a mean
free path l = 20 nm. The red short-dashed line corresponds to an
unscreened deformation-potential coupling, while the blue dasheddotted line corresponds to a screened deformation-potential coupling,with a Thomas-Fermi wave vector qTF = 3.2kF as appropriate
for graphene on SiO2 (Ref. 12). The dashed purple curve is thecontribution from the vector-potential coupling. Despite its muchsmaller bare coupling constant, we see that at low temperatures,
both the deformation-potential and vector-potential couplings makealmost equal contributions when both screening and disorder effects
are included.

G A W K m2

qTF = 4

105

2

0.1

DP unscreened
DP screened
VP

10 4
10 7
10 10
10 13
0.01

0.1

1

10

100

TK

FIG. 7. (Color online) Thermal conductance per unit area G/Aversus temperature T , showing the effects of electronic screening.
Parameters are identical to Fig. 6, except we have now taken themean free path to be l = 10 μm, meaning that we are effectively
in the clean limit. Unlike the disordered case shown in Fig. 6, we
now see that at low temperatures, the contribution of the vectorpotential coupling dominates that from the deformation-potential
coupling.

versus temperature for both mechanisms is presented, forboth the strongly and weakly screened cases. The fact thatboth mechanisms are comparable is markedly different fromwhat happens in the screened, disorder-free case, which isrealized when Tdis  T  Tsc . In this case, the VP heat ﬂux
will dominate the DP heat ﬂux at low temperatures, as the VP
heat ﬂux scales like T 4 , while the DP heat ﬂux scales likeT 6 (see Fig. 7). This behavior in the clean limit is similar
to expectations for e-phonon contribution to the electricalresistivity, where it has also been argued that the VP coupling
can dominate at low temperatures.26
However, we note that a recent experiment measuring
the e-phonon contribution to the electrical resistivity of asuspended graphene ﬂake suggests that screening does notseem to be playing a role even when T < Tsc ,27 as theresults are compatible with the predictions for an unscreened
deformation-potential interaction (see also Ref. 25). The
same conclusions could be drawn from recent experiments
measuring the hot electron heat ﬂux in graphene.2–4
Finally, we note that with our theory, the only way to obtaina T 3 power law in the heat ﬂux at low temperatures [i.e.,δ = 3 in Eq. (1)] is via an unscreened deformation potential.
Thus, measurements of the low-temperature heat ﬂux could
also serve as a diagnostic tool for assessing the importance ofscreening the deformation potential.

IV. CONCLUSIONS

We have presented a comprehensive theory showing how
electronic disorder modiﬁes the electron-phonon interaction ingraphene (both the vector-potential and deformation-potential
couplings), and how this in turn has observable consequences
for the heat ﬂux between the electrons and lattice (acoustic)phonons. We focused on the relatively simple situation where
the graphene is doped away from the Dirac point, and where
the impurity potential can be considered smooth on atomic

125443-10

ELECTRON-PHONON MEDIATED HEAT FLOW IN . . .

PHYSICAL REVIEW B 86, 125443 (2012)

scales (implying that the disorder potential preserves thepseudospin and valley symmetries of the graphene Hamiltonian). We found that the unusual diffusion dynamics ofelectrons in graphene that results from their chirality also
has implications for how disorder modiﬁes electron-phonon
physics. We also found that this modiﬁcation is quite different
for the deformation-potential coupling versus the effective
vector-potential coupling. In the absence of screening, thecontribution to the heat ﬂux in Eq. (1) from both couplings
has temperature dependence of T 4 in the clean limit T  s/ l,consistent with previous work. In the disorder limit T < s/ l,
however, disorder affects the two types of couplings differently
and the power law of heat ﬂux for the two couplings becomesdifferent. We found that the total effect of disorder enhances
the heat ﬂux from DP coupling; however, it suppresses theheat ﬂux from VP coupling. The power law in Eq. (1) for
DP coupling becomes T 3 in the disorder limit while for
VP coupling, the power law becomes T 5 in the same limit
without screening. Thus, without screening, the DP couplingis expected to dominate heat transport in both the clean limit
and disorder limit, given that DP coupling is believed to bemuch stronger than vector-potential coupling.
We also considered the effects of screening on the above
physics, adopting the Thomas-Fermi approximation. Only
the DP is affected by screening; it suppresses it, and thus
at low temperatures, its contribution to the heat ﬂux scales
like T 6 in the clean limit and T 5 in the disordered limit.
Without disorder, one would thus expect the VP to dominate
at low temperatures due to the screening of DP, similar
to expectations for the electron-phonon contribution to theresistivity in clean graphene sheets.26 However, when disordereffects are included along with screening, we ﬁnd that boththe DP and VP coupling mechanisms can make comparable
contributions at low temperature.
In the future, it would be interesting (though challenging)to extend these results to situations of lower doping (where
the Fermi energy is closer to the Diract point), and to include
a richer structure of disorder that can break the symmetriesof the clean graphene Hamiltonian (as has been done in, e.g.,
studies of weak localization20 ).

ACKNOWLEDGMENTS

We acknowledge useful discussions with Anton V. Andreev,
K. C. Fong, Tami Pereg-Barnea, and K. C. Schwab. This
work was supported by NSERC and the Canadian Institutefor Advanced Research (CIFAR).

APPENDIX A: KELDYSH FORMALISM OF THE
KINETIC EQUATIONThe Green’s function in the Keldysh space is a matrix as
Ĝ =

 K
Ĝ

ĜR

ĜA

0


,

(A1)

where ĜR and ĜA are the retarded and advanced Green’s
functions, respectively. The Keldysh component ĜK in general

could be parametrized as ĜK = ĜR ◦ F − F ◦ ĜA , where F is
the Hermitian distribution function matrix. The Green function
obeys the following Dyson equation:0 R −1
ˆR
Ĝ0
−

 
 A −1ˆA−Ĝ0ĜK ĜR= 1, (A2)
◦
ˆK
ĜA 0
−

where the bare retarded and advanced Green functions in
R/Agraphene are (Ĝ0 )−1 (x,x  ) = δ(x − x  )(i∂t  + i σˆ · ∂r  ) in
ˆ K are the retarded, advanced,
ˆ R/A and 
the sublattice basis. 
and Keldysh components of the electron self-energy, respectively. The circle ◦ indicates integration over intermediate
coordinates.
The equation for the Keldysh component reads −1  R −1
ˆK +F ◦ˆA−
ˆ R ◦ F.
F ◦ ĜA
− Ĝ0
◦F =
0
(A3)
The distribution function matrix F (x,x  ) in graphene could
be decomposed to four components,
F (x,x  ) =3
fi (x,x  )σ̂ i ,

(A4)

i=0

where σ̂ 0 is the 2 × 2 unit matrix and σ̂ i , i = 1,2,3 are thePauli matrices.
The Wigner transformation of Eq. (A3) gives the kinetic
equation of the distribution function as−3

i σ̂ ∂τ fi (t,ρ; ε,p) − 2ivF
i

i=0

3


ij k σ̂ k kj fi (t,ρ; ε,p)i,j,k=1

ˆ K (ε,p) + F (t,ρ; ε,p)
ˆ A (ε,p) − ˆ R (ε,p)F (t,ρ; ε,p)=
= I (ε,p; t),(A5)where the function f0 (t,ρ; ε,p) = 1 − 2n(t,ρ; ε,p) andn(t,ρ; ε,p) is the charge density distribution function,fi (t,ρ; ε,p), i = 1,2,3 represent the pseudospin density distribution, and ij k is the three-dimensional antisymmetrictensor. The parameters t and ρ are center-of-mass time andcoordinates, respectively. Since the translation symmetry inspace is restored after averaging over impurities, we drop theρ dependence in the text.
The right-hand side of this equation is the collision
integral I (ε,p; t) (times i) in the presence of interactions.
ˆ R, 
ˆ A , and ˆ K are, respectively, the retarded, advanced,
and Keldysh components of the electron self-energy due toe-phonon and impurity scatterings. From the left-hand sideof the above equation, one ﬁnds that the pseudospin density
distribution functions fi=1,2,3 are small in a factor of 1/EF τ
compared to f0 . To leading order, one can replace F (t,ρ; ε,p)on the right-hand side by the scalar charge component
f0 (ε,p; t).
The collision integral is a 2 × 2 matrix in the sublatticebasis, which can be decomposed to components of σ̂ 0 andσ̂ i (i = 1,2,3). The three components of σ̂ i (i = 1,2,3) give

125443-11

WEI CHEN AND AASHISH A. CLERKPHYSICAL REVIEW B 86, 125443 (2012)a measure of the pseudospin density induced by interaction,
while the σ̂ 0 component determines the collision integral forthe charge density distribution function and is the one ofinterest in this paper. From Eq. (A5), one gets the collision
integral for charge distribution function
I0 (ε,p) =Eq. (19) in such basis becomes
⎛
⎜

√

⎞T

1+q 2 l 2
q 2 l2⎟⎜ − i cos φ ⎟q ⎟mDP,diff (ω → 0,q) = iqξql ⎜⎜ qli⎟⎝ − ql sin φq ⎠1 ∂f0 (ε,p; t)i∂n(ε,p; t)ˆ K (ε,p)=−= − Tr{∂t2∂t4ˆ R )(ε,p)}
ˆA−
(A6)
+ [1 − 2n(ε,p; t)](

g1 , (B2)0where the superscript T means transpose of the column vector
to row vector and the same for mVP,diff below. Written in thesublattice basis, the renormalized deformation-potential vertex
becomes
√ 22
⎛
⎞
1+q l
ie−iφq
1
+
−
2
2
q l
√ql 2 2 ⎠ g1 .
M̂DP,diff (ω → 0,q) = iqξql ⎝
1+q l
ieiφq
− ql
1 + q 2 l2

as shown in Eq. (12).
APPENDIX B: FULL FORM OF THE DIFFUSIONPROPAGATOR AND RENORMALIZED e-PHONON
VERTEX

The full form of the diffusion propagator D(ω,q) crossing
the whole temperature regime vsF Tdis < T  TBG is quite
complicated, yet in the temperature regime T > vsF Tdis , thefrequency dependence of the diffusion propagator can be
dropped and the typical thermal phonon wave vector q is muchsmaller than the Fermi wave vector. In this case, the diffusion
propagator in Eq. (20) is simpliﬁed toD(ω → 0,q)
√ 22
⎛
1+q l
1 + q 2 l2⎜⎜⎜ i
⎜ − ql cos φq=⎜⎜⎜ i⎜ − ql sin φq
⎝
0

1+

(B3)
In the diffusive limit ωτ  1, ql  1, it reduces to Eq. (28)(dropping the frequency dependence there), while in the clean
limit ql  1, it reduces to the bare vertex.
The four-vector representation of the bare vector potential
for the LA phonon in Eq. (6) ismVP,0 = iqξql (0, sin 2φq , cos 2φq ,0)g2− qli cos φq− qli sin φq

1+ √

− √

1−cos 2φq2

− √1+q 2 l 2

sin 2φq2

sin 2φq2

1+q 2 l 2

1+ √1+cos 2φq1+q 2 l 2

02

1+q 2 l 2

0⎞

in the Pauli matrix basis. The renormalized vertex according
to Eq. (19) then becomes

0

⎟
⎟
⎟
0⎟
⎟
⎟
⎟
0⎟
⎠
1
(B1)

mVP,diff (ω → 0,q)⎛−i sin 3φq
ql

2

1+q l

−

i sin 3φq
ql

2

1+q l

⎞T

⎜
⎟

sin 4φ⎜ sin 2φ 1 + √ 1− √ q2 2 ⎟q⎟
22
2 1+q l2 1+q l ⎟ g2
cos 4φq ⎟⎜ cos 2φ 1 + √ 1+ √ 22⎠q⎝
2 2⎜= iqξql ⎜in the Pauli matrix basis.
The four-vector representation of the bare deformationpotential vertex in the Pauli matrix basis is mDP,0 =iqξql (1,0,0,0)g1 . The renormalized vertex according to
⎛

(B4)

2

1+q l

0
(B5)
in the same basis. Written in the sublattice basis, it reads


M̂VP,diff (ω → 0,q) = iqξql ⎝ −2iφq 
4iφq
1 + √ 1 2 2 + i √e 2 2
ie
2

1+q l

⎞


−4iφq
−ie2iφq 1 + √ 1 2 2 − i √e 2 2
2 1+q l2 1+q l ⎠g2 .
i sin 3φ
− ql q(B6)

In the diffusive limit, it reduces to Eq. (34) (again dropping the frequency dependence there), while in the clean limit ql  1, it
reduces to the bare vertex.
APPENDIX C: DETAILS OF THE MATRIX COLLISION INTEGRAL
1. Deformation potential

The renormalized deformation potential in the whole regime of ql is presented in Appendix B. Plugging in the renormalized
deformation potential to Eq. (21) and integrating over the phonon frequency, one gets the collision integral for deformationpotential coupling as






dq 2
1
1 + q 2l2∂n(ε,p; t)
i
2
ε
+
ω
(qξ
=i
1
+
g
)
+
IDP,0 (ε,p) =
q 
q
2
∂t
(2π )2 12τ
q 2l2
ε + ωq + 2τi − vF2 |p + q|2

qx px + qy py
vF
)
−
(ω
→
−ω
,q
→
−q)
,
(C1)
+i
+
h.c.
R(ε,ω
+ ivF l
q
q
q
l
q 2l2
where ωq = sq.
125443-12

ELECTRON-PHONON MEDIATED HEAT FLOW IN . . .

PHYSICAL REVIEW B 86, 125443 (2012)

The collision integral for deformation potential after average over the electron momentum becomes








1
1
1 + q 2l2
1
dq 2 2 21
dp¯

IDP,0 (ε) =
q
1
+
1
−
A(ε,p)I
(ε,p)
=
τ
g
ξ
−
DP,0
q
πν
(2π )2(2π )2 1
q 2l2
q 2l2
1 + q 2l2
1 + q 2l2
× R(ε,ωq ) − (ωq → −ωq ,q → −q) .
(C2)
The heat ﬂux is thenPDP (Te ,Tph ) = νdε εI¯DP,0 (ε) = FDP (Tph ) − FDP (Te ),(C3)

where the energy control function FDP (T ) is presented in Sec. III.

2. Vector potential

Plugging in the renormalized vector potential in Appendix B to the collision integral Eq. (21) and separating the component forthe charge distribution function, one gets the kinetic equation of the charge distribution function due to vector-potential coupling
after averaging over the angle of electron momentum as






∂n(ε; t)1
1
1
dq 2 2 2
1



q
1
+
1
−
R(ε,ωq )gξ−= I¯VP,0 (ε) = τq∂t(2π )2 2
2q 2 l 2
2 1 + q 2l2
1 + q 2l2
1 + q 2l2
(C4)
− (ω → −ωq ,q → −q) .
The heat ﬂux due to vector-potential coupling is
PVP (Te ,Tph ) = νdε εI¯VP,0 (ε) = FVP (Te ) − FVP (Tph ),(C5)where the energy control function FVP (T ) is presented in Sec. III.

1

J. C. W. Song, M. Y. Reizer, and L. S. Levitov, Phys. Rev. Lett. 109,
106602 (2012).
2
A. M. R. Baker, J. A. Alexander-Webber, T. Altebaeumer, and
R. J. Nicholas, Phys. Rev. B 85, 115403 (2012).
3
Kin Chung Fong and K. C. Schwab, Phys. Rev. X 2, 031006(2012).
4
A. C. Betz, F. Vialla, D. Brunel, C. Voisin, M. Picher, A. Cavanna,
A. Madouri, G. Fève, J.-M. Berroir, B. Plaçais et al., Phys. Rev.
Lett. 109, 056805 (2012).
5
J. Ziman, Principles of the Theory of Solids (Cambridge University
Press, Cambridge, UK, 1972).
6
F. Giazotto, T. T. Heikkilä, A. Luukanen, A. M. Savin, and J. P.
Pekola, Rev. Mod. Phys. 78, 217 (2006).
7
R. Bistritzer and A. H. MacDonald, Phys. Rev. Lett. 102, 206410(2009).
8
W.-K. Tse and S. Das Sarma, Phys. Rev. B 79, 235406 (2009).
9
M. L. Roukes, K. Schwab, E. A. Henriksen, and J. M. Worlock,
Nature (London) 404, 974 (2000).
10
S. S. Kubakaddi, Phys. Rev. B 79, 075417 (2009).
11
J. K. Viljas and T. T. Heikkilä, Phys. Rev. B 81, 245404 (2010).
12
S. Das Sarma, S. Adam, E. Hwang, and E. Rossi, Rev. Mod. Phys.
83, 407 (2011).
13
B. L. Al’tshuler, Sov. Phys. JETP 48, 670 (1978).
14
M. Y. Reizer and A. V. Sergeyev, Sov. Phys. JETP 63, 616 (1986).
15
A. Sergeev and V. Mitin, Phys. Rev. B 61, 6041 (2000).

16

A. Sergeev, M. Y. Reizer, and V. Mitin, Phys. Rev. Lett. 94, 136602(2005).
17
J. González, F. Guinea, and M. A. H. Vozmediano, Phys. Rev. Lett.
69, 172 (1992).
18
H. Suzuura and T. Ando, Phys. Rev. B 65, 235412 (2002).
19
A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, andA. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
20
K. Kechedzhi, E. McCann, and V. I. Fal’ko, Euro. Phys. J. 148, 39(2007).
21
K. Kechedzhi, O. Kashuba, and V. I. Fal’ko, Phys. Rev. B 77,
193403 (2008).
22
B. L. Altshuler, D. Khmel’nitzkii, A. I. Larkin, and P. A. Lee, Phys.
Rev. B 22, 5142 (1980).
23
S. Hikami, A. I. Larkin, and N. Nagaosa, Prog. Theor. Phys. 63,
707 (1980).
24
F. von Oppen, F. Guinea, and E. Mariani, Phys. Rev. B 80, 075420(2009).
25
H. Min, E. H. Hwang, and S. Das Sarma, Phys. Rev. B 83, 161404(2011).
26
E. Mariani and F. von Oppen, Phys. Rev. B 82, 195403 (2010).
27
D. K. Efetov and P. Kim, Phys. Rev. Lett. 105, 256805 (2010).
28
P. R. Wallace, Phys. Rev. 71, 622 (1947).
29
J. C. Slonczewski and P. R. Weiss, Phys. Rev. 109, 272 (1958).
30
E. H. Hwang and S. Das Sarma, Phys. Rev. B 77, 115449(2008).

125443-13

WEI CHEN AND AASHISH A. CLERK
31

PHYSICAL REVIEW B 86, 125443 (2012)

K. Kaasbjerg, K. S. Thygesen, and K. W. Jacobsen, Phys. Rev. B
85, 165440 (2012).
32
A. Kamenev and A. Levchenko, Adv. Phys. 58, 197 (2009).
33
J. Rammer, Quantum Field Theory of Non-Equilibrium States(Cambridge University Press, Cambridge, 2007).
34
A. A. Burkov, A. S. Núñez, and A. H. MacDonald, Phys. Rev. B
70, 155308 (2004).
35
A. A. Burkov and D. G. Hawthorn, Phys. Rev. Lett. 105, 066802(2010).
36
E. McCann, K. Kechedzhi, V. I. Fal’ko, H. Suzuura, T. Ando, and
B. L. Altshuler, Phys. Rev. Lett. 97, 146805 (2006).

37

N. H. Shon and T. Ando, J. Phys. Soc. Jpn. 67, 2421 (1998).
E. Chow, H. P. Wei, S. M. Girvin, and M. Shayegan, Phys. Rev.
Lett. 77, 1143 (1996).
39
R. Fletcher, Y. Feng, C. T. Foxon, and J. J. Harris, Phys. Rev. B 61,
2028 (2000).
40
D. V. Khveshchenko and M. Reizer, Phys. Rev. B 56, 15822 (1997).
41
E. Chow, H. P. Wei, S. M. Girvin, W. Jan, and J. E. Cunningham,
Phys. Rev. B 56, R1676 (1997).
42
D. V. Khveshchenko and M. Y. Reizer, Phys. Rev. Lett. 78, 3531
(1997).
43
M. Y. Reizer, Phys. Rev. B 40, 5411 (1989).
38

125443-14