# Magic distances for flat bands in twisted bilayer graphene

Antonio Palamara<EMAIL_ADDRESS>Dipartimento di Fisica,
Università della Calabria, Via P. Bucci, Cubo 30C, I-87036 Rende (CS), Italy
INFN, sezione LNF, Gruppo collegato di Cosenza, Via P. Bucci, Cubo 31C,
I-87036 Rende (CS), Italy  Michele Pisarra<EMAIL_ADDRESS>Dipartimento di Fisica, Università della Calabria, Via P. Bucci, Cubo 30C,
I-87036 Rende (CS), Italy  INFN, sezione LNF, Gruppo collegato di Cosenza, Via
P. Bucci, Cubo 31C, I-87036 Rende (CS), Italy  Antonello Sindona
<EMAIL_ADDRESS>Dipartimento di Fisica, Università della
Calabria, Via P. Bucci, Cubo 30C, I-87036 Rende (CS), Italy  INFN, sezione
LNF, Gruppo collegato di Cosenza, Via P. Bucci, Cubo 31C, I-87036 Rende (CS),
Italy

###### Abstract

Twisted bilayer graphene is known to host isolated and relatively flat bands
near charge neutrality, when tuned to specific _magic_ angles. Nonetheless,
these rotational misalignments, lying below $1.1$ degrees, result in long-
period moiré crystals, whose anomalous electronic properties are hardly
accessible to reliable atomistic simulations. Here, we present a map of
differently stacked graphene sheets, at arbitrary rotation angles
corresponding to precise interplanar distances, into an equivalence class
represented by _magic-angle_ twisted bilayer graphene. We determine the
equivalence relation in the class within a continuum model, and extend its
definition to a tight-binding approach. Then, we use density functional theory
to suggest that the _magic-angle_ physics may be characterized by costly
computational strategies on a twisted bilayer geometry, with conveniently
large stacking angles. Our results may pave the way for an ab initio
characterization of the unconventional topological phases and related
excitations, associated with currently observed low-energy quasi-flat bands.

##  I Introduction

Twisted bilayer graphene (TBG) refers to a unique structure formed by stacking
two graphene layers on top of each other, at given distance and rotational
mismatch, which offers an intriguing mixture of electronic and geometric
effects. In particular, the TBG lattice develops large moiré patterns for some
_magic_ twist angles, exhibiting strongly-correlated phenomena, such as
unconventional Mott insulator, Chern insulator and superconducting phases [1,
2, 3, 4, 5, 6]. These features are associated with the appearance of nearly
flat bands with vanishing group velocity at the Fermi level, as predicted by
theory [7, 8, 9, 10], and subsequently observed in a number of seminal
experiments [11, 12, 13, 14, 15, 16].

The discovery of socalled _magic-angle_ (MA) TBG has spawned the field of
twisttronics [17, 18]. Furthermore, MATBG has widened the already rich field
of Van der Waals (VdW) heterostructures [19], by introducing a new degree of
freedom in the ability to assemble different two-dimensional (2D) materials
through VdW interactions, with on-demand electronic, mechanical, and thermal
properties.

Nonetheless, the misalignment angle $\theta$ between the graphene sheets is
not the only geometric parameter responsible for the formation of quasi-flat
bands. The interlayer distance $d$ also plays an equally crucial role in the
development of these non-dispersive electronic states. Indeed, quasi-flat
bands can also be achieved by applying vertical pressures, of the order of the
gigapascal, to the layers of a non-magic-angle TBG system [20, 21, 22, 23,
24]. In the simplest picture [7], the emergence of quasi-flat bands arises
from the interplay of two energy scales. One involves the electronic states
over the moiré supelattice sites, and is determined by $\theta$. The other
reflects the coupling between the lattice sites of the TBG layers, and depends
on $d$.

Unfortunately, the unique features of MATBG are commonly observed at very
small twist angles $\theta{\sim}1.1^{\circ}$, corresponding to supercell sizes
of ${\sim}10^{4}$ atoms, which turns into a formidable challenge for the
modeling community. On these premises, we here discuss the existence of an
infinite set of crystal (and arguably quasi-crystal) structures, determined by
different values of $\theta$ associated to precise values of $d$, which belong
to an equivalence class, with MATBG serving as a representative member. The
map between elements of the class may allow a thorough investigation of the
magic-angle quasi-flat-band physics, by focusing on the most convenient
geometry within the class accessible to atomistic simulation methods. Other
mappings are available that allow to effectively reduce the unit cell size of
a TBG moiré crystal [25, 26, 27, 28]. However, they are defined in specific
tight binding (TB) models by operating on its free parameters. Our approach,
on the other hand, stems from an experimentally observable quantity, i.e., the
interlayer distance or the twist angle, being thus independent of the
underlying theory.

The reminder of this paper is organized as follows. In Sec. II, we introduce
the equivalence class within a continuum model (CM) approach [7, 8, 29, 30,
31, 32]. In Sec. III, we extend the map to a TB formulation [33]. In Sec. IV,
we suggest how the equivalence relation can be used in density-functional
theory (DFT) to access the Kohn-Sham (KS) spectrum of TBG at the smallest
magic-angle. Finally, in Sec. V, we outline a possible application of our
results on the time-dependent (TD) DFT treatment of MATBG.

##  II The equivalence class of MATBG

The CM provides a simple way to study and predict the emergence of quasi-flat
bands in a TBG of stacking angle $\theta$. Such an approach is based on the
single-valley Hamiltonian

$\mathcal{H}=\begin{bmatrix}\mathrm{v}_{\textsc{f}}\bm{\sigma}_{\theta/2}\cdot\mathbf{p}&\mathcal{T}(\mathbf{r})\\\
\mathcal{T}(\mathbf{r})^{\dagger}&\mathrm{v}_{\textsc{f}}\bm{\sigma}_{-\theta/2}\cdot\mathbf{p}\end{bmatrix},$
(1)

for a Dirac electron at position $\mathbf{r}$, with linear momentum
$\mathbf{p}=-i{\hbar}\bm{\nabla}$ [7, 8, 29, 30, 31, 32].

$\mathcal{H}$ acts on four-component Dirac bispinors, related to the top and
bottom graphene layers of the system, being respectively counter-clockwise
rotated by $\theta/2$ and $-\theta/2$, along the out-of-plane axis $z$. The
diagonal blocks in $\mathcal{H}$ account for the isolated graphene sheets, of
lattice constant $a=2.46$ Å, lying parallel to the $xy$ plane at distance $d$
from each other. These terms involve the Cartesian Pauli matrices
$\sigma_{x}$, $\sigma_{y}$, $\sigma_{z}$, through the rotated vector operators

$\bm{\sigma}_{\pm\theta/2}=e^{-i(\pm\theta/4)\sigma_{z}}(\sigma_{x},\sigma_{y})e^{i(\pm\theta/4)\sigma_{z}}.$
(2)

${\mathrm{v}_{\textsc{f}}}$ denotes the group velocity at the corners
(inequivalent $K$ points) of the hexagonal first Brillouin zone
($1^{\mathrm{st}}$BZ) of graphene, where the Dirac-cone bands meet. The off-
diagonal elements of $\mathcal{H}$ express the interlayer coupling, dictated
by the position-dependent hopping operator $\mathcal{T}(\mathbf{r})$, along
with its hermitian conjugate $\mathcal{T}(\mathbf{r})^{\dagger}$.

By the spatial and time-reversal symmetry groups of TBG at low energies and
small misalignment angles [30], the interlayer potential can be expanded into
a Fourier series over the wave vectors that connect equivalent Dirac-cone
vertices between the reciprocal sublattices of the top and bottom layers. The
series is dominated by nearest-neighbor contributions from the rotated
$1^{\mathrm{st}}$BZ’s of the two reciprocal sublattices, characterized by the
center-to-corner distance $K_{\textsc{d}}=4\pi/3a$.

