# A many-body perturbation approach to moiré bands in twisted bilayer graphene

Federico Escudero<EMAIL_ADDRESS>IFISUR, Departamento de
Física, CONICET, Universidad Nacional del Sur, Av. Alem 1253, B8000 Bahía
Blanca, Argentina

(August 31, 2024)

###### Abstract

We develop a many-body perturbation theory to account for the emergence of
moiré bands in the continuum model of twisted bilayer graphene. Our framework
is build upon treating the moiré potential as a perturbation that transfers
electrons from one layer to another through the exchange of the three wave
vectors that define the moiré Brillouin zone. By working in the two-band basis
of each monolayer, we analyze the one-particle Green’s function and introduce
a diagrammatic representation for the scattering processes. We then identify
the moiré-induced self-energy, relate it to the quasiparticle weight and
velocity of the moiré bands, and show how it can be obtained by summing
irreducible diagrams. We also connect the emergence of flat bands to the
behavior of the static self-energy at the magic angle. In particular, we show
that a vanishing Dirac velocity is a direct consequence of the relative
orientation of the momentum transfer vectors, suggesting that the origin of
magic angles in twisted bilayer graphene is intrinsically connected to its
geometrical properties. Our approach provides a many-body diagrammatic
framework that highlights the physical properties of the moiré bands.

##  I Introduction

In recent years, moiré heterostructures have emerged as a promising platform
to study strongly correlated phases of matter [1, 2, 3, 4, 5]. The most
prominent example is twisted bilayer graphene (TBG) [6, 7, 8], where
superconductivity and metal-insulator transitions have been experimentally
observed around the magic angle ($\theta\sim 1^{\circ}$) [9, 10]. These
correlated effects are thought to be induced by a sharp increase in the
density of states as the lowest moiré bands become flat at the magic angle [9,
11]. The ability to easily manipulate the parameters of the system, such as
the twist angle or the carrier density, provides a rich experimental platform
to test and understand the nature of the correlated phases [12, 7, 6]. Such
understanding may prove to be not only relevant for the phenomena observed in
moiré heterostructures, but also for elucidating the mechanisms responsible
for unconventional superconductivity in more complex systems, such as in
cuprates [13, 14].

A first step towards modeling the electronic properties in moiré
heterostructures is describing its band structure within the independent
electron approximation, from which correlations can then be, when possible,
treated perturbatively. There are, in general, two methods to obtain the moiré
bands: first principle calculations [15, 16, 17, 18, 19], such as DFT or tight
binding computations, and effective continuum models [9, 11, 20, 21]. First
principle calculations, although naturally more exact, require commensurate
superlattice structures, which for twist-induced heterostructures it only
occur at a discrete set of twist angles [22, 23]. In addition, the number of
atoms within the supercells increases as the twist angle decreases. Around the
magic angle the supercell can typically contain up to ten thousands atoms
[24], thus requiring a high computational capacity and time to perform first-
principle calculation [25, 26]. This has naturally led to the need of having
effective continuum models that can capture the low-twist angle behavior in
the system. One of the main advantages of the continuum models is that they
can be constructed even if the moiré pattern is incommensurate [11, 27]. With
only a few set of parameters, such as the hopping energies between different
sublattices, the continuum models can yield results in excellent agreement
with those obtained by first-principle calculations [21, 28].

There are different continuum models in the literature [9, 11, 27, 29, 30,
31]. The most sophisticated and realistic ones take into account lattice
relaxation effects [28, 29, 25, 31, 32], usually through the inclusion of
deformation-induced gauge fields. The rich physical properties arising from
the continuum model, especially around the magic angles, have further sparked
numerous studies on its mathematical properties [33, 34, 35]. For twisted
bilayer graphene, all of these continuum models share the essence of coupling
both layers through a moiré-induced potential, whose strength is modulated by
the twist angle [11]. In momentum space, the leading order Fourier expansion
of the moiré potential couples the Dirac points in each layer through the
exchange of three twist-dependent wave vectors [27], which are related by a
hexagonal symmetry that determines the moiré Brillouine zone [11]. This leads
to one the main features of the continuum models, namely that its Hamiltonian
has no inherent momentum cutoff, since any Dirac point can be always coupled
to another one. In practice, a sufficiently high momentum cutoff is determined
by when the resulting moiré bands of interest (usually the low-energy ones)
converge to a steady value. As the twist angle decreases the strength of the
moiré coupling, relative to the energy difference between the uncopled Dirac
points, increases. Consequently, around the magic angle one generally needs a
relatively large momentum cutoff in order to archive a diagonalization
convergence for the low-energy moiré bands, which translates to a continuum
model Hamiltonian of high matrix dimensions. Besides being less computational
appealing, this increase in the need of a higher matrix dimension hinders the
understanding of how the flat moiré bands emerge and behave. Numerous lower
dimension models, which still aim to capture the continuum model physics, have
been consequently developed, such as two-band models [9, 11], helical network
models [36, 37, 38], ten-band models [28], Kondo-lattice models [39, 40],
$\left(2+2\right)$ models [41], and topological heavy-fermions models [42,
43].

In this work, we propose and develop an alternative many-body perturbative
approach to study the moiré bands. Our motivation is twofold: on one hand the
perturbative approach provides a natural framework to study the increasing
incidence of the moiré coupling on the one-particle excitations of the system
as the twist angle is decreased; on the other hand, deploying many-body
diagrammatic techniques gives a direct physical picture of the moiré-induced
scattering events taking place, and how this may lead to flat bands. Our
framework is build upon the usual continuum model, but explicitly working in
the two-band basis of each uncoupled layer. In this basis the moiré potential
acts as a one-body operator that transfers electrons from one layer to another
by the exchange of the three wave vectors that determine the moiré Brillouin
zone. By treating this potential as a perturbation to the uncoupled layers, we
obtain the moiré-induced self-energy that renormalizes the linear massless
dispersions onto the moiré bands, and discuss its physical interpretation by
identifying the contributions of irreducible diagrams. From the self-energy we
further obtain the one-particle spectral density, the quasiparticle weight and
the velocity of the moiré bands. We also relate the properties of the lowest
moiré bands to the behavior of the static self-energy at the magic angle, and
show that the emergence of flat bands can be directly connected to the
particular geometrical properties of the moiré potential in twisted bilayer
graphene.

The rest of this work is organized as follows: in Sec. II we give a brief
recap of the continuum model of twisted bilayer graphene, which serves as a
basis for the subsequent perturbation approach. Next we describe the
formulation of the continuum model in the two-band basis of each layer. In
Sec. III we develop the proposed many-body perturbation approach to the moiré
bands. We introduce the one-particle Green’s function and give its exact
closed form through a matrix inversion of the continuum model Hamiltonian in
the band basis. We identify the moiré-induced self energy, compute it to
leading order in perturbation, and give its physical interpretation. We then
discuss the one-particle spectral density, the velocity renormalization of the
moiré bands, and the emergence and origin of flat bands around the magic
angle. We end with the conclusions in Sec. IV.

##  II Continuum model of twisted bilayer graphene

Figure 1: Reciprocal space representation of two twisted honeycomb Brillouin
zones, with corresponding reciprocal vectors $\mathbf{b}_{1,2}^{\pm\theta/2}$
and Dirac points $\mathbf{K}_{\pm},\mathbf{K}_{\pm}^{\prime}$ for the top
$\left(+\right)$ and bottom $\left(-\right)$ layer, each rotated by
$\pm\theta/2$. The moiré vectors are given by the difference between the
rotated reciprocal vectors,
$\mathbf{g}_{1,2}=\mathbf{b}_{1,2}^{\theta/2}-\mathbf{b}_{1,2}^{-\theta/2}$.
In the continuum model, the twisted Dirac points in different layers are
coupled, to leading order, by the exchange of three momentum transfer vectors
$\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}$ that determine the borders of
the moiré Brillouin zone (right). For momenta originating from the $K$ valley
of both layers, the moiré valleys are set by $K_{M}=K_{+}$ and
$K^{\prime}_{M}=K_{-}$, while around the $K^{\prime}$ valley
$K_{M}=K^{\prime}_{-}$ and $K^{\prime}_{M}=K^{\prime}_{+}$. The three momentum
transfer vectors relate equivalent valleys within the moiré Brillouin zone.

Two honeycomb lattices in a bilayer configuration, relatively rotated by an
angle $\theta$, give rise to a moiré pattern of characteristic length $L\sim
a/2\sin\left(\theta/2\right)$ [6], where $a$ is the lattice constant of the
layers ($a\simeq 0.246\,\mathrm{nm}$ in graphene). For arbitrary twist angles
the moiré pattern is in general incommensurate, i.e., it is not technically a
crystal structure [22, 23]. Nevertheless, one can geometrically describe the
moiré pattern by defining a set of two moiré vectors $\mathbf{g}_{i}$
($i=1,2$), constructed by taking the minimum difference between the rotated
reciprocal vectors in each layer [21, 27],
$\mathbf{g}_{i}=\mathbf{b}_{i}^{\theta/2}-\mathbf{b}_{i}^{-\theta/2}$, where
$\mathbf{b}_{i}^{\pm\theta/2}=R\left(\pm\theta/2\right)\mathbf{b}_{i}$ with
$R$ the rotation matrix and $\mathbf{b}_{i}$ the reciprocal vectors of a
honeycomb lattice. In real space, the moiré vectors so defined connect _all_
the AA stacking regimes of the moiré pattern. We choose the primitive lattice
vectors of the layers as $\mathbf{a}_{1}=a\left(1,0\right)$ and
$\mathbf{a}_{2}=a\left(1/2,\sqrt{3}/2\right)$, which upon rotation lead to the
reciprocal space orientation shown in Fig. 1.

The continuum model captures the low-energy electronic properties of the moiré
pattern, in the regime where the moiré length is much larger than the atomic
length [9]. While first-principle models require commensurate superlattice
structures, the continuum model holds even if the moiré pattern is
incommensurate [11, 27]. There are different formulations of the continuum
model in the literature; they all share, however, the essence of coupling the
Dirac points in each layer through a moiré-induced potential. In this work we
will, in particular, consider the Bistritzer–MacDonald model [11], for momenta
relative to the $K$ valley in both layers (see Fig. 1). The corresponding
Hamiltonian reads

$H=\left[\begin{array}[]{cc}h_{+}\left(\mathbf{k}\right)&T^{\dagger}\left(\mathbf{r}\right)\\\
T\left(\mathbf{r}\right)&h_{-}\left(\mathbf{k}\right)\end{array}\right],$ (1)

where $h_{\pm}\left(\mathbf{k}\right)=\hbar
v\bm{\sigma}_{\pm\theta/2}\cdot\mathbf{k}_{\pm}$ are the Dirac Hamiltonian in
the top $\left(+\right)$ and bottom $\left(-\right)$ layers, with
$v=\sqrt{3}ta/2\hbar$ the Fermi velocity in graphene ($t$ being the nearest
sublattice hopping energy), $\mathbf{k}_{\pm}$ the momenta relative to the
rotated $K$ points, and
$\bm{\sigma}_{\pm\theta/2}=\left(\sigma_{x,\pm\theta/2},\sigma_{y,\pm\theta/2}\right)$
the rotated Pauli matrix vector [44, 45]. The angle dependence in the Pauli
matrices leads to a small energy shift in the system; in what follow we
neglect this and set $\bm{\sigma}_{\pm\theta/2}\rightarrow\bm{\sigma}$ (this
restores the particle-hole symmetry in the system) [11, 30].

The coupling between both layers is dictated by the moiré potential
$T\left(\mathbf{r}\right)$. At small twist angles, its leading order Fourier
expansion is given by [11, 21, 27]

$T^{\alpha\beta}\left(\mathbf{r}\right)=w_{\alpha\beta}\sum_{j=1}^{3}e^{-i\mathbf{q}_{j}\cdot\mathbf{r}}T_{j}^{\alpha\beta},$
(2)

where $\alpha$ and $\beta$ are the sublattice indices in each layer,
$w_{\alpha\beta}$ is the corresponding coupling strength, and $\mathbf{q}_{j}$
are the three equal-length wave vectors arising from the difference between
the three equivalent rotated Dirac points, cf. Fig. 1. The wave vectors
$\mathbf{q}_{j}$ thus determine the moiré Brillouin zone (mBZ), and can be can
be obtained from the moiré vectors as
$\mathbf{q}_{1}=\left(2\mathbf{g}_{1}+\mathbf{g}_{2}\right)/3$,
$\mathbf{q}_{2}=R\left(2\pi/3\right)\mathbf{q}_{1}$ and
$\mathbf{q}_{3}=R\left(4\pi/3\right)\mathbf{q}_{1}$. For an initial AA bilayer
configuration, upon rotation the coupling matrices for each wave vector
$\mathbf{q}_{j}$ read [27]

$T_{j}=\left(\begin{array}[]{cc}w_{0}&w_{1}e^{-i\phi_{j}}\\\
w_{1}e^{i\phi_{j}}&w_{0}\end{array}\right),$ (3)

where $w_{0}=w_{AA}$ and $w_{1}=w_{AB}$ are the AA and AB hopping energies,
and $\phi_{j}=\left(j-1\right)2\pi/3$. Due to relaxation effects, the ratio
$w_{0}/w_{1}$ tends to decrease at low twist angles, effectively reducing the
energetically costly AA stacking regimes [46, 28, 47]. In this work, unless
otherwise stated, we consider $w_{1}=0.11$ eV, $w_{0}=0.8w_{1}$ and $t=2.8$
eV, which gives the first magic angle at $\theta_{M}\simeq 1.06^{\circ}$ [11,
48].

The moiré bands are obtained by diagonalizing the continuum model Hamiltonian
in momentum space. As each Dirac point couples to other three by the exchange
of the wave vectors $\mathbf{q}_{j}$, the continuum model Hamiltonian has no
bound in momentum space. However, only the low-energy moiré bands are
relevant, for which the eigenstates effectively converge up to a sufficiently
large momentum cutoff. In this work we adopt the original Tripod model
expansion of the continuum model Hamiltonian [11, 30], in which the momentum
is measured from a Dirac point and the cutoffs are set by shells of
equidistant Dirac points from the origin, cf. Fig. 2.

The continuum model Hamiltonian around the $K^{\prime}$ valley can be obtained
from the one in the $K$ valley by replacing
$\bm{\sigma}\rightarrow\bm{\sigma}^{*}=\left(\sigma_{x},-\sigma_{y}\right)$
and $\phi_{j}\rightarrow-\phi_{j}$ [21]. The moiré bands for momenta around
the layer valleys $K$ and $K^{\prime}$ (not to be confused with the moiré
valleys $K_{M}$) only differ slightly around the $\Gamma_{M}$ and $M_{M}$
moiré symmetry points [21]; their general properties, particularly as a
function of the twist angle, are nevertheless the same. As intervalley
scattering events within each layer are negligible at low energies, throughout
this paper we continue to work with the TBG Hamiltonian for the $K$ valley,
keeping in mind that the results for the $K^{\prime}$ valley can be trivially
obtained by the above correspondence.

Figure 2: (a) Tripod model of the TBG Hamiltonian expansion in momentum space,
in which the momentum is measured from the Dirac point in one layer. Empty
(filled) circles correspond to Dirac points in the top (bottom) layer. Circles
indicate the truncation (or cutoffs) at increasing _shells_ , corresponding to
equidistant Dirac points. (b) Schematic representation of successive momentum
exchanges $\mathbf{k}_{-}-\mathbf{k}_{+}=\mathbf{q}_{j}$ between adjacent
Dirac points, up to the third shell.

###  II.1 Band basis

