# Electron-$K$-Phonon Interaction In Twisted Bilayer Graphene

Chao-Xing Liu∗ Department of Physics, The Pennsylvania State University,
University Park, Pennsylvania 16802, USA  Department of Physics, Princeton
University, Princeton, NJ 08544  Yulin Chen  Department of Physics, University
of Oxford, Oxford, OX1 3PU, United Kingdom  Ali Yazdani  Department of
Physics, Princeton University, Princeton, NJ 08544  B. Andrei Bernevig∗
Department of Physics, Princeton University, Princeton, NJ 08544  Donostia
International Physics Center, P. Manuel de Lardizabal 4, 20018 Donostia-San
Sebastian, Spain  IKERBASQUE, Basque Foundation for Science, Bilbao, Spain

###### Abstract

We develop an analytic theory to describe the interaction between electrons
and $K$-phonons and study its influence on superconductivity in the _bare
bands_ of twisted bilayer graphene (TBG). We find that, due to symmetry and
the two-center approximation, only one optical $K$ phonon ($\sim 160meV$) of
graphene is responsible for inter-valley electron-phonon interaction. This
phonon has recently been found in angular-resolved photo-emission spectroscopy
to be responsible for replicas of the TBG flat bands. By projecting the
interaction to the TBG flat bands, we perform the full symmetry analysis of
phonon-mediated attractive interaction and pairing channels in the Chern
basis, and show that several channels are guaranteed to have gapless order
parameters. From the linearized gap equations, we find that the highest
$T_{c}$ pairing induced by this phonon is a singlet gapped s-wave inter-Chern-
band order parameter, followed closely by a gapless nematic d-wave intra-
Chern-band order parameter. We justify these results analytically, using the
topological heavy fermion mapping of TBG which has allowed us to obtain an
analytic form of phonon-mediated attractive interaction and to analytically
solve the linearized and $T=0$ gap equations. For the intra-Chern-band
channel, the nematic state with nodes is shown to be stabilized in the chiral
flat band limit. While the flat band Coulomb interaction can be screened
sufficiently enough - around Van-Hove singularities - to allow for electron-
phonon based superconductivity, it is unlikely that this effect can be
maintained in the lower density of states excitation bands around the
correlated insulator states.

###### pacs:

73.43.-f, 71.10.Fd, 03.65.Vf, 03.65.Ud

Introduction -  Superconductivity in twisted bilayer graphene (TBG) appears
within its phase diagram around the correlated insulator states [1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 4, 11, 12, 13, 14, 15, 16, 17, 18, 19]. Amongst the
mechanisms suggested for superconductivity are the phonons, spin fluctuations,
skyrmions, and others [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36, 29, 37, 38, 39, 40, 41, 42, 33, 43, 44, 23, 31]. Based on a recent
experiment that suggest a strong coupling between the graphene $K$-phonon and
the flat bands in TBG [45], we perform a comprehensive analysis of the
electron-$K$-phonon (e-K-ph) interaction and the resulting phonon-mediated
superconductivity on the bare flat bands of TBG. We develop an exhaustive
numerical, analytical, and symmetry based description of the e-K-ph
interaction in TBG and the symmetry classifications of the order parameter,
and find the competing singlet gapped inter-Chern-band channel and nematic
gapless intra-Chern-band channel. Armed with the heavy-fermion description of
TBG[46, 47, 48, 49, 50, 51, 52, 53, 54], the form factors of the $K$-phonon
induced attractive interaction can be analytically computed and matched well
to a full numerical calculations. An analysis of the Coulomb screening shows
that, due to the high density of states (DOS) of flat bands, the Coulomb
interaction might be strongly renormalized down near the Van Hove
singularities. However, it remains unclear if the Hartree-Fock bands of the
correlated insulator, with the lower DOS, can provide a similar result.

Model Hamiltonian for electron-phonon interaction in TBG - We consider the
deformation potential type of theory, described by a tight-binding (TB) model
for the electron Hamiltonian with the hopping parameters depending on the atom
positions $\tilde{\bf R}_{\alpha}^{l}={\bf R}^{l}+{\bf\tau}_{\alpha}^{l}+{\bf
u}^{l}({\bf R}^{l}_{\alpha})$ with a displacement field ${\bf u}^{l}({\bf
R}^{l}_{\alpha}={\bf R}^{l}+{\bf\tau}_{\alpha}^{l})$, where ${\bf R}^{l}$ and
${\bf\tau}_{\alpha}^{l}$ label the lattice vector and the sublattice atom
position ($\alpha=A,B$) at the layer $l$, respectively. By treating ${\bf u}$
as a perturbation, we expand the intra-layer Hamiltonian up to the linear
order in $\mathbf{u}$ (Supplementary Material (SM) Sec. II [55]). We only keep
$\mathbf{u}$-independent terms for the inter-layer Hamiltonian for TBG, thus
focusing on intra-layer electron-phonon (e-ph) interaction in this work. As
only the Dirac bands appear around the Fermi energy close to
$\pm\mathbf{K}_{D}=\pm\frac{4\pi}{3a_{0}}(1,0)$ in Brillouin zone (BZ) with
lattice constant $a_{0}$ in graphene, we also expand the Hamiltonian around
$\eta\mathbf{K}_{D}$ ($\eta=\pm$ labelling two valleys) and focus on Dirac
electrons around two valleys. Our full Hamiltonian consists of three parts

$\displaystyle H=H_{el}+H_{ph}+H_{eph}.$ (1)

Here $H_{el}$ describes the Dirac electrons located at valley $\eta=\pm$
momentum $\eta\mathbf{K}_{D}$ that are coupled through inter-layer tunnelings
and is given by the Bistritzer-MacDonald (BM) model [56] (SM Sec. V [55]),

${\hat{H}}_{\rm el}=\sum_{\eta s}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\alpha\alpha^{\prime}}\sum_{\mathbf{Q},\mathbf{Q}^{\prime}}h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(\eta)}(\mathbf{k})c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}c_{\mathbf{k},\mathbf{Q}^{\prime},\alpha^{\prime},\eta,s}$
(2)

where $c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}$ is the fermion annihilation
operator $\mathbf{k}$ is a momentum in the Moiré Brillouin zone (MBZ) (Fig.
1b), $\alpha$ is the sublattice index, and $s$ is spin. The vector
$\mathbf{Q}$ belongs to the lattice set
$\mathcal{Q}_{l\eta}=\\{l\eta\mathbf{q}_{2}+n_{1}\mathbf{b}_{M1}+n_{2}\mathbf{b}_{M2}\
|\ n_{1,2}\in\mathbb{Z}\\}$, where $l$ is the layer index,
$\mathbf{q}_{2}=k_{\theta}(\frac{\sqrt{3}}{2},\frac{1}{2})$,
$\mathbf{b}_{M1}=k_{\theta}(\frac{\sqrt{3}}{2},\frac{3}{2})$,
$\mathbf{b}_{M2}=k_{\theta}(-\frac{\sqrt{3}}{2},\frac{3}{2})$ and
$k_{\theta}=2|\mathbf{K}_{D}|\sin\frac{\theta}{2}$ with $\theta$ the twist
angle.
$h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(\eta)}(\mathbf{k})$
is given in SM Sec. V.A [55]. ${\hat{H}}_{\rm el}$ exhibits $C_{6v}$ and time-
reversal symmetries, generated by valley-switching $\pi/6$-rotation along z
axis ($\hat{C}_{6z}$), time reversal ($\hat{T}$), and $\pi$-rotation along y
axis ($C_{2y}$), and valley-preserving $\pi$-rotation along x axis
($\hat{C}_{2x}$), and the composite anti-unitary $C_{2z}T$. In addition,
${\hat{H}}_{\rm el}$ has a unitary particle-hole ($\hat{P}$) symmetry, as well
as a chiral symmetry $\hat{C}$ in the limit with vanishing $AA$ region hopping
($w_{0}=0$) [19]. A full discussion of symmetry of BM model [57, 58] is found
in SM Sec.V.B [55].

$H_{ph}$ describes the intra-layer in-plane phonon modes. Out-of-plane phonon
modes are decoupled from Dirac electrons for intra-layer e-ph interaction. The
dynamical matrix for a single layer graphene is derived in SM Sec.III [55]
based both on symmetry considerations and the microscopic model, up to the
next nearest neighbor interaction. The resulting in-plane phonon dispersion in
Fig. 1a reproduces that in literature [59, 60, 61, 62, 63] (SM Sec.III and
IV.B [55]). The phonon modes at $\Gamma$ and $\eta\mathbf{K}_{D}$ can induce
intra-valley and inter-valley e-ph interactions, respectively. In this paper
we focus on the $\eta\mathbf{K}_{D}$ phonons; the $\Gamma$ phonons will be
derived in [64]. At $\eta\mathbf{K}_{D}$, we have one $A_{1}$ ($\sim 160meV$),
one $A_{2}$ ($\sim 140meV$) and one 2D $E$ mode ($\sim 150meV$) of $C_{3v}$
group. Based on the deformation potential theory, we derive the e-ph
interaction $H_{eph}$ by expanding the TB Hamiltonian treating both the
momentum and phonon displacement field $\mathbf{u}$ as perturbations. For e-ph
interaction, we only keep the dominant zeroth-order in momentum for the
$\eta\mathbf{K}_{D}$ phonons. We find, due to both symmetry and the two-center
approximation (SM Sec.II.E [55]), that only the $A_{1}$ phonons at
$\mathbf{K}_{D}$ can scatter an electron from $\mathbf{K}_{D}$ to
$-\mathbf{K}_{D}$ [63]. The corresponding Hamiltonian reads

$\displaystyle H_{inter-
vall}^{op,A_{1}}\approx\frac{\gamma_{3}}{\sqrt{2N_{G}M\omega_{A1}}}\sum_{\tilde{\mathbf{k}},\tilde{\mathbf{k}}^{\prime},\eta,\alpha\beta}(b_{-\eta\mathbf{K}_{D}+\tilde{\mathbf{k}}-\tilde{\mathbf{k}}^{\prime},A_{1}}+b_{\eta\mathbf{K}_{D}-\tilde{\mathbf{k}}+\tilde{\mathbf{k}}^{\prime},A_{1}}^{\dagger})c^{\dagger}_{\tilde{\mathbf{k}}+\eta\mathbf{K}_{D},\alpha}(\sigma_{x})_{\alpha\beta}c_{\tilde{\mathbf{k}}^{\prime}-\eta\mathbf{K}_{D},\beta},$
(3)

where $\tilde{\mathbf{k}}$ is the electron momentum away from
$\eta\mathbf{K}_{D}$, $N_{G}$ is the number of atomic unit cells, $M$ is the
atomic mass, $\omega_{A1}$ is the $A_{1}$ phonon frequency, and $b$ and $c$
are phonon and electron annihilation operators. The material dependent
parameter $\gamma_{3}$ can be derived from the hopping potential as
$\gamma_{3}=2i\sum_{\mathbf{G}}e^{i(\tau_{A}-\tau_{B})\cdot\mathbf{G}}(\mathbf{G}+\mathbf{K}_{D})_{y}t(\mathbf{G}+\mathbf{K}_{D},0)\approx
17\text{eV/ \AA}$, where $\mathbf{G}$ is the reciprocal lattice vector and
$t(\mathbf{q})$ is the Fourier transform of the $\pi$-bond hopping function
between two carbon $p_{z}$ orbitals in graphene (Eq. 6 in SM Sec.I [55]). Our
next step is to re-write the electron momentum $\tilde{k}$ into the MBZ by
$\tilde{\mathbf{k}}=\mathbf{k}-\mathbf{Q}_{l\eta}$ with
$\mathbf{k}\in\text{MBZ}$, so that
$c_{\mathbf{k},\mathbf{Q}_{l\eta},\alpha,\eta,s}=c_{\eta\mathbf{K}_{D}^{l}+\tilde{\mathbf{k}},\alpha,l,s}$
and
$\sum_{\tilde{\mathbf{k}}}\rightarrow\sum_{\mathbf{k}\in\text{MBZ}}\sum_{\mathbf{Q}_{l\eta}}$,
where we have added the spin index $s$ and layer index $l$. Finally, we
project the e-ph interaction $H_{inter-vall}^{op,A_{1}}$ into the flat bands
of the BM Hamiltonian as

$\displaystyle H_{inter-vall}^{op,A_{1}}\approx\frac{1}{\sqrt{N_{G}}}\sum
G_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{Q}_{-l\eta}}^{\eta
nn^{\prime}l}\gamma_{\mathbf{k},n,\eta,s}^{\dagger}\gamma_{\mathbf{k}^{\prime},n^{\prime},-\eta,s}$
$\displaystyle(b_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime}-\mathbf{Q}_{-l\eta},l,A_{1}}+b_{\eta\mathbf{K}_{D}-\mathbf{k}+\mathbf{k}^{\prime}+\mathbf{Q}_{-l\eta},l,A_{1}}^{\dagger})$
(4)

where the summation includes
$\mathbf{k},\mathbf{k}^{\prime},n,n^{\prime},\eta,s,l,\mathbf{Q}_{-l\eta}$,
$\gamma_{\mathbf{k},n,\eta,s}^{\dagger}=\sum_{\mathbf{Q}\alpha}u_{\mathbf{Q}\alpha;n\eta}\left(\mathbf{k}\right)c_{\mathbf{k},\mathbf{Q},\eta,\alpha
s}^{\dagger}$ with
$u^{n}_{\mathbf{k},\mathbf{Q}^{\prime}_{l\eta},\alpha,\eta}$ the eigenstates
of
$h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(\eta)}(\mathbf{k})$.
The matrix element

$\displaystyle G_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{Q}_{-l\eta}}^{\eta
nn^{\prime}l}=\frac{\gamma_{3}}{\sqrt{2M\omega_{A_{1}}}}\sum_{\mathbf{Q}^{\prime}_{l\eta},\alpha\beta}$
$\displaystyle
u^{n\star}_{\mathbf{k},\mathbf{Q}^{\prime}_{l\eta},\alpha,\eta}\sigma^{x}_{\alpha\beta}u^{n^{\prime}}_{\mathbf{k}^{\prime},\mathbf{Q}^{\prime}_{l\eta}-\mathbf{Q}_{-l\eta},\beta,-\eta}$
(5)

characterizes the e-ph interaction strength for TBG and can be evaluated
numerically (and later analytically), as shown in SM Sec. VI.F [55]. We focus
on two flat bands (per valley per spin) of TBG, labelled by $n=\pm$. Instead
of the eigen-state basis, we work on the so-called ”Chern-band” basis, defined
by

$u^{e_{Y}}_{\mathbf{k},\mathbf{Q},\alpha,\eta}=\frac{1}{\sqrt{2}}(u^{n=+}_{\mathbf{k},\mathbf{Q},\alpha,\eta}+ie_{Y}u^{n=-}_{\mathbf{k},\mathbf{Q},\alpha,\eta})$
(6)

with $e_{Y}=\pm 1$. $u^{e_{Y}}_{\mathbf{k},\mathbf{Q},\alpha,\eta}$ carries
the Chern number $\pm 1$. On the Chern-band basis, the expressions for e-ph
interaction can be obtained by replacing the $n,n^{\prime}$ indices in Eqs.
(Electron-$K$-Phonon Interaction In Twisted Bilayer Graphene) and
(Electron-$K$-Phonon Interaction In Twisted Bilayer Graphene) with
$e_{Y},e^{\prime}_{Y}$ indices and $u^{n}_{\mathbf{k},\mathbf{Q},\alpha,\eta}$
in Eq. (Electron-$K$-Phonon Interaction In Twisted Bilayer Graphene) with
$u^{e_{Y}}_{\mathbf{k},\mathbf{Q},\alpha,\eta}$. Discrete symmetries can
constrain the form of the function
$G_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{Q}_{-l\eta}}^{\eta
e_{Y}e_{Y}^{\prime}l}$, as discussed in SM Sec.VI.D [55]. In particular, in
the chiral limit $w_{0}=0$ one can show that
$G_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{Q}_{-l\eta}}^{\eta
e_{Y}e_{Y}^{\prime}l}=\delta_{e_{Y},e_{Y}^{\prime}}G_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{Q}_{-l\eta}}^{\eta
e_{Y}e_{Y}l}$ has diagonal form on the Chern-band basis.

Phonon-mediated Electron-electron Interaction and Symmetry Classification of
Superconducting Pairing Channels - We next apply the Schrieffer-Wolff
transformation [65] to integrate out the phonon modes and obtain the phonon-
mediated electron-electron (el-el) interaction [25, 32]. We focus on the
Cooper pair channel of the attractive interaction, which takes the form

$\displaystyle
H_{ee}=-\frac{1}{N_{M}}\sum_{\mathbf{k},\mathbf{k}^{\prime},s,s_{1},e_{Y},e_{Y}^{\prime}}V^{\eta,e_{Y},e_{Y}^{\prime}}_{\mathbf{k},\mathbf{k}^{\prime}}$
$\displaystyle\gamma_{\mathbf{k}e_{Y}\eta
s}^{\dagger}\gamma_{-\mathbf{k}e^{\prime}_{Y},-\eta
s_{1}}^{\dagger}\gamma_{-\mathbf{k}^{\prime}e^{\prime}_{Y},\eta
s_{1}}\gamma_{\mathbf{k}^{\prime}e_{Y},-\eta s},$ (7)

where

$V^{\eta,e_{Y},e_{Y}^{\prime}}_{\mathbf{k},\mathbf{k}^{\prime}}=\frac{1}{N_{0}\omega_{A_{1}}}\sum_{\mathbf{G}_{M},l}G_{\mathbf{k},\mathbf{k}^{\prime},-l\eta\mathbf{q}_{2}+\mathbf{G}_{M}}^{\eta,e_{Y},l}G_{-\mathbf{k},-\mathbf{k}^{\prime},l\eta\mathbf{q}_{2}-\mathbf{G}_{M}}^{-\eta,e_{Y}^{\prime},l}$
(8)

with $\mathbf{G}_{M}$ the Moiré reciprocal lattice vectors, $N_{M}$ the number
of Moiré unit cells and $N_{0}$ the number of atomic unit cells in one Moiré
unit cell ($N_{G}=N_{0}\times N_{M}$). Discrete symmetries constrain the form
of the interaction parameter
$V^{\eta,e_{Y},e_{Y}^{\prime}}_{\mathbf{k},\mathbf{k}^{\prime}}$. The ones
leaving the momentum $(\mathbf{k},\mathbf{k}^{\prime})$ unchanged are: (1)
$\hat{C}_{2z}\hat{P}$:
$V^{\eta,e_{Y},e_{Y}^{\prime}}_{\mathbf{k},\mathbf{k}^{\prime}}=V^{-\eta,e_{Y},e_{Y}^{\prime}}_{\mathbf{k},\mathbf{k}^{\prime}}$;
(2) $\hat{C}_{2z}\hat{T}$:
$V^{\eta,e_{Y},e_{Y}^{\prime}}_{\mathbf{k},\mathbf{k}^{\prime}}=V^{\eta,-e_{Y},-e_{Y}^{\prime}\star}_{\mathbf{k},\mathbf{k}^{\prime}}$;
and (3) the combination of index reshuffling and $\hat{P}$ symmetry:
$V^{\eta,e_{Y}e_{Y}^{\prime}}_{\mathbf{k},\mathbf{k}^{\prime}}=V^{-\eta,e_{Y}^{\prime}e_{Y}}_{\mathbf{k},\mathbf{k}^{\prime}}$.
These three symmetry operations reduce the number of the independent
components of the $V$-function for a fixed $(\mathbf{k},\mathbf{k}^{\prime})$
from 8 complex parameters to 1 real
($V^{+,+-}_{\mathbf{k},\mathbf{k}^{\prime}}$) and 1 complex parameter
($V^{+,++}_{\mathbf{k},\mathbf{k}^{\prime}}$). Other discrete symmetries,
including $\hat{P}$, reshuffling, hermicity, $\hat{C}_{3z}$ and
$\hat{C}_{2z}$, relate the $V$-function at different
$(\mathbf{k},\mathbf{k}^{\prime})$. In particular, $\hat{C}_{3z}$ guarantees
$V^{\eta,e_{Y},e_{Y}}_{\mathbf{K}_{M},0}=0$ for the intra-Chern-band channels.
The projected Coulomb interaction into the flat bands of the BM model
possesses a large $U(4)\times U(4)$ spin-valley continuous symmetry [57, 18,
66]. the el-el interaction (Electron-$K$-Phonon Interaction In Twisted Bilayer
Graphene) breaks this symmetry down to the $U(2)_{e_{Y}=+}\times
U(2)_{e_{Y}=-}$ in the chiral limit and further to a total spin $SU(2)$
together with a valley charge $U(1)\otimes U(1)$ (SM Sec.VI.E [55]).

At the mean field level, the attractive interaction (Electron-$K$-Phonon
Interaction In Twisted Bilayer Graphene) is decomposed into the fermion
bilinear form $H_{\Delta}=\hat{\Delta}+\hat{\Delta}^{\dagger}$ with

$\hat{\Delta}=\sum\gamma_{\mathbf{k},e_{Y_{1}},\eta,s_{1}}^{\dagger}\Delta^{\eta}_{\mathbf{k};e_{Y_{1}}s_{1},e_{Y_{2}}s_{2}}\gamma_{-\mathbf{k},e_{Y_{2}},-\eta,s_{2}}^{\dagger},$
(9)

where the summation above includes the indices
$\mathbf{k},e_{Y_{1}},e_{Y_{2}},s_{1},s_{2},\eta$ and the gap function

$\Delta^{\eta}_{\mathbf{k};e_{Y_{1}}s_{1},e_{Y_{2}}s_{2}}=-\frac{1}{N_{M}}\sum_{\mathbf{k}^{\prime}}V_{\mathbf{k}\mathbf{k}^{\prime}}^{\eta
e_{Y_{1}}e_{Y_{2}}}\langle\gamma_{-\mathbf{k}^{\prime}e_{Y_{2}}\eta
s_{2}}\gamma_{\mathbf{k}^{\prime}e_{Y_{1}}-\eta s_{1}}\rangle.$ (10)

Since the interaction $V$-function does not involve spin, we can decompose
$\Delta^{\eta}_{\mathbf{k};e_{Y_{1}}s_{1},e_{Y_{2}}s_{2}}=\sum_{S,M}\Delta^{\eta,SM}_{\mathbf{k};e_{Y_{1}}e_{Y_{2}}}{\mathcal{S}}^{SM}_{s_{1}s_{2}}$,
where $S=0$ for spin singlet and $S=1$ ($M=-S,...,S$) for spin triplet (SM
Sec.VI.G.1 [55]).

