# Unconventional Flat Chern Bands and 2$e$ Charges in Skyrmionic Moiré
Superlattices

Yifei Guan  Institute of Physics, Ecole Polytechnique Fédérale de Lausanne
(EPFL), CH-1015 Lausanne, Switzerland  Oleg V. Yazyev  Institute of Physics,
Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Alexander Kruchkov  Institute of Physics, Ecole Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Department of Physics, Harvard
University, Cambridge, MA 02138, USA  Branco Weiss Society in Science, ETH
Zürich, CH-8092 Zürich, Switzerland

###### Abstract

The interplay of topological characteristics in real space and reciprocal
space can lead to the emergence of unconventional topological phases. In this
Letter, we implement a novel mechanism for generating higher-Chern flat bands
on the basis of twisted bilayer graphene (TBG) coupled to topological magnetic
structures in the form of the skyrmion lattice. In particular, we discover a
scenario for generating $|C|=2$ dispersionless electronic bands when the
skyrmion periodicity and the moiré periodicity are matched. Following the
Wilczek argument, the statistics of the charge-carrying excitations in this
case is bosonic, characterized by electronic charge $Q=2e$, that is even in
units of electron charge $e$. The required skyrmion coupling strength
triggering the topological phase transition is realistic, with its threshold
estimated as low as 4 meV. The Hofstadter butterfly spectrum of this phase is
different resulting in an unexpected quantum Hall conductance sequence
$\pm\frac{2e^{2}}{h},\ \pm\frac{4e^{2}}{h},...$ for TBG with skyrmion order.

The interplay between the real- and momentum-space topologies is a new
direction in exploring interacting topological phases of matter. Generically,
the band topology was introduced in condensed matter physics through the
quantum-Hall-like arguments,Thouless  _et al._ (1982); Haldane (1988) which
are typically insensitive to the real-space defects. Nevertheless, it later
became clear in the context of quantum Hall ferromagnetism that the
coexistence of real-space topology (skyrmions) and momentum-space topology
(Chern numbers) is possible for the flat bands.Sondhi  _et al._ (1993); Fertig
_et al._ (1994) More broadly, the question can be formulated in the form: Can
the real space topology influence the change in the momentum space topology?
The answer is positive;Lux  _et al._ (2021) and we provide a concrete example
of a realistic system.

The skyrmions have been called upon to explain possible pairing mechanisms in
twisted bilayer graphene (TBG) and similar heterostructures, Khalaf  _et al._
(2021); Chatterjee  _et al._ (2020) see also Refs. Zhang  _et al._ , 2019a;
Abanov and Wiegmann, 2001; Grover and Senthil, 2008. In particular, the
skyrmion-based theories are able to address not only the problem of possible
mechanism for unconventional superconductivity in TBG,Cao  _et al._ (2018) but
also provide reasonable explanationChatterjee  _et al._ (2020) for
nonmonotonic magnetoresistance observed in experiments.Sharpe  _et al._
(2019); Serlin  _et al._ (2020) Thus, the coexistence of the skyrmion order
and the moiré potentials, to one or another extent, seems feasible. Since this
is a consistent theoretical mechanism, we ask the following question: What is
the effect of the skyrmion order on the TBG flat bands themselves?

The topological flat bands can induce a non-collinear magnetic order, such as
skyrmions.Fertig  _et al._ (1997) In general, however, the skyrmion
excitations are neither bosonic, nor fermionic, and reshape themselves
continuously between the two opposite quantum statistics. Wilczek and Zee
(1983) It can be shown that excess electric charge density $\rho(\mathbf{r})$,
polarized by the presence of a skyrmion $\mathbf{S}(\mathbf{r})$ on a flat
Chern band is given by Sondhi  _et al._ (1993); Abanov and Wiegmann (2001);
Hsu and Chakravarty (2013)

$\displaystyle\rho(\mathbf{r})=Ce\,\frac{\mathbf{S}\cdot(\partial_{x}\mathbf{S}\times\partial_{y}\mathbf{S})}{4\pi},$
(1)

where $C$ is the Chern number of the underlying electronic band. Upon
integration over the unit cell one obtains

$\displaystyle Q=C\mathcal{W}e,$ (2)

where $\mathcal{W}$ is the skyrmion winding number (an integer). Thus,
skyrmions polarize a discrete electric charge when the underlying flat band is
topological.

The collective effect of the real-space and momentum-space topologies in
mechanism (2) plays an important role. Indeed, the effective magnetic field,
assigned to the skyrmion
$B(\mathbf{r})\propto\mathbf{S}\cdot(\partial_{x}\mathbf{S}\times\partial_{y}\mathbf{S})$
in (1) is essentially the Berry field.Nagaosa and Tokura (2013) In momentum
space, the gluing mechanism (2) can be rewritten as

$\displaystyle\rho(\mathbf{q})=e\sum_{\mathbf{k}}c^{{\dagger}}_{\mathbf{k}+\mathbf{q}}c_{\mathbf{k}}=\mathcal{W}\frac{F_{xy}}{2\pi},$
(3)

where $F_{xy}$ is the Berry curvature associated with the Chern band,
$C=\frac{1}{2\pi}\int_{\text{BZ}}F_{xy}dk_{x}dk_{y}$. Thus, according to Eq.
(3), it is possible to have an unconventional electronic pairing mechanism
induced by the skyrmion order when electrons are moving in the Chern band.
This has been recently pointed out in the models of skyrmionic
superconductivity in TBG, Chatterjee  _et al._ (2020); Khalaf  _et al._ (2021)
where unconventional pairing mechanism is resulting from the skyrmions in the
pseudospin space built upon reshuffling intrinsic $C=+1$ and $C=-1$
sectors.Khalaf  _et al._ (2021) However, a different situation is possible if
the moiré system develops a flat band of a higher Chern number. Namely, in the
$C=2N$ case, the electric charge polarized by the skyrmion is

$\displaystyle Q=2N\mathcal{W}e,$ (4)

that is bosonic. In other words, in case of $|C|=2$ electron pairs are
preformed independently of mechanism developed in Ref.Khalaf  _et al._ (2021).
In this Letter, we show that the flat bands in TBG acquire even Chern numbers
$|C|=2$ when electrons are firmly coupled to the underlying skyrmion lattice.

