# Interaction-driven Spontaneous Ferromagnetic Insulating States with Odd
Chern Numbers

Peizhi Mai  Department of Physics and Institute of Condensed Matter Theory,
University of Illinois Urbana-Champaign, Urbana, IL 61801, USA  Edwin W. Huang
Department of Physics and Institute of Condensed Matter Theory, University of
Illinois Urbana-Champaign, Urbana, IL 61801, USA  Jiachen Yu  Department of
Applied Physics, Stanford University, Stanford, CA 94305, USA  Geballe
Laboratory of Advanced Materials, Stanford, CA 94305, USA  Benjamin E. Feldman
Geballe Laboratory of Advanced Materials, Stanford, CA 94305, USA  Department
of Physics, Stanford University, Stanford, CA 94305, USA  Stanford Institute
for Materials and Energy Sciences, SLAC National Accelerator Laboratory, Menlo
Park, CA 94025, USA  Philip W. Phillips  Department of Physics and Institute
of Condensed Matter Theory, University of Illinois Urbana-Champaign, Urbana,
IL 61801, USA<EMAIL_ADDRESS>
###### Abstract

Motivated by recent experimental work on moiré systems in a strong magnetic
field, we compute the compressibility as well as the spin correlations and
Hofstadter spectrum of spinful electrons on a honeycomb lattice with Hubbard
interactions using the determinantal quantum Monte Carlo method. While the
interactions in general preserve quantum and anomalous Hall states, emergent
features arise corresponding to an antiferromagnetic insulator at half-filling
and other incompressible states following the Chern sequence $\pm(2N+1)$.
These odd integer Chern states exhibit strong ferromagnetic correlations and
arise spontaneously without any external mechanism for breaking the spin-
rotation symmetry. Analogs of these magnetic states should be observable in
general interacting quantum Hall systems. In addition, the interacting
Hofstadter spectrum is qualitatively similar to the experimental data at
intermediate values of the on-site interaction.

## Introduction

The hallmark of non-trivial topology of filled bands is a quantized Hall
conductance, an integer multiple of the quantum of conductance. The integer is
set by the Chern number[1]. Typically, deviations from integer Chern numbers
indicate that electron-electron interactions are important, as in the
fractional quantum Hall effect. However, there are several examples of
physical systems in which interactions dominate but the Chern number is still
an integer. One such mechanism that involves spin polarization and its
generalizations is quantum Hall ferromagnetism[2, 3, 4, 5, 6]. Moiré systems
in a magnetic field provide a second example in which symmetry-broken quantum
Hall insulators appear at high magnetic fields [7, 8, 9, 10, 11]. The extent
to which these phenomena are generic beyond graphene-based systems and
independent of lattice geometry is unknown.

Motivated by these phenomena, we report here a series of insulating states on
the honeycomb and square lattice which are driven by interactions. The series
we report has odd integer Chern numbers, $\pm 1,3,5,\cdots$. An analysis of
the spin correlations suggests that the spin rotation symmetry is
spontaneously broken resulting in ferromagnetism. Our observations here add
intrigue to the mixed role played by topology and interactions in 2D
materials. Our simulations reveal that explicit single-particle symmetry
breaking such as Zeeman splitting is not required and the full ferromagnetic
sequence arises spontaneously in a general bipartite lattice.

The evolution of electronic states in a perpendicular magnetic field has a
long history. Because a magnetic field preserves the crystal momentum, the
single-particle energy spectrum for non-interacting electrons is easily
obtained[12, 13, 14] by replacing the momentum ${\bf p}$ with ${\bf p}-e{\bf
A}/c$ where ${\bf A}$ is the magnetic vector potential, $e$ the electron
charge and $c$ the speed of light. In 2D, the resultant Hofstadter[14]
spectrum adequately describes the evolution of the tight-binding electronic
states as a function of the magnetic flux. Hidden in the wings of the
underlying butterfly spectrum are gapped states indexed by Chern numbers which
fix[1] the quantization of the Hall conductance. Moiré systems as in the case
of MATBG offer a new route to engineering gaps in the electronic spectrum
through the competition between band filling and the interaction energy[15,
16, 17, 18]. As the twist angle controls[19] the ratio of the kinetic to the
potential energy and leads to a complete quenching of the kinetic energy at
the magic angle, moiré systems in a magnetic field offer the ultimate
playground for studying the physics from the interplay between strong
correlation and magnetic field. With the kinetic energy quenched, moiré
systems encode the evolution of the Hofstadter spectrum in the presence of
strong interactions. This is currently an unsolved non-trivial problem.