The continuum model Hamiltonian, Eq. (1), is written in the sublattice basis
for each layer. To make a connection with perturbation theory, it will prove
more convenient to rather express it in the band basis of each layer, in which
the Dirac Hamiltonian is diagonalized. The matrix elements of the moiré
potential in the band basis are given by (cf. Appendix A)

$\left\langle\mathbf{k}_{-},s^{\prime}\right|\hat{T}_{j}\left|\mathbf{k}_{+},s\right\rangle=\mathcal{T}_{j,\mathbf{k}_{+},s,s^{\prime}}\delta_{\mathbf{k}_{-},\mathbf{k}_{+}+\mathbf{q}_{j}},$
(4)

where $s,s^{\prime}=\pm 1$ are the band indices and

$\mathcal{T}_{j,\mathbf{k},s,s^{\prime}}=\psi_{\mathbf{k}+\mathbf{q}_{j},s^{\prime}}^{\dagger}T_{j}\psi_{\mathbf{k},s}.$
(5)

Here $T_{j}$ are the moiré matrices given by Eq. (3) and

$\psi_{\mathbf{k},s}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
se^{i\theta_{\mathbf{k}}}\end{array}\right)$ (6)

are the pseudospinors eigenstates in graphene (accounting for the neglected
rotation of the Pauli matrices amounts to replace
$\theta_{\mathbf{k}}\rightarrow\theta_{\mathbf{k}}\pm\theta/2$ in
$\psi_{\mathbf{k},s}$). The overlap of these pseudospinors leads to the
momentum dependence of the moiré potential in the band basis. Specifically,
the AA hoppings lead to a contribution of the form $\mathcal{T}_{j}\sim
w_{0}\psi_{\mathbf{k}+\mathbf{q}_{j},s^{\prime}}^{\dagger}\psi_{\mathbf{k},s}$,
whereas the AB/BA hoppings lead to a contribution of the form
$\mathcal{T}_{j}\sim
w_{1}\psi_{\mathbf{k}+\mathbf{q}_{j},s^{\prime}}^{\dagger}\bm{\sigma}\cdot\bm{\phi}_{j}\psi_{\mathbf{k},s}$,
where $\bm{\phi}_{j}$ is the versor with an angle $\phi_{j}$ relative to the
$x$ axis. The latter overlap is similar to the one obtained for the
interaction between two chiral states mediated by a vector potential
$\mathbf{A}_{j}\sim w_{1}\bm{\phi}_{j}$.

The use of the band basis provides a natural framework from which one can
deploy perturbation methods into the problem. To this end, we first second
quantize the continuum model Hamiltonian and separate
$\hat{H}=\hat{H}_{0}+\hat{T}$, where

$\displaystyle\hat{H}_{0}$
$\displaystyle=\sum_{\mathbf{k},s}\sum_{\ell=\pm}\epsilon_{\mathbf{k},s}\hat{c}_{\ell,\mathbf{k},s}^{\dagger}\hat{c}_{\ell,\mathbf{k},s},$
(7) $\displaystyle\hat{T}$
$\displaystyle=\sum_{\mathbf{k},s,s^{\prime}}\sum_{j=1}^{3}\mathcal{T}_{j,\mathbf{k},s,s^{\prime}}\hat{c}_{-,\mathbf{k}+\mathbf{q}_{j},s^{\prime}}^{\dagger}\hat{c}_{+,\mathbf{k},s}+h.c.$
(8)

Here $\ell=\pm$ runs over the top and bottom layers. The operators
$\hat{c}_{\ell,\mathbf{k},s}^{\dagger}$ and $\hat{c}_{\ell,\mathbf{k},s}$
create and annihilate an electron in the $\ell$ layer, with momentum
$\mathbf{k}$ in the $s$ band and energy $\epsilon_{\mathbf{k},s}=s\hbar
v\left|\mathbf{k}\right|$. We adopt the convention that the momentum is
measured relative to a Dirac point in the top layer. The moiré potential
$\hat{T}$ can be thought as a one-body operator that exchanges electrons from
one layer to another by the transfer of a momentum $\mathbf{q}_{j}$, cf. Fig.
3.

Since the potential $\hat{T}$ is quadratic in the operators, the total
Hamiltonian can be, of course, exactly diagonalized:

$\hat{H}=\sum_{\mathbf{k}\in\mathrm{mBZ}}\sum_{\nu}E_{\mathbf{k},\nu}\hat{a}_{\mathbf{k},\nu}^{\dagger}\hat{a}_{\mathbf{k},\nu}.$
(9)

Here the operators $\hat{a}_{\mathbf{k},\nu}^{\dagger}$ and
$\hat{a}_{\mathbf{k},\nu}$ create and annihilate an electron in the $\nu$
moiré band with momentum $\mathbf{k}$ within the moiré BZ, and energy
$E_{\mathbf{k},\nu}$. Eq. (9) formally solves the problem, provided that the
energies $E_{\mathbf{k},\nu}$ are obtained by diagonalizing the continuum
model Hamiltonian. However, such matrix diagonalization implicitly requires a
truncation of the TBG Hamiltonian, on which the number of moiré bands in the
sum over $\nu$, as well as the corresponding expressions for their energies,
depend. The same holds for the moiré operators $\hat{a}_{\mathbf{k},\nu}$,
whose expansion in terms of the $\hat{c}$ operators generally depends on the
momentum cutoff. This sort of perturbation dependence in the continuum model
is hidden in Eq. (9). The effectiveness of the truncation, as measured by
diagonalization convergence, further depends on the twist angle: whereas at
relatively high angles it may be sufficient to consider a low momentum cutoff
to describe the lowest moiré bands, at lower twist angles one generally
requires a higher momentum cutoff for convergence. As a result, a physical
understanding of the emergent properties of the moiré bands becomes
increasingly harder to get as the twist angle decrease, since a larger
momentum cutoff translates to a TBG Hamiltonian of high dimensions.

##  III Many-body perturbation approach to moiré bands

In this work, rather than reducing the problem to a matrix diagonalization, we
propose and develop a many-body perturbation approach to study the nature of
the moiré bands. We treat the system as an effective many-body problem in
which the moiré potential $\hat{T}$ acts as a perturbation to the Hamiltonian
$\hat{H}_{0}$. Our interest will be the incidence of the perturbation on the
one-particle excitations in each layer. In particular, we consider the
imaginary time one-particle Green’s function [49]

$\mathcal{G}_{\gamma\gamma^{\prime}}\left(\tau\right)=-\left\langle
T_{\tau}\hat{c}_{\gamma}\left(\tau\right)\hat{c}_{\gamma^{\prime}}^{\dagger}\left(0\right)\right\rangle,$
(10)

where the index $\gamma=\left(\ell,\mathbf{k},s\right)$ accounts for the
layer, momentum and band indices of the $\hat{c}$ operators, $T_{\tau}$ is the
time-ordering operator, and $\left\langle\cdots\right\rangle$ is the ensemble
average over the interacting system. By treating the moiré potential $\hat{T}$
as a perturbation, the ensemble average $\left\langle\cdots\right\rangle$ can
be expanded as an ensemble average over the noninteracting system, cf.
Appendix B.

Since the problem is exactly solvable, the Green’s function can be obtained in
closed form. Indeed, in frequency-momentum space the Green’s function of the
system can be reduced to a matrix inversion [49]
$\mathcal{G}\left(i\omega_{n}\right)=\left(i\omega_{n}\mathbb{I}-H_{B}\right)^{-1}$,
where $\omega_{n}=\left(2n+1\right)\pi/\hbar\beta$ are the fermionic Matsubara
frequencies (with $n$ an integer and $\beta=1/k_{B}T$), $H_{B}$ is the
continuum model Hamiltonian in the band basis (cf. Appendix A), and
$\mathbb{I}$ is the identity matrix of the same dimension as $H_{B}$. The
Green’s function given by Eq. (10) can then be obtained as

$\mathcal{G}_{\gamma\gamma^{\prime}}\left(i\omega_{n}\right)=\left(i\omega_{n}\mathbb{I}-H_{B}\right)_{\gamma\gamma^{\prime}}^{-1},$
(11)

with $\gamma$ and $\gamma^{\prime}$ being matrix components such that
$\left(H_{B}\right)_{\gamma\gamma}=\epsilon_{\gamma}$ and
$\left(H_{B}\right)_{\gamma^{\prime}\gamma^{\prime}}=\epsilon_{\gamma^{\prime}}$.
The retarded Green’s function follows by analytical continuation
$i\omega_{n}\rightarrow\omega+i0^{+}$.

From a numerical point of view, Eq. (11) gives a straightforward method to
obtain the Green’s function of the system. However, as it stands, it does not
provide a direct physical picture of the moiré-induced scattering processes
taking place. To this end, in what follows we make a connection between Eq.
(11) and the usual many-body diagrammatic techniques. We identify the self-
energy, discuss its diagrammatic expansion, and relate it to the spectral
density, the quasiparticles velocity, and the emergence of flat bands around
the magic angle.

Figure 3: (a) Diagrammatic representation for the scattering between electrons
in different layers induced by the moiré potential. Solid lines represent the
_free_ Green’s function
$\mathcal{G}_{\mathbf{k},s}^{0}\left(i\omega\right)=1/\left(i\omega_{n}-\epsilon_{\mathbf{k},s}\right)$
for a Dirac electron with momentum $\mathbf{k}$ in the $s=\pm 1$ band and
energy $\epsilon_{\mathbf{k},s}=s\hbar v\left|\mathbf{k}\right|$. Empty
(filled) circles represent the moiré potential strength $\mathcal{T}$
$\left(\mathcal{T}^{*}\right)$ for the scattering of an electron in the top
(bottom) layer to the other layer. (b) A curved line connecting one empty and
one filled circle is used to indicate the same exchanged momentum
$\mathbf{q}_{j}$, regardless of possible in-between scattering events that may
modify the initial and final momenta.

###  III.1 Self-energy

Since the moiré potential, Eq. (8), exchanges electrons from one layer to
another, with different momenta and arbitrary band indices, in general the
Green’s function given by Eq. (11) will depend on the indices
$\gamma,\gamma^{\prime}$. Nevertheless, most of the information about the one-
particle excitations spectra can be deduced by analyzing the Green’s functions
with $\gamma=\gamma^{\prime}$. To simplify the notation we shall write
$\mathcal{G}_{\gamma\gamma}\equiv\mathcal{G}_{\gamma}$. The Dyson equation in
momentum-frequency space then leads to the exact expression [50] (from now on
we set $\hbar=1$)

$\mathcal{G}_{\gamma}\left(i\omega_{n}\right)=\frac{1}{i\omega_{n}-\bar{\epsilon}_{\gamma}-\Sigma_{\gamma}\left(i\omega_{n}\right)},$
(12)

where $\bar{\epsilon}_{\gamma}=\epsilon_{\gamma}-\mu$ are the unperturbed
energies measured relative to the chemical potential $\mu$ (henceforth we set,
for simplicity, the energy origin at $\mu=0$, and drop the bar notation), and
$\Sigma_{\gamma}\left(i\omega_{n}\right)$ is the proper self-energy
(henceforth just the _self-energy_). As with the Green’s function, the self-
energy can be directly obtained from the continuum model Hamiltonian through
Eq. (11), which implies that

$\Sigma_{\gamma}\left(i\omega_{n}\right)=i\omega_{n}-\epsilon_{\gamma}-\frac{1}{\left(i\omega_{n}\mathbb{I}-H_{B}\right)_{\gamma\gamma}^{-1}},$
(13)

where the component index $\gamma$ should be read as the index at which
$H_{\gamma\gamma}=\epsilon_{\gamma}$.

Figure 4: Top: Diagrammatic representation of the moiré-induced self-energy;
double lines indicate the full dressed Green’s function. Bottom: leading order
examples of irreducible diagrams that contribute to the self-energy. First row
is in terms of the Feynman diagrams detailed in Fig. 3. Internal momentum
transfers and band indices are to be summed over. Diagrams in which at a
midpoint the electron returns to the same initial momentum, as in the second
diagram, imply a change of band. Last row is a hexagonal picture
representation of the above diagrams, in the same order. Curved lines indicate
the path of an electron by coupling to successive Dirac points. For simplicity
only one possible equivalent path is shown.

The self-energy describes the renormalization of the uncouple dispersions
$\epsilon_{\gamma}$ due to the interaction with the moiré potential. One can
thus obtain valuable insights into the nature of the moiré bands by solely
analyzing the properties of the self-energy. One of the main advantages of the
developed approach is that the self-energy can be obtained diagrammatically by
summing irreducible diagrams. Some examples of those are shown in Fig. 4. The
diagrammatic approach gives not only a straightforward method to compute the
self-energy (in principle, to any perturbation order), but also a clear
physical picture of the involved moiré-induced scattering processes. This can
be very valuable for understanding the emergent properties of the moiré bands,
particularly when compared to the numerical diagonalization procedure. The
diagrammatic representation also allows one to deploy well-established many-
body perturbation techniques, which may be particularly useful for finding
effective theories for the moiré bands of interest.

It is worth noting that the perturbation description put forward here _does
not_ imply small corrections to the unperturbed energies, since those can
actually be strongly renormalized, especially at low twist angles. Indeed, the
diagonalization convergence requires an increasingly higher perturbation order
as the twist angle decreases (i.e., higher momentum cutoff in the continuum
model Hamiltonian), which naturally leads to a stronger energy
renormalization. As long as the dressed Green’s function converges (so that
Eq. (12) holds), such an increase in renormalization does not, however,
invalidate the perturbation approach.

####  III.1.1 Physical interpretation

The self-energy can be directly computed by summing irreducible diagrams, cf.
Figs. 3 and 4. Due to the strong coupling of the moiré potential at low twist
angles, and the nonnegligible contribution of interband transitions, one
generally needs to consider many irreducible diagrams in order to achieve a
convergence of the self-energy. A valuable guide to visualize the relevant
scattering events can be obtained by giving a diagrammatic representation to
Eq. (13). For conciseness, throughout this section we will focus on the self-
energy for the top layer, so that
$\gamma=\text{$\left(+,\mathbf{k},s\right)$}$. The results are readily
generalized for the bottom layer upon the replacement
$\mathbf{k}\rightarrow\mathbf{k}+\mathbf{q}_{1}$,
$\mathbf{q}_{j}\rightarrow-\mathbf{q}_{j}$ and
$\mathcal{T}_{j}\rightarrow\mathcal{T}_{j}^{*}$ (see Fig. 1).

We start by considering the lowest TBG Hamiltonian in the Tripod model, where
the momentum is truncated at the first shell (cf. Fig. 2). This leads to the
well known $8\times 8$ truncated TBG Hamiltonian, which in the band basis
takes the form

$H_{B}^{\left(1\right)}=\left(\begin{array}[]{cccc}h_{\mathbf{k}}^{B}&\mathcal{T}_{1,\mathbf{k}}^{\dagger}&\mathcal{T}_{2,\mathbf{k}}^{\dagger}&\mathcal{T}_{3,\mathbf{k}}^{\dagger}\\\
\mathcal{T}_{1,\mathbf{k}}&h_{\mathbf{k}+\mathbf{q}_{1}}^{B}&0&0\\\
\mathcal{T}_{2,\mathbf{k}}&0&h_{\mathbf{k}+\mathbf{q}_{2}}^{B}&0\\\
\mathcal{T}_{3,\mathbf{k}}&0&0&h_{\mathbf{k}+\mathbf{q}_{3}}^{B}\end{array}\right),$
(14)

where $h_{\mathbf{k}}^{B}=-v\left|\mathbf{k}\right|\sigma_{z}$ is the
diagonalized Dirac Hamiltonian and [cf. Eq. (5)]