The gap function can be classified according to the discrete symmetries. The
$C_{6v}$ group includes four 1D irreducible representations (irreps), e.g.
$A_{1,2}$ and $B_{1,2}$, and two 2D irreps, $E_{1,2}$. 1D irreps $A_{1,2}$ and
$B_{1,2}$ channels differ by their $\hat{C}_{2z}$ eigen-values,
$\lambda_{C_{2z}}=+1$ for $A_{1,2}$ and $\lambda_{C_{2z}}=-1$ for $B_{1,2}$.
Combining $\hat{C}_{2z}$ and reshuffling symmetries leads to
$\Delta^{\eta}_{\mathbf{k};e_{Y_{1}},e_{Y_{2}}}=\lambda_{C_{2z}}\Delta^{\eta}_{\mathbf{k};e_{Y_{2}},e_{Y_{1}}}$
for spin singlet and
$\Delta^{\eta}_{\mathbf{k};e_{Y_{1}},e_{Y_{2}}}=-\lambda_{C_{2z}}\Delta^{\eta}_{\mathbf{k};e_{Y_{2}},e_{Y_{1}}}$
for spin triplet. Thus, for intra-Chern-band pairing ($e_{Y_{1}}=e_{Y_{2}}$),
the $A_{1,2}$ channel must be spin singlet while the $B_{1,2}$ channel must be
spin triplet. Furthermore, the rotation $\hat{C}_{3z}$ ensures the existence
of nodes at $\mathbf{K}_{M}$ for the gap function of any 1D irrep intra-Chern-
band channel ($\Delta^{\eta}_{\mathbf{K}_{M};e_{Y},e_{Y}}=0$), while the
inter-Chern-band channel does not have such constraint. The 2D irreps $E_{1}$
and $E_{2}$ have different $\hat{C}_{2z}$ eigen-values, $\lambda_{C_{2z}}=+1$
for $E_{2}$ and $\lambda_{C_{2z}}=-1$ for $E_{1}$, similarly to the 1D irrep
case. Consequently, the $E_{2}$ channel must be spin singlet while the $E_{1}$
channel must be spin triplet for intra-Chern-band pairings. $\hat{C}_{3z}$
guarantees nodes at $\Gamma_{M}$ for both intra- and inter-Chern-band
channels, and it requires additional nodes at $\mathbf{K}_{M}$ for the inter-
Chern-band channels for both 2D $E_{1,2}$ pairings. Besides discrete
symmetries, the continuous $U(2)_{e_{Y}=1}\times U(2)_{e_{Y}=-1}$ spin
symmetry in the chiral limit guarantees the singlet and triplet pairings of
inter-Chern-band channel to be degenerate in the chiral flat band limit. The
full symmetry analysis of the gap functions can be found in SM Sec.VI.G [55].

Figure 1:  (a) Phonon dispersion of graphene. The irreps for phonon modes at
$\Gamma$ and $\mathbf{K}_{D}$ are labelled. Inset: BZ of graphene. (b) MBZ of
TBG. (c) and (d) shows the momentum dependence of the normalized gap function
$|\Delta_{\mathbf{k}}|$ for the inter-Chern-band $A_{1}$ singlet (or the
$A_{2}$ triplet) channel and intra-Chern-band pairing 2D $E_{2}$ singlet
channel, respectively. The inset in (c) shows the $|\Delta_{\mathbf{k}}|$
along the dashed line $k_{y}=0$ for both the inter-Chern-band (red) and intra-
Chern-band (blue) channels. The momenta ${\bf\Gamma}_{M},\mathbf{K}_{M},{\bf
M}_{M}$ are labelled in MBZ in (b) and (d).

Gap Equations and Self-consistent Solution of Pairing Channels - The
linearized gap equation (LGE) for the attractive interaction
(Electron-$K$-Phonon Interaction In Twisted Bilayer Graphene) can be derived
by evaluating $\langle\gamma_{-\mathbf{k}^{\prime}e_{Y_{2}}\eta
s_{2}}\gamma_{\mathbf{k}^{\prime}e_{Y_{1}}-\eta s_{1}}\rangle$ in Eq. (10) and
expanding it to linear order of the gap function. In the chiral flat band
limit, e.g. the band width is much smaller than the critical temperature
$T_{c}$, the LGE is derived as

$2k_{B}T\Delta^{\eta,SM}_{\mathbf{k};e_{Y_{1}}e_{Y_{2}}}=\frac{1}{N_{M}}\sum_{\mathbf{k}^{\prime}}V^{\eta
e_{Y_{1}}e_{Y_{2}}}_{\mathbf{k},\mathbf{k}^{\prime}}\Delta^{-\eta,SM}_{\mathbf{k}^{\prime};e_{Y_{1}}e_{Y_{2}}}.$
(11)

This is eigen-equation problem for the matrix $V^{\eta
e_{Y_{1}}e_{Y_{2}}}_{\mathbf{k},\mathbf{k}^{\prime}}$: the $T_{c}$ is
determined by the largest eigen-value and the symmetry of the gap function is
determined by that of its eigenvector. As mentioned, the only two independent
components of the $V$-function (complex
$V^{+++}_{\mathbf{k},\mathbf{k}^{\prime}}$ and real
$V^{++-}_{\mathbf{k},\mathbf{k}^{\prime}}$) leads to two independent LGEs for
the intra- and inter-Chern-band channels, respectively. The form of the LGE
suggests that all the gap functions are doubly degenerate at $T_{c}$ in the
flat-band limit. They belong either to two degenerate 1D irreps or one 2D
irrep. We first numerically solve these two LGEs from Eq. (359), and find the
forms of the gap functions with the largest eigen-values, as shown in Fig. 1.
Our numerical calculations show $k_{B}T_{c}\sim 0.21meV$ for the inter-Chern-
band channel, slightly larger than $k_{B}T_{c}\sim 0.16meV$ for the intra-
Chern-band channel. For the inter-Chern-band channels, the gap function is
almost a constant in Fig. 1a, featuring a fully gapped s-wave pairing with
even $\hat{C}_{2z}$-parity ($A_{1}$ or $A_{2}$ irrep). In the chiral flat-band
limit, spin singlet and triplet pairings are degenerate, as required by the
continuous $U(2)\times U(2)$ spin symmetry (SM Sec.VI.E [55]). Including
kinetic energy splits this degeneracy and makes the spin singlet $A_{1}$ irrep
channel to have the highest $T_{c}$. For the intra-Chern-band channel, one can
see nodes appear at the $\Gamma_{M}$ in Fig. 1b. As our previous symmetry
analysis shows that the gap function should have nodes at $\mathbf{K}_{M}$ for
the 1D irrep ($A_{1,2},B_{1,2}$) and $\Gamma_{M}$ for the 2D irrep
($E_{1,2}$), numerical results should correspond to a 2D irrep. Numerically
analyzing the symmetry property of the gap function suggests that the intra-
Chern-band channel belongs to the 2D $E_{2}$ irrep with spin singlet. Full
numerical results are discussed in SM Sec.VI.H.2 and 3 [55].

Our results for the intra-Chern-band channels reveal a d-wave character of the
gap. Using the heavy fermion formalism of TBG [46, 64], we analytically obtain
$V^{\eta,e_{Y},e_{Y}}_{\mathbf{k},\mathbf{k}^{\prime}}$:

$\displaystyle
V^{\eta,e_{Y},e_{Y}}_{\mathbf{k},\mathbf{k}^{\prime}}=U^{*}_{e_{Y},\mathbf{k}}U_{e_{Y},\mathbf{k}^{\prime}};\;\;U_{e_{Y},\mathbf{k}}=\frac{\sqrt{V_{0}}}{k^{2}+b^{2}}k_{e_{Y}}^{2},$
(12)

with $k_{e_{Y}}=k_{x}+ie_{Y}k_{y}$ ($e_{Y}=\pm$). This interaction allows us
to solve the LGE analytically to obtain the $T_{c}$,

$k_{B}T_{c}=\frac{\tilde{V}_{0}}{2},\;\;\tilde{V}_{0}=\frac{1}{N_{M}}\sum_{\mathbf{k}}V_{0}\frac{k^{4}}{(k^{2}+b^{2})^{2}},$
(13)

where $V_{0}$ and $b$ are material dependent parameters. The corresponding
self-consistent gap function takes the d-wave form

$\displaystyle\left(\begin{array}[]{c}\Delta^{+,00}_{\mathbf{k},e_{Y}e_{Y}}\\\
\Delta^{-,00}_{\mathbf{k},e_{Y}e_{Y}}\end{array}\right)=\Delta_{e_{Y}}\frac{k_{-e_{Y}}^{2}}{k^{2}+b^{2}}\left(\begin{array}[]{c}1\\\
1\end{array}\right)$ (18)

with $e_{Y}=\pm$ and $\Delta_{e_{Y}}$ a parameter to be determined. Time
reversal, if exists, requires
$\left(\Delta^{-,00}_{\mathbf{k};--},\Delta^{+,00}_{\mathbf{k};--}\right)=\left(\Delta^{+,00}_{\mathbf{k};++},\Delta^{-,00}_{\mathbf{k};++}\right)^{*}$.
The d-wave nature of the gap function suggests the possibility of the nodal
superconductivity. However, one should note that the single-particle
Hamiltonian is not diagonal in the Chern-band basis. The Bogoliubov-de Gennes
(BdG) spectrum must be checked with kinetic energy added. The BdG Hamiltonian
for the intra-Chern-band pairing is block diagonal and one block
$\mathcal{H}_{BdG}^{+,+}$ on the basis
$(\gamma^{\dagger}_{\mathbf{k},e_{Y}=\pm,+,s=\uparrow},\gamma_{-\mathbf{k},e_{Y}=\pm,-,s=\downarrow})$
reads

$\displaystyle
H_{BdG}^{+,+}(\mathbf{k})=\begin{pmatrix}h_{+}(\mathbf{k})&\Delta^{+}_{\mathbf{k}}\\\
(\Delta^{+}_{\mathbf{k}})^{\dagger}&-h_{-}^{\star}(-\mathbf{k})\end{pmatrix}$
(19)

with
$h_{\eta}(\mathbf{k})=(d_{0,\eta}(\mathbf{k})-\mu)\zeta^{0}+d_{x,\eta}(\mathbf{k})\zeta^{x}$
and
$\Delta^{+}_{\mathbf{k}}=Diag[\Delta^{+,00}_{\mathbf{k},++},\Delta^{+,00}_{\mathbf{k},--}]$.
Here
$d_{0,\eta}(\mathbf{k})=(\epsilon_{+,\eta}(\mathbf{k})+\epsilon_{-,\eta}(\mathbf{k}))/2$
and
$d_{x,\eta}(\mathbf{k})=(\epsilon_{+,\eta}(\mathbf{k})-\epsilon_{-,\eta}(\mathbf{k}))/2$,
where $\epsilon_{\pm,\eta}(\mathbf{k})$ are the eigen-energies for the two
low-energy flat bands (per valley per spin) of the BM model ${\hat{H}}_{\rm
el}$ (2). The corresponding energy spectrum can possess nodes when the pairing
amplitudes of two Chern-band channels are equal,
$|\Delta_{e_{Y}=+}|=|\Delta_{e_{Y}=-}|=\Delta_{0}$, which corresponds to the
Euler pairing discussed in Ref.[67, 32]. Point nodes appear at the location
defined by two conditions (1)
$\cos((\Phi_{\mathbf{k},-}-\Phi_{\mathbf{k},+})/2)=0$, where
$\Phi_{\mathbf{k},e_{Y}}=\varphi_{e_{Y}}-2e_{Y}\theta_{\mathbf{k}}$ with
$\Delta_{e_{Y}}=\Delta_{0}e^{i\varphi_{e_{Y}}}$ and
$k_{e_{Y}}=ke^{ie_{Y}\theta_{\mathbf{k}}}$; and (2)
$d_{x,\mathbf{k}}^{2}=(d_{0,\mathbf{k}}-\mu)^{2}+\Delta_{0,\mathbf{k}}^{2}$
with $\Delta_{0,\mathbf{k}}=\Delta_{0}\frac{k^{2}}{k^{2}+b^{2}}$, as discussed
in SM Sec.VI.H.5 [55]. The first condition determines the momentum angle for
the nodes while the second gives the momentum amplitude, thus together fixing
the location of point nodes in the 2D momentum space. We next solve the self-
consistent gap equation at zero temperature for the interaction form (415).
With the gap function ansatz
$\Delta_{\mathbf{k};e_{Y}}=\Delta_{e_{Y}}\frac{k^{2}_{-e_{Y}}}{k^{2}+b^{2}}$,
we find a self-consistent gap equation

$\Delta_{e_{Y}}=\frac{V_{0}}{N_{M}}\sum_{\mathbf{k}^{\prime},e_{Y_{1}}}\frac{{k^{\prime}}^{2}_{e_{Y}}}{k^{\prime
2}+b^{2}}u_{-\mathbf{k}^{\prime},e_{Y}e_{Y_{1}}}w^{*}_{-\mathbf{k}^{\prime},e_{Y}e_{Y_{1}}},$
(20)

where
$\psi_{\mathbf{k},e_{Y_{1}}}=(u_{\mathbf{k},\pm,e_{Y_{1}}},w_{\mathbf{k},\pm
e_{Y_{1}}})$ ($e_{Y_{1}}=\pm$) are the eigen-wave functions with the positive
eigen-energies of the BdG Hamiltonian $H_{BdG}^{+,+}(\mathbf{k})$ (19). Fig.
2A shows the chemical potential dependence of the gap functions and the
condensation energy. The Euler pairing $|\Delta_{+}|=|\Delta_{-}|$ is always
energetically favored for a non-flat band width $\sim 0.3meV$, quite different
from chiral d-wave pairing in doped graphene [68, 69, 26]. For the chemical
potential $\mu$ below $0.1meV$, a nodal superconductor phase with four point
nodes (Fig. 2B) located at the positions determined by two conditions
discussed above [32]. With increasing $\mu$, four nodes move towards
$\Gamma_{M}$ and eventually a gapped superconductor phase (Fig. 2C) appears
for $\mu>0.1meV$.

Figure 2:  (A) The superconductor order parameter amplitudes $|\Delta_{\pm}|$
(red circles and blue crosses) and the ground state energy (black dots) as a
function of $\mu$. The superconducting phase has nodes in the shadowed regime.
(B) and (C) show the BdG spectrum with and without nodes at $\mu=0.04meV$ and
$0.14meV$, respectively. The single-particle bandwidth is set around $0.3meV$.

The energy scale of the Coulomb interaction in TBG is $\sim 24meV$ [57], much
larger than the estimated energy scale of e-K-ph mediated attractive
interaction $\sim 0.3meV$ [70]. Near the Van-Hove singularities of flat bands,
the screening can significantly reduce the Coulomb interaction to a similar
order as e-K-ph mediated interaction due to the large DOS (SM Sec.VI.H.7
[55]), thus making superconductivity from this mechanism possible. If,
however, the DOS is that of the Hartree-Fock bands of correlated insulators,
the screening might not be enough to reduce the Coulomb interaction. Hence,
superconductivity from this $K$-phonon flat bare band mechanism could appear
only when the correlated insulator states are suppressed.

Conclusion -  In conclusion, we develop a theory for the projected e-K-ph
interaction of the flat bands and the resulting superconductor pairing
channels in TBG. We find the inter-Chern-band s-wave singlet pairing and the
intra-Chern-band d-wave nematic singlet pairing have highest $T_{c}$, and the
$T_{c}$ of inter-Chern-band channel is slightly higher than the intra-Chen-
band channel. The intra-Chen-band channel can have nodes in a large parameter
regime. From the estimate of the screened Coulomb interaction, we argue that
this mechanism requires the correlated insulators to be suppressed.

Acknowledgement – We would like to acknowledge Biao Lian, Xi Dai and Zhida
Song for the helpful discussion. B.A.B.’s heavy fermion in twisted bilayer
research was supported by the DOE Grant No. DE-SC0016239. BAB’s sabbatical
support also comes from the Simons Investigator Grant No. 404513, the Gordon
and Betty Moore Foundation through Grant No.GBMF8685 towards the Princeton
theory program, the Gordon and Betty Moore Foundation’s EPiQS Initiative
(Grant No. GBMF11070), Office of Naval Research (ONR Grant No.
N00014-20-1-2303), and the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program (grant agreement no.
101020833). CXL also acknowledges the support through the Penn State
MRSEC–Center for Nanoscale Science via NSF award DMR-2011839. AY acknowledges
support from the Gordon and Betty Moore Foundation’s EPiQS initiative grant
GBMF9469, DOE- BES grant DE-FG02-07ER46419, NSF-DMR-1904442, ARO MURI
(W911NF-21-2-0147), and ONR N00012- 21-1-2592. CXL and AY also acknowledges
the support from the NSF-MERSEC (Grant No. MER- SEC DMR 2011750). YLC
acknowledges the support from the Oxford-ShanghaiTech collaboration project
and the Shanghai Municipal Science and Technology Major Project (grant
2018SHZDZX02).

## References

  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, _et al._ , Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, _et al._ , Nature 574, 653 (2019).
  * Codecido  _et al._ [2019] E. Codecido, Q. Wang, R. Koester, S. Che, H. Tian, R. Lv, S. Tran, K. Watanabe, T. Taniguchi, F. Zhang, _et al._ , Science Advances 5, eaaw9770 (2019).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Cao  _et al._ [2021] Y. Cao, D. Rodan-Legrain, J. M. Park, N. F. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, science 372, 264 (2021).
  * Oh  _et al._ [2021] M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, X. Liu, K. Watanabe, T. Taniguchi, and A. Yazdani, Nature 600, 240 (2021).
  * Liu  _et al._ [2021] X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. Li, Science 371, 1261 (2021).
  * Saito  _et al._ [2020] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, Nature Physics 16, 926 (2020).
  * Arora  _et al._ [2020] H. S. Arora, R. Polski, Y. Zhang, A. Thomson, Y. Choi, H. Kim, Z. Lin, I. Z. Wilson, X. Xu, J.-H. Chu, _et al._ , Nature 583, 379 (2020).
  * Andrei and MacDonald [2020] E. Y. Andrei and A. H. MacDonald, Nature materials 19, 1265 (2020).
  * Liu and Dai [2021] J. Liu and X. Dai, Nature Reviews Physics 3, 367 (2021).
  * Kwan  _et al._ [2021] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck, Physical Review X 11, 041063 (2021).
  * Nuckolls  _et al._ [2023] K. P. Nuckolls, R. L. Lee, M. Oh, D. Wong, T. Soejima, J. P. Hong, D. Călugăru, J. Herzog-Arbeitman, B. A. Bernevig, K. Watanabe, _et al._ , arXiv preprint arXiv:2303.00024 (2023).
  * Călugăru  _et al._ [2022] D. Călugăru, N. Regnault, M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, A. Yazdani, O. Vafek, and B. A. Bernevig, Physical review letters 129, 117602 (2022).
  * Kang  _et al._ [2021] J. Kang, B. A. Bernevig, and O. Vafek, Physical review letters 127, 266402 (2021).
  * Kang and Vafek [2018] J. Kang and O. Vafek, Physical Review X 8, 031088 (2018).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Physical review letters 122, 246401 (2019).
  * Vafek and Kang [2020] O. Vafek and J. Kang, Physical Review Letters 125, 257602 (2020).
  * Khalaf  _et al._ [2021] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Science advances 7, eabf5299 (2021).
  * Khalaf  _et al._ [2022] E. Khalaf, P. Ledwith, and A. Vishwanath, Physical Review B 105, 224508 (2022).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Physical Review X 8, 031089 (2018).
  * You and Vishwanath [2019] Y.-Z. You and A. Vishwanath, npj Quantum Materials 4, 1 (2019).
  * Lian  _et al._ [2019] B. Lian, Z. Wang, and B. A. Bernevig, Physical review letters 122, 257002 (2019).
  * Wu  _et al._ [2018] F. Wu, A. H. MacDonald, and I. Martin, Physical review letters 121, 257001 (2018).
  * Wu [2019] F. Wu, Physical Review B 99, 195114 (2019).
  * Wu  _et al._ [2019] F. Wu, E. Hwang, and S. D. Sarma, Physical Review B 99, 165112 (2019).
  * Wu and Sarma [2019] F. Wu and S. D. Sarma, Physical Review B 99, 220507 (2019).
  * Chichinadze  _et al._ [2020] D. V. Chichinadze, L. Classen, and A. V. Chubukov, Physical Review B 101, 224513 (2020).
  * Fernandes and Fu [2021] R. M. Fernandes and L. Fu, Physical review letters 127, 047001 (2021).
  * Wang  _et al._ [2021] Y. Wang, J. Kang, and R. M. Fernandes, Physical Review B 103, 024506 (2021).
  * Yu  _et al._ [2022a] J. Yu, M. Xie, F. Wu, and S. D. Sarma, arXiv preprint arXiv:2202.02353 (2022a).
  * Sharma  _et al._ [2020] G. Sharma, M. Trushin, O. P. Sushkov, G. Vignale, and S. Adam, Physical Review Research 2, 022040 (2020).
  * Isobe  _et al._ [2018] H. Isobe, N. F. Yuan, and L. Fu, Physical Review X 8, 041041 (2018).
  * Roy and Juričić [2019] B. Roy and V. Juričić, Physical Review B 99, 121407 (2019).
  * Peltonen  _et al._ [2018] T. J. Peltonen, R. Ojajärvi, and T. T. Heikkilä, Physical Review B 98, 220504 (2018).
  * Kennes  _et al._ [2018] D. M. Kennes, J. Lischner, and C. Karrasch, Physical Review B 98, 241407 (2018).
  * Gonzalez and Stauber [2019] J. Gonzalez and T. Stauber, Physical review letters 122, 026801 (2019).
  * Cea and Guinea [2021] T. Cea and F. Guinea, Proceedings of the National Academy of Sciences 118, e2107874118 (2021).
  * Shavit  _et al._ [2021] G. Shavit, E. Berg, A. Stern, and Y. Oreg, Physical Review Letters 127, 247703 (2021).
  * Liu  _et al._ [2018] C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Physical review letters 121, 217001 (2018).
  * Fidrysiak  _et al._ [2018] M. Fidrysiak, M. Zegrodnik, and J. Spałek, Physical Review B 98, 085436 (2018).
  * Chou  _et al._ [2019] Y.-Z. Chou, Y.-P. Lin, S. D. Sarma, and R. M. Nandkishore, Physical Review B 100, 115128 (2019).
  * Kozii  _et al._ [2019] V. Kozii, H. Isobe, J. W. Venderbos, and L. Fu, Physical Review B 99, 144507 (2019).
  * [45] Cheng Chen, Kevin P. Nuckolls, Shuhan Ding, Wangqian Miao, Dillon Wong, Myungchul Oh, Ryan L. Lee, Shanmei He, Cheng Peng, Ding Pei, Yiwei Li, Shihao Zhang, Jianpeng Liu, Zhongkai Liu, Chris Jozwiak, Aaron Bostwick, Eli Rotenberg, Xi Dai, Chaoxing Liu, B. Andrei Bernevig, Yao Wang, Ali Yazdani, Yulin Chen, Strong Intervalley Electron-Phonon Coupling in Magic-Angle Twisted Bilayer Graphene, submitted.
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Physical Review Letters 129, 047601 (2022).
  * Călugăru  _et al._ [2023] D. Călugăru, M. Borovkov, L. L. Lau, P. Coleman, Z.-D. Song, and B. A. Bernevig, arXiv preprint arXiv:2303.03429 (2023).
  * Hu  _et al._ [2023a] H. Hu, G. Rai, L. Crippa, J. Herzog-Arbeitman, D. Călugăru, T. Wehling, G. Sangiovanni, R. Valenti, A. Tsvelik, and B. Bernevig, arXiv preprint arXiv:2301.04673 (2023a).
  * Hu  _et al._ [2023b] H. Hu, B. A. Bernevig, and A. M. Tsvelik, arXiv preprint arXiv:2301.04669 (2023b).
  * Lau and Coleman [2023] L. L. Lau and P. Coleman, arXiv preprint arXiv:2303.02670 (2023).
  * Datta  _et al._ [2023] A. Datta, M. J. Calderon, A. Camjayi, and E. Bascones, arXiv preprint arXiv:2301.13024 (2023).
  * Yu  _et al._ [2023] J. Yu, M. Xie, B. A. Bernevig, and S. D. Sarma, arXiv preprint arXiv:2301.04171 (2023).
  * Chou and Sarma [2022] Y.-Z. Chou and S. D. Sarma, arXiv preprint arXiv:2211.15682 (2022).
  * Zhou and Song [2023] G.-D. Zhou and Z.-D. Song, arXiv preprint arXiv:2301.04661 (2023).
  * sm [2] See Supplemental Material at [URL] for details of the derivations of electron-phonon interaction in graphene and its projection to the flat bands of twisted bilayer graphene, symmetry analysis of phonon-mediated electron-electron interaction and irreducible superconductor pairing channels, the linearized and T=0 gap equations, the superconductor ground state energy and the Coulomb screening effect.
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Bernevig  _et al._ [2021] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Physical Review B 103, 205413 (2021).
  * Song  _et al._ [2021] Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Physical Review B 103, 205412 (2021).
  * Neto  _et al._ [2009] A. C. Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim, Reviews of modern physics 81, 109 (2009).
  * Sahoo and Mishra [2012] R. Sahoo and R. R. Mishra, Journal of Experimental and Theoretical Physics 114, 805 (2012).
  * Thingstad  _et al._ [2020] E. Thingstad, A. Kamra, J. W. Wells, and A. Sudbø, Physical Review B 101, 214513 (2020).
  * Maultzsch  _et al._ [2004] J. Maultzsch, S. Reich, C. Thomsen, H. Requardt, and P. Ordejón, Physical review letters 92, 075501 (2004).
  * Mohr  _et al._ [2007] M. Mohr, J. Maultzsch, E. Dobardžić, S. Reich, I. Milošević, M. Damnjanović, A. Bosak, M. Krisch, and C. Thomsen, Physical Review B 76, 035439 (2007).
  * [64] Liu, Chao-Xing and Bernevig, B Andrei, in preparation.
  * Schrieffer and Wolff [1966] J. R. Schrieffer and P. A. Wolff, Physical Review 149, 491 (1966).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Physical Review X 10, 031034 (2020).
  * Yu  _et al._ [2022b] J. Yu, Y.-A. Chen, and S. D. Sarma, Physical Review B 105, 104515 (2022b).
  * Black-Schaffer and Honerkamp [2014] A. M. Black-Schaffer and C. Honerkamp, Journal of Physics: Condensed Matter 26, 423201 (2014).
  * Nandkishore  _et al._ [2012] R. Nandkishore, L. S. Levitov, and A. V. Chubukov, Nature Physics 8, 158 (2012).
  * [70] Note that our calculation is performed in the chiral limit, so we may underestimate the electron-phonon coupling strength.
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, Physical review letters 123, 036401 (2019).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Physical review letters 122, 106405 (2019).
  * Ahn  _et al._ [2019] J. Ahn, S. Park, and B.-J. Yang, Physical Review X 9, 021013 (2019).
  * Po  _et al._ [2019] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Physical Review B 99, 195455 (2019).
  * Savary  _et al._ [2017] L. Savary, J. Ruhman, J. W. Venderbos, L. Fu, and P. A. Lee, Physical Review B 96, 214514 (2017).
  * Sigrist and Ueda [1991] M. Sigrist and K. Ueda, Reviews of Modern physics 63, 239 (1991).