This problem is complicated by the fact that the simple replacement of the
momentum by ${\bf p}-e{\bf A}/c$ fails in the presence of interactions because
interactions in general mix crystal momenta as in the case of the Hubbard
interaction. Consequently, while theoretical efforts have addressed certain
limits of the interacting Hofstadter problem[20, 21, 22, 23, 24, 25, 26, 27,
28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39], no analytical method exists
to determine the complete spectrum in a magnetic field in the presence of
interactions. Nonetheless, this is an urgent problem in condensed matter
physics given the plethora of experiments on MATBG and related systems that
are focused on revealing the low-energy physics resulting from the interplay
between a magnetic field and strong correlation.

Theoretically, one has three options: 1) phenomenology, 2) some type of mean-
field theory, dynamical[32, 33] or otherwise[34, 35, 36, 37] or 3) serious
numerics which so far have been limited to exact diagonalization[38, 39] on
few-particle systems. We pursue the last option in this paper as no benchmarks
have been established for even the simplest model of interacting electrons on
any of the lattices relevant to either MATBG or the transition metal
dichalcogenide systems. We focus on spinful fermions primarily on a honeycomb
lattice including only nearest neighbor hopping and Hubbard interactions under
an external magnetic field, and perform a determinantal quantum Monte Carlo
(DQMC) simulation for all densities and magnetic fluxes. DQMC is an unbiased
and numerically exact method to capture the full quantum fluctuations for
correlated systems. In a prior work with Hubbard interactions, features such
as the local compressibility and other thermodynamic quantities were
calculated using DQMC as a function of the magnetic flux for the square
lattice[27] and no ferromagnetism was reported. We focus here on a honeycomb
lattice as it is closer to the underlying geometry of most existing moiré
systems. In general, we find that the interactions preserve the integer
quantum and anomalous Hall states of the non-interacting system. However, the
interactions do generate an antiferromagnetic insulating state at half-
filling, as expected, and also emergent interaction-driven insulating states
in both of the honeycomb and square lattices. The Chern sequence for these
states is $\pm(2N+1)$. All such states exhibit strong ferromagnetic
correlations. This represents a numerically exact evidence for such
interaction-driven states in full density region based on the Hubbard
interaction.

## Results

### Non-interacting quantum Hall effects

To begin with, we present the non-interacting charge compressibility
$\chi=\partial\langle n\rangle/\partial\mu$ in Fig. 1 as a function of
magnetic flux and electron density and compare with the results for the square
lattice at $\beta=20/t$. In both the square and honeycomb lattices, particle-
hole symmetry obtains as both are bipartite and the model contains only
nearest neighbor hopping. The straight lines in Fig. 1 correspond to solutions
to the Diophantine equation [1]

$\displaystyle\langle n\rangle=r\frac{\phi}{\phi_{0}}+s,$ (1)

in which $\langle n\rangle=\langle N_{\text{e}}\rangle/N_{c}$ ($N_{c}$ is the
number of unit cells), $r$ is an integer given by the inverse slope of the
straight lines and $s$ is the offset given by the intercept. $r$ defines the
Chern number. We have chosen to plot the filling from $[0,4]$ to take into
account the spin and sublattice degeneracy in the honeycomb lattice but from
$[0,2]$ for the square lattice in which only a spin degeneracy exists. Hence,
there is only a factor of $2$ in translating the densities between the two
systems. In both Fig. 1 panels (a) and (b), the Diophantine lines[1] starting
from the bottom left (right) corners have $r=\pm 2N$ (the factor of 2 accounts
for spin degeneracy) with $N=1,2\cdots$ corresponding to spin-unpolarized
quantum Hall states, while in only Fig. 1(a), the lines that start from half-
filling ($\langle n\rangle=2$) at zero-field and have $r=\pm 4(N+1/2)$ (the
factor of 4 accounts for spin and sub-lattice degeneracy) with $N=0,1,2\cdots$
indicating the anomalous quantum Hall effect[40, 41, 42]. Thus, the honeycomb
lattice is more closely aligned to the physics observed in MATBG than is the
square lattice.

### Turning on interactions

Next, we explore how interactions change this pattern in the honeycomb
lattice. We use DQMC to calculate the compressibility,

$\displaystyle\chi=\beta\chi_{c}=\frac{\beta}{N}\sum_{{\bf i},{\bf
j}}\left[\langle n_{\bf i}n_{\bf j}\rangle-\langle n_{\bf i}\rangle\langle
n_{\bf j}\rangle\right],$ (2)