Figure 1: Unconventional Chern flat bands of $|C|$$=$$2$ in twisted bilayer
graphene matched with the skyrmion lattice. (a) Band structure and density of
states (DOS) of magic-angle TBG without skyrmion order. The associated Chern
number is zero. (b) The inclusion skyrmion order ($m_{0}=4$ meV, $\mu=8$ meV)
results in flat Chern bands with $|C|=2$.

Underlying skyrmion lattice.—The calculations below are independent of a
particular origin of the skyrmion lattice (SkL), and can be applied to both
(i) the case of proximity to a substrate hosting skyrmion lattice, or (ii) the
case of spontaneous formation of skyrmion lattice in the moiré system. In the
former scenario (i), the periodicity of SkL $\lambda_{\text{SkL}}=J/2D$
dictated by the interplay of Heisenberg ($J$) and Dzialoshinsky-Moriya ($D$)
interactions can be engineeredTokura and Kanazawa (2020) to match the
periodicity of the moiré superlattice of TBG ($a_{\text{M}}=12.9$ nm at the
magic angle $\theta=1.08^{\circ}$). The later scenario (ii), represented by
the SkL forming as a result of instability associated with the flat band, is
supported by recent theoretical calculations.Wu and Sarma (2020); Bömerich
_et al._ (2020) At $\nu=3/4$ filling of the conduction band, TBG has been
reported to exhibit giant spontaneous magnetism, resulting from internal Berry
fields of the band structure, that is orbital in nature.Sharpe  _et al._
(2019); Serlin  _et al._ (2020) It is interesting that the TBG flat band can
be viewed as a Landau level in non-homogeneous magnetic field,Tarnopolsky  _et
al._ (2019); Popov and Milekhin (2021); San-Jose  _et al._ (2012) and this
real-space Berry field has strong inhomogeneities of the same periodicity as
the patern moiré potential.Ledwith  _et al._ (2020) Thus, the resulting SkL
instability, if allowed by the mechanism of Ref. Wu and Sarma, 2020, will
naturally inherit the same periodicity as the pristine TBG moiré superlattice.
In what follows below, we consider the skyrmion lattice matched with the
periodicity of TBG, but similar principles can be applied to other moiré
heterostructures.

Figure 2: (a)Phase diagram of flat Chern band phases in twisted bilayer
graphene with the skyrmion order. We see that at the low skyrmion field the
topologically–trivial phases are energetically preferred. However, at higher
skyrmion fields, the phase diagram shows a remarkable new phase of $C$$=$$2$
flat composite bands. (b) The calculated electric polarization of the
different TBG phases with skyrmion lattice, see Eq. (6). The numerical
parameters used: $\mu=8$ meV, $m_{0}=4$ mev ($C=2$ cylinder) and $\mu=2$ meV,
$m_{0}=4$ mev ($C=0$ cylinder). (c) Sketched phase volume depicting the
polarized charge 2$e$ (in units of Fig.2a).

Band structure and unconventional Chern flat bands.—We start by addressing the
effect of skyrmion order on the band structure of TBG. We demonstrate that the
skyrmion lattice, when being commensurate with the moiré periodicity
$a_{\text{M}}$ at the first magic angle, changes the electronic response of
the system beyond recognition. We define the underlying SkL as Mühlbauer  _et
al._ (2009); Nagaosa and Tokura (2013)

$\displaystyle\mathbf{m}(\mathbf{r})=\mathbf{m}_{0}+\mu\sum^{3}_{j=1}(\mathbf{S}_{\mathbf{b}_{j}}e^{i\mathbf{b}_{j}\mathbf{r}}+\mathbf{S}_{-\mathbf{b}_{j}}e^{-i\mathbf{b}_{j}\mathbf{r}}),$
(5)

where $m_{0}$ is the uniform magnetization component, $\mu$ the SkL strength
and reciprocal vectors $\mathbf{b}_{1}=\mathbf{q}_{2}-\mathbf{q}_{1}$,
$\mathbf{b}_{2}=\mathbf{q}_{3}-\mathbf{q}_{1}$,
$\mathbf{b}_{3}=\mathbf{q}_{3}-\mathbf{q}_{2}$ are the same for the SkL and
the moiré superlattice (here $\mathbf{q}_{1}=k_{\theta}(0,-1)$,
$\mathbf{q}_{2,3}=k_{\theta}(\pm\sqrt{3}/2,1/2)$, with moiré wave number
$k_{\theta}=2k_{D}\sin\frac{\theta}{2}$, $k_{D}$ being Dirac momentum of the
monolayer graphene and $\theta$ being the twist angle). For the considered
ranges of parameters $m_{0},\mu=0..10$ meV, the effective local magnetic
fields are sufficient to modify the band structure (Fig. 1) and promote
unconventional Chern phases (Fig. 2).

We focus our attention on the effect of the underlying skyrmion order (5) on
the electronic band structure of TBG at the magic angle $\theta=1.08^{\circ}$.
We have first confirmed this effect within the continuum model (see
Supplemental Materials), and then performed the tight-binding model
calculations for spinful electrons on a magic-angle TBG configuration that
takes into account lattice relaxation effects. The effect of skyrmion spin
texture is introduced via the on-site terms
$H_{s}(\mathbf{r}_{i})=\boldsymbol{\sigma}\cdot\mathbf{m}(\mathbf{r}_{i})$ for
each of the orbitals. Further details of our calculations are presented in the
Methods section.