Supplementary Material for ”Electron-$K$-Phonon Interaction In Twisted Bilayer
Graphene”

###### Contents

  1. I Notation
  2. II Electron-phonon interaction in a single layer graphene
    1. II.1 Zero Displacement
    2. II.2 Finite Out of-Plane Displacement
    3. II.3 Finite In-Plane Displacement
    4. II.4 Inter-valley Electron-Phonon Interaction with In-plane Phonons
      1. II.4.1 Physical Meaning of the phonon Modes at $\mathbf{K}_{D}$
    5. II.5 Symmetry properties of the Electron-Phonon Interaction
  3. III Phonons in a Single Layer Graphene
    1. III.1 Dynamical Matrix
    2. III.2 Dynamical Matrix and Phonon Modes at $\mathbf{K}_{D}$
  4. IV Electron and Phonon-Band Projected Electron-Phonon Hamiltonian
    1. IV.1 Projected Inter-Valley Electron-Phonon Interaction
    2. IV.2 Numerical Calculations of Electron-Phonon Interaction in a Single-Layer Graphene
  5. V Review of Bistritzer-MacDonald model for Twisted Bilayer Graphene
    1. V.1 Basis and Single-particle Hamiltonian
    2. V.2 Symmetries of the BM model
  6. VI Intra-Layer Inter-Valley Electron-phonon Interaction in Twisted Bilayer Graphene
    1. VI.1 Projected Intra-Layer Inter-Valley Electron-Phonon Interaction for the Flat Bands
    2. VI.2 Phonon-Mediated Electron-Electron Interaction for the Flat Bands
    3. VI.3 Gauge Fixing for Eigen-state Basis and Chern Band Basis
      1. VI.3.1 Eigen-state Basis
      2. VI.3.2 Gauge Fixing for Eigen-state Basis
      3. VI.3.3 Gauge Fixing for Chern band basis
      4. VI.3.4 The chiral symmetry at $w_{0}=0$
    4. VI.4 Discrete Symmetry Analysis of the Intra-Layer Inter-valley Electron-Phonon Interaction and the phonon-induced electron-electron interaction
    5. VI.5 Continuous Symmetry Analysis of the Intra-layer and Inter-valley Electron-Phonon Interaction and the phonon-induced electron-electron interaction
    6. VI.6 Numerical results for the phonon-induced electron-electron interaction
    7. VI.7 Symmetry Classification of Irreducible Pairing Channels
      1. VI.7.1 Discrete Symmetry analysis of the Gap Function
      2. VI.7.2 Continuous Symmetry Analysis of Irreducible Pairing Channels
    8. VI.8 Gap equation and Possible Pairing Channels for the Flat Bands
      1. VI.8.1 Linearized gap equation
      2. VI.8.2 Intra-Chern-band channels from linearized gap equation
      3. VI.8.3 Inter-Chern-band channels from linearized gap equation
      4. VI.8.4 Bogoliubov-de Gennes Hamiltonian, full self-consistent gap equation and ground state energy
      5. VI.8.5 Bogoliubov-de Gennes Hamiltonian for the Intra-Chern-band Channel
      6. VI.8.6 Bogoliubov-de Gennes Hamiltonian for the Inter-Chern-band Channel
      7. VI.8.7 Estimate of screened Coulomb Interaction

##  I Notation

In this section, we will first describe our notations for the whole paper. We
denote ${\bf a}_{1}^{l},{\bf a}_{2}^{l}$ as the un-strained, microscopic
lattice vectors of the graphene layer $l$ ($l=\pm$ for the top and bottom
layers), while ${\bf G}_{1,2}^{l}$ are the corresponding reciprocal lattice
vectors. ${\bf G}_{a}^{l}\cdot{\bf a}_{b}^{l}=2\pi\delta_{ab}$. In the layer
$l$, there are two sub-lattice positions
${\bf\tau}_{A}^{l},{\bf\tau}_{B}^{l}$. We assume the graphene layer $l$ is
rotated by an in-plane angle $l\theta/2$ followed by a shift
$l\frac{l_{0}}{2}\hat{z}$ along the z direction, so any vector ${\bf v}^{l}$
is related to the un-rotated vector ${\bf v}$ by ${\bf
v}^{l}=\hat{R}^{l}_{\theta/2}{\bf
v}+l\frac{l_{0}}{2}\hat{z}=l\frac{\theta}{2}\hat{z}\times{\bf
v}+l\frac{l_{0}}{2}\hat{z}$, where $\hat{R}^{l}_{\theta/2}$ labels the
rotation operator, $\hat{z}$ is the out-of-plane unit vector and $l_{0}$ is
the equilibrium distance between two graphene layers. Thus, we only need to
give the values of all the un-rotated vectors below. Two primitive lattice
vectors are ${\bf a}_{1}=a_{0}(\frac{-1}{2},\frac{\sqrt{3}}{2},0)$ and ${\bf
a}_{2}=-a_{0}(\frac{1}{2},\frac{\sqrt{3}}{2},0)$, where $a_{0}=2.46\AA$ is the
lattice constant. The two primitive reciprocal lattice vectors spanning the
Brillouin zone (BZ) are

${\bf
G}_{1}=\frac{4\pi}{\sqrt{3}a_{0}}(-\frac{\sqrt{3}}{2},\frac{1}{2});\;\;{\bf
G}_{2}=\frac{4\pi}{\sqrt{3}a_{0}}(-\frac{\sqrt{3}}{2},-\frac{1}{2}).$ (21)

Since the Fermi surface in graphene is around two $K$-valley, we also label
them as $\eta\mathbf{K}_{D}$ with $\mathbf{K}_{D}=\frac{4\pi}{3a_{0}}(1,0)$
and $\eta=\pm$. The atoms are at positions

${\bf\tau}_{A}=\frac{a_{0}}{2}(1,-\frac{1}{\sqrt{3}}),\;\;{\bf\tau}_{B}=\frac{a_{0}}{2}(1,\frac{1}{\sqrt{3}}).$
(22)

For two layers $l=\pm$, the atoms on position $\alpha=A,B$ in the top/bottom
later reside at

${\bf R}^{l}_{\alpha}=\hat{R}_{\theta/2}^{l}{\bf
R}+l\frac{l_{0}}{2}\hat{z}+{\bf\tau}_{\alpha}^{l}={\bf
R}^{l}+{\bf\tau}_{\alpha}^{l}$ (23)

where ${\bf R}=m{\bf a}_{1}+n{\bf a}_{2}$ is a unit cell vector.

The hopping between two atomic centers at a distance $r$ is

$\displaystyle t({\bf r}_{\parallel},z)=V_{pp\pi}(1-(\frac{{\bf
r}\cdot\hat{z}}{r})^{2})e^{-(r-\frac{a}{\sqrt{3}})/r_{0}}+V_{pp\sigma}(\frac{{\bf
r}\cdot\hat{z}}{r})^{2}e^{-\frac{r-l_{0}}{r_{0}}}$ (24)

where ${\bf r}_{\parallel}=(x,y)$, $r_{0}$ characterizes the decaying length
of atomic orbitals of carbon atoms, $a=a_{0}/\sqrt{3}$ is the distance between
the $A$ and $B$ atoms, and $V_{pp\pi}$ and $V_{pp\sigma}$ are the
hybridization parameters of $\pi$ and $\sigma$ bonds. We performed the Fourier
transform to the in-plane coordinate of $t({\bf r}_{\parallel},z)$,

$t_{{\bf q},z}=\frac{1}{S_{0}}\int d^{2}{\bf r}_{\parallel}e^{-i{\bf
q}\cdot{\bf r}_{\parallel}}t({\bf r}_{\parallel},z)$ (25)

where $S_{0}$ is the unit cell area.

Figure 3: Lattice and Brillouin zone for graphene.

In the discussion below, we mainly focus on the in-plane hopping term, and
direct calculation gives

$t_{{\bf
q},0}=\frac{V_{pp\pi}e^{\frac{a_{0}}{\sqrt{3}r_{0}}}}{S_{0}}\frac{2\pi
r_{0}^{2}}{(1+(qr_{0})^{2})^{3/2}}$ (26)

In the calculation we have used

$r_{0}=0.184a_{0};\;V_{pp\pi}=-2.7eV;\;S_{0}=\frac{\sqrt{3}a_{0}^{2}}{2};\;|\mathbf{K}_{D}|=1.7028\AA^{-1}$
(27)

In the presence of a phonon field, the atom position should include the
displacement field ${\bf u}^{l}({\bf R}^{l}_{\alpha})$ as

${\bf R}_{\alpha}^{l}={\bf R}^{l}+{\bf\tau}_{\alpha}^{l}\rightarrow\tilde{\bf
R}_{\alpha}^{l}={\bf R}^{l}+{\bf\tau}_{\alpha}^{l}+{\bf u}^{l}({\bf
R}^{l}_{\alpha})$ (28)

where

${\bf u}^{l}({\bf R}^{l}+\tau_{\alpha}^{l})=\frac{1}{\sqrt{N_{G}}}\sum_{\bf
p}{\bf u}_{\bf p}^{l,\alpha}e^{i{\bf p}\cdot({\bf R}^{l}+\tau_{\alpha}^{l})}$
(29)

where the sum over ${\bf p}$ extends of the graphene BZ and the $N_{G}$ are
the total number of atomic unit cells in graphene. Another important property
of the phonon field is its periodic gauge. With the form of inversion Fourier
transform

$\mathbf{u}_{\mathbf{q}}^{l,\alpha}=\frac{1}{\sqrt{N_{G}}}\sum_{\mathbf{R}}\mathbf{u}^{l}(\mathbf{R}^{l}+\tau_{\alpha}^{l})e^{-i\mathbf{q}\cdot(\mathbf{R}^{l}+\tau_{\alpha}^{l})},$
(30)

we can see

$\mathbf{u}^{\alpha}_{\mathbf{q}+\mathbf{G}}=e^{-i\mathbf{G}\cdot\tau_{\alpha}}\mathbf{u}^{l\alpha}_{\mathbf{q}}$
(31)

which is called the ”periodic gauge”.

Here we only introduce our notation for the single layer graphene, and the
notation for the Moiré lattice for twisted bilayer graphene (TBG) will be
introduced in Sec. V.

##  II Electron-phonon interaction in a single layer graphene

On the basis states

$\left|\psi_{{\bf p},\alpha}^{l}\right\rangle=\frac{1}{\sqrt{N_{G}}}\sum_{{\bf
R}^{l}}e^{i{\bf p}\cdot({\bf R}^{l}+{\bf\tau}_{\alpha}^{l})}\left|\psi({\bf
R}_{\alpha}^{l}+{\bf u}^{l}({\bf R}^{l}+{\bf\tau}_{\alpha}^{l})\right\rangle,$
(32)

the matrix Hamiltonian between the momentum ${\bf p},{\bf p}^{\prime}$ in
layers $a,b=1,2$ and orbitals $\alpha,\beta=A,B$ is given by

$\displaystyle H_{ll^{\prime}}^{\alpha\alpha^{\prime}}({\bf
pp}^{\prime})=\frac{1}{N_{G}}\sum_{{\bf R}^{l},{\bf
R}^{l^{\prime}}}t(\tilde{\bf R}_{\alpha}^{l}-\tilde{\bf
R}_{\alpha^{\prime}}^{l^{\prime}})e^{-i({\bf R}_{\alpha}^{l}\cdot{\bf p}-{\bf
R}_{\alpha^{\prime}}^{l^{\prime}}\cdot{\bf p}^{\prime})}$ (33)

or

$H_{ll^{\prime}}^{\alpha\alpha^{\prime}}({\bf
pp}^{\prime})=\frac{S_{0}}{(2\pi)^{2}N_{G}}\int d^{2}{\bf q}\left(\sum_{{\bf
R}^{l},{\bf R}^{l^{\prime}}}e^{-i({\bf R}_{\alpha}^{l}\cdot({\bf p-q})-{\bf
R}_{\alpha^{\prime}}^{l^{\prime}}\cdot({\bf p}^{\prime}-{\bf q}))}e^{i{\bf
q}\cdot({\bf u}({\bf R}_{\alpha}^{l})-{\bf u}({\bf
R}_{\alpha^{\prime}}^{l^{\prime}}))}\right)t({\bf q},z_{ll^{\prime}})$ (34)

where $z_{ll^{\prime}}$ is the height difference between the layers $l$ and
$l^{\prime}$. In this section, we first discuss the electron-phonon (e-ph)
interaction in the same layer $l=l^{\prime}$, and may drop the layer index $l$
in some expressions for simplicity.

###  II.1 Zero Displacement

We have

$\frac{1}{N_{G}}\sum_{{\bf R}^{l}}e^{-i{\bf R}^{l}\cdot({\bf p-q})}=\sum_{{\bf
G}^{l}}\delta_{{\bf p-q},{\bf G}^{l}}$ (35)

where ${\bf G}^{l}$ is a reciprocal lattice vector in the layer $l$.

For the layer $l=l^{\prime}$, we have (since $z_{ll}=0$)

$\displaystyle H_{0,ll}^{\alpha\alpha^{\prime}}({\bf
pp}^{\prime})=\frac{S_{0}N_{G}}{(2\pi)^{2}}\int d^{2}{\bf q}\sum_{{\bf
G}_{a}^{l},{\bf G}_{b}^{l}}\left(\delta_{{\bf p-q},-{\bf
G}_{a}^{l}}\delta_{{\bf p^{\prime}-q},-{\bf G}_{b}^{l}}e^{-i{\bf
G}_{b}^{l}\cdot{\bf\tau}_{\alpha^{\prime}}^{l}+i{\bf
G}_{a}^{l}\cdot{\bf\tau}_{\alpha}^{l}}\right)t({\bf q},0)$
$\displaystyle=\sum_{{\bf G}_{a}^{l},{\bf G}_{b}^{l}}\delta_{{\bf p}+{\bf
G}_{a}^{l},{\bf p}^{\prime}+{\bf G}_{b}^{l}}e^{-i{\bf
G}_{b}^{l}\cdot{\bf\tau}_{\alpha^{\prime}}^{l}+i{\bf
G}_{a}^{l}\cdot{\bf\tau}_{\alpha}^{l}}t({\bf p+G}_{a}^{l},0)$ (36)

where ${\bf G}_{a}^{l}$, ${\bf G}_{b}^{l}$ are two sets of reciprocal lattice
vectors in the BZ of the layer $l$. The summation in Eq. (36) respects the
symmetries of the problem. When we restrict ${\bf p,p}^{\prime}$ to the same
(first) BZ of graphene, we have $\delta_{{\bf p+G}_{a}^{l},{\bf
p^{\prime}+G}_{b}^{l}}=\delta_{\bf p,p^{\prime}}\delta_{{\bf G}_{a}^{l},{\bf
G}_{b}^{l}}$ and hence

$H_{0ll}^{\alpha\alpha^{\prime}}({\bf pp}^{\prime})=\delta_{\bf
p,p^{\prime}}\sum_{{\bf G}_{a}^{l}}e^{-i{\bf
G}_{a}^{l}\cdot({\bf\tau}_{\alpha^{\prime}}^{l}-{\bf\tau}_{\alpha}^{l})}t({\bf
p+G}_{a}^{l},0)$ (37)

We then have the tight-binding Hamiltonian for the single layer graphene

$t_{\bf p}^{\alpha\alpha^{\prime}}=\sum_{{\bf G}_{a}^{l}}e^{-i{\bf
G}_{a}^{l}\cdot({\bf\tau}_{\alpha^{\prime}}^{l}-{\bf\tau}_{\alpha}^{l})}t({\bf
p+G}_{a}^{l},0),$ (38)

which has a node at the ${\bf K}_{D}$ point.

###  II.2 Finite Out of-Plane Displacement

For the out-of-plane displacement, we can use hopping $t({\bf q},z=u^{z}({\bf
R}_{\alpha}^{l})-u^{z}({\bf R}_{\alpha}^{\prime l}))$. With the small
displacement approximation, we can linearize the expression in small $u$, and
sum over graphene lattice, leading to the intra-layer Hamiltonian

$\displaystyle H_{\perp ll}^{\alpha\alpha^{\prime}}({\bf
pp}^{\prime})=\frac{S_{0}}{(2\pi)^{2}N_{G}}\int d^{2}{\bf q}\left(\sum_{{\bf
R}_{1}^{l},{\bf R}_{2}^{l}}e^{-i({\bf R}_{1\alpha}^{l}({\bf p-q})-{\bf
R}_{2\alpha^{\prime}}^{l}({\bf p^{\prime}-q}))}(u^{z}({\bf
R}_{\alpha}^{l})-u^{z}({\bf R}_{\alpha}^{\prime l}))\partial_{z}t({\bf
q},z)|_{z=0}\right)$ $\displaystyle=\frac{1}{N^{1/2}_{G}}\sum_{{\bf p}_{1}\in
BZ}\sum_{{\bf G}_{a}^{l},{\bf G}_{b}^{l}}e^{i(\tau_{\alpha}^{l}\cdot{\bf
G}_{a}^{l}-\tau_{\alpha^{\prime}}^{l}\cdot{\bf G}_{b}^{l})}\delta_{{\bf
p^{\prime}+p_{1}+G}_{b}^{l},{\bf p+G}_{a}^{l}}$ $\displaystyle\left(u_{{\bf
p}_{1}}^{zl,\alpha}\partial_{z}t({\bf G}_{a}^{l}+{\bf
p-p_{1}},z)|_{z=0}-u_{\bf p_{1}}^{zl,\alpha^{\prime}}\partial_{z}t({\bf
G}_{b}^{l}+{\bf p^{\prime}+p_{1}},z)|_{z=0}\right).$ (39)

One can show that

$\frac{\partial t({\bf r},t)}{\partial z}|_{z=0}=0$ (40)

and hence there is no out-of-plane phonon coupling to the graphene electrons
(this will also be derived by symmetry), and thus

$H_{\perp ll}^{\alpha\alpha^{\prime}}({\bf pp}^{\prime})=0$ (41)

Thus, we only focus on in-plane phonon coupling below.

###  II.3 Finite In-Plane Displacement

For the in-plane displacement, we can use the Fourier transform form of the
hopping term, $t({\bf q},z=0)$. With the small displacement approximation, we
can linearize ${\bf u}$, and sum over graphene lattice, leading to the intra-
layer Hamiltonian

$\displaystyle H_{\parallel ll}^{\alpha\alpha^{\prime}}({\bf
pp}^{\prime})=\frac{S_{0}}{(2\pi)^{2}N_{G}}\int d^{2}{\bf q}\left(\sum_{{\bf
R}_{1}^{l},{\bf R}_{2}^{l}}e^{-i({\bf R}_{1\alpha}^{l}\cdot({\bf p-q})-{\bf
R}_{2\alpha^{\prime}}^{l}({\bf p^{\prime}-q}))}i{\bf q}\cdot({\bf u}({\bf
R}_{1\alpha}^{l})-{\bf u}({\bf R}_{2\alpha^{\prime}}^{l}))\right)t({\bf q},0)$
$\displaystyle=\frac{S_{0}N^{1/2}_{G}}{(2\pi)^{2}}\int d^{2}{\bf q}\times
t({\bf q},0)i{\bf q}\cdot$ $\displaystyle\sum_{{\bf p}_{1}\in BZ}\sum_{{\bf
G}_{a}^{l},{\bf G}_{b}^{l}}\left(\delta_{{\bf-p_{1}+p-q},-{\bf
G}_{a}^{l}}\delta_{{\bf p^{\prime}-q},-{\bf
G}_{b}^{l}}e^{i({\bf\tau}_{\alpha}^{l}\cdot{\bf
G}_{a}^{l}-{\bf\tau}_{\alpha^{\prime}}^{l}\cdot{\bf G}_{b}^{l}))}{\bf u}_{{\bf
p}_{1}}^{l,\alpha}-\delta_{{\bf p_{1}+p^{\prime}-q},-{\bf
G}_{b}^{l}}\delta_{{\bf p-q},-{\bf
G}_{a}^{l}}e^{i({\bf\tau}_{\alpha}^{l}\cdot{\bf
G}_{a}^{l}-{\bf\tau}_{\alpha^{\prime}}^{l}\cdot{\bf G}_{b}^{l})}{\bf u}_{{\bf
p}_{1}}^{l,\alpha^{\prime}}\right)$ (42)