in the presence of Hubbard interactions, where $\chi_{c}$ is the charge
correlation function. Due to the Fermionic sign problem, we are only able to
calculate the compressibility for the full density and flux region for a
system size $N_{\text{site}}=6\times 6\times 2$ with an interaction strength
up to $U/t=4$ and temperature as low as $T/t=0.125$ (or $\beta=8t^{-1}$). In
the first row of Fig. 2, as $U$ increases, the non-interacting lines in the
compressibility are softened and a middle vertical line appears as a single-
particle gap develops. At the largest $U/t=4$, we can still observe the
dominant and sub-dominant lines, indicating the resilience of the non-
interacting quantum Hall effect against interactions. Note that the lines not
merging at $\phi/\phi_{0}=0,0.5$ in Fig. 2(b,c) is due to strong finite size
effects (see Supplementary Figure 1 and 2) and thus disregarded in Fig.
2(e,f). Also of note is the emergence of a new feature at $\langle n\rangle=2$
for $U/t=4$. This corresponds to a dip in the density of states, a precursor
to the Mott gap[43, 44]. The second row of these figures corresponds to a
simulation at the lower temperature of $\beta=20/t$ for $U=0$ and $U/t=2$.
Fig. 2(d) shows a sharpening of the Diophantine features as the temperature is
lowered by more than a factor of two to $\beta=20/t$. At $U/t=2$, the vertical
line at $\langle n\rangle=2$ possibly indicates a gap opening for finite
field. At this temperature, even for such a modest value of $U$, the
suppression of the density of states is evident and the corresponding
insulator is antiferromagnetic (see Supplementary Figure 6). Also of note is
the state indicated by the solid red line in Fig. 2(c,e). This line evolves as
a function of the magnetic flux with a slope of unity. Such a state is absent
from the non-interacting sequence as it has a Chern number of $\pm 1$. Notably
this state is visible in the map of the spin correlation (Fig. 2(f)),

$\displaystyle\chi_{s}=\sum_{\bf r}S({\bf r})=\frac{1}{N}\sum_{{\bf i},{\bf
r}}\langle S^{z}_{\bf i}S^{z}_{{\bf i}+{\bf r}}\rangle.$ (3)

In fact, the light features in Fig. 2(f) reveal that these possibly
incompressible states all have a markedly enhanced spin correlation with
distinct slopes as a function of magnetic flux and density which differ from
those of the non-interacting system even at much lower temperatures (see
Supplementary Figure 4). It is the origin of these states that is the
principal focus of this paper.

We gain further insight into the possible emergence of incompressible states
by taking slices through the charge correlation, $\chi_{c}=\chi/\beta$, at
particular values of the magnetic flux. In Fig. 3(a-c), we show the charge
correlation explicitly at $\phi/\phi_{0}=11/36$, chosen to avoid finite-size
effect (see Supplementary Figure 1 and 2) for the honeycomb lattice and
$\phi/\phi_{0}=2/9$ for the square lattice. For the honeycomb lattice, Fig.
3(a,b) illustrates that as the temperature is lowered from $\beta=8/t$ to
$\beta=20/t$, the dual-dip feature for $U=0$ in the vicinity of $\langle
n\rangle=2$ gives rise to a full quantum Hall state. With the interaction
increasing to $U=4t$, this dual-dip feature contains a depression precisely at
$\langle n\rangle=2$. This is the incompressible Mott gap at $\langle
n\rangle=2$ [43, 44]. Panel Fig. 3(c) displays the analogous trend for the
square lattice. Away from half-filling, we find several emergent states which
have a suppressed charge correlation indicating the possible onset of a gap.
The states occur around fillings of $\langle n\rangle=11/36$ and $11/12$. The
same is true of their particle-hole equivalents. It is precisely the first of
such states that is highlighted in red in Fig. 2(e). The analogous states are
also present for the square lattice in Fig. 3(c). All such dips in $\chi_{c}$
are enhanced as the interaction strength increases. Further, such behavior
persists even as the system size increases (see Supplementary Figure 3). To
uncover the possible cause of these states, we focus on the spin
susceptibility $\chi_{s}$. For $U=0$, $\chi_{c}=4\chi_{s}$. However, the
second row of Fig. 3 shows that whenever the charge correlation exhibits an
interaction-driven dip, the spin correlation shows a peaked structure. Fig.
3(e) shows that as the temperature is lowered, the peak of the spin
correlation increases at the fillings where the charge correlation develops a
dip.