$\mathcal{T}_{j,\mathbf{k}}=\left(\begin{array}[]{cc}\mathcal{T}_{j,\mathbf{k},-,-}&\mathcal{T}_{j,\mathbf{k},+,-}\\\
\mathcal{T}_{j,\mathbf{k},-,+}&\mathcal{T}_{j,\mathbf{k},+,+}\end{array}\right).$
(15)

The self-energy given by Eq. (13), for the top layer, can then be recasted as

$\Sigma_{+,\mathbf{k},s}^{\left(1\right)}\left(i\omega_{n}\right)=\Sigma_{\mathbf{k},s}^{\mathrm{intra}}\left(i\omega_{n}\right)+\Sigma_{\mathbf{k},s}^{\mathrm{inter}}\left(i\omega_{n}\right),$
(16)

where

$\displaystyle\Sigma_{\mathbf{k},s}^{\mathrm{intra}}\left(i\omega_{n}\right)$
$\displaystyle=\sum_{j,r}\frac{\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}},$
(17)
$\displaystyle\Sigma_{\mathbf{k},s}^{\mathrm{inter}}\left(i\omega_{n}\right)$
$\displaystyle=\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},-s,r}\mathcal{T}_{j,\mathbf{k},s,r}^{*}}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}$
$\displaystyle\quad\times\left(i\omega_{n}-\epsilon_{\mathbf{k},-s}-\sum_{j,r}\frac{\left|\mathcal{T}_{j,\mathbf{k},-s,r}\right|^{2}}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)^{-1}.$
(18)

To give a diagrammatic interpretation of these expressions, we first recall
the Feynman rules of the system shown in Fig. 3. From it one readily sees that
the two irreducible contributions to $\Sigma^{\left(1\right)}$ have a simple
physical interpretation:

  * •

The term $\Sigma_{\mathbf{k},s}^{\mathrm{intra}}$ accounts for intraband
transitions in which an electron in the state
$\left|\mathbf{k},s\right\rangle$ scatters back and forth to the opposite
layer only once, returning to the same initial state
$\left|\mathbf{k},s\right\rangle$.

  * •

The term $\Sigma_{\mathbf{k},s}^{\mathrm{inter}}$ accounts for all interband
transitions in which an electron in the state
$\left|\mathbf{k},s\right\rangle$ transitions to the opposite band state
$\left|\mathbf{k},-s\right\rangle$ by scattering back and forth to the
opposite layer, before finally scattering again back to the initial state
$\left|\mathbf{k},s\right\rangle$.

The corresponding diagrammatic representation is shown in Fig. 5. The full
perturbation calculation that effectively leads to Eq. (16) to leading order
can be found in the Appendix B.

Figure 5: Diagrammatic representation of the moiré-induced self-energy
$\Sigma_{+,\mathbf{k},s}^{\left(1\right)}=\Sigma_{\mathbf{k},s}^{\mathrm{intra}}+\Sigma_{\mathbf{k},s}^{\mathrm{inter}}$
corresponding to the first shell truncation of the continuum model
Hamiltonian. The first diagram $\Sigma_{\mathbf{k},s}^{\mathrm{intra}}$ gives
the intraband contribution, where an electron scatters back and forth to the
opposite layer, returning to the initial state in the same band. The second
diagram $\Sigma_{\mathbf{k},s}^{\mathrm{inter}}$ gives the interband
contribution, where an electron transitions to the opposite band by scattering
back and forth to the opposite layer, before scattering back again to the
initial state.

The leading order diagrammatic representation of Eq. (13) suggest that, in
general, the self-energy for a shell truncation of the continuum model
Hamiltonian is given by the sum of all irreducible scattering processes with
momentum transfer up to the corresponding cutoff. It is important to note that
the perturbation order in the computation of $\Sigma$ (cf. Appendix B) does
not generally equal number of shells in the truncated Hamiltonian. This is due
to the incidence of interband scattering events, which imply higher
perturbation orders than intraband ones.

Consider the next truncation of the TBG Hamiltonian at the second shell. The
corresponding self-energy $\Sigma{}^{\left(2\right)}$ can be obtained directly
from Eq. (13) by inverting the $20\times 20$ truncated Hamiltonian. The
resulting expression, however, is highly cumbersome and hard to work with,
hindering the inference of physical insight from it. Instead, one can obtain a
relatively simple expression for $\Sigma^{\left(2\right)}$ by using the
diagrammatic approach, taking into account all the possible irreducible
scattering events up to second order in momentum transfer. The full
calculation is large and left to the Appendix C (see Fig. 11 for the
diagrammatic representation).

As the number of shells increase, the exact diagrammatic computation of the
self-energy get increasingly harder very rapidly. The situation actually
becomes quite complicated beyond the second shell, since then one has to take
into account multiple connected paths that an electron can take (see, for
instance, the last diagram of $\Sigma$ in Fig. 4). This complexity in the
computation of the self-energy, by diagrammatic methods, can be in a way
related to the nonlinear matrix dimension increase of the continuum model
Hamiltonian as the momentum cutoff increases.

###  III.2 Spectral density

The one-particle spectral density is given by
$A_{\gamma}=-\frac{1}{\pi}\mathrm{Im}\mathcal{G}_{\gamma}^{R}$, where
$\mathcal{G}_{\gamma}^{R}=\mathcal{G}_{\gamma}\left(i\omega_{n}\rightarrow\omega+i0^{+}\right)$
is the retarded Green’s function (cf. Eq. (11); we continue to consider the
case $\gamma^{\prime}=\gamma$). The renormalized quasiparticle energies follow
from the peaks in the spectral density, which occur when [51]

$\omega-\epsilon_{\gamma}-\Sigma_{\gamma}^{\prime}\left(\omega\right)=0,$ (19)

where
$\Sigma^{\prime}\left(\omega\right)=\mathrm{Re}\left[\Sigma\left(i\omega_{n}\rightarrow\omega+i0^{+}\right)\right]$
is the real part of the retarded self-energy. The quasiparticle energies
$\omega$ that satisfy the above equation are the same eigenvalues that
diagonalize the continuum model Hamiltonian, i.e., the moiré bands. However,
each state carries different quasiparticle weight [51, 52]. To see this, we
first note that since each quasiparticle state diagonalizes the TBG
Hamiltonian, in perturbation terms their lifetime is _infinite_. Consequently,
the imaginary part of the self-energy _vanishes_ when
$\omega=E_{\mathbf{k},\nu}$ [this can be proved rigorously from Eq. (11)].
Thus
$A_{\gamma}=\delta\left[\omega-\epsilon_{\gamma}-\Sigma_{\gamma}^{\prime}\left(\omega\right)\right]$,
and resolving the delta we obtain

$A_{\gamma}=\sum_{\nu}Z_{\nu,\gamma}\delta\left(\omega-
E_{\mathbf{k},\nu}\right),$ (20)

where

$Z_{\nu,\gamma}^{-1}=1-\left(\frac{\partial\Sigma_{\gamma}^{\prime}}{\partial\omega}\right)_{\omega=E_{\mathbf{k},\nu}}$
(21)

is the quasiparticle weight of the moiré band with energy
$E_{\mathbf{k},\nu}$. It can be proved that the normalization condition
$\sum_{\nu}Z_{\nu,\gamma}=1$ is satisfied. The quasiparticle weight
$Z_{\nu,\gamma}$ _measures_ the amplitude of creating a single-particle state
with energy $E_{\mathbf{k},\nu}$ from the interaction of the bare electrons
with the moiré potential:

$Z_{\nu,\gamma}\sim\left|\left\langle\mathbf{k},\nu\right|\hat{c}_{\gamma}^{\dagger}\left|\Omega\right\rangle\right|^{2},$
(22)

where $\left|\Omega\right\rangle$ is the total (interacting) ground state and
$\left|\mathbf{k},\nu\right\rangle$ is the excited single-particle state with
energy $E_{\mathbf{k},\nu}$.

Figure 6: Quasiparticle weight (QW) of the moiré bands, in the top and bottom
layers, for twist angles $\theta=5^{\circ}$ and $\theta=1.06^{\circ}$ (magic
angle). The QWs are calculated with Eq. (21) summing the contribution of
electrons in both Dirac bands. At low twist angles the QW is mainly on the
lowest moiré bands, while at lower twist angles the QW disperses to remote
moiré bands.

Fig. 6 shows the quasiparticle weight (QW) in the top and bottom layers,
calculated by summing Eq. (21) over both Dirac bands in each case. The
strongest QW is always on the lowest (in energy) moiré bands, and around the
$K_{M}$ ($K^{\prime}_{M}$) for the top (bottom) layer, consistent with Fig. 1.
At relatively high twist angles only low-momentum couplings contribute to the
self-energy, so the QW is almost entirely contained in those lowest moiré
bands. This means that the quasiparticle renormalization around the Dirac
point in a layer is only influenced by the coupling to the nearest Dirac
points in the other layer. At lower twist angles the QW disperses to higher
moiré bands, resulting in a significant reduction in the weight of the lowest
moiré bands. Around the magic angle ($\theta_{M}\sim 1.06^{\circ}$) the QW of
the lowest moiré bands are effectively reduced by more than a half, with
higher bands having an increasingly higher weight. In addition, the QW of the
lowest moiré bands in both layers becomes more uniformly distributed in
momentum space, reflecting the high hybridization of the Dirac points.
Breaking the layer degeneracy, for instance by an effective mass in one layer
arising from a hBN substrate effect [53], can thus result in a larger split of
the narrow moiré bands as the twist angle decreases [54].

Figure 7: Difference between the QW contribution of each Dirac band in the top
layer, for the low-energy moiré bands around the magic angle. Color red (blue)
indicate a QW coming mainly from electrons in the top (bottom) Dirac band. An
inversion behavior in the lowest moiré bands takes place around the magic
angle $\theta=1.06^{\circ}$.

By differentiating in each layer the QW contribution of electrons in each
Dirac band, one observes that an interesting band-inversion behavior takes
place around the magic angle [55]. This is seen in Fig. 7, where we show the
QW difference $\Delta Z=Z_{\nu,\mathbf{k},+}-Z_{\nu,\mathbf{k},-}$ in the top
layer, for the low-energy moiré bands. For $\theta>\theta_{M}$ the QW of the
top (bottom) lowest moiré band comes mainly from electrons in the top (bottom)
Dirac band. This behavior is expected for a low moiré coupling, or when
interband transitions are nondominant. However, as the moiré bands flatten
around the magic angle, a band-inversion behavior takes place, with the QW of
the top lowest moiré band now coming mainly from electrons in the bottom Dirac
band, and vice versa for the other lowest moiré band. As the twist angle is
further decreased, one sees a non-uniform, momentum-dependent QW behavior in
the lowest moiré bands, where the weight of electrons in the top Dirac band
may dominate around the $\Gamma_{M}$ point, while that of electrons in the
bottom Dirac band may dominate around the $K_{M}$ point (in the particular
case of $\theta=0.95^{\circ}$ in Fig. 7, such behavior arises because the
remote bands touch the lowest bands at the $\Gamma_{M}$ point). In general, a
band-inversion takes place when two moiré bands touch at a particular momentum
$\mathbf{k}$, or when at such point the quasiparticle velocity vanish. This
band-inversion behavior has been noted before by projecting the moiré
eigenstates on the eigenstates of the decoupled layers [55].

The quasiparticle weight in each layer can be experimentally relevant in
angle-resolved photoemission spectroscopy (ARPES) of moiré heterostructures
[56, 57, 58, 59], particularly because then, by design, the topmost layer
typically contributes the most to the signal. Within the simplest sudden
approximation, the resolved photoemission intensity $I$ can be directly linked
to the one-particle spectral density [60, 52]

$I\sim\mathcal{M}A_{\mathbf{k}}\left(\omega\right)f\left(\omega\right),$ (23)

where $f\left(\omega\right)$ is the Fermi-Dirac distribution and $\mathcal{M}$
are the matrix elements accounting for the signal dependence on the incident
light intensity, polarization and frequency (see, e.g., Refs. [61, 62] for a
detailed account). The above equation gives a direct intensity dependence on
the spectral density and thus a signal pattern similar to that in Fig. 6,
which is in relatively well agreement with recent experiments [63, 64, 65]. We
note that photoemission spectroscopy measurements have been previously
explained by a band-unfolding approach [57, 63], which in essence is similar
to what Eq. (22) measures.

###  III.3 Renormalized velocity and flat bands

The emergence of flat moiré bands is one of the most remarkable and intriguing
aspects of the continuum model, prompting many works devoted to study their
nature and origin [48, 30, 66, 32]. Within the perturbation scheme developed
here, the moiré-induced self-energy provides a direct account for the
quasiparticle velocity renormalization and the origin of moiré flat bands. In
line with the previous sections, in what follows we will consider the
perturbation on electrons in the top Dirac band of the top layer, and for
brevity omit those indices in the self-energy. Thus throughout this section it
should be understood that
$\Sigma{}_{+,\mathbf{k},s}\rightarrow\Sigma{}_{\mathbf{k}}$. The analysis
done, however, is valid for both layers and Dirac bands.

From Eq. (19) it follows that the quasiparticle velocity
$\mathbf{v}_{\mathbf{k},\nu}^{\star}=\nabla_{\mathbf{k}}E_{\mathbf{k},\nu}$
vector of the moiré band $\nu$ is given by

$\mathbf{v}_{\mathbf{k},\nu}^{\star}=Z_{\mathbf{k},\nu}\left[v\hat{\mathbf{k}}+\left(\nabla_{\mathbf{k}}\Sigma^{\prime}_{\mathbf{k}}\right)_{\omega=E_{\mathbf{k},\nu}}\right],$
(24)

where $Z_{\mathbf{k},\nu}$ is the quasiparticle weight given by Eq. (21), and
$\nabla_{\mathbf{k}}\epsilon_{\mathbf{k}}=v\hat{\mathbf{k}}$ is the velocity
vector in graphene. We are, in particular, interested in the renormalized
velocity $\mathbf{v}^{\star}$ of the lowest (in energy) moiré bands at the
Dirac point. To leading order in the perturbation, Eq. (16) gives
$\Sigma^{\prime}_{\mathbf{k}=0}\left(\omega=0\right)=0$ and
$\nabla_{\mathbf{k}}\Sigma^{\prime}_{\mathbf{k}}\left(k=0=\omega\right)\propto\hat{\mathbf{k}}$,
and it can be numerically checked from Eq. (13) that the same holds at higher
orders. Consequently, $\mathbf{v}^{\star}=v^{\star}\hat{\mathbf{k}}$ with

$\frac{v^{\star}}{v}=\frac{1+\left(\partial\Sigma^{\prime}_{\mathbf{k}}/\partial
k\right)_{0}/v}{1-\left(\partial\Sigma^{\prime}_{\mathbf{k}}/\partial\omega\right)_{0}},$
(25)

where the subindex “0” implies evaluating the derivatives at $k=0=\omega$. To
leading order, Eq. (16) gives
$\left(\partial\Sigma^{\prime}_{\mathbf{k}}/\partial
k\right)_{0}=-3v\alpha_{1}^{2}$ and
$\left(\partial\Sigma^{\prime}_{\mathbf{k}}/\partial\omega\right)_{0}=-3\left(\alpha_{0}^{2}+\alpha_{1}^{2}\right)$,
where $\alpha_{i}=w_{i}/vk_{\theta}$ are the moiré coupling strengths.
Replacing in Eq. (25) then leads to the well-known result
$v^{\star}/v=\left(1-3\alpha_{1}^{2}\right)/\left(1+3\alpha_{0}^{2}+3\alpha_{1}^{2}\right)$
when the TBG Hamiltonian is truncated at the first shell [11, 48, 30]. The
computation of $v^{\star}$, using the self-energy, can be carried out at
higher perturbation orders, cf. Appendix D.