where ${\bf G}_{a}^{l},{\bf G}_{b}^{l}$ are two sets of reciprocal lattice
vectors in the the BZ of layer $l$. Transforming the $\delta$ functions

$\displaystyle\delta_{{\bf-p_{1}+p-q},-{\bf G}_{a}^{l}}\delta_{{\bf
p^{\prime}-q},-{\bf G}_{b}^{l}}=\delta_{{\bf q},\frac{{\bf G}_{a}^{l}+{\bf
G}_{b}^{l}+{\bf p+p^{\prime}-p}_{1}}{2}}\delta_{{\bf p^{\prime}+p}_{1}+{\bf
G}_{b}^{l},{\bf p+G}_{a}^{l}}$ $\displaystyle\delta_{{\bf
p_{1}+p^{\prime}-q},-{\bf G}_{b}^{l}}\delta_{{\bf p-q},-{\bf
G}_{a}^{l}}=\delta_{{\bf q},\frac{{\bf G}_{a}^{l}+{\bf G}_{b}^{l}+{\bf
p+p^{\prime}+p_{1}}}{2}}\delta_{{\bf p^{\prime}+p_{1}+G}_{b}^{l},{\bf
p+G}_{a}^{l}},$ (43)

we find the matrix elements of the intra-layer Hamiltonian,

$\displaystyle H_{\parallel l}^{\alpha\alpha^{\prime}}({\bf
pp}^{\prime})=\frac{1}{N^{1/2}_{G}}\sum_{{\bf p}_{1}\in BZ}\sum_{{\bf
G}_{a}^{l},{\bf G}_{b}^{l}}e^{i({\bf\tau}_{\alpha}^{l}\cdot{\bf
G}_{a}^{l}-{\bf\tau}_{\alpha^{\prime}}^{l}\cdot{\bf G}_{b}^{l})}\delta_{{\bf
p^{\prime}+p_{1}+G}_{b}^{l},{\bf p+G}_{a}^{l}}$
$\displaystyle\left(i\frac{{\bf G}_{a}^{l}+{\bf G}_{b}^{l}+{\bf
p+p^{\prime}-p_{1}}}{2}\cdot{\bf u}_{{\bf p}_{1}}^{l,\alpha}t(\frac{{\bf
G}_{a}^{l}+{\bf G}_{b}^{l}+{\bf p+p^{\prime}-p_{1}}}{2},0)-i\frac{{\bf
G}_{a}^{l}+{\bf G}_{b}^{l}+{\bf p+p^{\prime}+p_{1}}}{2}\cdot{\bf u}_{{\bf
p}_{1}}^{l,\alpha^{\prime}}t(\frac{{\bf G}_{a}^{l}+{\bf G}_{b}^{l}+{\bf
p+p^{\prime}+p_{1}}}{2},0)\right)=$
$\displaystyle=\frac{1}{N^{1/2}_{G}}\sum_{{\bf p}_{1}\in BZ}\sum_{{\bf
G}_{a}^{l},{\bf G}_{b}^{l}}e^{i({\bf\tau}_{\alpha}^{l}\cdot{\bf
G}_{a}^{l}-{\bf\tau}_{\alpha^{\prime}}^{l}\cdot{\bf G}_{b}^{l})}\delta_{{\bf
p^{\prime}+p_{1}+G}_{b}^{l},{\bf p+G}_{a}^{l}}$ $\displaystyle\left(i({\bf
G}_{a}^{l}+{\bf p-p_{1}})\cdot{\bf u}_{{\bf p}_{1}}^{l,\alpha}t({\bf
G}_{a}^{l}+{\bf p-p_{1}},0)-i({\bf G}_{b}^{l}+{\bf p^{\prime}+p_{1}})\cdot{\bf
u}_{{\bf p}_{1}}^{l,\alpha^{\prime}}t({\bf G}_{b}^{l}+{\bf
p^{\prime}+p_{1}},0)\right)$ (44)

It should be noted that ${\bf p}_{1}$ summation in the above expression is
over one BZ. For convenience, we can choose it within the first BZ. From the
$\delta$-function, we have ${\bf p_{1}}={\bf p+G}_{a}^{l}-{\bf
p^{\prime}-G}_{b}^{l}$, so if we keep the summation over ${\bf G}_{b}^{l}$,
only one ${\bf G}_{a}^{l}$ is allowed to make ${\bf p}_{1}$ in the first BZ.
Actually the choice of ${\bf G}_{a}^{l}$ will not influence the results. To
see that, let’s assume the choice of ${\bf G}_{a0}^{l}$ can make ${\bf
p}_{10}={\bf p+G}_{a0}^{l}-{\bf p^{\prime}-G}_{b}^{l}$ within the first BZ,
and now we consider another choice ${\bf G}_{a}^{l}={\bf
G}_{a0}^{l}+\Delta{\bf G}^{l}$, so ${\bf p_{1}=p+G}_{a}^{l}-{\bf
p^{\prime}-G}_{b}^{l}={\bf p}_{10}+\Delta{\bf G}^{l}$ is outside the first BZ.
For the first term in Eq. (44), we have

$\displaystyle e^{i(\tau_{\alpha}^{l}\cdot{\bf
G}_{a}^{l}-\tau_{\alpha^{\prime}}^{l}\cdot{\bf G}_{b}^{l})}i({\bf
G}_{a}^{l}+{\bf p-p_{1}})\cdot{\bf u}_{{\bf p}_{1}}^{l,\alpha}t({\bf
G}_{a}^{l}+{\bf p-p_{1}},0)$ $\displaystyle=e^{i(\tau_{\alpha}^{l}\cdot({\bf
G}_{a0}^{l}+\Delta{\bf G}^{l})-\tau_{\alpha^{\prime}}^{l}\cdot{\bf
G}_{b}^{l})}i({\bf G}_{a0}^{l}+{\bf p-p_{10}})\cdot{\bf u}_{{\bf
p}_{10}+\Delta{\bf G}^{l}}^{l,\alpha}t({\bf G}_{a0}^{l}+{\bf p-p_{10}},0)$
$\displaystyle=e^{i(\tau_{\alpha}^{l}\cdot{\bf
G}_{a0}^{l}-\tau_{\alpha^{\prime}}^{l}\cdot{\bf G}_{b}^{l})}i({\bf
G}_{a0}^{l}+{\bf p-p}_{10})\cdot{\bf u}_{{\bf p}_{10}}^{l,\alpha}t({\bf
G}_{a0}^{l}+{\bf p-p}_{10},0)$ (45)

where we have used the periodic gauge
$\mathbf{u}^{l,\alpha}_{\mathbf{p}_{10}+\Delta\mathbf{G}^{l}}=e^{-i\tau_{\alpha}^{l}\cdot\Delta{\bf
G}^{l}}{\bf u}_{\mathbf{p}_{10}}^{l,\alpha}$. From the above expression, one
can see that the first term in Eq. (44) remains the same if one changes from
the summation over ${\bf p}_{10}$ to that over ${\bf p}_{1}$. This means one
just needs to pick one ${\bf G}_{a}^{l}$ (e.g. one possible choice ${\bf
G}_{a}^{l}={\bf G}_{b}^{l}$), while which choice will not influence the
results. Similar analysis can also be applied to the second term in Eq. (44).

Both in graphene and TBG, the electronic states with the energies around the
Fermi energy is only close to the $\eta{\bf K}_{D}$ point where $\eta=\pm$ is
the valley. Thus, one needs to expand the electron-phonon Hamiltonian around
$\eta{\bf K}_{D}$. The electron momenta $\mathbf{p}$ and $\mathbf{p}^{\prime}$
in Eq. (II.3) can be expanded around the same valley ($\eta{\bf K}_{D}$ for a
fixed $\eta$) or opposite valleys ($\eta{\bf K}_{D}$ and $-\eta{\bf K}_{D}$),
corresponding to intra-valley and inter-valley e-ph interaction, respectively.
In this work, we focus on inter-valley e-ph interaction.

###  II.4 Inter-valley Electron-Phonon Interaction with In-plane Phonons

In this section, we choose $\mathbf{p}=\mathbf{k}+\eta\mathbf{K}_{D}$,
$\mathbf{p}^{\prime}=\mathbf{k}^{\prime}-\eta\mathbf{K}_{D}$ where
$k,k^{\prime}\ll|\mathbf{K}_{D}|$, and correspondingly, we also need to
decompose the phonon momentum into a large value (the difference between the
two valleys) and a small value around that,

$\mathbf{p}_{1}=\mathbf{p}+\mathbf{P}_{1},;\;\;p\ll P_{1}$ (46)

where $\mathbf{P}_{1}$ determines the major part of phonon momentum. The
momentum $\delta$-function in $H_{\parallel
ll}^{\alpha\alpha^{\prime}}(\mathbf{p}\mathbf{p}^{\prime})$ becomes

$\delta_{\mathbf{k}+\eta\mathbf{K}_{D}+\mathbf{G}_{a}^{l},\mathbf{k}^{\prime}-\eta\mathbf{K}_{D}+\mathbf{p}_{1}+\mathbf{G}_{b}^{l}}=\delta_{\mathbf{k},\mathbf{k}^{\prime}+\mathbf{p}}\delta_{2\eta
K_{D}+\mathbf{G}_{a}^{l}-\mathbf{G}_{b}^{l},\mathbf{P}_{1}}.$ (47)

For the momentum $\mathbf{p}$ in the first BZ, there are $3$ possible values
for $\mathbf{P}_{1}$ momenta

$\mathbf{P}_{1}=-\eta\mathbf{K}_{D};-\eta(\mathbf{K}_{D}+\mathbf{G}_{1});-\eta(\mathbf{K}_{D}+\mathbf{G}_{2}).$
(48)

However, as discussed in Sec. II.3, the summation over $\mathbf{p}_{1}$ in the
Hamiltonian (44) is within one BZ and the choice of different BZs will not
influence the result due to the periodic gauge of the phonon modes. Thus, we
make the choice of $\mathbf{P}_{1}=-\eta\mathbf{K}_{D}$ below, which
corresponds to the K-phonon, and the Hamiltonian reads

$\displaystyle H_{\parallel
l}^{\eta\alpha,-\eta\alpha^{\prime}}(\mathbf{k}\mathbf{k}^{\prime})=\frac{i}{N^{1/2}_{G}}\sum_{\mathbf{G}_{a}^{l}}e^{i(\tau_{\alpha}^{l}-\tau_{\alpha^{\prime}}^{l})\cdot\mathbf{G}_{a}^{l}}e^{i\tau_{\alpha^{\prime}}^{l}\cdot(\mathbf{P}_{1}-2\eta\mathbf{K}_{D})}$
$\displaystyle\left((\mathbf{G}_{a}^{l}+\eta\mathbf{K}_{D}+\mathbf{k}^{\prime}-\mathbf{P}_{1})\cdot\mathbf{u}_{\mathbf{P}_{1}+\mathbf{k}-\mathbf{k}^{\prime}}^{l,\alpha}t(\mathbf{G}_{a}^{l}+\eta\mathbf{K}_{D}+\mathbf{k}^{\prime}-\mathbf{P}_{1},0)-(\mathbf{G}_{a}^{l}+\eta\mathbf{K}_{D}+\mathbf{k})\cdot\mathbf{u}_{\mathbf{P}_{1}+\mathbf{k}-\mathbf{k}^{\prime}}^{l,\alpha^{\prime}}t(\mathbf{G}_{a}^{l}+\eta\mathbf{K}_{D}+\mathbf{k},0)\right)$
$\displaystyle=\frac{i}{N^{1/2}_{G}}\sum_{\mathbf{G}_{a}^{l}}e^{i(\tau_{\alpha}^{l}-\tau_{\alpha^{\prime}}^{l})\cdot\mathbf{G}_{a}^{l}}\left((\mathbf{G}_{a}^{l}+2\eta\mathbf{K}_{D}+\mathbf{k}^{\prime})\cdot\mathbf{u}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime}}^{l,\alpha}t(\mathbf{G}_{a}^{l}+2\eta\mathbf{K}_{D}+\mathbf{k}^{\prime},0)\right.$
$\displaystyle\left.-(\mathbf{G}_{a}^{l}+\eta\mathbf{K}_{D}+\mathbf{k})\cdot\mathbf{u}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime}}^{l,\alpha^{\prime}}t(\mathbf{G}_{a}^{l}+\eta\mathbf{K}_{D}+\mathbf{k},0)\right),$
(49)

where we have used

$\displaystyle 3\mathbf{K}_{D}=\mathbf{G}_{1}+\mathbf{G}_{2}$ $\displaystyle
e^{i\tau_{A}\cdot(\mathbf{G}_{1}+\mathbf{G}_{2})}=e^{i\tau_{B}\cdot(\mathbf{G}_{1}+\mathbf{G}_{2})}=1.$
(50)

We next parameterize the above expression and define

$\displaystyle\gamma_{0}=\sum_{\mathbf{G}_{a}}t(\mathbf{G}_{a}+\eta\mathbf{K}_{D},0),$
$\displaystyle\gamma_{1}=\sum_{\mathbf{G}_{a}}(\widehat{\mathbf{G}_{a}+\eta\mathbf{K}_{D}})_{x}(\mathbf{G}_{a}+\eta\mathbf{K}_{D})_{x}\frac{\partial{t(\mathbf{q},0)}}{\partial
q}|_{q=|\mathbf{K}_{D}+\eta\mathbf{G}_{a}|},$ (51)
$\gamma_{2}=\sum_{\mathbf{G}_{a}}e^{i(\tau_{A}-\tau_{B})\cdot\mathbf{G}_{a}}(\widehat{\mathbf{G}_{a}+\mathbf{K}_{D}})_{x}(\mathbf{G}_{a}+\mathbf{K}_{D})_{x}\frac{\partial{t(\mathbf{q},0)}}{\partial
q}|_{q=|\mathbf{K}_{D}+\mathbf{G}_{a}|}$ (52)

and

$\gamma_{3}=2i\sum_{\mathbf{G}_{a}}e^{i(\tau_{A}-\tau_{B})\cdot\mathbf{G}_{a}}(\mathbf{G}_{a}+\mathbf{K}_{D})_{y}t(\mathbf{G}_{a}+\mathbf{K}_{D},0).$
(53)

We also define the acoustic and optical basis for the displacement field as

$\mathbf{u}_{\mathbf{p}}^{l,\alpha}=\mathbf{u}_{\mathbf{p}}^{l,ac}+(-1)^{\alpha}\mathbf{u}_{\mathbf{p}}^{l,op}$
(54)

where we define $(-1)^{\alpha=A(B)}=-1(+1)$.

With the above definitions, we can simplify the inter-valley e-ph Hamiltonian
as

$\displaystyle H_{\parallel
l}^{\eta,-\eta}(\mathbf{k}\mathbf{k}^{\prime})=-\frac{i}{\sqrt{N_{G}}}(\gamma_{0}+\gamma_{1})(\mathbf{k}-\mathbf{k}^{\prime})\cdot(\mathbf{u}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime}}^{l,ac}\sigma_{0}-\mathbf{u}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime}}^{l,op}\sigma_{z})$
$\displaystyle+\frac{1}{\sqrt{N_{G}}}\gamma_{3}\sigma_{x}(-u^{l,op}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},y}-i\eta
u^{l,ac}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},x})$
$\displaystyle+\frac{i}{\sqrt{N_{G}}}\gamma_{2}[(-(k_{x}-k^{\prime}_{x})\sigma_{x}+\eta(k_{y}+k^{\prime}_{y})\sigma_{y})u^{l,ac}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},x}+((k_{y}-k^{\prime}_{y})\sigma_{x}+\eta(k_{x}+k^{\prime}_{x})\sigma_{y})u^{l,ac}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},y}]$
$\displaystyle+\frac{1}{\sqrt{N_{G}}}\gamma_{2}[((k_{x}+k^{\prime}_{x})\sigma_{y}+\eta(k_{y}-k_{y}^{\prime})\sigma_{x})u^{l,op}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},x}+(-(k_{y}+k^{\prime}_{y})\sigma_{y}+\eta(k_{x}-k_{x}^{\prime})\sigma_{x})u^{l,op}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},y}],$
(55)

where the Pauli matrix $\sigma$ is for the sublattice basis
($(\sigma_{i})_{\alpha,\alpha^{\prime}}$ with $\alpha,\alpha^{\prime}=A,B$).

In addition, we also define

$v_{F}=\left|\sum_{\mathbf{G}_{a}}e^{i\mathbf{G}_{a}\cdot(\tau_{A}-\tau_{B})}\widehat{(\mathbf{K}_{D}+\mathbf{G}_{a})}_{x}\frac{\partial{t(\mathbf{q},0)}}{\partial
q}|_{q=|\mathbf{K}_{D}+\mathbf{G}_{a}|}\right|$ (56)

for the Fermi velocity of graphene.

All these parameters can be numerically evaluated and their values are listed
in Table (1) for both the $3$ $\mathbf{G}_{a}$ and for 2000 (infinite)
$\mathbf{G}_{a}$’s. One can see a large difference for the obtained values for
all the parameters when keeping only $3$ $\mathbf{G}_{a}$, as compared to
keeping 2000 (infinite) $\mathbf{G}_{a}$’s. For example, when keeping of order
2000 $G_{a}$’s, we find $\gamma_{1}=63.506$ and $\gamma_{0}=-61.4306$, which
almost cancel each other, and hence $\gamma=\gamma_{0}+\gamma_{1}=2.07565$.
For the first $3$ $\mathbf{G}_{a}$’s approximation, $\gamma_{0}=-22.7898$ and
$\gamma_{1}=12.7393$, and hence $\gamma=-10.0505$, which is larger and of
opposite sign as compared to keeping all $\mathbf{G}_{a}$. Here
$\gamma_{0},\gamma_{1}$ and $\gamma_{2}$ is in the unit of eV while
$\gamma_{3}$ is in the unit of $eV/\AA$.

| Infinite | 3 Shells | 3/Infinite
---|---|---|---
$v_{f}$ | $5.229$ | $7.3206$ | $1.4$
$\gamma_{0}$ | $-61.4306$ | $-22.7898$ | $0.371$
$\gamma_{1}$ | $63.5062$ | $12.7393$ | $0.2006$
$\gamma=\gamma_{0}+\gamma_{1}$ | $2.07565$ | $-10.0505$ | $-4.8421$
$\gamma_{2}$ | $6.7281$ | $6.3917$ | $0.95$
$\gamma_{3}$ | $-17.088$ | $-37.5936$ | $2.2$
Table 1: Table for all the parameters in the Hamiltonian (II.4). Here
$\gamma_{0},\gamma_{1}$ and $\gamma_{2}$ is in the unit of eV while
$\gamma_{3}$ is in the unit of $eV/\AA$.

In the above Hamiltonian (II.4), the term in the second line ($\gamma_{3}$
term) is zeroth order in $\mathbf{k}$ while all the other terms linearly
depend on $\mathbf{k}$. Thus, we below only keep the zeroth-order term for
e-ph interaction for K-phonons.

####  II.4.1 Physical Meaning of the phonon Modes at $\mathbf{K}_{D}$

We now focus on the leading term of the inter-valley electron-phonon
Hamiltonian

$\displaystyle
H=\frac{\gamma_{3}}{\sqrt{N_{G}}}\sum_{\mathbf{k},\mathbf{k}^{\prime},\eta,\alpha\beta}(-u^{l,op}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},y}-i\eta
u^{l,ac}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},x})c^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}(\sigma_{x})_{\alpha\beta}c_{\mathbf{k}^{\prime}-\eta\mathbf{K}_{D},\beta}$
(57)

where the Hamiltonian is hermitian once the summation over $\eta$ is
performed. We take the Fourier transform form of the phonon field

$\displaystyle\mathbf{u}_{\alpha}(\mathbf{k}-\mathbf{k}^{\prime}-\eta\mathbf{K}_{D})=\frac{1}{\sqrt{N_{G}}}\sum_{\mathbf{R}}e^{i\eta\mathbf{K}_{D}\cdot(\mathbf{R}+\tau_{\alpha})}e^{-i(\mathbf{k}-\mathbf{k}^{\prime})\cdot(\mathbf{R}+\tau_{\alpha})}\mathbf{u}_{\alpha}(\mathbf{R}+\tau_{\alpha}).$
(58)

