#  Magnetic properties of low-angle twisted bilayer graphene at three-quarters
filling

Kevin J. U. Vidarte  Instituto de Física, Universidade Federal Rio de Janeiro,
21941-972 Rio de Janeiro, RJ, Brazil  Caio Lewenkopf  Instituto de Física,
Universidade Federal Fluminense, 24210-346 Niterói, RJ, Brazil

###### Abstract

We present a theoretical investigation of the magnetic properties exhibited by
twisted bilayer graphene systems with small twist angles, where the appearance
of flat minibands strongly enhances electron-electron interactions effects. We
consider a tight-binding Hamiltonian combined with a Hubbard mean-field
interaction term and employ a real-space recursion technique to self-
consistently calculate the system’s local density of states. The $O({\cal N})$
efficiency of the recursion method makes it possible to address superlattices
of very large size by means of a full real-space analysis. Our procedure is
validated by comparison with literature momentum-space calculations that find
a magnetic phase in charge-neutral twisted bilayer graphene. We use our method
to investigate the properties of small-angle twisted bilayer graphene at
three-quarters filling of the conduction miniband. Our calculations indicate
the emergence of a ferromagnetic phase that can be understood in terms of the
Stoner mechanism, in line with recent experimental observations.

##  I Introduction

Twisting two layers of graphene by a specific angle has led to the discovery
of a wealth of unexpected phenomena [1, 2]. Particularly intriguing is the
experimental discovery of insulating phases and superconducting states [3, 4]
in small-angle twisted bilayer graphene (TBG) systems has opened new paths for
the investigation of graphene systems [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15],
playing a key role in a line of investigation called twistronics [16], and has
rapidly attracted a large and active community of researchers in condensed
matter physics and materials science [2, 17].

The origin of these remarkable phenomena in low-angle TBG is attributed to the
appearance of flat electronic bands near the charge neutrality point (CNT)
[18, 19], as explained by the continuum model [20, 21, 22, 23]. The
corresponding electronic states are localized, leading to enhanced electron-
electron interaction effects, which are believed to originate the rich
properties observed in low-angle TBG. From a theoretical point of view, the
challenge is to treat a system with a unit cell that has thousands of atoms
with strongly interacting electrons.

The focus of our paper are the magnetic properties of low-angle TBG systems at
finite doping. The motivation is a recent experiment [24] that observed the
emergence of a ferromagnetic phase in a TBG with a twist angle $\theta\approx
1.16^{\circ}$ at 3/4 filling of the conduction miniband and the debate about
the underlying mechanism responsible for this experimental finding [24, 25].

The onset of magnetism in graphene-based systems has been intensively studied,
both experimentally and theoretically [26]. The latter span a wide variety of
settings, such as zero-dimensional graphene nanofragments [27, 28, 29], one-
dimensional graphene nanoribbons [30, 31, 32, 33, 34, 35], defect-induced
magnetism in bulk graphene [36, 37]. The description of magnetic moments
induced by edge terminations and vacancies has been addressed by both ab
initio calculations and the tight-binding approximation with an on-site
Hubbard electron-electron interaction term. The magnetic properties of these
systems are nicely described by a mean-field theory, with the exception of
vacancy-induced Kondo correlations [38, 39, 40].

In this paper, we employ a tight-binding Hamiltonian with a mean-field Hubbard
on-site interaction term to compute the low-energy electronic and magnetic
properties of low-angle TGBs. We use the Haydock-Heine-Kelly (HHK) recursive
technique [41, 42, 43, 44] to calculate the spin-resolved LDOS with spin
$\sigma=\\{\uparrow,\downarrow\\}$ of TBG as a function of the bilayer twist
angle. Being an order $\mathcal{N}$ method, the HHK method makes possible to
compute the Green’s functions of TBG with very large primitive unit cells. The
latter combined with a self-consistent mean-field calculation allows us to
investigate the electron localization properties and the magnetization of TBGs
at arbitrary filling factors of the flat minibands. Motivated by a recent
experimental paper [24], our focus is the study of emergent ferromagnetism
close to $3/4$-filling at the conduction miniband of a TBG with $\theta\approx
1.16^{\circ}$ [24].

The paper is organized as follows. In Sec. II we review the geometric
properties of commensurate TBG moiré structures, present the Hamiltonian
model, and the numerical method developed for this study. In Sec. III we
discuss the LDOS of low-angle TBG systems. First, we show our results for non-
interacting electrons. Next, we calculate the magnetic properties of TBGs at
the charge neutrality point (CNP). The excellent agreement between our results
and those of Ref. [45] serves to validate our method. Finally, we investigate
the magnetic properties of a TBG system with $\theta=1.16^{\circ}$ at 3/4
filling and find the emergence of a ferromagnetic phase, in line with recent
experiments [24]. We present our conclusion and an outline of future research
in Sec. IV.

##  II Theory and methods

###  II.1 TBG geometric and symmetry properties

The stacking geometry of TBGs is characterized, as standard [20, 46, 47], by
the twist of one graphene layer with respect to the other around a given site
starting from the AA-stacked bilayer, forming a moiré pattern. We describe a
TBG system by twisting the upper layer at an angle $\theta$ as described in
Ref. [20, 21]. The primitive lattice vectors of the bottom layer are
$\textbf{a}^{b}_{1}=\sqrt{3}a_{0}\hat{\bf e}_{x}$ and
$\textbf{a}^{b}_{2}=\sqrt{3}a_{0}/2\left(\hat{\bf e}_{x}+\sqrt{3}\hat{\bf
e}_{y}\right)$ with the carbon-carbon bond length $a_{0}=1.42$ Å, and those of
top layer as $\textbf{a}^{t}_{i}=R(\theta)\textbf{a}^{b}_{i}$, where
$R(\theta)$ is the rotation matrix. In this study, we take the spacing between
graphene layers as $d_{0}=3.35$ Å [19] and do not consider lattice relaxation
effects.

The lattice structure of a TBG system is periodic if the periods of the two
graphene layers match, giving a finite unit cell. Hence, the periodicity
condition requires the lattice translation vector
$m\textbf{a}_{1}^{b}+n\textbf{a}_{2}^{b}$ of the bottom (unrotated) layer and
$n\textbf{a}_{1}^{t}+m\textbf{a}_{2}^{t}$ in the top (rotated) layer, with $m$
and $n$ integers, to coincide. Hence, the twist angle for a commensurate
structure is related to $(m,n)$ by [48, 49, 50]