Figure 1a compares the electronic band structure of TBG without the skyrmion
order with that (Fig. 1b) for a representative case characterized by $m_{0}=4$
meV and $\mu=8$ meV. Importantly, these parameters result in the $Q=2e$ phase
as shown in phase diagram in Fig. 2a. In absence of SkL ($m_{0}=0$ meV,
$\mu=0$ meV) magic-angle TBG has 8 flat bands near the Fermi level isolated
from remote bands (Fig.1a). At the $K$ point of mini Brillouin zone, the
flatband manifold has degenerate Dirac-like band crossing with renormalized
Fermi velocity, while the degeneracy at the $\Gamma$ point due to the spin and
valley degrees of freedom. Introducing the uniform magnetization ($m_{0}\neq
0$) lifts the spin degeneracy splitting the flat band manifold into 2 groups,
while leaving intact the degeneracy at the Dirac points (see the Supplemental
Material). In the presence of skyrmion order ($\mu\neq 0$), the spin-orbit
terms couple the spin with the bands opening gaps at points $K$ and
$K^{\prime}$ and contributing to the hybridization of flat bands at the
$\Gamma$ point (Fig. 1b). In this case, the 8 flat bands decompose into 4
doubly degenerate sub-bands, which remain significantly flat. Employing the
band flatness criterion of Bistritzer and MacDonaldBistritzer and MacDonald
(2011) (renormalized Fermi velocity $v_{F}$ the at the $K$ point), we observe
that these bands become very flat in the $K$, $K^{\prime}$ valleys
($v_{F}=3.3\times 10^{4}$ m/s without SkL vs. $v_{F}=2\times 10^{2}$ m/s with
SkL). For comparison, the overall bandwidths in TBG with SkL are significantly
lower than in the pristine TBG case (5.3 mev, 13 meV in pristine TBG and 2.7
meV, 5 meV, 8.5 meV, 9.3 meV in TBG with SkL, given in ascending bandwidth
order). Importantly, the density of states (DOS) remains significant upon
adding SkL term (Fig. 1), thus enabling strong correlations. To quantify the
the band flatness further, we investigate the band flatness criterion in terms
of wave functions.Kruchkov (2021) We thus confirm that upon adding the SkL
term, the flat bands remain well-defined with clear real-space localization
resulting into significant electronic density of states (Fig. 1b), with two of
the resulting flat Chern bands having the high Chern number $|C|=2$. Thus,
non-collinear spin texture of skyrmion lattice coupled to charge-carrier in
twisted bilayer graphene modifies its original band structure, but preserves
the flatness in the atomistic model with lattice relaxations We now proceed to
the effect on band topology.

Figure 3: Hofstadter and Wannier diagrams of TBG without and with the skyrmion
order Hofstadter (a) and Wannier (b) diagrams of twisted bilayer graphene
without the skyrmion order. Hofstadter (c) and Wannier (d) diagrams of twisted
bilayer graphene with the skyrmion order. We clearly see a different pattern
of gaps and Chern numbers, attributed to the magnetic-field responds of a flat
Chern band (see main text). All calculations are performed with the tight-
binding model accounting for atomic relaxations of TBG. Here $n_{S}=8$ (total
bands) and $\Phi_{0}=h/e$.

Skyrmion-induced band topology and the phase diagram.— The pristine TBG has
‘band structure with hidden Chern numbers $C=\pm 1$ per valley that sum up to
zero Chern number in total. Thus, in absence of skyrmion lattice or explicit
symmetry breaking, TBG is not expected to demonstrate Chern indices of the
composite flat bands (i.e. $m_{0}=\mu=0$ meV in phase diagram shown in Fig.2a
is topologically trivial). In what follows below, we discuss the half-filling
of the conduction and valence bands ($\nu=3/4$ and $\nu=1/4$ fillings of the 8
moiré flat bands, respectively), motivated by the experimentally relevant
filling $\nu=3/4$ (see e.g. Ref.Sharpe  _et al._ (2019)). The formalism of
hybrid Wannier functions (HWFs)Soluyanov and Vanderbilt (2011) is employed to
calculate the positions of the Wannier charge centers (WCC). This formalism
gives us advantage of both addressing the Chern numbers as the winding of WCC
on the cylinder (see Fig. 2b), and linking the nontrivial winding to electric
polarization: The sum of all Wannier charge centers $\langle
W^{h}_{n}|\mathbf{r}|W^{h}_{n}\rangle$ is gauge-invariant (mod lattice
constant), and relates to the electric polarization of the systemKing-Smith
and Vanderbilt (1993); Resta (1993); Coh and Vanderbilt (2009)

$\displaystyle P_{e}=e\sum_{n}\langle
W^{(h)}_{n}|\mathbf{r}|W^{(h)}_{n}\rangle.$ (6)

Upon including the skyrmion order the $C=2$ phase exhibits $\Delta P_{e}=2ea$,
witnessed through the polarized charge $2e$. For definiteness, in Fig. 2a and
the text below, the phase diagram traces the total Chern number at $\nu=3/4$,
that is of the lowest three composite flat bands, which is the same in
absolute value as the Chern number of the upper composite flat band, in the
range from 0 to 10 meV for both the uniform component $m_{0}$ and the skyrmion
strength $\mu$. This range is dictated by the typical bandwidth scale in the
system. Larger values of $m_{0}$ (over 20 meV) drives the flat bands to
overlap with remote bands and destroys the separation condition of the flat
bands, while $\mu$ does not induce the overlap explicitly.

The phase diagram in Figure 2a can be divided into three parts: two
topologically-trivial phases (gapped $C=0$ phase and gapless phase denoted
”gapless” for which the Chern numbers cannot defined). Below the critical
value $\mu\approx 4$ meV, the system remains in one of the trivial phases.
Upon increasing the skyrmion strength above $\mu=4$ meV, the new topological
phase with $C=2$ emerges. This phase, however is possible if only both $\mu$
and $m_{0}$ are finite: As shown in Fig. 2a, the WCC winding is zero when
$m_{0}=0$. One can find a tiny band gap near the $\Gamma$ point at $m_{0}=0$,
in which the Chern number is zero. However, a weak ferromagnetic component
$m_{0}\ll\mu$ is sufficient to drive the system into the $C=2$ phase, provided
the skyrmion strength with $\mu>4$ meV. On the other hand, the ferromagnetic
component $m_{0}$ does not contribute to the finite Chern number itself, since
in the region $\mu<4$ meV the total Chern number is always zero (or undefined
as a consequence of gaplessness). We highlight that the skyrmion strength of 4
meV is a realistic value, comparably smaller than the bandwidth itself, and a
small uniform magnetization turns the system into a higher-Chern topological
phase.

Hofstadter spectra.— Another qualitative difference of the system with
skyrmionic order can be traced in the response to the magnetic field. In this
connection, the established apparatus involves the Hofstadter butterfly
spectrum, which depicts the evolution of Landau levels in the system versus
magnetic flux.Hofstadter (1976) The Hofstadter butterfly spectra, together
with the Wannier diagrams,Wannier (1978) provide a comprehensive and sensitive
map of the Chern numbers, and is directly linked to the observable quantum
Hall conductance sequences.Dean  _et al._ (2013); Hunt  _et al._ (2013);
Ponomarenko  _et al._ (2013) In Figure 3, we show the difference between the
magnetic field response of magic-angle TBG with and without the skyrmion
order. For the skyrmion-free TBG, the Landau levels at small magnetic fluxes
are clearly recognizable, and their field-dependence in the limit of low field
is consistent with the behavior of Dirac electrons scaling as $\sqrt{B}$ (Fig.
3a). In this limit, the Hofstadter spectrum is energetically bound to the
bandwidth of the original flat band, which is expected for the trivial
phase.Herzog-Arbeitman  _et al._ (2020) Upon increasing the magnetic flux, a
moderate broadening of the Landau levels is observed, while sufficiently large
gaps separate low-index Landau levels.