Figure 1:  (a) Rotated $1^{\mathrm{st}}$BZ’s of the graphene sublattices in
TBG. $K_{\pm\theta/2}$ denote the equivalent Dirac-cone vertices, connected by
wave vectors $\mathbf{q}_{1}$, $\mathbf{q}_{2}$ and $\mathbf{q}_{3}$.
$K_{\textsc{d}}$ is the distance of the $K_{\pm\theta/2}$ points from the
$1^{\mathrm{st}}$BZ-center, or $\Gamma$ point. $\mathbf{b}_{1}$ and
$\mathbf{b}_{2}$ label the basis vectors for the unrotated reciprocal lattice
of graphene. (b) Reciprocal space of TBG. $\mathbf{q}_{1}$, $\mathbf{q}_{2}$
and $\mathbf{q}_{3}$ generate a triangular arrangement of equivalent
$\tilde{K}$ points, which, together with their inequivalent counterparts, make
the periodically replicated $1^{\mathrm{st}}$BZ of the moiré, with center
$\tilde{\Gamma}$ and mid-point $\tilde{M}$.

Hereinafter, we use $\mathbf{a}_{1}{=}a(1,\sqrt{3})/2$,
$\mathbf{a}_{2}{=}a(-1,\sqrt{3})/2$ as basis of real-space lattice vectors,
and $\mathbf{b}_{1}{=}2\pi(1,1/\sqrt{3})/a$,
$\mathbf{b}_{2}{=}2\pi(-1,1/\sqrt{3})/a$ as basis of reciprocal-space lattice
vectors, for unrotated graphene. With this choice, reported in Fig. 1a, the
leading terms in the Fourier expansion of $\mathcal{T}(\mathbf{r})$ are
associated to the wave vectors $\mathbf{q}_{1}{=}k_{\theta}(0,-1)$,
$\mathbf{q}_{2}{=}k_{\theta}(\sqrt{3}/2,1/2)$ and
$\mathbf{q}_{3}{=}k_{\theta}(-\sqrt{3}/2,1/2)$, of equal modulus
$k_{\theta}{=}2{K_{\textsc{d}}}\sin({\theta}/{2})$, which generate the
reciprocal lattice of TBG, sketched in Fig. 1b. Accordingly,
$\mathcal{T}(\mathbf{r})$ is approximated to the three-element sum

$\mathcal{T}(\mathbf{r})=\sum^{3}_{l=1}T_{l}e^{i\mathbf{q}_{l}\cdot\mathbf{r}},$
(3)

whose Fourier coefficients

$\displaystyle T_{1}=$
$\displaystyle\omega_{\textsc{AA}}\sigma_{0}+\omega_{\textsc{AB}}\sigma_{x},$
$\displaystyle T_{2}=$
$\displaystyle\omega_{\textsc{AA}}\sigma_{0}-\frac{\omega_{\textsc{AB}}}{2}\left(\sigma_{x}-\sqrt{3}\sigma_{y}\right),$
$\displaystyle T_{3}=$
$\displaystyle\omega_{\textsc{AA}}\sigma_{0}-\frac{\omega_{\textsc{AB}}}{2}\left(\sigma_{x}+\sqrt{3}\sigma_{y}\right),$
(4)

depend on the Cartesian Pauli matrices, the $2{\times}2$ identity matrix
$\sigma_{0}$, plus two interlayer hopping energies, $\omega_{\textsc{AA}}$ and
$\omega_{\textsc{AB}}$, being respectively related to AA and AB (or BA)
stacking centers [34].

Working under the limiting condition
$\omega_{\textsc{AA}}{=}\omega_{\textsc{AB}}{=}\omega$, the interlayer Fourier
coefficients of Eq. (4) become $T_{l}{=}\omega\tilde{T}_{l}$, and the
reference Hamiltonian of Eq. (1) reduces to the Bistritzer-MacDonald (BM)
Hamiltonian [7]. The latter establishes a minimal model capable of predicting
the emergence of quasi-flat bands with vanishing group velocity at the Fermi
level.

We now rely on the dimensionless position operator
$\tilde{\mathbf{r}}=k_{\theta}\mathbf{r}$ and associated linear-momentum
operator $\tilde{\mathbf{p}}=-i\bm{\nabla}_{\tilde{\mathbf{r}}}$, which allow
the BM Hamiltonian to be written as

$\mathcal{H}=\begin{bmatrix}{\hbar}{\mathrm{v}_{\textsc{f}}}k_{\theta}\bm{\sigma}_{\theta/2}\cdot\tilde{\mathbf{p}}&\omega\tilde{\mathcal{T}}(\tilde{\mathbf{r}})\\\
\omega\tilde{\mathcal{T}}(\tilde{\mathbf{r}})^{\dagger}&{\hbar}{\mathrm{v}_{\textsc{f}}}k_{\theta}\bm{\sigma}_{-\theta/2}\cdot\tilde{\mathbf{p}}\end{bmatrix}.$
(5)

Eq. (5) depends on the dimensionless hopping operator
$\tilde{\mathcal{T}}(\tilde{\mathbf{r}})$, obtained by replacing $T_{l}$ in
Eq. (3) with the $\omega$-independent Fourier coefficients $\tilde{T}_{l}$.
Two characteristic energies are outlined. The Dirac-cone energy
${\hbar}\mathrm{v}_{\textsc{f}}k_{\theta}$, associated to the (dimensionless)
unperturbed Hamiltonian

$\tilde{\mathcal{H}}_{0}=\begin{bmatrix}\bm{\sigma}_{\theta/2}\cdot\tilde{\mathbf{p}}&0\\\
0&\bm{\sigma}_{-\theta/2}\cdot\tilde{\mathbf{p}}\end{bmatrix},$ (6)

and the interlayer coupling energy $\omega$, associated to the (dimensionless)
perturbing potential

$\tilde{\mathcal{H}}_{\textsc{t}}=\begin{bmatrix}0&\tilde{\mathcal{T}}(\tilde{\mathbf{r}})\\\
\tilde{\mathcal{T}}(\tilde{\mathbf{r}})^{\dagger}&0\end{bmatrix}.$ (7)

The BM Hamiltonian may be thus reexpressed as

$\mathcal{H}=\hbar
k_{\theta}\mathrm{v}_{\textsc{f}}\left(\tilde{\mathcal{H}}_{0}+\gamma\tilde{\mathcal{H}}_{\textsc{t}}\right),\qquad\gamma=\dfrac{\omega}{\hbar
k_{\theta}\mathrm{v}_{\textsc{f}}}.$ (8)

Notably, this operator depends on the single parameter $\gamma$ that controls
the physics of the system. This implies that varying $\omega$ and $\theta$,
while keeping $\gamma$ constant, leaves the eigenstates of $\mathcal{H}$
unchanged, while the corresponding eigenvalues are scaled by the common factor
$\hbar k_{\theta}\mathrm{v}_{\textsc{f}}$.

Numerical calculations show the appearance of quasi-flat energy bands with
vanishing Fermi velocity, for suitable choices of $\gamma$ and $\omega$. In
particular, $\gamma{=}3^{-1/2}$ and $\omega{=}\omega_{\textsc{ma}}{=}0.11$ eV
lead to a sequence of magic angles with maximum value
$\theta_{\textsc{ma}}{\sim}1.05^{\circ}$ [7].

Let us now denote by $\mathcal{O}_{\textsc{bm}}$ the set of all Hamiltonians
of the BM type given by Eq. (8), depending on $\omega$ and $\theta$. Let us
consider two generic elements
$\mathcal{H}_{1}=\mathcal{H}(\omega_{1},\theta_{1})$ and
$\mathcal{H}_{2}=\mathcal{H}(\omega_{2},\theta_{2})$ in
$\mathcal{O}_{\textsc{bm}}$. We may introduce the equivalence relation

$\mathcal{H}_{1}\underset{\gamma}{\sim}\mathcal{H}_{2}\iff\frac{\omega_{1}}{\hbar
k_{\theta_{1}}\mathrm{v}_{\textsc{f}}}=\frac{\omega_{2}}{\hbar
k_{\theta_{2}}\mathrm{v}_{\textsc{f}}}=\mathcal{\gamma}\in\mathbb{R}$ (9)

By doing so, we have partitioned $\mathcal{O}_{\textsc{bm}}$ into infinite
equivalence classes of the form

$[\mathcal{H}]_{\gamma}\equiv\\{\mathcal{H}\in\mathcal{O}_{\textsc{bm}}:\mathcal{H}_{\gamma}\underset{\gamma}{\sim}\mathcal{H}\\},$
(10)