$\theta(m,n)=\arccos{\left(\dfrac{1}{2}\dfrac{m^{2}+n^{2}+4mn}{m^{2}+n^{2}+mn}\right)},$
(1)

with a lattice constant

$L=a_{0}\sqrt{3(m^{2}+n^{2}+mn)}=\dfrac{|m-n|\sqrt{3}a_{0}}{2\sin{\theta/2}}.$
(2)

The commensurate unit cell contains $N=4(m^{2}+n^{2}+mn)$ atoms.

Due to the symmetry of the honeycomb lattice, when the translation vector of
the top layer is fixed to $\boldsymbol{\delta}^{t}=0$, a twist at an angle
$\theta=120^{\circ}$ transforms the AA-stacked bilayer into itself. In turn, a
twist by $\theta=60^{\circ}$ transforms the bilayer from AA- to AB-stacking.
TBGs with $\theta$ and $-\theta$ are mirror images that share equivalent band
structures [48].

Each graphene layer is constituted by two sublattices, ${\rm A}^{b({\rm
or}\;t)}$ and ${\rm B}^{b({\rm or}\;t)}$. The commensurate structures occur in
two different forms distinguished by their sublattice parities [49, 50].
Figure 1(a) shows an example of an even commensurate structure under
sublattice exchange. It is characterized by having 3 three-fold symmetric
positions that correspond to the stacking of the ${\rm A}^{b}{\rm A}^{t}$ and
${\rm B}^{b}{\rm B}^{t}$ sites and by hexagonal centers that overlap, ${\rm
H}^{b}{\rm H}^{t}$. In contrast, Fig. 1(b) shows an odd commensurate
structure. Here, the top and bottom coinciding sites correspond to ${\rm
A}^{b}{\rm A}^{t}$ at the origin, and the two remaining three-fold symmetric
positions are occupied by ${\rm B}^{b}$(or ${\rm B}^{t}$)-sublattice sites of
one layer aligned with the hexagon centers ${\rm H}^{t}$(or ${\rm H}^{b}$) of
its neighbor layer. Therefore, the angles $60^{\circ}-\theta$ and $-\theta$
followed by a relative translation of the upper layer by
$\boldsymbol{\delta}^{t}=(\textbf{a}_{1}^{t}+\textbf{a}_{2}^{t})/3$ form
commensurate partners with unit cells of equal areas but opposite sublattice
parities.

Figure 1: Commensurate structure partners with $\theta\approx 21.8^{\circ}$
that corresponds $(1,2)$. Purple (orange) circles indicate the carbon sites of
the top (bottom) layer. The black lines indicate the primitive unit cell. The
3 blue disks correspond to the three-fold symmetry positions. (a) The twist
angle $60^{\circ}-\theta$ corresponds to an even symmetry under sublattice
exchange at three-fold symmetric positions. (b) The twist angle $-\theta$
followed by $\boldsymbol{\delta}^{t}$ translation represents an odd symmetry
under sublattice exchange. The green disks indicate the equivalent points. The
green dashed line corresponds to the two-fold rotation axis.

A TBG with $\theta=30^{\circ}$ is a special case, since its crystal structure
is its own commensurate partner and corresponds to an elementary two-
dimensional quasicrystalline lattice [51, 52].

Here, we work with odd commensurate structures. In Fig. 1(b), the carbon atom
sites, indicated by the green disks, are equivalent due to the three-fold
symmetric positions and two-fold rotation axis (green dashed line). These
symmetry properties allow us to reduce the numerical computation by
identifying the equivalent sites within the primitive unit cell.

###  II.2 Model Hamiltonian

We model the electronic properties of low-angle TBG systems using a tight-
binding Hamiltonian with a Hubbard interaction term [26, 53], namely

$H=H_{\rm TB}+H_{\rm U},$ (3)

where $H_{\rm TB}$ stands for the TBG tight-binding Hamiltonian and $H_{\rm
U}$ for a Hubbard on-site electron-electron interaction term.

The tight-binding Hamiltonian that describes the low-energy electronic
structure of TBG reads

$H_{\rm TB}=\sum_{i,j\sigma}\left(t_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}+{\rm
H.c}\right),$ (4)

where the operators $c_{i\sigma}$ and $c^{\dagger}_{i\sigma}$ annihilate and
create an electron with spin $\sigma=\left\\{\uparrow,\downarrow\right\\}$ at
site $i$, respectively. $t_{ij}$ is the transfer integral between the Wannier
electronic orbitals centered at the carbon sites $i$ and $j$. The transfer
integral $t_{ij}$ depends on the interatomic distance and the relative
orientation between $p_{z}$ orbitals at each site and is usually parameterized
as [19, 48, 54]

$t_{ij}=V_{pp\pi}\left[1-\left(\dfrac{d_{0}}{d}\right)^{2}\right]+V_{pp\sigma}\left(\dfrac{d_{0}}{d}\right)^{2},$
(5)

with

$\displaystyle
V_{pp\pi}=V^{0}_{pp\pi}e^{-(d-a_{0})/\delta},\;\;\;V_{pp\sigma}=V^{0}_{pp\sigma}e^{-(d-a_{0})/\delta},$
(6)

where $d=|\textbf{r}_{i}-\textbf{r}_{j}|^{1/2}$. Here, $V^{0}_{pp\pi}=-2.7$ eV
is the transfer integral between the nearest-neighbor atoms in the monolayer
graphene, $V^{0}_{pp\sigma}=0.48$ eV is the interlayer transfer integral
between vertically located atoms, and $\delta=0.319a_{0}$ is the decay length
of the transfer integral [54, 19]. The transfer integral for $d>5$ Å is
exponentially small and can be safely neglected. The intra- and interlayer
hopping matrix elements $t_{ij}$ have been determined by fitting the
dispersion relations of graphene monolayer and graphene AB-stacked bilayer
obtained from ab initio calculations [19].

The Hubbard term reads

$H_{\rm U}=U\sum_{i}n_{i\uparrow}n_{i\downarrow},$ (7)

where $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ is the spin-resolved
electron number operator at site $i$. The parameter $U>0$ gives the magnitude
of the on-site Coulomb repulsion [53, 26]. The magnitude of $U$ in graphene
systems is has been extensively discussed [55, 56] with estimates that vary
from $0.5\cdots 2.0V^{0}_{pp\pi}$. In this study we consider, unless otherwise
stated, $U=V^{0}_{pp\pi}$.