The situation changes dramatically upon including the skyrmion order (Fig. 3c;
for $m_{0}=4$ meV and $\mu=8$ meV). First of all, the low-field behavior is no
longer recognizable as a simple $\sqrt{B}$ behavior since the Dirac physics is
no longer relevant in this case. Instead, the behavior is qualitatively
consistent with universal response of a Chern band to magnetic fields.Herzog-
Arbeitman  _et al._ (2020) A different pattern of pronounced gaps emerges upon
increasing the flux. The broadening and hybridization between sub-bands is
promoted, leading to a more complex butterfly structure. We observe that the
Hofstadter spectrum of the flat band is not bounded by the bandwidth of the
original flat band, and can hybridize with other bands, which in turn
trivializes the band topology leading to $C=0$ gaps. Importantly, the
particle-hole asymmetry is pronounced, which results in the particle-hole
asymmetric Wannier diagrams and asymmetric Hall response. Since the TBG with
SkL host a Chern band characterized by $C=2$, we observe a different pattern
of in-gap Chern numbers, as indicated in Fig. 3c. In order to draw the
quantitative difference between these two cases, we proceed to the analysis of
the Wannier diagrams.

Table 1: Summary of quantum Hall conductance sequences at neutrality.  TBG ($C_{3}$ symmetric) |  $\pm 4\frac{e^{2}}{h}$, $\pm 12\frac{e^{2}}{h}$, …
---|---
TBG+SkL ($C$$=$$0$ phase) |  $\pm 4\frac{e^{2}}{h}$ , $\pm 8\frac{e^{2}}{h}$ , …
TBG+SkL ($C$$=$$2$ phase) |  $\pm 2\frac{e^{2}}{h}$ , $\pm 4\frac{e^{2}}{h}$, …
TBG experiments |  $\pm 2\frac{e^{2}}{h}$ , $\pm 4\frac{e^{2}}{h}$, … + further sample-dependent40

Hall conductance sequences.—The Wannier diagram Wannier (1978) is obtained
from the Hofstadter butterfly spectrum by presenting the statistical weight of
the states below a given gap versus the magnetic flux. Each line in the
Wannier diagram depicts a gap in the Hofstadter butterfly spectrum, while its
slope characterizes the in-gap Chern numbers ($dn/dB=Ce/h$). In Figs. 3b,d we
plot the Wannier diagrams that correspond to the respective Hofstadter spectra
shown in Figs. 3a,c. To facilitate the analysis, we overlay the lines that
trace the minimal intensities in the Wannier diagrams (Fig. 3c,3d). We find
that for pristine TBG with the preserved $C_{3}$ symmetry, the relevant Hall
conductance sequence is $\pm 4\frac{e^{2}}{h}$, $\pm 12\frac{e^{2}}{h}$, …,
that is in increments of 8 conductance quanta.Zhang  _et al._ (2019b) On
contrary, the presence of skyrmion order changes the Hall sequence
dramatically. For the TBG with SkL in the $C=2$ phase, the Hall conductance
”fingerprint” (Fig. 3) at neutrality reads

$\displaystyle G_{\text{TBG}}^{\text{SkL}}=\pm\frac{2e^{2}}{h},\
\pm\frac{4e^{2}}{h},...\ \ \ $ (7)

One might be tempted to interpret these results as a trivial lifting of the
spin degeneracy by the uniform component $m_{0}$ of the skyrmion phase.
However, we observe that by changing the skyrmion parameters to $m_{0}=4$ meV
and $\mu=2$ meV, for which the system is in the $C=0$ phase (see Fig. 2a), a
different Hall sequence $\pm 4\frac{e^{2}}{h},-8\frac{e^{2}}{h}$ is obtained,
even though the system still has the same uniform magnetization $m_{0}$.
Instead, we believe that this effect is connected to the fundamentally
different Hosftadter spectrum for the trivial and Chern bands, and hence,
Wannier diagrams. The obtained results are summarized in Table I. Wannier
diagrams can be directly probed in experiment by performing Shubnikov–de Haas
measurements on high-quality magic-angle TBG samples.Lu  _et al._ (2019) We
remark that the experimentally observed QH sequences are strongly sample-
dependent,111The question of the Landau fan (LF) in TBG is a hard one. The
original studiesCao  _et al._ (2018) reported the Hall sequence, in units
$e^{2}/h$: $\pm 4$, $\pm 8$, $\pm 12$ at neutrality, instead of expected $\pm
4$, $\pm 12$…Moon and Koshino (2012), thus puzzling theorists. The consequent
experimentsLu  _et al._ (2019) on more homogeneous samples brought even
further surprise with the QH sequence $\pm 2$, $\pm 4$, $\pm 8$. Moreover, LF
with $\pm 2$, $\pm 4$ was reported in Ref.Sharpe  _et al._ (2019); LF with
$\pm 2$, $\pm 4$, $\pm 6$ was reported in TBG stabilized by WSe2;Arora  _et
al._ (2020) and LF with $-2$, $-4$ in the sample of Ref.Pierce  _et al._
(2021). For the high-field behavior of Landau fans, we refer readers to the
recent observation of correlated Chern insulators in TBG e.g. in Refs. Pierce
_et al._ (2021); Xie  _et al._ (2021) however the two leading contributions to
the low-field Landau fans can be considered as $\pm 2\frac{e^{2}}{h}$, $\pm
4\frac{e^{2}}{h}$. It is certainly interesting that TBG with SkL in the $C=2$
phase (but not the $C=0$ phase) gives a reasonable Hall sequence, without
demanding $C_{3z}$ or $C_{2z}$ breaking.Zhang  _et al._ (2019b)