Now we look at the full density- and magnetic-flux-dependent spin correlation
$\chi_{s}$ for $U/t=2$ and $U/t=4$ in Fig. 4(a) and (b) respectively at their
lowest temperatures. Straight lines with inverse slope corresponding to Chern
number $C=\pm 1,3,5$ are plotted and found to be aligned with the ridges of
the spin correlation. We choose one representative point (not the brightest)
at each line and study its temperature evolution at different $U$, presented
in Fig. 4 (c-e). In all cases, when $U=0$, the spin correlation deceases along
with temperature. However, for a finite $U$, the spin correlation blows up as
temperature decreases (below a critical temperature for Fig. 4 (d) and (e)).
The ultimate spin state is revealed from a spatial map of the real-space
static spin susceptibility:

$\displaystyle S({\bf r},\omega=0)=\frac{1}{N}\int_{0}^{\beta}\sum_{{\bf
i}}\langle S^{z}_{{\bf i}+{\bf r}}(\tau)S^{z}_{\bf i}(0)\rangle d\tau,$ (4)

presented in Fig. 4 (f-h), at $U/t=2$ and the lowest temperature ($\beta=30/t$
as circled in Fig. 4 (c-e) respectively). This quantity is more sensitive in
detecting fluctuating order at finite temperature than the zero-time spin
correlation[45]. The color map signifies positive spin correlation across the
lattice relative to the site at the origin. Such same-sign correlations are
indicative of ferromagnetism. Fig. 4(f) for a Chen number $C=1$ exhibits
strong ferromagnetic susceptibility. Fig. 4(g) with a Chen number $C=3$ also
displays a clear ferromagnetic pattern. Fig. 4(h) corresponding to Chen number
$C=5$ reveals an evident tendency towards ferromagnetism at lower $T$, though
not fully ferromagnetic as in the other cases. In addition, we also expect the
$C=\pm 1$ ferromagnetic states to exist in the middle of Fig. 4 (a,b) (as
depicted by the dashed line) with extrapolation to $\langle n\rangle=2$ at
zero flux. But limited by the sign problem, we have yet been not able to
investigate low enough temperature to unearth a clear ferromagnetic pattern at
these densities.

The full picture is now apparent. The charge dips and enhanced spin
correlations, which have no counterpart in the non-interacting system in Fig.
2, correspond to ferromagnetic insulators with odd integer Chern numbers. Both
the spin correlations and magnitude of the charge gap are enhanced as the
temperature is lowered. The same trend holds for the square lattice as is
evident from Fig. 3(c,f). This behavior matches expectations from quantum Hall
ferromagnetism, but here we show that local interactions are sufficient to
induce odd Chern integer states on both the honeycomb and square lattices. We
are led to the conclusion that such insulating states are generically present
in bipartite lattices with interactions with no need for fine-tuning or
single-particle splitting. While there is some indication of charge ordering
(see Supplementary Figure 5), it is not compelling at this level of study and
hence we leave this for a future publication.

We finally display the benchmark calculation of the Hofstadter spectrum as
defined by the local density of states, the quantity directly measured
experimentally. Our focus in Fig. 5 is at half-filling and $U/t=2,4$, obtained
from constructing an analytic continuation with DEAC on the DQMC local Green
function. The comparison between DEAC and the analytical result at $U=0$ (see
Supplementary Figure 7) gives us some idea about the resolution of DEAC and
offers a guide as to how to interpret the interacting system results. Since
there is no sign problem at half-filling, we are able to conduct the
calculation at a low temperature ($\beta=30/t$). Panels Fig. 5(a,b) show how
the antiferromagnetic gap comes into full view by $U/t=4$ and some hint of it
appears already at the modest value of $U/t=2$ only with finite magnetic
field. The gap at $U/t=2$ is most likely to be of the Slater type[46, 47]
because the interaction strength is only around $1/3$ of the bare bandwidth
and the insulating state appears at much lower temperature than that required
for the formation of antiferromagnetic correlation. On the other hand, the gap
at $U/t=4$ is closer to a Mott gap because it is established at a much high
temperature (shown in Fig. 2(c)), consistent with previous studies on the
Hubbard model in honeycomb lattice[43, 44]. While the corresponding
experimental figure is at variable filling[18], which is inaccessible because
of the sign problem, the overall features qualitatively reproduce the
experimental results for moderate values of $U/t=2$.

## Discussion