We solve the TBG Hamiltonian for a given filling factor in the mean-field
approximation. As standard, we write $n_{i\sigma}\equiv\langle
n_{i\sigma}\rangle+\delta n_{i\sigma}$ and neglect the quadratic terms in
$\delta n_{i\sigma}$ that are responsible for electronic correlations. The
mean field Hubbard Hamiltonian reads

$H_{\rm U}^{\rm MF}=U\sum_{i}\left[n_{i\uparrow}\langle
n_{i\downarrow}\rangle+n_{i\downarrow}\langle n_{i\uparrow}\rangle-\langle
n_{i\uparrow}\rangle\langle n_{i\downarrow}\rangle\right].$ (8)

The last term is a constant for a given local magnetic configuration and can
be absorbed in the chemical potential.

###  II.3 Numerical method

Here, we describe how we calculate the ground state charge density and
magnetic properties of TBG systems modeled by the Hamiltonian $H=H_{\rm
TB}+H_{\rm U}^{\rm MF}$.

Due to the large number of atoms in a low-angle TBG primitive cell, even mean-
field calculations can be computationally very costly. Several effective
Hamiltonian models have been developed to reduce the numerical effort [57, 45,
58, 59]. Here, we perform a full real-space calculation using the HHK
recursive method, which is very efficient to compute the single-particle LDOS
of large systems. The HHK method [42, 43, 44, 41] puts forward a very
efficient Lanczos-like $O({\cal N})$ recursive procedure that transforms an
arbitrary sparse Hamiltonian matrix in a tridiagonal one. Next, it evaluates
the diagonal Green’s function in real space by a continued fraction expansion,
which is much more numerically amenable than a full diagonalization.

By a suitable choice of the chemical potential, our approach allows us to
consider charge neutral as well as systems with a finite doping, namely,
$n_{\rm dop}=N_{e}/A_{\rm PUC}$, where $N_{e}$ is the number of electrons in
excess to the CNP and $A_{\rm PUC}=3\sqrt{3}a_{0}^{2}/8\sin^{2}(\theta/2)$ is
the area of the TBG moiré cell with $N$ atoms.

We implement the self-consistent mean-field calculation as standard:

($i$) We start with an initial set of occupation numbers $\langle
n_{i\sigma}\rangle$, with the constraint $\sum_{i\sigma}\langle
n_{i\sigma}\rangle=N+N_{e}$, where the sum runs over all sites of the moiré
unit cell, set by the considered doping. The occupations can be chosen
randomly or respecting some given symmetry condition.

($ii$) Using the HHK recursion technique [41, 42, 43, 44], we compute the LDOS
of the electronic system defined by the Hamiltonian, Eq. (3). The spin-
resolved LDOS at a given site $i$ can be written as

$\displaystyle\nu_{j}(\epsilon)$
$\displaystyle=-\frac{1}{\pi}\lim_{\eta\rightarrow 0^{+}}{\rm
Im}\;G_{jj}(\epsilon+i\eta).$ (9)

In practice, the regularization parameter $\eta$ is considered as finite and
its magnitude can be attributed the self-energy correction due to disorder,
that is ubiquitous in graphene systems. The HHK method enables us to compute
the LDOS with $O(\mathcal{N})$ operations and the DOS with
$O(\mathcal{N}^{2})$. Hence, it is much more efficient that a direct
diagonalization, which demands $O(\mathcal{N}^{3})$ operations.

($iii$) Next, we determine $\langle n_{i\sigma}\rangle$, the average electron
occupation number with spin $\sigma$ at the site $i$. In general, at a given
temperature $T$, $\langle n_{i\sigma}\rangle$ reads

$\langle n_{i\sigma}\rangle=\int_{-\infty}^{\infty}{\rm d}\epsilon
f_{\mu}(\epsilon)\nu_{i\sigma}(\epsilon),$ (10)

with
$f_{\mu}(\epsilon)=\left\\{\exp{\left[\beta(\epsilon-\mu)\right]}+1\right\\}^{-1}$,
where $\beta=1/k_{\rm B}T$, $k_{\rm B}$ is Boltzmann constant, and $\mu$ is
chemical potential. At zero absolute temperature, $\mu$ is equal to the Fermi
energy $\varepsilon_{F}$ and the Fermi distribution becomes a step function.
For simplicity, here we consider $T=0$. As standard, the Fermi energy
$\varepsilon_{F}$ (or the chemical potential $\mu$) is determined by the
doping.

($iv$) We examine whether the output occupation numbers $\langle
n_{i\sigma}\rangle$ coincide with the input ones, within a tolerance of
$10^{-6}$. If the answer is positive, we end the self-consistent loop. If not,
we return to ($ii$). The computation of the updated spin densities is then
repeated iteratively until all values of $\langle n_{i\sigma}\rangle$ are
converged.

The self-consistent solution provides the spin polarization at the $i$th site,

$\displaystyle p_{z,i}$ $\displaystyle=\dfrac{\langle
n_{i\uparrow}\rangle-\langle n_{i\downarrow}\rangle}{2}.$ (11)

The $p_{z,i}$ is translated in a local magnetization $m_{z,i}=-g\mu_{\rm
B}p_{z,i}$, where $\mu_{\rm B}$ is the Bohr magneton and the electron
$g$-factor is $g=2$. Hence, the total magnetization per moiré cell reads

$\displaystyle M_{\rm PUC}$ $\displaystyle=\sum^{N}_{i=1}m_{z,i}.$ (12)

We study two cases: $(i)$ $N_{e}=0$, corresponding to the CNP, and $(ii)$
$N_{e}=3$ corresponding to $3/4$-filling of the conduction miniband [24]. For
the CNP case, a previous theoretical study working in reciprocal space has
predicted that low-angle TBG systems have an antiferromagnetic ground state
[45]. In this case, we find computationally convenient to start the self-
consistent loop with the configuration: $\langle
n^{A}_{i\uparrow}\rangle=1/2+\Delta n^{A}$, $\langle
n^{A}_{i\downarrow}\rangle=1/2-\Delta n^{A}$, $\langle
n^{B}_{i\uparrow}\rangle=1/2+\Delta n^{B}$ and $\langle
n^{B}_{i\downarrow}\rangle=1/2-\Delta n^{B}$, where $\Delta n^{B}=-\Delta
n^{A}$. For the $3/4$-filling case, we start with a random set of $\langle
n_{i\sigma}\rangle$ with the constraint $\sum_{i\sigma}\langle
n_{i\sigma}\rangle=N+N_{e}$, where $N_{e}=3$.