Conclusion.—In this Letter, we considered the influence of the real-space
skyrmionic order on the dispersion and topology of flat bands. For the
particular example of magic-angle twisted bilayer graphene, we find that the
commensurate skyrmionic order itself redefines both the band structure and
topology of the system. Surprisingly, the real-space skyrmion order influences
the electronic bands topology, giving the rise to an unexpected $C=2$ phase.
To our knowledge, this is the first realistic system which provides a robust
and affirmative answer to the question of the interplay between topologies in
real and reciprocal spaces.Lux  _et al._ (2021) According to the field-
theoretical arguments,Sondhi  _et al._ (1993); Abanov and Wiegmann (2001); Hsu
and Chakravarty (2013) the elementary excitations in this case are $Q=2e$, a
consequence important for the transport phenomena in twisted bilayer graphene.
Experimentally, the reported effect manifests in magnetotransport observables
with the leading quantum Hall sequence $\pm 2\frac{e^{2}}{h}$, $\pm
4\frac{e^{2}}{h}$. This Landau fan fingerprint resonates well with
experimental reports, and arises naturally in our model as a direct
consequence of the emergent skyrmion order.

METHODS

Tight binding model for twisted bilayer graphene.—The electronic structure
calculations are performed with the tight-binding Hamiltonian capturing the
atomistic physics of twisted bilayer graphene in free-electron formalism,

$\displaystyle H=\sum_{ij}t_{ij}c^{{\dagger}}_{i}c_{j},$ (8)

where $c^{{\dagger}},c$ are creation and annihilation operators and the
numerical values of hopping terms $t_{ij}$ are computed within the Slater-
Koster formalism by taking into account hopping of $p_{z}$ orbital on each
atom. The hopping term between two $p_{z}$ orbitals is expressed as a sum of
$\pi$ and $\sigma$ type overlapping, $t_{ij}$$=$$V_{i}^{\pi}+V_{ij}^{\sigma}$,
conventionally approximated in the formNam and Koshino (2017)

$\displaystyle
t_{ij}=V_{0}^{\pi}\exp(-\frac{r-a_{0}}{r_{0}})\sin^{2}\theta+V_{0}^{\sigma}\exp(-\frac{r-d_{0}}{r_{0}})\cos^{2}\theta,$
(9)

where $r$$=$$|\mathbf{r}_{ij}|$ and $\theta$ is the angle between
$\mathbf{r}_{ij}$ and the $\hat{e}_{z}$,
$\theta$$=$$\angle\\{\mathbf{r}_{ij},(0,0,1)\\}$. The numerical parameters
have been adapted from Ref.Nam and Koshino (2017), hence
$V_{0}^{\pi}$$=$$-2.7$ eV is the hopping integral between nearest-neighbor
(NN) atoms in the same layer; $V_{0}^{\sigma}$=0.48 eV is the overlap between
the aligned atoms in different layers; here $a_{0}$$=$$0.142$ nm and
$d_{0}$=0.334 nm correspond to the nearest neighbor distance and the
interlayer distance; $r_{0}$=0.184$a$ (we use $a$$=$$\sqrt{3}a_{0}$ as the
lattice constant of graphene) is the hopping strength decay length. The cutoff
real-space is chosen at 9.84 Å; in the outer region the coupling is
negligible. To the main approximation, the hopping amplitudes are barely
affected by the real-space magnetic order. Hence on the tight-binding level
the effect of local magnetization enters the Hamiltonian as on-site spin-orbit
coupling terms on each atom. In the basis
$(\psi_{i}^{\uparrow},\psi_{i}^{\downarrow})$ the spin-orbit coupling reads as
$H_{\text{SOC}}$$=$$\boldsymbol{\sigma}\cdot\mathbf{m}(\mathbf{r}_{i})$, where
$\boldsymbol{\sigma}$$=$$(\sigma_{x},\sigma_{y},\sigma_{z})$ are the
conventional 2$\times$2 Pauli matrices. At atomic position $i$, the on-site
Hamiltonian term is

$\displaystyle H_{i}^{\text{os}}$ $\displaystyle=\begin{pmatrix}1&0\\\
0&1\end{pmatrix}V_{i}+\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\mu_{x}+\begin{pmatrix}0&\text{-}i\\\
i&0\end{pmatrix}\mu_{y}+\begin{pmatrix}1&0\\\
0&\text{-}1\end{pmatrix}[\mu_{z}+m_{0}],$ (10)

where we distinguish between the uniform magnetization
($m_{0}$$=$$\langle\mathbf{m}(\mathbf{r})\rangle$) and the skyrmion order
itself $\boldsymbol{\mu}$$=$$\mathbf{m}(\mathbf{r})$$-$$\mathbf{m}_{0}$. The
overall real-space magnetic phase is approximated by

$\displaystyle\mathbf{m}(\mathbf{r})=\mathbf{m}_{0}+\mu\sum_{j=1}(\mathbf{S}_{\mathbf{q}_{j}}e^{i\mathbf{q}_{j}\mathbf{r}}+\mathbf{S}_{-\mathbf{q}_{j}}e^{-i\mathbf{q}_{j}\mathbf{r}}),$
(11)

where for the skyrmion lattice it is sufficient to cut the Fourier terms by
the first triade $\mathbf{q}_{1}+\mathbf{q}_{2}+\mathbf{q}_{3}$$=$$0$, with
$|\mathbf{q}_{i}|$ commensurate with moiré periodicity. Adding such on-site
terms on each orbital, the full spinful Hamiltonian with the underlying
skyrmion order reads

$\displaystyle H$ $\displaystyle=H_{\text{hop}}+H_{\text{os}}=\sum_{i\neq
j}t_{ij}\sigma_{0}+\sum_{i}\sigma_{0}V_{i}$ (12)
$\displaystyle+\sum_{i}\sigma_{x}\mu_{x}(r_{i})+\sigma_{y}\mu_{y}(r_{i})+\sigma_{z}[\mu_{z}(r_{i})+m_{0}].$
(13)

Since the atomic relaxation has an effect on the interlayer coupling in AA and
AB stacking region, and has dramatic effects on band flatness in TBG, we build
the tight-binding hamiltonian based on the relaxed atomic structures. The
atomic structure relaxation is done with classical potential on the LAMMPS
package.Plimpton (1995) The classical force field contains the short range
terms for covalent bonds and long-range terms which describes the Van der
Waals interaction. The implemented long-range potentials are obtained from the
parameterized version of the Kolmogorov-Crespi potential and confirmed to
match the DFT results in Ref. Gargiulo and Yazyev (2017). The short range
interaction of long-range carbon bond order potential (LCBOP) was adopted from
Los and Fasolino (2003) without modification.