We have studied here the evolution of the excitation spectrum of a Hubbard-
interacting electron gas in the presence of a strong perpendicular magnetic
field on bipartite lattices. We have shown that while the interactions
preserve the non-interacting integer quantum and anomalous Hall states, new
states do emerge from the interactions. In addition to the antiferromagnetic
gap at half-filling, we have discovered a series of odd-integer Chern
insulating ferromagnetic states which exhibit enhanced positive spin
correlations as the temperature is lowered. In light of ferromagnetism as the
underlying cause of the insulating states, that the Chern number is odd is
easily understood. Our work suggests ferromagnetism is generic requiring just
modest magnetic fields and strong interactions to generate the full sequence
of odd-integer states. Note while the sequence we observe departs from the
middle or the edge of the band, an odd-integer sequence emanating from the
$\langle n\rangle=1$ or $\langle n\rangle=3$ fillings on the honeycomb lattice
would require spin-orbit coupling in the Hamiltonian, thereby generalizing the
utility of this work.

## Methods

We study the Hofstadter-Hubbard model on a honeycomb lattice,

$\displaystyle\begin{aligned} H&=-t\sum_{\langle{\bf i}{\bf
j}\rangle\sigma}\exp(i\phi_{{\bf i},{\bf j}})c^{\dagger}_{{\bf
i}\sigma}c^{\phantom{\dagger}}_{{\bf j}\sigma}-\mu\sum_{{\bf i},\sigma}n_{{\bf
i}\sigma}\\\ &+U\sum_{{\bf i}}(n_{{\bf i}\uparrow}-\frac{1}{2})(n_{{\bf
i}\downarrow}-\frac{1}{2}),\end{aligned}$ (5)

where $t$ represents the nearest neighbor hopping; $c^{\dagger}_{{\bf
i}\sigma}$ ($c_{{\bf i}\sigma}$) creates (annihilates) an electron with spin
$\sigma$ at site ${\bf i}$, $\mu$ is the chemical potential, $U$ is the on-
site interaction. Due to the presence of a uniform magnetic field, we use the
Peierls substitution[12] to introduce the phase through the flux threading,

$\displaystyle\phi_{{\bf i},{\bf j}}=\frac{2\pi}{\phi_{0}}\int_{r_{\bf
i}}^{r_{\bf j}}{\bf A}\cdot d{\bf l},$ (6)

where $\phi_{0}=h/e$ in the hopping term is a result of the quantized magnetic
field and the integration is over the straight line path from site ${\bf i}$
to ${\bf j}$.

We simulate this Hamiltonian Eq. (5) on a finite cluster
$N_{\text{site}}=2L^{2}$. The honeycomb lattice contains two sub-lattices,
which explains the factor of $2$, with lattice constant $a=1$ and $L$ the
number of site along either lattice basis respectively for each sub-lattice.
We adjust the modified periodic boundary conditions in Ref. [ 48] to the
honeycomb lattice. To obtain a single-value wave function requires the flux
quantization condition $\phi/\phi_{0}=n_{f}/N_{c}$ with $n_{f}$ an integer and
$N_{c}=L^{2}$ the number of unit cells. The symmetric gauge ${\bf
A}=(x\hat{y}-y\hat{x})B/2$ is chosen for this calculation.

We apply DQMC[49, 50, 51] to this model Eq. (5) and calculate the
compressibility and Green function. The jackknife resampling is used to
estimate the standard error of the mean as error bars in DQMC results. With
the local Green function, we compute the local density of states using the
recently developed Differential Evolution for Analytic Continuation (DEAC)
[52].

## Data Availability

The data for this study is available at
https://zenodo.org/record/7608167#.Y-AxyezML6g.

## Code Availability

The DQMC code used for this project can be obtained at
https://github.com/edwnh.

## Acknowledgements

BEF, PM, PWP, and JY acknowledge support for the computation and conception of
this project from Quantum Sensing and Quantum Materials (QSQM), an Energy
Frontier Research Center funded by the US Department of Energy (DOE), Office
of Science, Basic Energy Sciences (BES), under award no. DE-SC0021238. E.W.H.
was supported by the Gordon and Betty Moore Foundation EPiQS Initiative
through the grants GBMF 4305 and GBMF 8691. This work used an analysis of
insulating states in MATBG funded through DMR-2111379 for his work on MATBG.
PM thanks Nathan Nichols for the help on DEAC simulation. The DQMC calculation
of this work used the Advanced Cyberinfrastructure Coordination Ecosystem:
Services & Support (ACCESS) Expanse supercomputer through the research
allocation TG-PHY220042, which is supported by National Science Foundation
grant number ACI-1548562[53].