##  III Results

We begin this section by presenting a study of the single-particle LDOS of TBG
systems with a focus on the formation of low-energy minibands at small twist
angles $\theta$. Next, we use the obtained LDOS to calculate the magnetization
of low-angle TBG systems using the self-consistent procedure outlined in Sec.
II.3. We address on two cases: ($i$) TBGs at the CNP and ($ii$) TBGs at
$3/4$-filling of the conduction band.

###  III.1 Non-interacting electronic densities

We compute the local electronic properties of TBGs using the HHK recursion
method [41, 42, 43, 44]. Being a real space approach, the HHK calculations
take advantage of the symmetry properties of the odd commensurate TBG
structures discussed above. We have numerically verified that the predicted
equivalent sites indeed display the same LDOS. This is used to reduce the
numerical effort by a factor of 6.

The regularization factor is set to be at most $\eta\approx 25$ meV, a value
that can be attributed to a small disorder concentration. We return to this
point later on.

To highlight the sites with the most prominent enhancement of the LDOS in
TBGs, we consider the moiré region where the AA stacking is placed at the
center of the TBG primitive unit cell. Figure 2(a) shows the LDOS for a moiré
unit cell with $(m,n)=(22,23)$ corresponding to a twist angle $\theta\approx
1.47^{\circ}$ and containing 6076 carbon atoms. The AB- (or BA-) stacking
region is zoomed in the inset. Close to the CNP, well-localized states are
found in the AA-stacking region, leading to an enhanced LDOS on atoms around
the AA stacking, with much smaller LDOS at AB- and BA-stacking regions, in
agreement with previous studies [6, 19].

Figure 2: Non-interacting electronic densities: (a) LDOS at the CNP,
$\nu_{i}({\rm CNP})$ (in arbitrary units), for the $(m,n)=(22,23)$
commensurate lattice, that corresponds to $\theta\approx 1.47^{\circ}$ and
6076 atoms within the moiré unit cell (black line). The areas of the purple
disks are proportional to $\nu_{i}({\rm CNP})$. (b) LDOS at the AA dimer
sites, $\nu_{\rm AA}$, as a function of the energy (in electron volts) for odd
commensurate $n-m=1$ TBG systems with $\theta\approx 1.61^{\circ}$,
$1.02^{\circ}$ and $0.82^{\circ}$. (c) $\nu_{\rm AA}$ at the CNP as a function
of $\theta$ (in degrees) corresponding to odd commensurate moiré structures
with $n-m=1$.

Figure 2(b) shows the LDOS at the AA dimer site as a function of the energy
$\varepsilon$ around the CNP Fermi energy for $\theta\approx 1.16^{\circ}$,
$1.02^{\circ}$ and $0.82^{\circ}$, that belong to the family of odd
commensurate lattices with $n-m=1$. The appearance of the central peak can be
interpreted in terms of the continuum model [21], which associates the
enhancement of the LDOS at the vicinity of the magic angle to the formation of
a flat band at the CNP. Figure 2(c) displays the evolution of the AA dimer
site LDOS at the CNP as a function of the small twist angles $\theta$ for the
$n-m=1$ family of odd commensurate lattice structures.

Figures 2(b) and (c) show that when the twist angle decreases, the LDOS of the
AA-stacking region increases significantly at the CNP. The maximum of the peak
appears at the twist angle of $\theta\approx 1.16^{\circ}$ corresponding to
$(m,n)=(28,29)$ with 9748 carbon atoms within the moiré cell. For the $n-m=1$
structures considered here, this twist angle is the closest to the magic angle
$\theta\approx 1.1^{\circ}$ predicted by the continuum model [21, 23].

As the twist angle decreases ($\theta<1.16^{\circ}$), the height of the
central DOS peak becomes smaller, but the satellite peaks around the main one,
both in the valence and conduction bands, see Fig. 2(b), become increasingly
intense, and move closer to the CNP. Eventually, these peaks begin to overlap,
showing another peak maximum at the twist angle close to $0.53^{\circ}$, see
Fig. 2(c), corresponding to $(m,n)=(62,63)$ with $46876$ carbon atoms within
the moiré cell. This behavior is in good agreement with theoretical
predictions based on the continuum model [21].

We note that in Fig. 2(c) we compute $\nu_{i}(\epsilon)$ using $\eta=5$ meV.
This choice guarantees that $\eta$ is smaller than the flat-band broadening in
a magic angle twisted bilayer graphene (MATBG) that is about 10 meV [4, 21]
and helps to accurately determine the peak maxima at the CNP as a function of
the twist angle $\theta$. In contrast, in Figs. 2(a)-(b) we use $\eta=25$ meV.
For this parameter choice, the flat-band broadening is greater than 10 meV,
but we stress that this does not change the possibility of filling just the
flat band with 4 electrons per primitive unit cell [4, 3, 7, 5], since the
satellite peaks are separated from the central one by gaps of approximately
0.17 eV. Considering that the physical origin of a finite $\eta$ can be
attributed to disorder, the latter can be the reason for the lack of
experimental observations of magic angles $\theta_{n}=1.1^{\circ}/n$ with
$n>1$.

Let us now depart from the CNP and discuss finite doping. Figure 3(a) shows
the DOS of the single-particle miniband of a TBG with $\theta\approx
1.47^{\circ}$. Anticipating the appearance of a gap at the CNP when the
interaction is switched on, we define a valence and a conduction miniband
separated at the CNP. The red area represents the half-band filling of a TBG
system where the Fermi energy is at the CNP. This filling corresponds to 4
electrons per moiré cell in the valence miniband. The green area represents
the electron doping at the $3/4$ filling of the conduction miniband. This
doping corresponds to 3 electrons per moiré cell in excess of the CNP at a
Fermi energy $\varepsilon_{F}=872$ meV.

Figure 3: Electron doping of the conduction miniband: (a) DOS (in arbitrary
units) of the flat miniband as a function of the energy (in electron volts).
The CNP corresponds to the DOS peak. The green area represents 3/4-filling of
the conduction miniband. (b) Carrier distribution within the moiré cell (black
line) in real space for $\theta\approx 1.47^{\circ}$. The areas of the green
disks are proportional to $\Delta n_{i}$, the occupation in excess to the CNP
filling.