Hybrid Wannier functions formalism—The Chern numbers of the moiré flat bands
are calculated through the formalism of hybrid Wannier function methods
introduced by Soluyanov and Vanderbilt,Soluyanov and Vanderbilt (2011) see
also Refs.Taherinejad  _et al._ (2014); Wu  _et al._ (2018) The Wannier
functions $W_{n}(\mathbf{R})$ are conventionally defined as localized orbitals
which are build through integral representation through the Bloch wave
functions $u_{k}$,Wannier (1962)

$\displaystyle
W_{n}(\boldsymbol{R})=\frac{1}{(2\pi)^{2}}\int_{\text{BZ}}e^{i\boldsymbol{k}\cdot(\boldsymbol{r}-\boldsymbol{R})}|u_{n\boldsymbol{k}}\rangle.$
(14)

On contrary, the hybrid Wannier functions (HWFs) are wannierized only along
one of the two spatial dimensions (e.g. $x$), thus remaining a Bloch-periodic
function in momentum along the other dimension (e.g. $k_{y}$). HWFs had proved
to be useful to track the topological properties, in particular Chern numbers
of flat electronic bands.Gresch  _et al._ (2017) The underlying idea is to
track the evolution of the hybrid Wannier functions $|W^{(h)}_{n}\rangle$
along a $k_{y}$-path, and the Chern number is given by the total winding of
the Wannier centers. This approach is based on the observation that every
nontrivial Chern number presents an obstruction for constructing maximally-
localized Wannier functions.Marzari  _et al._ (2012) The numerical advantage
we implement is: instead of finding the full 2D Wannier functions, we
calculate the Wannier charge centers(WCCs) in one of the dimensions (e.g.
$x$), and track it along the momentum on the other dimension($k_{y}$). Here we
give an example of the calculation of hybrid Wannier functions. For the 2D
electronic system with its eigenstates in BZ being
$|u_{\mathbf{k}}(k_{x},k_{y})\rangle$ in the conventional Bloch band
framework, the Wannier charge centers (WCC) are defined through the phase
accumulation,

$\displaystyle|u_{\mathbf{k}}(k_{x}+2\pi,k_{y})\rangle=e^{i\phi_{x}}|u_{\mathbf{k}}(k_{x},k_{y})\rangle,\
\ \ \phi_{x}(k_{y})\equiv\text{WCC}(k_{y}).$ (15)

The numerical calculations are performed on the discretized BZ. The phase
factor $\phi_{x}$ is calculated by evolving from $k_{x}=0$ to $k_{x}=2\pi$,
and tracing the accumulated phase in each increment

$\displaystyle\phi_{x}(k_{y})=\arg\left[\prod_{i=1}^{n-1}|u_{\mathbf{k}}(k_{x}^{i},k_{y})\rangle\,\langle
u_{\mathbf{k}}(k_{x}^{i+1},k_{y})|\right]$ (16)

in order to preserve the smooth gauge. Thus defined, the WCC winding on the
cylinder determines the Chern numbers by counting the number of times the WCC
is crossing the boundaries.Gresch  _et al._ (2017) The winding of the Wannier
center is giving the Chern number of a single or a composite band. In
multiple-band systems, the total Chern number is determined by the addition of
the winding of all the bands. This method becomes particularly handy for
evaluating the Chern numbers of the flat electronic bands.

Effective continuum model.—We consider the twisted bilayer graphene tuned to
the first magic angle, and coupled to the skyrmion magnetic order with the
same modulation vectors as the Moire pattern in TBG. At such small angles and
large supercell scales, the effective physics of both the TBG and SkL can be
described within the effective continuum models (see e.g. Ref.Tarnopolsky  _et
al._ (2019)). We first consider the SkL with uniform (ferromagnetic) component
to be small. The effective Hamiltonian of the TBG coupled to the magnetic
skyrmion Lattice $\mathbf{m}(\mathbf{r})$ can be written as:

$\displaystyle\mathcal{H}=\begin{pmatrix}-iv\boldsymbol{\sigma}_{-\theta/2}\nabla-\mathbf{m}(\mathbf{r})\,\mathbf{s}&T(\mathbf{r})\\\
T^{{\dagger}}(\mathbf{r})&-iv\boldsymbol{\sigma}_{+\theta/2}\nabla-\mathbf{m}(\mathbf{r})\mathbf{s}\end{pmatrix},$
(17)

where $\mathbf{s}$ is electron spin. Here $v$ is the input fermi velocity and
rotated Pauli matrices for each layer defined as
$\boldsymbol{\sigma}_{\theta}$$=$$e^{-i\theta\sigma_{z}/2}(\sigma_{x},\sigma_{y})e^{+i\theta\sigma_{z}/2}$.
The interlaying coupling is given by the 2x2 matrix

$\displaystyle
T(\mathbf{r})=\sum_{j=1}^{3}T_{j}e^{-i\mathbf{q}_{j}\mathbf{r}},$ (18)

where vectors $\mathbf{q}_{1}$, $\mathbf{q}_{2}$, $\mathbf{q}_{3}$ constitute
the three-fold rotation with
$\mathbf{q}_{1}+\mathbf{q}_{2}+\mathbf{q}_{3}$$=$$0$,
$|\mathbf{q}_{j}|$$=$$k_{\theta}$$=$$2k_{D}\sin(\theta/2)$, responsible for
the Moire lattice formation. We here use
$\mathbf{q}_{1}$$=$$k_{\theta}(0,-1)$,
$\mathbf{q}_{2,3}$$=$$k_{\theta}(\pm\sqrt{3}/2,1/2)$. In general the
symmetries of the TBG allow coupling matrices $T_{j}$ of the form,

$\displaystyle
T_{j+1}=w_{\text{AA}}\,\sigma_{0}+w_{\text{AB}}\,(\sigma_{x}\cos
n\phi+\sigma_{y}\sin n\phi),\ \ \ \phi=\pi/3$

We consider a magnetic skyrmion Latice (SkL) matched to the twisted graphene
bilayer in the way that both the supercell periodicities coincide precisely,
and thus the skyrmion lattice is build on the
$\mathbf{Q}_{1}$$=$$\mathbf{b}_{1}$, $\mathbf{Q}_{2}$$=$$-\mathbf{b}_{2}$,
$\mathbf{Q}_{3}$$=$$\mathbf{b}_{3}$ commensurate with the Moire lattice. The
skyrmion lattice order parameter to a good accuracy can be written as