## Competing interests

The authors declare no competing financial or non-financial interests.

## Author Contributions

E. W. H. and P. M. developed the DQMC code. P. M. carried out the
calculations. All authors analyzed the results and wrote the manuscript. B. E.
F. and P. W. P. supervised the project.

## References

  * [1] Thouless, D. J., Kohmoto, M., Nightingale, M. P. & den Nijs, M.  Quantized hall conductance in a two-dimensional periodic potential.  _Phys. Rev. Lett._ 49, 405–408 (1982).
  * [2] Sondhi, S. L., Karlhede, A., Kivelson, S. A. & Rezayi, E. H.  Skyrmions and the crossover from the integer to fractional quantum hall effect at small zeeman energies.  _Phys. Rev. B_ 47, 16419–16426 (1993).
  * [3] Manfra, M., Goldberg, B., Pfeiffer, L. & West, K.  Skyrmions and the $\nu=1$ quantum hall ferromagnet.  _Acta Physica Polonica A_ 92, 621 (1997).
  * [4] Arovas, D. P., Karlhede, A. & Lilliehöök, D.  $\mathrm{SU}(n)$ quantum hall skyrmions.  _Phys. Rev. B_ 59, 13147–13150 (1999).
  * [5] Ezawa, Z. F. & Hasebe, K.  Interlayer exchange interactions, $\mathrm{SU}(4)$ soft waves, and skyrmions in bilayer quantum hall ferromagnets.  _Phys. Rev. B_ 65, 075311 (2002).
  * [6] Nomura, K. & MacDonald, A. H.  Quantum hall ferromagnetism in graphene.  _Phys. Rev. Lett._ 96, 256602 (2006).
  * [7] Spanton, E. M. _et al._ Observation of fractional chern insulators in a van der waals heterostructure.  _Science_ 360, 62–66 (2018).
  * [8] Wang, L. _et al._ Evidence for a fractional fractal quantum hall effect in graphene superlattices.  _Science_ 350, 1231–1234 (2015).
  * [9] Dean, C. R. _et al._ Hofstadter’s butterfly and the fractal quantum hall effect in moirésuperlattices.  _Nature_ 497, 598–602 (2013).
  * [10] Yu, G. L. _et al._ Hierarchy of hofstadter states and replica quantum hall ferromagnetism in graphene superlattices.  _Nature Physics_ 10, 525–529 (2014).
  * [11] Hunt, B. _et al._ Massive dirac fermions and hofstadter butterfly in a van der waals heterostructure.  _Science_ 340, 1427–1430 (2013).
  * [12] Peierls, R.  Zur theorie des diamagnetismus von leitungselektronen.  _Z. Physik_ 80, 763–791 (1933).
  * [13] Harper, P. G.  Single band motion of conduction electrons in a uniform magnetic field.  _Proceedings of the Physical Society. Section A_ 68, 874–878 (1955).
  * [14] Hofstadter, D. R.  Energy levels and wave functions of bloch electrons in rational and irrational magnetic fields.  _Phys. Rev. B_ 14, 2239–2249 (1976).
  * [15] Saito, Y. _et al._ Hofstadter subband ferromagnetism and symmetry-broken chern insulators in twisted bilayer graphene.  _Nature Physics_ 17, 478–481 (2021).
  * [16] Park, J. M., Cao, Y., Watanabe, K., Taniguchi, T. & Jarillo-Herrero, P.  Flavour hund’s coupling, chern gaps and charge diffusivity in moirégraphene.  _Nature_ 592, 43–48 (2021).
  * [17] Choi, Y. _et al._ Correlation-driven topological phases in magic-angle twisted bilayer graphene.  _Nature_ 589, 536–541 (2021).
  * [18] Yu, J. _et al._ Correlated hofstadter spectrum and flavour phase diagram in magic-angle twisted bilayer graphene.  _Nature Physics_ 18, 825–831 (2022).
  * [19] Bistritzer, R. & MacDonald, A. H.  Moiré bands in twisted double-layer graphene.  _Proceedings of the National Academy of Sciences_ 108, 12233–12237 (2011).
  * [20] Gudmundsson, V. & Gerhardts, R. R.  Effects of screening on the Hofstadter butterfly.  _Physical Review B_ 52, 16744–16752 (1995).
  * [21] Pfannkuche, D. & MacDonald, A. H.  Quantum Hall effect of interacting electrons in a periodic potential.  _Physical Review B_ 56, R7100–R7103 (1997).
  * [22] Doh, H. & Salk, S.-H. S.  Effects of electron correlations on the Hofstadter spectrum.  _Physical Review B_ 57, 1312–1315 (1998).
  * [23] Czajka, K., Gorczyca, A., Maśka, M. M. & Mierzejewski, M.  Hofstadter butterfly for a finite correlated system.  _Physical Review B_ 74, 125116 (2006).
  * [24] Apalkov, V. M. & Chakraborty, T.  Gap Structure of the Hofstadter System of Interacting Dirac Fermions in Graphene.  _Physical Review Letters_ 112, 176401 (2014).
  * [25] Mishra, A., Hassan, S. R. & Shankar, R.  Effects of interaction in the Hofstadter regime of the honeycomb lattice.  _Physical Review B_ 93, 125134 (2016).
  * [26] Andrews, B. & Soluyanov, A.  Fractional quantum Hall states for moir\’e superstructures in the Hofstadter regime.  _Physical Review B_ 101, 235312 (2020).
  * [27] Ding, J. K. _et al._ Thermodynamics of correlated electrons in a magnetic field.  _Communications Physics_ 5, 204 (2022).
  * [28] Shaffer, D., Wang, J. & Santos, L. H.  Theory of hofstadter superconductors.  _Phys. Rev. B_ 104, 184501 (2021).
  * [29] Shaffer, D., Wang, J. & Santos, L. H.  Unconventional self-similar hofstadter superconductivity from repulsive interactions.  _Nature Communications_ 13, 7785 (2022).
  * [30] Wang, X. & Vafek, O.  Narrow bands in magnetic field and strong-coupling hofstadter spectra.  _Phys. Rev. B_ 106, L121111 (2022).
  * [31] Herzog-Arbeitman, J., Chew, A., Efetov, D. K. & Bernevig, B. A.  Reentrant correlated insulators in twisted bilayer graphene at 25 t ($2\pi$ flux).  _Phys. Rev. Lett._ 129, 076401 (2022).
  * [32] Acheche, S., Arsenault, L.-F. & Tremblay, A.-M. S.  Orbital effect of the magnetic field in dynamical mean-field theory.  _Phys. Rev. B_ 96, 235135 (2017).
  * [33] Markov, A. A., Rohringer, G. & Rubtsov, A. N.  Robustness of the topological quantization of the hall conductivity for correlated lattice electrons at finite temperatures.  _Phys. Rev. B_ 100, 115102 (2019).
  * [34] Gudmundsson, V. & Gerhardts, R. R.  Effects of screening on the hofstadter butterfly.  _Phys. Rev. B_ 52, 16744–16752 (1995).
  * [35] Doh, H. & Salk, S.-H. S.  Effects of electron correlations on the hofstadter spectrum.  _Phys. Rev. B_ 57, 1312–1315 (1998).
  * [36] Tu, W.-L., Schindler, F., Neupert, T. & Poilblanc, D.  Competing orders in the hofstadter $t--j$ model.  _Phys. Rev. B_ 97, 035154 (2018).
  * [37] Zang, J., Wang, J., Cano, J. & Millis, A. J.  Hartree-fock study of the moiré hubbard model for twisted bilayer transition metal dichalcogenides.  _Phys. Rev. B_ 104, 075150 (2021).
  * [38] Barelli, A., Bellissard, J., Jacquod, P. & Shepelyansky, D. L.  Double butterfly spectrum for two interacting particles in the harper model.  _Phys. Rev. Lett._ 77, 4752–4755 (1996).
  * [39] Czajka, K., Gorczyca, A., Maśka, M. M. & Mierzejewski, M.  Hofstadter butterfly for a finite correlated system.  _Phys. Rev. B_ 74, 125116 (2006).
  * [40] Semenoff, G. W.  Condensed-matter simulation of a three-dimensional anomaly.  _Phys. Rev. Lett._ 53, 2449–2452 (1984).
  * [41] Novoselov, K. S. _et al._ Two-dimensional gas of massless dirac fermions in graphene.  _Nature_ 438, 197–200 (2005).
  * [42] Zhang, Y., Tan, Y.-W., Stormer, H. L. & Kim, P.  Experimental observation of the quantum hall effect and berry’s phase in graphene.  _Nature_ 438, 201–204 (2005).
  * [43] Assaad, F. F. & Herbut, I. F.  Pinning the order: The nature of quantum criticality in the hubbard model on honeycomb lattice.  _Phys. Rev. X_ 3, 031010 (2013).
  * [44] Ostmeyer, J. _et al._ Semimetal–mott insulator quantum phase transition of the hubbard model on the honeycomb lattice.  _Phys. Rev. B_ 102, 245105 (2020).
  * [45] Mai, P., Karakuzu, S., Balduzzi, G., Johnston, S. & Maier, T. A.  Intertwined spin, charge, and pair correlations in the two-dimensional hubbard model in the thermodynamic limit.  _Proceedings of the National Academy of Sciences_ 119, e2112806119 (2022).
  * [46] Gull, E., Werner, P., Wang, X., Troyer, M. & Millis, A. J.  Local order and the gapped phase of the hubbard model: A plaquette dynamical mean-field investigation.  _Europhysics Letters_ 84, 37009 (2008).
  * [47] Schäfer, T. _et al._ Fate of the false mott-hubbard transition in two dimensions.  _Phys. Rev. B_ 91, 125109 (2015).
  * [48] Assaad, F. F.  Depleted kondo lattices: Quantum monte carlo and mean-field calculations.  _Phys. Rev. B_ 65, 115104 (2002).
  * [49] Blankenbecler, R., Scalapino, D. J. & Sugar, R. L.  Monte carlo calculations of coupled boson-fermion systems. i.  _Phys. Rev. D_ 24, 2278–2286 (1981).
  * [50] Hirsch, J. E.  Two-dimensional hubbard model: Numerical simulation study.  _Phys. Rev. B_ 31, 4403–4419 (1985).
  * [51] White, S. R. _et al._ Numerical study of the two-dimensional hubbard model.  _Phys. Rev. B_ 40, 506–516 (1989).
  * [52] Nichols, N. S., Sokol, P. & Del Maestro, A.  Parameter-free differential evolution algorithm for the analytic continuation of imaginary time correlation functions.  _Phys. Rev. E_ 106, 025312 (2022).
  * [53] Towns, J. _et al._ Xsede: Accelerating scientific discovery.  _Computing in Science Engineering_ 16, 62–74 (2014).