We can see that the term $\eta\mathbf{K}_{D}\cdot\mathbf{R}$ is not periodic
in $\mathbf{R}$, and this problem can be fixed by enlarging the unit cell. By
choosing $\mathbf{R}=n\mathbf{a}_{1}+m\mathbf{a}_{2}$, we have

$\mathbf{K}_{D}\cdot\mathbf{R}=-\frac{2\pi}{3}(m+n).$ (59)

By choosing

$m+n=3n_{1}+s\;\;\;$ (60)

with $s=0,1,2$ and

$\sum_{\mathbf{R}}=\sum_{m,n}=\sum_{n_{1},m,s},$ (61)

We have

$\mathbf{R}^{\prime}=3n_{1}\mathbf{a}_{1}+m(\mathbf{a}_{2}-\mathbf{a}_{1}),\;\;\mathbf{R}=\mathbf{R}^{\prime}+s\mathbf{a}_{1},$
(62)

so

$\displaystyle\mathbf{u}_{\alpha}(\mathbf{k}-\mathbf{k}^{\prime}-\eta\mathbf{K}_{D})=\frac{1}{\sqrt{N_{G}}}\sum_{s,\mathbf{R}^{\prime}}e^{i\eta\mathbf{K}_{D}\cdot(s\mathbf{a}_{1}+\tau_{\alpha})}e^{-i(\mathbf{k}-\mathbf{k}^{\prime})\cdot(\mathbf{R}^{\prime}+s\mathbf{a}_{1}+\tau_{\alpha})}\mathbf{u}(\mathbf{R}^{\prime}+s\mathbf{a}_{1}+\tau_{\alpha})$
$\displaystyle=\frac{1}{\sqrt{3}}\sum_{s}e^{i\frac{2\pi}{3}\eta(1-s)}\mathbf{u}_{\alpha
s}(\mathbf{k}-\mathbf{k}^{\prime})$ (63)

where we have defined the $\mathbf{u}$ modes in the larger unit cell as

$\mathbf{u}_{\alpha
s}(\mathbf{q})=\sqrt{\frac{3}{N_{G}}}\sum_{\mathbf{R}^{\prime}}e^{-i\mathbf{q}\cdot(\mathbf{R}^{\prime}+s\mathbf{a}_{1}+\tau_{\alpha})}\mathbf{u}(\mathbf{R}^{\prime}+s\mathbf{a}_{1}+\tau_{\alpha}).$
(64)

For our choices of vectors, we note that

$\tau_{A}\cdot\mathbf{K}_{D}=\tau_{B}\cdot\mathbf{K}_{D}=\frac{2\pi}{3}$ (65)

The fact that they have the same value allows us to separate them once again
into optical and acoustic combinations

$u_{\alpha s}=u_{ac,s}+(-1)^{\alpha}u_{op,s}$ (66)

and to have

$\displaystyle\mathbf{u}_{ac/op}(\mathbf{k}-\mathbf{k}^{\prime}-\eta\mathbf{K}_{D})=\frac{1}{\sqrt{3}}\sum_{s}e^{i\frac{2\pi}{3}\eta(1-s)}\mathbf{u}_{ac/op,s}(\mathbf{k}-\mathbf{k}^{\prime})$
(67)

And the Hamiltonian becomes

$\displaystyle
H=\frac{\gamma_{3}}{\sqrt{3N_{G}}}\sum_{\mathbf{k},\mathbf{k}^{\prime},\eta,\alpha\beta}\sum_{s=1,2,3}e^{i\frac{2\pi}{3}\eta(1-s)}(-u_{op,s,y}(\mathbf{k}-\mathbf{k}^{\prime})-i\eta
u_{ac,s,x}(\mathbf{k}-\mathbf{k}^{\prime}))c^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}(\sigma_{x})_{\alpha\beta}c_{\mathbf{k}^{\prime}-\eta\mathbf{K}_{D},\beta}$
(68)

We may separate the valley index from electron momentum and introduce the
field operator

$\hat{\Psi}_{\mathbf{k},\alpha}=(c_{\mathbf{k}+\mathbf{K}_{D},\alpha},c_{\mathbf{k}-\mathbf{K}_{D},\alpha})^{T}$
(69)

and the Pauli matrix $\tau$ for the valley index. The Hamiltonian is then
written as

$\displaystyle
H=\frac{\gamma_{3}}{2\sqrt{3N_{G}}}\sum_{\mathbf{k},\mathbf{k}^{\prime}}\hat{\Psi}^{\dagger}_{\mathbf{k}}\sigma_{x}\left(\begin{array}[]{cc}0&\sum_{\alpha,s}(\mathbf{d}^{r}_{\alpha
s}\cdot\mathbf{u}_{\alpha,s}-i\mathbf{d}^{i}_{\alpha
s}\cdot\mathbf{u}_{\alpha,s})\\\ \sum_{\alpha,s}(\mathbf{d}^{r}_{\alpha
s}\cdot\mathbf{u}_{\alpha s}+i\mathbf{d}^{i}_{\alpha
s}\cdot\mathbf{u}_{\alpha,s})&0\end{array}\right)\hat{\Psi}_{\mathbf{k}^{\prime}}$
(72)
$\displaystyle=\frac{\gamma_{3}}{2\sqrt{3N_{G}}}\sum_{\mathbf{k},\mathbf{k}^{\prime},\alpha,s}\hat{\Psi}^{\dagger}_{\mathbf{k}}\sigma_{x}\left((\mathbf{d}^{r}_{\alpha
s}\cdot\mathbf{u}_{\alpha,s})\tau_{x}+(\mathbf{d}^{i}_{\alpha
s}\cdot\mathbf{u}_{\alpha,s})\tau_{y}\right)\hat{\Psi}_{\mathbf{k}^{\prime}},$
(73)

where we have used

$\displaystyle\sum_{s}e^{i\eta\theta_{s}}(-u_{op,s,y}-i\eta
u_{ac,s,x}))=\frac{1}{2}\sum_{s}e^{i\eta\theta_{s}}(-(u_{B,s,y}-u_{A,s,y})-i\eta(u_{A,s,x}+u_{B,s,x}))$
$\displaystyle=\frac{1}{2}\sum_{\alpha,s}(\mathbf{d}^{r}_{\alpha
s}\cdot\mathbf{u}_{\alpha s}+i\eta\mathbf{d}^{i}_{\alpha
s}\cdot\mathbf{u}_{\alpha,s})$ (74)

with

$\displaystyle\mathbf{d}^{r}_{A,s}=\cos{\theta_{s}}\hat{e}_{y}+\sin{\theta_{s}}\hat{e}_{x},\;\;\mathbf{d}^{r}_{B,s}=-\cos{\theta_{s}}\hat{e}_{y}+\sin{\theta_{s}}\hat{e}_{x},$
$\displaystyle\mathbf{d}^{i}_{A,s}=\sin{\theta_{s}}\hat{e}_{y}-\cos{\theta_{s}}\hat{e}_{x},\;\;\mathbf{d}^{i}_{B,s}=-\sin{\theta_{s}}\hat{e}_{y}-\cos{\theta_{s}}\hat{e}_{x},$
(75)

and $\theta_{s}=\frac{2\pi}{3}(1-s)$. Here $\mathbf{u}_{\alpha,s}$ is a
function of $\mathbf{k}-\mathbf{k}^{\prime}$. The above e-ph Hamiltonian takes
the same form as in Ref. [25].

Now let’s focus on the vibration mode at $\mathbf{K}_{D}$
($\mathbf{k}\sim\mathbf{k}^{\prime}\sim 0$), so
$\mathbf{u}_{\alpha,s}(0)=\sqrt{3/N_{G}}\sum_{\mathbf{R}^{\prime}}\mathbf{u}(\mathbf{R}^{\prime}+s\mathbf{a}_{1}+\tau_{\alpha})$.
We can view $\mathbf{d}^{r(i)}_{\alpha s}$ as the vector, along which the
$(\alpha,s)$-atom vibrates. Explicitly, the $\mathbf{d}^{r}$-vectors are given
by
$\mathbf{d}^{r}_{A,0}=-\frac{1}{2}\hat{e}_{y}+\frac{\sqrt{3}}{2}\hat{e}_{x}$,
$\mathbf{d}^{r}_{B,0}=\frac{1}{2}\hat{e}_{y}+\frac{\sqrt{3}}{2}\hat{e}_{x}$,
$\mathbf{d}^{r}_{A,1}=\hat{e}_{y}$, $\mathbf{d}^{r}_{B,1}=-\hat{e}_{y}$,
$\mathbf{d}^{r}_{A,2}=-\frac{1}{2}\hat{e}_{y}-\frac{\sqrt{3}}{2}\hat{e}_{x}$
and
$\mathbf{d}^{r}_{B,2}=\frac{1}{2}\hat{e}_{y}-\frac{\sqrt{3}}{2}\hat{e}_{x}$,
which can be depicted by the green arrows in Fig. 4b and correspond to the
$A_{1}$ phonon mode. The vibration configuration of each atoms for the
$\mathbf{d}^{i}$ phonon mode is shown in Fig. 4c and corresponds to the
$A_{2}$ mode (Here $A_{1}$ and $A_{2}$ corresponds to the irreducible
representations of the $C_{3v}$ group).

Figure 4:  (a) Enlarged unit cell for phonon modes at $\mathbf{K}_{D}$. (b)
and (c) are for the $A_{1}$ and $A_{2}$ phonon modes at $\mathbf{K}_{D}$.

###  II.5 Symmetry properties of the Electron-Phonon Interaction

In this section, we will present a construction of the electron-phonon
Hamiltonian from the symmetry principle. We now write all the symmetry
operators on the electron and phonon modes (we consider the phonon as a field
and expand it around $\eta\mathbf{K}_{D}$):

$\displaystyle\hat{C}_{3}\hat{c}^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}\hat{C}_{3}^{-1}=e^{i\eta\tau_{\alpha}\cdot\mathbf{G}_{1}}\hat{c}^{\dagger}_{\hat{C}_{3}\mathbf{k}+\eta\mathbf{K}_{D},\alpha}$
$\displaystyle\hat{M}_{y}\hat{c}^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}\hat{M}_{y}^{-1}=\sum_{\beta}(\sigma_{x})_{\beta\alpha}\hat{c}^{\dagger}_{\hat{M}_{y}\mathbf{k}+\eta\mathbf{K}_{D},\beta}$
$\displaystyle\hat{M}_{x}\hat{c}^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}\hat{M}_{x}^{-1}=\hat{c}^{\dagger}_{\hat{M}_{x}\mathbf{k}-\eta\mathbf{K}_{D},\alpha}$
$\displaystyle\hat{T}\hat{c}^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}\hat{T}^{-1}=\hat{c}^{\dagger}_{-\mathbf{k}-\eta\mathbf{K}_{D},\alpha}$
$\displaystyle\hat{c}^{\dagger}_{\mathbf{k}+\mathbf{G},\alpha}=e^{i\tau_{\alpha}\cdot\mathbf{G}}\hat{c}^{\dagger}_{\mathbf{k},\alpha}$
$\displaystyle u_{\alpha
i}(\mathbf{p}+\mathbf{G})=e^{-i\tau_{\alpha}\cdot\mathbf{G}}u_{\alpha
i}(\mathbf{p})$ $\displaystyle\hat{C}_{3}u_{\alpha
i}(\mathbf{k}+\eta\mathbf{K}_{D})\hat{C}_{3}^{-1}=e^{-i\eta\tau_{\alpha}\cdot\mathbf{G}_{1}}\sum_{i^{\prime}}u_{\alpha
i^{\prime}}(\hat{C}_{3}\mathbf{k}+\eta\mathbf{K}_{D})C_{3i^{\prime}i}$
$\displaystyle\hat{M}_{y}u_{\alpha
i}(\mathbf{k}+\eta\mathbf{K}_{D})\hat{M}_{y}^{-1}=\sum_{i^{\prime}\beta}(\sigma_{x})_{\beta\alpha}u_{\beta
i^{\prime}}(\hat{M}_{y}\mathbf{k}+\eta\mathbf{K}_{D})M_{yi^{\prime}i}$
$\displaystyle\hat{M}_{x}u_{\alpha
i}(\mathbf{k}+\eta\mathbf{K}_{D})\hat{M}_{x}^{-1}=\sum_{i^{\prime}}u_{\alpha
i^{\prime}}(\hat{M}_{x}\mathbf{k}-\eta\mathbf{K}_{D})M_{xi^{\prime}i}$
$\displaystyle\hat{T}u_{\alpha
i}(\mathbf{k}+\eta\mathbf{K}_{D})\hat{T}^{-1}=u_{\alpha
i}(-\mathbf{k}-\eta\mathbf{K}_{D})$ (76)

In the above, the matrix representations are defined as

$\displaystyle C_{3}=\left(\begin{array}[]{cc}-1/2&-\sqrt{3}/2\\\
\sqrt{3}/2&-1/2\end{array}\right)$ (79) $\displaystyle
M_{x}=\left(\begin{array}[]{cc}-1&0\\\ 0&1\end{array}\right)$ (82)
$\displaystyle M_{y}=\left(\begin{array}[]{cc}1&0\\\ 0&-1\end{array}\right).$
(85)

We now use the above symmetry transformation to construct the inter-valley
electron-phonon Hamiltonian.

The generic form of the intra-valley electron-phonon Hamiltonian reads

$\displaystyle
H_{2}=\frac{1}{\sqrt{N_{G}}}\sum_{\mathbf{k},\mathbf{k}^{\prime},\eta;\alpha\beta\theta;j}F_{\alpha\beta\theta;j}^{\eta}(\mathbf{k},\mathbf{k}^{\prime})c^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}u_{\theta
j}(\mathbf{k}-\mathbf{k}^{\prime}-\eta\mathbf{K}_{D})c_{\mathbf{k}^{\prime}-\eta\mathbf{K}_{D},\beta}$
(86)

where $F_{\alpha\beta\theta;j}^{\eta}(\mathbf{k},\mathbf{k}^{\prime})$ is the
generic inter-valley electron-phonon interaction parameter. We consider the
case of $F_{\alpha\beta\theta;j}^{\eta}(\mathbf{k},\mathbf{k}^{\prime})$
independent of $\mathbf{k},\mathbf{k}^{\prime}$ in order to obtain the leading
term. However, we first obtain the symmetry conditions considering the entire
functional dependence. Under the generic $\hat{g}$-symmetry transformations,
we have

$\displaystyle\hat{g}c^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}\hat{g}^{-1}=\sum_{\eta^{\prime},\alpha^{\prime}}c^{\dagger}_{\hat{g}\mathbf{k}+\eta^{\prime}\mathbf{K}_{D},\alpha^{\prime}}D_{c}(\hat{g})_{\eta^{\prime}\alpha^{\prime};\eta\alpha}$
$\displaystyle\hat{g}u_{\alpha
i}(\mathbf{p}+\eta\mathbf{K}_{D})\hat{g}^{-1}=\sum_{\alpha^{\prime}i^{\prime}\eta^{\prime}}u_{\alpha^{\prime}i^{\prime}}(\hat{g}\mathbf{p}+\eta^{\prime}\mathbf{K}_{D})S_{i^{\prime}i}^{\alpha^{\prime}\eta^{\prime};\alpha\eta}(\hat{g})$
(87)

where $S_{i^{\prime}\eta^{\prime};i\eta}^{\alpha^{\prime}\alpha}$ is a
combination of the space rotation and sublattice and valley rotation. The
equation $\hat{g}H\hat{g}^{-1}=H$ gives

$\displaystyle\sum_{\eta\alpha\beta\theta
i}D_{c}(\hat{g})_{\eta^{\prime}\alpha^{\prime};\eta\alpha}S_{i^{\prime},i}^{\theta^{\prime}-\eta^{\prime};\theta-\eta}D^{\star}_{c}(\hat{g})_{-\eta^{\prime}\beta^{\prime};-\eta\beta}F_{\alpha\beta\theta;i}^{\eta}(\mathbf{k},\mathbf{k}^{\prime})=F_{\alpha^{\prime}\beta^{\prime}\theta^{\prime};i^{\prime}}^{\eta^{\prime}}(\hat{g}\mathbf{k},\hat{g}\mathbf{k}^{\prime})$
(88)

We now write the explicit form of the matrices and the constraints

$\displaystyle\hat{M}_{y}:D_{c\eta^{\prime}\alpha^{\prime};\eta\alpha}=\delta_{\eta^{\prime}\eta}\sigma^{x}_{\alpha^{\prime}\alpha},\;\;S_{i^{\prime}i}^{\theta^{\prime}\eta^{\prime};\theta\eta}=\sigma^{x}_{\theta^{\prime}\theta}(-1)^{i+1}\delta_{i^{\prime}i}\delta_{\eta,\eta^{\prime}}$
$\displaystyle\implies
F_{\bar{\alpha}\bar{\beta}\bar{\theta};i}^{\eta}(\hat{M}_{y}\mathbf{k},\hat{M}_{y}\mathbf{k}^{\prime})=(-1)^{i+1}F_{\alpha\beta\theta;i}^{\eta}(\mathbf{k},\mathbf{k}^{\prime})$
$\displaystyle\hat{M}_{x}:D_{c\eta^{\prime}\alpha^{\prime};\eta\alpha}=\delta_{\eta^{\prime}\bar{\eta}}\delta_{\alpha^{\prime}\alpha},\;\;S_{i^{\prime}i}^{\theta^{\prime}\eta^{\prime};\theta\eta}=\delta_{\theta^{\prime}\theta}(-1)^{i}\delta_{i^{\prime}i}\delta_{\eta^{\prime}\bar{\eta}}$
$\displaystyle\implies
F_{\alpha\beta\theta;i}^{\bar{\eta}}(\hat{M}_{x}\mathbf{k},\hat{M}_{x}\mathbf{k}^{\prime})=(-1)^{i}F_{\alpha\beta\theta;i}^{\eta}(\mathbf{k},\mathbf{k}^{\prime})$
$\displaystyle\hat{C}_{3}:D_{c\eta^{\prime}\alpha^{\prime};\eta\alpha}=\delta_{\eta^{\prime}\eta}\delta_{\alpha^{\prime}\alpha}e^{i\eta\tau_{\alpha}\cdot\mathbf{G}_{1}};\;\;S_{i^{\prime}i}^{\theta^{\prime}\eta^{\prime};\theta\eta}=e^{-i\eta\tau_{\theta}\cdot\mathbf{G}_{1}}\delta_{\eta^{\prime}\eta}\delta_{\theta^{\prime}\theta}C_{3i^{\prime}i}$
$\displaystyle\implies
F_{\alpha\beta\theta;i^{\prime}}^{\eta}(\hat{C}_{3}\mathbf{k},\hat{C}_{3}\mathbf{k}^{\prime})=e^{i\eta(\tau_{\alpha}+\tau_{\beta}+\tau_{\theta})\cdot\mathbf{G}_{1}}\times\sum_{i}F_{\alpha\beta\theta;i}^{\eta}(\mathbf{k},\mathbf{k}^{\prime})C_{3i^{\prime}i}$
$\displaystyle\hat{T}:D_{c\eta^{\prime}\alpha^{\prime};\eta\alpha}=\delta_{\eta^{\prime}\bar{\eta}}\delta_{\alpha^{\prime}\alpha},\;\;S_{i^{\prime}i}^{\theta^{\prime}\eta^{\prime};\theta\eta}=\delta_{\theta^{\prime}\theta}\delta_{i^{\prime}i}\delta_{\eta^{\prime},\bar{\eta}}$
$\displaystyle\implies
F_{\alpha\beta\theta;i}^{\bar{\eta}}(-\mathbf{k},-\mathbf{k}^{\prime})=F_{\alpha\beta\theta;i}^{\eta\star}(\mathbf{k},\mathbf{k}^{\prime})$
(89)

The Hermiticity of the $H_{2}$ Hamiltonian requires

$\text{Hermiticity:}\;\;\;F_{\alpha\beta\theta;i}^{\eta}(\mathbf{k},\mathbf{k}^{\prime})=F_{\beta\alpha\theta;i}^{\bar{\eta}\star}(\mathbf{k}^{\prime},\mathbf{k})$
(90)

We are aided by the following identity

$e^{i\eta(\tau_{\alpha}+\tau_{\beta}+\tau_{\theta})\cdot\mathbf{G}_{1}}=e^{i(-1)^{\eta}((-1)^{\alpha}+(-1)^{\beta}+(-1)^{\theta})\frac{2\pi}{3}}$
(91)

This is a simple system of linear equations to solve, giving rise to the _two_
possible solutions. The first solution is given by

$\displaystyle F^{\eta}_{\alpha\beta\theta
x}=-\frac{\gamma_{3}}{2}\eta\sigma^{x}_{\alpha\beta},\;\;F^{\eta}_{\alpha\beta\theta
y}=-\frac{\gamma_{3}}{2}\sigma^{x}_{\alpha\beta}(-1)^{\theta}$ (92)

which corresponds to the inter-valley phonon contribution that we have
obtained, $-\gamma_{3}(i\eta u^{ac}_{x}+u^{op}_{y})\sigma_{x}$.

The second solution is the missing term whose only nonzero components are

$F_{\alpha\alpha\bar{\alpha};x}^{\eta}=i\eta;\;\;F_{\alpha\alpha\bar{\alpha};y}^{\eta}=(-1)^{\alpha}$
(93)

which would correspond to a new term in the Hamiltonian

$\displaystyle\sum_{\mathbf{k}\mathbf{k}^{\prime}\eta}c^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}c_{\mathbf{k}^{\prime}-\eta\mathbf{K}_{D}\alpha}\times(i\eta
u_{\bar{\alpha}x}(\mathbf{k}-\mathbf{k}^{\prime}-\eta\mathbf{K}_{D})+(-1)^{\alpha}u_{\bar{\alpha}y}(\mathbf{k}-\mathbf{k}^{\prime}-\eta\mathbf{K}_{D}))$

from which we can see two crucial properties of the extra terms

  * •

First, it is diagonal in the electron sublattice index, hence it can only come
from second nearest neighbor, or larger.

  * •

Crucially, the sublattice index of the electron $\alpha$ and that of the
phonon $\bar{\alpha}$ are opposite. It means that the phonon on one sublattice
couples to the electron annihilation and creation operators on the other
sublattice. Hence it is now clear why this term does not appear in our
Hamiltonian or _in any other two-center approximation_. The two center
approximation reads

$t(R_{\alpha}+u_{\alpha}-R_{\beta}-u_{\beta})c_{R_{\alpha}}^{\dagger}c_{R_{\beta}}$
(94)

and hence the phonon sublattice index has to be either that of $c^{\dagger}$
or that of $c$. It cannot be neither that of $c^{\dagger}$ nor that of $c$. We
hence see that the absence of this term is due to the two-center approximation
of our Hamiltonian (and of most tight-binding Hamiltonians).