$\displaystyle\mathbf{m}(\mathbf{r})=\mathbf{m}_{0}+\mu\sum_{j=1}(\mathbf{S}_{\mathbf{b}_{j}}e^{i\mathbf{b}_{j}\mathbf{r}}+\mathbf{S}_{-\mathbf{b}_{j}}e^{-i\mathbf{b}_{j}\mathbf{r}}),$
(19)

where we are primarily building the skyrmion lattice on the moire modulation
vectors $\mathbf{b}_{1}$$=$$\mathbf{q}_{2}-\mathbf{q}_{1}$,
$\mathbf{b}_{2}$$=$$\mathbf{q}_{3}-\mathbf{q}_{1}$,
$\mathbf{b}_{3}$$=$$\mathbf{q}_{3}-\mathbf{q}_{2}$, (thus
$\mathbf{Q}_{1}+\mathbf{Q}_{2}+\mathbf{Q}_{3}$$=$$0$ form the SkL
crystalline), plus further momenta terms if needed. Here
$\mathbf{m}_{0}||\mathbf{B}$ is a uniform component proportional to the
strength of external magnetic field. The typical band structures are listed in
Supplementary Material, and qualitatively coincide with the tight binding
calculations. Importatly, the continuum model shows that the bands cannot be
made perfectly flat once the SkL is added.