## Figure Legend

Figure 1: Non-interacting compressibility for square and honeycomb lattice.
Non-interacting compressibility as a function of magnetic flux $\phi$ and
electron density $\langle n\rangle$ for (a) honeycomb lattice
($N_{\text{site}}=36\times 36\times 2$) and (b) square lattice
($N_{\text{site}}=40\times 40$) at $\beta/t=20$. The leading quantum Hall
states are labeled with the corresponding Chern number.  Figure 2:
Compressibility and spin correlation in the presence of interactions.
Compressibility as a function of magnetic flux and electron density in an
$N_{\text{site}}=36\times 36\times 2$ cluster for (a,d) $U/t=0$; an
$N_{\text{site}}=6\times 6\times 2$ cluster for (b,e) $U/t=2$ and (c) $U/t=4$.
The first row is for the inverse temperature $\beta=8/t$ and the second row
for $\beta=20/t$. Panel (f) shows the corresponding spin correlation at
$U/t=2$ and $\beta=20/t$. It has a reverse color bar compared to those of the
charge compressibilty. Figure 3: Compressibility and spin correlation at a
specific magnetic flux. The first row shows the charge correlation at varying
interaction strength with different fixed parameter sets. Panel (a) and (b)
are for the honeycomb lattice at $\phi/\phi_{0}=11/36$, while panel (c) is for
the square lattice at $\phi/\phi_{0}=2/9$. The corresponding temperatures are
(a) $\beta=8/t$, (b) $\beta=20/t$, and (c) $\beta=5/t_{\text{sq}}$ ($t_{sq}$
is the nearest neighbor hopping for the square lattice). The second row shows
the corresponding spin correlation and shares the same legend as does the
first row. Figure 4: Spin correlation and ferromagnetic states. The first row
shows the spin correlation at (a) $U/t=2,\beta=20/t$ and (b)
$U/t=4,\beta=8/t$, with the green lines depicting the odd integer Chern states
aligned with the ridge of the spin correlation. The second row shows the spin
correlation for selected points (marked at each Chern state in panels (a) and
(b)) as a function of temperature under different interaction strengths. The
third row presents the real-space zero-frequency spin susceptibility pattern
for the circled points in the corresponding second row (at $U/t=2$ and the
lowest temperature $\beta=30/t$).  Figure 5: Density of states at half-
filling. The density of states at half-filling ($\langle n\rangle=2$) and
$\beta=30/t$ for (a) $U/t=2$ and (b) $U/t=4$.