Figure 8: (a) Static self-energy at the magic angle $\theta=1.06^{\circ}$ for
$w_{0}=0$ (chiral model) and $w_{0}=0.8w_{1}$. (b) Density plots for the
uncoupled dispersion $-\epsilon_{\mathbf{k}}/vk_{\theta}$ (top), and the
static self-energy $\Sigma_{\mathbf{k}}\left(0\right)/vk_{\theta}$ at the
magic angle (bottom), centered at the Dirac point $K_{M}$. Except for small
deviations around the $\Gamma_{M}$ point, the static self-energy converges to
$\Sigma_{\mathbf{k}}\left(0\right)\sim-vk$, leading to zero-energy flat bands.

When $v^{\star}$ vanishes the lowest moiré bands actually flatten across the
whole moiré BZ [11, 30]. To understand this overall behavior one needs to go
beyond the analysis of the quasiparticle velocity at the Dirac point. The
self-energy again provides a natural framework for this. Indeed, from Eq. (19)
it follows that low-energy ($\omega\sim 0$) flat bands emerge when the
condition

$\epsilon_{\mathbf{k}}+\Sigma_{\mathbf{k}}\left(0\right)\sim 0$ (26)

is satisfied, which can only occur when the static self-energy behaves as
$\Sigma\propto-k$ to leading order in $\mathbf{k}$.

Considering the self-energy given by Eq. (16), expanding in orders of
$k/k_{\theta}$ (where $k_{\theta}=\left|\mathbf{q}_{j}\right|$), at $\omega=0$
one has

$\displaystyle\Sigma_{\mathbf{k}}^{\left(1\right)}\left(0\right)$
$\displaystyle\simeq-3vk\left[\alpha_{1}^{2}-\alpha_{0}^{2}\sin\left(3\theta_{\mathbf{k}}\right)\left(k/k_{\theta}\right)\right.$
$\displaystyle\quad\left.-2\alpha_{0}\alpha_{1}\cos\left(3\theta_{\mathbf{k}}\right)\left(k/k_{\theta}\right)^{2}+\mathcal{O}\left(k/k_{\theta}\right)^{3}\right].$
(27)

Up to order $\sim\left(k/k_{\theta}\right)$, the flat band condition
$\epsilon_{\mathbf{k}}+\Sigma_{\mathbf{k}}^{\left(1\right)}\left(0\right)\sim
0$ is then satisfied at the first magic angle where $\alpha_{1}^{2}\sim 1/3$
and $v^{\star}\sim 0$. Since the two leading order corrections in Eq. (27)
scale with $\sim\alpha_{0}$, the moiré bands become more flat as the AA
hopping energies decrease.

Eq. (27) only gives a partial picture at the simplest first shell
approximation. At higher perturbation orders (i.e., higher number of shells)
one actually sees that at exactly the magic angle the self-energy converges to
$\Sigma_{\mathbf{k}}\left(0\right)\sim-vk$, only up to small deviations around
the $\Gamma_{M}$ point, see Fig. 8. As a result, the uncoupled dispersion
$\epsilon_{\mathbf{k}}=vk$ is renormalized to zero-energy moiré bands that
satisfy the flat band condition
$\epsilon_{\mathbf{k}}+\Sigma_{\mathbf{k}}\left(0\right)\sim 0$ over almost
the whole moiré Brillouin zone. In the chiral limit one has
$\Sigma_{\mathbf{k}}\left(0\right)\sim-vk$ in _all_ the mBZ, leading to the
formation of _absolute_ flat bands [48].

###  III.4 Origin of flat bands

Figure 9: Renormalized Dirac velocity at the first magic angle in the chiral
limit, as a function of the orientation $\varphi$ of the momentum transfer
vectors $\mathbf{q}_{j}$ relative to the phases $\phi_{j}$ in the moiré
potential (see inset). The solid line is the converged numerical result, while
the dashed line is the analytical result at the first shell given by Eq. (30).
The velocity only vanishes in the TBG case $\varphi=90^{\circ}$.

We now ask what properties of the moiré potential lead to those in the self
energy that, in particular, allow the flat band condition (26) to be
satisfied. A first point to note is that the behavior of the moiré potential,
as given by Eq. (5), depends on both the phases $\phi_{j}$ of the moiré
matrices and the momentum transfer vectors $\mathbf{q}_{j}$. In the continuum
model the phases $\phi_{j}$ arise from the underlying hexagonal symmetry of
the layers (technically, from phase factors $\sim
e^{i\mathbf{b}\cdot\bm{\delta}}$ relating the reciprocal lattice vectors
$\mathbf{b}$ and the basis vectors $\bm{\delta}$ of a honeycomb lattice [9,
11, 21]), and are therefore _independent_ of the geometrical properties of the
moiré pattern. The momentum transfer vectors $\mathbf{q}_{j}$, on the other
hand, depend on the orientation and shape of the moiré BZ [67, 68, 69, 70],
which are determined by the difference between the reciprocal lattice vectors
in each layer.

To understand how both $\phi_{j}$ and $\mathbf{q}_{j}$ may influence the
emergence of flat bands, we consider the leading order self-energy given by
Eq. (16), but keeping now explicitly its dependence on those parameters. At
the origin point $k=0=\omega$ one has

$\Sigma_{\mathbf{0}}^{\left(1\right)}\left(0\right)=-6vk_{\theta}\alpha_{0}\alpha_{1}\cos\varphi,$
(28)

where

$\varphi=\theta_{\mathbf{q}_{j}}-\phi_{j}$ (29)

is the $j$-independent difference between the angles $\theta_{\mathbf{q}_{j}}$
of the momentum transfer vectors and the phases $\phi_{j}$. Therefore in
general $\Sigma_{\mathbf{0}}\left(0\right)$ only vanishes in the TBG case
$\varphi=90^{\circ}$ (see Fig. 1), which implies that only then there is a
zero energy solution to Eq. (19) at $k=0$; otherwise the Dirac points are
shifted in energy.

The energy shift of the Dirac points modifies the computation of the
corresponding quasiparticle velocity $v^{\star}$, since then the derivatives
of the self-energy must be evaluated at $\omega$ such that
$\omega-\Sigma_{\mathbf{0}}\left(\omega\right)=0$, which in general depends on
$\varphi$. The situation simplifies in the particular chiral limit
$\alpha_{0}=0$, where Eq. (28) gives $\Sigma_{\mathbf{0}}\left(0\right)=0$ for
all $\varphi$. In that case, the quasiparticle velocity in terms of $\varphi$
reads (cf. Appendix E)

$\frac{v^{\star}}{v}=\frac{\sqrt{1+6\alpha_{1}^{2}\cos\left(2\varphi\right)+9\alpha_{1}^{4}}}{1+3\alpha_{1}^{2}}.$
(30)

Crucially, we see that $v^{\star}$ can only vanish in the TBG case where
$\varphi=90^{\circ}$. The same behavior is numerically observed at higher
perturbation orders, see Fig. 9.

Figure 10: Bottom lowest moiré band for the magic angle $\theta=1.06^{\circ}$
and different orientations $\varphi$ of the momentum transfer vectors in the
continuum model, in steps of $10^{\circ}$ from top to bottom (see inset in
Fig. 9), for the cases (a) $w_{0}=0$ (chiral model) and (b) $w_{0}=0.8w_{1}$.
Flat bands across the mBZ only occur in the TBG case $\varphi=90^{\circ}$.

Eqs. (28) and (30) put the value of $\varphi$ as a key factor in the origin of
magic angles. This observation can be tested by studying the behavior of the
moiré bands for different orientations of the momentum transfer vectors
$\mathbf{q}_{j}$ in the continuum model. The results for the bottom lowest
moiré band is shown in Fig. 10, for both the chiral and nonchiral cases. In
line with Eq. (30), at the magic angle only the TBG case $\varphi=90^{\circ}$
leads to the formation of flat bands across the whole mBZ. This behavior also
occurs in the general nonchiral case, where in addition the Dirac points are
effectively shifted in energy, leading to a stronger suppression of the flat
bands as $\varphi$ changes.

Fig. 10, and in particular Eq. (30), incidentally give a simple explanation as
to why flat bands in TBG seem to disappear in the presence of strain [71, 67,
70], whose effect generally modifies the momentum transfer vectors. Even in
those cases in which the hexagonal symmetry of the moiré pattern in preserved,
as e.g. with a combination of twist and biaxial strain [70], the rotation of
the wave vectors $\mathbf{q}_{j}$ would generally suppress the formation of
flat bands.

##  IV Conclusions

We have developed a many-body perturbation approach to study the properties of
moiré bands in twisted bilayer graphene. This approach differs from previous
studies on several aspects. First, it is build upon a perturbation of the two
uncopled Dirac bands in each monolayer. This provides a natural framework to
deploy perturbation methods and gives a direct physical picture of the moiré-
induced scattering processes. Secondly, our approach treats the perturbation
by the moiré potential within the full framework of a many-body problem. In
this scheme, the moiré potential acts as an effective one-body potential that
transfers electrons from one layer to another, in all possible Dirac bands, by
the exchange of the three wave vectors that determine the moiré Brillouin
zone. The strength of this coupling is naturally accounted by the value of the
twist angle; namely, the decreasing length of the momentum transfer vectors as
the twist angle is lowered translates to an increase in the coupling between
the Dirac points, and therefore in the perturbation strength.

The developed perturbation approach leads to a direct diagrammatic
representation of the scattering processes taking place. By a set of simple
rules, one can obtain, describe or approximate, the moiré-induced self-energy
of the electrons, and from it the moiré bands of interest. This not only
brings the full power of Feynman diagrams into the problem, but also, and
perhaps more importantly, provides a physical picture of the relevant
scattering events. Furthermore, since the perturbation problem is quadratic in
operators, one can naturally obtain, through a straightforward matrix
diagonalization, a closed expression for the Green’s function. As a result,
there is a direct relation between the moiré bands properties at a given shell
cutoff of the continuum model, and the diagrammatic expansion of the Dyson
equation.

The moiré-induced self-energy can be related to different one-particle
properties of interest. In particular, we have obtained the quasiparticle
weight of the moiré bands for electrons in each layer and Dirac band. The
layer distinction can be particularly useful for interpreting photoemission
spectroscopy experiments in moiré heterostructures, where usually the topmost
layer gives the larger contribution to the measured signal. In line with
previous studies, the perturbation approach naturally yields an increasingly
weight of higher moiré bands as the twist angle decreases, reflecting the
contribution of multiple high-momentum scatterings events at successive Dirac
points. Around the magic angle, an interesting band-inversion behavior takes
place: the contribution of electrons in the top and bottom Dirac bands to the
quasiparticle weight of the flat bands is inverted.

The self-energy also gives a natural account for the emergence of flat bands
in twisted bilayer graphene. We have shown that exactly at the magic angle the
static self-energy practically cancels out the uncoupled linear dispersion in
each layer, except only for small deviations around the $\Gamma_{M}$ point in
the nonchiral case, leading to zero-energy flat bands over almost the whole
moiré Brillouin zone. By analyzing how the moiré potential influences such
behavior in the self-energy, we have further shown that flat bands can emerge
only under a particular orientation of the momentum transfer vectors, which is
indeed satisfied in twisted bilayer graphene. For different orientations of
the momentum transfer vectors the Dirac velocity can never vanish. The
geometrical properties of twisted bilayer graphene thus seem to play a crucial
role in the formation of flat bands.

###### Acknowledgements.

I thank Pierre A. Pantaleón, Juan Sebastián Ardenghi and Francisco Guinea for
discussions. This paper was partially supported by the CONICET (Argentina
National Research Council) and Universidad Nacional del Sur (UNS). F. E.
acknowledges a research fellowship from CONICET.

