# Supplemental Material for “Twisted Trilayer Graphene: a Precisely Tunable
Platform for Correlated Electrons”

Ziyan Zhu  Department of Physics, Harvard University, Cambridge, Massachusetts
02138, USA  Stephen Carr  Department of Physics, Harvard University,
Cambridge, Massachusetts 02138, USA  Daniel Massatt  Department of Statistics,
The University of Chicago, Chicago, Illinois 60637, USA  Mitchell Luskin
School of Mathematics, University of Minnesota - Twin Cities, Minneapolis,
Minnesota 55455, USA  Efthimios Kaxiras  Department of Physics, Harvard
University, Cambridge, Massachusetts 02138, USA  John A. Paulson School of
Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts
02138, USA

The Supplemental Material includes four sections. In Section I, we discuss the
geometry of the twisted trilayer graphene (tTLG) and calculate the higher-
order moiré of moiré lengths. In Section II, we present a detailed derivation
of the momentum-space model and test its convergence. In Section III, we
derive analytically the magic angles in tTLG. Finally, in Section IV, we
compare and contrast our results with a simplified model proposed by Mora  _et
al._ (2019) as well as results obtained with a full model without the low-
energy expansion proposed by Amorim and Castro (2018).

##  I Calculation of moiré of moiré lengths

The atomic and reciprocal space geometry of tTLG with two independent twist
angles are shown in Fig. S1(a). The monolayer lattice vectors are defined as
the column vectors of the following matrix:

$\displaystyle A_{0}=a_{G}\begin{bmatrix}1&1/2\\\
0&\sqrt{3}/2\end{bmatrix}=\begin{bmatrix}\bm{a}_{1}&\bm{a}_{2}\end{bmatrix},$
(S1)

where $a_{G}=2.4768\,$Å is the graphene lattice constant (as obtained from
DFT). The $l$-th layer will be referred to as L$\ell$. We assume that L2 is
unrotated, with L1 rotated clockwise by $\theta_{12}$ and L3 rotated
counterclockwise by $\theta_{23}$. Defining the counterclockwise rotation
matrix