where $\mathcal{H}_{\gamma}$ is a representative of the $\gamma$ class. If
$\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ belong to $[\mathcal{H}]_{\gamma}$,
then they have the same eigenstates and present _homothetic spectra_ , which
differ by a multiplication factor. In particular, the spectrum of
$\mathcal{H}_{1}$ is obtained from the spectrum of $\mathcal{H}_{2}$
multiplied by

$s^{1}_{2}=\frac{k_{\theta_{1}}}{k_{\theta_{2}}}=\frac{\sin(\theta_{1}/2)}{\sin(\theta_{2}/2)}.$
(11)

Among the different $\gamma$ classes in $\mathcal{O}_{\textsc{bm}}$, we can
identify a specific equivalent class, denoted $\gamma_{\textsc{ma}}$, for
which a representative element is just one of the above outlined MA
Hamiltonians, for example the one set by $\omega{=}\omega_{\textsc{ma}}$ and
$\theta{=}\theta_{\textsc{ma}}$. Accordingly, all the Hamiltonians in
$[\mathcal{H}]_{\gamma_{\textsc{ma}}}$ present an equal number of energy bands
with vanishing group velocity at the Fermi level.

At this point, we use the fact that the interlayer hopping energy $\omega$
depends on the interlayer distance $d$ [35, 33]. We may then assume that a
bijective function $F$ of $d$ exist, being such that

$\omega(d)=\omega(d_{0})F(d),$ (12)

where $F(d_{0})=1$ at a reference distance $d_{0}$.

Let us consider a specific equivalence class $\gamma$, spanned by the
representative twist angle $\theta_{0}$ at the interlayer distance $d_{0}$. By
the defining relation (9) and the bijectivity of $F$, there exists a unique
interlayer distance

$d_{\gamma}(\theta)=F^{-1}\left(\frac{k_{\theta}}{k_{\theta_{0}}}\right),$
(13)

for any other twist angle $\theta$, which determines a BM Hamiltonian of the
same class.

Let us now focus back on the equivalence class of MATBG, represented by the
Hamiltonian $\mathcal{H}_{\gamma_{\textsc{ma}}}$ with geometric parameters
$\theta_{\textsc{ma}}$ and $\omega_{\textsc{ma}}$. We may set
$\theta_{0}{=}\theta_{\textsc{ma}}$ and associate $\omega_{\textsc{ma}}$ to
the well-known interlayer distance $d_{0}{=}3.349$ Å of bilayer graphene or
graphite. Then, we can take any twist angle $\theta$ and its corresponding
_magic_ interlayer distance $d_{\gamma_{\textsc{ma}}}(\theta)$. These two
parameters generate a BM Hamiltonian in
$[\mathcal{H}]_{\gamma_{\textsc{ma}}}$, which hosts the same number of bands,
with vanishing Fermi velocity, as the MATBG Hamiltonian. However, the tiny
dispersive widths observed in the quasi-flat bands of MATBG on the meV scale,
are amplified by the $k_{\theta}/k_{\theta_{\textsc{ma}}}$ factor, which may
pose a limitation on the Hamiltonian of the $\gamma_{\textsc{ma}}$ class with
large stacking angles.

The construction carried out so far is valid at the CM level, within the
limits outlined in appendix B. Inclusion of relaxation effects is possible by
a straightforward redefinition of the Fourier coefficients $\tilde{T}_{l}$, as
mentioned in appendix C. Nonetheless, $d$ and $\theta$ are two experimental
geometric parameters, which uniquely determine, up to symmetry operations, the
configuration of a TBG. For this reason, it may be possible to construct an
equivalence class for the electronic states near the Fermi level, regardless
of the underlying theory.

Such a conjecture leads to two remarkable consequences. (i) In experiments
involving TBG samples at small twist angles, the magic quasi-flat bands at the
Fermi level should be achieved by tuning an external pressure, as reported in
some previous studies [20, 21, 22, 23, 24]. (ii) The possibility to have
quasi-flat bands at the Fermi level, in a TBG of twist angles relatively
larger than the magic angle, may allow the characterization of MATBG, using an
equivalent auxiliary system made of relatively small numbers of atoms per unit
cell. While other scaling procedures with the same aim do exist[25, 26, 27,
28] (see also AppendixA), connecting them to the defining parameters of TBG,
the stacking angle and distance, may open to atomistic approaches, beyond the
CM and TB approaches.

##  III An equivalence relation for tight-binding Hamiltonians

The present section is devoted to numerically assess the existence of the
MATBG equivalence class, and derive a suitable functional form for Eq. (13).
As a validation method, we employ a tight-binding approach, focusing on the
electronic properties of TBG structures at different values of $\theta$ and
$d$. The TB Hamiltonian takes the quadratic form

$\mathcal{H}=\sum_{i}\epsilon_{i}c^{\dagger}_{i}c_{i}-\sum_{i\neq
j}t_{ij}c^{\dagger}_{i}c_{j},$ (14)

expressed in terms of the creation operators $c^{\dagger}_{i}$, and
annihilation operators $c_{i}$, for an electron with energy $\epsilon_{i}$ at
the lattice site $i$. The $t_{ij}$ matrix elements define the hopping
amplitudes between lattice sites $i$ and $j$, depending on their relative
position vector $\mathbf{r}_{ij}$, of modulus $r_{ij}$ and perpendicular
component $z_{ij}$ to the graphene sheets.

A minimal description of TBG is expressed in terms of $p_{z}$ orbitals
centered at the lattice sites [36], within the Slater-Koster (SK) formalism
[37]. The on-site energies $\epsilon_{i}$ have the $p$-orbital symmetry. The
hopping amplitudes are split into intralayer ($pp\pi$) and interlayer
($pp\sigma$) parts,

$t_{ij}=n_{ij}^{2}V_{pp\sigma}(r_{ij})+(1-n_{ij}^{2})V_{pp\pi}(r_{ij}),$ (15)

modulated by the out-of-plane directional cosines $n_{ij}=z_{ij}/r_{ij}$. The
non-standard SK parameterization

$\displaystyle V_{pp\pi}(r_{ij})=$ $\displaystyle-
t_{0}e^{-\delta(r_{ij}-a_{c})}c(r_{ij})$ (16) $\displaystyle
V_{pp\sigma}(r_{ij})=$ $\displaystyle
t_{1}e^{-\eta(d-d_{0})}e^{-\delta(r_{ij}-d)}c(r_{ij})$ (17)

provides a reasonably accurate model, which accounts for structural
deformation and out-of-plane atomic corrugation effects [24]. Here,
$a_{c}{=}a/\sqrt{3}{=}1.42$ Å is the in-plane nearest-neighbor distance
between two atomic sites. $d_{0}{=}3.349$ Å labels the above-introduced
equilibrium distance between the TBG planes, without applied pressure.
$t_{0}{=}2.8$ eV and $t_{1}{=}0.44$ eV respectively denote the intralayer and
interlayer hopping energies. $c(r){=}(1+e^{(r-r_{c})/l_{c}})^{-1}$ is a smooth
function in which $r_{c}{=}0.265$ Å and $l_{c}{=}5.0$ Å. The ranges of
$V_{pp\pi}$ and $V_{pp\sigma}$ are controlled by $\delta{=}2.18$ Å-1. Vertical
compression effects are encompassed in $\eta{=}0.58$ Å-1.

TB calculations are performed using the TBPLas package [38], by retaining all
the hoppings with $r_{ij}{\lesssim}7.5$ Å. Using the standard SK
parameterization [39], i.e., plugging $d{=}d_{0}$ in Eq. (17), the TB
Hamiltonian presents four nearly flat bands with vanishing Fermi velocity at
$\theta{=}1.050^{\circ}$. This fact is reported in Fig. 2a and 2c, where the
band structure of MATBG is plotted along the
$\tilde{\Gamma}\tilde{K}\tilde{M}\tilde{\Gamma}$ border of the irreducible
$1^{\mathrm{st}}$BZ of the system, depicted in Fig. 1b.

To further simplify the treatment, we select TBG structures in which the
mismatch angle between the two graphene layers leads to commensurate
arrangements [29, 40, 41, 8]. In this specific scenario, using AA-stacked
bilayer graphene as reference, the angle of mismatch between the two TBG
layers reads