## References

  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018b).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Maximized electron interactions at the magic angle in twisted bilayer graphene, Nature 572, 95 (2019).
  * Oh  _et al._ [2021] M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, X. Liu, K. Watanabe, T. Taniguchi, and A. Yazdani, Evidence for unconventional superconductivity in twisted bilayer graphene, Nature 600, 240 (2021).
  * Jiang  _et al._ [2019] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J. Mao, and E. Y. Andrei, Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene, Nature 573, 91 (2019).
  * Andrei and MacDonald [2020] E. Y. Andrei and A. H. MacDonald, Graphene bilayers with a twist, Nature Materials 19, 1265 (2020).
  * Nimbalkar and Kim [2020] A. Nimbalkar and H. Kim, Opportunities and Challenges in Twisted Bilayer Graphene: A Review, Nano-Micro Letters 12, 126 (2020).
  * Aggarwal  _et al._ [2023] D. Aggarwal, R. Narula, and S. Ghosh, A primer on twistronics: a massless Dirac fermion’s journey to moiré patterns and flat bands in twisted bilayer graphene, Journal of Physics: Condensed Matter 35, 143001 (2023).
  * Lopes dos Santos  _et al._ [2007] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Graphene Bilayer with a Twist: Electronic Structure, Physical Review Letters 99, 256802 (2007).
  * Andrei  _et al._ [2021] E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H. MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young, The marvels of moiré materials, Nature Reviews Materials 6, 201 (2021).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * 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, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019).
  * Taillefer [2010] L. Taillefer, Scattering and pairing in cuprate superconductors, Annual Review of Condensed Matter Physics 1, 51 (2010).
  * Singh [2021] N. Singh, Leading theories of the cuprate superconductivity: A critique, Physica C: Superconductivity and its Applications 580, 1353782 (2021).
  * Reich  _et al._ [2002] S. Reich, J. Maultzsch, C. Thomsen, and P. Ordejón, Tight-binding description of graphene, Physical Review B 66, 035412 (2002).
  * Suárez Morell  _et al._ [2010] E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Flat bands in slightly twisted bilayer graphene: Tight-binding calculations, Physical Review B 82, 121407 (2010).
  * Sboychakov  _et al._ [2015] A. O. Sboychakov, A. L. Rakhmanov, A. V. Rozhkov, and F. Nori, Electronic spectrum of twisted bilayer graphene, Physical Review B 92, 075402 (2015).
  * Lin and Tománek [2018] X. Lin and D. Tománek, Minimum model for the electronic structure of twisted bilayer graphene and related structures, Physical Review B 98, 081410 (2018).
  * Carr  _et al._ [2019a] S. Carr, S. Fang, H. C. Po, A. Vishwanath, and E. Kaxiras, Derivation of Wannier orbitals and minimal-basis tight-binding Hamiltonians for twisted bilayer graphene: First-principles approach, Physical Review Research 1, 033072 (2019a).
  * Lopes dos Santos  _et al._ [2012] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Continuum model of the twisted graphene bilayer, Physical Review B 86, 155449 (2012).
  * Moon and Koshino [2013] P. Moon and M. Koshino, Optical absorption in twisted bilayer graphene, Physical Review B 87, 205404 (2013).
  * Mele [2010] E. J. Mele, Commensuration and interlayer coherence in twisted bilayer graphene, Physical Review B 81, 161405 (2010).
  * Shallcross  _et al._ [2010] S. Shallcross, S. Sharma, E. Kandelaki, and O. A. Pankratov, Electronic structure of turbostratic graphene, Physical Review B 81, 165105 (2010).
  * Carr  _et al._ [2020] S. Carr, S. Fang, and E. Kaxiras, Electronic-structure methods for twisted moiré layers, Nature Reviews Materials 5, 748 (2020).
  * Leconte  _et al._ [2022] N. Leconte, S. Javvaji, J. An, A. Samudrala, and J. Jung, Relaxation effects in twisted bilayer graphene: A multiscale approach, Physical Review B 106, 115410 (2022).
  * Haddadi  _et al._ [2020] F. Haddadi, Q. Wu, A. J. Kruchkov, and O. V. Yazyev, Moiré Flat Bands in Twisted Double Bilayer Graphene, Nano Letters 20, 2410 (2020).
  * Koshino [2015] M. Koshino, Interlayer interaction in general incommensurate atomic layers, New Journal of Physics 17, 015014 (2015).
  * Carr  _et al._ [2019b] S. Carr, S. Fang, Z. Zhu, and E. Kaxiras, Exact continuum model for low-energy electronic states of twisted bilayer graphene, Physical Review Research 1, 013001 (2019b).
  * Guinea and Walet [2019] F. Guinea and N. R. Walet, Continuum models for twisted bilayer graphene: Effect of lattice deformation and hopping parameters, Physical Review B 99, 205134 (2019).
  * Bernevig  _et al._ [2021] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Twisted bilayer graphene. I. Matrix elements, approximations, perturbation theory, and a $k\ifmmode\cdot\else\textperiodcentered\fi{}p$ two-band model, Physical Review B 103, 205411 (2021).
  * Kang and Vafek [2023] J. Kang and O. Vafek, Pseudomagnetic fields, particle-hole asymmetry, and microscopic effective continuum Hamiltonians of twisted bilayer graphene, Physical Review B 107, 075408 (2023).
  * Ceferino and Guinea [2023] A. Ceferino and F. Guinea, Pseudomagnetic fields in fully relaxed twisted bilayer and trilayer graphene 10.48550/arXiv.2311.15052 (2023).
  * Becker  _et al._ [2022] S. Becker, M. Embree, J. Wittsten, and M. Zworski, Mathematics of magic angles in a model of twisted bilayer graphene, Probability and Mathematical Physics 3, 69 (2022).
  * Becker  _et al._ [2023] S. Becker, T. Humbert, and M. Zworski, Integrability in the Chiral Model of Magic Angles, Communications in Mathematical Physics 403, 1153 (2023).
  * Watson  _et al._ [2023] A. B. Watson, T. Kong, A. H. MacDonald, and M. Luskin, Bistritzer–MacDonald dynamics in twisted bilayer graphene, Journal of Mathematical Physics 64, 031502 (2023).
  * San-Jose and Prada [2013] P. San-Jose and E. Prada, Helical networks in twisted bilayer graphene under interlayer bias, Physical Review B 88, 121408 (2013).
  * Efimkin and MacDonald [2018] D. K. Efimkin and A. H. MacDonald, Helical network model for twisted bilayer graphene, Physical Review B 98, 035404 (2018).
  * De Beule  _et al._ [2021] C. De Beule, F. Dominguez, and P. Recher, Network model and four-terminal transport in minimally twisted bilayer graphene, Physical Review B 104, 195410 (2021).
  * Chou and Das Sarma [2023] Y.-Z. Chou and S. Das Sarma, Kondo Lattice Model in Magic-Angle Twisted Bilayer Graphene, Physical Review Letters 131, 026501 (2023).
  * Hu  _et al._ [2023] H. Hu, B. A. Bernevig, and A. M. Tsvelik, Kondo Lattice Model of Magic-Angle Twisted-Bilayer Graphene: Hund’s Rule, Local-Moment Fluctuations, and Low-Energy Effective Theory, Physical Review Letters 131, 026502 (2023).
  * Bennett  _et al._ [2023] D. Bennett, D. T. Larson, L. Sharma, S. Carr, and E. Kaxiras, Twisted bilayer graphene revisited: minimal two-band model for low-energy bands 10.48550/arXiv.2310.12308 (2023).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Magic-Angle Twisted Bilayer Graphene as a Topological Heavy Fermion Problem, Physical Review Letters 129, 047601 (2022).
  * Shi and Dai [2022] H. Shi and X. Dai, Heavy-fermion representation for twisted bilayer graphene systems, Physical Review B 106, 245129 (2022).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Reviews of Modern Physics 81, 109 (2009).
  * Das Sarma  _et al._ [2011] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Electronic transport in two-dimensional graphene, Reviews of Modern Physics 83, 407 (2011).
  * Jain  _et al._ [2016] S. K. Jain, V. Juričić, and G. T. Barkema, Structure of twisted and buckled bilayer graphene, 2D Materials 4, 015018 (2016).
  * Koshino and Nam [2020] M. Koshino and N. N. T. Nam, Effective continuum model for relaxed twisted bilayer graphene and moir\’e electron-phonon interaction, Physical Review B 101, 195425 (2020).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Origin of Magic Angles in Twisted Bilayer Graphene, Physical Review Letters 122, 106405 (2019).
  * Mahan [1990] G. D. Mahan,  _Many-Particle Physics_ (Springer US, 1990).
  * Jishi [2013] R. A. Jishi,  _Feynman Diagram Techniques in Condensed Matter Physics_ (Cambridge University Press, 2013).
  * Coleman [2015] P. Coleman,  _Introduction to Many-Body Physics_ (Cambridge University Press, 2015).
  * Berthod [2018] C. Berthod,  _Spectroscopic Probes of Quantum Matter_ (IOP Publishing, 2018).
  * Long  _et al._ [2023] M. Long, Z. Zhan, P. A. Pantaleón, J. Silva-Guillén, F. Guinea, and S. Yuan, Electronic properties of twisted bilayer graphene suspended and encapsulated with hexagonal boron nitride, Physical Review B 107, 115140 (2023).
  * Long  _et al._ [2022] M. Long, P. A. Pantaleón, Z. Zhan, F. Guinea, J. Silva-Guillén, and S. Yuan, An atomistic approach for the structural and electronic properties of twisted bilayer graphene-boron nitride heterostructures, npj Computational Materials 8, 1 (2022).
  * Yu  _et al._ [2023] G. Yu, Y. Wang, M. I. Katsnelson, and S. Yuan, Origin of the magic angle in twisted bilayer graphene from hybridization of valence and conduction bands, Physical Review B 108, 045138 (2023).
  * Ohta  _et al._ [2012] T. Ohta, J. T. Robinson, P. J. Feibelman, A. Bostwick, E. Rotenberg, and T. E. Beechem, Evidence for Interlayer Coupling and Moiré Periodic Potentials in Twisted Bilayer Graphene, Physical Review Letters 109, 186807 (2012).
  * Nishi  _et al._ [2017] H. Nishi, Y.-i. Matsushita, and A. Oshiyama, Band-unfolding approach to moiré-induced band-gap opening and Fermi level velocity reduction in twisted bilayer graphene, Physical Review B 95, 085420 (2017).
  * Jones  _et al._ [2020] A. J. H. Jones, R. Muzzio, P. Majchrzak, S. Pakdel, D. Curcio, K. Volckaert, D. Biswas, J. Gobbo, S. Singh, J. T. Robinson, K. Watanabe, T. Taniguchi, T. K. Kim, C. Cacho, N. Lanata, J. A. Miwa, P. Hofmann, J. Katoch, and S. Ulstrup, Observation of Electrically Tunable van Hove Singularities in Twisted Bilayer Graphene from NanoARPES, Advanced Materials 32, 2001656 (2020).
  * Lisi  _et al._ [2021] S. Lisi, X. Lu, T. Benschop, T. A. de Jong, P. Stepanov, J. R. Duran, F. Margot, I. Cucchi, E. Cappelli, A. Hunter, A. Tamai, V. Kandyba, A. Giampietri, A. Barinov, J. Jobst, V. Stalman, M. Leeuwenhoek, K. Watanabe, T. Taniguchi, L. Rademaker, S. J. van der Molen, M. P. Allan, D. K. Efetov, and F. Baumberger, Observation of flat bands in twisted bilayer graphene, Nature Physics 17, 189 (2021).
  * Lu  _et al._ [2012] D. Lu, I. M. Vishik, M. Yi, Y. Chen, R. G. Moore, and Z.-X. Shen, Angle-Resolved Photoemission Studies of Quantum Materials, Annual Review of Condensed Matter Physics 3, 129 (2012).
  * Pal and Mele [2013] A. Pal and E. J. Mele, Nodal surfaces in photoemission from twisted bilayer graphene, Physical Review B 87, 205444 (2013).
  * Amorim [2018] B. Amorim, General theoretical description of angle-resolved photoemission spectroscopy of van der Waals structures, Physical Review B 97, 165414 (2018).
  * Li  _et al._ [2022] Y. Li, S. Zhang, F. Chen, L. Wei, Z. Zhang, H. Xiao, H. Gao, M. Chen, S. Liang, D. Pei, L. Xu, K. Watanabe, T. Taniguchi, L. Yang, F. Miao, J. Liu, B. Cheng, M. Wang, Y. Chen, and Z. Liu, Observation of coexisting dirac bands and moiré flat bands in magic‐angle twisted trilayer graphene, Advanced Materials 34, 10.1002/adma.202205996 (2022).
  * Nunn  _et al._ [2023] J. E. Nunn, A. McEllistrim, A. Weston, A. Garcia-Ruiz, M. D. Watson, M. Mucha-Kruczynski, C. Cacho, R. V. Gorbachev, V. I. Fal’ko, and N. R. Wilson, ARPES Signatures of Few-Layer Twistronic Graphenes, Nano Letters 23, 5201 (2023).
  * Jiang  _et al._ [2023] Z. Jiang, K. Hsieh, A. J. H. Jones, P. Majchrzak, C. Sahoo, K. Watanabe, T. Taniguchi, J. A. Miwa, Y. P. Chen, and S. Ulstrup, Revealing flat bands and hybridization gaps in a twisted bilayer graphene device with microARPES, 2D Materials 10, 045027 (2023).
  * Wang  _et al._ [2023] W.-C. Wang, F.-W. Chen, K.-S. Lin, J. T. Hou, H.-C. Lin, and M.-Y. Chou, Origin of magic angles in twisted bilayer graphene: The magic ring 10.48550/arXiv.2309.10026 (2023).
  * Bi  _et al._ [2019] Z. Bi, N. F. Q. Yuan, and L. Fu, Designing flat bands by strain, Physical Review B 100, 035448 (2019).
  * Sinner  _et al._ [2023] A. Sinner, P. A. Pantaleón, and F. Guinea, Strain-Induced Quasi-1D Channels in Twisted Moiré Lattices, Physical Review Letters 131, 166402 (2023).
  * Kögl  _et al._ [2023] M. Kögl, P. Soubelet, M. Brotons-Gisbert, A. V. Stier, B. D. Gerardot, and J. J. Finley, Moiré straintronics: a universal platform for reconfigurable quantum materials, npj 2D Materials and Applications 7, 1 (2023).
  * Escudero  _et al._ [2023] F. Escudero, A. Sinner, Z. Zhan, P. A. Pantaleón, and F. Guinea, Designing moiré patterns by strain 10.48550/arXiv.2309.08671 (2023).
  * Huder  _et al._ [2018] L. Huder, A. Artaud, T. Le Quang, G. T. de Laissardière, A. G. Jansen, G. Lapertot, C. Chapelier, and V. T. Renard, Electronic Spectrum of Twisted Graphene Layers under Heterostrain, Physical Review Letters 120, 156405 (2018).

##  Appendix A Continuum model Hamiltonian in the band basis

Here we give details on the construction of the TBG Hamiltonian in the band
basis. We adopt the original Tripod model [11, 30], with the origin set a
Dirac point (cf. Fig. 2). Following the convention in the main text, we take
the origin Dirac point belonging to the top layer.

The continuum model Hamiltonian in the band basis, which we shall denote as
$H_{B}$, can be directly obtained from the usual one in the sublattice basis
$H$, Eq. (1), by the similarity transformation

$H_{B}=\mathcal{C}^{-1}H_{\mathrm{}}\mathcal{C},$ (31)

where $\mathcal{C}$ is the matrix that diagonalizes the Dirac Hamiltonian
$H_{0}$ (i.e., that without the moiré potential coupling). In momentum space
$\mathcal{C}$ is a matrix whose columns are thus given by the eigenstates of
$H_{0}$,

$\mathcal{C}=\left(\begin{array}[]{ccccc}\Psi_{\mathbf{k}}&0&0&0&\cdots\\\
0&\Psi_{\mathbf{k}+\mathbf{q}_{1}}&0&0&\cdots\\\
0&0&\Psi_{\mathbf{k}+\mathbf{q}_{2}}&0&\cdots\\\
0&0&0&\Psi_{\mathbf{k}+\mathbf{q}_{3}}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right),$ (32)

where [cf. Eq. (6)]

$\Psi_{\mathbf{k}}=\left(\begin{array}[]{cc}\psi_{\mathbf{k},-}&\psi_{\mathbf{k},+}\end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&1\\\
-e^{i\theta_{\mathbf{k}}}&e^{i\theta_{\mathbf{k}}}\end{array}\right).$ (33)

Note that $\mathcal{C}^{-1}=\mathcal{C}^{\dagger}$. The similarity
transformation (31) leads to the moiré matrices in the band basis:

$\left\langle\mathbf{k}_{-}\left|\hat{T}_{j}\right|\mathbf{k}_{+}\right\rangle_{B}=\Psi_{\mathbf{k}_{-}}^{\dagger}T_{j}\Psi_{\mathbf{k}_{+}}\delta_{\mathbf{k}_{-},\mathbf{k}_{+}+\mathbf{q}_{j}}$
(34)

where $T_{j}$ is given by Eq. (3). The matrix components in the above equation
lead to Eq. (5) in the main text. The above result can be also obtained by
making a direct basis change of the moiré matrix elements, using the relation
$\left|\mathbf{k},s\right\rangle=\left(\left|A\right\rangle+se^{i\theta_{\mathbf{k}}}\left|B\right\rangle\right)/\sqrt{2}$
between the band basis and the sublattice basis.

The computation of the Green’s function through Eq. (11) is done by a matrix
inversion of the continuum model Hamiltonian in the band basis. Given the
relation (31), this inversion can be also directly obtained from the original
TBG Hamiltonian $H$, Eq. (1) (in momentum space), and a basis change:

$\mathcal{G}\left(i\omega_{n}\right)=\mathcal{C}^{\dagger}\left(i\omega_{n}\mathbb{I}-H\right)^{-1}\mathcal{C},$
(35)

In particular, for the two cases
$\gamma=\gamma^{\prime}=\left(+,\mathbf{k},s=\pm\right)$, the matrix
components of $\left(i\omega_{n}\mathbb{I}-H\right)_{\gamma\gamma}^{-1}$,
given our convention, are $\gamma\gamma=\left(1,1\right)$ and
$\gamma\gamma=\left(2,2\right)$ for the bottom $\left(s=-1\right)$ and top
$\left(s=1\right)$ Dirac band, respectively. Defining

$\tilde{H}\left(i\omega_{n}\right)=i\omega_{n}\mathbb{I}-H,$ (36)

we thus have

$\displaystyle\mathcal{G}_{+,\mathbf{k},-}\left(i\omega_{n}\right)$
$\displaystyle=\mathcal{C}_{1\mu}^{\dagger}\tilde{H}_{\mu\nu}^{-1}\left(i\omega_{n}\right)\mathcal{C}_{\nu
1},$ (37) $\displaystyle\mathcal{G}_{+,\mathbf{k},+}\left(i\omega_{n}\right)$
$\displaystyle=\mathcal{C}_{2\mu}^{\dagger}\tilde{H}_{\mu\nu}^{-1}\left(i\omega_{n}\right)\mathcal{C}_{\nu
2},$ (38)