Figure 3(b) shows the carrier distribution $\Delta n_{i}$ in excess to the CNP
distribution within the moiré unit cell for $\theta\approx 1.47^{\circ}$. The
areas of the green disks are proportional to the carrier occupation number at
the atomic site $i$, $\Delta n_{i}$. Owing to the localized LDOS, see Fig. 2,
the carrier distribution is also concentrated at the AA-stacking region.

###  III.2 Magnetism at the CNP

Using the mean-field Hubbard Hamiltonian, we implement the self-consistent
scheme described in Sec. II.3 to study the emergence of antiferromagnetic
states in low twist angle TBGs at the CNP.

Figure 4(a) displays the local spin polarization, $p_{z,i}$, within the moiré
unit cell for a $(m,n)=(22,23)$ lattice corresponding to a rotation angle
$\theta\approx 1.47^{\circ}$ with $U/V^{0}_{pp\pi}=1$. The regions with the
largest local magnetic moments correspond to those of enhanced LDOS, see Fig.
2, strongly suggesting that the emergence of magnetism can be attributed to
the Stoner mechanism [53]. In the AA stacking region, the imbalance between
$\langle n_{i\uparrow}\rangle$ and $\langle n_{i\downarrow}\rangle$ leads to
the emergence of an antiferromagnetic ground state at the CNP, as can be seen
in the zoom. We note in passing that the converged solution is the same if the
we start the self-consistent loop with a distribution of random or
antiferromagnetic occupations $\langle n_{i\sigma}\rangle$.

Figure 4: Antiferromagnetic phase at the CNP: (a) Local spin polarization
$p_{z,i}$ within the moiré cell (black line) for a TBG with $\theta\approx
1.47^{\circ}$ and $U/V^{0}_{pp\pi}=1$. The disk areas are proportional to
$|p_{z,i}|$, where the red and blue disks correspond to positive and negative
spin polarizations, respectively. (b) Local spin polarization $p_{z,i}$ as a
function of the radial distance with respect to the AA dimer site (or AB non-
dimer site) for $U/V^{0}_{pp\pi}=1.2$, $1.4$ and $1.5$ with $\theta\approx
1.47^{\circ}$.

Figure 4(b) shows the local spin polarization for a TBG system with
$\theta\approx 1.47^{\circ}$ for different values of the on-site Coulomb
interaction, $U/V^{0}_{pp\pi}=1.2$, $1.4$ and $1.5$, as a function of the
radial distance to the AA dimer site. Here, the distance between the AA dimer
site and the AB nondimer one is $D=55.3$ Å. In line with Fig. 4(a), one
observes that the spin polarization shows a maximum at the AA-stacking region,
and gradually decreases and eventually vanishes for the sites that are closer
to the AB stacking. The magnitude of local spin polarization becomes
increasingly larger with increasing on-site Coulomb interaction $U$, and the
overall dependence with the radial distance to the AA dimer site is similar
for all $U$ values considered. Note that since the ground state is
antiferromagnetic, the total spin polarization per moiré cell is zero for all
values of $U$. Figure 4 shows excellent agreement with the reciprocal space
calculation reported in Ref. [45].

Figure 5 displays the LDOS of the flat miniband calculated for AA dimer site,
namely, ${\rm LDOS}({\rm AA},\epsilon)=\nu_{{\rm
AA},\uparrow}(\epsilon)+\nu_{{\rm AA},\downarrow}(\epsilon)$ for a twist angle
of $\theta\approx 1.47^{\circ}$ with $U/V^{0}_{pp\pi}=1$. Due to the electron-
electron interaction, the LDOS of the flat miniband in the AA dimer site is
split into two peaks almost symmetric around the Fermi energy, here set to
$\epsilon_{F}=0$, and a small gap is opened. The separation between the two
LDOS peaks is approximately $0.031$ eV. Note that since the charge density in
not homogeneous, one expects ${\rm LDOS}({\rm AA},\epsilon)$ to differ from
the DOS($\varepsilon$). However, since the two low-energy LDOS peaks are
absent in the atomic sites at the AB-stacking region, the ${\rm LDOS}({\rm
AA},\epsilon)$ captures the main features of ${\rm DOS}(\epsilon)$ at low
energies.

Figure 5:  LDOS at the AA dimer site, $\nu_{{\rm AA},\uparrow}+\nu_{{\rm
AA},\downarrow}$ (in arbitrary units) as a function of energy $\varepsilon$
(in eV) for $\theta\approx 1.47^{\circ}$ with $U/V^{0}_{pp\pi}=1$. For
convenience, the Fermi energy at the CNP is set to zero.

Figure 6(a) shows the maximum spin polarization that corresponds to the AA
dimer site, $p_{z,{\rm AA}}$, as a function of the on-site Coulomb repulsion
$U/V^{0}_{pp\pi}$ for different twist angles $\theta\approx 1.47^{\circ}$,
$1.30^{\circ}$ and $1.08^{\circ}$. We observe that the local spin polarization
at the AA dimer site increases monotonically with increasing $U$. Furthermore,
$p_{z,{\rm AA}}$ versus $U/V^{0}_{pp\pi}$ shows a similar (almost linear)
slope for $\theta\approx 1.08^{\circ}$ and $\theta\approx 1.30^{\circ}$, while
for $\theta\approx 1.47^{\circ}$ both $p_{z,{\rm AA}}$ and its slope with
respect to $U/V^{0}_{pp\pi}$ are much smaller. This suggests that
$\theta\approx 1.47^{\circ}$ is close to a critical angle where the magnetic
phase ceases to exist. A study of the critical values of $U_{c}/V^{0}_{pp\pi}$
for the zero temperature normal antiferromagnetic transition of TBG systems at
the CNP can be found in Ref. [45], which estimates the critical values of
$U_{c}/V^{0}_{pp\pi}$.

Figure 6:  (a) Local spin polarization at the AA dimer site versus
$U/V^{0}_{pp\pi}$ for $\theta\approx 1.47^{\circ}$, $1.30^{\circ}$ and
$1.08^{\circ}$. (b) Local spin polarization at the AA dimer site as a function
of the twist angle $\theta$ (in degrees) for $U/V^{0}_{pp\pi}=1.0$ and $1.5$.