$\cos(\theta_{m})=\frac{3m^{2}+3m+1/2}{3m^{2}+3m+1},$ (18)

with $m$ being an integer. Moving the coordinate system so that it aligns with
the bottom graphene sheet, the moiré super-structure for commensurate TBG,
with twist angle $\theta_{m}$, is generated by the primitive vectors

$\begin{bmatrix}\mathbf{t}_{1}(m)\\\
\mathbf{t}_{2}(m)\end{bmatrix}=\begin{bmatrix}m&m+1\\\
-(m+1)&2m+1\end{bmatrix}\begin{bmatrix}\mathbf{a}_{1}\\\
\mathbf{a}_{2}\end{bmatrix}.$ (19)

The latter are obtained as linear combinations of the chosen basis vectors
$\mathbf{a}_{1}$, $\mathbf{a}_{2}$ for unrotated graphene, see discussion in
Sec. II and Fig. 1a.

As discussed above and shown in Fig. 2a, with the TB parameters at hand, the
Hamiltonian of Eq. (14) satisfies the MA condition, at the geometric
parameters $\theta_{31}{=}1.050^{\circ}$, obtained by letting $m{=}31$ in Eq.
(18), and $d_{31}{=}d_{0}$. The next step is constructing an equivalence class
for this Hamiltonian. To do so, we have sorted the $d$ parameter to different
interlayer distances, at fixed twist angle $\theta_{m}$. In the procedure, we
have looked for a band structure hosting nearly non-dispersive band states at
the Fermi level, which most closely resemble the MA bands when compressed by
the scaling factor $s^{31}_{m}$, cfr. Eq. (11). We found a well-converged
_magic distance_ that fulfills this conditions, for each tested $\theta_{m}$.

To better illustrate the process, we consider two random twist angles greater
than $\theta_{31}$, namely, $\theta_{9}{=}3.481^{\circ}$ and
$\theta_{11}{=}2.876^{\circ}$. The sorting procedure returns the magic
distances $d_{9}{=}2.7329$ Å, for $\theta_{9}$, and $d_{11}{=}2.8304$ Å, for
$\theta_{11}$. Let us inspect the energy spectra of the two TBG systems.

Figure 2:  (a) Band structure of MATBG, as obtained with the geometric
parameters $\theta_{31}{=}1.050^{\circ}$ and $d_{31}=3.349$ Å, within the TB
approach illustrated in the main text. (b) Band dispersions of two TBG
structures with $\theta_{9}{=}3.481^{\circ}$ and $d_{9}{=}2.7329$ Å (green
lines), $\theta_{11}{=}2.876^{\circ}$ and $d_{11}{=}2.8304$ Å (blue lines).
(c) Compressed band dispersions of panel (b), obtained with the scaling
factors $s^{31}_{9}=0.3017$ (green), $s^{31}_{11}=0.3652$ (blue), and
superimposed to the MATBG bands (red). The three spectra, reported against the
$\tilde{\Gamma}\tilde{K}\tilde{M}\tilde{\Gamma}$ contour of the moiré
$1^{\mathrm{st}}$BZ, shown in Fig. 1b, coincide within ${\sim}0.04$ eV around
the Fermi level (gray shaded area). (d)-(f) Projected local density at the
Fermi level of (d) MATBG, (e) TBG with $\theta_{9}$-$d_{9}$, and (f) TBG with
$\theta_{11}$-$d_{11}$, calculated by applying the TB propagation method. The
in-plane coordinates are provided both in nm and moiré lattice units (mlu).
The black parallelograms denote the moiré unit cells of the three TBG systems.
Density peaks occur around the AA staking regions, while density dips
characterize the AB stacking regions.

Fig. 2b reports the band structures obtained with $\theta_{9}$-$d_{9}$ (green
lines) and $\theta_{11}$-$d_{11}$ (blue lines), which, at a rough look, appear
to differ by an overall scaling factor. Furthermore, the MATBG spectrum in
Fig. 2a is extraordinarily similar to the TBG spectra in Fig. 2b, over a
different energy scale. We may then apply the scaling factors
$s^{31}_{9}{=}k_{\theta_{31}}/k_{\theta_{9}}{=}0.3017$ and
$s^{31}_{11}{=}k_{\theta_{31}}/k_{\theta_{11}}{=}0.3652$ to the bands in Fig.
2b, and compare them with the MATBG band system.

As shown in Fig. 2c, the scaled dispersions from the two TBG samples
practically coincide with the MATBG dispersions in their low-energy parts,
within ${\sim}0.04$ eV around the Fermi level, which establishes the range of
validity of the CM, see also appendix B. Strictly speaking, the equivalence
class $\gamma_{\textsc{ma}}$ is well-defined within this range. Nonetheless,
we may extend the definition (9) at the TB level, by saying that two TB
Hamiltonian are equivalent if their energy bands differ by an overall scaling
factor in the CM limit.

We have thus reasonably proved that the low-energy spectra from TBG structures
with $\theta_{9}$-$d_{9}$ and $\theta_{11}$-$d_{11}$ are both mapped onto the
magic angle spectrum. As a further step towards the validation of the
equivalence class, we look at the space distribution of the quasi-flat band
states. In particular, we focus on the local electronic density, projected
onto the quasi-flat-band range. Fig. 2d-f reports the projected density
profiles for the three TBG systems, which appear to be strongly localized
around the AA stacking regions of the corresponding moiré lattices, while they
are practically vanishing in the AB stacking regions. This result is consisted
with a number of previous studies [33, 42, 14, 43, 44]. More importantly, it
allows us to conclude that the eigenstates corresponding to quasi-flat bands,
in the same equivalence class, share the same localization properties.

Figure 3: Magic distance $d_{\textsc{m}}$ vs stacking angle $\theta$, as
obtained through the TB method of the main text. Also reported on the top
horizontal axis are the values of $m$, cfr. Eq. (18). The red dashed line
corresponds to the merit function of Eq. (21), with the best-fit parameter
$\beta=1.957$ Å-1.

The just outlined procedure was applied to a variety of TBG systems, with
twist angles falling in the range of validity of the CM, see appendix B. In
particular, Fig. 3 reports the interlayer magic distances $d_{\textsc{m}}$ for
which nearly non-dispersive bands were obtained as function of the twist angle
$\theta$, in the range of $\theta_{31}{=}1.050^{\circ}$ to
$\theta_{5}{=}6.009^{\circ}$. As already pointed out, the main effect of
reducing the distance between the TBG layers is to modify the interlayer
interaction. In the TB approach, this results in a variation of the interlayer
hopping integral. Turning to the CM limit, we want to derive a smooth and
bijective dependence of the hopping parameter $\omega$ on the interlayer
distance $d$, thus providing an explicit expression to the $F$ function in Eq.
(12). We may assume the simple exponential scaling law

$\omega(d)=\omega(d_{0})e^{-\beta(d-d_{0})},\quad d\geq d_{0},$ (20)

parametrized in terms of the inverse damping length $\beta$. Then, placing
ourselves in the equivalence class of MATBG, we obtain that the data of Fig. 3
are correlated by the following merit function

$d_{\textsc{m}}(\theta)=d_{0}-\frac{1}{\beta}\ln\left[\frac{\sin({\theta}/{2})}{\sin({\theta_{0}}/{2})}\right],$
(21)

with $d_{0}=d_{31}$ and $\theta_{0}=\theta_{31}$. Eq. (21) is indeed in
excellent agreement with the data at the best-fit value $\beta=1.957$ Å-1, as
also reported in Fig. 3.

We have thus provided reasonable arguments that an equivalence class for the
MA condition exists and is well-defined at the TB level.

##  IV Equivalence-class approach and DFT