##  III Phonons in a Single Layer Graphene

In this section, we will discuss the phonon dynamical matrix and phonon eigen-
modes. In real-space, in any dimensions, the dynamica matrix $D$ reads

$(D_{\mathbf{R}\alpha,\mathbf{R}^{\prime},\alpha^{\prime}})_{ij}$ (95)

where $\mathbf{R},\mathbf{R}^{\prime}$ are the unit cells, the
$\alpha,\alpha^{\prime}$ are the atoms, and the $i,j\in x,y,z$ are the
components in 3D. The presence of stability (the free energy is invariant to
adding any constant vector $\vec{a}$ to the displacements:
$\mathbf{u}_{r\alpha}\rightarrow\mathbf{u}_{r\alpha}+\mathbf{a}$ ) gives

$\sum_{\mathbf{R}^{\prime},\alpha^{\prime}}(D_{\mathbf{R}\alpha,\mathbf{R}^{\prime},\alpha^{\prime}})_{ij}=0$
(96)

Which gives, in a spring model,

$\sum_{\beta}D_{\alpha\beta}(\mathbf{q}=0)=0,\forall\alpha;\;\sum_{\alpha}D_{\alpha\beta}(\mathbf{q}=0)=0,\forall\beta$
(97)

Which guarantees the presence of an acoustic phonon mode.

###  III.1 Dynamical Matrix

The dynamical matrix can be obtained from the potential $V$ of the system

$(D_{\mathbf{R}\alpha,\mathbf{R}^{\prime}\alpha^{\prime}})_{ii^{\prime}}=\frac{\partial^{2}V}{\partial
u_{\mathbf{R}\alpha i}\partial
u_{\mathbf{R}^{\prime}\alpha^{\prime}i^{\prime}}}$ (98)

As such, if there exists a symmetry $\hat{g}$ of the system with
representation matrix $G_{ij}$, we have the relation

$D_{\hat{g}\mathbf{R},\hat{g}\alpha;\hat{g}\mathbf{R}^{\prime},\hat{g}\alpha^{\prime}}=GD_{\mathbf{R}\alpha;\mathbf{R}^{\prime}\alpha^{\prime}}G^{-1}$
(99)

Next let’s consider the constraint on the dynamical matrix from the symmetry
operators, $\hat{C}_{3},\hat{M}_{x}$ $\hat{M}_{y}$ and $\hat{T}$, with their
representation matrices given in the above. Due to translation symmetry,
$(D_{\mathbf{R}\alpha,\mathbf{R}^{\prime}\alpha^{\prime}})_{ii^{\prime}}=(D_{\alpha\alpha^{\prime}}(\mathbf{R}_{\alpha}-\mathbf{R}^{\prime}_{\alpha^{\prime}}))_{ii^{\prime}}$
only depends $\mathbf{R}_{\alpha}-\mathbf{R}^{\prime}_{\alpha^{\prime}}$. We
only consider the nearest neighbor (nn) and next nearest neighbor (nnn)
interaction. Three nn vectors are labelled as
$\delta^{nn}_{AB,\mu=1,2,3}=(a_{0}/2,a_{0}/2\sqrt{3}),(-a_{0}/2,a_{0}/2\sqrt{3}),(0,-a_{0}/\sqrt{3})$
that satisfies $\hat{C}_{3}\delta^{nn}_{i}=\delta^{nn}_{i+1}$, and six nnn
vectors are
$\delta^{nnn}_{\alpha\alpha,\mu=1,2,3}=-\mathbf{a}_{1}-\mathbf{a}_{2},\mathbf{a}_{1},\mathbf{a}_{2}$
and $-\delta^{nnn}_{\alpha\alpha,\mu=1,2,3}$ with
$\hat{C}_{3}\delta^{nnn}_{i}=\delta^{nnn}_{i+1}$. In the above expressions, we
use the convention $\delta^{nn}_{4}=\delta^{nn}_{1}$ and
$\delta^{nnn}_{4}=\delta^{nnn}_{1}$.

For the nn interaction that couples the A and B atoms,
$\delta^{nn}_{AB,3}=(0,-a_{0}/\sqrt{3})$ remains the same under $\hat{M}_{x}$
and thus

$M_{x}D_{AB}(\delta^{nn}_{3})M_{x}^{-1}=D_{AB}(\delta^{nn}_{3}),$ (100)

which requires the diagonal form of $D_{AB}(\delta^{nn}_{3})$, given by

$D_{AB}(\delta^{nn}_{3})=-\left(\begin{array}[]{cc}K_{0t}&0\\\
0&K_{0r}\end{array}\right).$ (101)

Other nn interaction terms can be related to $D_{AB}(\delta^{nn}_{3})$ through
$C_{3}$ rotation,

$D_{AB}(\delta^{nn}_{i+1})=C_{3}D_{AB}(\delta^{nn}_{i})C_{3}^{-1}.$ (102)

$\hat{M}_{y}$ reverses the sign of $\delta^{nn}_{AB,\mu}$ and interchanges A
and B sublattice, and thus

$D_{BA}(-\delta^{nn}_{i})=M_{y}D_{AB}(\delta^{nn}_{i})M_{y}^{-1}.$ (103)

For the nnn interaction, we first consider the matrix
$D_{AA}(\delta^{nnn}_{1})$. Mirror $\hat{M}_{x}$ requires

$D_{AA}(-\delta^{nnn}_{1})=M_{x}D_{AA}(\delta^{nnn}_{1})M_{x}^{-1}$ (104)

. Furthermore, the symmetric condition for the dynamical matrix gives rise to

$D_{AA}(\delta^{nnn}_{1})=(D_{AA}(-\delta^{nnn}_{1}))^{T}.$ (105)

Thus, we have the condition

$(D_{AA}(\delta^{nnn}_{1}))^{T}=M_{x}D_{AA}(\delta^{nnn}_{1})M_{x}^{-1},$
(106)

which leads to the form

$D_{AA}(\delta^{nnn}_{1})=-\left(\begin{array}[]{cc}K_{r1}&K_{2}\\\
-K_{2}&K_{t1}\end{array}\right).$ (107)

$D_{AB}(\delta^{nnn}_{2,3})$ can be related to $D_{AB}(\delta^{nnn}_{1})$ from
$C_{3}$ rotation,

$D_{AB}(\delta^{nnn}_{i+1})=C_{3}D_{AB}(\delta^{nnn}_{i})C_{3}^{-1},$ (108)

while $D_{AB}(-\delta^{nnn}_{i})$ can be related to $D_{AB}(\delta^{nnn}_{i})$
from mirror $\hat{M}_{x}$,

$D_{AB}(-\delta^{nnn}_{i})=M_{x}D_{AB}(\delta^{nnn}_{i})M_{x}^{-1}.$ (109)

Finally, time reversal $\hat{T}$ requires all the five parameters
$K_{t0},K_{r0},K_{t1},K_{r1},K_{2}$ to be real. With all the components in the
real space, we can transform the dynamical matrix to the momentum space, which
is given by

$D_{\alpha\alpha^{\prime}}(\mathbf{k})=\sum_{\mathbf{R}}e^{-i\mathbf{k}\cdot(\mathbf{R}+\tau_{\alpha}-\tau_{\alpha^{\prime}})}D_{\alpha\alpha^{\prime}}(\mathbf{R}+\tau_{\alpha}-\tau_{\alpha^{\prime}})$
(110)

The explicit form of the dynamical matrix under the basis
$(u_{Bx},u_{By},u_{Ax},u_{Ay})^{T}$ is given by

$\displaystyle D(\mathbf{k})=d_{ij}(k_{x},k_{y})\sigma_{i}\otimes\xi_{j}$
$\displaystyle
d_{00}(k_{x},k_{y})=\left(K_{\text{r1}}+K_{\text{t1}}\right)\left(-2\cos\left(\frac{k_{x}}{2}\right)\cos\left(\frac{\sqrt{3}k_{y}}{2}\right)-\cos\left(k_{x}\right)+3\right)+\frac{1}{2}\left(3K_{\text{r0}}+3K_{\text{t0}}\right)$
$\displaystyle
d_{01}(k_{x},k_{y})=\sqrt{3}\sin\left(\frac{k_{x}}{2}\right)\sin\left(\frac{\sqrt{3}k_{y}}{2}\right)\left(K_{\text{r1}}-K_{\text{t1}}\right)$
$\displaystyle
d_{03}(k_{x},k_{y})=-\left(K_{\text{r1}}-K_{\text{t1}}\right)\left(\cos\left(k_{x}\right)-\cos\left(\frac{k_{x}}{2}\right)\cos\left(\frac{\sqrt{3}k_{y}}{2}\right)\right)$
$\displaystyle
d_{11}(k_{x},k_{y})=\frac{1}{2}\sqrt{3}\left(K_{\text{r0}}-K_{\text{t0}}\right)\sin\left(\frac{k_{x}}{2}\right)\sin\left(\frac{k_{y}}{2\sqrt{3}}\right)$
$\displaystyle
d_{21}(k_{x},k_{y})=\frac{1}{2}\sqrt{3}\left(K_{\text{r0}}-K_{\text{t0}}\right)\sin\left(\frac{k_{x}}{2}\right)\cos\left(\frac{k_{y}}{2\sqrt{3}}\right)$
$\displaystyle
d_{10}(k_{x},k_{y})=-\frac{1}{2}\left(K_{\text{r0}}+K_{\text{t0}}\right)\left(2\cos\left(\frac{k_{x}}{2}\right)\cos\left(\frac{k_{y}}{2\sqrt{3}}\right)+\cos\left(\frac{k_{y}}{\sqrt{3}}\right)\right)$
$\displaystyle
d_{20}(k_{x},k_{y})=\frac{1}{2}\left(K_{\text{r0}}+K_{\text{t0}}\right)\left(2\cos\left(\frac{k_{x}}{2}\right)\sin\left(\frac{k_{y}}{2\sqrt{3}}\right)-\sin\left(\frac{k_{y}}{\sqrt{3}}\right)\right)$
$\displaystyle
d_{13}(k_{x},k_{y})=-\frac{1}{2}\left(K_{\text{r0}}-K_{\text{t0}}\right)\left(\cos\left(\frac{k_{x}}{2}\right)\cos\left(\frac{k_{y}}{2\sqrt{3}}\right)-\cos\left(\frac{k_{y}}{\sqrt{3}}\right)\right)$
$\displaystyle
d_{23}(k_{x},k_{y})=\frac{1}{2}\left(K_{\text{r0}}-K_{\text{t0}}\right)\sin\left(\frac{k_{y}}{2\sqrt{3}}\right)\left(\cos\left(\frac{k_{x}}{2}\right)+2\cos\left(\frac{k_{y}}{2\sqrt{3}}\right)\right)$
$\displaystyle
d_{32}(k_{x},k_{y})=4K_{2}\sin\left(\frac{k_{x}}{2}\right)\left(\cos\left(\frac{k_{x}}{2}\right)-\cos\left(\frac{\sqrt{3}k_{y}}{2}\right)\right)$
(111)

where $\sigma$ operates in the $A,B$ sub-lattices while $\xi$ operates in the
$x,y$ components. The eigen-equation of dynamical matrix is given by

$D(\mathbf{k})\epsilon_{\mathbf{k}}=M\omega^{2}\epsilon_{\mathbf{k}}.$ (112)

Below we focus on the dynamical matrix and the corresponding phonon modes at
$\mathbf{K}_{D}$.

###  III.2 Dynamical Matrix and Phonon Modes at $\mathbf{K}_{D}$

Close to the $K_{D}$ point, the dynamical matrix is not solvable, so we first
consider the eigen-states at the $\mathbf{K}_{D}$ point

$\displaystyle
D(\eta\mathbf{K}_{D})=\frac{3}{2}\left(K_{\text{r0}}+3K_{\text{r1}}+K_{\text{t0}}+3K_{\text{t1}}\right)\sigma_{0}\otimes\sigma_{0}+$
$\displaystyle+\frac{3}{4}\left(K_{\text{r0}}-K_{\text{t0}}\right)(\sigma_{1}\otimes\sigma_{3}+\eta\sigma_{2}\otimes\sigma_{1})-$
$\displaystyle-\eta 3\sqrt{3}K_{2}\sigma_{3}\otimes\sigma_{3}$ (113)

Let’s define

$\displaystyle\alpha_{1}=\frac{3}{2}\left(K_{\text{r0}}+3K_{\text{r1}}+K_{\text{t0}}+3K_{\text{t1}}\right)$
$\displaystyle\alpha_{2}=\frac{3}{2}\left(K_{\text{r0}}-K_{\text{t0}}\right)$
$\displaystyle\alpha_{3}=3\sqrt{3}K_{2}$ (114)

and the eigenvalues and eigenstates are

$\displaystyle
E_{E_{-}}=\alpha_{1}-\alpha_{3};\;E_{E_{+}}=\alpha_{1}-\alpha_{3}$ (115)
$\displaystyle E_{A_{1}}=\alpha_{1}+\alpha_{2}+\alpha_{3}$ $\displaystyle
E_{A_{2}}=\alpha_{1}-\alpha_{2}+\alpha_{3}$ $\displaystyle\epsilon_{\eta
K_{D}}^{E_{-}}=\frac{1}{\sqrt{2}}(0,0,\eta i,1)^{T},\;\;\epsilon_{\eta
K_{D}}^{E_{+}}=\frac{1}{\sqrt{2}}(-\eta i,1,0,0)^{T},$
$\displaystyle\epsilon_{\eta K_{D}}^{A_{1}}=\frac{1}{2}(-\eta i,-1,-\eta
i,1)^{T},\;\;\epsilon_{\eta K_{D}}^{A_{2}}=\frac{1}{2}(\eta i,1,-\eta
i,1)^{T}$

and the eigen-frequencies are related to eigenvalues by

$\omega_{r}=\sqrt{E_{r}/M},\;\;r=E_{\pm},A_{2},A_{1}.$ (116)

The projected Dynamical matrix into the basis of the doublet $E_{\pm}$ is

$\displaystyle\left\langle\epsilon_{\eta\mathbf{K}_{D}}^{E_{\mp}}\right|D(\eta\mathbf{K}_{D})\left|\epsilon_{\eta\mathbf{K}_{D}}^{E_{\mp}}\right\rangle=$
$\displaystyle=E_{E_{-}}\sigma_{0}+\frac{\alpha_{2}}{2\sqrt{3}}(k_{y}\sigma_{y}-\eta
k_{x}\sigma_{x})$ (117)

The symmetry properties of these states are:

$\displaystyle\hat{M}_{y}\epsilon_{\eta\mathbf{K}_{D}}^{A_{1}}=\epsilon_{\eta\mathbf{K}_{D}}^{A_{1}};\;\;\hat{M}_{y}\epsilon_{\eta\mathbf{K}_{D}}^{A_{2}}=-\epsilon_{\eta\mathbf{K}_{D}}^{A_{2}}$
$\displaystyle\hat{M}_{y}\epsilon_{\eta\mathbf{K}_{D}}^{E_{-}}=-\epsilon_{\eta\mathbf{K}_{D}}^{E_{+}},\;\;\;\hat{M}_{y}\epsilon_{\eta\mathbf{K}_{D}}^{E_{+}}=-\epsilon_{\eta\mathbf{K}_{D}}^{E_{-}}$
(118)

The energies of the phonon bands away from $\eta\mathbf{K}_{D}$ are

$\displaystyle\left\langle\epsilon_{\eta\mathbf{K}_{D}}^{A_{1,2}}\right|D(\eta\mathbf{K}_{D})\left|\epsilon_{\eta\mathbf{K}_{D}}^{A_{1,2}}\right\rangle=$
$\displaystyle=E_{A_{1,2}}+\frac{k^{2}}{8}(-6\sqrt{3}K_{2}-K_{r0}-3K_{r1}+K_{t0}-3K_{t1})$
(119)

Hence the effective masses of phonons away from $\eta\mathbf{K}_{D}$ are the
same for the $A_{1,2}$ modes.

##  IV Electron and Phonon-Band Projected Electron-Phonon Hamiltonian

Since the inter-valley phonons always have a finite frequency, we keep only
the zeroth-order terms of momenta and neglect the rest for the inter-valley
electron-phonon Hamiltonian

$\displaystyle H_{\parallel
l}^{\eta,-\eta}(\mathbf{k},\mathbf{k}^{\prime})=\frac{1}{\sqrt{N_{G}}}\gamma_{3}\sigma_{x}(-u^{l,op}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},y}-i\eta
u^{l,ac}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},x})$ (120)

For the electron Dirac Hamiltonian $v_{f}(\eta
k_{x}\sigma_{x}+k_{y}\sigma_{y})$, we have the energies $E_{\pm}=\pm v_{f}k$
and the corresponding eigen-states

$\displaystyle\psi^{\pm}_{\mathbf{k},\eta}=\frac{1}{\sqrt{2}}\left[\begin{array}[]{c}\pm\eta
e^{-i\eta\phi}\\\ 1\end{array}\right],\;\;\;e^{i\phi}=\frac{k_{x}+ik_{y}}{k}$
(123)

The electron operator can also be expanded on the eigen-states with

$c_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}=\sum_{m=1,2}\psi^{m}_{\mathbf{k},\eta\alpha}\gamma_{\mathbf{k},\eta}^{m}.$
(124)

For the phonon part, with the eigenstates of the dynamical matrix
$\epsilon^{n}_{\mathbf{q}\alpha}$, we have that the displacements
$b_{\mathbf{q}n}$ become decomposed in phonon band basis

$u_{\mathbf{q}t}=\sum_{n=1,2,3,4}\epsilon_{\mathbf{q}t}^{n}\frac{b_{\mathbf{q}n}+b_{-\mathbf{q}n}^{\dagger}}{\sqrt{2M\omega_{\mathbf{q}n}}}$
(125)

where $t=1,2,3,4$ corresponds to $u_{Bx},u_{By},u_{Ax},u_{Ay}$, respectively,
and $n$ labels the energy eigenstates (from lower energy to higher energy).

Below we will project the electron-phonon Hamiltonian into the eigen-state
basis of both electrons and phonons.

###  IV.1 Projected Inter-Valley Electron-Phonon Interaction

We next consider the inter-valley e-ph interaction and focus on the
$k$-independent terms.

$\displaystyle
H=\frac{1}{\sqrt{N_{G}}}\gamma_{3}\sum_{\mathbf{k},\mathbf{k}^{\prime},\eta,\alpha\beta}(-u^{l,op}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},y}-i\eta
u^{l,ac}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},x})c^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}(\sigma_{x})_{\alpha\beta}c_{\mathbf{k}^{\prime}-\eta\mathbf{K}_{D},\beta}$
$\displaystyle=-\frac{1}{2\sqrt{N_{G}}}\gamma_{3}\sum_{m,n=\pm}\sum_{r=1}^{4}\sum_{\mathbf{k},\mathbf{k}^{\prime},\eta,\alpha\beta}\frac{b_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},r}+b_{\eta\mathbf{K}_{D}-\mathbf{k}+\mathbf{k}^{\prime},r}^{\dagger}}{\sqrt{2M\omega_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},r}}}\gamma_{\mathbf{k},\eta}^{n\dagger}\gamma_{\mathbf{k}^{\prime},-\eta}^{m}$
(126)
$\displaystyle((\epsilon^{r}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},2}-\epsilon^{r}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},4})+i\eta(\epsilon^{r}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},1}+\epsilon^{r}_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},3}))\psi_{\mathbf{k},\eta,\alpha}^{n\star}(\sigma_{x})_{\alpha\beta}\psi_{\mathbf{k}^{\prime},-\eta,\beta}^{m}$

Using the phonon eigen-states obtained at $\mathbf{K}_{D}$, we can obtain the
approximate form of the couplings for $\mathbf{k}=\mathbf{k}^{\prime}$ as

$\displaystyle
r=E_{-}:\;\;((\epsilon^{r}_{-\eta\mathbf{K}_{D},2}-\epsilon^{r}_{-\eta\mathbf{K}_{D},4})+i\eta(\epsilon^{r}_{-\eta\mathbf{K}_{D},1}+\epsilon^{r}_{-\eta\mathbf{K}_{D},3}))=0$
$\displaystyle
r=E_{+}:\;\;((\epsilon^{r}_{-\eta\mathbf{K}_{D},2}-\epsilon^{r}_{-\eta\mathbf{K}_{D},4})+i\eta(\epsilon^{r}_{-\eta\mathbf{K}_{D},1}+\epsilon^{r}_{-\eta\mathbf{K}_{D},3}))=0$
$\displaystyle
r=A_{2}:\;\;((\epsilon^{r}_{-\eta\mathbf{K}_{D},2}-\epsilon^{r}_{-\eta\mathbf{K}_{D},4})+i\eta(\epsilon^{r}_{-\eta\mathbf{K}_{D},1}+\epsilon^{r}_{-\eta\mathbf{K}_{D},3}))=0$
$\displaystyle
r=A_{1}:\;\;((\epsilon^{r}_{-\eta\mathbf{K}_{D},2}-\epsilon^{r}_{-\eta\mathbf{K}_{D},4})+i\eta(\epsilon^{r}_{-\eta\mathbf{K}_{D},1}+\epsilon^{r}_{-\eta\mathbf{K}_{D},3}))=-2$
(127)

Hence the electrons at the $\mathbf{K}_{D}$ point only couple with the $A_{1}$
phonons, while the couplings to all other 3 phonon modes are proportional to
the momentum $\mathbf{k}-\mathbf{k}^{\prime}$, being zero when
$\mathbf{k}=\mathbf{k}^{\prime}$. With

$\sum_{\alpha\beta}\psi^{n\star}_{\mathbf{k}\eta\alpha}\sigma^{x}_{\alpha\beta}\psi^{m}_{\mathbf{k}^{\prime}-\eta\beta}=\frac{\eta}{2}(ne^{i\eta\phi}-me^{i\eta\phi^{\prime}})$
(128)

We then have the inter-valley electron-phonon coupling