Figure 6(b) shows the maximum spin polarization, $p_{z,{\rm AA}}$, as a
function of the magic angle $\theta$ for $U/V^{0}_{pp\pi}=1.0$ and $1.5$. We
note that $p_{z,{\rm AA}}$ shows a maximum close to $\theta=1.30^{\circ}$ for
different values of $U/V^{0}_{pp\pi}$, consistent with Fig. 6(a). This value
of $\theta$ is slightly shifted with respect to the magic angle, $\theta_{\rm
MA}\approx 1.1^{\circ}$, which can be attributed to small differences in the
corresponding single-particle model Hamiltonian parameterizations.

###  III.3 Magnetism at $3/4$-filling of the conduction miniband

Let us now examine the ferromagnetic phase in a TBG system with a twist angle
$\theta\approx 1.16^{\circ}$ at $3/4$-filling of the conduction miniband whose
origin has been the subject of discussion in the literature.

We consider an odd commensurate lattice with $(m,n)=(28,29)$ corresponding to
9748 carbon atoms. Our results are obtained by following the self-consistent
procedure described in Sec. II.3 for a finite doping. Here, we set
$U/V^{0}_{pp\pi}=1$.

Figure 7(a) displays the local spin polarization, $p_{z,i}$, within a circular
area with a radius of $21$ Å centered at the AA dimer site. These results show
that the ground state is ferromagnetic as characterized by the imbalance
between $\langle n_{i\uparrow}\rangle$ and $\langle n_{i\downarrow}\rangle$ in
the AA stacking region. Since the spin polarization is largest at the regions
of enhanced LDOS, we attribute the emergence of ferromagnetism to the Stoner
instability mechanism.

Figure 7: Ferromagnetic phase at $3/4$-filling of the conduction miniband: (a)
Local spin polarization of a TBG with $\theta\approx 1.16^{\circ}$ and
$U/V^{0}_{pp\pi}=1$ within the AA-stacking region. The areas of the red disks
are proportional to the spin polarization $p_{z,i}$. (b) Local spin
polarization $p_{z,i}$ versus the distance $D$ from the $i$th site to the AA
dimer site up to the AB nondimer site.

Figure 7(b) shows the local spin polarization, $p_{z,i}$, as a function of the
distance $D$ from the AA dimer site. Similarly to the CNP case, the spin
polarization is largest at the AA dimer site and becomes smaller, as $D$
increases and eventually vanishes as the $i$th site approaches the AB nondimer
position. On the other hand, in distinction to the previous case, since the
ground state is ferromagnetic, the total magnetization is not zero. We find
that the total spin polarization per moiré cell is $\sim 0.19$ and $M_{\rm
PUC}=-0.22$ ${\rm eV}\cdot{\rm T}^{-1}$.

Figure 8(a) presents the LDOS calculated at the AA dimer site, ${\rm
LDOS}({\rm AA},\epsilon)$, within an energy window that contains the flat
minibands near the Fermi energy (set to $\varepsilon_{F}=0$ for clarity).
Similar to the previous case, we argue that the LDOS(AA,$\varepsilon$)
reproduces the behavior of DOS($\varepsilon$) in the vicinity of
$\varepsilon_{F}$. The separation of the two LDOS peaks around the Fermi
energy is approximately 0.029 eV. The existence of two LDOS peaks in the
valence miniband suggests the existence of two minibands that overlap in the
CNP, as indicated in Ref. [3].

Figure 8:  (a) LDOS in the AA dimer site, $\nu_{{\rm AA},\uparrow}+\nu_{{\rm
AA},\downarrow}$ (in arbitrary units) as a function of the energy (in eV). (b)
LDOS for spin up $\nu_{{\rm AA},\downarrow}$ (red line) and spin down
$\nu_{{\rm AA},\uparrow}$ (blue line) at the AA dimer site (in arbitrary
units) as a function of the energy (in eV). Both (a) and (b) correspond to
$\theta\approx 1.16^{\circ}$ with $U/V^{0}_{pp\pi}=1$. For convenience, the
Fermi energy at $3/4$-filling of the conduction miniband is set to zero.

Figure 8(b) depicts the spin up and spin down LDOS at the AA dimer site,
denoted as $\nu_{{\rm AA},\uparrow}(\epsilon)$ and $\nu_{{\rm
AA},\downarrow}(\epsilon)$, respectively. A clear spin-imbalance is evident
between the spin-up and spin-down occupation numbers in the filled miniband
(fmb). Quantitatively, the occupation numbers are $\langle n^{\rm fmb}_{{\rm
AA},\uparrow}\rangle=9.7\times 10^{-4}$ and $\langle n^{\rm fmb}_{{\rm
AA},\downarrow}\rangle=6.1\times 10^{-4}$. This imbalance translates to a spin
polarization $p_{z,{\rm AA}}=1.8\times 10^{-4}$ at the AA dimer site,
corresponding to a local magnetic moment $m_{z,{\rm AA}}=-0.21$ ${\rm
meV}\cdot{\rm T}^{-1}$.

##  IV Conclusions

In this work, we investigate the emergence of magnetism in low-angle TBG
systems using a numerical real-space approach. We have considered a tight-
binding Hamiltonian with a mean-field Hubbard term and obtained the ground
state of TBG systems by a self-consistent iteration procedure. Notably, owing
to the HHK recursive technique [41, 42, 43, 44], our approach allows
calculations to be efficiently performed for very large moiré cells.

To validate our method, we compare our results for low-angle TBGs at the CNP
with those previously reported in the literature [45]. Specifically, we have
found the emergence of an antiferromagnetic phase for low-angle TBGs with a
maximum local spin polarization $p_{z,i}$ near the magic angle. For all
computed twist angles $\theta$, the agreement of $p_{z,i}$ with Ref. [45] is
excellent, demonstrating the accuracy of our computational procedure.

The main finding of this study is the emergence of a ferromagnetic phase in
low-angle TBGs at $3/4$-filling of the conduction miniband. Motivated by
recent experiments [24, 60] we have investigated the DOS, the local and total
magnetization of TBGs with $\theta\approx 1.16^{\circ}$. Our calculations for
the case of $3/4$-filling of the conduction miniband (3 electrons in excess to
the CNP per moiré cell) show that the system ground state is a ferromagnetic
insulator with a gap of $\Delta\approx 0.029$ eV, in agreement with the
experimental findings [24]. We attribute this behavior to the large LDOS at
the AA-stacking region, as predicted by the continuum model [21, 23], that
causes a strong enhancement of the electron-electron interaction, a key
element for the Stoner mechanism. Our results provide an alternative scenario
to the substrate-induced gap proposed in Ref. [25].

