Supplemental Material for:
First-principles electron-phonon interactions and electronic
transport in large-angle twisted bilayer graphene
Shiyuan Gao,
1
,
2
Jin-Jian Zhou,
3
Yao Luo,
1
Marco Bernardi
1
,
∗
1
Department of Applied Physics and Materials Science, and Department of Physics,
California Institute of Technology, Pasadena, California 91125, USA
2
Department of Physics and Astronomy, Institute for Quantum Matter,
Johns Hopkins University, Baltimore, Maryland 21218, USA
3
School of Physics, Beijing Institute of Technology, Beijing 100081, China
∗
To whom correspondence should be addressed. E-mail: bmarco@caltech.edu.
1
I. ELECTRONIC BAND STRUCTURE OF MLG AND LARGE-ANGLE TBLG
Figure S1 shows the electronic band structure of MLG and large-angle tBLG with 21.8
◦
and 13.2
◦
twist angles obtained from Wannier interpolation. For the two twist angles con-
sidered here, the tBLG systems have doubly degenerate Dirac cones with nearly the same
Fermi velocity as MLG. Therefore, the density of states (DOS) in our large-angle tBLG
systems vanishes at charge neutrality.
FIG. S1. Electronic band structure of MLG, 21.8
◦
tBLG and 13.2
◦
tBLG, calculated from DFT
(cyan solid line) and Wannier interpolation (blue dashed line).
2
II. E-PH COUPLING AND CARRIER SCATTERING RATES IN MLG
For an easier comparison with previous results in the literature, we show the carrier
scattering rates and
e
-ph coupling matrix elements in MLG in Figure S2. Figure S2 (a-b)
show carrier scattering rates in linear and log scale respectively, and the energy range on
the x-axis is slightly extended compared to Figure 3 in the main text. Under linear scale, an
almost-no-scattering window similar to those in Ref. [1] appears in our results. Although
such a window may seem to suggest that the important phonons are the optical phonons,
it is important to note that the scattering channels most relevant to transport are those
within a few
k
B
T
of the Fermi energy (weighted by the transport distribution function, i.e.
the
∂f
n
k
(
T
)
∂ε
n
k
term in Eq. (3)), which fall almost completely within this almost-no-scattering
window. This means this window needs to be finely resolved for the calculation of transport
properties, and it turns out that the acoustic phonon scattering would dominate here and
become almost the singular contributor to phonon-limited transport at temperatures below
∼
250 K.
The scattering rates due to acoustic phonons vanishes linearly at the Dirac point, as
shown by Figure S2 (a) (note the acoustic phonon scattering rates are multiplied by 5). The
reason for this behavior is that the energy and momentum matching for the acoustic phonon
would be determined by the intersection between the carrier dispersion (Dirac cone) and the
phonon dispersion (a much flatter cone in the energy-momentum space with a slope equal
to the sound velocity and its tip sitting at some point on the Dirac cone representing the
energy and momentum of the carrier). Because the sound velocity is much smaller than the
Fermi velocity in graphene, the length of the intersection arc is almost identical to the DOS
at the carrier energy. The scattering rates due to optical phonon is instead determined by
the DOS
±
ℏ
ω
ν
q
away from the carrier energy. This causes the dip in the optical phonon
scattering rates around -0.2 eV.
Figure S2 (c) shows the phonon dispersion in MLG color-coded according to e-ph coupling
strength. Because of the lack of interlayer interaction and layer-breathing mode, the phonon
branches that correspond to atomic motions perpendicular to the layer, as shown by the two
blue-colored phonon branches in the plot, does not contribute to e-ph coupling in MLG.
Therefore, the intervalley phonon scattering in tBLG shown in Figure. 4(d)-(f) of the main
text is an emergent phenomenon in tBLG without a monolayer counterpart.
3
FIG. S2. (a-b) Electron-phonon scattering rate as a function of electron energy in MLG, shown
separately for the acoustic and optical modes, in (a) linear scale and (b) log scale. The temperature
is set at 300 K and the Fermi energy, indicated by the dotted black line, is set at 0.1 eV above
the Dirac point. In plot (a) the acoustic phonon scattering rates are multiplied by 5. (c) Phonon
dispersions in MLG from DFPT and color-coded according to the log of the
e
-ph coupling strength,
log
|
g
ν
(
q
)
|
.
III. PHONONS, E-PH COUPLING, AND RESISTIVITY IN TBLG USING PHONONS
COMPUTED AT DIFFERENT TEMPERATURES
In the main text, we use phonon eigenvalues and eigenvectors from the temperature-
dependent effective potential (TDEP) method to compute
e
-ph coupling matrix elements
and transport properties. The results presented in the main text use phonons obtained for
a fixed temperature of 30 K for 21.8
◦
tBLG, and 300 K for 13.2
◦
tBLG. For completeness,
in Fig. S3 below we compare the phonon dispersions,
e
-ph matrix elements, and resistivity
in 21.8
◦
tBLG obtained using phonons computed at two different temperatures (30 K and
300 K) in TDEP. For a 300 K phonon temperature, the acoustic phonons have a slightly
higher energy and weaker
e
-ph coupling, which results in a slightly lower resistivity compared
to the 30 K phonon results shown in the main text. Overall, the temperature chosen for the
phonon calculation does not change the results appreciably
−
for example, the resistivity
differs by about 5
−
30% in the temperature range studied here when using phonons computed
at 30 K versus 300 K. Therefore, the error resulting from using a fixed TDEP phonon
temperature for
e
-ph coupling and resistivity calculations over the entire temperature range
studied here does not impact our conclusions.
4
FIG. S3. Results for 21.8
◦
tBLG using two diffeent TDEP phonon temperatures. (a) Low-energy
phonon dispersions in 21.8
◦
tBLG, calculated with TDEP at 30 K (cyan solid line) and 300 K
(magenta dashed line). (b) The
e
-ph coupling matrix elements,
|
g
ν
(
k
,
q
)
|
versus
q
, for the six
low-energy modes in 21.8
◦
tBLG, computed with phonons from TDEP at 30 K (blue dots) and
300 K (magenta crosses). (c) Temperature dependent resistivity of 21.8
◦
tBLG, calculated with
phonons from TDEP at 30 K (cyan solid line) and 300 K (magenta dashed line).
IV. EFFECT OF THE RELAXATION TIME APPROXIMATION
ON THE RESISTIVITY OF MLG AND 21.8
◦
TBLG
In the main text, we use the relaxation time approximation (RTA) to compute the
phonon-limited resistivity. This can be improved by solving the full Boltzmann transport
equation using an iterative approach (ITA) [2]. In Fig. S4 below we compare the two ap-
proaches for both MLG and 21.8
◦
tBLG. In both systems, the difference between the RTA
and ITA resistivity is less than 3% below 300 K, and less than 5
−
10% at higher tempera-
tures. Due to the somewhat higher computational cost of ITA, especially for 13.2
◦
tBLG,
we choose to use the RTA throughout our work.
5
FIG. S4. Phonon-limited resistivity versus temperature, calculated using the RTA (solid lines) and
ITA (dashed lines), in MLG and 21.8
◦
tBLG. The Fermi level is 0.1 eV above the Dirac point.
V. SCALING OF E-PH MATRIX ELEMENTS WITH UNIT CELL SIZE
When comparing the
e
-ph matrix elements of MLG with tBLG at different twist angles,
it is important to know how the
e
-ph matrix elements
g
mnν
(
k
,
q
) depend on the size of the
simulation cell. In a supercell that is
N
times the size of the primitive unit cell, there will
be
N
times more electronic bands (indices (
m,n
)) and phonon branches (index
ν
) relative
to the unit cell calculation, and the Brillouin zone (BZ) will become 1
/N
times the unit-cell
BZ size (the BZ is spanned by the wave-vectors
k
,
q
). Therefore, there will be a total of
N
times more matrix elements
g
mnν
(
k
,
q
). However, we now show that only 1
/N
of them
are non-zero due to an additional phase condition that needs to be satisfied, and thus the
number of non-zero elements in
g
mnν
(
k
,
q
) is the same as for the primitive unit cell. In
addition, the non-zero elements in the supercell are 1
/
√
N
times the corresponding values
in the unit cell due to wave-function normalization. For this reason, in Fig. 2(c) of the main
text the
e
-ph matrix elements in 21.8
◦
tBLG are compared to the matrix elements in MLG
multiplied by 1
/
√
14, where the factor 1
/
√
7 comes from the wave-function normalization
in the supercell and an additional factor of 1
/
√
2 from the different normalization of the
phonon eigenvector as we go from a monolayer to a bilayer. These scaling relations ensure
that physical quantities such as the scattering rate remain unchanged (see below).
6
To give a practical example, we show how this scaling works for a simple 1D model of
e
-ph interactions. The scattering rate at
T
= 0, where only the phonon emission term is
present, can be calculated using Fermi’s golden rule:
Γ
nk
=
2
π
ℏ
L
BZ
X
mν
Z
BZ
dq
|
g
mnν
(
k,q
)
|
2
δ
(
ε
nk
−
ℏ
ω
νq
−
ε
mk
+
q
)
,
(1)
where
L
BZ
= 2
π/a
is the size of the BZ, and
a
is the length of the primitive unit cell. Consider
a toy model with band dispersion
ε
(
k
) =
−
t/
2
×
cos(
ka
), where
t >
0 is a hopping parameter,
coupled to a flat optical phonon branch with energy
ℏ
ω
O
through an
e
-ph coupling constant
g
. For the electronic state at
k
=
π/a
, the scattering rate in this model is
Γ
k
=
π/a
=
a
ℏ
Z
2
π/a
0
dq
|
g
|
2
δ
(
t/
2
−
(
t/
2) cos(
qa
)
−
ℏ
ω
O
)
.
(2)
If we further assume
ℏ
ω
O
=
t/
2, then the argument inside the delta function has two zeros
at
qa
=
π
±
π/
2. Therefore we have
Γ
k
=
π/a
=
2
a
ℏ
g
2
1
|
t
2
a
sin(
qa
)
|
qa
=
π
±
π/
2
=
4
g
2
ℏ
t
.
(3)
Let us consider how this picture changes if we describe the same system with a supercell
(SC) of size 2
a
. We first examine how
g
changes due to the folding of the BZ and the change
in the wave functions in the SC. Both the electronic bands and phonon dispersions are folded
onto the smaller BZ of size 2
π/
(2
a
) =
π/a
. Therefore, we introduce band indices (
n,m
)
spanning the two bands, and phonon mode indices
ν
spanning the two phonon branches,
and thus in the SC we write
g
as
g
mnν
. In the primitive unit cell, the Bloch wave function
is
φ
k
(
r
) =
e
ikr
u
k
, where
u
k
is the lattice-periodic part satisfying
u
k
(
r
) =
u
k
(
r
+
a
). In
the SC, using a reduced BZ scheme with
|
k
|
< π/
2
a
, the
|
k
|
> π/
2
a
region is folded
into the SC BZ, and the wave functions of the folded electronic states can be written as
φ
−
k
±
π/a
(
r
) =
e
i
(
k
±
π/a
)
r
(
e
∓
iπr/a
u
k
)
/
√
2, while the states with
|
k
|
< π/
2
a
remain unchanged
and can be written as
φ
+
k
(
r
) =
e
ikr
u
k
/
√
2. The 1
/
√
2 factor comes from the wave function
normalization in the SC, and we use the ’
−
’ and ’+’ symbols to label the two bands. The
SC-lattice periodic part of the wave function changes sign when we move to an adjacent
unit cell for the ’
−
’ band, and does not change sign for the ’+’ band. (Therefore, the SC-
lattice periodic parts for the ’+’ and ’
−
’ bands are orthogonal to each other for the same
k
).
Similarly, in the SC we have two phonon branches in which the atoms located one unit cell
apart either move in phase (+) or with opposite phase (
−
) relative to each other. Therefore,
7
in the SC we have a total of 8 different
e
-ph matrix elements:
g
+++
,
g
++
−
,
g
+
−
+
, etc.,
describing coupling between different electronic bands and phonon branches.
The
e
-ph matrix elements are calculated as
g
mnν
(
k,q
)
∝ ⟨
ψ
mk
+
q
|
∂
qν
V
|
ψ
nk
⟩
, which con-
tains an integral, over the SC of size 2
a
, of a function with periodicity of
a
multiplied by a
phase factor of either 1 or
e
±
iπr/a
, depending on whether there is an even or odd number of
’
−
’ indices in
g
. Therefore, the odd matrix elements
g
++
−
,
g
+
−
+
, etc. vanish, while the even
ones
g
+++
,
g
+
−−
, etc. are equal to
g/
√
2 due to the normalization factor, where
g
is the
e
-ph matrix element in the original primitive unit cell. Another quantity of interest is the
band-averaged
e
-ph matrix element
|
g
ν
|
=
p
P
mn
|
g
mnν
|
2
/N
b
, which is also equal to
g/
√
2.
We are now ready to check whether the SC description of the system gives the same
scattering rate, thus keeping the physical picture consistent. For simplicity, we choose a SC
BZ from
k
=
−
π/
2
a
+
ε
to
k
=
π/
2
a
+
ε
(where
ε
is a very small wave-vector) so that our
initial and final states in the example avoid the zone boundary. The initial state
k
=
π/a
is folded onto the ’
−
’ branch at
k
= 0, and in the SC BZ there are two final states (same
as in the primitive unit cell case) matching the energy-conservation condition given by the
delta function for
ℏ
ω
O
=
t/
2: these two final states are
k
=
π/
2
a
in the ’+’ branch, and
k
=
π/
2
a
in the ’
−
’ branch. The matrix elements connecting the initial and final states are
thus
g
+
−
ν
(
k
= 0
,q
=
π/
2
a
) and
g
−−
ν
(
k
= 0
,q
=
π/
2
a
), respectively. Plugging into Eq. (1),
we obtain
Γ
−
,k
=0
=
2
π
ℏ
L
SC
−
BZ
X
ν
Z
SC
−
BZ
dq
[
|
g
+
−
ν
|
2
δ
(
−
t/
2 cos
qa
) +
|
g
−−
ν
|
2
δ
(
t/
2 cos
qa
)]
.
(4)
Note that the size of the SC BZ
L
SC
−
BZ
is
π/a
, and thus we have
Γ
−
,k
=0
=
4
ℏ
t
X
ν
(
|
g
+
−
ν
|
2
+
|
g
−−
ν
|
2
)
(5)
Based on the analysis above,
g
+
−−
and
g
−−
+
are
g/
√
2, and
g
+
−
+
and
g
−−−
are 0, and thus
we arrive at the same scattering rate as for the primitive unit cell (see Eq. (3)), Γ
−
,k
=0
=
4
g
2
ℏ
t
.
Taken together, this derivation shows the correct way to transform the first-principles
e
-
ph matrix elements of a primitive unit cell calculation to compare with results from a larger
supercell calculation. To summarize, a calculation on a SC that is
N
times the size of the
primitive unit cell will give
e
-ph matrix elements smaller, by a factor 1
/
√
N
, than those of
the primitive cell. This scaling factor, derived here for the first time to our knowledge, is
essential to compare
e
-ph interactions in tBLG and MLG.
8
VI. INTERLAYER VERSUS INTRALAYER E-PH COUPLING IN 21.8
◦
TBLG
To separate the interlayer and intralayer contributions to
e
-ph interactions in tBLG,
we sum the
e
-ph matrix elements in the Wannier basis, ̃
g
ij
(
R
e
,R
p
), over pairs of Wannier
centers (
R
e
,R
p
) on the same layer for intralayer
e
-ph interactions, and on opposite layers
for interlayer
e
-ph interactions. The resulting
e
-ph matrix elements in momentum space are
given in Fig. S5(a) for
q
near the Γ point, which mediates intravalley scattering, and in
Fig. S5(b) for
q
near the K point, which mediates intervalley scattering. Figure S5(a) shows
that intralayer
e
-ph interactions give the primary contribution to intravalley scattering. Yet,
the interlayer contribution is non-negligible, and can be as large as 20% of the intralayer
contribution. Therefore, intravalley scattering in tBLG is greater than in MLG due to
the interlayer contribution, which comes primarily from the LB mode (see Fig. S5(a)).
This result still holds true at low temperature, where intervalley scattering is suppressed.
Conversely, Fig. S5(b) shows that intervalley scattering is dominated by
interlayer
e
-ph
interactions, and is thus significantly enhanced in tBLG compared to MLG.
FIG. S5. Electron-phonon coupling strength in 21.8
◦
tBLG for (a) intravalley scattering, with
q
near the Γ point, and (b) intervalley scattering, with
q
near the K point. The intralayer and
interlayer contributions, obtained by analyzing the
e
-ph matrix elements in the Wannier basis, are
plotted separately in both panels.
9
VII. PHONON-LIMITED RESISTIVITY OF 21.8
◦
TBLG FOR DIFFERENT
DOPING LEVELS
Figure S6 shows the phonon-limited resistivity of 21.8
◦
tBLG for Fermi energies of up
to 0.4 eV (which corresponds to a doping density of 3
.
57
×
10
13
cm
−
2
). Similar to the case
of acoustic phonon-limited resistivity of MLG [3], the main effect of changing the doping
level is a change of the Bloch-Gr ̈uneisen temperature
T
BG
, which is proportional to the
Fermi energy. In the power-law region between
T
BG
and 250 K, the resistivity curves for
different Fermi levels are very close to each other [see Fig. S6(a)], and by fitting the curve
for
E
f
=0.2 eV (see below) we find that the power law
ρ
∝
T
1
.
193
±
0
.
007
accurately describes
the resistivity vs. temperature trend. The curves being so close to each other also means
that the dependence of the resistivity on Fermi energy or doping concentration is very weak
at intermediate temperatures above
T
BG
, as confirmed by the results in Fig. S6(b).
FIG. S6. (a) Resistivity vs. temperature for four different Fermi levels above the Dirac point in
21.8
◦
tBLG. (b) The same resistivity as in (a) shown as a function of Fermi energy at three different
temperatures.
The exponents in the
ρ
versus
T
curves are obtained by fitting log(
ρ
) vs. log(
T
) with
linear regression, which also provides the standard error (computed under the assumption
of residual normality). In Table I below, we show the fitted exponents using data in three
different temperature ranges above
T
BG
. (The fitting results point to very small differences
10
in the exponents depending on which temperature region is used for fitting.) The results
show a clear deviation from a
T
-linear scaling in large-angle tBLG.
Temperature
range (K)
MLG
E
f
= 0
.
1 eV
21.8
◦
tBLG
E
f
= 0
.
1 eV
21.8
◦
tBLG
E
f
= 0
.
2 eV
13.2
◦
tBLG
E
f
= 0
.
1 eV
30
−
160
1.001
±
0.009
1.173
±
0.002
1.263
±
0.022
1.288
±
0.016
50
−
160
0.985
±
0.003
1.176
±
0.003
1.224
±
0.017
1.264
±
0.020
70
−
160
0.981
±
0.006
1.171
±
0.003
1.193
±
0.007
1.229
±
0.016
TABLE I. Slope in the log(
ρ
)-log(
T
) curve from linear regression in three different temperature
ranges. Results are given for MLG, 21.8
◦
tBLG at two different Fermi energies, and 13.2
◦
tBLG.
11
VIII. ELIASHBERG FUNCTION AND EFFECTIVE COUPLING CONSTANT
λ
FOR MLG, 21.8
◦
TBLG AND 13.2
◦
TBLG
From the first-principles
e
-ph coupling matrix elements, we also compute the Eliashberg
function using approaches similar to [4, 5]:
α
2
F
(
ω
) =
1
N
f
N
k
N
q
X
mnν
kq
|
g
mnν
(
k
,
q
)
|
2
δ
(
ε
n
k
−
E
f
)
δ
(
ε
m
k
+
q
−
E
f
)
δ
(
ℏ
ω
−
ℏ
ω
ν
q
)
,
(6)
where
N
f
is the DOS at the Fermi energy
E
f
, while
N
k
and
N
q
are the number of
k
- and
q
-points respectively. From the Eliashberg function, we calculate the effective coupling con-
stant as
λ
= 2
R
α
2
F
(
ω
)
/ω dω
. We compute
α
2
F
(
ω
) using
q
-point grids of 60
×
60 for MLG,
24
×
24 for 21.8 tBLG
◦
and 12
×
12 for 13.2
◦
tBLG, and
k
-point grids of 240
×
240 for MLG,
90
×
90 for 21.8 tBLG
◦
and 54
×
54 for 13.2
◦
tBLG. The
δ
-functions in Eq. 6 are replaced
by Gaussian functions with respective broadening values of 60 meV for the delta functions
involving electronic energies and 4 meV for the delta function involving phonon energies.
Figure S7 shows the Eliashberg function and effective coupling constant
λ
for MLG, 21.8
tBLG
◦
and 13.2
◦
tBLG. Panels (a)-(c) on the left show the Eliashberg function decomposed
into contributions from different phonon modes (for a Fermi energy of 0.2 eV) and panels
(d)-(f) on the right show the mode-decomposed and total
λ
as a function of Fermi energy.
These results show that the Eliashberg function is dominated by
e
-ph coupling with optical
phonons in MLG and large-angle tBLG. The contribution from acoustic low-energy inter-
valley phonon modes increases for decreasing twist angles. The effective coupling constant
λ
is nearly proportional to the Fermi energy, and the slope of
λ
vs. Fermi energy in 21.8
◦
tBLG is slightly higher than in MLG, and slightly lower than in 13.2
◦
tBLG.
12
FIG. S7. (a)-(c) Eliashberg function
α
2
F
(
ω
) for different families of phonon modes in (a) MLG,
(b) 21.8
◦
tBLG, and (c) 13.2
◦
tBLG. The Fermi energy is set to 0.2 eV for calculations in (a)-(c).
(d)-(f) Effective coupling constant
λ
vs. Fermi energy for the same systems as in (a)-(c) respectively.
For better visualization, the contribution from acoustic phonon modes is multiplied by 50 times in
panels (a)-(c) and by 3 times in (d)-(f).
[1] C.-H. Park, F. Giustino, M. L. Cohen, and S. G. Louie, Velocity renormalization and carrier
lifetime in graphene from the electron-phonon interaction, Phys. Rev. Lett.
99
, 086804 (2007).
[2] J.-J. Zhou, J. Park, I.-T. Lu, I. Maliyov, X. Tong, and M. Bernardi, Perturbo: A software
package for ab initio electron–phonon interactions, charge transport and ultrafast dynamics,
13
Comput. Phys. Commun.
264
, 107970 (2021).
[3] E. H. Hwang and S. Das Sarma, Acoustic phonon scattering limited carrier mobility in two-
dimensional extrinsic graphene, Phys. Rev. B
77
, 115449 (2008).
[4] E. R. Margine and F. Giustino, Two-gap superconductivity in heavily
n
-doped graphene: Ab
initio migdal-eliashberg theory, Phys. Rev. B
90
, 014518 (2014).
[5] Y. W. Choi and H. J. Choi, Dichotomy of electron-phonon coupling in graphene moir ́e flat
bands, Phys. Rev. Lett.
127
, 167001 (2021).
14