$\displaystyle\mathcal{R}(\theta)=\begin{bmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{bmatrix},$ (S2)

the lattice vectors of the three layers can be written as
$A_{1}=\mathcal{R}(-\theta_{12})A_{0}$, $A_{2}=A_{0}$, and
$A_{3}=\mathcal{R}(\theta_{23})A_{0}$ respectively, with the column vectors
denoted as $\bm{a}^{(\ell)}_{i}$ for $\ell=1,2,3$ and $i=1,2$. The monolayer
reciprocal lattice vectors are given by the columns of $G_{\ell}=2\pi
A_{\ell}^{-T}$. For example, the reciprocal lattice vectors of L2 are
$\bm{b}^{(2)}_{1}=2\pi/a_{G}(1,-\sqrt{3}/3)$ and
$\bm{b}^{(2)}_{2}=2\pi/a_{G}(0,2/\sqrt{3})$. The K point of L2 is given as
$K_{\mathrm{L}{2}}=(2\bm{b}^{(2)}+\bm{b}^{(1)})/3=(4\pi/(3a_{G}),0).$ The
reciprocal lattice vectors of layers L1 and L3 can be obtained by acting
$\mathcal{R}(-\theta_{12})$ and $\mathcal{R}(\theta_{23})$ on $\bm{b}^{(i)}$
for $i=1,2$. We also denote the monolayer unit cell of layer $\ell$ to be
$\Gamma^{(\ell)}$ and the reciprocal space to be $\Gamma^{(\ell)*}$.

Figure S1: Lattice geometries of the tTLG system. (a) The twisted trilayer
graphene system in real space (left) and momentum space with their original
monolayer reciprocal lattice vectors (right). (b) Bilayer moiré Brillouin zone
for $\theta_{12}=\theta_{23}$ with high symmetry points

The twisted trilayer system exhibits higher order moiré of moiré patterns due
to the interference between the two bilayer moiré patterns. To the lowest
order, the bilayer moiré length between layers $i$ and $j$ is given by
$a_{M}^{ij}=a_{G}/\sin\theta_{ij}$. We denote the bilayer moiré superlattice
between layers $ij$ to be $\Gamma^{ij}$, spanned by the column vectors of
matrix $A_{ij}=(\bm{a}^{(ij)}_{1}\quad\bm{a}^{(ij)}_{2})$. The bilayer moiré
Brillouin zone between layers $i$ and $j$ are given by the column vectors of
$G_{ij}=G_{j}-G_{i}=2\pi(A_{j}^{-T}-A_{i}^{-T}).$ The lattice vector of the
moiré superlattice is the column vectors of $A_{ij}=2\pi G_{ij}^{-T}$. After
some algebra, we obtain the lattice vectors for the two bilayer supercells,
$\bm{a}_{i}^{12},\bm{a}_{i}^{23}$ for $i=1,2$:

$\displaystyle\bm{a}^{12}_{1}$
$\displaystyle=\frac{a_{G}}{2(1-\cos\theta_{12})}\left[(1-\cos\theta_{12})\hat{x}-\sin\theta_{12}\hat{y}\right],$
$\displaystyle\bm{a}^{12}_{2}$
$\displaystyle=\frac{a_{G}}{2(1-\cos\theta_{12})}\Bigg{[}\left(\frac{1-\cos\theta_{12}}{2}+\frac{\sqrt{3}}{2}\sin\theta_{12}\right)\hat{x}+\left(-\frac{1}{2}\sin\theta_{12}+\frac{\sqrt{3}}{2}(1-\cos\theta_{12})\right)\hat{y}\Bigg{]},$
$\displaystyle\bm{a}^{23}_{1}$
$\displaystyle=\frac{a_{G}}{2(1-\cos\theta_{23})}\left[(1-\cos\theta_{23})\hat{x}+\sin\theta_{23}\hat{y}\right],$
$\displaystyle\bm{a}^{23}_{2}$
$\displaystyle=\frac{a_{G}}{2(1-\cos\theta_{23})}\Bigg{[}\left(\frac{1-\cos\theta_{23}}{2}-\frac{\sqrt{3}}{2}\sin\theta_{23}\right)\hat{x}+\left(\frac{1}{2}\sin\theta_{23}+\frac{\sqrt{3}}{2}(1-\cos\theta_{23})\right)\hat{y}\Bigg{]}.$
(S3)

Note that there is a small twist angle between the two bilayer moiré
superlattices. Moreover, for general twist angles
$\theta_{12}\neq\theta_{23}$, the two bilayer moiré cells have a lattice
mismatch. The twist angle $\phi$ and the lattice mismatch $\delta$ between the
two bilayer moiré patterns give rise to higher-order moiré of moiré lengths.
The primitive reciprocal lattice vectors of a given harmonic $(m,n)$ are given
as the column vectors of $G_{mn}^{\mathrm{H}}=mG_{12}-nG_{23}$. Inverting
$G_{mn}^{\mathrm{H}}$, we obtain the moiré of moiré supercell in real space
$A_{mn}^{\mathrm{H}}=\frac{1}{2\pi}(G_{mn}^{\mathrm{H}})^{-T}$. The norm of
the column vectors are the moiré of moiré lengths, denoted as
$\lambda_{mn}^{\mathrm{H}}$. For $m=n=1$, the moiré of moiré length
$\lambda_{11}^{\mathrm{H}}$ is explicitly given as

$\lambda_{11}^{\mathrm{H}}=\frac{(1+\delta)a_{M}^{12}}{\sqrt{2(1+\delta)(1-\cos\phi)+\delta^{2}}},$
(S4)

where
$\phi=\cos^{-1}\left(\frac{\bm{a}^{(12)}_{1}\cdot\bm{a}^{23}_{1}}{|\bm{a}^{(12)}_{1}||\bm{a}^{(23)}_{1}|}\right)$
is the twist angle between the bilayer moiré supercells and
$\delta=\frac{\sin\theta_{23}}{\sin\theta_{12}}-1$ is the lattice mismatch
between the two bilayer moiré supercells such that
$a^{23}_{M}=(1+\delta)a^{12}_{M}$. Equation (S4) agrees with the first-order
approximation for the moiré length for a twisted bilayer with a lattice
mismatch Yankowitz  _et al._ (2012), with the lattice constant being the
bilayer moiré length between L1 and L2.

A dominant moiré of moiré length does not necessarily exist nor evolve
smoothly under the continuous change of the twist angle. To see this, we will
consider different harmonics of the higher-order moiré patterns. To find the
dominant harmonic for an arbitrary pair of twist angles, we calculate
$A_{mn}^{\mathrm{H}}$ for $|m|,|n|\leq 15$ numerically and find the $(m,n)$
such that the norm of $G_{mn}^{\mathrm{H}}$ is the smallest, or, equivalently,
that the moiré of moiré length $\lambda_{mn}$ is largest. We are neglecting
the cases where higher order harmonics dominate, such as the cases where
$\theta_{12}$ and $\theta_{23}$ are different by more than a factor of 15. In
those cases, the two bilayers moiré supercells have very different sizes and
become essentially decoupled, which is not the focus of our study. Figure
S2(a) shows the moiré of moiré harmonics for varying $\theta_{12}$ at a fixed
$\theta_{23}=2.8^{\circ}$, indicating the non-smooth dependence of the
dominant moiré of moiré length on the twist angle. In tTLG, there exists a
supercell approximation when there is a clear dominant harmonic, that is when
$\theta_{12}\approx N\theta_{23}$ or $\theta_{23}\approx N\theta_{12}$ for
$N\in\mathbb{Z}$. For example, at
$\theta_{12}=2.6^{\circ},\theta_{23}=2.8^{\circ}$, the dominant harmonic is
$(1,1)$ [see Fig. S2(b)] and at
$\theta_{12}=1.35^{\circ},\theta_{23}=2.8^{\circ}$, the dominant harmonic is
$(2,1)$ [see Fig. S2(d)]. However, there are cases where there is no clear
dominant moiré of moiré. For example, in Fig. S2(c), it is difficult to
visually discern a large repeating pattern and the estimated moiré of moiré
lattice vectors fail to capture the relevant length scale. This is because
near $\theta_{12}=1.8^{\circ}$, many harmonics, such as (3,2), (5,4), and
(5,3), all have comparable lengths [see Fig. S2(a) the corresponding point].

Figure S2: (a) Moiré of moiré of lengths $\lambda_{mn}^{\mathrm{H}}$ as a
function of $\theta_{12}$ for $\theta_{23}=2.8^{\circ}$. Each color
corresponds to a different set of $(m,n)$. The thick black line indicates the
dominant length. (b)-(d) Examples moiré of moiré geometries, corresponding to
the red crosses in (a). Top: red and blue scattered points are the lattice
points of the bilayer moiré supercells between L1-L2 and L2-L3 respectively.
Due to the different twist angles, the moiré lattice vectors are slightly
rotated and have different lattice constants. Black vectors indicate estimated
dominant moiré of moiré supercell lattice vectors. A blowup of the small boxed
area is shown below, with points representing the atomic positions of each
monolayer graphene, for L1-L2 on the left half and for L2-L3 on the right
half. Red and blue vectors are the bilayer moiré lattice vectors of L1-L2 and
L2-L3 respectively.

##  II Momentum-space model

In this section, we offer a detailed derivation of our momentum-space model
and density of states formalism, and study the convergence as a function of
the momentum space cutoff radius.

###  II.1 Detailed derivation of momentum-space model

To model the electronic structure of the tTLG system, we start from a tight-
binding approximation for each individual layer; we take into account the
interlayer hopping in a transverse tight-binding approximation between nearest
neighbors. We start by writing the Hamiltonian for the trilayer as a sum of
the following terms

$H=\sum_{l=1}^{3}H^{\ell}+\sum_{l=1,2}\left(V^{\ell,\ell+1}+V^{\ell+1,\ell}\right),$
(S5)

where $H^{\ell}$ is the Hamiltonian for the $\ell$-th layer and $V^{ij}$
describes the interlayer hopping. For simplicity, we only consider the
interlayer couplings between adjacent layers. DFT calculations predict that
the interlayer coupling between L1 and L3 is roughly 10 times smaller than the
coupling between adjacent layers [e.g., between L1 and L2] Carr  _et al._
(2020). In a second quantized notation, $H^{\ell}$ can be written as

$\displaystyle
H^{\ell}=-t\sum_{\bm{R}^{(\ell)}}c_{\ell,A}^{\dagger}(\bm{R}^{(\ell)})[c_{\ell,B}(\bm{R}^{(\ell)})+c_{\ell,B}(\bm{R}^{(\ell)}-\bm{a}^{(\ell)}_{1})+c_{\ell,B}(\bm{R}^{(\ell)}-\bm{a}^{(\ell)}_{2})]+h.c.,$
(S6)

where $c_{\ell,\alpha}^{\dagger}$ and $c_{\ell,\alpha}$ are the creation and
annihilation fermionic operators of the orbital $\alpha$ in layer $l$,
$\bm{a}^{(\ell)}_{1,2}$ are the lattice vectors of layer $l$, and $t$ is the
hopping parameter between nearest neighbors. As for the interlayer coupling,
we define the following overlap matrix element in the tight-binding basis

$t_{\alpha\beta}^{ij}(\bm{R}^{(i)},\bm{R}^{(j)})=\bra{i,\bm{R}^{(i)},\alpha}H\ket{j,\bm{R}^{(j)},\beta},$
(S7)

where $\alpha$ and $\beta$ denotes the sublattice degree of freedom. The
interlayer Hamiltonian in the second-quantized notation is

$V^{ij}=\sum_{\bm{R}^{(i)},\alpha,\bm{R}^{(j)},\beta}c_{i,\alpha}^{\dagger}(\bm{R}^{(i)})t_{\alpha\beta}^{ij}(\bm{R}^{(i)},\bm{R}^{(j)})c_{j,\beta}(\bm{R}^{(j)}).$
(S8)

We obtain the Hamiltonian in the momentum basis at a center site momentum
$\bm{k}$. Defining $\bm{Q}^{(\ell)}=\bm{k}^{(\ell)}+\bm{k}$ for
$\bm{k}^{(\ell)}\in\Gamma^{(\ell)*}$. We perform the Fourier transform as
follows

$\displaystyle
c^{\dagger}_{\ell,\alpha}(\bm{R}^{(\ell)})=\frac{1}{\sqrt{|\Gamma^{(\ell)*}|}}\int_{\Gamma^{(\ell)*}}\mathrm{d}\bm{k}^{(\ell)}e^{i\bm{Q}^{(\ell)}\cdot(\bm{R}^{(\ell)}+\bm{\tau}^{(\ell)}_{\alpha})}c^{\dagger}_{\ell,\bm{k}^{(\ell)},\alpha},$
$\displaystyle
c_{\ell,\alpha}(\bm{R}^{(\ell)})=\frac{1}{\sqrt{|\Gamma^{(\ell)*}|}}\int_{\Gamma^{(\ell)*}}\mathrm{d}\bm{k}^{(\ell)}e^{-i\bm{Q}^{(\ell)}\cdot(\bm{R}^{(\ell)}+\bm{\tau}^{(\ell)}_{\alpha})}c_{\ell,\bm{k}^{(\ell)},\alpha},$
(S9)

where the integral is over the Brillouin zone of the $\ell$-th layer,
$\bm{\tau}^{(\ell)}_{A}=\bm{0}$,
$\bm{\tau}^{(\ell)}_{B}=1/3(\bm{a}^{(\ell)}_{1}+\bm{a}^{(\ell)}_{2})$. The
inverse of the transform in Eq. (S9) is

$\displaystyle
c^{\dagger}_{\ell,\bm{k}^{(\ell)},\alpha}=\frac{1}{\sqrt{|\Gamma^{(\ell)*}|}}\sum_{\bm{R}^{(\ell)}}e^{-i\bm{Q}^{(\ell)}\cdot(\bm{R}^{(\ell)}+\bm{\tau}^{(\ell)}_{\alpha})}c^{\dagger}_{\ell,\alpha}(\bm{R}^{(\ell)}),$
$\displaystyle
c_{\ell,\bm{k}^{(\ell)},\alpha}=\frac{1}{\sqrt{|\Gamma^{(\ell)*}|}}\sum_{\bm{R}^{(\ell)}}e^{i\bm{Q}^{(\ell)}\cdot(\bm{R}^{(\ell)}+\bm{\tau}^{(\ell)}_{\alpha})}c_{\ell,\alpha}(\bm{R}^{(\ell)}),$
(S10)

where $|\Gamma^{(\ell)*}|$ is the area of the Brillouin zone in the $\ell$-th
layer. The intralayer Hamiltonian in the Bloch basis can now be written as
follows

$\displaystyle H^{\ell}$
$\displaystyle=-\frac{t}{|\Gamma^{(\ell)*}|}\sum_{\bm{R}^{(\ell)}}\int_{\Gamma^{(\ell)*}}\mathrm{d}\bm{k}^{(\ell)}\int_{\Gamma^{(\ell)*}}\mathrm{d}\bm{k}^{\prime(\ell)}e^{i(\bm{k}^{(\ell)}-\bm{k}^{\prime(\ell)})\cdot\bm{R}^{(\ell)}}\sum_{\bm{s}^{(\ell)}_{i}}e^{i\bm{k}\cdot\bm{s}^{(\ell)}_{i}}c^{\dagger}_{\ell,\bm{k},A}c_{\ell,\bm{k}^{\prime(\ell)},B}$
$\displaystyle=-t\int_{\Gamma^{(\ell)*}}\mathrm{d}\bm{k}^{(\ell)}\sum_{\bm{s}^{(\ell)}}e^{i\bm{Q}^{(\ell)}\cdot\bm{s}^{(\ell)}_{i}}c^{\dagger}_{\ell,\bm{k}^{(\ell)},A}c_{\ell,\bm{k}^{(\ell)},B},$
(S11)

where we use the Poisson summation formula,
$\sum_{\bm{R}^{(\ell)}}e^{i\bm{k}^{(\ell)}\cdot\bm{R}^{(\ell)}}=|\Gamma^{(\ell)*}|\sum_{\bm{G}^{(\ell)}}\delta_{\bm{k}^{(\ell)},\bm{G}^{(\ell)}}$.
We also define $\bm{s}^{(\ell)}_{i}$ to describe the nearest neighbor
separation between $A$ and $B$ sublattices, which are given as
$\bm{s}^{(\ell)}_{1}=1/3(\bm{a}^{(\ell)}_{1}+\bm{a}^{(\ell)}_{2}),\bm{s}^{(\ell)}_{2}=1/3(-2\bm{a}^{(\ell)}_{1}+\bm{a}^{(\ell)}_{2}),\bm{s}^{(\ell)}_{3}=1/3(\bm{a}^{(\ell)}_{1}-2\bm{a}^{(\ell)}_{2})$.
The intralayer Hamiltonian in the basis of $c_{\ell,\bm{k}^{(\ell)},\alpha}$
can then be written as

$H^{\ell}(\bm{Q}^{(\ell)})=-t\begin{bmatrix}0&f_{\ell}(\bm{Q}^{(\ell)})\\\
f_{\ell}^{*}(\bm{Q}^{(\ell)})&0\end{bmatrix},$ (S12)

where
$f_{\ell}(\bm{Q}^{(\ell)})=\sum_{\bm{s}^{(\ell)}_{i}}e^{i\bm{Q}^{(\ell)}\cdot\bm{s}^{(\ell)}_{i}}$.
The Hamiltonian is equivalent to the monolayer graphene tight-binding model at
a given momentum $\bm{Q}^{(\ell)}$ Castro Neto  _et al._ (2009). For the
intralayer Hamiltonian, there is no constraint on $\bm{Q}^{(\ell)}$.

Similarly, we write the interlayer Hamiltonian in the
$c^{\dagger}_{\ell,\bm{k}^{(\ell)},\alpha}$ basis

$\displaystyle
V^{ij}=\int_{\Gamma^{(i)*}}\mathrm{d}\bm{k}^{(i)}\int_{\Gamma^{(j)*}}\mathrm{d}\bm{k}^{(j)}\sum_{\alpha\beta}c^{\dagger}_{i,\bm{k}^{(i)},\alpha}T^{\alpha\beta}_{ij}(\bm{k}^{(i)},\bm{k}^{(j)})c_{j,\bm{k}^{(j)},\beta},$
(S13)

where we use Eq. (S9) and

$\displaystyle
T_{\alpha\beta}^{ij}(\bm{k}^{(i)},\bm{k}^{(j)})=\frac{1}{\sqrt{|\Gamma^{(i)*}||\Gamma^{(j)*}|}}\sum_{\bm{R}^{(i)},\bm{R}^{(j)}}e^{i\bm{Q}^{(i)}\cdot(\bm{R}^{(i)}+\bm{\tau}^{(i)}_{\alpha})}t_{\alpha\beta}^{ij}(\bm{R}^{(i)},\bm{R}^{(j)})e^{-i\bm{Q}^{(j)}\cdot(\bm{R}^{(j)}+\bm{\tau}^{(j)}_{\beta})}.$
(S14)

We now apply the two center approximation

$t_{\alpha\beta}^{ij}(\bm{R}^{(i)},\bm{R}^{(j)})=t_{\alpha\beta}^{ij}(\bm{R}^{(i)}+\bm{\tau}^{(i)}_{\alpha}-\bm{R}^{(j)}-\bm{\tau}^{(j)}_{\beta}),$
(S15)

and write the interlayer coupling in terms of a two-dimensional Fourier
Transform

$\displaystyle
t_{\alpha\beta}^{ij}(\bm{R}^{(i)},\bm{R}^{(j)})=t^{ij}_{\alpha\beta}(\bm{R}^{(i)}+\bm{\tau}^{(i)}_{\alpha}-\bm{R}^{(j)}-\bm{\tau}^{(j)}_{\beta})$
$\displaystyle=\int\,\frac{\mathrm{d}\bm{p}}{(2\pi)^{2}}e^{i\bm{p}\cdot(\bm{R}^{(i)}+\bm{\tau}^{(i)}_{\alpha}-\bm{R}^{(j)}-\bm{\tau}^{(j)}_{\beta})}\tilde{t}^{ij}_{\alpha\beta}(\bm{p}).$
(S16)

Plugging Eq. (S16) into Eq. (S14), the interlayer coupling matrix element in
momentum space is

$\displaystyle T_{\alpha\beta}^{ij}(\bm{k}^{(i)},\bm{k}^{(j)})$
$\displaystyle=\frac{1}{\sqrt{|\Gamma^{(i)*}||\Gamma^{(j)*}|}}\int\,\frac{\mathrm{d}\bm{p}}{(2\pi)^{2}}\sum_{\bm{R}^{(i)},\bm{R}^{(j)}}e^{i(\bm{Q}^{(i)}+\bm{p})\cdot(\bm{R}^{(i)}+\bm{\tau}^{(i)}_{\alpha})}t_{\alpha\beta}^{ij}(\bm{p})e^{-i(\bm{Q}^{(j)}+\bm{p})\cdot(\bm{R}^{(j)}+\bm{\tau}^{(j)}_{\beta})}$
$\displaystyle=\sqrt{|\Gamma^{(i)*}||\Gamma^{(j)*}|}\sum_{\bm{G}^{(i)},\bm{G}^{(j)}}\int\,\frac{\mathrm{d}\bm{p}}{(2\pi)^{2}}e^{i\bm{G}^{(i)}\cdot\bm{\tau}^{(i)}_{\alpha}}\tilde{t}^{ij}_{\alpha\beta}(\bm{p})e^{-i\bm{G}^{(j)}\cdot\bm{\tau}^{(i)}_{\beta}}\delta_{\bm{k}+\bm{k}^{(i)}-\bm{p},\bm{G}^{(i)}}\delta_{\bm{k}+\bm{k}^{(j)}-\bm{p},\bm{G}^{(j)}}$
$\displaystyle=\frac{1}{|\Gamma|}\sum_{\bm{G}^{(i)},\bm{G}^{(j)}}e^{i\bm{G}^{(i)}\cdot\bm{\tau}^{(i)}_{\alpha}}\tilde{t}^{ij}_{\alpha\beta}(\bm{k}^{(i)}+\bm{k}-\bm{G}^{(i)})e^{-i\bm{G}^{(j)}\cdot\bm{\tau}^{(j)}_{\beta}}\delta_{\bm{k}^{(i)}-\bm{G}^{(i)},\bm{k}^{(j)}-\bm{G}^{(j)}}.$
(S17)

In the last step, we use the Possion summation rule and
$|\Gamma^{(\ell)*}|=4\pi^{2}|\Gamma|^{-1}$, where $|\Gamma|$ is the monolayer
unit cell area. We have obtained the scattering selection rule
$\bm{k}^{(i)}-\bm{G}^{(i)}=\bm{k}^{(j)}-\bm{G}^{(j)}$ for $i=j\pm 1$, which
imposes the constraint on the values of allowed $\bm{k}^{(\ell)}$.

Combining the intralayer and interlayer terms, the Hamiltonian in the
$c^{\dagger}_{\ell,\bm{k}^{(\ell)},\alpha}$ basis can be represented as a
$3\times 3$ block given in Eq. (1) of the main text.

###  II.2 Low-energy limit

We can greatly simplify the model by taking the low-energy limit. Each
$H^{\ell}$ can be expanded around its Dirac point,
$\bm{k}^{(\ell)}=K_{\mathrm{L}\ell}+\bm{q}^{(\ell)}$, as a rotated Dirac
Hamlitonian $H_{D}^{\ell}(\bm{q})$ for
$\bm{q}=\bm{k}+\bm{k}^{(\ell)}-K_{\mathrm{L\ell}}$:

$\displaystyle H^{1}(\bm{k})\approx
H_{D}^{1}(\bm{q})=v_{F}\begin{bmatrix}0&e^{i\theta_{12}}q_{+}\\\
e^{-i\theta_{12}}q_{-}&0\end{bmatrix},$ $\displaystyle H^{2}(\bm{k})\approx
H_{D}^{2}(\bm{q})=v_{F}\begin{bmatrix}0&q_{+}\\\ q_{-}&0\end{bmatrix},$
$\displaystyle H^{3}(\bm{k})\approx
H_{D}^{3}(\bm{q})=v_{F}\begin{bmatrix}0&e^{-i\theta_{23}}q_{+}\\\
e^{i\theta_{23}}q_{-}&0\end{bmatrix},$ (S18)

where $q_{\pm}=q_{x}\pm iq_{y}$. For the interlayer coupling, we substitute
$\bm{k}^{(\ell)}=\bm{q}^{(\ell)}+K_{\mathrm{L}{\ell}}$ into Eq. (S17),

$\displaystyle
T_{\alpha\beta}^{ij}(\bm{q}^{(i)},\bm{q}^{(j)})=\frac{1}{|\Gamma|}\sum_{\bm{G}^{(i)},\bm{G}^{(j)}}e^{i\bm{G}^{(i)}\cdot\bm{\tau}^{(i)}_{\alpha}}\tilde{t}_{\alpha\beta}^{ij}(\bm{k}+K_{\mathrm{L}{i}}+\bm{q}^{(i)}+\bm{G}^{(i)})e^{-i\bm{G}^{(j)}\cdot\bm{\tau}^{(j)}_{\alpha}}\delta_{\bm{q}^{(i)}+K_{\mathrm{L}{i}}-\bm{G}^{(i)},\bm{q}^{(j)}+K_{\mathrm{L}{j}}-\bm{G}^{(j)}}.$
(S19)

For momenta near the Dirac point, since
$|\bm{q}^{(i)}|,|\bm{k}|\ll|K_{\mathrm{L}{i}}|,|\bm{G}^{(i)}|$ we can
approximate
$\tilde{t}_{\alpha\beta}^{ij}(\bm{k}+K_{\mathrm{L}{i}}+\bm{q}^{(i)}+\bm{G}^{(i)})\approx\tilde{t}_{\alpha\beta}^{ij}(K_{\mathrm{L}{i}}+\bm{G}^{(i)}).$
This approximation can lead to the suppression of particle-hole asymmetry in
the tight-binding model Carr  _et al._ (2019); Fang  _et al._ (2019). Due to
the rapid decay of the hopping parameter $\tilde{t}(\bm{p})$ as $\bm{p}$
increases Bistritzer and MacDonald (2011), we keep only the first shell in the
summation in Eq. (S19):

$\displaystyle
T_{\alpha\beta}^{ij}(\bm{q}^{(i)},\bm{q}^{(j)})=\sum_{n=1}^{3}T_{n,\alpha\beta}^{ij}\delta_{\bm{q}^{(i)}-\bm{q}^{(j)},-\bm{q}_{n}^{ij}},$
(S20)

where $\bm{q}_{1}^{ij}=K_{\mathrm{L}{i}}-K_{\mathrm{L}{j}}$,
$\bm{q}_{2}^{ij}=\mathcal{R}^{-1}(2\pi/3)\bm{q}_{1}^{ij}$, and
$\bm{q}_{3}^{ij}=\mathcal{R}(2\pi/3)\bm{q}_{1}^{ij}$ (see Fig. S1b). We
include out-of-plane relaxation by letting
$t^{ij}_{AA}=t^{ij}_{BB}=\omega_{0}=0.07\,\mathrm{eV}$ and
$t^{ij}_{AB}=t^{ij}_{BA}=\omega_{1}=0.11\,\mathrm{eV},$ which matches with the
interlayer coupling in Nam and Koshino (2017) and Carr  _et al._ (2019). In
matrix form,

$\displaystyle T^{ij}_{1}=\begin{bmatrix}\omega_{0}&\omega_{1}\\\
\omega_{1}&\omega_{0}\end{bmatrix},T^{ij}_{2}=\begin{bmatrix}\omega_{0}&\omega_{1}\bar{\phi}\\\
\omega_{1}\phi&\omega_{0}\end{bmatrix},T^{ij}_{3}=\begin{bmatrix}\omega_{0}&\omega_{1}\phi\\\
\omega_{1}\bar{\phi}&\omega_{0}\end{bmatrix},$ (S21)

where $\phi=\exp(i\frac{2\pi}{3})$, $\bar{\phi}=\exp(-i\frac{2\pi}{3})$.

Figure S3: Top: DOS obtained using different sizes of momentum-space basis,
$n_{k}$. Bottom: errors in the DOS corresponding to the vertical lines on the
top with the same color. The error is defined as the difference between the
DOS value at the given energy and the DOS at the largest cutoff shown on the
top. The dashed lines are guides to the eyes showing power law scaling of the
error as a function of $n_{k}$.

###  II.3 Density of States

The DOS at a given energy $\epsilon$, $\mathcal{D}(\epsilon)$, for an
incommensurate tight-binding model is defined as Carr  _et al._ (2017)

$\displaystyle\mathcal{D}(\epsilon)=\sum_{r}\frac{1}{N}\sum_{n=1}^{N}\delta(\epsilon-\epsilon_{n})|\psi_{n}(r)|^{2},$
(S22)

where the $r$ sum is over all real space lattice positions, $n$ is the band
index, and $\psi_{n}(r)$ is the corresponding eigenfunction. To obtain the DOS
numerically, we use a Gaussian function
$\phi_{\epsilon,\kappa}(x)=\frac{2\sqrt{\ln 2}}{\sqrt{\pi}\kappa}\exp[-4\ln
2\frac{(x-\epsilon)^{2}}{\kappa^{2}}]$ to approximate the $\delta$ function,
and $\kappa$ is the full-width-half-maximum of the Gaussian, which determines
the energy resolution of the DOS Massatt  _et al._ (2018). We can transform
the DOS equation to momentum space:

$\displaystyle\mathcal{D}(\epsilon)=\frac{\mathcal{N}}{2}\sum_{\alpha=A,B}\sum_{\ell=1,2}\int_{\Gamma^{(\ell,\ell+1)*}}\sum_{n}\phi_{\epsilon,\kappa}(\epsilon_{n,\bm{k}})|\psi_{n,\bm{k}}|^{2}\
\mathrm{d}\bm{k},$ (S23)

where $\mathcal{N}$ is a normalization constant, $\epsilon_{n,\bm{k}}$ is an
energy within the energy window
$[\epsilon-\Delta\epsilon/2,\epsilon+\Delta\epsilon/2),$ $\Delta\epsilon$ is
the energy interval, $\psi_{n,\bm{k}}$ and $\epsilon_{n,\bm{k}}$ is an eigen-
pair of the Hamiltonian $\mathcal{H}(\bm{k})$ in Eq. (1) of the main text
associated with the center site $\bm{k}$ and band $n$. The integral is
evaluated over the bilayer moiré Brillouin zone between layers $\ell$ and
$\ell+1$, $\Gamma^{(\ell,\ell+1)*}$, and we discretize
$\Gamma^{(\ell,\ell+1)*}$ using a $22\times 22$ grid to evaluate the integral.
We adapt $\kappa$ based on the area of the integration domain
$\Gamma^{(\ell,\ell+1)*}$ as $\theta_{\ell,\ell+1}$ changes.

In order to make a direct comparison between the DOS at different twist
angles, we need to properly normalize the DOS. For a given cutoff radius, we
first calculate the DOS of the intralayer Hamiltonian only, which reduces to
three independent copies of monolayer graphene. Near the charge-neutrality
point, the DOS per $\mathrm{eV}$ per $\mathrm{nm}^{2}$ is given by Castro Neto
_et al._ (2009)

$\mathcal{D}(\epsilon)=\frac{6}{\pi}\frac{|\epsilon|}{v_{F}^{2}},$ (S24)

where the prefactor includes a factor 3 from the number of layers as well as a
factor of 4 from spin and valley degeneracies. We then obtain a normalization
constant by fixing the prefactor to the expected slope given in Eq. (S24) and
use the same constant for the DOS of the full Hamiltonian.

###  II.4 Convergence

The incommensurability of the tTLG system leads to an infinite number of
coupled momenta within any finite cutoff radius. Due to the additional
constraints we impose on the magnitude of $\bm{G}^{(\ell)}$, we neglect
degrees of freedom that can contribute to the low energy states. As a result,
there is no guaranteed convergence. Figure S3(a)-(d) shows the DOS and the
corresponding errors for different numbers of momentum degrees of freedom for
tTLG with two different sets of twist angles. In both cases, as the cutoff
increases, the error does not decay significantly. Note that in the case of
$\theta_{12}=2^{\circ},\theta_{23}=2.4^{\circ}$, the drop in error is most
likely a numerical artifact and further increasing the cutoff will not likely
to reduce the error. However, the physically relevant features, such as the
magnitude of the DOS maximum and the positions of the VHS, are relatively
stable as the cutoff increases. In contrast, Fig. S3(e), (f) shows the fast
convergence of the DOS in tBLG as a function of cutoff radius. This is because
in tBLG, increasing the cutoff radius does not change the number of relevant
low-energy degrees of freedom. In this work, we choose a cutoff at the
4${}^{\textrm{th}}$ honeycomb shell (i.e., $k_{c}=4|\bm{b}^{(\ell)}|$
corresponding to $\sim 5\,600$ momenta). This choice was made by considering
both computational efficiency and the accuracy of physical properties of
interest.

##  III Effective Hamiltonian and renormalized Fermi velocity

We examine the limit in which the momentum-space is truncated at the first
honeycomb shell. The truncation gives rise to the following $14\times 14$
Hamiltonian:

$\mathcal{H}(\bm{q})=\begin{bmatrix}H^{1}_{D}(\bm{q}+\bm{q}^{12}_{1})&&&(T_{1}^{12})^{\dagger}&&&\\\
&H^{1}_{D}(\bm{q}+\bm{q}^{12}_{2})&&(T_{2}^{12})^{\dagger}&&&\\\
&&H^{1}_{D}(\bm{q}+\bm{q}^{12}_{3})&(T_{3}^{12})^{\dagger}&&&\\\
T_{1}^{12}&T_{2}^{12}&T_{3}^{12}&H^{2}_{D}(\bm{q})&T_{1}^{23}&T_{2}^{23}&T_{3}^{23}\\\
&&&(T_{1}^{23})^{\dagger}&H^{3}_{D}(\bm{q}+\bm{q}^{23}_{1})\\\
&&&(T_{2}^{23})^{\dagger}&&H^{3}_{D}(\bm{q}+\bm{q}^{23}_{2})\\\
&&&(T_{3}^{23})^{\dagger}&&&H^{3}_{D}(\bm{q}+\bm{q}^{23}_{3})\end{bmatrix}.$
(S25)

This Hamiltonian acts on seven two-component spinors
$\Psi=(\psi_{11},\psi_{12},\psi_{13},\psi_{20},\psi_{31},\psi_{32},\psi_{33})$,
where $\ell$ and $j$ in $\psi_{\ell j}$ denote the layer and the momentum
basis index respectively. Using this Hamiltonian, we can derive an expression
for the renormalized Fermi velocity $v_{F}^{*}$. We first define the
dimensionless quantities $\alpha_{12}=\omega/v_{F}k_{\theta_{12}}$ and
$\alpha_{23}=\omega/v_{F}k_{\theta_{23}}$, where
$k_{\theta_{ij}}=\frac{8\pi\sin{(\theta_{ij}/2)}}{3a_{G}}$. For simplicity, we
assume $\omega_{0}=\omega_{1}=\omega$ and neglect the angular dependence in
$H_{D}^{\ell}$ by letting $H_{D}^{\ell}$ be an unrotated Dirac Hamiltonian:
$H_{D}^{\ell}(\bm{q})=v_{F}\bm{\sigma}\cdot\bm{q}$, where
$\bm{\sigma}=(\sigma_{x},-\sigma_{y})$ is the Pauli matrix. The zero-energy
state of the Hamiltonian satisfies
$\mathcal{H}\Psi=\sum_{j=1}^{7}c_{j}\Psi_{j}=\bm{0}$, where $c_{j}$ is the
column vectors of $\mathcal{H}$, and $\Psi_{j}$ is the $j$-th component of the
spinor $\Psi$. Therefore, we obtain the following relation between components
of $\Psi$

$\Psi_{j}=-(H_{D}^{\ell})^{-1}T^{\dagger}\psi_{20},$ (S26)

where $j\neq 4$ ($\Psi_{j}$ is not a state on L2 or $\psi_{20}$). Using this,
the effective Hamiltonian to the leading order in $\bm{q}$ is

$\displaystyle\bra{\Psi}\mathcal{H}^{2}_{D}(\bm{q})\ket{\Psi}$
$\displaystyle=\frac{v_{F}}{1+6(\alpha_{12}^{2}+\alpha_{23}^{2})}\psi_{20}^{\dagger}\Bigg{\\{}\bm{\sigma}\cdot\bm{q}+\sum_{n=1}^{3}\Big{[}T_{n}^{12}(H_{D}^{1}(\bm{q}+\bm{q}_{n}^{12}))^{-1}(\bm{\sigma}\cdot\bm{q})(H_{D}^{1}(\bm{q}+\bm{q}_{n}^{12}))^{-1}T^{12\dagger}_{n}$
$\displaystyle\quad+T_{n}^{23}(H_{D}^{3}(\bm{q}+\bm{q}_{n}^{23}))^{-1}(\bm{\sigma}\cdot\bm{q})(H_{D}^{3}(\bm{q}+\bm{q}^{23}_{n}))^{-1}T^{23\dagger}_{n}\Big{]}\Bigg{\\}}\psi_{20}$
$\displaystyle=v_{F}^{*}\psi_{20}^{\dagger}\bm{\sigma}\cdot\bm{q}\psi_{20},$
(S27)

where the renormalized Fermi velocity $v^{\star}_{F}$ is

$\frac{v_{F}^{*}}{v_{F}}=\frac{1-3(\alpha_{12}^{2}+\alpha_{23}^{2})}{1+6(\alpha_{12}^{2}+\alpha_{23}^{2})}.$
(S28)

Figure S4(a) shows the $v_{F}^{*}$ to $v_{F}$ ratio as a function of
$\theta_{12}$ at a few values of $\theta_{23}$. As $\theta_{23}$ increases,
the $v_{F}^{*}/v_{F}$ ratio approaches the tBLG curve. Figure S4(b) shows
$v_{F}^{*}/v_{F}$ for equal twist angles, which shows that perturbation theory
predicts that $v_{F}^{*}$ can still go to zero at $1.72^{\circ}$. However, in
our numerical calculation using the full Hamiltonian, we do not observe a
complete flattening of bands at this twist angle.

Figure S4: (a) The ratio of renormalized Fermi velocity $v_{F}^{*}$ to the
monolayer Fermi velocity $v_{F}$ as a function of $\theta_{12}$ for given
values of $\theta_{23}$. Black dashed line shows the tBLG $v_{F}^{*}/v_{F}$
ratio. (b) $v_{F}^{*}/v_{F}$ ratio as a function of twist angle for
$\theta_{12}=\theta_{23}$. Figure S5: Comparison between the renormalized
Fermi velocity $v_{F}^{*}$ of the Hamiltonian in Eq. (S25) calculated
analytically (solid line) and numerically (dashed line).

Finally, we show that our assumption in the analytic calculation of an
unrotated Dirac Hamiltonian for the intralayer Hamiltonian and
$\omega_{0}=\omega_{1}$ does not significantly change the magic angle
estimate. Figure S5 compares the $v_{F}^{*}$ obtained analytically and
numerically and show that the two curves and the magic angle do not differ
significantly. In the numerical calculation, we diagonalize the $14\times 14$
Hamiltonian with rotated Dirac equation for the intralayer terms and
$\omega_{0}=0.07\,\mathrm{eV},\omega_{1}=0.11\,\mathrm{eV}$ for the interlayer
terms. At $\theta_{23}=2.5^{\circ}$, the magic angle obtained analytically and
numerically differ by 1.1%.

##  IV Comparison to other models

Figure S6: Comparison of the bilayer moiré Brillouin zone geometry between our
model and Mora  _et al._ (2019) model with $\theta_{12}=\theta_{23}$. Left:
two bilayer moiré Brillouin zones are misaligned by a small twist angle;
right: the two bilayer moiré Brillouin zones are approximated to be aligned.
Figure S7: The momentum degrees of freedom in the low-energy limit of the
simplified model for (a) $\theta_{12}=\theta_{23}$ and (b)
$2\theta_{12}=\theta_{23}$ with high symmetry points.

In this section, we compare our results to two other works Mora  _et al._
(2019); Amorim and Castro (2018). We first compare our results with the model
proposed by Mora et al. Mora  _et al._ (2019) and use it to gain further
insights into our findings. In this alternate model, a different momentum-
space basis is used by aligning the two bilayer moiré Brillouin zones [Fig.
S6]. This approximation ignores the incommensurability of the system, making a
two-dimensional momentum space crystal with the periodicity of the bilayer
moiré Brillouin zone. As a result, the problem’s complexity reduces to that of
a bilayer. Formally, the Hamiltonian can still be written as the $3\times 3$
block as in Eq. (1) in the main text, but the size of the basis is reduced to
be on the same order as tBLG. We implemented two cases: (1)
$\theta_{12}=\theta_{23}$ and (2) $2\theta_{12}=\theta_{23}$. Figure S7 shows
the momentum-space basis for these two cases. In case (2), the larger bilayer
Brillouin zone (L1-L2) is folded onto the smaller Brillouin zone (L2-L3) in
momentum space. This model essentially describes a system consisted of
$2\times 2$ L1-L2 moiré supercell and a L2-L3 moiré supercell. Figure S9 shows
a comparison between the DOS obtained from the two models. We keep the values
of $\omega_{0}$ and $\omega_{1}$ the same as our model and use the same
approach to normalize the DOS for a direct comparison. We cut off the basis at
the $4^{\mathrm{th}}$ shell and use a grid size $22\times 22$ for the density
of states. The Gaussian FWHM we use is 5 meV for $\theta<2^{\circ}$ and 8 meV
for $\theta\geq 2^{\circ}$, where $\theta$ is the twist angle that determines
the size of the Brillouin zone.

For $\theta_{12}=\theta_{23}$, Fig. S8(a) shows the DOS obtained with the
simplified model, which agrees qualitatively with the DOS from our model [Fig.
3(a) of the main text]. However, here the DOS has the sharpest peak between
$1.7^{\circ}$ and $1.8^{\circ}$, and at $2.1^{\circ}$ the VHS have a larger
width compared to our model. Figure S8(b)-(d) shows that the location of peaks
away from the CNP are also very different from our model.

Figure S8: DOS obtained with the Mora  _et al._ (2019) model DOS on a
logarithmic color scale at equal twist angles (same color scale as Fig. 3(a)
in the main text for a direct comparison) (b)-(d) DOS states along the black
dashed line in (a) for (b) $\theta=1.8^{\circ}$, (c) $\theta=2.1^{\circ}$, (d)
$\theta=2.5^{\circ}$, where the black solid lines are obtained using the Mora
_et al._ (2019) model, and the blue dashed lines are obtained using our full
model.

For $2\theta_{12}=\theta_{23}$, the two models predict similar trend for the
VHS evolution, and the simplified model makes the right prediction for the
magic angle. This is expected from perturbation theory, since the magic angle
condition does not rely on the existence of a moiré of moiré cell [as was
shown in Section III]. However, the magnitude of the DOS differs significantly
between the two models. This is because there are two flat bands near the CNP
in the simplified model, whereas in our full model, there is a large number of
nearly overlapping flat bands due to incommensurability [Fig. S10]. Figure S10
compares the band structure from our model and the simplified model. The two
band structures are qualitatively similar but our model shows a large number
of bands due to the lack of a periodic Brillouin zone. Furthermore, the
aligned-bilayer approximation will exclude correlated phases that depend on
band-hybridization or symmetries from the moiré of moiré length scale. Note
that we do not plot the relative layer weights (color) of the band structure
in the simplified model because of the way that the Brillouin zone is wrapped
– the L1 degrees of freedom are wrapped on top of the L3 degrees of freedom.
Therefore, the wavefunction weights from the two models are not directly
comparable for this particular high symmetry line cut.

Figure S9: DOS as a function of $\theta_{12}$ with $2\theta_{12}=\theta_{23}$
using (a) our full model and (b) the Mora  _et al._ (2019) model, both on a
logarithmic color scale. Figure S10: Comparison of band structures and DOS at
$\theta_{12}=1.4^{\circ},\theta_{23}=2.8^{\circ}$ between our model (top)
along the green dashed line in Fig. S1(b) and the Mora  _et al._ (2019) model
(bottom) along the black dashed line in Fig. S7(b). In (a), colors represent
the weight of the wavefunctions at the center site. Red, blue, and green
represents weights purely from L1, L2, L3 respectively, and colors in between
represent hybridization between different layers. A colormap is provided on
the top left corner. The DOS from the two models are shown on the same scale.

We can use these results to further support our argument of bilayer moiré
hybridization at equal twist angles. In this simplified model, sharpest VHS
occur between $1.7^{\circ}$ and $1.8^{\circ}$, which is in better agreement
with the magic angle prediction from perturbation theory. In our model, the
sharpest peak and the narrowest width occurs at a larger angle
($2.1^{\circ}$). If this phenomenon is caused by moiré hybridization, the
simplified model would not have it since it does not have the moiré of moiré
scale. Indeed, the DOS from the two models differ most significantly at
$2.1^{\circ}$ [see Fig. S8(c)].

As we argue in the main text, adding electrons from the CNP at a low carrier
concentration on the order of the tTLG moiré of moiré cells fills one flat
band near the CNP in Fig. S10(a) at a time. Injecting electrons at a carrier
concentration comparable to the bilayer moiré cell density would fill all
these flat bands near the CNP. The simplified model can again be used to
understand this argument. The model also predicts some band flattening at
certain twist angles, but there are only two flat bands near the CNP [Fig.
S10(c)]. Filling electrons to these two bands is equivalent to filling the
bilayer moiré cell, since their momentum-space basis has the periodicity of
bilayer moiré Brillouin zone and there is no moiré of moiré length in this
model. These two flat bands near the CNP can be qualitatively considered as
the limit where all the flat bands from our model overlap exactly on top of
each other. Therefore, in terms of filling the supercell, filling the two flat
bands from the simplified model is equivalent to filling all flat bands in the
full model.

In addition to its inability to make predictions about electronic behaviors at
the moiré of moiré scale, another major limitation of the model is its
difficulty to generalize to arbitrary twist angles. For each set of twist
angles on a different $(m,n)$ harmonic, it requires the derivation of a new
basis by folding the bilayer moiré Brillouin zone, while our model’s basis is
insensitive to the choice of angles and overcomes this limitation.

We can also use our model to study the case where L1 and L3 are twisted in the
same direction (when $\theta_{12}$ and $\theta_{23}$ take opposite signs).
This case has been studied theoretically by Amorim and Castro (2018) and its
spectral properties have been investigated experimentally by Zuo  _et al._
(2018). Unlike our model, Amorim and Castro (2018) does not take the low-
energy limit [see Section II.2]. Figure S11 shows the band structure and the
corresponding DOS of $\theta_{12}=-2.81^{\circ},\theta_{23}=2.1^{\circ}$
obtained with our model, which is the same case as Figs. 1(a) and 2 presented
in Amorim and Castro (2018). The results from the two models show an
agreement, with the same VHS positions. The difference in the band structure
can be most likely attributed to the different ways of truncating the
momentum-space bases between the two models.

Figure S11: (a) Band structure along a high symmetry line that connects the
Dirac points of the three layers and (b) density of states at
$\theta_{12}=-2.81^{\circ},\theta_{23}=2.1^{\circ}$. The colormap in (a) is
the same as in Fig. S10(a).

## References

  * Mora  _et al._ (2019) Christophe Mora, Nicolas Regnault, and B. Andrei Bernevig, “Flatbands and Perfect Metal in Trilayer Moiré Graphene,” Phys. Rev. Lett. 123, 026402 (2019).
  * Amorim and Castro (2018) B. Amorim and Eduardo V. Castro, “Electronic spectral properties of incommensurate twisted trilayer graphene,” arXiv e-prints , arXiv:1807.11909 (2018), arXiv:1807.11909 [cond-mat.mes-hall] .
  * Yankowitz  _et al._ (2012) Matthew Yankowitz, Jiamin Xue, Daniel Cormode, Javier D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, Pablo Jarillo-Herrero, Philippe Jacquod, and Brian J. Leroy, “Emergence of superlattice Dirac points in graphene on hexagonal boron nitride,” Nature Physics 8, 382–386 (2012).
  * Carr  _et al._ (2020) Stephen Carr, Chenyuan Li, Ziyan Zhu, Efthimios Kaxiras, Subir Sachdev, and Alexander Kruchkov, “Ultraheavy and ultrarelativistic dirac quasiparticles in sandwiched graphenes,” Nano Letters 20, 3030–3038 (2020).
  * 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–162 (2009).
  * Carr  _et al._ (2019) Stephen Carr, Shiang Fang, Ziyan Zhu, and Efthimios Kaxiras, “Exact continuum model for low-energy electronic states of twisted bilayer graphene,” Physical Review Research 1, 013001 (2019).
  * Fang  _et al._ (2019) Shiang Fang, Stephen Carr, Ziyan Zhu, Daniel Massatt, and Efthimios Kaxiras, “Angle-Dependent Ab initio Low-Energy Hamiltonians for a Relaxed Twisted Bilayer Graphene Heterostructure,” arXiv e-prints , arXiv:1908.00058 (2019), arXiv:1908.00058 [cond-mat.mes-hall] .
  * Bistritzer and MacDonald (2011) Rafi Bistritzer and Allan H. MacDonald, “Moiré bands in twisted double-layer graphene,” Proceedings of the National Academy of Science 108, 12233–12237 (2011).
  * Nam and Koshino (2017) Nguyen N. T. Nam and Mikito Koshino, “Lattice relaxation and energy band modulation in twisted bilayer graphene,” Phys. Rev. B 96, 075311 (2017).
  * Carr  _et al._ (2017) Stephen Carr, Daniel Massatt, Shiang Fang, Paul Cazeaux, Mitchell Luskin, and Efthimios Kaxiras, “Twistronics: Manipulating the electronic properties of two-dimensional layered structures through their twist angle,” Phys. Rev. B 95, 075420 (2017).
  * Massatt  _et al._ (2018) Daniel Massatt, Stephen Carr, Mitchell Luskin, and Christoph Ortner, “Incommensurate heterostructures in momentum space,” Multiscale Modeling & Simulation 16, 429–451 (2018).
  * Zuo  _et al._ (2018) Wei-Jie Zuo, Jia-Bin Qiao, Dong-Lin Ma, Long-Jing Yin, Gan Sun, Jun-Yang Zhang, Li-Yang Guan, and Lin He, “Scanning tunneling microscopy and spectroscopy of twisted trilayer graphene,” Physical Review B 97, 035440 (2018).