We now apply the equivalence class concept to DFT, suggesting how the low-
energy features of some selected TBG systems, arranged in the commensurate
lattices, set by Eqs. (18) and (19), can be linked to the MA condition. DFT
calculations were run using the Quantum Espresso (QE) package [45, 46], within
the projector augmented wave (PAW) method [47] provided by the PseudoDojo PAW
datasets [48]. The Perdew-Burke-Ernzerhof (PBE) exchange-correlation
functional [49] was adopted, in conjunction with a plane-wave cutoff of $30$
Ry, and an energy convergence criterion of $10^{-6}$ Ry for the self-
consistent runs. The TBG slab was given an out-of-plane vacuum region of $20$
Å. The $1^{\mathrm{st}}$BZ was sampled over Monkhorst-Pack grids [50] of
$2{\times}2$ to $6{\times}6$ points, equivalent to average spacings between
$0.0004$ Å-1 and $0.007$ Å-1. A Marzari-Vanderbilt smearing function [51] of
$0.01$ Ry in width was included, to fasten convergence.

We focused on three specific stacking angles, namely,
$\theta_{5}{=}6.009^{\circ}$, $\theta_{6}{=}5.086^{\circ}$ and
$\theta_{7}{=}4.408^{\circ}$, being sufficiently larger than the magic angle,
while falling around the range of applicability of the CM model. Then, we
adopted a similar sorting procedure as the TB method, that is, we carried out
several DFT calculations by varying the interlayer distance until nearly non-
dispersive KS bands at the Fermi level were detected. The DFT search returned
the magic distances $\tilde{d}_{5}{=}2.450$ Å, $\tilde{d}_{6}{=}2.540$ Å and
$\tilde{d}_{7}{=}2.616$ Å, which gave the least dispersive band dispersions at
the Fermi level, within numerical uncertainties below $0.05$ Å. Interestingly,
the TB predictions for the magic distances are $d_{5}{=}2.470$ Å,
$d_{6}{=}2.546$ Å and $d_{7}{=}2.615$ Å, as extracted from Fig. 3 with the
same stacking angles. This very good agreement suggests using the TB magic
distances as a starting guess to speed up the DFT search process.

In Fig. 4, we report the DFT spectra obtained with
$\theta_{5}$-$\tilde{d}_{5}$, $\theta_{6}$-$\tilde{d}_{6}$ and
$\theta_{7}$-$\tilde{d}_{7}$ by respectively processing $364$, $508$ and $676$
atoms per unit-cell. Looking at Fig. 4a-c, we see the main features of the
band structure, given by the presence of four nearly flat bands around the
Fermi level, with vanishing group velocity at the corners of the moiré
$1^{\mathrm{st}}$BZ, namely the $\tilde{K}$ points and their inequivalent
counterparts shown in Fig. 1b. For the angle $\theta_{5}$, even though we
manage to extract a suitable magic distance, the four flat bands, compressed
by $s^{31}_{5}$, exhibit a greater bandwidth as compared to the
$\theta_{6}$-$\tilde{d}_{6}$ and $\theta_{7}$-$\tilde{d}_{7}$ cases. A wider
bandwidth was also detected at the TB level, which demonstrates that
$\theta{\sim}6^{\circ}$ is the upper limit of applicability of our results. In
Fig. 4d-f, we show the local electron density, projected onto the quasi-non-
dispersive states, in the range of $-0.05$ eV to $0.05$ eV, around the Fermi
energy. This quantity shows that the wave functions of the four bands are
strongly localized in the AA stacking regions.

Figure 4: (a)-(c) DFT band structure of a TBG with the geometric parameters
(a) $\theta_{7}{=}4.408^{\circ}$-$\tilde{d}_{7}{=}2.616$ Å, (b)
$\theta_{6}{=}5.086^{\circ}$-$\tilde{d}_{6}{=}2.540$ Å and (c)
$\theta_{5}{=}6.009^{\circ}$-$\tilde{d}_{5}{=}2.450$ Å. The energy ticks in
the secondary vertical axis (red) are scaled to the width of the MATBG energy
spectrum of Fig. 2a, by (a) $s^{31}_{7}{=}0.2383$, (b) $s^{31}_{6}{=}0.2065$
and (c) $s^{31}_{5}{=}0.1748$. (d)-(f) Localized density associated to the
four quasi-flat bands for (d) $\theta_{7}$-$\tilde{d}_{7}$, (e)
$\theta_{6}$-$\tilde{d}_{6}$ and (f) $\theta_{5}$-$\tilde{d}_{5}$, near charge
neutrality. All other settings are as in Fig. 2.

We should also point out that, due to the large unit cell size of the magic-
angle twisted bilayer graphene, we do not have direct access to its
fundamental ground-state properties, through DFT. Indeed, the twist angle
$\theta_{31}$ produces a commensurate TBG lattice with $11908$ atoms per unit
cell. An atomistic self-consistent computation for such a system would be
unfeasible even by largest, currently available, high-performance computing
resources. Furthermore, it is not guaranteed that the MA condition is
fulfilled for the $\theta_{31}{\sim}1.050^{\circ}$ value, at the DFT level.
Therefore, the $\theta_{7}$-$\tilde{d}_{7}$, $\theta_{6}$-$\tilde{d}_{6}$ and
$\theta_{5}$-$\tilde{d}_{5}$ conditions may not precisely fall in the same
equivalence class of MATBG. However, they both share with MATBG the two main
features, that is four nearly flat bands at the Fermi level with vanishing
group velocity at the $1^{\mathrm{st}}$BZ corners, and eigenfunctions strongly
localized in the AA stacking region of the moiré lattice. In addition, the
agreement of TB and DFT approaches on the magic distances, suggests that the
inaccessible MA spectrum of DFT at $\theta_{31}$ may be close to the ones
obtained by applying the TB compression factors $s^{31}_{7}{=}0.2383$,
$s^{31}_{6}{=}0.2065$ and $s^{31}_{5}{=}0.1748$ to our DFT spectra. The
compressed DFT band structures of TBG, for the $\theta_{7}$-$\tilde{d}_{7}$,
$\theta_{6}$-$\tilde{d}_{6}$ and $\theta_{5}$-$\tilde{d}_{5}$ cases, can be
respectively read from Fig. 4a-c, by looking at the secondary vertical axis
with red ticks. The three scaled spectra show very similar dispersions within
$100$ meV around the Fermi level, however, they should present identical
features. We here interpret these differences as a measure of error in finding
the DFT features of MATBG by the equivalent class approach. We may accordingly
infer that the three considered TBG systems are in a _nearby_ equivalence
class to that of MATBG, being characterized by the interaction parameter
$\gamma{=}\gamma_{\textsc{ma}}+\delta\gamma$, in the CM, with
${\delta\gamma/\gamma_{\textsc{ma}}}{\ll}1$.

To further support this statement, in Fig. 5, we compare the results for
MATBG, as obtained with the TB approach outlined in Sec. III, and the above
discussed $\theta_{7}$-$\tilde{d}_{7}$ and $\theta_{6}$-$\tilde{d}_{6}$ cases,
along with compression factors $s^{31}_{7}$ and $s^{31}_{6}$. We see that the
two DFT compressed spectra are practically overlapping, see Fig. 5a-c, within
an error window of less than ${\sim}5$ meV, which occurs in the quasi-flat
band region, see Fig. 5c.

Figure 5: (a)-(c) DFT band structure of TBG, with geometric parameters
$\theta_{6}$-$\tilde{d}_{6}$, $\theta_{7}$-$\tilde{d}_{7}$, and compression
factors $s^{31}_{6}$, $s^{31}_{7}$, vs. TB band structure of MATBG. (a)
Occupied dispersive levels, (b) unoccupied dispersive levels, (c) quasi-flat
bands. (d), (e) Projected electron density onto the quasi-flat band states for
(d) MATBG and (e) TBG structures of (a)-(c).

In addition, the DFT and TB occupied dispersive levels are also very close to
one another, see Fig. 5a. Similar conclusions can be raised on the unoccupied
dispersive parts of the spectra, given in Fig. 5b, apart from a global shift
of about $20$ meV. The DFT and TB quasi-flat bands, on the other hand, appear
to have rather different behaviors, see Fig. 5c. Nonetheless, the tiny
dispersive ranges of these bands present similar widths, namely, ${\sim}16$
meV with $\theta_{6}$-$\tilde{d}_{6}$ and ${\sim}10$ meV with either MATBG or
$\theta_{7}$-$\tilde{d}_{7}$. In addition, the projected electron densities
onto these bands display consistent peaks and dips, when reported on the same
moiré-lattice-unit scale, see Fig. 5c,e. We also need to point out that TB is
a parametric approach, adjusted to experiments or consolidated theory, whereas
DFT is an ab initio approach. Thus, the recorded differences between the two
methods cannot be surprising. Indeed, TB and DFT calculations present even
more marked differences when applied to known systems, such as graphene,
leading to $15$-$20$ % inaccuracies in the estimate of the Fermi velocity [52]
and unoccupied band levels [53].