Hofstadter diagrams: Lanczos algorithm.— For computation of the Hofstadter
diagrams the diagonalization is performed on a magnetic supercell retaining
the periodicity of the boundary conditions. Following the Landau gauge,
$\mathbf{A}$$=$$(-By,0,0)$, each tight-binding calculation on such magnetic
supercells contain diagonalization of large matrices, usually up to dimension
$\sim$ $10^{7}$. Regular diagonalization methods(with complexity $O(N^{3})$
fail quickly in such a large job due to the consumption of time and memory. In
order to perform the required calculation of the Hofstadter and Wannier
diagrams, we make use of the Lanczos algorithm implemented in the WannierTools
package.Wu  _et al._ (2018) The Lanczos algorithm solves the eigenvalue or
density of states problems by iteration, requiring much less memory than the
direct diagonalization. The idea of the algorithm is to iteratively apply the
Hamiltonian matrix to the set of initial vectors, and subsequently transform
the Hamiltonian transformed into a tridiagonal form during the process. To
diagonalize a Hamiltonian $H$, the algorithm starts with an initial vector
$|v_{1}\rangle$ whose norm is 1, each iteration is done as following:

1\. The main diagonal element is given by $\alpha_{i}=\langle
v_{i}|H|v_{i}\rangle$.

2\. Get the second vector
$|v_{i}+1\rangle=H|v_{i}\rangle-\alpha_{i}|v_{i}\rangle$.

3\. Find the sub-leading diagonal $\beta_{i}=||(|Hv_{i}+1\rangle)||$.

4\. Rescale $|v_{i+1}\rangle=|v_{i+1}\rangle/\beta_{i}$.

The operations are further performed iteratively until the convergence is
achieved. By the above iteration, the Hamiltonian is transformed into a
tridiagonal form:

$\displaystyle H_{\text{tri}}=\begin{pmatrix}\alpha_{1}&\beta_{1}&0&0&0\\\
\beta_{1}&\alpha_{2}&\beta_{2}&0&0\\\ 0&\beta_{2}&...&\beta_{n-1}&0\\\
0&0&\beta_{n-1}&\alpha_{n}&\beta_{n}\\\ 0&0&0&\beta_{n}&...\end{pmatrix}.$
(20)

Having the tridiagonalized Hamiltonian, the density of states is thereafter
evaluated through the Green’s functions,

$\displaystyle\rho(E+i\varepsilon)=-\text{Im}\Tr([G(E+i\varepsilon)])=-\text{Im}\Tr[(E+i\varepsilon)I-H]^{-1}.$
(21)

One further applies the continued fraction method to calculate
$G$$=$$[(E+i\varepsilon)I-H]^{-1}$. For the complex frequency
$\Omega$$=$$E+i\varepsilon$, the continued fraction representation is

$\displaystyle G(\Omega)=(\Omega-
H_{\text{tr}})^{-1}=\frac{1}{\Omega-\alpha_{1}-\frac{\beta_{1}^{2}}{\Omega-\alpha_{2}-\frac{\beta_{2}^{2}}{\Omega-...}}}.$
(22)

Further-on, the DOS in the selected energy range is obtained by evaluating the
Green’s function iteratively.

Acknowledgments. The authors thank Bertrand Halperin and Henrik Rønnow for
useful discussion. A.K. thanks Ashvin Vishwanath for multiple discussions on
the topic. The project was supported by the Branco Weiss Society in Science,
ETH Zurich, through the research grant on flat bands, strong interactions and
SYK physics, and by the Swiss National Science Foundation, grants No. 172543
and P2ELP2_175278. Computations have been performed at the Swiss National
Supercomputing Centre (CSCS) under project s1008 and the facilities of
Scientific IT and Application Support Center of EPFL.

## References

  * Thouless  _et al._ (1982) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Physical Review Letters 49, 405 (1982).
  * Haldane (1988) F. D. M. Haldane, Physical Review Letters 61, 2015 (1988).
  * Sondhi  _et al._ (1993) S. L. Sondhi, A. Karlhede, S. Kivelson, and E. Rezayi, Physical Review B 47, 16419 (1993).
  * Fertig  _et al._ (1994) H. Fertig, L. Brey, R. Côté, and A. MacDonald, Physical Review B 50, 11018 (1994).
  * Lux  _et al._ (2021) F. R. Lux, S. Ghosh, P. Prass, E. Prodan, and Y. Mokrousov, preprint arXiv:2103.01047 (2021).
  * Khalaf  _et al._ (2021) E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Science Advances 7, eabf5299 (2021).
  * Chatterjee  _et al._ (2020) S. Chatterjee, N. Bultinck, and M. P. Zaletel, Physical Review B 101, 165141 (2020).
  * Zhang  _et al._ (2019a) Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T. Senthil, Physical Review B 99, 075127 (2019a).
  * Abanov and Wiegmann (2001) A. Abanov and P. B. Wiegmann, Physical Review Letters 86, 1319 (2001).
  * Grover and Senthil (2008) T. Grover and T. Senthil, Physical Review Letters 100, 156804 (2008).
  * Cao  _et al._ (2018) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018).
  * 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, Science 365, 605 (2019).
  * Serlin  _et al._ (2020) M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Science 367, 900 (2020).
  * Fertig  _et al._ (1997) H. Fertig, L. Brey, R. Côté, A. MacDonald, A. Karlhede, and S. Sondhi, Physical Review B 55, 10671 (1997).
  * Wilczek and Zee (1983) F. Wilczek and A. Zee, Phys. Rev. Lett. 51, 2250 (1983).
  * Hsu and Chakravarty (2013) C.-H. Hsu and S. Chakravarty, Physical Review B 87, 085114 (2013).
  * Nagaosa and Tokura (2013) N. Nagaosa and Y. Tokura, Nature Nanotechnology 8, 899 (2013).
  * Tokura and Kanazawa (2020) Y. Tokura and N. Kanazawa, Chemical Reviews 121, 2857 (2020).
  * Wu and Sarma (2020) F. Wu and S. D. Sarma, Physical Review B 102, 165118 (2020).
  * Bömerich  _et al._ (2020) T. Bömerich, L. Heinen, and A. Rosch, Physical Review B 102, 100408 (2020).
  * Tarnopolsky  _et al._ (2019) G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Physical Review Letters 122, 106405 (2019).
  * Popov and Milekhin (2021) F. K. Popov and A. Milekhin, Physical Review B 103, 155150 (2021).
  * San-Jose  _et al._ (2012) P. San-Jose, J. González, and F. Guinea, Physical Review Letters 108, 216802 (2012).
  * Ledwith  _et al._ (2020) P. J. Ledwith, G. Tarnopolsky, E. Khalaf, and A. Vishwanath, Physical Review Research 2, 023237 (2020).
  * Mühlbauer  _et al._ (2009) S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Science 323, 915 (2009).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Kruchkov (2021) A. Kruchkov, preprint arXiv:2105.14672 (2021).
  * Soluyanov and Vanderbilt (2011) A. A. Soluyanov and D. Vanderbilt, Phys. Rev. B 83, 035108 (2011).
  * King-Smith and Vanderbilt (1993) R. King-Smith and D. Vanderbilt, Physical Review B 47, 1651 (1993).
  * Resta (1993) R. Resta, Europhysics Letters 22, 133 (1993).
  * Coh and Vanderbilt (2009) S. Coh and D. Vanderbilt, Physical review letters 102, 107603 (2009).
  * Hofstadter (1976) D. R. Hofstadter, Physical Review B 14, 2239 (1976).
  * Wannier (1978) G. Wannier, physica status solidi (b) 88, 757 (1978).
  * Dean  _et al._ (2013) C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, _et al._ , Nature 497, 598 (2013).
  * Hunt  _et al._ (2013) B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, _et al._ , Science 340, 1427 (2013).
  * Ponomarenko  _et al._ (2013) L. Ponomarenko, R. Gorbachev, G. Yu, D. Elias, R. Jalil, A. Patel, A. Mishchenko, A. Mayorov, C. Woods, J. Wallbank, _et al._ , Nature 497, 594 (2013).
  * Herzog-Arbeitman  _et al._ (2020) J. Herzog-Arbeitman, Z.-D. Song, N. Regnault, and B. A. Bernevig, Physical Review Letters 125, 236804 (2020).
  * Zhang  _et al._ (2019b) Y.-H. Zhang, H. C. Po, and T. Senthil, Physical Review B 100, 125104 (2019b).
  * Lu  _et al._ (2019) X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, _et al._ , Nature 574, 653 (2019).
  * Note (1) The question of the Landau fan (LF) in TBG is a hard one. The original studiesCao  _et al._ (2018) reported the Hall sequence, in units $e^{2}/h$: $\pm 4$, $\pm 8$, $\pm 12$ at neutrality, instead of expected $\pm 4$, $\pm 12$…Moon and Koshino (2012), thus puzzling theorists. The consequent experimentsLu  _et al._ (2019) on more homogeneous samples brought even further surprise with the QH sequence $\pm 2$, $\pm 4$, $\pm 8$. Moreover, LF with $\pm 2$, $\pm 4$ was reported in Ref.Sharpe  _et al._ (2019); LF with $\pm 2$, $\pm 4$, $\pm 6$ was reported in TBG stabilized by WSe2;Arora  _et al._ (2020) and LF with $-2$, $-4$ in the sample of Ref.Pierce  _et al._ (2021). For the high-field behavior of Landau fans, we refer readers to the recent observation of correlated Chern insulators in TBG e.g. in Refs. Pierce  _et al._ (2021); Xie  _et al._ (2021).
  * Nam and Koshino (2017) N. N. T. Nam and M. Koshino, Physical Review B 96, 075311 (2017).
  * Plimpton (1995) S. Plimpton, Journal of Computational Physics 117, 1 (1995).
  * Gargiulo and Yazyev (2017) F. Gargiulo and O. V. Yazyev, 2D Materials 5, 015019 (2017).
  * Los and Fasolino (2003) J. H. Los and A. Fasolino, Physical Review B 68, 024107 (2003).
  * Taherinejad  _et al._ (2014) M. Taherinejad, K. F. Garrity, and D. Vanderbilt, Physical Review B 89, 115102 (2014).
  * Wu  _et al._ (2018) Q. Wu, S. Zhang, H.-F. Song, M. Troyer, and A. A. Soluyanov, Computer Physics Communications 224, 405 (2018).
  * Wannier (1962) G. H. Wannier, Reviews of Modern Physics 34, 645 (1962).
  * Gresch  _et al._ (2017) D. Gresch, G. Autes, O. V. Yazyev, M. Troyer, D. Vanderbilt, B. A. Bernevig, and A. A. Soluyanov, Physical Review B 95, 075146 (2017).
  * Marzari  _et al._ (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Reviews of Modern Physics 84, 1419 (2012).
  * Moon and Koshino (2012) P. Moon and M. Koshino, Physical Review B 85, 195458 (2012).
  * Arora  _et al._ (2020) H. S. Arora, R. Polski, Y. Zhang, A. Thomson, Y. Choi, H. Kim, Z. Lin, I. Z. Wilson, X. Xu, J.-H. Chu, _et al._ , Nature 583, 379 (2020).
  * Pierce  _et al._ (2021) A. T. Pierce, Y. Xie, J. M. Park, E. Khalaf, S. H. Lee, Y. Cao, D. E. Parker, P. R. Forrester, S. Chen, K. Watanabe, _et al._ , preprint arXiv:2101.04123 (2021).
  * Xie  _et al._ (2021) Y. Xie, A. T. Pierce, J. M. Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, _et al._ , preprint arXiv:2107.10854 (2021).