Ferromagnetism in narrow bands of moir\'e superlattices

TL;DR

This paper demonstrates that nearly flat narrow bands in various moiré superlattices can host robust ferromagnetic insulators at integer fillings, with stability influenced by band flatness, interactions, and quantum geometry.

Contribution

It provides analytical and numerical evidence for ferromagnetism in flat bands of moiré materials, extending results to finite bandwidth and larger systems using DMRG and exact diagonalization.

Findings

Maximally spin and valley polarized insulators at integer fillings in flat bands.

Interactions can destabilize ferromagnetic states by inducing effective dispersive terms.

Results applicable even when valley Chern number is zero, emphasizing the role of quantum geometry.

Abstract

Many graphene moir\'e superlattices host narrow bands with non-zero valley Chern numbers. We provide analytical and numerical evidence for a robust spin and/or valley polarized insulator at total integer band filling in nearly flat bands of several different moir\'e materials. In the limit of a perfectly flat band, we present analytical arguments in favor of the ferromagnetic state substantiated by numerical calculations. Further, we numerically evaluate its stability for a finite bandwidth. We provide exact diagonalization results for models appropriate for ABC trilayer graphene aligned with hBN, twisted double bilayer graphene, and twisted bilayer graphene aligned with hBN. We also provide DMRG results for a honeycomb lattice with a quasi-flat band and non-zero Chern number, which extend our results to larger system sizes. We find a maximally spin and valley polarized insulator at all integer fillings when the band is sufficiently flat. We also show that interactions may induce effective dispersive terms strong enough to destabilize this state. These results still hold in the case of zero valley Chern number (for example, trivial side of TLG/hBN). We give an intuitive picture based on extended Wannier orbitals, and emphasize the role of the quantum geometry of the band, whose microscopic details may enhance or weaken ferromagnetism in moir\'e materials.

Equations

$$H_V=\sum _{\mathbf{q}}\delta \tilde{\rho }(\mathbf{q})V(\mathbf{q})\delta \tilde{\rho }(-\mathbf{q})$$

$$\tilde{\rho }(\mathbf{q})$$

$${\tilde{\rho }}_s(\mathbf{q})$$

$${\lambda }_s(\mathbf{k}+\mathbf{q},\mathbf{k})=⟨u_{s,\mathbf{k}+\mathbf{q}}|u_{s,\mathbf{k}}⟩$$

$$\delta \tilde{\rho }(\mathbf{q})|\psi ⟩=0$$