##  V Discussion

The significance of elucidating the existence of an equivalence class lies in
the fact that MATBG, due to its very large moiré unit cell, represents a
computationally inaccessible system, even to its ground state properties, when
ab initio type methods are employed. The existence of equivalence classes, and
the MATBG equivalence class in particular, may open the way to computations of
ground-state and excited-state features, using more convenient geometric
structures that are computationally affordable.

As a paradigmatic example, we suggest concentrating on the the dynamical
density-density response function of the TBG electrons around the quasi-flat
band region [54]. The latter can be put in a form that depends on the DFT
energies and wavefunctions of the system, plus, the electronic temperature
embedded in the Fermi-Dirac distribution. At least in the random phase
approximation (RPA), the interaction can be reduced to a smooth Coulomb-like
form. By the equivalence-class approach developed here, it may be possible to
map one density-density response function onto another for structures
belonging to the same equivalence class of MATBG, by carrying out TDDFT
computations in moiré lattice units and suitably adjusting temperature, poles,
and amplitudes using overall scaling factors similar to $s^{1}_{2}$ of Eq.
(11). Following preliminary calculations to be run at the CM and/or TB levels,
we may thus explore the nature of plasmon excitations in MATBG, setting up a
computationally accessible strategy on an auxiliary TBG system, with a more
convenient geometry. Indeed, a systematic study of light-trapping mechanisms
related to plasmonic effects in this area is still lacking. Recent
experimental studies seem to rule out a purely electronic coupling [55, 56],
whereas, from a theoretical point of view, it has been proposed that
collective excitations play a crucial role in the emergence of unconventional
superconductivity phenomena [57, 58, 59].

Building upon the same considerations, it might be further possible to
investigate the superconducting-normal state transition of a TBG with quasi-
flat bands at a specific magic distance. In this case, the critical
temperature should be re-normalized with a suitable scaling law, because we
expect that critical temperature increases with increasing the stacking angle.

We finally observe that results discussed here are not limited to TBG. An
equivalence relation could be defined on other twisted materials, made of two
or multiple atomic layers that host Dirac-cone states in their isolated
phases. Indeed, the weak interlayer interaction should always lead, under
suitable limiting conditions, to a reference Hamiltonian of the BM type, which
should be expressed in a form similar to Eq. (8). A crucial point of the
treatment, particularly challenging in heterogeneous multilayer structures,
would be to construct a smooth, bijective function of the hopping parameters
to the interlayer distances.

##  VI Conclusion

In conclusion, we have made explicit that MATBG belongs to a specific
equivalence class, being well-defined in the framework of a continuum theory.
All TBG structures of this class are generated by the twist angle and a
_magic_ interlayer distance. They exhibit homothetic spectra, fixed number of
bands with vanishing Fermi velocity, and a common basis of eigenstates, which
can be one-to-one mapped to MATBG through a dimensionless scaling factor
associated to the twist angle. The equivalence relation can be extended to
tight-binding models, by running a search algorithm on the interlayer distance
at fixed twist angle, until the best homothetic match in the quasi-flat bands
of the selected TBG with those of MATBG is found. Restricting our observations
to commensurate lattices, we have identified a one-to-one smooth
correspondence between TBG structures with twist angles below
${\sim}6^{\circ}$ and MATBG. Next, we have analyzed the DFT ground-state of
superimposed graphene monolayers, twisted over the same limited angular range.
We have set up an ab initio procedure to determine the optimal interlayer
distance, for each tested case, which leads to quasi-flat bands with vanishing
group velocity at the Fermi level. Finally, we have forced the validity of the
equivalent-class concept to tentatively reconstruct the inaccessible DFT
spectrum of MATBG. Our results may open the possibility to study, by using
costly ab initio techniques, the electronic properties of systems with
affordable numbers of atoms per unit cell, but with magic-angle-like
properties. Additionally, they may guide the _experimental mapping_ of a TBG
system onto MATBG by suggesting the correct magic distance to be achieved in
compressing the system with an external pressure.

##  VII Acknowledgements

We thank N. Lo Gullo, J. Settino, and F. Plastina for useful discussion. This
research was partially supported _Centro Nazionale di Ricerca in High-
Performance Computing, Big Data and Quantum Computing_ , PNRR 4 2 1.4, CI
CN00000013, CUP H23C22000360005. We acknowledge the _Marconi_ , _Marconi_ 100,
_Galileo_ 100 and Leonardo high performance computing resources, provided by
the CINECA consortium (Italy), within the INF16_npqcd project, under the
CINECA-INFN agreement. We also acknowledge the _Newton_ high performance
computing cluster, provided by the University of Calabria.