We expect our methodology to be useful for the study of interaction effects of
other moiré 2D systems with large primitive unit cells. For instance, with
modest adaptations, our approach can be used in the analysis of anisotropic
ferromagnetism dominated by the orbital magnetic moment in 3/4-filling at the
conduction miniband of MATBG systems [60]. By combining our calculations with
with a molecular dynamics procedure, we plan to analyze lattice relaxation
effects as a function of twist angle. In addition, the presence of an external
magnetic field, that typically requires the use of large supercells in
standard approaches, in our method is trivial to account.

###### Acknowledgements.

We thank J. Vahedi for the discussions and for providing the data of his
computations of TBG systems. This work was partially supported by the
Brazilian Institute of Science and Technology (INCT) in Carbon Nanomaterials
and the Brazilian agencies CAPES, CNPq, FAPEMIG, and FAPERJ.

## References

  * MacDonald [2019] A. H. MacDonald, Bilayer graphene’s wicked, twisted road, Physics 12, 12 (2019).
  * Andrei and MacDonald [2020] E. Andrei and A. MacDonald, Graphene bilayers with a twist, Nat. Mater. 19, 1265–1275 (2020).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] 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 (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).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Maximized electron interactions at the magic angle in twisted bilayer graphene, Nature 572, 95 (2019).
  * 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, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Serlin  _et al._ [2020] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, Intrinsic quantized anomalous Hall effect in a moiré; heterostructure, Science 367, 900 (2020).
  * Saito  _et al._ [2021] Y. Saito, F. Yang, J. Ge, X. Liu, T. Taniguchi, K. Watanabe, J. I. A. Li, E. Berg, and A. F. Young, Isospin pomeranchuk effect in twisted bilayer graphene, Nature 592, 220 (2021).
  * Utama  _et al._ [2021] M. I. B. Utama, R. J. Koch, K. Lee, N. Leconte, H. Li, S. Zhao, L. Jiang, J. Zhu, K. Watanabe, T. Taniguchi, P. D. Ashby, A. Weber-Bargioni, A. Zettl, C. Jozwiak, J. Jung, E. Rotenberg, A. Bostwick, and F. Wang, Visualization of the flat electronic band in twisted bilayer graphene near the magic angle twist, Nat. Phys. 17, 184 (2021).
  * Lisi  _et al._ [2021] S. Lisi, X. Lu, T. Benschop, T. A. de Jong, P. Stepanov, J. R. Duran, F. Margot, I. Cucchi, E. Cappelli, A. Hunter, A. Tamai, V. Kandyba, A. Giampietri, A. Barinov, J. Jobst, V. Stalman, M. Leeuwenhoek, K. Watanabe, T. Taniguchi, L. Rademaker, S. J. van der Molen, M. P. Allan, D. K. Efetov, and F. Baumberger, Observation of flat bands in twisted bilayer graphene, Nat. Phys. 17, 189 (2021).
  * Tseng  _et al._ [2022] C.-C. Tseng, X. Ma, Z. Liu, K. Watanabe, T. Taniguchi, J.-H. Chu, and M. Yankowitz, Anomalous Hall effect at half filling in twisted bilayer graphene, Nat. Phys. 18, 1745 (2022).
  * Liu  _et al._ [2022] C. Liu, Z. Li, R. Qiao, Q. Wang, Z. Zhang, F. Liu, Z. Zhou, N. Shang, H. Fang, M. Wang, Z. Liu, Z. Feng, Y. Cheng, H. Wu, D. Gong, S. Liu, Z. Zhang, D. Zou, Y. Fu, J. He, H. Hong, M. Wu, P. Gao, P.-H. Tan, X. Wang, D. Yu, E. Wang, Z.-J. Wang, and K. Liu, Designed growth of large bilayer graphene with arbitrary twist angles, Nat. Mater. 21, 1263 (2022).
  * Lin  _et al._ [2022] J.-X. Lin, Y.-H. Zhang, E. Morissette, Z. Wang, S. Liu, D. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, and J. I. A. Li, Spin-orbit–driven ferromagnetism at half moiré filling in magic-angle twisted bilayer graphene, Science 375, 437 (2022).
  * Jaoui  _et al._ [2022] A. Jaoui, I. Das, G. Di Battista, J. Díez-Mérida, X. Lu, K. Watanabe, T. Taniguchi, H. Ishizuka, L. Levitov, and D. K. Efetov, Quantum critical behaviour in magic-angle twisted bilayer graphene, Nat. Phys. 18, 633 (2022).
  * 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).
  * Wang  _et al._ [2019] J. Wang, X. Mu, L. Wang, and M. Sun, Properties and applications of new superlattice: twisted bilayer graphene, Mater. Today Phys. 9, 100099 (2019).
  * 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).
  * Trambly de Laissardière  _et al._ [2010] G. Trambly de Laissardière, D. Mayou, and L. Magaud, Localization of Dirac electrons in rotated graphene bilayers, Nano Lett. 10, 804 (2010).
  * 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).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proc. Natl. Acad. Sci. U.S.A. 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).
  * 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).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Zhang  _et al._ [2019] Y.-H. Zhang, D. Mao, and T. Senthil, Twisted bilayer graphene aligned with hexagonal boron nitride: Anomalous Hall effect and a lattice model, Phys. Rev. Res. 1, 033126 (2019).
  * Yazyev [2010] O. V. Yazyev, Emergence of magnetism in graphene materials and nanostructures, Rep. Prog. Phys. 73, 056501 (2010).
  * Fernández-Rossier and Palacios [2007] J. Fernández-Rossier and J. J. Palacios, Magnetism in graphene nanoislands, Phys. Rev. Lett. 99, 177204 (2007).
  * Wang  _et al._ [2009] W. L. Wang, O. V. Yazyev, S. Meng, and E. Kaxiras, Topological frustration in graphene nanoflakes: Magnetic order and spin logic devices, Phys. Rev. Lett. 102, 157201 (2009).
  * Feldner  _et al._ [2010] H. Feldner, Z. Y. Meng, A. Honecker, D. Cabra, S. Wessel, and F. F. Assaad, Magnetism of finite graphene samples: Mean-field theory compared with exact diagonalization and quantum Monte Carlo simulations, Phys. Rev. B 81, 115416 (2010).
  * Fujita  _et al._ [1996] M. Fujita, K. Wakabayashi, K. Nakada, and K. Kusakabe, Peculiar localized state at zigzag graphite edge, J. Phys. Soc. Japan 65, 1920 (1996).
  * Son  _et al._ [2006] Y.-W. Son, M. L. Cohen, and S. G. Louie, Energy Gaps in Graphene Nanoribbons, Phys. Rev. Lett. 97, 216803 (2006).
  * Fernández-Rossier [2008] J. Fernández-Rossier, Prediction of hidden multiferroic order in graphene zigzag ribbons, Phys. Rev. B 77, 075430 (2008).
  * Pisani  _et al._ [2007] L. Pisani, J. A. Chan, B. Montanari, and N. M. Harrison, Electronic structure and magnetic properties of graphitic ribbons, Phys. Rev. B 75, 064418 (2007).
  * Tao  _et al._ [2011] C. Tao, L. Jiao, O. V. Yazyev, Y.-C. Chen, J. Feng, X. Zhang, R. B. Capaz, J. M. Tour, A. Zettl, S. G. Louie, H. Dai, and M. F. Crommie, Spatially resolving edge states of chiral graphene nanoribbons,, Nat. Phys. 7, 616 (2011).
  * Carvalho  _et al._ [2014] A. R. Carvalho, J. H. Warnes, and C. H. Lewenkopf, Edge magnetization and local density of states in chiral graphene nanoribbons, Phys. Rev. B 89, 245444 (2014).
  * Yazyev and Helm [2007] O. V. Yazyev and L. Helm, Defect-induced magnetism in graphene, Phys. Rev. B 75, 125408 (2007).
  * Miranda  _et al._ [2016] V. G. Miranda, L. G. G. V. Dias da Silva, and C. H. Lewenkopf, Coulomb charging energy of vacancy-induced states in graphene, Phys. Rev. B 94, 075114 (2016).
  * Chen  _et al._ [2011] J.-H. Chen, L. Li, W. G. Cullen, E. D. Williams, and M. S. Fuhrer, Tunable Kondo effect in graphene with defects, Nat. Phys. 7, 535 (2011).
  * Miranda  _et al._ [2014] V. G. Miranda, L. G. G. V. Dias da Silva, and C. H. Lewenkopf, Disorder-mediated Kondo effect in graphene, Phys. Rev. B 90, 201101 (2014).
  * Jiang  _et al._ [2018] Y. Jiang, P.-W. Lo, D. May, G. Li, G.-Y. Guo, F. B. Anders, T. Tanigushi, K. Watanabe, J. Mao, and E. Y. Andrei, Inducing kondo screening of vacancy magnetic moments in graphene with gating and local curvature, Nat. Commun. 9, 2349 (2018).
  * Vidarte and Lewenkopf [2022] K. J. U. Vidarte and C. Lewenkopf, High-energy Landau levels in graphene beyond nearest-neighbor hopping processes: Corrections to the effective Dirac Hamiltonian, Phys. Rev. B 106, 155414 (2022).
  * Haydock  _et al._ [1972] R. Haydock, V. Heine, and M. J. Kelly, Electronic structure based on the local atomic environment for tight-binding bands, J. Phys. C: Solid State Phys. 5, 2845 (1972).
  * Haydock  _et al._ [1975] R. Haydock, V. Heine, and M. J. Kelly, Electronic structure based on the local atomic environment for tight-binding bands: II, J. Phys. C: Solid State Phys. 8, 2591 (1975).
  * Haydock [1980] R. Haydock, Comput. Phys. Commun. 20, 11 (1980).
  * 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).
  * Rozhkov  _et al._ [2016] A. Rozhkov, A. Sboychakov, A. Rakhmanov, and F. Nori, Electronic properties of graphene-based bilayer systems, Phys. Rep. 648, 1 (2016).
  * Katsnelson [2020] M. Katsnelson,  _The Physics of Graphene_, Nanoscience and technology (Cambridge University Press, 2020).
  * Moon and Koshino [2013] P. Moon and M. Koshino, Optical absorption in twisted bilayer graphene, Phys. Rev. B 87, 205404 (2013).
  * Mele [2010] E. J. Mele, Commensuration and interlayer coherence in twisted bilayer graphene, Phys. Rev. B 81, 161405 (2010).
  * Mele [2012] E. J. Mele, Interlayer coupling in rotationally faulted multilayer graphenes, Journal of Physics D: Applied Physics 45, 154004 (2012).
  * Yao  _et al._ [2018] W. Yao, E. Wang, C. Bao, Y. Zhang, K. Zhang, K. Bao, C. K. Chan, C. Chen, J. Avila, M. C. Asensio, J. Zhu, and S. Zhou, Quasicrystalline 30∘ twisted bilayer graphene as an incommensurate superlattice with strong interlayer coupling, Proc. Natl. Acad. Sci. U.S.A. 115, 6928 (2018).
  * Ahn  _et al._ [2018] S. J. Ahn, P. Moon, T.-H. Kim, H.-W. Kim, H.-C. Shin, E. H. Kim, H. W. Cha, S.-J. Kahng, P. Kim, M. Koshino, Y.-W. Son, C.-W. Yang, and J. R. Ahn, Dirac electrons in a dodecagonal graphene quasicrystal, Science 361, 782 (2018).
  * Auerbach [1998] A. Auerbach,  _Interacting Electrons and Quantum Magnetism_, Graduate Texts in Contemporary Physics (Springer New York, 1998).
  * Uryu [2004] S. Uryu, Electronic states and quantum transport in double-wall carbon nanotubes, Phys. Rev. B 69, 075402 (2004).
  * Hancock  _et al._ [2010] Y. Hancock, A. Uppstu, K. Saloriutta, A. Harju, and M. J. Puska, Generalized tight-binding transport model for graphene nanoribbon-based systems, Phys. Rev. B 81, 245402 (2010).
  * Schüler  _et al._ [2013] M. Schüler, M. Rösner, T. O. Wehling, A. I. Lichtenstein, and M. I. Katsnelson, Optimal hubbard models for materials with nonlocal coulomb interactions: Graphene, silicene, and benzene, Phys. Rev. Lett. 111, 036601 (2013).
  * 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).
  * 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).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Magic-angle twisted bilayer graphene as a topological heavy fermion problem, Phys. Rev. Lett. 129, 047601 (2022).
  * Sharpe  _et al._ [2021] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Evidence of orbital ferromagnetism in twisted bilayer graphene aligned to hexagonal boron nitride, Nano Lett. 21, 4299 (2021).