$$[\delta {\tilde{\rho }}_{\pm }(\mathbf{q}),\delta {\tilde{\rho }}_{\pm }({\mathbf{q}}^{\prime })=\pm 2i\sin \left(\frac{(\mathbf{q}\times {\mathbf{q}}^{\prime })l_B^2}{2}\right)\delta {\tilde{\rho }}_{\pm }(\mathbf{q}+{\mathbf{q}}^{\prime })$$

$$[\delta \tilde{\rho }(\mathbf{q}),\delta \tilde{\rho }({\mathbf{q}}^{\prime })]=2i\sin \left(\frac{(\mathbf{q}\times {\mathbf{q}}^{\prime })l_B^2}{2}\right)I^z(\mathbf{q}+{\mathbf{q}}^{\prime })$$

$$I^z(\mathbf{q})|\psi ⟩=0,\forall \mathbf{q}\ne 0$$

$$⟨\psi |I^z(-\mathbf{q})I^z(\mathbf{q})|\psi ⟩\underset{\mathbf{q}\to 0}{=}\frac{\sqrt{{\kappa }_v{\rho }_{sv}}}{2}|\mathbf{q}|$$

$$V(q)=U\frac{1}{q}\left(1-e^{-q{r}_0}\right)$$

$$H=H_V-\frac{W_{\mathrm{bare}}}{2}\sum _{\mathbf{k}}\cos \left(\mathbf{k}\cdot {\mathbf{a}}_1+\mathbf{k}\cdot {\mathbf{a}}_2\right)c_{\mathbf{k}}^{†}{c}_{\mathbf{k}}$$

$$\rho_{s}=\frac{\partial E(\theta_{s})}{\partial\theta_{s}}\biggr{\rvert}_{\theta_{s}=0}$$

$$H_1$$

$$f_{\alpha ,I}(\mathbf{k})=\sum _m{u}_{\alpha ,m}^I(\mathbf{k})c_{\alpha ,m}(\mathbf{k})$$

$$\sum _m{u}_{\alpha ,m}^{I\ast }{u}_{\alpha ,m}^J={\delta }_{IJ}$$

$$\sum _I{u}_{\alpha ,m}^{I\ast }{u}_{\alpha ,n}^I={\delta }_{mn}$$

$$H_{\mathrm{int}}=\frac{1}{2}\sum _{{\mathbf{k}}_1,{\mathbf{k}}_2,\mathbf{q}}\sum _{\alpha ,\beta }\sum _{I,J}V(\mathbf{q})f_{\alpha ,I}^{†}({\mathbf{k}}_1+\mathbf{q})f_{\beta ,J}^{†}({\mathbf{k}}_2-\mathbf{q})f_{\beta ,J}({\mathbf{k}}_2)f_{\alpha ,I}({\mathbf{k}}_1)$$

$$H_{\mathrm{int}}$$

$$⟨u_{\alpha ,m_1}({\mathbf{k}}_1+\mathbf{q})⟩u_{\alpha ,m_2}({\mathbf{k}}_1⟨u_{\beta ,n_1}({\mathbf{k}}_2-\mathbf{q})⟩u_{\beta ,n_2}({\mathbf{k}}_2)$$

$${\rho }_{\alpha }(\mathbf{q})=\sum _{m_1,m_2}\sum _{\mathbf{k}}{c}_{\alpha ,m_1}^{†}(\mathbf{k}+\mathbf{q})c_{\alpha ,m_2}(\mathbf{k})⟨u_{\alpha ,m_1}(\mathbf{k}+\mathbf{q})⟩u_{\alpha ,m_2}(\mathbf{k})$$

$$H_{\mathrm{int}}=$$

$$-\frac{1}{2}\sum _{\mathbf{k},\mathbf{q}}\sum _{\alpha }\sum _{m_1,m_2}\sum _{n_2}V(\mathbf{q})c_{\alpha ,m_1}^{†}(\mathbf{k})c_{\alpha ;n_2}(\mathbf{k})⟨u_{\alpha ;m_1}(\mathbf{k})⟩u_{\alpha ;m_2}(\mathbf{k}-\mathbf{q}⟨u_{\alpha ;m_2}(\mathbf{k}-\mathbf{q})⟩u_{\alpha ,n_2}(\mathbf{k})$$

$$=$$

$$H_{\mathrm{int}}^{(1)}$$

$$=\frac{1}{2}\sum _{\mathbf{q}}\sum _{\alpha \beta }V(\mathbf{q}){\tilde{\rho }}_{\alpha }(\mathbf{q}){\tilde{\rho }}_{\beta }(-\mathbf{q})-\sum _{\alpha }\sum _{m,n\in M}\sum _{\mathbf{k}}{\tilde{\xi }}_{mn}^{\alpha }(\mathbf{k})c_{\alpha ;m}^{†}(\mathbf{k})c_{\alpha ;n}(\mathbf{k})$$

$${\tilde{\rho }}_{\alpha }(\mathbf{q})=\sum _{m_1,m_2\in M}\sum _{\mathbf{k}}{c}_{\alpha ,m_1}^{†}(\mathbf{k}+\mathbf{q})c_{\alpha ,m_2}(\mathbf{k})⟨u_{\alpha ,m_1}(\mathbf{k}+\mathbf{q})⟩u_{\alpha ,m_2}(\mathbf{k}$$

$${\tilde{\xi }}_{mn}^{\alpha }(\mathbf{k})$$

$$H_{\mathrm{int}}^{(2)}$$

$$⟨u_{\alpha ;m_1}({\mathbf{k}}_1+\mathbf{q})⟩u_{\alpha ;m_2}({\mathbf{k}}_1⟨u_{\beta ;n_1}({\mathbf{k}}_2-\mathbf{q})⟩u_{\beta ;n_2}({\mathbf{k}}_2)$$

$$H_{\mathrm{int}}^{(2)}=-\frac{1}{2}\sum _{\alpha }\sum _{\mathbf{k},\mathbf{q}}\sum _{m,n\in M}\sum _{m_2\in O}V(\mathbf{q})c_{\alpha ;m}^{†}(\mathbf{k})c_{\alpha ;n}(\mathbf{k})⟨u_{\alpha ;m}(\mathbf{k})⟩u_{\alpha ;m_2}(\mathbf{k}-\mathbf{q}⟨u_{\alpha ;m_2}(\mathbf{k}-\mathbf{q})⟩u_{\alpha ;n}(\mathbf{k})$$

$$H_{\mathrm{int}}^{(3)}$$

$$⟨u_{\alpha ,m_1}({\mathbf{k}}_1+\mathbf{q})⟩u_{\alpha ,m_2}({\mathbf{k}}_1⟨u_{\beta ,n_1}({\mathbf{k}}_2-\mathbf{q})⟩u_{\beta ,n_2}({\mathbf{k}}_2)$$

$$H_{\mathrm{int}}^{(3)}$$

$$H_{\mathrm{int}}^{(4)}$$

Full text

Following the remarkable discovery of correlated insulators and superconductivity in Magic Angle Twisted Bilayer Graphene (TBG) Cao et al. (2018a, b), a great deal of attention has been lavished on various “moiré materials". In TBG, the correlated insulator/superconductor has since been observed by other groups Yankowitz et al. (2019); Lu et al. (2019), leading to a wealth of new information. When the TBG is further aligned with a hexagonal Boron Nitride substrate (TBG/hBN), emergent ferromagnetism and a large Sharpe et al. (2019) or even quantized Serlin et al. (2019) anomalous Hall effect is observed at $3/4$ filling of the conduction band.

In other experiments, the moiré bands formed when ABC stacked trilayer graphene is aligned with one of the hBN (TLG/hBN) substrates display interesting strong correlation physicsChen et al. (2019a, , b). Applying a perpendicular displacement field $D$ enables control of the bandwidth leading to a gate-tunable correlated insulator at half-fillingChen et al. (2019a). Furthermore the sign of $D$ enables changing the band topologyZhang et al. (2019a); Chittari et al. (2019); Chen et al. (2019b): for one sign of $D$ the bands in each of the two valleys have equal and opposite Chern number while for the other sign of $D$, the Chern number in either valley is zero. Remarkably in the topologically non-trivial side, the $1/4$ filled state in the valence band shows ferromagnetism and a quantized anomalous Hall effect with Chern number 2Chen et al. (2019b). In yet other experiments, in twisted double bilayer graphene systems (TDBG, i.e bilayer graphene twisted relative to another bilayer graphene to a magic angle), good evidence is found for a spin polarized insulator Liu et al. ; Shen et al. (2019); Cao et al. at half-filling of the conduction band which gives way Liu et al. ; Shen et al. (2019), upon doping, to a superconductor which is likely also spin polarized. Theoretically, the conduction band in the TDBG system also has a non-zero Chern number which is equal and opposite for the two valleys Zhang et al. (2019a); Lee et al. (2019); Koshino (2019); Liu and Dai (2019).

In this paper, we study the physics of narrow moiré bands in the strong interaction limit, i.e when the Coulomb interaction is much larger than the bandwidth. Our focus is on the variety of moiré systems where the conduction and valence bands are separated from each other by energy gaps 111This is however not the case in the original TBG where there are Dirac points connecting conduction and valence bands, and our results do not directly apply to this system. We call $\nu_{T}$ the filling fraction including the spin and valley degrees of freedom, i.e the number of electrons per moiré unit cell in the active band. We restrict our attention to total integer fillings $\nu_{T}=1,2,3$, and study the nature of the correlated insulators in this strong interaction limit. Previous papers have presented physical arguments, and supporting Hartree-Fock calculations, that the natural fate of the system in this limit is a spin-valley ferromagnetic insulator Zhang et al. (2019a); Lee et al. (2019); Bultinck et al. (2019); Zhang and Senthil (2019); Zhang et al. (2019b); Xie and MacDonald (2018). However Hartree-Fock theory typically overestimates the stability of ferromagnetic states. Thus it is important to substantiate the physical arguments for spin-valley ferromagnetism through other less biased calculations. In the context of strained graphene, Ref. Ghaemi et al., 2012 provided numerical evidence for a valley-polarized insulator at fractional filling. Here we present analytical arguments and numerical calculations - both exact diagonalization (ED) and density matrix renormalization group (DMRG) - for strongly interacting nearly flat band models pertinent to many moiré materials. We show that indeed spin-valley ferromagnetic states are stabilized in the flat-band limit. Our results highlight the robustness of ferromagnetism in a narrow band even when it is topologically trivial Zhang and Senthil (2019). Moreover, we provide a quantitative estimate of the stability range of the ferromagnetic states in terms of the interaction to bandwidth ratio. We also show that intervalley coherent order is always disfavored compared to ferromagnetism in the limit where the flat band is a Landau level. We emphasize that the interactions renormalize the bare (non-interacting) dispersion of the band through various mechanisms detailed in the paper. The strong interaction limit is thus defined as the regime where residual interactions far exceed the renormalized bandwidth.

We consider a two dimensional material with spinful electrons occupying bands in each of two disconnected valleys (denoted $+$ and $-$) such that time-reversal maps one valley to the other. As a result, the $+$ and $-$ valleys of the active band have equal and opposite Chern number $C$. We consider density-density interactions of typical strength $U$, and focus on the strong coupling limit $U/W\gg 1$, where $W$ is the bandwidth of the active band.

Physical argument for ferromagnetism First consider the case of a topologically trivial band with non-zero Berry flux distribution. Such a situation arises, in TLG/h-BN for one sign of the displacement field $D$ Zhang et al. (2019a). Using a Wannier basis Brouder et al. (2007); Marzari et al. (2012), we can build an effective tight-binding model for the active band  Zhang and Senthil (2019). The Wannier functions have a finite extension which is needed to capture the Berry flux density in the band. Projecting the Coulomb interaction onto the Wannier basis leads to an on-site Hubbard interaction of order $U$ (as well as smaller terms between neighboring sites) but also to an inter-site ferromagnetic Hund’s interaction $J_{F}=g_{s}U$ Zhang and Senthil (2019). The coefficient $g_{s}$ depends on the overlap of Wannier densities at neighboring sites but stays finite even when $W\rightarrow 0$. In the large-$U$ limit, the ground state has a fixed number $\nu_{T}\in\mathbb{Z}$ of electrons at each site. The active degrees of freedom are local moments in the spin-valley space. In the strict limit $W\rightarrow 0$, the only coupling that survives between these local spin-valley degrees of freedom is $J_{F}$, giving rise to a spin-valley ferromagnet Zhang and Senthil (2019) (see also Ref. Kang and Vafek, 2019). As $W$ increases, there will also be antiferromagnetic inter-site superexchange $\frac{\alpha W^{2}}{U}$ with $\alpha$ a constant independent of $W$ and $U$. This antiferromagnetic exchange can dominate over the inter-site Hund’s exchange only when $W>\sqrt{\frac{g_{s}}{\alpha}}U$. Thus so long as $W/U$ is small enough we get a spin-valley ferromagnetic Mott insulator. In this mechanism the larger the extension of the Wannier functions, the larger the coefficient $g_{s}$, and thus the stronger the ferromagnetism.

Turning next to the case of topologically non-trivial $\pm$ Chern bands, symmetric Wannier functions cannot be localized Brouder et al. (2007); Marzari et al. (2012); we may view this as the limit of Wannier functions with infinite extension. Intuitively, the Hund’s effect will only be stronger than in the topologically trivial case, and hence a spin-valley ferromagnetic insulator is the likely ground state at all integer fillings.

Analytical considerations For a perfectly flat band separated by a large gap from other bands, the effective Hamiltonian takes the form of the Coulomb interaction projected onto the active band

[TABLE]

where $\delta\tilde{\rho}(\mathbf{q})=\tilde{\rho}(\mathbf{q})-\rho_{0}\delta^{2}(\mathbf{q)}$ is the deviation of the projected density from the average density $\rho_{0}$. The total projected density operator $\tilde{\rho}$ is summed over spin and valley indices and can be written

[TABLE]

Here $\tilde{\rho}_{\pm}$ are the projected densities in each valley $s$. $\sigma=\uparrow,\downarrow$ is the spin index. The $\lambda_{s}$ are valley dependent form factors which are defined in terms of the Bloch eigenstates $|u_{s,\mathbf{k}}\rangle$ through

[TABLE]

Because of the Berry flux distribution, $\lambda_{s}$ is a non-trivial function of $\mathbf{k}$ and $\mathbf{q}$. It is readily verified that $\tilde{\rho}(\mathbf{-q})=\tilde{\rho}(\mathbf{q})^{\dagger}$ as befits the total density operator. The Hamiltonian $H_{V}$ is invariant under a $U(2)\times U(2)$ rotation corresponding to independent charge and spin conservation within each valley. It is not $SU(4)$ invariant in spin-valley space due to the form factors. With the further assumption that the interaction $V(\mathbf{q})\geq 0$ for all $\mathbf{q}$ (satisfied for Coulomb and for short range repulsive potentials), it becomes clear that $H_{V}$ is positive semi-definite. Thus any state $\ket{\psi}$ that satisfies

[TABLE]

for all $\mathbf{q}$ is an exact ground state.

At $\nu_{T}$ integer, consider the state obtained by filling up $\nu_{T}$ bands to form an insulator. At $\nu_{T}=2$ this can be spin-polarized or valley-polarized. At $\nu_{T}=1,3$ this must be both spin and valley polarized. These states satisfy Eq. 5, hence they are exact eigenstates of $H_{V}$ with eigenvalue [math]. It follows that any ground state of $H_{V}$ must satisfy Eq. (5) for all $\mathbf{q}$. The remaining question is whether the spin-valley polarized states are the unique ground states. Indeed the same argument applied to a half-filled Hubbard model in the flat band limit (zero hopping) would yield ferromagnetic ground states; but these are degenerate with all other spin configurations and hence are not unique. In contrast in the flat Chern band, spin flips in the ferromagnetic state generically cost energy. This is well known in the idealized case of a band with uniform Berry curvature, i.e quantum Hall ferromagnetism Sondhi et al. (1993); Ezawa and Tsitsishvili (2009). In this case, the spin stiffness $\rho_{s}$ can be calculated exactly Zhang et al. (2019a) and is proportional to the square of the Chern number. This suggests that these ferromagnetic states may indeed be the unique ground states for any flat Chern band.

Besides spin-valley polarized states, inter-valley coherent (IVC) states are a plausible ground state candidate for moiré systems at integer filling $\nu_{T}$ Jung et al. (2015); Bultinck et al. (2019); Zhang et al. (2019b, a); Lee et al. (2019) (see also Ref. Po et al., 2018). To address these, it is instructive to consider a toy-model where the active band is the lowest Landau level of a system with opposite magnetic fields $\pm B$ for each valley. Then the projected density operators satisfy the Girvin-MacDonald-Platzman (GMP) algebra Girvin et al. (1986) ($l_{B}$ is the magnetic length):

[TABLE]

Since the commutator has the opposite sign for the two opposite valleys, the total density satisfies

[TABLE]

where $I^{z}(\mathbf{q})=\delta\tilde{\rho}_{+}(\mathbf{q})-\delta\tilde{\rho}_{-}(\mathbf{q})$ is the Fourier transform of the valley charge density. Thus in any ground state, by applying Eq. (5) to the left hand side of Eq. (7), we find that

[TABLE]

Thus the valley charge density cannot have fluctuations at any non-zero $\mathbf{q}$ in a ground state. If however there is IVC order then the breaking of valley $U_{v}(1)$ symmetry will lead to non-zero fluctuations of $I^{z}(\mathbf{q})$. For instance from the expected Goldstone fluctuations of the phase of the IVC order parameter the ground state correlator

[TABLE]

(Here $\kappa_{v}$ is the valley charge susceptibility,and $\rho_{sv}$ is the phase stiffness of the IVC order parameter). The exact result in Eq. (8) is clearly in conflict with this expectation. We conclude therefore that the IVC ordered state is not a ground state in this model. Moving away from the toy model above, introducing a small bandwidth and Berry curvature fluctuations will clearly not change this result 222In a generic Chern band, the GMP algebra is indeed approximated in the $\mathbf{q}\rightarrow 0$ limit Parameswaran et al. (2013); Repellin et al. (2014), and Eq. (8) still holds in this limit. We thus expect that, as suggested by Hartree-Fock, the IVC state is disfavored relative to the spin-valley polarized states in the nearly flat band limit for generic $\pm C$ Chern bands. In what follows we will provide numerical results supporting this expectation, and find that spin-valley ferromagnetism is indeed the ground state of nearly flat $\pm$ Chern bands at integer total filling.

Note that in realistic models, various effects are responsible for the dispersion of the active band. Besides the bare bandwidth $W_{\mathrm{bare}}$ of the non-interacting model, two bilinear terms of order $U$ contribute to the overall bandwidth. The first term is the Fock term stemming from the interaction between electrons in the active band and in the lower (fully occupied) bands. The second bilinear term is the difference between $H_{V}$ and its normal-ordered couterpart, and is proportional to the fluctuations of the squared form factor $|\lambda_{s}(\mathbf{k}+\mathbf{q},\mathbf{k})|^{2}$. See the supplementary material for more details. In our analytical considerations, we considered the ideal case where the renormalized bandwidth vanishes.

Exact diagonalization results for a single moiré band We consider the continuum momentum-space models Bistritzer and MacDonald (2011) of three iconic moiré systems (TBG/hBN Zhang et al. (2019b), TDBG Zhang et al. (2019a); Lee et al. (2019), TLG/hBN Zhang et al. (2019a)). Unless otherwise noted, for concreteness we respectively choose the twist angles $\theta=1.05^{\circ}$, $1.2^{\circ}$ and [math], and the displacement field $D=0$, $40\mathrm{mV}$ and $50\mathrm{mV}$. The active band is the valence (TLG/hBN) or the conduction (TBG/hBN and TDBG) band and has Chern number $C=\pm 1$ (TBG/hBN), $C=\pm 2$ (TDBG) and $C=\pm 3$ (TLG/hBN)333In TLG/hBN, a quantized anomalous Hall effect was in fact observed at $C=2$ Chen et al. (2019b). Here, we use the continuum non-interacting model which predicts a Chern $C=\pm 3$ band. A scenario involving interactions and supported by Hartree-Fock calculations was proposed in Ref. Chen et al., 2019b to explain this discrepancy.. We also consider the trivial ($C=0$) band obtained by switching the sign of the coupling to the top hBN layer in TBG/hBN or the sign of the displacement field $D$ in TLG/hBN. We take the limit where the active band is separated from other bands by a gap much larger than its bandwidth. The Hamiltonian is obtained by normal-ordering the projected Hamiltonian $H_{V}$ of Eq. (1), where the screened Coulomb interaction takes the form

[TABLE]

$r_{0}$ is the screening length (we choose $r_{0}=5.0$ in units of the moiré lattice constant). We also consider the addition of a kinetic term which gives the active band a width $W$.

[TABLE]

where $\mathbf{a_{1}},\mathbf{a_{2}}$ are the moiré lattice vectors, and we have taken a simplified dispersion (not the realistic one). We study this model using exact diagonalization at integer filling $\nu_{T}=1,2$ ($\nu_{T}=3$ is related to $\nu_{T}=1$ in our model through a particle-hole transformation). We call $N_{s}$ the number of moiré unit cells, and choose the aspect ratio of the finite cluster to be close to $1$.

We start by investigating the nature of the ground state in the limit $W_{\mathrm{bare}}=0$. In spite of the band flatness, the normal ordering of the interaction induces a finite dispersion (see supplementary materials), such that ferromagnetism is not guaranteed. Nevertheless, for $\nu_{T}=1$, we find that the ground state is always fully spin and valley polarized; the resulting state has a quantum anomalous Hall effect $\sigma_{xy}=Ce^{2}/h$. At $\nu_{T}=2$, maximal polarization of spin or valley leads to several correlated insulators all related by $U(2)\times U(2)$ symmetry in our model. For example, one with full spin polarization (but $I_{z}=0$), which is a valley-Hall insulator if $C\neq 0$, and one with full valley polarization ($S_{z}=0$), with anomalous Hall effect $\sigma_{xy}=2Ce^{2}/h$.444Including the weak inter-valley Hund’s interaction would break the $U(2)\times U(2)$ symmetry and lift the degeneracy between these two states. Numerically, we find that these $\nu_{T}=2$ polarized insulators indeed have the lowest energy in all models except for one: the $C=0$ TBG/hBN. In this case, the ground state is partially polarized, but the important finite-size effects prevent us from identifying its nature.

Adding a finite bandwidth $W$, we find that the ferromagnetic phase survives up to $W_{c}\simeq 0.05U$ in some models, and up to $W_{c}=0.2U$ in others at filling $\nu_{T}=1$. At $\nu_{T}=2$, the ferromagnet is relatively less stable, with critical values ranging from $W_{c}=0.04U$ to $W_{c}=0.1U$. The critical bandwidth is extracted from a finite-size extrapolation of exact diagonalization results which is detailed in Fig. 1.

The spin stiffness measures the energy change from twisting the spin boundary conditions

[TABLE]

To evaluate $\rho_{s}$, we calculated the spin-wave dispersion exactly for large systems (hundreds of moiré unit cells) by restricting the calculation to the valley-polarized, $S_{z}=N/2-1$ sector. Fig. 2a) shows this dispersion in the case of the $C=3$ model, and displays a remarkable data collapse for $36<N_{s}<256$. For these larger systems, calculating the spectrum of all spin and valley polarization sectors is impossible within exact diagonalization, but we rely on our previous result and assume that the other polarization sectors do not affect the low-energy properties. Fig. 2b) illustrates the influence of the Berry curvature distribution $F(\mathbf{k})$ on the spin stiffness ($F(\mathbf{k})$ is changed by tuning the microscopic parameters such as the twist angle $\theta$ or the displacement field $D$). While $F(\mathbf{k})$ does not uniquely determine the spin stiffness 555Even in the case of uniform Berry curvature, the real part of the quantum geometric tensor also enters the expression of $\rho_{s}$ Zhang et al. (2019a), for a given model (a given color in Fig. 2b), larger Berry curvature fluctuations appear to enhance $\rho_{s}$.

Ferromagnetism in the spin-valley Haldane model We now turn to DMRG calculations on infinitely long cylinders, which help us circumvent the size limitations of exact diagonalization. We use a tight-binding model, which facilitates DMRG 666While it is in principle possible to perform DMRG using a hybrid momentum-real space basis Motruk et al. (2016); Ehlers et al. (2017), it comes with the cost of longer-range interactions and more complex algorithms., but comes with an additional cost: due to the non-trivial Chern number of the active band, we must consider a two-band model. Additionally, the two-band model permits the consideration of the band gap energy scale, which we have supposed to be infinite until now. Our toy model is a tight-binding model on the honeycomb lattice based on the Haldane model Haldane (1988) with on-site Hubbard interaction of strength $U$.

[TABLE]

The hopping amplitudes $t_{ij}$ are non-zero for first ($t$) and second ($t_{2}$) neighbor and realize the Haldane model. We also consider third ($t_{3}$) and fourth ($t_{4}$) neighbor hopping. We tune these parameters to obtain a narrow conduction band with Chern number $C=\pm 1$ and $C=\pm 2$; we use the parameters $t_{2}/t=0.315e^{0.209\pi i}$ , $t_{3}=t_{4}=0$ for $C=\pm 1$ and $t_{2}/t=1.312i$, $t_{3}/t=1.312$, $t_{4}/t=0.524$.

We numerically obtained the ground state of $H_{1}$ using DMRG in the infinite cylinder geometry with $L_{y}=2$ unit cells along the perimeter. We extracted the spin and valley polarization of the ground state for several values of $U$. In spite of a large bond dimension dependence, we find that the onset of ferromagnetism $U_{c}/W$ is always smaller than the band gap (see Fig. 3c). Working in the valley-polarized or spin-polarized limits, we can better approach convergence in our DMRG calculations, and simulate wider cylinders (up to $L_{y}=3$). Our results are shown in Fig. 3 and confirm the ferromagnetic nature of the ground state for $U>U_{c}$ where $U_{c}\simeq 3W$ (the band gap is respectively $\Delta/W=6.0$ and $5.45$ for the $C=\pm 1$ and $C=\pm 2$ parameters). $U_{c}$ seemingly decreases with increasing system size $L_{y}$ (our complementary exact diagonalization results on this model show the same trend, see supplementary material for details), giving us confidence that $U_{c}<\Delta$ in the thermodynamic limit.

Discussion In this paper we have shown both analytical and numerical evidence for spin and valley polarization in nearly flat bands, which naturally emerge in several graphene moiré superlattices. Our results demonstrate a valley and spin polarized quantum anomalous Hall insulator at $\nu_{T}=1$ or $\nu_{T}=3$. At $\nu_{T}=2$, a spin-polarized valley Hall insulator and a valley-polarized quantum anomalous Hall insulator are both possible. Indeed recent experiments have already observed signatures for spin polarization in twisted double bilayer graphene and anomalous Hall effect in TBG/hBNSharpe et al. (2019) and TLG/hBNChen et al. (2019b).

The phenomenon we described is reminiscent of quantum Hall ferromagnetism, but there are important differences. Most naively, the finite bandwidth may destroy ferromagnetism, and we have quantitatively evaluated the position of this transition. A less expected difference is that flat-band ferromagnetism appears even when the Chern number is zero, due to the non-zero Berry flux. Deviating from the Landau level situation through large Berry curvature fluctuations has two opposite effects, which respectively destabilize and stabilize ferromagnetism: it may increase the strength of interaction-induced dispersive terms, but it may also enhance the spin stiffness. A natural future direction is to study the possibility of fractional quantum Hall effect from similar spontaneous time reversal breaking at fractional filling.

Acknowledgement We are thankful to Yin-Chen He and Kun Yang for insightful discussions. DMRG calculations were performed using the TeNPy Library (version 0.4.0)Hauschild and Pollmann (2018). This work was supported by NSF grant DMR-1911666, and partially through a Simons Investigator Award from the Simons Foundation to Senthil Todadri. C.R. is supported by the Marie Sklodowska-Curie program under EC Grant agreement No. 751859. Part of this work was performed during a visit of CR and TS at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.

Author contributions CR performed the ED simulations, ZD performed the DMRG simulations. YHZ provided the momentum space and lattice models. CR, YHZ and TS supervised the work of ZD. All authors contributed to the discussions and analysis of results.

I Renormalization of the bare bandwidth by interactions

In this section, we discuss the effective band broadening originating from the projection of the interaction onto the active band.

We first recall the description of the density interaction in the case of a periodic system without any band projection. We define $f(\mathbf{k})$ as operators for the microscopic electrons. $c(\mathbf{k})$ is the band operator. We have the following relation

[TABLE]

where $\alpha$ is the spin-valley index, $m$ is the band index and $I$ is the orbital index at the microscopic level (for example, in the Haldane model $I=A,B$ is the sublattice index). $u_{\alpha,m}(\mathbf{k})$ is the Bloch wavefunction of the band specified by $\alpha,m$. It is a vector for which $I$ is the index. Therefore there is a completeness relation:

[TABLE]

and an orthogonality relation

[TABLE]

In terms of $f$ the interaction term is written

[TABLE]

Substituting Eq. 14 we get

[TABLE]

We define the density operator

[TABLE]

With simple algebra, we can write $H_{\mathrm{int}}$ in the form $\rho(\mathbf{q})\rho(-\mathbf{q})$

[TABLE]

where we have used the completeness relation (15) for $\ket{u_{\alpha,m_{2}}(\mathbf{k-q})}$. Clearly the second term is just a chemical potential term. For a fixed density, the normal ordering does not matter.

Next we project the Hamiltonian to the active bands. The Hilbert space is a tensor product ${\cal H}=\prod_{m}{\cal H}_{m}$, where ${\cal H}_{m}$ is the Hilbert space of band m. Suppose $M$ is the set of active bands and $\bar{M}$ is the set of all the other (fully filled or fully empty) bands. We have ${\cal H}={\cal H}_{M}\otimes{\cal H}_{\bar{M}}$. We want to project to the subspace ${\cal H}_{S}={\cal H}_{M}\otimes\ket{\Psi_{\bar{M}}}$, where $\ket{\Psi_{\bar{M}}}$ is a product state in $H_{\bar{M}}$; $\ket{\Psi_{\bar{M}}}$ is the Slater determinant such that each momentum is either occupied or empty depending on the specific band.

We assume $\ket{\Psi_{\bar{M}}}$ is fixed and then we can trace it out to obtain an effective Hamiltonian in ${\cal H}_{M}$. The projection to ${\cal H}_{S}$ consists in calculating the expectation value of operators in ${\cal H}_{\bar{M}}$ for the state $\ket{\Psi_{\bar{M}}}$. The kinetic term is just a constant. We focus on the four fermion interaction in Eq. 18.

Let us consider $c^{\dagger}_{\alpha,m_{1}}(\mathbf{k_{1}+q})c^{\dagger}_{\beta,n_{1}}(\mathbf{k_{2}-q})c_{\beta,n_{2}}(\mathbf{k_{2}})c_{\alpha,m_{2}}(\mathbf{k_{1}})$. We can group the terms using the number of indices belonging to $M$. If there is an odd number of indices belonging to $M$, the term vanishes in ${\cal H}_{S}$. Therefore, we only need to consider two groups of terms: (1)$m_{1},m_{2},n_{1},n_{2}\in M$ (2)$m_{1},n_{2}\in M,m_{2},n_{1}\in\bar{M}$; (3)$m_{1},m_{2}\in M,n_{1},n_{2}\in\bar{M}$. (4) $n_{1},m_{2}\in M,m_{1},n_{2}\in\bar{M}$; (5)$n_{1},n_{2}\in M,m_{1},m_{2}\in\bar{M}$. Note that $m_{1},n_{1}\in M,m_{2},n_{2}\in\bar{M}$ vanishes.

Next we calculate the above five terms one by one. The first term gives:

[TABLE]

where

[TABLE]

[TABLE]

We find that if we do not use normal ordering, we need to add a bilinear term, which effectively behaves as a kinetic term proportional to the interaction strength.

The second term $H^{(2)}_{\mathrm{int}}$ is

[TABLE]

Projecting to $H_{S}$ is equivalent to substituting $c^{\dagger}_{\beta;n_{1}}(\mathbf{k_{2}-q})c_{\alpha;m_{2}}(\mathbf{k_{1}})$ with its expectation values under $\ket{\Psi_{\bar{M}}}$. Defining $O\subset\bar{M}$ as the set of occupied bands. It is easy to find that

[TABLE]

The term (3) is

[TABLE]

$\langle c^{\dagger}_{\beta;n_{1}}(\mathbf{k_{2}-q})c_{\beta;n_{2}}(\mathbf{k_{2}})\rangle=0$ for $\mathbf{q}\neq 0$. As a result

[TABLE]

where ${\cal N}_{O}$ is the number of occupied bands and $N_{s}$ is the number of $k$ points, and $V(0)=Ur_{0}$ for the screened Coulomb interaction Eq. (10) . This term is identical to the term (5). If there is only one active band ($m=n$), which is the case considered in this paper, then these terms are just a chemical potential which does not contribute to the bandwidth.

Term (4) is the same as the term (2) and gives

[TABLE]

One can easily calculate the expectation value and get

[TABLE]

which is equal to $H_{\mathrm{int}}^{(2)}$.

Summing all five terms $H_{\mathrm{int}}^{(i)}$, we get the final form of projected interaction

[TABLE]

or equivalently

[TABLE]

We have the additional bilinear terms $\tilde{\xi}^{\alpha}_{mn}(\mathbf{k})$ (defined in Eq. (23)) and

[TABLE]

Finally, we have the full Hamiltonian in two equivalent forms.

[TABLE]

or

[TABLE]

where we denote $H_{K}^{0}$ as the bare dispersion, and the effective kinetic terms are

[TABLE]

and

[TABLE]

In $H^{1}_{K}$, we add Hartree-Fock terms (actually only the Fock term) originating from occupied bands. For the graphene moiré terms, we expect that the bare band structure from the continuum model will be greatly renormalized by this effect. In $H^{2}_{K}$, we further include a dispersion term from the active band itself.

Using $H^{0}_{K}$, $H^{1}_{K}$ and $H^{2}_{K}$ we can define three bandwidths $W_{\mathrm{bare}}$ (the bare bandwidth), $W_{1}$ and $W_{2}$. In the main text, we have used $W_{2}=0$ for our analytical argument; in this limit, we have shown that the spin and valley polarized state is an exact ground state of the interaction. In the numerical exact diagonalizations, we have used $W_{1}=0$ as the flat-band limit; in this case there is no guarantee that the spin-valley polarized state is a ground state, and in fact, our results show examples where it is an excited state.

In a different perspective, we can calculate the effective dispersion at a certain filling $\nu_{T}$. In a quantum Hall ferromagnetic state, we can assume a flavor $\alpha$ is fully filled. The corresponding dispersion from Hartree-Fock calculation for the band $\alpha$ is:

[TABLE]

We can again define a bandwidth $W_{\mathrm{eff}}$ for the fully filled state. Note that there is an additional factor of 2. Therefore for the active band $W_{\mathrm{eff}}$ is again different from $W_{2}$. Compared to $W_{1}$, $W_{\mathrm{eff}}$ includes the Hatree-Fock term from the active band itself. However, $W_{2}$ only includes one half of this term. In that sense, it is surprising that $W_{2}=0$ is the condition for flat band ferromagnetism.

II Additional numerical results

In this section, we give additional numerical results on the two-band (Honeycomb) model $H_{1}$ Eq. (13) of the main text.

II.1 Exact diagonalization results

First, we complement our DMRG study with results of exact diagonalization (ED). We limit ourselves to the spin-polarized or valley-polarized versions of $H_{1}$, since only very few system sizes are accessible within the full (non-polarized) model using ED. Our results are summarized in Fig. 4. We find that the onset of ferromagnetism is the same in these two limits in finite size, and thus show only one plot. The results are consistent with the results of our DMRG study, with small differences which can be accounted for by invoking finite-size effects.

II.2 Bond dimension convergence of the DMRG results

Due to the Wannier obstruction of topological bands, it is necessary for a real space simulation to include both conduction and valence bands. This is highly challenging even in the parameter regimes where insulating ground states are expected, because a fully filled Chern band still carries large entanglement. More specifically, even though the insulating (fully-polarized) ground state is a product state in the (momentum-space) basis of Bloch wave functions, the entanglement entropy associated with a real-space cut is large. As a result, the critical value of the interaction $U_{c}$ still displays relatively large variations with bond dimension in our results (see Fig. 3c of the main text).

In the main text we have relied on energetics to detect the onset of ferromagnetism; we ran independent DMRG simulations for each spin and valley polarization sector, and compared all these energies to infer the spin and valley polarization of the ground state. Alternatively, we can use the fact that the correlation length associated with spin-spin correlations $\xi_{S\cdot S}$ diverges in the ferromagnetic state. Since $\xi_{S\cdot S}$ is very sensitive to the bond dimension, this method provides a more conservative estimate for $U_{c}/W$. The evolution of $\xi_{S\cdot S}$ with the interaction strength is shown in Fig. 5 for various values of the bond dimension. It shows us that finite bond dimension effects are still important at the bond dimensions we used. Nevertheless, even with this more conservative method, $U_{c}/W$ is still significantly smaller than the band gap.