## 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, _et al._ , Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018b).
  * 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).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, _et al._ , Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Flavour hund’s coupling, chern gaps and charge diffusivity in moiré graphene, Nature 592, 43–48 (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).
  * 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, Phys. Rev. B 86, 155449 (2012).
  * 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, Phys. Rev. B 82, 121407 (2010).
  * Rosendo López  _et al._ [2020] M. Rosendo López, F. Peñaranda, J. Christensen, and P. San-Jose, Flat bands in magic-angle vibrating plates, Phys. Rev. Lett. 125, 214301 (2020).
  * Xie  _et al._ [2019] Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene, Nature 572, 101 (2019).
  * 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, _et al._ , Maximized electron interactions at the magic angle in twisted bilayer graphene, Nature 572, 95 (2019).
  * 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).
  * Li  _et al._ [2010] G. Li, A. Luican, J. Lopes dos Santos, A. Castro Neto, A. Reina, J. Kong, and E. Andrei, Observation of van hove singularities in twisted graphene layers, Nature physics 6, 109 (2010).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, _et al._ , Electronic correlations in twisted bilayer graphene near the magic angle, Nature physics 15, 1174 (2019).
  * 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, _et al._ , Observation of flat bands in twisted bilayer graphene, Nature Physics 17, 189 (2021).
  * Carr  _et al._ [2017] S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, and E. Kaxiras, Twistronics: Manipulating the electronic properties of two-dimensional layered structures through their twist angle, Phys. Rev. B 95, 075420 (2017).
  * Elizabeth [2019] G. Elizabeth, How magic angle graphene is stirring up physics, Nature 565, 15 (2019).
  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Van der waals heterostructures, Nature 499, 419 (2013).
  * Chittari  _et al._ [2018] B. L. Chittari, N. Leconte, S. Javvaji, and J. Jung, Pressure induced compression of flatbands in twisted bilayer graphene, Electronic Structure 1, 015001 (2018).
  * Yndurain [2019] F. Yndurain, Pressure-induced magnetism in rotated graphene bilayers, Phys. Rev. B 99, 045423 (2019).
  * Carr  _et al._ [2018] S. Carr, S. Fang, P. Jarillo-Herrero, and E. Kaxiras, Pressure dependence of the magic twist angle in graphene superlattices, Phys. Rev. B 98, 085144 (2018).
  * Ge  _et al._ [2021a] L. Ge, K. Ni, X. Wu, Z. Fu, Y. Lu, and Y. Zhu, Emerging flat bands in large-angle twisted bi-layer graphene under pressure, Nanoscale 13, 9264 (2021a).
  * Lin  _et al._ [2020] X. Lin, H. Zhu, and J. Ni, Pressure-induced gap modulation and topological transitions in twisted bilayer and twisted double bilayer graphene, Phys. Rev. B 101, 155405 (2020).
  * Sainz-Cruz  _et al._ [2021] H. Sainz-Cruz, T. Cea, P. A. Pantaleón, and F. Guinea, High transmission in twisted bilayer graphene with angle disorder, Phys. Rev. B 104, 075144 (2021).
  * Vahedi  _et al._ [2021] J. Vahedi, R. Peters, A. Missaoui, A. Honecker, and G. T. de Laissardière, Magnetism of magic-angle twisted bilayer graphene, SciPost Phys. 11, 083 (2021).
  * Gonzalez-Arraga  _et al._ [2017] L. A. Gonzalez-Arraga, J. L. Lado, F. Guinea, and P. San-Jose, Electrically controllable magnetism in twisted bilayer graphene, Phys. Rev. Lett. 119, 107201 (2017).
  * Andrade  _et al._ [2023] E. Andrade, P. A. Pantaleón, F. Guinea, and G. G. Naumis, Flat bands and electronic localization in twisted bilayer graphene nanoribbons, Phys. Rev. B 108, 235418 (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, Phys. Rev. Lett. 99, 256802 (2007).
  * Balents [2019] L. Balents, General continuum model for twisted bilayer graphene and arbitrary smooth deformations, SciPost Phys. 7, 048 (2019).
  * Carr  _et al._ [2019] S. Carr, S. Fang, Z. Zhu, and E. Kaxiras, Exact continuum model for low-energy electronic states of twisted bilayer graphene, Phys. Rev. Res. 1, 013001 (2019).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Origin of magic angles in twisted bilayer graphene, Phys. Rev. Lett. 122, 106405 (2019).
  * Trambly de Laissardière  _et al._ [2012] G. Trambly de Laissardière, D. Mayou, and L. Magaud, Numerical studies of confined states in rotated bilayers of graphene, Phys. Rev. B 86, 125413 (2012).
  * Bernevig  _et al._ [2021] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Twisted bilayer graphene. iii. interacting hamiltonian and exact symmetries, Phys. Rev. B 103, 205413 (2021).
  * Bistritzer and MacDonald [2010] R. Bistritzer and A. H. MacDonald, Transport between twisted graphene layers, Phys. Rev. B 81, 245412 (2010).
  * Lin and Tománek [2018] X. Lin and D. Tománek, Minimum model for the electronic structure of twisted bilayer graphene and related structures, Phys. Rev. B 98, 081410 (2018).
  * Slater and Koster [1954] J. C. Slater and G. F. Koster, Simplified lcao method for the periodic potential problem, Phys. Rev. 94, 1498 (1954).
  * Li  _et al._ [2023] Y. Li, Z. Zhan, X. Kuang, Y. Li, and S. Yuan, Tbplas: A tight-binding package for large-scale simulation, Computer Physics Communications 285, 108632 (2023).
  * Shi  _et al._ [2020] H. Shi, Z. Zhan, Z. Qi, K. Huang, E. v. Veen, J. A. Silva-Guillén, R. Zhang, P. Li, K. Xie, H. Ji, M. I. Katsnelson, S. Yuan, S. Qin, and Z. Zhang, Large-area, periodic, and tunable intrinsic pseudo-magnetic fields in low-angle twisted bilayer graphene, Nature Communications 11, 371 (2020).
  * Mele [2010] E. J. Mele, Commensuration and interlayer coherence in twisted bilayer graphene, Phys. Rev. B 81, 161405 (2010).
  * Shallcross  _et al._ [2010] S. Shallcross, S. Sharma, E. Kandelaki, and O. A. Pankratov, Electronic structure of turbostratic graphene, Phys. Rev. B 81, 165105 (2010).
  * Kim  _et al._ [2017] K. Kim, A. DaSilva, S. Huang, B. Fallahazad, S. Larentis, T. Taniguchi, K. Watanabe, B. J. LeRoy, A. H. MacDonald, and E. Tutuc, Tunable moiré bands and strong correlations in small-twist-angle bilayer graphene, Proceedings of the National Academy of Sciences 114, 3364 (2017).
  * Brihuega  _et al._ [2012] I. Brihuega, P. Mallet, H. González-Herrero, G. Trambly de Laissardière, M. M. Ugeda, L. Magaud, J. M. Gómez-Rodríguez, F. Ynduráin, and J.-Y. Veuillen, Unraveling the intrinsic and robust nature of van hove singularities in twisted bilayer graphene by scanning tunneling microscopy and theoretical analysis, Phys. Rev. Lett. 109, 196802 (2012).
  * Tong  _et al._ [2022] L.-H. Tong, Q. Tong, L.-Z. Yang, Y.-Y. Zhou, Q. Wu, Y. Tian, L. Zhang, L. Zhang, Z. Qin, and L.-J. Yin, Spectroscopic visualization of flat bands in magic-angle twisted monolayer-bilayer graphene: Coexistence of localization and delocalization, Phys. Rev. Lett. 128, 126401 (2022).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Quantum espresso: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with quantum espresso, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Jollet  _et al._ [2014] F. Jollet, M. Torrent, and N. Holzwarth, Generation of Projector Augmented-Wave atomic data: A 71 element validated table in the XML format, Comput. Phys. Commun. 185, 1246 (2014).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B 13, 5188 (1976).
  * Marzari  _et al._ [1999] N. Marzari, D. Vanderbilt, A. De Vita, and M. C. Payne, Thermal contraction and disordering of the al(110) surface, Phys. Rev. Lett. 82, 3296 (1999).
  * Trevisanutto  _et al._ [2008] P. E. Trevisanutto, C. Giorgetti, L. Reining, M. Ladisa, and V. Olevano, Ab initio $gw$ many-body effects in graphene, Phys. Rev. Lett. 101, 226405 (2008).
  * Sindona  _et al._ [2017] A. Sindona, M. Pisarra, C. Vacacela Gomez, P. Riccardi, G. Falcone, and S. Bellucci, Calibration of the fine-structure constant of graphene by time-dependent density-functional theory, Phys. Rev. B 96, 201408 (2017).
  * Giuliani and Vignale [2005] G. F. Giuliani and G. Vignale,  _Quantum theory of the electron liquid_ (Cambridge Univ. Press, Cambridge, 2005).
  * Stepanov  _et al._ [2020] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Untying the insulating and superconducting orders in magic-angle graphene, Nature 583, 375 (2020).
  * Ge  _et al._ [2021b] J. Ge, Y. Saito, T. Taniguchi, K. Watanabe, and A. Young, Independent superconductors and correlated insulators in twisted bilayer graphene, in  _APS March Meeting Abstracts_ (2021) pp. R42–005.
  * Sharma  _et al._ [2020] G. Sharma, M. Trushin, O. P. Sushkov, G. Vignale, and S. Adam, Superconductivity from collective excitations in magic-angle twisted bilayer graphene, Phys. Rev. Res. 2, 022040 (2020).
  * Peng  _et al._ [2023] L. Peng, I. Yudhistira, G. Vignale, and S. Adam, Theoretical determination of the effect of a screening gate on plasmon-induced superconductivity in twisted bilayer graphene (2023), arXiv:2309.14767 [cond-mat.supr-con] .
  * Cea and Guinea [2021] T. Cea and F. Guinea, Coulomb interaction, phonons, and superconductivity in twisted bilayer graphene, Proceedings of the National Academy of Sciences 118, e2107874118 (2021).
  * Uchida  _et al._ [2014] K. Uchida, S. Furuya, J.-I. Iwata, and A. Oshiyama, Atomic corrugation and electron localization due to moiré patterns in twisted bilayer graphenes, Phys. Rev. B 90, 155451 (2014).
  * van Wijk  _et al._ [2015] M. M. van Wijk, A. Schuring, M. I. Katsnelson, and A. Fasolino, Relaxation of moiré patterns for slightly misaligned identical lattices: graphene on graphite, 2D Materials 2, 034010 (2015).
  * Dai  _et al._ [2016] S. Dai, Y. Xiang, and D. J. Srolovitz, Twisted bilayer graphene: Moiré with a twist, Nano Letters 16, 5923 (2016).
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Lattice relaxation and energy band modulation in twisted bilayer graphene, Phys. Rev. B 96, 075311 (2017).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene, Phys. Rev. X 8, 031087 (2018).

##  Appendix A Alternative equivalence relations for tight-binding
Hamiltonians

In the main text, we put much effort in defining an equivalence relation for
MATBG by varying the interlayer distance. In particular, the TB calculations
presented in Sec. III were based on a Slater-Koster parametrization of the
interlayer hopping potential $V_{pp\sigma}$, cfr. Eq. (17), which allowed us
to operate on the $d$ parameter of the compression factor
$e^{-\eta(d-d_{0})}$. By doing so, we established a search algorithm to
achieve homothetic, nearly non-dispersive states, being equivalent to those of
the MATBG at angles greater than $\theta_{\textsc{ma}}$.

Given the equivalence relation of Eq. (9), we may explore alternative
approaches to derive Hamiltonians of the same equivalence class at the TB
level, by modifying the interlayer hopping in Eq. (15). For example, we can
adopt the following parametrization

$V_{pp\sigma}(r_{ij})=t_{1}\frac{\sin(\theta/2)}{\sin(\theta_{\textsc{ma}}/2)}e^{\delta(d_{0}-r_{ij})}c(r_{ij}),$
(22)

holding close to the magic angle, which corresponds to replace
$e^{-\eta(d-d_{0})}$ in Eq. (17) with the scaling factor that sets the
equivalence relation of MATBG, cfr. Eq. (11). Accordingly, we can operate on
$\theta_{\textsc{ma}}$ at constant $d_{0}$ to find an optimal magic value,
associated to $\theta$, which establishes the equivalence relation for the
corresponding TBG Hamiltonian. In other words, we can _move_ within the
equivalence class of MATBG at constant interlayer distance and look for a
renormalized magic angle $\tilde{\theta}_{M}$ to build up the equivalence
relation.

Further scaling procedures are available by modifying the other _free_
parameters in the equivalence relation (9). If one leaves the interlayer
coupling energy $\omega$ untouched, a renormalization of the angle $\theta$ is
obtained by modifying the Fermi velocity $\mathrm{v}_{\textsc{f}}$. This can
be achieved at the TB level by modifying the intralayer hopping in the
following way:

$V_{pp\sigma}(r_{ij})=-t_{0}\frac{\sin(\theta_{\textsc{ma}}/2)}{\sin(\theta/2)}e^{\delta(a_{c}-r_{ij})}c(r_{ij}).$
(23)

In this case, two _equivalent_ Hamiltonian present identical spectra in the CM
limit, because $\hbar\mathrm{v}_{\textsc{f}}{k_{\tilde{\theta}}}$ in Eq. (8)
is kept constant. A similar procedure, with additional spacial scaling to
preserve the Dirac-cone value, has been employed to reduce the number of sites
per unit cell in the TB modelling of TBG [25, 26, 27] and TBG-based systems
[28].

##  Appendix B Limits of the equivalence-class approach

The equivalence class of MATBG was robustly defined within the framework of
the single-valley BM model, as reported in Sec. II. Nonetheless, the BM
Hamiltonian of Eq. (8) reasonably accounts for the ground-state electronic
properties of misaligned graphene planes at stacking angles
${\lesssim}10^{\circ}$, regardless of whether or not the bilayer structure is
periodic. Indeed, a twist of ${\sim}10^{\circ}$ entails a reciprocal space
displacement of ${\sim}0.3$Å-1 from the Dirac point, where the Dirac-cone
bands cease to be isotropic [7]. Furthermore, the Dirac-cone approximation in
graphene holds for electronic energies below ${\sim}0.8$ eV, relative to the
Fermi level. The Tight-binding simulations discussed in Sec. III suggest that
more strict conditions should be imposed to guarantee an equivalence relation
for MATBG. Looking at Fig. 2c, we see that two TBG structures belonging to the
same equivalence class present homothetic spectral properties over a limited
energy range of about ${\sim}40$ meV around the Fermi level. In addition, we
have verified that a one-to-one mapping of TBG onto MATBG is reasonably
achieved for twist angles ${\lesssim}6^{\circ}$.

On more general grounds, a TBG system of twist angle $\theta_{m}$, cfr. Eq.
(18), is well-represented in terms of Bloch wave functions
$\psi^{m}_{n\mathbf{k}}(\mathbf{r})=e^{i\mathbf{k}\cdot\mathbf{r}}u^{m}_{n\mathbf{k}}(\mathbf{r})$
and energy levels $\epsilon^{m}_{n\mathbf{k}}$, measured from the neutrality
point, where $n$ labels the band-index, $\hbar\mathbf{k}$ denotes the crystal
momentum, and the functions $u^{m}_{n\mathbf{k}}$ have the periodicity of the
commensurate TBG lattice. Roughly speaking, the weak moiré potential couples
massless Dirac fermions between the single graphene layers with a mean
interaction energy being much less than the unperturbed Dirac-cone energy
$\hbar\mathrm{v}_{\textsc{f}}{k_{\theta_{m}}}$. Now, given an equivalence
class $\gamma$ and its representative $\mathcal{H}_{\gamma}$, the properties
of the equivalence class must be valid for
$|\epsilon^{m}_{n\mathbf{k}}|<\frac{k_{\theta_{m}}}{k_{\theta_{\gamma}}}\hbar\mathrm{v}_{\textsc{f}}k_{\theta_{\gamma}}$.
Accordingly, in either TB or DFT atomistic approaches, the equivalence class
of MATBG is expressed by the transformation relations

$\epsilon^{\textsc{ma}}_{n\mathbf{k}}=\frac{k_{\theta_{\textsc{ma}}}}{k_{\theta_{m}}}\epsilon^{m}_{n\mathbf{k}}\qquad
u^{\textsc{ma}}_{n\mathbf{k}}(k_{\theta_{\textsc{ma}}}\mathbf{r})=u^{m}_{n\mathbf{k}}(k_{\theta_{m}}\mathbf{r}),$
(24)

which hold under the limiting condition

$|\epsilon^{m}_{n\mathbf{k}}|<\frac{k_{\theta_{m}}}{k_{\theta_{\textsc{ma}}}}\hbar\mathrm{v}_{\textsc{f}}k_{\theta_{\textsc{ma}}}.$
(25)

##  Appendix C The equivalence class of MATBG with lattice relaxation

In a real TBG system, the structure undergoes a spontaneous relaxation [60,
61, 62, 63], which reflects the different interlayer interactions across the
moiré unit cell. This effect was neglected throughout the main text, because
our main goal was to identify geometric configurations exhibiting quasi-flat
bands and gather them into equivalence classes.

Nonetheless, based on the arguments presented in Sec. II, and following some
recent studies [64, 31, 63], we can express a CM formulation of the relaxed
TBG lattice in terms of the dimensionless parameter
$\alpha{=}{\omega_{\textsc{AA}}}/{\omega_{\textsc{AB}}}{\neq}1$. Accordingly,
Eq. (5) can be rewritten as

$\mathcal{H}=\begin{bmatrix}{\hbar}{\mathrm{v}_{\textsc{f}}}k_{\theta}\bm{\sigma}_{\theta/2}\cdot\tilde{\mathbf{p}}&{\omega_{\textsc{AB}}}\tilde{\mathcal{T}}(\tilde{\mathbf{r}})\\\
{\omega_{\textsc{AB}}}\tilde{\mathcal{T}}(\tilde{\mathbf{r}})^{\dagger}&{\hbar}{\mathrm{v}_{\textsc{f}}}k_{\theta}\bm{\sigma}_{-\theta/2}\cdot\tilde{\mathbf{p}}\end{bmatrix},$
(26)

and the Fourier coefficients of the intelayer potential acquire an additional
dependence on $\alpha$, namely,

$\displaystyle\tilde{T}_{1}=$ $\displaystyle\alpha\sigma_{0}+\sigma_{x},$
$\displaystyle\tilde{T}_{2}=$
$\displaystyle\alpha\sigma_{0}-\frac{1}{2}\left(\sigma_{x}-\sqrt{3}\sigma_{y}\right),$
$\displaystyle\tilde{T}_{3}=$
$\displaystyle\alpha\sigma_{0}-\frac{1}{2}\left(\sigma_{x}+\sqrt{3}\sigma_{y}\right).$
(27)

Then, the equivalence relation set in Eqs. (8)-(10) can be reformulated using
$\gamma=\dfrac{\omega_{\textsc{AB}}}{\hbar
k_{\theta}\mathrm{v}_{\textsc{f}}}$. Turning to atomistic TB or DFT
approaches, the treatment of relaxation effects in the search algorithm for
the magic distance requires a more careful scrutiny.