$\displaystyle H^{op}_{inter-
vall}=\frac{1}{\sqrt{N_{G}}}\gamma_{3}\sum_{m,n=\pm}\sum_{\mathbf{k},\mathbf{k}^{\prime},\eta,\alpha\beta}\frac{b_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},A_{1}}+b_{\eta\mathbf{K}_{D}-\mathbf{k}+\mathbf{k}^{\prime},A_{1}}^{\dagger}}{\sqrt{2M\omega_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},A_{1}}}}\gamma_{\mathbf{k},\eta}^{n\dagger}\gamma_{\mathbf{k}^{\prime},-\eta}^{m}\psi_{\mathbf{k},\eta,\alpha}^{n\star}(\sigma_{x})_{\alpha\beta}\psi_{\mathbf{k}^{\prime},-\eta,\beta}^{m}$
$\displaystyle\approx\frac{1}{2\sqrt{N_{G}}}\gamma_{3}\sum_{m,n=\pm}\sum_{\mathbf{k},\mathbf{k}^{\prime},\eta}\frac{b_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},A_{1}}+b_{\eta\mathbf{K}_{D}-\mathbf{k}+\mathbf{k}^{\prime},A_{1}}^{\dagger}}{\sqrt{2\sqrt{M(\alpha_{1}+\alpha_{2}+\alpha_{3})}}}\gamma_{\mathbf{k},\eta}^{n\dagger}\gamma_{\mathbf{k}^{\prime},-\eta}^{m}\eta(ne^{i\eta\phi}-me^{i\eta\phi^{\prime}})$
(129)

where $\phi,\phi^{\prime}$ are the angles of the momenta
$\mathbf{k},\mathbf{k}^{\prime}$.

  * •

Intra-band ($m=n$) electrons have the matrix element

$m=n:\eta n(e^{i\eta\phi}-e^{i\eta\phi^{\prime}})$ (130)

this is zero for $\phi=\phi^{\prime}$, i.e. for
$\mathbf{k},\mathbf{k}^{\prime}$ being parallel, and is maximum for
$\phi=\pi+\phi^{\prime}$, i.e. for $\mathbf{k},\mathbf{k}^{\prime}$ being
anti-parallel.

  * •

Inter-band ($m\neq n$) electrons have the matrix element

$m=-n:\eta n(e^{i\eta\phi}+e^{i\eta\phi^{\prime}})$ (131)

this is zero for $\phi=\pi+\phi^{\prime}$, i.e. for
$\mathbf{k},\mathbf{k}^{\prime}$ antiparallel and is maximum for
$\phi=\phi^{\prime}$, i.e. for $\mathbf{k},\mathbf{k}^{\prime}$ parallel.

For further reference, we will also need the electron-phonon interaction on
the sublattice basis, _not_ projected into the graphene bands - as we will
need to project it into the TBG bands, which are different from the graphene
bands. This form of the inter-valley electron-phonon interaction reads

$\displaystyle H_{inter-
vall}^{op}\approx\frac{1}{\sqrt{N_{G}}}\gamma_{3}\sum_{\mathbf{k},\mathbf{k}^{\prime},\eta,\alpha\beta}\frac{b_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},A_{1}}+b_{\eta\mathbf{K}_{D}-\mathbf{k}+\mathbf{k}^{\prime},A_{1}}^{\dagger}}{\sqrt{2\sqrt{M(\alpha_{1}+\alpha_{2}+\alpha_{3})}}}c^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},\alpha}(\sigma_{x})_{\alpha\beta}c_{\mathbf{k}^{\prime}-\eta\mathbf{K}_{D},\beta}.$
(132)

###  IV.2 Numerical Calculations of Electron-Phonon Interaction in a Single-
Layer Graphene

In this section, we will describe our numerical calculations of energy
spectrum for both electrons and phonons in a single layer graphene. The energy
spectrum for electrons and phonons are shown in Fig. 5(a) and (b),
respectively, which are consistent with the literature [59, 60, 61, 62, 63].
The parameters for the calculation of energy spectrum of electrons and phonons
are summarized in table 2. For the phonon dispersion, we label four in-plane
phonon modes as the modes 1, 2, 3 and 4, with their irreducible
representations (irreps) at high symmetry momenta shown in table 3. Here we
only consider the $C_{6v}$ symmetry group for the in-plane modes of graphene,
as the z-direction symmetry will be lost when considering TBG. The phonon
modes 1 and 2 are degenerate have zero energy at ${\bf\Gamma}$, so they
describe the acoustic phonon modes with $E_{1}$ irrep of $C_{6v}$, while the
phonon modes 3 and 4 are degenerate at ${\bf\Gamma}$ and describe two optical
phonon modes with $E_{2}$ irrep. These degeneracies are split away from
$\Gamma$. At $\mathbf{K}_{1}$ (the same as $\mathbf{K}_{D}$), the modes 2 and
3 become degenerate, forming the $E$ irrep, while the mode 1 belongs to
$A_{2}$ irrep and the mode 4 to $A_{1}$ irrep, as shown in Fig. 5 (b).

Figure 5: (a) Energy dispersion for electrons in single-layer graphene. Here the blue dashed lines are for the nearest neighbor approximation, while the red solid lines are for the full long-range hopping. The inset shows the Brillouin zone of graphene. Here $K_{1}$ is just $K_{D}$ defined in the main text. (b) Phonon dispersion in single-layer graphene.  Quantity | Value
---|---
$a_{0}$ |  $2.46$Å
$r_{0}$ | $0.148a_{0}$
$V_{pp\pi}$ | $-2.7eV$
$K_{t0}$ | $2.94\times 10^{5}cm^{-2}$
$K_{r0}$ | $5.15\times 10^{5}cm^{-2}$
$K_{t1}$ | $-8.58\times 10^{4}cm^{-2}$
$K_{r1}$ | $1.33\times 10^{5}cm^{-2}$
$K_{2}$ | 0
Table 2: Table for parameters used for calculating energy spectrum of electrons and phonons in single-layer graphene.  Phonon modes | 1 | 2 | 3 | 4
---|---|---|---|---
$\Gamma$ ($C_{6v}$) | $E_{1}$ | $E_{1}$ | $E_{2}$ | $E_{2}$
$K$ ($C_{3v}$) | $A_{2}$ | $E$ | $E$ | $A_{1}$
Table 3: Symmetry representations for the in-plane phonons at high symmetry
momenta $\Gamma$ and $\mathbf{K}$ in single-layer graphene.

For the numerical evaluation of electron-phonon interaction, we introduce the
interaction parameter $g$ as (Since we only consider a single-layer graphene
here, we have dropped the layer index)

$\displaystyle
H_{ep}=\frac{1}{\sqrt{N_{G}}}\sum_{\mathbf{p}\mathbf{p}^{\prime},\alpha\beta
r}g^{r}_{\alpha\beta}(\mathbf{p},\mathbf{p}^{\prime})c^{\dagger}_{\mathbf{p},\alpha}c_{\mathbf{p}^{\prime},\beta}(b_{\mathbf{p}-\mathbf{p}^{\prime},r}+b^{\dagger}_{\mathbf{p}^{\prime}-\mathbf{p},r})$
$\displaystyle=\frac{1}{\sqrt{N_{G}}}\sum_{\mathbf{p}\mathbf{p}^{\prime},nm,r}g^{s}_{nm}(\mathbf{p},\mathbf{p}^{\prime})\gamma^{n\dagger}_{\mathbf{p}}\gamma^{m}_{\mathbf{p}^{\prime}}(b_{\mathbf{p}-\mathbf{p}^{\prime},r}+b^{\dagger}_{\mathbf{p}^{\prime}-\mathbf{p},r})$
(133)

where $N_{G}$ is the number of unit cells, $\alpha,\beta=A,B$, $n,m$ label the
eigen-states of electrons and $s$ labels the eigen-state of phonon modes. Here
$g^{r}_{\alpha\beta}(\mathbf{p},\mathbf{p}^{\prime})$ is in the sub-lattice
basis for electrons while $g^{r}_{nm}(\mathbf{p},\mathbf{p}^{\prime})$ is in
the eigen-state basis. Explicitly,
$g^{r}_{\alpha\beta}(\mathbf{p},\mathbf{p}^{\prime})$ is given by

$\displaystyle
g^{r}_{\alpha\beta}(\mathbf{p},\mathbf{p}^{\prime})=i\sqrt{\frac{\hbar}{2M\omega_{\mathbf{p}-\mathbf{p}^{\prime},r}}}\sum_{G}e^{i\mathbf{G}\cdot(\tau_{\alpha}-\tau_{\beta})}\left(t(\mathbf{p}^{\prime}+\mathbf{G})(\mathbf{p}^{\prime}+\mathbf{G})\cdot{\bf\epsilon}^{r}_{\mathbf{p}-\mathbf{p}^{\prime},\alpha}-t(\mathbf{p}+\mathbf{G})(\mathbf{p}+\mathbf{G})\cdot{\bf\epsilon}^{r}_{\mathbf{p}-\mathbf{p}^{\prime},\beta}\right)$
(134)

which can be directly obtained from substituting the phonon eigen-mode
expansion Eq. (125) into Eq. (44). Here we have recovered the constant $\hbar$
for the correct unit.

$\displaystyle
g^{r}_{nm}(\mathbf{p},\mathbf{p}^{\prime})=\sum_{\alpha\beta}\psi^{n*}_{\alpha,\mathbf{p}}g^{r}_{\alpha\beta}(\mathbf{p},\mathbf{p}^{\prime})\psi^{m}_{\beta,\mathbf{p}^{\prime}}.$
(135)

With the expansion of the momentum around $\eta\mathbf{K}_{D}$, one can
directly relate the function $g$ to the parameters $\gamma_{1,2,3}$ defined in
the early sections. We can give an estimate of the magnitude of $g$ from the
value of $\gamma_{3}$ as
$|\gamma_{3}|\sqrt{\frac{\hbar}{2M\omega_{op,K}}}=|\gamma_{3}|a\sqrt{\frac{\hbar^{2}}{2m_{e}a_{B}^{2}}\frac{1}{\hbar\omega_{op,K}}\frac{m_{e}}{M}\frac{a_{B}^{2}}{a^{2}}}$,
where $\frac{\hbar^{2}}{2m_{e}a_{B}^{2}}$ is Rydberg energy,
$\hbar\omega_{op,K}$ is optical phonon energy at $K$, $m_{e}/M$ is the mass
ration between electron and carbon atom, $a_{B}$ is the Bohr radius and
$a=a_{0}/\sqrt{3}$ is the NN distance between two carbon atoms in graphene.
With $|\gamma_{3}|\sim 17eV/$Å, $a\sim 1.42\AA$ and the dimensionless quantity
$\sqrt{\frac{\hbar^{2}}{2m_{e}a_{B}^{2}}\frac{1}{\hbar\omega_{op,\Gamma}}\frac{m_{e}}{M}\frac{a_{B}^{2}}{a^{2}}}\sim
0.02$, we find the magnitude of $g$ is around $0.56$ eV. To evaluate $g$
numerically, it is more efficient to perform the calculation in the real space
instead of the momentum space since the hopping function $t$ decays
exponentially in the real space, but only in power-law in the momentum space.
To get a real space expression, we can use the identity

$\displaystyle\sum_{\mathbf{G}}e^{i\mathbf{G}\cdot(\tau_{\alpha}-\tau_{\beta})}i(\mathbf{p}+\mathbf{G})t(\mathbf{p}+\mathbf{G})=\sum_{\mathbf{R}}e^{-i\mathbf{p}\cdot(\mathbf{R}+\tau_{\alpha}-\tau_{\beta})}\nabla
t(\mathbf{r})|_{\mathbf{r}=\mathbf{R}+\tau_{\alpha}-\tau_{\beta}},$ (136)

and obtain

$\displaystyle
g^{r}_{\alpha\beta}(\mathbf{p},\mathbf{p}^{\prime})=\sqrt{\frac{\hbar}{2M\omega_{r,\mathbf{p}-\mathbf{p}^{\prime}}}}\sum_{\mathbf{R}}\left(e^{-i\mathbf{p}^{\prime}\cdot\delta_{\alpha\beta}}{\bf\epsilon}^{r}_{\alpha,\mathbf{p}-\mathbf{p}^{\prime}}-e^{-i\mathbf{p}\cdot\delta_{\alpha\beta}}{\bf\epsilon}^{r}_{\beta,\mathbf{p}-\mathbf{p}^{\prime}}\right)\cdot\left.\nabla
t(\mathbf{r})\right|_{\mathbf{r}=\delta_{\alpha\beta}}$
$\displaystyle=-\sqrt{\frac{\hbar}{2M\omega_{r,\mathbf{p}-\mathbf{p}^{\prime}}}}\sum_{\mathbf{R}}\frac{t(\delta_{\alpha\beta})}{r_{0}}\left(e^{-i\mathbf{p}^{\prime}\cdot(\mathbf{R}+\tau_{\alpha}-\tau_{\beta})}{\bf\epsilon}^{r}_{\alpha,\mathbf{p}-\mathbf{p}^{\prime}}-e^{-i\mathbf{p}\cdot(\mathbf{R}+\tau_{\alpha}-\tau_{\beta})}{\bf\epsilon}^{r}_{\beta,\mathbf{p}-\mathbf{p}^{\prime}}\right)\cdot\frac{\delta_{\alpha\beta}}{|\delta_{\alpha\beta}|},$
(137)

where $\delta_{\alpha\beta}=\mathbf{R}+\tau_{\alpha}-\tau_{\beta}$. This is
the expression for our numerical simulations.

Figure 6:  The electron-phonon interaction parameter
$g^{r}_{nm}(\mathbf{p},\mathbf{p}^{\prime})$ as a function of $p^{\prime}$ for
a fixed $\mathbf{p}=\mathbf{K}_{1}+\delta\mathbf{p}$, where
$\delta\mathbf{p}=\frac{1}{a_{0}}(0.001,0)$. (a), (b), (c) and (d) are for
four phonon modes $r=1,2,3,4$, respectively, and $n,m=\pm$. The red hexagon
shows the Brillouin zone. Here the momenta $\mathbf{p}$ and
$\mathbf{p}^{\prime}$ are in the unit of $1/a_{0}$ with $a_{0}=2.46$Å.  Figure
7:  Zoom in for the electron-phonon parameter
$g(\mathbf{p},\mathbf{p}^{\prime})$ as a function of $\mathbf{p}^{\prime}$
with a fixed $\mathbf{p}=\mathbf{K}_{1}+\delta\mathbf{p}$
($\delta\mathbf{p}=\frac{1}{a_{0}}(0.001,0)$) for the phonon mode $r=4$
(inter-valley scattering). (a) and (b) show the zoom in of
$g^{r=4}_{\alpha\beta}(\mathbf{p},\mathbf{p}^{\prime})$ ($\alpha,\beta=A,B$)
and $g^{r=4}_{nm}(\mathbf{p},\mathbf{p}^{\prime})$ ($n,m=\pm$) for
$\mathbf{p}^{\prime}$ around $\mathbf{K}_{2}^{\prime}$. (c) shows the
analytical solutions for the electron-phonon parameters
$g^{r=4}_{nm}(\mathbf{p},\mathbf{p}^{\prime})$ ($n,m=\pm$) for
$\mathbf{p}^{\prime}$ around $\mathbf{K}_{2}^{\prime}$ for the phonon mode 4
from Eq. (IV.1). Here the momenta $\mathbf{p}$ and $\mathbf{p}^{\prime}$ are
in the unit of $1/a_{0}$ with $a_{0}=2.46$Å.

Next let’s discuss our numerical results in Figs. 6 and 7 based on the full
expression Eq. (137). Fig. 6 shows
$g^{r}_{nm}(\mathbf{p}=\mathbf{K}_{1}+\delta\mathbf{p},\mathbf{p}^{\prime})$
($\mathbf{K}_{1}=\mathbf{K}_{D}$) as a function of $\mathbf{p}^{\prime}$ in
the whole Brillouin zone for different phonon modes (a small
$\delta\mathbf{p}=\frac{1}{a_{0}}(0.001,0)$ is chosen to get rid of the
degeneracy of electron states at the Dirac point). Since electron Fermi
surface is only around $\pm\mathbf{K}_{D}$, we can focus on the corners of the
Brillouin zone (BZ), which is shown by the red hexagon in Fig. 6. Furthermore,
we focus on the inter-valley e-ph interaction, which means
$\mathbf{p}^{\prime}$ around
$\mathbf{K}_{1}^{\prime},\mathbf{K}_{2}^{\prime},\mathbf{K}_{3}^{\prime}$, as
we fix $\mathbf{p}=\mathbf{K}_{1}+\delta\mathbf{p}$ and the phonon momentum is
$\mathbf{q}=\mathbf{p}-\mathbf{p}^{\prime}$. We first notice that the e-ph
interaction for the phonon modes 1, 2 and 3 are approaching zero when the
momentum $\mathbf{p}^{\prime}$ is close to
$\mathbf{K}_{1}^{\prime},\mathbf{K}_{2}^{\prime},\mathbf{K}_{3}^{\prime}$,
which is consistent with our analytical solution (See Eq. IV.1) in Sec. IV.1.
Thus, only phonon modes 4 couple to electron for the $K$-phonons. The phonon
mode 4 belongs to $E_{2}$ irrep of $C_{6v}$ at ${\bf\Gamma}$ and $A_{1}$ irrep
of $C_{3v}$ at $\mathbf{K}_{D}$, and we show the zoom-in for
$\mathbf{p}^{\prime}$ around $\mathbf{K}_{2}^{\prime}$ (inter-valley
scattering) in Fig. 7(a) and (b). Here (a) and (b) correspond to the parameter
$g$ in the sublattice basis and the eigen-state basis, respectively, for the
electron part. On the sublattice basis, we find $g$ is non-zero for the inter-
sublattice components, while its value approaches zero for the intra-valley
components, consistent with the analytical results shown in Eq. (132). On the
eigen-state basis, we find $g$ shows a singular behavior at
$\mathbf{K}_{2}^{\prime}$ for the inter-valley scattering. In Fig. 7(c), we
also show the analytical results of inter-valley scattering from Eq. (IV.1),
which is consistent with the numerical results in Fig. 7 (b) on the eigen-
state basis.

##  V Review of Bistritzer-MacDonald model for Twisted Bilayer Graphene

Here we first briefly review the Bistritzer-MacDonald (BM) model [56] of
magic-angle twisted bilayer graphene (MATBG) by introducing our notations.
Readers may refer to the supplementary materials of Ref. [71] for more
details.

###  V.1 Basis and Single-particle Hamiltonian

As shown in Fig. 8(a), the low-energy states around the Dirac points at
$\mathbf{K}_{D}$ of the two graphene layers form the Moiré Brillouin zone
(MBZ) in the valley $\mathbf{K}_{D}$. Similarly, those at the
$\mathbf{K}_{D}^{\prime}$ momenta of the two layers form another MBZ in the
valley $\mathbf{K}_{D}^{\prime}$. Thus, in the BM model there are two
independent sets of basis from two valleys. We use the index $\eta$ ($=+$ for
$\mathbf{K}_{D}$ and $-$ for $\mathbf{K}_{D}^{\prime}$) to label the two
graphene valleys. We denote the basis as
$c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}$, where $\mathbf{k}$ is a momentum in
the MBZ, $\mathbf{Q}$ takes values in the lattice shown in Fig. 8(b)
(illustrated below), $\alpha=1,2$ represents the graphene sublattice,
$\eta=\pm$ represents the graphene valley, and $s=\uparrow,\downarrow$ is the
spin index. To construct the MBZ and the $\mathbf{Q}$ lattices in Fig. 8(b),
we define ${\bf q}_{1}={\bf K}_{D}^{-}-{\bf K}_{D}^{+}=k_{\theta}(0,-1),{\bf
q}_{2}=\hat{C}_{3z}{\bf q}_{1}={\bf K}_{D}^{-}+{\bf G}_{1}^{2}-{\bf
K}_{D}^{+}-{\bf G}_{1}^{1}=k_{\theta}(\frac{\sqrt{3}}{2},\frac{1}{2}),{\bf
q}_{3}=C_{3z}^{2}{\bf q}_{1}={\bf K}_{D}^{-}+{\bf G}_{2}^{2}-{\bf
K}_{D}^{+}-{\bf G}_{2}^{1}=k_{\theta}(-\frac{\sqrt{3}}{2},\frac{1}{2})$, where
$k_{\theta}=2|{\bf K}_{D}|\sin\frac{\theta}{2}$ and $\theta$ is the twist
angle between two layers. Here $\mathbf{K}_{D}^{+}$ and $\mathbf{K}_{D}^{-}$
are the $\mathbf{K}_{D}$ momentum of the top and bottom layers There are two
types of $\mathbf{Q}$ lattices in Fig. 8(b): the blue lattice
$\mathcal{Q}_{+}=\\{\mathbf{q}_{2}+n_{1}\mathbf{b}_{M1}+n_{2}\mathbf{b}_{M2}\
|\ n_{1,2}\in\mathbb{Z}\\}$ and the red lattice
$\mathcal{Q}_{-}=\\{-\mathbf{q}_{2}+n_{1}\mathbf{b}_{M1}+n_{2}\mathbf{b}_{M2}\
|\ n_{1,2}\in\mathbb{Z}\\}$, where $\mathbf{b}_{M1},\mathbf{b}_{M2}$ are Moiré
reciprocal lattice basis

$\mathbf{b}_{M1}=\mathbf{q}_{2}-\mathbf{q}_{1}=k_{\theta}(\frac{\sqrt{3}}{2},\frac{3}{2}),\qquad\mathbf{b}_{M2}=\mathbf{q}_{3}-\mathbf{q}_{1}=k_{\theta}(-\frac{\sqrt{3}}{2},\frac{3}{2}).$
(138)

For $\mathbf{Q}\in\mathcal{Q}_{+}$, the basis is defined as

$(\mathbf{Q}\in\mathcal{Q}_{+})\qquad
c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}=\begin{cases}\displaystyle\frac{1}{\sqrt{N_{\rm
G}}}\sum_{\mathbf{R}\in\mathrm{top}}e^{i(\mathbf{K}_{D}^{+}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}+\tau_{\alpha})}c_{\mathbf{R},\alpha,s}^{\dagger}\qquad&\text{if}\;\eta=+\\\
\displaystyle\frac{1}{\sqrt{N_{\rm
G}}}\sum_{\mathbf{R}^{\prime}\in\mathrm{bottom}}e^{i(-\mathbf{K}_{D}^{-}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}^{\prime}+\tau_{\alpha}^{\prime})}c_{\mathbf{R}^{\prime},\alpha,s}^{\dagger}\qquad&\text{if}\;\eta=-\end{cases},$
(139)

and for $\mathbf{Q}\in\mathcal{Q}_{-}$, the basis is defined as

$(\mathbf{Q}\in\mathcal{Q}_{-})\qquad
c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}=\begin{cases}\displaystyle\frac{1}{\sqrt{N_{\rm
G}}}\sum_{\mathbf{R}^{\prime}\in\mathrm{bottom}}e^{i(\mathbf{K}_{D}^{-}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}^{\prime}+\tau_{\alpha}^{\prime})}c_{\mathbf{R}^{\prime},\alpha,s}^{\dagger}\qquad&\text{if}\;\eta=+\\\
\displaystyle\frac{1}{\sqrt{N_{\rm
G}}}\sum_{\mathbf{R}\in\mathrm{top}}e^{i(-\mathbf{K}_{D}^{+}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}+\tau_{\alpha})}c_{\mathbf{R},\alpha,s}^{\dagger}\qquad&\text{if}\;\eta=-\end{cases}\
.$ (140)