where repeated indices are to be summed over. Here it should be understand
that $\tilde{H}_{\mu\nu}^{-1}\left(i\omega_{n}\right)$ are the
$\left(\mu,\nu\right)$ components of $\tilde{H}^{-1}\left(i\omega_{n}\right)$.
From Eqs. (32) and (33) it follows that

$\displaystyle\mathcal{G}_{+,\mathbf{k},\pm}\left(i\omega_{n}\right)$
$\displaystyle=\frac{1}{2}\left(\tilde{H}_{11}^{-1}+\tilde{H}_{22}^{-1}\pm\tilde{H}_{12}^{-1}e^{i\theta_{\mathbf{k}}}\pm\tilde{H}_{21}^{-1}e^{-i\theta_{\mathbf{k}}}\right).$
(39)

This provides direct way to obtain the Green’s function by working entirely in
the sublattice basis. From it we readily see that

$\displaystyle\mathcal{G}_{+,\mathbf{k},+}\left(i\omega_{n}\right)+\mathcal{G}_{+,\mathbf{k},-}\left(i\omega_{n}\right)$
$\displaystyle=\tilde{H}_{11}^{-1}+\tilde{H}_{22}^{-1},$ (40)
$\displaystyle\mathcal{G}_{+,\mathbf{k},+}\left(i\omega_{n}\right)-\mathcal{G}_{+,\mathbf{k},-}\left(i\omega_{n}\right)$
$\displaystyle=\tilde{H}_{12}^{-1}e^{i\theta_{\mathbf{k}}}+\tilde{H}_{21}^{-1}e^{-i\theta_{\mathbf{k}}}.$
(41)

It should be noted that dependence of $\mathcal{G}$ on the moiré potential
$\mathcal{T}_{j,\mathbf{k},s,s^{\prime}}$ in the band basis, when computed
using Eq. (39), would not be explicit; it would rather be implicitly contained
in the couplings between the Hamiltonian $\sim\bm{\sigma}\cdot\mathbf{k}$ of
different Dirac points. From a diagrammatic point of view, and physical
interpretation, it is thus more convenient to compute $\mathcal{G}$ by
directly working in the band representation.

##  Appendix B Perturbation expansion of the moiré potential

By treating the moiré potential $\hat{T}$, Eq. (8), as a perturbation, the
Green’s function given by Eq. (10) can be expanded as [49, 50, 51] (we set
$\hbar=1$)

$\mathcal{G}_{\gamma\gamma^{\prime}}\left(\tau\right)=-\left\langle
T_{\tau}\hat{c}_{\gamma}\left(\tau\right)\hat{c}_{\gamma^{\prime}}^{\dagger}\left(0\right)\hat{U}\left(\beta,0\right)\right\rangle_{0,c},$
(42)

where $\beta=1/k_{B}T$ and $\left\langle\cdots\right\rangle_{0,c}$ is the
ensemble average over the noninteracting system, taking into account only the
contribution of connected diagrams. The evolution operator reads

$\hat{U}\left(\tau,\tau^{\prime}\right)=\sum_{l=0}^{\infty}\frac{1}{l!}\left(-1\right)^{l}\int_{\tau^{\prime}}^{\tau}d\tau_{1}\ldots\int_{\tau^{\prime}}^{\tau}d\tau_{l}T_{\tau}\left[\hat{T}\left(\tau_{1}\right)\ldots\hat{T}\left(\tau_{l}\right)\right].$
(43)

The perturbation expansion of the evolution operators naturally accounts for
the shell truncation of the continuum model Hamiltonian. Following Sec. III,
we are interested in the case $\gamma=\gamma^{\prime}$ in Eq. (42); we note
$\mathcal{G}_{\gamma\gamma}\left(\tau\right)=\mathcal{G}_{\gamma}\left(\tau\right)$.
In that situation, since the moiré potential $\hat{T}$ always exchanges one
electron from one layer to another, only even orders in $l$ contribute to the
expansion of the evolution operator. For conciseness we shall consider the
particular case $\gamma=\left(+,\mathbf{k},s\right)$. The leading order
correction to the Green’s function (second order in $l$) reads

$\mathcal{G}_{+,\mathbf{k},s}^{\left(2\right)}\left(\tau\right)=-\frac{1}{2!}\left(-1\right)^{2}\int_{0}^{\beta}d\tau_{1}\int_{0}^{\beta}d\tau_{2}\left\langle
T_{\tau}\hat{c}_{+,\mathbf{k},s}\left(\tau\right)c_{+,\mathbf{k},s}^{\dagger}\left(0\right)\hat{T}\left(\tau_{1}\right)\hat{T}\left(\tau_{2}\right)\right\rangle_{0,c}.$
(44)

Replacing Eq. (8) leads to the following four kinds of ensemble averages:

$\displaystyle C_{1}$ $\displaystyle=\left\langle
T_{\tau}\hat{c}_{+,\mathbf{k},s}\left(\tau\right)\hat{c}_{+,\mathbf{k},s}^{\dagger}\left(0\right)\hat{c}_{-,\mathbf{k}_{1}+\mathbf{q}_{j_{1}},s_{1}^{\prime}}^{\dagger}\left(\tau_{1}\right)\hat{c}_{+,\mathbf{k}_{1},s_{1}}\left(\tau_{1}\right)\hat{c}_{-,\mathbf{k}_{2}+\mathbf{q}_{j_{2}},s_{2}^{\prime}}^{\dagger}\left(\tau_{2}\right)\hat{c}_{+,\mathbf{k}_{2},s_{2}}\left(\tau_{2}\right)\right\rangle_{0,c},$
(45) $\displaystyle C_{2}$ $\displaystyle=\left\langle
T_{\tau}\hat{c}_{+,\mathbf{k},s}\left(\tau\right)\hat{c}_{+,\mathbf{k},s}^{\dagger}\left(0\right)\hat{c}_{-,\mathbf{k}_{1}+\mathbf{q}_{j_{1}},s_{1}^{\prime}}^{\dagger}\left(\tau_{1}\right)\hat{c}_{+,\mathbf{k}_{1},s_{1}}\left(\tau_{1}\right)\hat{c}_{+,\mathbf{k}_{2},s_{2}}^{\dagger}\left(\tau_{2}\right)\hat{c}_{-,\mathbf{k}_{2}+\mathbf{q}_{j_{2}},s_{2}^{\prime}}\left(\tau_{2}\right)\right\rangle_{0,c},$
(46) $\displaystyle C_{3}$ $\displaystyle=\left\langle
T_{\tau}\hat{c}_{+,\mathbf{k},s}\left(\tau\right)\hat{c}_{+,\mathbf{k},s}^{\dagger}\left(0\right)\hat{c}_{+,\mathbf{k}_{1},s_{1}}^{\dagger}\left(\tau_{1}\right)\hat{c}_{-,\mathbf{k}_{1}+\mathbf{q}_{j_{1}},s_{1}^{\prime}}\left(\tau_{1}\right)\hat{c}_{-,\mathbf{k}_{2}+\mathbf{q}_{j_{2}},s_{2}^{\prime}}^{\dagger}\left(\tau_{2}\right)\hat{c}_{+,\mathbf{k}_{2},s_{2}}\left(\tau_{2}\right)\right\rangle_{0,c},$
(47) $\displaystyle C_{4}$ $\displaystyle=\left\langle
T_{\tau}\hat{c}_{+,\mathbf{k},s}\left(\tau\right)\hat{c}_{+,\mathbf{k},s}^{\dagger}\left(0\right)\hat{c}_{+,\mathbf{k}_{1},s_{1}}^{\dagger}\left(\tau_{1}\right)\hat{c}_{-,\mathbf{k}_{1}+\mathbf{q}_{j_{1}},s_{1}^{\prime}}\left(\tau_{1}\right)\hat{c}_{+,\mathbf{k}_{2},s_{2}}^{\dagger}\left(\tau_{2}\right)\hat{c}_{-,\mathbf{k}_{2}+\mathbf{q}_{j_{2}},s_{2}^{\prime}}\left(\tau_{2}\right)\right\rangle_{0,c}.$
(48)

Applying Wick’s theorem, considering the anticommutation relation
$\left\\{\hat{c}_{\ell,\mathbf{k},s},\hat{c}_{\ell^{\prime},\mathbf{k}^{\prime},s^{\prime}}^{\dagger}\right\\}=\delta_{\ell,\ell^{\prime}}\delta_{\mathbf{k},\mathbf{k}^{\prime}}\delta_{s,s^{\prime}}$
(all others anticommutators vanish), and taking into account only those
contractions that represent connected diagrams, yields $C_{1}=0=C_{4}$ and

$\displaystyle C_{2}$
$\displaystyle=-\mathcal{G}_{\mathbf{k},s}^{0}\left(\tau-\tau_{2}\right)\delta_{\mathbf{k},\mathbf{k}_{2}}\delta_{s,s_{2}}\mathcal{G}_{\mathbf{k},s}^{0}\left(\tau_{1}\right)\delta_{\mathbf{k},\mathbf{k}_{1}}\delta_{s,s_{1}}\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j_{1}},s_{1}^{\prime}}^{0}\left(\tau_{2}-\tau_{1}\right)\delta_{\mathbf{k}_{2}+\mathbf{q}_{j_{2}},\mathbf{k}_{1}+\mathbf{q}_{j_{1}}}\delta_{s^{\prime}_{2},s^{\prime}_{1}},$
(49) $\displaystyle C_{3}$
$\displaystyle=-\mathcal{G}_{\mathbf{k},s}^{0}\left(\tau-\tau_{1}\right)\delta_{\mathbf{k},\mathbf{k}_{1}}\delta_{s,s_{1}}\mathcal{G}_{\mathbf{k},s}^{0}\left(\tau_{2}\right)\delta_{\mathbf{k},\mathbf{k}_{2}}\delta_{s,s_{2}}\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j_{2}},s_{2}^{\prime}}^{0}\left(\tau_{1}-\tau_{2}\right)\delta_{\mathbf{k}_{2}+\mathbf{q}_{j_{2}},\mathbf{k}_{1}+\mathbf{q}_{j_{1}}}\delta_{s^{\prime}_{2},s^{\prime}_{1}},$
(50)

where

$\left\langle
T_{\tau}\hat{c}_{\ell,\mathbf{k},s}\left(\tau\right)\hat{c}_{\ell^{\prime},\mathbf{k}^{\prime},s^{\prime}}^{\dagger}\left(\tau^{\prime}\right)\right\rangle_{0}=-\delta_{\ell,\ell^{\prime}}\delta_{\mathbf{k},\mathbf{k}^{\prime}}\delta_{s,s^{\prime}}\mathcal{G}_{\mathbf{k},s}^{0}\left(\tau-\tau^{\prime}\right)$
(51)

is the noninteracting Green’s function. Both $C_{2}$ and $C_{3}$ give
topologically equivalent diagrams. Thus

$\mathcal{G}_{+,\mathbf{k},s}^{\left(2\right)}\left(\tau\right)=\sum_{j,s^{\prime}}\left|\mathcal{T}_{j,\mathbf{k},s,s^{\prime}}\right|^{2}\int_{0}^{\beta}d\tau_{1}\int_{0}^{\beta}d\tau_{2}\mathcal{G}_{\mathbf{k},s}^{0}\left(\tau-\tau_{2}\right)\mathcal{G}_{\mathbf{k},s}^{0}\left(\tau_{1}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j},s^{\prime}}^{0}\left(\tau_{2}-\tau_{1}\right).$
(52)

Using the Fourier transform

$\mathcal{G}_{\mathbf{k},s}^{0}\left(\tau\right)=\frac{1}{\beta}\sum_{n}\mathcal{G}_{\mathbf{k},s}^{0}\left(i\omega_{n}\right)e^{-i\omega_{n}\tau},\qquad\omega_{n}=\left(2n+1\right)\pi/\beta,$
(53)

and the relation

$\int_{0}^{\beta}d\tau\exp\left[i\left(\omega_{n_{j}}-\omega_{n_{j^{\prime}}}\right)\tau\right]=\beta\delta_{\omega_{n_{j}},\omega_{n_{j^{\prime}}}},$
(54)

we get

$\mathcal{G}_{+,\mathbf{k},s}^{\left(2\right)}\left(i\omega_{n}\right)=\mathcal{G}_{\mathbf{k},s}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k},s}^{0}\left(i\omega_{n}\right)\sum_{j,s^{\prime}}\left|\mathcal{T}_{j,\mathbf{k},s,s^{\prime}}\right|^{2}\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j},s^{\prime}}^{0}\left(i\omega_{n}\right),$
(55)

where

$\mathcal{G}_{\mathbf{k},s}^{0}\left(i\omega_{n}\right)=\frac{1}{i\omega_{n}-\epsilon_{\mathbf{k},s}}.$
(56)

Since the perturbation is quadratic in the operators, there is no dependence
on the statistic of the particles. One readily sees that Eq. (55) represents
the correction to the Green’s function coming from the _intraband_
contribution to the leading order self-energy, cf. Eq. (17) (the other
interband contribution comes from higher orders in the perturbation
expansion).

From the above derivation one can easily infer the Feynman rules for the
system, see Fig. 3. As an example of applying the Feynman rules, the
contribution to the self-energy of the five diagrams in Fig. 4 (in the same
order as shown) read

$\displaystyle\Sigma_{\mathbf{k},s}^{\left(l=2\right)}\left(i\omega_{n}\right)$
$\displaystyle=\sum_{j}\sum_{r}\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},s,r}^{*}\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j},r}^{0}\left(i\omega_{n}\right),$
(57)
$\displaystyle\Sigma_{\mathbf{k},s}^{\left(l=4,a\right)}\left(i\omega_{n}\right)$
$\displaystyle=\sum_{j}\sum_{j^{\prime}}\sum_{r_{1},r_{2}}\mathcal{T}_{j,\mathbf{k},s,r_{1}}\mathcal{T}_{j,\mathbf{k},-s,r_{1}}^{*}\mathcal{T}_{j^{\prime},\mathbf{k},-s,r_{2}}\mathcal{T}_{j^{\prime},\mathbf{k},s,r_{2}}^{*}\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j},r_{1}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k},-s}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j^{\prime}},r_{2}}^{0}\left(i\omega_{n}\right),$
(58)
$\displaystyle\Sigma_{\mathbf{k},s}^{\left(l=4,b\right)}\left(i\omega_{n}\right)$
$\displaystyle=\sum_{j}\sum_{j^{\prime}\neq
j}\sum_{r_{1},r_{2}}\sum_{s_{1}}\mathcal{T}_{j,\mathbf{k},s,r_{1}}\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{1},r_{1}}^{*}\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{1},r_{2}}\mathcal{T}_{j,\mathbf{k},s,r_{2}}^{*}$
$\displaystyle\qquad\times\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j},r_{1}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{1}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j},r_{2}}^{0}\left(i\omega_{n}\right),$
(59)
$\displaystyle\Sigma_{\mathbf{k},s}^{\left(l=6,a\right)}\left(i\omega_{n}\right)$
$\displaystyle=\sum_{j}\sum_{j^{\prime}\neq j}\sum_{j^{\prime\prime}\neq
j^{\prime}}\sum_{r_{1},r_{2},r_{3}}\sum_{s_{1},s_{2}}\mathcal{T}_{j,\mathbf{k},s,r_{1}}\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{1},r_{1}}^{*}\mathcal{T}_{j^{\prime\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{1},r_{2}}$
$\displaystyle\times\mathcal{T}_{j^{\prime\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{2},r_{2}}^{*}\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{2},r_{3}}\mathcal{T}_{j,\mathbf{k},s,r_{3}}^{*}$
$\displaystyle\times\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j},r_{1}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{1}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}}+\mathbf{q}_{j^{\prime\prime}},r_{2}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{2}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j},r_{3}}^{0}\left(i\omega_{n}\right),$
(60)
$\displaystyle\Sigma_{\mathbf{k},s}^{\left(l=6,b\right)}\left(i\omega_{n}\right)$
$\displaystyle=\sum_{j}\sum_{j^{\prime}\neq j}\sum_{j^{\prime\prime}\neq
j,j^{\prime}}\sum_{r_{1},r_{2},r_{3}}\sum_{s_{1},s_{2}}\mathcal{T}_{j,\mathbf{k},s,r_{1}}\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{1},r_{1}}^{*}\mathcal{T}_{j^{\prime\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{1},r_{2}}$
$\displaystyle\times\mathcal{T}_{j,\mathbf{k}-\mathbf{q}_{j^{\prime}}+\mathbf{q}_{j^{\prime\prime}},s_{2},r_{2}}^{*}\mathcal{T}_{j^{\prime},\mathbf{k}-\mathbf{q}_{j^{\prime}}+\mathbf{q}_{j^{\prime\prime}},s_{2},r_{3}}\mathcal{T}_{j^{\prime\prime},\mathbf{k},s,r_{3}}^{*}$
$\displaystyle\times\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j},r_{1}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},s_{1}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}}+\mathbf{q}_{j^{\prime\prime}},r_{2}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}-\mathbf{q}_{j^{\prime}}+\mathbf{q}_{j^{\prime\prime}},s_{2}}^{0}\left(i\omega_{n}\right)\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j^{\prime\prime}},r_{3}}^{0}\left(i\omega_{n}\right).$
(61)

##  Appendix C Diagrammatic computation of the self-energy at second shell

As noted in the main text, the self-energy for a momentum truncation at the
second shell (see Fig. 2) can be directly computed by inverting the
corresponding $20\times 20$ TBG Hamiltonian, but this leads to a highly
cumbersome expression, which is hard to work with and gives little physical
insight. For this reason, here we rather compute the self-energy by means of
only diagrammatic methods. This approach highlights the involved scattering
processes and provides a workable analytical expression from which one can
further study different properties of the moiré bands, such as the
quasiparticle velocity renormalization (Appendix D).

Figure 11: Left: Diagrammatic representation of the moiré-induced self-energy
in the top layer, up to a momentum truncation at the second shell (cf. Fig.
2). Dot-dashed Green’s functions represent the dressed propagator for an
electron in the bottom layer, with momentum $\mathbf{k}+\mathbf{q}_{j}$, by
the coupling to an electron in the top layer with momentum
$\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}}$, where
$\mathbf{q}_{j^{\prime}}\neq\mathbf{q}_{j^{\prime}}$. Right: schematic
hexagonal picture of all the scattering processes up to a second shell
truncation (top panel), and the dressed dot-dashed Green’s function (bottom
panel).

We consider, as in the main text, the self-energy for the top layer. The
diagrammatic representation that accounts for all scattering processes with
momentum exchange up to
$\Delta\mathbf{k}\sim\left|\mathbf{q}_{1}-\mathbf{q}_{2}\right|$ (second shell
in Fig. 2) is shown in Fig. 11. Applying the Feynman rules we get

$\displaystyle\Sigma_{+,\mathbf{k},s}^{\left(2\right)}\left(i\omega_{n}\right)$
$\displaystyle=\sum_{j,r_{1},r_{2}}\mathcal{T}_{j,\mathbf{k},s,r_{1}}\mathcal{T}_{j,\mathbf{k},-s,r_{2}}^{*}\tilde{\mathcal{G}}_{\mathbf{k}+\mathbf{q}_{j},r_{1},r_{2}}\left(i\omega_{n}\right)\sum_{j,r_{1},r_{2}}\mathcal{T}_{j,\mathbf{k},-s,r_{2}}\mathcal{T}_{j,\mathbf{k},s,r_{1}}^{*}\tilde{\mathcal{G}}_{\mathbf{k}+\mathbf{q}_{j},r_{1},r_{2}}\left(i\omega_{n}\right)$
$\displaystyle\qquad\times\left(i\omega_{n}-\epsilon_{\mathbf{k},-s}-\sum_{j,r_{1},r_{2}}\mathcal{T}_{j,\mathbf{k},-s,r_{1}}\mathcal{T}_{j,\mathbf{k},-s,r_{2}}^{*}\tilde{\mathcal{G}}_{\mathbf{k}+\mathbf{q}_{j},r_{1},r_{2}}\left(i\omega_{n}\right)\right)^{-1}$
$\displaystyle\qquad\qquad+\sum_{j,r_{1},r_{2}}\mathcal{T}_{j,\mathbf{k},s,r_{1}}\mathcal{T}_{j,\mathbf{k},s,r_{2}}^{*}\tilde{\mathcal{G}}_{\mathbf{k}+\mathbf{q}_{j},r_{1},r_{2}}\left(i\omega_{n}\right).$
(62)

Here we defined
$\tilde{\mathcal{G}}_{\mathbf{k}+\mathbf{q}_{j},r_{1},r_{2}}\left(i\omega_{n}\right)$
as the Green’s function for the propagation of
$\left|-,\mathbf{k}+\mathbf{q}_{j},r_{1}\right\rangle$ to
$\left|-,\mathbf{k}+\mathbf{q}_{j},r_{2}\right\rangle$, taking into account
only the interactions with the other layer with
$\mathbf{q}_{j^{\prime}}\neq\mathbf{q}_{j}$. It can be obtained
diagrammatically by differentiating the cases $r_{1}=r_{2}$ and $r_{1}=-r_{2}$
(cf. Fig. 11), with the result

$\displaystyle\tilde{\mathcal{G}}_{\mathbf{k}+\mathbf{q}_{j},r,r}\left(i\omega_{n}\right)$
$\displaystyle=\frac{1}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}-\Sigma^{\prime\prime}_{\mathbf{k}+\mathbf{q}_{j},r,r}\left(i\omega_{n}\right)},$
(63)
$\displaystyle\tilde{\mathcal{G}}_{\mathbf{k}+\mathbf{q}_{j},r,-r}\left(i\omega_{n}\right)$
$\displaystyle=\frac{1}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j},-r}-\Sigma^{\prime\prime}_{\mathbf{k}+\mathbf{q}_{j},-r,-r}\left(i\omega_{n}\right)}\frac{\Sigma^{\prime\prime}_{\mathbf{k}+\mathbf{q}_{j},r,-r}\left(i\omega_{n}\right)}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}.$
(64)

where

$\displaystyle\Sigma^{\prime\prime}_{\mathbf{k}+\mathbf{q}_{j},r_{1},r_{2}}\left(i\omega_{n}\right)$
$\displaystyle=\sum_{j^{\prime}\neq
j,r^{\prime}}\frac{\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime},r_{1}}^{*}\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime},-r_{2}}}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime}}}\sum_{j^{\prime}\neq
j,r^{\prime}}\frac{\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime},-r_{2}}^{*}\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime},r_{2}}}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime}}}$
$\displaystyle\times\left(i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j},-r_{2}}-\sum_{j^{\prime}\neq
j,r^{\prime}}\frac{\left|\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime},-r_{2}}\right|^{2}}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime}}}\right)^{-1}+\sum_{j^{\prime}\neq
j,r^{\prime}}\frac{\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime},r_{1}}^{*}\mathcal{T}_{j^{\prime},\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime},r_{2}}}{i\omega_{n}-\epsilon_{\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}},r^{\prime}}}.$
(65)

It is important to note that $j^{\prime}\neq j$ in the summation over
$j^{\prime}$, in order to not overcount scattering processes. We have checked
numerically that the self-energy given by Eq. (62) agrees with the one
obtained from the general expression given by Eq. (13), for a momenta cutoff
at the second shell.

Note that Eq. (62) reduces to Eq. (16) when the momentum is truncated at the
first shell and therefore one does not consider the dressing of the propagator
with momentum $\mathbf{k}+\mathbf{q}_{j}$ by the interaction with higher order
Dirac points, i.e., when
$\tilde{\mathcal{G}}_{\mathbf{k}+\mathbf{q}_{j},r_{1},r_{2}}\rightarrow\mathcal{G}_{\mathbf{k}+\mathbf{q}_{j},r_{1}}^{0}\delta_{r_{1},r_{2}}$.

##  Appendix D Dirac velocity renormalization

The renormalized velocity $v^{\star}$ of the lowest moiré bands at the Dirac
point can be directly computed from the self-energy (Sec. III.3). Here we
give, for completeness and reference, details of the analytical calculation of
$v^{\star}$ up to the second shell momentum cutoff, using Eqs. (16) and (62).
Of course, this approach complements and yields the same results as usually
obtained by traditional perturbation methods applied to a truncated TBG
Hamiltonian [11, 48, 30]. As we are interested in the $k$ and $\omega$
derivatives, around the zero point $k=0=\omega$, of the real part of the
retarded self-energy, we shall directly replace $i\omega_{n}\rightarrow\omega$
when doing the analytical continuation.

Following Sec. III.3 in the main text, we consider the self-energy for the top
layer, but keeping here in general both bands $s$. The generalization of Eq.
(25) for both bands is

$\frac{v^{\star}}{v}=\frac{1+s\left(\partial\Sigma^{\prime}_{\mathbf{k},s}/\partial
k\right)_{0}/v}{1-\left(\partial\Sigma^{\prime}_{\mathbf{k},s}/\partial\omega\right)_{0}}.$
(66)

To condense the notation we will note $\Sigma^{\prime}_{\mathbf{k},s}=\Sigma$.

###  D.1 First shell

The self-energy when the momentum is truncated at the first shell is given by
Eq. (16). From Eq. (5) one has

$\displaystyle\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}$
$\displaystyle=w_{0}^{2}\frac{1+sr\cos\left(\theta_{\mathbf{k}}-\theta_{\mathbf{k}+\mathbf{q}_{j}}\right)}{2}+w_{1}^{2}\frac{1+sr\cos\left(\theta_{\mathbf{k}}+\theta_{\mathbf{k}+\mathbf{q}_{j}}-2\phi_{j}\right)}{2}$
$\displaystyle\qquad+sw_{0}w_{1}\left[\cos\left(\theta_{\mathbf{k}}-\phi_{j}\right)+sr\cos\left(\theta_{\mathbf{k}+\mathbf{q}_{j}}-\phi_{j}\right)\right],$
(67)
$\displaystyle\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}$
$\displaystyle=iw_{0}^{2}\frac{sr\sin\left(\theta_{\mathbf{k}}-\theta_{\mathbf{k}+\mathbf{q}_{j}}\right)}{2}+iw_{1}^{2}\frac{sr\sin\left(\theta_{\mathbf{k}}+\theta_{\mathbf{k}+\mathbf{q}_{j}}-2\phi_{j}\right)}{2}$
$\displaystyle\qquad+isw_{0}w_{1}\sin\left(\theta_{\mathbf{k}}-\phi_{j}\right).$
(68)

It follows that

$\left(\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}=\left(\frac{\partial}{\partial
k}\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}=\left(\frac{\partial}{\partial\omega}\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}=0,$
(69)

where the subindex “0” implies evaluation at $k=0=\omega$. Consequently the
interband term given by Eq. (18) makes no contribution to the first order
derivatives $\left(\partial\Sigma/\partial k\right)_{0}$ and
$\left(\partial\Sigma/\partial\omega\right)_{0}$. We thus continue considering
only the intraband contribution given by Eq. (17), which yields