Here $N_{\rm G}$ is the number of graphene unit cells in each layer,
$\mathbf{R}$ and $\mathbf{R}^{\prime}$ index graphene lattices in the the top
and bottom layers, respectively, $\tau_{\alpha}$ and $\tau_{\alpha}^{\prime}$
are the sublattice vectors of the two layers, respectively,
$\mathbf{K}_{D}^{+}$ and $\mathbf{K}_{D}^{-}$ are the $\mathbf{K}_{D}$
momentum of the top and bottom layers, respectively, and
$c_{\mathbf{R},\alpha,s}$ ($c_{\mathbf{R}^{\prime},\alpha,s}$) is the fermion
operator with spin $s$ at the atom site $\mathbf{R}+\tau_{\alpha}$
($\mathbf{R}^{\prime}+\tau_{\alpha}^{\prime}$). The $\mathcal{Q}_{+}$
($\mathcal{Q}_{-}$) lattice is defined in such a way that
$\eta\mathbf{K}_{D}^{+}+\mathbf{k}-\mathbf{Q}$
($\eta\mathbf{K}_{D}^{-}+\mathbf{k}-\mathbf{Q}$) with
$\mathbf{Q}\in\mathcal{Q}_{+}$ ($\mathbf{Q}\in\mathcal{Q}_{-}$) is the Dirac
point position $\eta\mathbf{K}_{D}^{+}$ ($\eta\mathbf{K}_{D}^{-}$) when
$\mathbf{k}$ equals to the high symmetry point $\eta\mathbf{K}_{M}$
($\eta\mathbf{K}_{M}^{\prime}$) of the MBZ, as shown in Fig. 8(b).

Figure 8: The moiré BZ and the $\mathbf{Q}$ lattice. (a) The blue and red
hexagons represent the BZ’s of the top and bottom graphene layers,
respectively. The blue and red dots represent the Diract points
($\mathbf{K}_{D}$ momentum) from the two layers. The vectors
$\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}$ connecting the Dirac points of
the two layers span the moiré BZ in the valley $\mathbf{K}_{D}$ (shaded with
light green). Similarly, the Dirac points at the $\mathbf{K}_{D}^{\prime}$
momentum form another Moiré BZ in the valley $\mathbf{K}_{D}^{\prime}$ (shaded
with light yellow). (b) We denote the blue and red $\mathbf{Q}$ lattices as
$\mathcal{Q}_{+}$ and $\mathcal{Q}_{-}$ respectively. In the BM model, the
states in the valley $\mathbf{K}_{D}$ ($\mathbf{K}_{D}^{\prime}$) of the top
layer and in the valley $\mathbf{K}_{D}^{\prime}$ ($\mathbf{K}_{D}$) of the
bottom layer contribute to the lattice $\mathcal{Q}_{+}$ ($\mathcal{Q}_{-}$).

The BM model is given by

${\hat{H}}_{\rm 0}=\sum_{\eta s}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\alpha\alpha^{\prime}}\sum_{\mathbf{Q},\mathbf{Q}^{\prime}}h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(\eta)}(\mathbf{k})c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}c_{\mathbf{k},\mathbf{Q}^{\prime},\alpha^{\prime},\eta,s}$
(141)

where the single particle Hamiltonian reads

$h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(+)}(\mathbf{k})=v_{F}(\mathbf{k}-\mathbf{Q})\cdot\bm{\sigma}\delta_{\mathbf{Q},\mathbf{Q}^{\prime}}+\sum_{j=1}^{3}[T_{j}]_{\alpha\alpha^{\prime}}\delta_{\mathbf{Q},\mathbf{Q}^{\prime}\pm\mathbf{q}_{j}},\qquad
h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(-)}(\mathbf{k})=h_{-\mathbf{Q}\alpha,-\mathbf{Q}^{\prime}\alpha^{\prime}}^{(+)*}(-\mathbf{k}),$
(142)
$T_{j}=w_{0}\sigma_{0}+w_{1}\sigma_{x}\cos\frac{2\pi(j-1)}{3}+w_{1}\sigma_{y}\sin\frac{2\pi(j-1)}{3}\
.$ (143)

Here $w_{0}$ and $w_{1}$ are the interlayer couplings in the AA-stacking and
AB-stacking regions, respectively. In this work, we adopt the parameters
$v_{F}=5.944\mathrm{eV\cdot\mathring{A}}$,
$k_{\theta}=0.0312\mathrm{\mathring{A}^{-1}}$, $w_{1}=110\mathrm{meV}$. The
relation
$h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(-)}(\mathbf{k})=h_{-\mathbf{Q}\alpha,-\mathbf{Q}^{\prime}\alpha^{\prime}}^{(+)*}(-\mathbf{k})$
is due to the time-reversal symmetry that transform the two valleys to each
other. The single-particle Hamiltonian (upon to a unitary transformation [71])
is periodic with respect to the reciprocal lattice vectors $\mathbf{b}_{M1}$,
$\mathbf{b}_{M2}$. The MBZ and high symmetry momenta are defined in Fig. 8(b).
The corresponding real space unit cell and maximal Wyckoff positions are shown
in Fig. 8(c). The $1a$ and $2c$ positions correspond to the AA-stacking and
AB-stacking regions, respectively. We denote the real space lattice basis as
$\mathbf{a}_{M1}$, $\mathbf{a}_{M2}$, they satisfy
$\mathbf{a}_{Mi}\cdot\mathbf{b}_{Mj}=2\pi\delta_{ij}$. We plot the band
structures in the valley $\eta=+$ for different twisted angles around the
magic angle in the chiral limit $w_{0}=0$ in Fig. 9 .

The fact that the band structure is labeled by $\mathbf{k}$ in the MBZ implies
that $\mathbf{k}$ labels the eigenvalues of translation operators of the Moiré
lattice. We hence define the translation operator $T_{\mathbf{R}}$ as

$T_{\mathbf{R}}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}T_{\mathbf{R}}^{-1}=e^{-i\mathbf{k}\cdot\mathbf{R}}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}\
,$ (144)

where $\mathbf{R}=n_{1}\mathbf{a}_{M1}+n_{2}\mathbf{a}_{M2}$, with
$n_{1,2}\in\mathbb{Z}$, is a Moiré lattice vector. One should not confuse
$T_{\mathbf{R}}$ with the time-reversal symmetry ($\hat{T}$) defined in the
early section. We now verify Eq. (144) at the commensurate twist angles, where
$\mathbf{a}_{M1}$ and $\mathbf{a}_{M2}$ are integer linear combinations of the
microscopic graphene lattice vectors $\mathbf{a}_{1},\mathbf{a}_{2}$.
$T_{\mathbf{R}}$ can be defined as translation acting on the atomic orbitals:
$\mathbf{R}^{\prime}\to\mathbf{R}^{\prime}+\mathbf{R}$, leading to

$\displaystyle
T_{\mathbf{R}}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}T_{\mathbf{R}}^{-1}=\frac{1}{\sqrt{N_{\rm
G}}}\sum_{\mathbf{R}^{\prime}\in
l}e^{i(\eta\mathbf{K}_{D}^{l}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}^{\prime}+\tau_{\alpha})}c_{\mathbf{R}^{\prime}+\mathbf{R},\alpha,s}^{\dagger}=\frac{1}{\sqrt{N_{\rm
G}}}\sum_{\mathbf{R}^{\prime}\in
l}e^{-i(\eta\mathbf{K}_{D}^{l}+\mathbf{k}-\mathbf{Q})\cdot\mathbf{R}}e^{i(\eta\mathbf{K}_{D}^{l}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}^{\prime}+\tau_{\alpha})}c_{\mathbf{R}^{\prime},\alpha,s}^{\dagger}\
,$ (145)

where $l=\eta$ ($-\eta$) for $\mathbf{Q}\in\mathcal{Q}_{+}$
($\mathcal{Q}_{-}$). For commensurate angle,
$\eta\mathbf{K}_{D}^{l}-\mathbf{Q}$ is a Moiré reciprocal lattice and hence
$e^{-i(\eta\mathbf{K}_{D}^{l}+\mathbf{k}-\mathbf{Q})\cdot\mathbf{R}}=e^{-i\mathbf{k}\cdot\mathbf{R}}$.
Then Eq. (144) is justified at commensurate angles. We emphasize that Eq.
(144) is also well-defined even at non-commensurate angles.

Figure 9:  The band structure in the valley $\eta=+$ of TBG near the magic
angle in the chiral limit ($w_{0}=0$). (a), (b) and (c) are for three twist
angles $\theta=1.063152^{\circ}$, $\theta=1.06^{\circ}$ and
$\theta=1.05^{\circ}$. (d), (e) and (f) show the zoom in for the dispersion of
the flat bands for these three angles. One can see the profiles of the
dispersion look similar for these three angles, while the band width changes
$0.02$ meV to 2 meV. The other parameters of the BM model are given by
$v_{F}=5.944\mathrm{eV}\cdot\mathring{A}$,
$|\mathbf{K}|=1.703\mathrm{\mathring{A}^{-1}}$, $w_{1}=110\mathrm{meV}$.

###  V.2 Symmetries of the BM model

The single-particle Hamiltonian in each valley has the symmetry of the
magnetic space group $P6^{\prime}2^{\prime}2$ (# 177.151 in the BNS setting),
which is generated by $\hat{C}_{2z}\hat{T}$, $\hat{C}_{3z}$, $\hat{C}_{2x}$
and translation symmetries. Since the two valleys are related by time-reversal
symmetry $\hat{T}$, the total system also has the $\hat{C}_{2z}$ symmetry
(product of $\hat{C}_{2z}T$ and $\hat{T}$). The full crystalline symmetries of
the two valleys form the space group $P622$ (#177), which is generated by
$\hat{C}_{6z}$, $\hat{C}_{2x}$, and translation symmetries. Thus, these
symmetry operators satisfy

$[\hat{C}_{3z},H_{0}]=[\hat{C}_{2z},H_{0}]=[\hat{C}_{2x},H_{0}]=[\hat{T},H_{0}]=0,$
(146)

and we write the symmetry action on the fermion operators as

$\hat{g}c_{\mathbf{k},\mathbf{Q},\eta,\alpha,s}^{\dagger}\hat{g}^{-1}=\sum_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime}}c_{\hat{g}\mathbf{k},\mathbf{Q}^{\prime},\eta^{\prime},\alpha^{\prime},s}^{\dagger}[D(\hat{g})]_{\mathbf{Q}^{\prime}\eta^{\prime}\alpha^{\prime},\mathbf{Q}\eta\alpha}\
.$ (147)

The $D$ matrices for the above discrete symmetries are given by

$[D(C_{3z})]_{\mathbf{Q}^{\prime}\eta^{\prime}\beta,\mathbf{Q}\eta\alpha}=\delta_{\mathbf{Q}^{\prime},C_{3z}\mathbf{Q}}\delta_{\eta^{\prime},\eta}(e^{i\eta\frac{2\pi}{3}\sigma_{z}})_{\beta\alpha}\;\;\;[D(C_{2x})]_{\mathbf{Q}^{\prime}\eta^{\prime}\beta,\mathbf{Q}\eta\alpha}=\delta_{\mathbf{Q}^{\prime},C_{2x}\mathbf{Q}}\delta_{\eta^{\prime},\eta}(\sigma_{x})_{\beta\alpha},$
(148)
$[D(C_{2z})]_{\mathbf{Q}^{\prime}\eta^{\prime}\beta,\mathbf{Q}\eta\alpha}=\delta_{\mathbf{Q}^{\prime},-\mathbf{Q}}\delta_{\eta^{\prime},-\eta}(\sigma_{x})_{\beta\alpha},\;\;\;[D(T)]_{\mathbf{Q}^{\prime}\eta^{\prime}\beta,\mathbf{Q}\eta\alpha}=\delta_{\mathbf{Q}^{\prime},-\mathbf{Q}}\delta_{\eta^{\prime},-\eta}\delta_{\beta,\alpha,}$
(149)
$[D(C_{2z}T)]_{\mathbf{Q}^{\prime}\eta^{\prime}\beta,\mathbf{Q}\eta\alpha}=[D(C_{2z})D(T)]_{\mathbf{Q}^{\prime}\eta^{\prime}\beta,\mathbf{Q}\eta\alpha}=\delta_{\mathbf{Q}^{\prime},\mathbf{Q}}\delta_{\eta^{\prime},\eta}(\sigma_{x})_{\beta,\alpha}.$
(150)

Here and below we will use $\sigma_{x,y,z}$ ($\sigma_{0}$) and $\tau_{x,y,z}$
($\tau_{0}$) for the Pauli (identity) matrices in the sub-lattice and valley
spaces. For the sublattice, we choose $\alpha=A,B$ to be $1,2$ for $\sigma$
matrix.

Besides the crystalline symmetries, the BM model also has a unitary particle-
hole (PH) symmetry $\hat{P}$ that anti-commutes with the single-particle
Hamiltonian, $\\{\hat{P},H_{0}\\}=0$. $\hat{P}$ transforms $\mathbf{k}$ to
$-\mathbf{k}$ and the corresponding $D$ matrix reads

$Pc^{\dagger}_{\mathbf{k},\mathbf{Q},\eta,\alpha,s}P^{-1}=\sum_{\mathbf{Q}^{\prime}\eta^{\prime}\beta}[D(P)]_{\mathbf{Q}^{\prime}\eta^{\prime}\beta,\mathbf{Q}\eta\alpha}c^{\dagger}_{-\mathbf{k},\mathbf{Q}^{\prime},\eta^{\prime},\beta,s},\;\;\;\quad[D(P)]_{\mathbf{Q}^{\prime}\eta^{\prime}\beta,\mathbf{Q}\eta\alpha}=\delta_{\mathbf{Q}^{\prime},-\mathbf{Q}}\delta_{\eta^{\prime},\eta}\delta_{\beta,\alpha}\zeta_{\mathbf{Q}}\
,$ (151)

where $\zeta_{\mathbf{Q}}=\pm 1$ for $\mathbf{Q}\in\mathcal{Q}_{\pm}$,
respectively. Note that $P$ transforms creation operators to creation
operators (rather than annihilation operators), and maps sites
$\mathbf{Q}\in\mathcal{Q}_{\pm}$ into $-\mathbf{Q}\in\mathcal{Q}_{\mp}$. The
PH transformation $P$ satisfies

$P^{2}=-1,\qquad[P,C_{3z}]=0,\qquad\\{P,C_{2x}\\}=0,\qquad\\{P,C_{2z}\\}=0,\qquad\\{P,T\\}=0,\qquad[P,C_{2z}T]=0.$
(152)

The particle-hole symmetry will be broken if the $\theta$-dependence of the
single-layer Hamiltonian or quadratic terms in $\mathbf{k}$ of the single-
layer Hamiltonian are taken into account.

Besides the above discrete symmetry, the BM model also possesses the
continuous U(2)$\times$U(2) spin-charge rotation symmetry, and the
corresponding symmetry operators are given by

$\hat{S}^{ab}=\sum_{\mathbf{k}}(\tau^{a})_{\eta\eta^{\prime}}(s^{b})_{ss^{\prime}}c_{\mathbf{k},\mathbf{Q},\eta,\alpha,s}^{\dagger}c_{\mathbf{k},\mathbf{Q},\eta^{\prime},\alpha,s^{\prime}}\
,\qquad(a=0,z,\quad b=0,x,y,z)\ ,$ (153)

where $\tau$ and $s$ are Pauli matrices for valley and spin degree of freedom.

If, furthermore, $w_{0}=0$, the model acquires an effective chiral symmetry
$\hat{C}$ that anti-commute with the single-particle Hamiltonian [72].
$\hat{C}$ is referred to as the first chiral symmetry in Refs. [58, 57].
$\hat{C}$ leaves $\mathbf{k}$ invariant and its $D$ matrix reads

$[D(\hat{C})]_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime},\mathbf{Q}\alpha\eta}=\delta_{\mathbf{Q}^{\prime},\mathbf{Q}}[\sigma_{z}]_{\alpha^{\prime}\alpha}[\tau_{0}]_{\eta^{\prime}\eta}\
.$ (154)

In the first-quantized formalism the algebra between $C$ and other symmetries
are given by

$C^{2}=1,\qquad[C,T]=0,\qquad[C,C_{3z}]=0,\qquad\\{C,C_{2x}\\}=0,\qquad\\{C,C_{2z}T\\}=0,\qquad[C,P]=0\
.$ (155)

##  VI Intra-Layer Inter-Valley Electron-phonon Interaction in Twisted Bilayer
Graphene

###  VI.1 Projected Intra-Layer Inter-Valley Electron-Phonon Interaction for
the Flat Bands

We define the short-hand notation for the $\mathbf{Q}$ lattice in the layer
$l=\pm$ (top/bottom) and valley $\eta=\pm$ as

$\mathcal{Q}_{l\eta}=\\{l\eta\mathbf{q}_{2}+n_{1}\mathbf{b}_{M1}+n_{2}\mathbf{b}_{M2}\
|\ n_{1,2}\in\mathbb{Z}\\},$ (156)

which is consistent with the definition in the previous Sec. V.1. We have the
equivalence between the operators $c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}$ in
Eq. (139) in the MBZ and the ones $c_{\mathbf{k},\alpha,l,s}$ in a single
layer graphene,

$c_{\mathbf{k},\mathbf{Q}_{l\eta},\alpha,\eta,s}=c_{\eta\mathbf{K}_{D}^{l}+\mathbf{k}-\mathbf{Q}_{l\eta},\alpha,l,s},\;\;\;k\in
MBZ$ (157)

where $\mathbf{K}_{D}^{l}$ is the $K$ momentum in the layer $l$.

With the electron-phonon Hamiltonian now summed over the two layers $l=\pm$,
and with $\mathbf{K}_{D}^{l}$ being the graphene Dirac momentum in the layer
$l$, we have

$\displaystyle H_{inter-
vall}^{op}\approx\frac{\gamma_{3}}{\sqrt{2N_{G}{\sqrt{M(\alpha_{1}+\alpha_{2}+\alpha_{3})}}}}\sum_{\mathbf{k},\mathbf{k}^{\prime},\eta,\alpha\beta,l,s}(b_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime},l,A_{1}}+b_{\eta\mathbf{K}_{D}-\mathbf{k}+\mathbf{k}^{\prime},l,A_{1}}^{\dagger})c^{\dagger}_{\mathbf{k}+\eta\mathbf{K}_{D},l,\alpha,s}(\sigma_{x})_{\alpha\beta}c_{\mathbf{k}^{\prime}-\eta\mathbf{K}_{D},l,\beta,s}$

where we have added the layer $l=\pm$ index as well as the spin index in the
electrons. In this summation, $\mathbf{k}$ is summed over the vicinity of the
Dirac node, over a radius much larger than the first MBZ (but much smaller
smaller than $\mathbf{K}_{D}$). We hence break the sum over the momentum into

$\sum_{\mathbf{k}}\rightarrow\sum_{\mathbf{k}\in
MBZ}\sum_{\mathbf{Q}_{l\eta}}$ (159)

where the lattice $\mathbf{Q}_{l,\eta}$ depends on $\eta,l$ for electrons
around the $\eta\mathbf{K}_{D}$ point of the $l$ layer, and find

$\displaystyle H_{inter-
vall}^{op}\approx\frac{\gamma_{3}}{\sqrt{2N_{G}{\sqrt{M(\alpha_{1}+\alpha_{2}+\alpha_{3})}}}}\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
MBZ}\sum_{l,\eta}\sum_{\mathbf{Q}_{l\eta},\mathbf{Q}^{\prime}_{-l\eta}}\sum_{\alpha\beta,s}$
$\displaystyle(b_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime}-\mathbf{Q}_{l\eta}+\mathbf{Q}^{\prime}_{-l\eta},l,A_{1}}+b_{\eta\mathbf{K}_{D}-\mathbf{k}+\mathbf{k}^{\prime}+\mathbf{Q}_{l\eta}-\mathbf{Q}^{\prime}_{-l\eta},l,A_{1}}^{\dagger})c^{\dagger}_{\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{Q}_{l\eta},l,\alpha
s}(\sigma_{x})_{\alpha\beta}c_{-\eta\mathbf{K}_{D}+\mathbf{k}^{\prime}-\mathbf{Q}^{\prime}_{-l\eta},l,\beta,s}$
$\displaystyle=\frac{\gamma_{3}}{\sqrt{2N_{G}{\sqrt{M(\alpha_{1}+\alpha_{2}+\alpha_{3})}}}}\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
MBZ}\sum_{l,\eta}\sum_{\mathbf{Q}_{l\eta},\mathbf{Q}^{\prime}_{-l\eta}}\sum_{\alpha\beta,s}$
$\displaystyle(b_{-\eta\mathbf{K}_{D}+\mathbf{k}-\mathbf{k}^{\prime}-\mathbf{Q}_{l\eta}+\mathbf{Q}^{\prime}_{-l\eta},l,A_{1}}+b_{\eta\mathbf{K}_{D}-\mathbf{k}+\mathbf{k}^{\prime}+\mathbf{Q}_{l\eta}-\mathbf{Q}^{\prime}_{-l\eta},l,A_{1}}^{\dagger})c^{\dagger}_{\mathbf{k},\mathbf{Q}_{l\eta},\alpha,\eta,s}(\sigma_{x})_{\alpha\beta}c_{\mathbf{k}^{\prime},\mathbf{Q}^{\prime}_{-l\eta},\beta,-\eta,s}$
(160)

We next project the above Hamiltonian into the flat bands of the BM model with
the expansion

$c_{\mathbf{k},\mathbf{Q}_{l\eta},\alpha,\eta,s}=\sum_{n=1,2}u^{n}_{\mathbf{k},\mathbf{Q}_{l\eta},\alpha,\eta}\gamma_{\mathbf{k}n\eta
s}$ (161)

where $u_{\mathbf{k},\mathbf{Q}_{l\eta},\alpha,\eta}$ is the eigenstate of the
BM model. This gives

$\displaystyle H_{inter-
vall}^{op}\approx\frac{\gamma_{3}}{\sqrt{2N_{G}{\sqrt{M(\alpha_{1}+\alpha_{2}+\alpha_{3})}}}}\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
MBZ}\sum_{n,n^{\prime}=1,2}\sum_{s}\gamma_{\mathbf{k},n,\eta,s}^{\dagger}\gamma_{\mathbf{k}^{\prime},n^{\prime},-\eta,s}$