$\displaystyle\left(\frac{\partial\Sigma}{\partial k}\right)_{0}$
$\displaystyle=-\sum_{j=1}^{3}\sum_{r=\pm
1}\left(\frac{\partial\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{\partial
k}\right)_{0}\frac{1}{rvk_{\theta}}+\sum_{j=1}^{3}\sum_{r=\pm
1}\frac{\left|\mathcal{T}_{j,\mathbf{0},s,r}\right|^{2}}{\left(vk_{\theta}\right)^{2}}\left(\frac{\partial\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}{\partial
k}\right)_{0}.$ (70)

where
$\mathcal{T}_{j,\mathbf{0},s,r}=\mathcal{T}_{j,\mathbf{k}=\mathbf{0},s,r}$ and
we replaced $\epsilon_{\mathbf{q}_{j}}=vk_{\theta}$, which is independent of
$j$ ($\left|\mathbf{q}_{j}\right|=k_{\theta}$ being the moiré BZ length). We
get

$\displaystyle\sum_{j=1}^{3}\sum_{r=\pm
1}\left(\frac{\partial\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{\partial
k}\right)_{0}\frac{1}{rvk_{\theta}}$
$\displaystyle=\frac{3}{2}sv\left(\alpha_{0}^{2}+\alpha_{1}^{2}\right),$ (71)
$\displaystyle\sum_{j=1}^{3}\sum_{r=\pm
1}\frac{\left|\mathcal{T}_{j,\mathbf{0},s,r}\right|^{2}}{\left(vk_{\theta}\right)^{2}}\left(\frac{\partial\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}{\partial
k}\right)_{0}$
$\displaystyle=\frac{3}{2}sv\left(\alpha_{0}^{2}-\alpha_{1}^{2}\right),$ (72)

where $\alpha_{i}=w_{i}/vk_{\theta}$. Thus,

$\displaystyle\left(\frac{\partial\Sigma}{\partial k}\right)_{0}$
$\displaystyle=-3sv\alpha_{1}^{2}.$ (73)

For the $\omega$-derivative we get

$\left(\frac{\partial\Sigma}{\partial\omega}\right)_{0}=-\sum_{j=1}^{3}\sum_{r=\pm
1}\frac{\left|\mathcal{T}_{j,\mathbf{0},s,r}\right|^{2}}{\left(vk_{\theta}\right)^{2}}=-3\left(\alpha_{0}^{2}+\alpha_{1}^{2}\right).$
(74)

Replacing in Eq. (66) leads to the well-known leading order renormalization

$\frac{v^{\star}}{v}=\frac{1-3\alpha_{1}^{2}}{1+3\alpha_{0}^{2}+3\alpha_{1}^{2}}.$
(75)

The case $w_{0}=w_{1}$ gives the original result first obtained in Ref. [11].

###  D.2 Second shell

Truncating the momentum at the second shell yields the self-energy obtained in
Appendix C. The moiré potentials in Eq. (65) are given by Eqs. (67) and (68),
upon the replacement $j\rightarrow j^{\prime}$, $s\rightarrow r^{\prime}$ and
$\mathbf{k}\rightarrow\mathbf{k}+\mathbf{q}_{j}-\mathbf{q}_{j^{\prime}}$. From
Eq. (65) it follows that, at $k=0=\omega$,

$\displaystyle\Sigma^{\prime\prime}_{\mathbf{q}_{j},r,r}\left(0\right)$
$\displaystyle=-rvk_{\theta}\alpha_{0}^{2},$ (76)
$\displaystyle\Sigma^{\prime\prime}_{\mathbf{q}_{j},r,-r}\left(0\right)$
$\displaystyle=0,$ (77)

and therefore, at the same point, Eqs. (63) and (64) reduce to

$\displaystyle\tilde{\mathcal{G}}_{\mathbf{q}_{j},r,r}\left(0\right)$
$\displaystyle=-\frac{r}{vk_{\theta}}\frac{1}{1-\alpha_{0}^{2}},$ (78)
$\displaystyle\tilde{\mathcal{G}}_{\mathbf{q}_{j},r,-r}\left(0\right)$
$\displaystyle=0.$ (79)

Moreover, one has

$\displaystyle\sum_{j,r}\mathcal{T}_{j,\mathbf{0},s,r}\mathcal{T}_{j,\mathbf{0},-s,r}^{*}\tilde{\mathcal{G}}_{\mathbf{q}_{j},r,r}\left(0\right)$
$\displaystyle=0.$ (80)

Given these relations, from Eq. (62) we have

$\displaystyle\left(\frac{\partial\Sigma}{\partial k}\right)_{0}$
$\displaystyle=\sum_{j,r}\left(\frac{\partial\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{\partial
k}\right)_{0}\tilde{\mathcal{G}}_{\mathbf{q}_{j},r,r}\left(0\right)+\sum_{j,r_{1},r_{2}}\mathcal{T}_{j,\mathbf{0},s,r_{1}}\mathcal{T}_{j,\mathbf{0},s,r_{2}}^{*}\left(\frac{\partial\tilde{\mathcal{G}}_{\mathbf{k}+\mathbf{q}_{j},r_{1},r_{2}}}{\partial
k}\right)_{0}.$ (81)

The first term is the same as the one appearing in the first shell case, Eq.
(71), but with the inclusion of the factor
$\sim\left(1-\alpha_{0}^{2}\right)^{-1}$ in Eq. (78), so

$\displaystyle\sum_{j,r}\left(\frac{\partial\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{\partial
k}\right)_{0}\tilde{\mathcal{G}}_{\mathbf{q}_{j},r,r}\left(0\right)$
$\displaystyle=-\frac{3}{2}sv\frac{\alpha_{0}^{2}+\alpha_{1}^{2}}{1-\alpha_{0}^{2}}.$
(82)

For the second term in Eq. (81), using Eqs. (63) and (64) we get

$\displaystyle\sum_{j,r_{1},r_{2}}\mathcal{T}_{j,\mathbf{0},s,r_{1}}\mathcal{T}_{j,\mathbf{0},s,r_{2}}^{*}\left(\frac{\partial\tilde{\mathcal{G}}_{\mathbf{k}+\mathbf{q}_{j},r_{1},r_{2}}}{\partial
k}\right)_{0}$
$\displaystyle=\frac{1}{\left(1-\alpha_{0}^{2}\right)^{2}}\sum_{j,r}\frac{\left|\mathcal{T}_{j,\mathbf{0},s,r}\right|^{2}}{\left(vk_{\theta}\right)^{2}}\left[\left(\frac{\partial\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}{\partial
k}\right)_{0}+\left(\frac{\partial\Sigma^{\prime\prime}_{\mathbf{k}+\mathbf{q}_{j},r,r}}{\partial
k}\right)_{0}\right]$
$\displaystyle\qquad-\frac{1}{1-\alpha_{0}^{2}}\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{0},s,r}\mathcal{T}_{j,\mathbf{0},s,-r}^{*}}{\left(vk_{\theta}\right)^{2}}\left(\frac{\partial\Sigma^{\prime\prime}_{\mathbf{k}+\mathbf{q}_{j},r,-r}}{\partial
k}\right)_{0}.$ (83)

The first summation with the term proportional to
$\left(\partial\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}/\partial k\right)_{0}$
is given by Eq. (72). For the other terms, a long but straightforward
calculation using Eq. (65) yields

$\displaystyle\sum_{j,r}\frac{\left|\mathcal{T}_{j,\mathbf{0},s,r}\right|^{2}}{\left(vk_{\theta}\right)^{2}}\left(\frac{\partial\Sigma^{\prime\prime}_{\mathbf{k}+\mathbf{q}_{j},r,r}}{\partial
k}\right)_{0}$
$\displaystyle=\frac{1}{2}sv\left[\left(\alpha_{0}^{2}-\alpha_{1}^{2}\right)^{2}+8\alpha_{0}^{2}\alpha_{1}^{2}\right],$
(84)
$\displaystyle\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{0},s,r}\mathcal{T}_{j,\mathbf{0},s,-r}^{*}}{\left(vk_{\theta}\right)^{2}}\left(\frac{\partial\Sigma^{\prime\prime}_{\mathbf{k}+\mathbf{q}_{j},r,-r}}{\partial
k}\right)_{0}$
$\displaystyle=-\frac{1}{2}sv\frac{\left(\alpha_{1}^{2}-2\alpha_{0}^{2}\right)\left(\alpha_{0}^{2}+\alpha_{1}^{2}\right)}{1-\alpha_{0}^{2}}.$
(85)

Combining the previous results and simplifying, for the $k$-derivative we have

$\left(\frac{\partial\Sigma}{\partial
k}\right)_{0}=-sv\frac{3\alpha_{1}^{2}-\alpha_{0}^{4}-\alpha_{1}^{4}-4\alpha_{0}^{2}\alpha_{1}^{2}}{\left(1-\alpha_{0}^{2}\right)^{2}}.$
(86)

For the other derivative over $\omega$, proceeding as with the momentum
derivatives we get

$\left(\frac{\partial\Sigma}{\partial\omega}\right)_{0}=-\frac{3\left(\alpha_{0}^{2}+\alpha_{1}^{2}\right)+2\left(\alpha_{0}^{4}+\alpha_{1}^{4}+4\alpha_{1}^{2}\alpha_{0}^{2}\right)}{\left(1-\alpha_{0}^{2}\right)^{2}}.$
(87)

Finally, replacing in Eq. (66) yields

$\frac{v^{\star}}{v}=\frac{\left(1-\alpha_{0}^{2}\right)^{2}-3\alpha_{1}^{2}+\alpha_{0}^{4}+\alpha_{1}^{4}+4\alpha_{0}^{2}\alpha_{1}^{2}}{\left(1-\alpha_{0}^{2}\right)^{2}+3\left(\alpha_{0}^{2}+\alpha_{1}^{2}\right)+2\left(\alpha_{0}^{4}+\alpha_{1}^{4}+4\alpha_{0}^{2}\alpha_{1}^{2}\right)}.$
(88)

In the chiral case $\alpha_{0}=0$, the above reduces to
$v^{\star}/v=\left(1-3\alpha_{1}^{2}+\alpha_{1}^{4}\right)/\left(1+3\alpha_{1}^{2}+2\alpha_{1}^{4}\right)$,
in agreement with Ref. [48] up to $\sim\alpha_{1}^{4}$.

##  Appendix E First-shell velocity $v^{\star}$ as function of $\phi_{j}$ and
$\mathbf{q}_{j}$

Here we give details of the calculation of $v^{\star}$ at the first shell
using Eq. (16), but keeping explicitly the dependence on the phases $\phi_{j}$
and momentum transfer vectors $\mathbf{q}_{j}$ as they enter in the moiré
potential given by Eq. (5); cf. Sec. III.4. As before, we focus on the
perturbation on the top layer, and for simplicity note throughout
$\Sigma_{+,\mathbf{k},s}^{\left(1\right)}\rightarrow\Sigma_{\mathbf{k},s}$. We
will, in particular, consider the chiral limit $\alpha_{0}=0$, which
simplifies the calculation because then the Dirac points are never shifted in
energy [cf. Eq. (28)], and thus $v^{\star}$ can be obtained, as before, by
simply evaluating the derivatives of the self-energy at $k=\omega=0$.

The calculation of $v^{\star}$ can be done by following the same procedure as
in the particular TBG case in which
$\theta_{\mathbf{q}_{j}}-\phi_{j}=90^{\circ}$, where $\theta_{\mathbf{q}_{j}}$
is the angle of the momentum transfer vectors $\mathbf{q}_{j}$ (see Fig. 1).
However, since for $\theta_{\mathbf{q}_{j}}-\phi_{j}\neq 90^{\circ}$ the
derivatives of the interband contribution to the self-energy do not vanish,
one must be careful with the divergences of the last term in Eq. (18) when
evaluated at $k=0=\omega$. One way to circumvent that is to rather start from
Eq. (19),

$\omega-\epsilon_{\mathbf{k},s}-\left(\Sigma_{\mathbf{k},s}^{+}+\Sigma_{\mathbf{k},s}^{\mathrm{-}}\right)=0,$
(89)

where we have explicitly separated the intraband and interband contributions
to the self-energy, Eqs. (17) and (18), and to simplify the notation we have
noted

$\Sigma_{\mathbf{k},s}^{+}=\Sigma_{\mathbf{k},s}^{\mathrm{intra}},\quad\Sigma_{\mathbf{k},s}^{-}=\Sigma_{\mathbf{k},s}^{\mathrm{inter}}.$
(90)

Now we multiply Eq. (89) by the last factor in Eq. (18), so that

$\left[\omega-\epsilon_{\mathbf{k},s}-\left(\Sigma_{\mathbf{k},s}^{+}+\Sigma_{\mathbf{k},s}^{\mathrm{-}}\right)\right]\left(\omega-\epsilon_{\mathbf{k},-s}-\Sigma_{\mathbf{k},-s}^{+}\right)=0.$
(91)

Using $\epsilon_{\mathbf{k},-s}=-\epsilon_{\mathbf{k},s}$ we get

$\omega^{2}-\epsilon_{\mathbf{k},s}^{2}-\left(\omega-\epsilon_{\mathbf{k},s}\right)\Sigma_{\mathbf{k},-s}^{+}-\left(\omega+\epsilon_{\mathbf{k},s}\right)\Sigma_{\mathbf{k},s}^{+}+\Sigma_{\mathbf{k},s}^{+}\Sigma_{\mathbf{k},-s}^{+}-\Sigma_{\mathbf{k},s}^{-}\left(\omega+\epsilon_{\mathbf{k},s}-\Sigma_{\mathbf{k},-s}^{+}\right)=0.$
(92)

The last factor
$\Sigma_{\mathbf{k},s}^{-}\left(\omega+\epsilon_{\mathbf{k},s}-\Sigma_{\mathbf{k},-s}^{+}\right)$
gets effectively rid of the divergences of $\Sigma_{\mathbf{k},s}^{-}$ at
$k=0=\omega$. It is convenient to define

$\tilde{\Sigma}_{\mathbf{k},s}=\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}},$
(93)

so that
$\Sigma_{\mathbf{k},s}^{-}\left(\omega+\epsilon_{\mathbf{k},s}-\Sigma_{\mathbf{k},-s}^{+}\right)=\tilde{\Sigma}_{\mathbf{k},s}\tilde{\Sigma}_{\mathbf{k},s}^{*}$,
cf. Eq. (18). To obtain the quasiparticle velocity
$v^{\star}=\left(\partial\omega/\partial k\right)_{0}$ at the Dirac point, we
now derive Eq. (92) twice in momentum, evaluate all at $k=0=\omega$, and use
the chiral limit results

$\displaystyle\left(\sum_{j,r}\frac{\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}$
$\displaystyle=\left(\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}=0,$
(94) $\displaystyle\left(\frac{\partial}{\partial
k}\sum_{j,r}\frac{\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}$
$\displaystyle=3sv\alpha_{1}^{2}\cos\left(2\varphi\right),\qquad\left(\frac{\partial}{\partial
k}\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}=i3sv\alpha_{1}^{2}\sin\left(2\varphi\right),$
(95)
$\displaystyle\left(\frac{\partial}{\partial\omega}\sum_{j,r}\frac{\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}$
$\displaystyle=-3\alpha_{1}^{2},\qquad\left(\frac{\partial}{\partial\omega}\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}=0,$
(96) $\displaystyle\left(\frac{\partial^{2}}{\partial
k^{2}}\sum_{j,r}\frac{\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}$
$\displaystyle=\left(\frac{\partial^{2}}{\partial\omega^{2}}\sum_{j,r}\frac{\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}=\left(\frac{\partial^{2}}{\partial
k\partial\omega}\sum_{j,r}\frac{\left|\mathcal{T}_{j,\mathbf{k},s,r}\right|^{2}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}=0,$
(97) $\displaystyle\left(\frac{\partial^{2}}{\partial
k^{2}}\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}$
$\displaystyle=\left(\frac{\partial^{2}}{\partial\omega^{2}}\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}=\left(\frac{\partial^{2}}{\partial
k\partial\omega}\sum_{j,r}\frac{\mathcal{T}_{j,\mathbf{k},s,r}\mathcal{T}_{j,\mathbf{k},-s,r}^{*}}{\omega-\epsilon_{\mathbf{k}+\mathbf{q}_{j},r}}\right)_{0}=0,$
(98)

where $\varphi=\theta_{\mathbf{q}_{j}}-\phi_{j}$ (which is $j$-independent).
The last two rows imply, by Eq. (89), that also
$\left(\partial^{2}\omega/\partial k^{2}\right)_{0}=0$. Then we get

$\displaystyle 0$
$\displaystyle=\left(v^{\star}\right)^{2}\left(1-\partial_{\omega}^{0}\Sigma_{\mathbf{k},s}^{+}-\partial_{\omega}^{0}\Sigma_{\mathbf{k},-s}^{+}+\partial_{\omega}^{0}\Sigma_{\mathbf{k},s}^{+}\partial_{\omega}^{0}\Sigma_{\mathbf{k},-s}^{\mathrm{+}}\right)$
$\displaystyle\quad-v^{\star}\left(\partial_{k}^{0}\Sigma_{\mathbf{k},s}^{\mathrm{+}}+\partial_{k}^{0}\Sigma_{\mathbf{k},-s}^{\mathrm{+}}+v\partial_{\omega}^{0}\Sigma_{\mathbf{k},s}^{+}-v\partial_{\omega}^{0}\Sigma_{\mathbf{k},-s}^{+}-\partial_{k}^{0}\Sigma_{\mathbf{k},s}^{\mathrm{+}}\partial_{\omega}^{0}\Sigma_{\mathbf{k},-s}^{\mathrm{+}}-\partial_{k}^{0}\Sigma_{\mathbf{k},-s}^{\mathrm{+}}\partial_{\omega}^{0}\Sigma_{\mathbf{k},s}^{\mathrm{+}}\right)$
$\displaystyle\quad-v^{2}-v\left(\partial_{k}^{0}\Sigma_{\mathbf{k},s}^{\mathrm{+}}-\partial_{k}^{0}\Sigma_{\mathbf{k},-s}^{+}\right)+\partial_{k}^{0}\Sigma_{\mathbf{k},s}^{+}\partial_{k}^{0}\Sigma_{\mathbf{k},-s}^{+}-\partial_{k}^{0}\tilde{\Sigma}_{\mathbf{k},s}\partial_{k}^{0}\tilde{\Sigma}_{\mathbf{k},s}^{*},$
(99)

where, to condense the notation, we have noted
$\partial_{k}^{0}\rightarrow\left(\partial/\partial k\right)_{0}$ and
$\partial_{\omega}^{0}\rightarrow\left(\partial/\partial\omega\right)_{0}$.
Replacing the derivatives yields

$0=\left(v^{\star}\right)^{2}\left(1+6\alpha_{1}^{2}+9\alpha_{1}^{4}\right)-v^{2}-sv\left[6sv\alpha_{1}^{2}\cos\left(2\varphi\right)\right]-\left[3v\alpha_{1}^{2}\cos\left(2\varphi\right)\right]^{2}-\left[3v\alpha_{1}^{2}\sin\left(2\varphi\right)\right]^{2},$
(100)

which leads to Eq. (30) in the main text. Note that the interband contribution
comes only from the last term in Eq. (99), which vanishes when
$\varphi=90^{\circ}$. This means that moiré-induced interband scatterings,
within each layer, do not contribute to $v^{\star}$ only for the specific
orientation of the momentum transfer vectors in TBG; otherwise they influence
$v^{\star}$ and can actually prevent it for vanishing.