A mechanism for anomalous Hall ferromagnetism in twisted bilayer graphene

TL;DR

This paper models the anomalous Hall effect in twisted bilayer graphene using a Landau level approach, revealing phase transitions and unique Chern band properties that differ from traditional quantum Hall ferromagnetism.

Contribution

It introduces a lowest Landau level model capturing the anomalous Hall effect and predicts phase transitions and excitonic vortex lattice formation in twisted bilayer graphene.

Findings

Phase transition from spin-valley polarized insulator to metal with increased bandwidth.

Chern band structure prevents inter-valley coherent phases.

Potential formation of excitonic vortex lattice at high interaction anisotropy.

Abstract

We use a lowest Landau level model to study the recent observation of an anomalous Hall effect in twisted bilayer graphene. This effective model is rooted in the occurrence of Chern bands which arise due to the coupling between the graphene device and its encapsulating substrate. Our model exhibits a phase transition from a spin-valley polarized insulator to a partial or fully valley unpolarized metal as the bandwidth is increased relative to the interaction strength, consistent with experimental observations. In sharp contrast to standard quantum Hall ferromagnetism, the Chern number structure of the flat bands precludes an instability to an inter-valley coherent phase, but allows for an excitonic vortex lattice at large interaction anisotropy.

Equations

$$c_{\pm ,\text{k}}^{†}=\frac{1}{\sqrt{N_x}}\sum _{n=0}^{N_x-1}{e}^{\pm i{k}_x(k_y+nQ)l_B^2}{c}_{\pm ,k_y+nQ}^{†},$$

$$n_{\pm }(\text{q})$$

$$H^i$$

$${\Delta }_{\mathbf{k}}={\Delta }_0\sum _{j=-\infty }^{\infty }{e}^{-i\frac{\pi }{2}{j}^2}{e}^{-\frac{1}{4}(2k_y+jQ{)}^2{l}_B^2-i{k}_x(2k_y+jQ)l_B^2}$$

$$\displaystyle e^{HF}(1,\Delta_{0})-e^{HF}(1,0)=\frac{1}{8h^{2}}\bigg{[}\int_{{\mathbf{k}},{\mathbf{q}}}u_{0}({\mathbf{q}})|\Delta_{+}-\Delta_{-}|^{2}$$

$$\displaystyle+\int_{{\mathbf{k}},{\mathbf{q}}}u_{1}({\mathbf{q}})|\Delta_{+}+\Delta_{-}|^{2}-4u_{1}(\mathbf{0})\int_{{\mathbf{k}}}|\Delta_{\mathbf{k}}|^{2}\bigg{]}\,,~{}~{}~{}$$

$$H(\text{k})=\sum _{{\text{g}}_1,{\text{g}}_2}\left(h^{++}(R(\theta /2)(\text{k}+\text{X}+{\text{g}}_1)){\delta }_{{\text{g}}_1,{\text{g}}_2}+h^{--}(R(-\theta /2)(\text{k}+\text{X}+{\text{g}}_1)){\delta }_{{\text{g}}_1,{\text{g}}_2}+\sum _{\tilde{\text{g}}}\left[T_{\tilde{\text{g}}}^{+-}{\delta }_{{\text{g}}_1,{\text{g}}_2+\tilde{\text{g}}}+T_{\tilde{\text{g}}}^{-+}{\delta }_{{\text{g}}_1+\tilde{\text{g}},{\text{g}}_2}\right]\right)$$

$$h(\text{k})=\left(\begin{array}{cc}0& e^{i\text{k}\cdot {\text{R}}_A}+e^{i\text{k}\cdot ({\text{R}}_A-{\text{R}}_1)}+e^{i\text{k}\cdot ({\text{R}}_A-{\text{R}}_2)}\\ e^{-i\text{k}\cdot {\text{R}}_A}+e^{-i\text{k}\cdot ({\text{R}}_A-{\text{R}}_1)}+e^{-i\text{k}\cdot ({\text{R}}_A-{\text{R}}_2)}& 0\end{array}\right),$$

$$T_{\text{0}}$$

$$T_{{\text{g}}_1}$$

$$T_{{\text{g}}_2}$$

$$H^{\tau }(\text{k})$$

$$c_{\text{k},\tau ,s}^{†}$$

$$=$$

$$\equiv$$

$${\psi }_{\text{r},s}^{†}$$

$$n_{[\text{r}]}(\text{q})$$

$$=$$

$$H^i$$

$$=$$

$$V_{\text{q}}([\text{r}]-[{\text{r}}^{\prime }])=\sum _{\text{t}}{e}^{i\text{q}\cdot \text{t}}V(\text{t}+[\text{r}]-[{\text{r}}^{\prime }]),$$

$$V_{{\tau }_1,{\tau }_1^{\prime },{\tau }_2,{\tau }_2^{\prime }}^{\text{k},{\text{k}}^{\prime }}(\text{q})=\sum _{[\text{r}],[{\text{r}}^{\prime }]}{\varphi }_{{\tau }_1,\text{k}}^{\ast }([\text{r}]){\varphi }_{{\tau }_1^{\prime },\text{k}+\text{q}}([\text{r}])e^{i({\tau }_2^{\prime }-{\tau }_2)\text{X}\cdot [\text{r}]}{e}^{-i\text{q}\cdot ([\text{r}]-[{\text{r}}^{\prime }])}{V}_{\text{q}}([\text{r}]-[{\text{r}}^{\prime }]){\varphi }_{{\tau }_2,{\text{k}}^{\prime }}^{\ast }([{\text{r}}^{\prime }]){\varphi }_{{\tau }_2^{\prime },{\text{k}}^{\prime }-\text{q}}([{\text{r}}^{\prime }])e^{i({\tau }_1^{\prime }-{\tau }_1)\text{X}\cdot [\text{r}{]}^{\prime }}$$

$$d_{x_0,k_y,\tau ,s}^{†}$$

$$=$$

$$\equiv$$

$$\Delta (\mathbf{r})=\sum _k{C}_k{e}^{iky}{e}^{-\frac{1}{2{\xi }^2}(x-k{\xi }^2{)}^2},\xi =\frac{l_B}{\sqrt{2}}$$

$${\mathcal{T}}_{\mathbf{R}}{\mathcal{T}}_{{\mathbf{R}}^{\prime }}=e^{iq\mathbf{B}\cdot (\mathbf{R}\times {\mathbf{R}}^{\prime })/\hslash }{\mathcal{T}}_{{\mathbf{R}}^{\prime }}{\mathcal{T}}_{\mathbf{R}}$$

$${\mathcal{T}}_1=e^{-ia{p}_y/\hslash };{\mathcal{T}}_2=e^{i\pi /2}{e}^{iay/l_B^2}{e}^{-ia(p_x+p_y)/(2\hslash )}$$

$$\mathcal{T}(a\hat{y})\Delta (\mathbf{r})=\sum _k{C}_k{e}^{-ika}{e}^{iky}{e}^{-\frac{1}{2{\xi }^2}(x-k{\xi }^2{)}^2}=\Delta (\mathbf{r})⟹k=k_j=\frac{2\pi j}{a}=jQ\text{ for }j\in \mathbb{Z}$$

$$\Delta (\mathbf{r})=\sum _{j=-\infty }^{\infty }{C}_j{e}^{i{k}_{j}y}{e}^{-\frac{1}{2{\xi }^2}(x-k_j{\xi }^2{)}^2}$$

$${\mathcal{T}}_2\Delta (\mathbf{r})$$

$$\Delta (\mathbf{r})=C_0\sum _{j=-\infty }^{\infty }{e}^{i\frac{\pi }{2}{j}^2}{e}^{i{k}_{j}y}{e}^{-\frac{1}{2{\xi }^2}(x-k_j{\xi }^2{)}^2}$$

$${\Delta }_{00}(\mathbf{k},{\mathbf{k}}^{\prime })=\int d\mathbf{r}\Delta (\mathbf{r}){\varphi }_{+,\mathbf{k}}^{\ast }(\mathbf{r}){\varphi }_{-,{\mathbf{k}}^{\prime },}(\mathbf{r})$$

$${\Delta }_{\mathbf{k}}={\Delta }_0\sum _{j=-\infty }^{\infty }{e}^{-i\frac{\pi }{2}{j}^2}{e}^{-\frac{1}{4}(2k_y+jQ{)}^2{l}_B^2}{e}^{-i{k}_x(2k_y+jQ)l_B^2}$$

Full text

Moiré graphene systems are a class of simple van der Waals heterostructures Geim and Grigorieva (2013) hosting interaction driven low-energy physics, making them an exciting platform to advance our understanding of correlated quantum matter. In twisted bilayer graphene (TBG) with a small twist angle between adjacent layers, interaction effects are enhanced by van Hove singularities coming from 8 nearly flat bands around charge neutrality (CN) in the Moiré- or mini-Brillouin zone (mBZ) Bistritzer and MacDonald (2011); Mele (2010); Lopes dos Santos et al. (2012, 2007); Trambly de Laissardière et al. (2010); Shallcross et al. (2010); Suárez Morell et al. (2010); Moon and Koshino (2012); Li et al. (2009); Luican et al. (2011); Yan et al. (2012); Brihuega et al. (2012); Ohta et al. (2012); Havener et al. (2014); Suárez Morell et al. (2010); de Gail et al. (2011); Uchida et al. (2014); Sboychakov et al. (2015); Jung et al. (2014); Wong et al. (2015); Fang and Kaxiras (2016). Observation of correlated insulating states when 2 or 6 of the 8 TBG flat bands are filled confirms the importance of interactions Cao et al. (2016); Kim et al. (2017); Cao et al. (2018); Yankowitz et al. (2018a); Kerelsky et al. (2018); Choi et al. (2019); Chen et al. (2018).

Recent experiments indicate that certain magic angle graphene devices have large resistance peaks at $\nu=0,3$, with the latter featuring an anomalous Hall (AH) effect detected via hysteresis in the Hall conductance as a function of the out-of-plane magnetic field Sharpe et al. (2019). The Hall conductance is of order $e^{2}/h$ but not yet quantized. Some have detected an meV-scale gap at CN, and a hysteretic behaviour of the Hall conductance with applied field at $\nu=-1$ Efetov (2019). In this work we discuss how the breaking of the 180-degree rotational symmetry ($C_{2z}$) by a partially aligned hexagonal boron-nitride (h-BN) substrate could explain these observations. A variety of works Jung et al. (2015); Hunt et al. (2013); Amet et al. (2013); Lee et al. (2016); Yankowitz et al. (2018b); Zibrov et al. (2018); Kim et al. (2018) have found that h-BN opens up a band gap at the Dirac points of graphene whose magnitude depends on the graphene / h-BN alignment angle, reaching $\Delta_{AB}\sim 17$meV Kim et al. (2018) to $\sim 30$meV Yankowitz et al. (2018b); Zibrov et al. (2018) at perfect alignment. Notably, even in seemingly unaligned devices with little or no observable h-BN induced Moiré potential, band gaps of several meV are still observed Kim et al. (2018); Zibrov et al. (2018). In TBG, the substrate can likewise gap out the band Dirac points at the $K_{\pm}$ points of the mBZ, splitting the bands as $8=4+4$ to create a gap at CN. We find that for certain sublattice splittings the resulting flat bands have Chern number $C=\pm 1$. This makes the TBG case similar to ABC stacked trilayer graphene, where under an appropriately directed electric field the flat bands have Chern numbers $\pm 3$ Zhang et al. (2018).

Accounting for the $C_{2z}$-breaking substrate, the basic structure of the problem is as follows. The gap at CN allows us to focus only on the four nearly degenerate conduction (valence) bands for fillings above (below) CN, i.e, $\nu>0$ ($\nu<0$). These four Chern bands are uniquely labeled by their valley $\tau=+,-$ and spin $s=\uparrow,\downarrow$; time-reversal switches the valley index and enforces opposite Chern numbers for bands from opposite valleys. Since a $|C|=1$ band is topologically equivalent to a Landau level (LL), the problem is roughly analogous to a spinful bilayer quantum Hall problem with one flux quanta per unit cell, but with opposite layers (valleys) experiencing opposite magnetic fields. The LLs are degenerate, but as in a quantum Hall ferromagnet (QHFM)Sondhi et al. (1993) at integer filling the electrons may open a gap by spontaneously polarizing into a subset of these LLs, or a coherent superposition of them. In conventional quantum Hall bilayers at filling $\nu=1$, interactions generically drive inter-layer coherence, e.g., the exciton condensate Eisenstein (2014); Tutuc et al. (2004). But the twist here is the opposing Chern numbers of the two valleys. We find that the Chern number structure provides a topological reason for penalizing a coherent state: an exciton condensate between $C=1,-1$ bands is analogous to a superconductor in a strong magnetic field, which forces vortices into the order parameter, reducing the gain in the correlation energy. Hence, a spontaneously valley-polarized (VP) state is stable and exhibits AH effect with Hall resistance $\sim h/e^{2}$ (QAH if completely spin and valley polarized). Further, pinning of valley-polarization by an out-of-plane $B^{z}$ due to a large orbital g-factor explains the presence of the $R_{xy}$ hysteresis loop observed in Ref. Sharpe et al. (2019).

The possibility of spin and valley polarization and/or quantum anomalous Hall physics and chiral edge states in TBG has been discussed previously in Refs. Zhang et al. (2018); Dodaro et al. (2018); Thomson et al. (2018); Kang and Vafek (2018); Ochi et al. (2018); Xie and MacDonald (2018); Xu and Balents (2018); Lin and Nandkishore (2019); Liu et al. (2018); Zhang and Senthil (2018), albeit from a different perspective. We also note that a recent self-consistent Hartree-Fock (HF) treatment of the continuum model exhibits spontaneous $C_{2z}T$ breaking at CN, though the resulting Chern numbers were $C=\pm 2$ Xie and MacDonald (2018).

Substrate-induced Dirac mass and Chern numbers– We model the effect of the h-BN substrate Jung et al. (2015) by including in our band calculations a uniform but $C_{2z}$ breaking A-B sublattice splitting $\Delta_{t}$ and $\Delta_{b}$ on the top and bottom layer respectively (see sup for details). While h-BN may also introduce a Moiré potential, its magnitude falls off much more rapidly with alignment angle than $\Delta_{t/b}$ Kim et al. (2018). For our calculations we used a twist angle $\theta\approx 1.05^{\circ}$, and have taken a phenomenological corrugation effect into account by using a larger AB/BA inter-layer hopping $w_{1}$ as compared to the AA/BB inter-layer hopping $w_{0}$. Taking $w_{0}/w_{1}=0.85$ results in flat bands separated from the dispersing bands by an energy gap of approximately $20$ meV (for zero sublattice splittings).

With sublattice splitting, the phases of the $\tau=+$ valley (or $K$-valley of monolayer graphene) Moiré Hamiltonian for different parameter regimes of $\Delta_{t}$ and $\Delta_{b}$ are shown in Fig. 1. We find four different regions where both Dirac cones in the mBZ are gapped because of the sublattice splittings. In these regions, there are two isolated flat bands. We find that these four regions have bands with Chern numbers Fukui et al. (2005) $C=\pm 1$ or $C=0$, and are separated from each other by a Dirac point at either the $K_{-}$ or $K_{+}$ point in the mBZ. In Fig. 1 we show the Chern number of the flat band for the $\tau=+$ valley above (below) CN in green (orange). The Chern number for the flat bands from the $\tau=-$ valley can be obtained by time-reversal.

The location of the $C=\pm 1$ phases can be understood from the fact that for small $\Delta_{t}=\Delta_{b}>0$ or $\Delta_{t}=\Delta_{b}<0$, the leading order effect of the sublattice potentials is to generate Dirac masses with the same sign at both the $K_{-}$ and $K_{+}$ points of the mBZ. Because both Dirac cones in a single valley have the same chirality, this leads to bands with Chern number $\pm 1$, a feature earlier work dubbed a “flipped Haldane model”Zou et al. (2018) (see also Po et al. (2018, 2018); Song et al. (2018)). From Fig. 1 we see that even if only one of the layers has a non-zero sublattice splitting, the strong inter-layer coupling ensures that both Dirac cones at the mBZ $K$-points acquire a mass.These findings can also be inferred analytically within the “chiral” approximation of tBLG San-Jose et al. (2012); Tarnopolsky et al. (2018), in which all bands are sub-lattice polarized and carry Chern number $C=\sigma\tau$, where $\sigma$ denotes sublattice.

Metal - valley polarization competition– In this work, we focus only on the four flat conduction bands above the CNP (the highlighted band in Fig. 1 and its valley and spin counterparts). In the supplement, we numerically justify this for TBG, showing that $\Delta_{t}\sim 15$ meV ($\Delta_{b}=0$) creates a $30$ meV gap between valence and conduction bands sup . To phenomenologically model the effect of interactions in this set of bands we adopt a lowest Landau level (LLL) description. We can map the Chern bands to a LLL by constructing the Wannier-Qi states Marzari et al. (2012); Qi (2011); sup . In the following, we use an approximation where the Wannier-Qi states of the flat bands are replaced by the continuum LLL wave functions of a two-dimensional electron gas. Physically, this amounts to neglecting the inhomogeneous Berry curvature in the Chern bands. The AH effect and edge transport reported in Ref. Sharpe et al., 2019 can be explained if there is one VP hole per Moiré unit cell. From the data in Ref. Sharpe et al. (2019) is not possible to exclude a spin-unpolarized, gapless phase. If the spins do polarize however, the underlying mechanism is expected to be the same as in conventional QHFM Sondhi et al. (1993), and is not sensitive to the opposite Chern numbers of the two valleys. Therefore, in the analysis below we ignore spin and focus on the mechanism of valley polarization. Considering the uniform repulsive nature of the projected Coulomb interaction and the numerical evidence against stripes in the LLL Shibata and Yoshioka (2004), we disregard the possibility of interaction-induced charge density waves, and focus on the competition between valley-polarized, inter-valley coherent and metallic phases. For this we need to introduce two parameters in our LLL toy model: the bandwidth and the interaction anisotropy. To achieve a non-zero bandwidth we use a square lattice potential, that sidesteps the complexities of a hexagonal lattice and allows analytical progress.

We consider a torus of length $L_{x}$ ($L_{y}$) in the $x$ ($y$) direction, with a magnetic field perpendicular to the surface. We choose units in which $L_{x}L_{y}=2\pi N_{\phi}l_{B}^{2}\equiv N_{\phi}a^{2}$, where $N_{\phi}$ is the number of flux quanta piercing the torus, and $l_{B}=(\hbar/eB)^{-1/2}$ is the magnetic length. In particular, we will take $L_{x}=N_{x}a$ and $L_{y}=N_{y}a$, with $N_{\phi}=N_{x}N_{y}$. Next to the magnetic field, we also add a periodic potential $V_{P}(x,y)=w(\cos(2\pi x/a)+\cos(2\pi y/a))$, such that there is exactly $2\pi$ flux in each unit cell. The potential is invariant under translations over $a$ in both the $x$ and $y$-direction, which means that the momenta $k_{x}=n\frac{2\pi}{N_{x}a}$ and $k_{y}=n\frac{2\pi}{N_{y}a}$ ($n\in\mathbb{Z}$) are good quantum numbers.

We are interested in the physics in the LLL with Chern numbers $C=1,-1$. The electron creation operator projected in these subspaces takes the form $\psi^{\dagger}_{\pm}(x,y)=\frac{1}{\sqrt{L_{y}l_{B}\sqrt{\pi}}}\sum_{k}e^{iky-\frac{1}{2l_{B}^{2}}(x\mp kl_{B}^{2})^{2}}c_{\pm,k}^{\dagger}$, where we have chosen the Landau gauge which explicitly preserves (continuous) translation symmetry in the $y$-direction, such that $k=2\pi n/L_{y}=2\pi n/N_{y}a$ with $n\in\{0,1,\dots,N_{x}N_{y}\}$. We now proceed in analogy to Ref. Mishmash et al., 2018, and define the Bloch states $c^{\dagger}_{\pm,(k_{x},k_{y})}=c^{\dagger}_{\pm,\textbf{k}}$ as

[TABLE]

where $Q=\sqrt{2\pi}/l_{B}=2\pi/a$. The density operator in the LLL $n_{\pm}(\textbf{q})=\int\mathrm{d}\textbf{r}\,e^{-i\textbf{q}\cdot\textbf{r}}\psi^{\dagger}_{\pm}(\textbf{r})\psi_{\pm}(\textbf{r})$ takes the form

[TABLE]

where the form factor is given by $F(\textbf{q})=e^{-\textbf{q}^{2}l_{B}^{2}/4}$. In the Bloch basis, the Hamiltonian term associated with the periodic potential takes the diagonal form $H^{p}=\sum_{\textbf{k}}\varepsilon_{\textbf{k}}(c^{\dagger}_{+,\textbf{k}}c_{+,\textbf{k}}+c^{\dagger}_{-,\textbf{k}}c_{-,\textbf{k}})$, with $\varepsilon_{\textbf{k}}=-we^{-\pi/2}[\cos(k_{x}a)+\cos(k_{y}a)]$.

We are interested in the effect of density-density interactions on the LLL electrons moving in the periodic potential, described by the following Hamiltonian:

[TABLE]

where we neglect the small inter-valley scattering terms sup . We will consider a general repulsive interaction of the form $V({\mathbf{q}})F^{2}({\mathbf{q}})=u_{0}({\mathbf{q}})(\mathds{1}+\tau^{x})+u_{1}({\mathbf{q}})(\mathds{1}-\tau^{x})$. In analogy to quantum Hall ferromagnetism Sondhi et al. (1993); Eisenstein (2014); Ezawa and Tsitsishvili (2009) and related strongly coupled systems Neupert et al. (2012a, 2011), at half-filling of the two bands we expect that the main effect of $H^{i}$ is to introduce a valley Hund’s coupling between the electrons resulting in an insulating ground state. On the other hand, the kinetic term $H^{p}$ coming from the periodic potential favors a metal over the VP insulator. To study the competition between these two phases, we perform a HF analysis using Slater determinants with correlation matrix $\langle c^{\dagger}_{\tau,\textbf{k}}c_{\tau^{\prime},\textbf{k}^{\prime}}\rangle=\delta_{\tau,\tau^{\prime}}\delta_{\textbf{k},\textbf{k}^{\prime}}\Theta(\epsilon^{\tau}_{F}-\epsilon_{\textbf{k}})$, such that $\sum_{\tau}\sum_{\textbf{k}}\Theta(\epsilon^{\tau}_{F}-\epsilon_{\textbf{k}})=N_{\phi}$. The possibility of inter-valley coherent states is addressed in the next section. For each Slater determinant, we define the corresponding valley polarization $P_{v}$ as $P_{v}=(N_{+}-N_{-})/N_{\phi}$, where $N_{+}$ ($N_{-}$) is the number of electrons in the $+$ ($-$) valley. Without loss of generality, we restrict to $P_{v}>0$.

We first consider an isotropic ($u_{1}(\textbf{q})=0$) dual-gate screened Coulomb potential with LLL form factors $u_{0}(\textbf{q})=2\pi Ue^{-\textbf{q}^{2}l_{B}^{2}/2}\tanh{(d|\textbf{q}|)}/|\textbf{q}|$, and screening length $d=a$. Using this interaction potential, we calculated the HF energy $E^{HF}$ sup . We find that for $W/U\lesssim 0.6$, where $W\equiv 4we^{-\pi/2}$ is the bandwidth, the completely VP state indeed has the lowest energy. When $W/U\approx 0.6$, the valley polarization $P_{v}$ of the optimal Slater determinant jumps and starts decreasing continuously, indicating a first-order Mott transition from the VP insulator to an itinerant valley-ferromagnet. Around $W/U\approx 2.0$, $P_{v}$ continuously goes to zero and a conventional metallic phase sets in sup .

Inter-valley coherence and exciton vortex lattice– In bilayer QH ferromagnets, the insulating layer-polarized state is unstable to a uniform exciton condensate or inter-layer coherent state in presence of infinitesimal interaction anisotropy $u_{1}(\textbf{q})>0$ Eisenstein (2014). The situation here is different as even with $u_{1}(\textbf{q})=0$, there is no SU$(2)$ valley symmetry because of the Chern number mismatch. The VP state therefore only breaks discrete symmetries, indicating there will be no instability of this insulating state. Another, more physical, way to understand the absence of an exciton condensation instability is to use an analogy with type II superconductors. Because electrons in bands with an opposite Chern numbers effectively see opposite magnetic fields, an electron-hole condensate $\Delta({\mathbf{r}})=\langle c^{\dagger}_{+,{\mathbf{r}}}c_{-,{\mathbf{r}}}\rangle$ will behave like a charge 2e superconducting order parameter in a perpendicular magnetic field. However, in our scenario a Meissner-like effect, corresponding to uniform amplitude of the exciton order parameter, is ruled out from the outset. Rather, the magnetic field must leak through vortices in the exciton order parameter, leading to an excitonic vortex lattice phase. In this section, we show that both the VP insulator and the unpolarized metal are energetically favorable to the exciton vortex lattice, for sufficiently small interaction anisotropy $u_{1}(\textbf{q})$.

For our LLL model, we can derive an exact expression for the exciton vortex lattice order parameter $\Delta(\textbf{r})$. To respect all symmetries of the square lattice, we expect $\Delta(\textbf{r})$ to have vortices at both the lattice sites and the plaquette centers, leading to a $4\pi$ vorticity in each unit cell. In the analytically tractable limit, we can uniquely determine $\Delta({\mathbf{r}})$ up to a translation by demanding its invariance under the magnetic translations ${\mathcal{T}}(a\hat{x})$ and ${\mathcal{T}}\left(\frac{a}{2}(\hat{x}+\hat{y})\right)$, connecting the anticipated vortices sup . In Fig. 2 we plot the magnitude of $\Delta(\textbf{r})$ thus obtained, from which we clearly see the expected Abrikosov vortex lattice. Projecting $\Delta({\mathbf{r}})$ to the LLL Bloch basis wavefunctions $\phi_{\pm,{\mathbf{k}}}({\mathbf{r}})$ leads to a diagonal order parameter

[TABLE]

where $\Delta_{0}$ represents the overall strength of the exciton condensate. $\Delta_{\mathbf{k}}$ has two nodes with identical phase winding at ${\mathbf{k}}=\pm(\pi/2,-\pi/2)$, as shown in Fig. 2 sup .

The presence of two zeros in the BZ with the same phase winding is a topological requirement for the exciton order parameter, and is not an artifact of our effective LLL model. In an isolated band $a$ with non-zero Chern number $C_{a}$, the phase of the electron creation operator $c^{\dagger}_{a,{\mathbf{k}}}$ cannot be chosen to be both continuous and single-valued over the BZ. In particular, it must wind $2\pi C_{a}$ times along the boundary of the BZ in a continuous gauge choice. This implies that the phase of $\Delta_{\mathbf{k}}=\langle c^{\dagger}_{+,{\mathbf{k}}}c_{-,{\mathbf{k}}}\rangle$ winds $2\pi(C_{a}-C_{b})=4\pi$ times along the BZ boundary for bands from opposite valleys with $C_{a}=1$ and $C_{b}=-1$, which precisely corresponds to winding around two zeros with identical chirality.

We now demonstrate that variational states with an exciton vortex lattice have higher energy than the VP state or the metal for small anisotropy $u_{1}$ in the interaction $H^{i}$. We consider the Slater determinant ground state $\ket{\psi_{MF}}$ of the mean-field Hamiltonian $H_{MF}=\sum_{{\mathbf{k}},\tau,\tau^{\prime}}c^{\dagger}_{{\mathbf{k}},\tau}h_{\tau,\tau^{\prime}}({\mathbf{k}})c_{{\mathbf{k}},\tau^{\prime}}$, where $h_{\tau,\tau^{\prime}}(\textbf{k})=\epsilon_{{\mathbf{k}}}\mathds{1}+h\tau^{z}+\text{Re}(\Delta_{\textbf{k}})\tau^{x}+\text{Im}(\Delta_{\textbf{k}})\tau^{y}$. $\ket{\psi_{MF}}$ is characterized by the valley polarization $P_{v}$ (determined by $h$) and an exciton vortex lattice of strength $\Delta_{0}$, to be treated as variational parameters. The correlation matrix evaluated in this state takes the form of the projector $\langle c^{\dagger}_{\tau,{\mathbf{k}}}c_{\tau^{\prime}{\mathbf{k}}^{\prime}}\rangle=P_{\tau,\tau^{\prime}}({\mathbf{k}})\delta_{{\mathbf{k}},{\mathbf{k}}^{\prime}}$, which can be used to evaluate the regularized HF energy density $e^{HF}(P_{v},\Delta_{0})$ of the variational state for a given microscopic interaction at a fixed filling $\nu=1$. We find that the global minimum of $e^{HF}$ lies at $|P_{v}|=1$ and $\Delta_{0}=0$ for the insulator in the limit of flat bands and isotropic interaction ($u_{1}=0$) sup . We next show that the states of interest, with a fixed valley polarization $P_{v}$ at filling $\nu=1$, are stable to the formation of an vortex lattice in presence of small interaction anisotropy. To do this, we consider the difference in energy density $e^{HF}(P_{v},\Delta_{0})-e^{HF}(P_{v},0)$ perturbatively in $|\Delta_{0}|$ for arbitrary repulsive interaction parametrized by $u_{0}$ and $u_{1}$; a positive difference would indicate that $\Delta_{0}=0$ corresponds to an energy minimum. For the polarized phase, we find

[TABLE]

where $\Delta_{\pm}\equiv\Delta_{{\mathbf{k}}\pm{\mathbf{q}}/2}$ sup . For a uniform exciton condensate, $\Delta_{\mathbf{k}}=\Delta_{0}$ and this energy difference is negative sup . However, for an exciton order parameter formed with electrons and holes from opposite Chern bands, $\nabla_{\mathbf{k}}\Delta_{\mathbf{k}}\neq 0$. Therefore, when $u_{1}$ is sufficiently small compared to $u_{0}$ the energy of the state with non-zero $\Delta_{\mathbf{k}}$ is higher. So the VP state with $\Delta_{0}=0$, previously shown to be the ground state with an isotropic interaction for small $W/u_{0}$, is indeed robust to small interaction anisotropy. Analogous computations sup show that the unpolarized metal ($P_{v}=0=\Delta_{0}$) is stable to the vortex lattice as well. An approximate phase diagram of our model for a short-range (LLL-projected) interaction anisotropy $u_{1}({\mathbf{q}})=u_{1}e^{-\textbf{q}^{2}l_{B}^{2}/2}$ is presented in Fig. 3. For TBG, we expect $W/U\lesssim 0.2$ from the ratio of the bandwidth to the Coulomb interaction, and the anisotropy $u_{1}/U\lesssim 0.01$ to be small sup ; Chatterjee et al. (2019), indicating a VP phase consistent with experiments Sharpe et al. (2019); Serlin et al. (2019). In the supplement, we numerically solve the mean-field equations for TBG on hBN at $\nu=3$ and confirm that the spin and VP QAH state is indeed the ground state.

Valley Zeeman effect– Having argued in favor of a VP state at $\nu=3$, we turn to the observed hysteresis in the $\nu=3$ Hall conductance as a function of out-of-plane magnetic field $B^{z}$  Sharpe et al. (2019). To this end, we compute the orbital $g_{v}$-factor for the TBG conduction bands. In a band $\tau$ without time-reversal electrons can carry a momentum-dependent orbital moment $m_{\tau,\textbf{k}}$ Xiao et al. (2010); Chang and Niu (1996). Time-reversal ensures that $m_{\tau,\textbf{k}}=-m_{-\tau,-\textbf{k}}$, which averaged over the mBZ produces a valley-Zeeman splitting $E=-g_{v}\frac{\tau^{z}}{2}\mu_{B}B^{z}$. We find that for $\Delta_{b}=0,\Delta_{t}\sim 10-30$ meV, $g_{v}$ ranges from approximately -2 to -6 sup . Note that for $B^{z}>0$, the $C=1$ band comes down in energy. The sign of this effect is in agreement with the Landau fans of Refs. Sharpe et al. (2019); Serlin et al. (2019).

Conclusion– We showed that broken inversion symmetry in TBG due to substrate (h-BN) coupling leads to two Chern bands per valley. Spontaneous polarization of holes in spin and valley space then leads to an AH state at $\nu=3$. Using a LLL model, a HF analysis establishes a stable VP state as the ground state when the bandwidth is small compared to the interaction strength. The opposite Chern numbers for the two valleys precludes uniform inter-valley coherence. The resultant exciton vortex lattice structure reduces correlation energy gain and stabilizes valley-polarization. This result agrees with numerical work on a Hubbard model Neupert et al. (2012b).

Note added– Recently, a quantized AHE with net Chern number $C=1$ has been observed for a gapped insulator at $\nu=3$ in TBG aligned with h-BN Serlin et al. (2019), consistent with our theoretical results. Quantized AHE arising from valley-Chern bands have also been observed Chen et al. (2019) and proposed Liu et al. (2019); Lee et al. (2019) in other Moiré heterostructures, in accordance with our phenomenological picture of interaction in nearly flat bands with opposite Chern numbers.

Acknowledgements– We thank Aaron Sharpe, Eli Fox and David Goldhaber-Gordon for discussions about their data and sharing their insights. We also thank Ryan Mishmash for explaining the details of Ref. Mishmash et al., 2018 to two of us (SC and NB), and Ehud Altman, Senthil Todadri and Andrea Young for inspiring discussions. Our work overlaps with concurrent work by Y. Zhang, D. Mao and T. Senthil Zhang et al. (2019). MZ and NB were supported by the DOE, office of Basic Energy Sciences under contract no. DE-AC02-05-CH11231. SC acknowledges support from the ERC synergy grant UQUAM via E. Altman.

I Flat bands with sublattice splitting

As was shown in Ref. Jung et al. (2015); Hunt et al. (2013); Amet et al. (2013), the h-BN substrate generates a substantial Dirac mass when it is nearly aligned with the graphene sheet. We model this by introducing a $C_{2z}$ symmetry breaking sublattice-staggered potential $\Delta_{t}$ and $\Delta_{b}$ in respectively the top and bottom layer graphene sheet.

For the Moiré Hamiltonian, we consider a commensurate Moiré pattern, obtained from an AA stacked bilayer where the top and bottom graphene layers are rotated relative to each other along an out-of-plane rotation axis centered at a hexagon by an angle $\theta$. This gives a Moiré super lattice with microscopic $C_{6z}$ symmetry, which is found to be a very good approximate low-energy symmetry even for microscopically less symmetric Moiré patterns obtained from different initial stacking alignments or different rotation axis Zou et al. (2018). We choose to work with a commensurate pattern in order to use sharply defined Moiré lattice and reciprocal lattice vectors. However, the relevant properties of the electronic band structure around charge neutrality do not rely on the assumption of commensurability. In figure 4(a) we show the mono-layer graphene Brillouin zone with our convention for the reciprocal lattice basis vectors and the high symmetry points $K_{+}$ and $K_{-}$.

We now consider following spinless (spin-orbit coupling is negligible) single-valley Moiré Hamiltonian

[TABLE]

Here, $\textbf{g}_{1}$ and $\textbf{g}_{2}$ lie on the Moiré reciprocal lattice, $R(\theta)$ is a rotation matrix over angle $\theta$, $h^{++}(\textbf{k})=th(\textbf{k})+\Delta_{t}\sigma^{z}$ ($h^{--}(\textbf{k})=th(\textbf{k})+\Delta_{b}\sigma^{z}$) is the mono-layer graphene Hamiltonian of the top (bottom) layer with hopping strength $t=2.61$ eV and a sublattice splitting $\Delta_{t}\sigma^{z}$ ($\Delta_{b}\sigma^{z}$). The mono-layer graphene Hamiltonian is given by

[TABLE]

where $\textbf{R}_{1},\textbf{R}_{2}$ are the graphene Bravais lattice vectors and $\textbf{R}_{A},\textbf{R}_{B}$ are the sublattice vectors. X is the position of the center of the mBZ at the mono-layer $K_{+}$-points as shown in Fig.4(b). In the commensurate case we are considering here, X lies on the Moiré reciprocal lattice. The inter-layer coupling is given by the matrices

[TABLE]

where $\omega=e^{i2\pi/3}$, $\textbf{g}_{1}=(R(\theta/2)-R(-\theta/2))\textbf{G}_{1}$ and $\textbf{g}_{2}=(R(\theta/2)-R(-\theta/2))\textbf{G}_{2}$. The AB inter-layer hopping strength is $w_{1}=98$ meV. To phenomenologically incorporate the corrugation of the bilayer system we have used an AA-AB inter-layer hopping ratio $w_{0}/w_{1}=0.85$.

In Fig. 5 we show the resulting flat bands in the mBZ around charge neutrality of the single-valley Moiré Hamiltonian along high-symmetry paths, for different strengths sublattice splittings $\Delta_{t}$ and $\Delta_{b}$. The twist angle in these calculations was $\theta\approx 1.05^{\circ}$. When $\Delta_{t}=\Delta_{b}=0$, the flat bands have Dirac cones at $K_{+}$ and $K_{-}$ and are separated from the dispersive bands by an energy gap of approximately $20$ meV. If one of the sublattice splittings is non-zero, both Dirac cones acquire a mass because of the strong inter-layer coupling. In Fig. 5 we also show different plots with $\Delta_{t}=15$ meV constant and decreasing negative $\Delta_{b}$ to show the two Chern number changing transitions where a Dirac cone occurs at either $K_{+}$ or $K_{-}$.

II Suppression of inter-valley scattering

We write the single valley Moiré Hamiltonian schematically as

[TABLE]

where again the vectors $\textbf{g}_{i}$ lie on the reciprocal lattice of the Moiré super lattice. Here, we introduced the notation that $\tau\in\{+,-\}$ represents the different Dirac valleys, located at the high symmetry $K-$points of the mono-layer graphene BZ. The sublattice degree of freedom is denoted by $\sigma\in\{A,B\}$, and the two graphene layers are labeled by $\xi\in\{+,-\}$. The carbon atoms are located at positions r, such that r is of the form $\textbf{r}=R(\xi\theta/2)(m\textbf{R}_{1}+n\textbf{R}_{2}+\textbf{R}_{\sigma})$, where $m,n\in\mathbb{Z}$, $\textbf{R}_{1},\textbf{R}_{2}$ are the graphene Bravais lattice vectors and $\textbf{R}_{A},\textbf{R}_{B}$ the sublattice vectors.

Importantly, for $H^{\tau}(\textbf{k})$ we define the momentum k relative to the center of the mini-Brillouin zone located at the mono-layer $K_{\tau}$-points of top and bottom layer. In the second line we diagonalized the Moiré Hamiltonian using the unitary matrices $U$. Because we are interested in one band per valley, we drop the $\mu$ band index and associate $\tau$ with the band label. In this notation, we write the flat band states in each valley as

[TABLE]

Here we have introduced the notation $\textbf{r}=[\textbf{r}]+\textbf{t}$, where $[\textbf{r}]$ is the part of r lying in the Moire unit cell centered at the origin and t are Moire lattice vectors. To clarify the notation in going from Eq. (12) to Eq. (13), recall that the position label r contains the information contained in the labels $\xi$ and $\sigma$ via the relation $\textbf{r}=R(\xi\theta/2)(m\textbf{R}_{1}+n\textbf{R}_{2}+\textbf{R}_{\sigma})$. The index $s$ refers to spin. Inverting the above expression now gives us the electron raising operator at position r projected in the two flat bands:

[TABLE]

With the projected electron raising operator, which we simply rewrite as $\psi^{\dagger}_{\textbf{t},[\textbf{r}],s}$, we now define

[TABLE]

Using the above expression we can write the microscopic Coulomb interaction projected into the flat bands as

[TABLE]

where

[TABLE]

and $V(\textbf{r}-\textbf{r}^{\prime})$ is the microscopic Coulomb interaction. The projected interaction coefficients are given by

[TABLE]

Now it is important to remember that $\phi_{\tau,\textbf{k}}(\textbf{r})=\sum_{\textbf{g}}U_{\tau,\textbf{k}}(\xi,\sigma,\textbf{g})e^{i\textbf{g}\cdot\textbf{r}}$. Because for the flat band states $U_{\tau,\textbf{k}}(\xi,\sigma,\textbf{g})$ decays fast with $|\textbf{g}|$, $\phi_{\tau,\textbf{k}}(\textbf{r})$ varies slowly within the Moiré unit cell. So if $V(\textbf{r}-\textbf{r}^{\prime})$ is sufficiently long-range (like Coulomb), then the sums over $[\textbf{r}]$ and $[\textbf{r}^{\prime}]$ will suppress the terms with $\tau_{1}\neq\tau^{\prime}_{1}$ and $\tau_{2}\neq\tau^{\prime}_{2}$. For this reason, we restrict to the dominant density-density terms in our effective Landau-level problem.

III Construction of Wannier-Qi states

Because the single-valley flat bands split by the one-sided staggered potential have Chern number $\pm 1$, one cannot construct exponentially localized Wannier functions for these bands Thouless (1984). However, using the right gauge choice it is possible to construct Wannier functions that are exponentially localized along one direction. Using these quasi-one dimensional Wannier states there exists a natural mapping from the lattice system to a Landau-level system, as pointed out by Qi in Ref.Qi (2011). Here we review this mapping in the context of TBG.

Consider a system with periodic boundary conditions along two directions, which we refer to as the $x$ and $y$-directions. Using the flat band states as defined in the previous appendix, we construct the superlattice Wannier-Qi functions as follows:

[TABLE]

Here, $e^{i\alpha_{\tau}(\textbf{k})}$ ensures an optimal gauge choice such that the functions $W_{\tau,x_{0},k_{y}}(\textbf{r})$ are exponentially localized in the $x$-direction around the lattice position $x_{0}$. We now imagine adiabatically threading $2\pi$ flux through the hole of the torus, such that the flux is felt by a particle moving on closed path in the $y$-direction. Because of the Chern number $|C|=1$, this adiabatic process will change the polarization in the $x$-direction by one ‘polarization quantum’ Resta (1992); King-Smith and Vanderbilt (1993), which means that the centers of the Wannier functions all shift by one Moire lattice vector along the $x$-axis (the direction in which they shift depends on the sign of the Chern number). This implies that $W_{\tau,x_{0},k_{y}+g}(\textbf{r})=W_{\tau,x_{0}+\tau t_{x},k_{y}}(\textbf{r})$, where $g$ is the norm of the Moiré reciprocal basis vectors. Therefore, we can use $k\equiv k_{y}+\tau g$ as a single label for our Wannier-Qi states $W_{\tau,k}(\textbf{k})$ (for each $k_{y}\in[0,g]$, there is one Wannier function with a particular value for $x_{0}$ in each Moiré unit cell). One can now straightforwardly map the Chern band to a LLL, by replacing each Wannier-Qi function $W_{\tau,k}(\textbf{r})$ by the corresponding LLL Gaussian wave function. One of the main approximations in using the LLL states instead of the Wannier-Qi states of the twisted bilayer is that we ignore any Berry-curvature inhomogeneity.

IV Exciton vortex lattice in the lowest Landau level

The perpendicular magnetic field seen by the exciton order parameter $\Delta({\mathbf{r}})=\langle c^{\dagger}_{+,{\mathbf{r}}}c_{-,{\mathbf{r}}}\rangle$ induces a vortex lattice in the order parameter. Since $\Delta({\mathbf{r}})$ essentially behaves like a charge $q=2e$ object in a magnetic field, the solution to this vortex lattice may be obtained by solving the Ginzburg Landau (GL) equation for $\Delta({\mathbf{r}})$. For analytical tractability, we focus on vortex lattice solution of the linearized GL equation. Our solution is exact only at the upper critical field $H_{c2}$ of the corresponding superconductor, but we expect our results to be valid more generally. In this limit, the problem reduces to the solving the Schrodinger equation for a single particle of charge $2e$. In the Landau gauge ${\mathbf{A}}=Bx\hat{y}$, this solution is given by

[TABLE]

The exciton vortex lattice we consider has $2\pi$ flux through each plaquette of the square lattice of side $a$, i.e, $a^{2}=2\pi l_{B}^{2}$. Since each elementary vortex carries a flux of $\pi$, we therefore expect two elementary vortices within a plaquette. Inspired by the computation of similar vortex patterns for the superconducting order parameter in Ref. Mishmash et al. (2018), we choose the vortex lattice wavefunction to be symmetric under magnetic translations ${\mathcal{T}}_{1}={\mathcal{T}}(a\hat{y})$ and ${\mathcal{T}}_{2}={\mathcal{T}}\left(\frac{a}{2}(\hat{x}+\hat{y})\right)$. Note that the magnetic translation operators for a particle of charge $q$ in a magnetic field $\mathbf{B}$ satisfy the following algebra:

[TABLE]

Since $\Delta({\mathbf{r}})$ is a charge $q=2e$ order parameter, we have $q\mathbf{B}\cdot({\mathbf{R}}_{1}\times{\mathbf{R}}_{2})/\hbar=(2e)B(a^{2}/2)/\hbar=2\pi$ implying that ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$ commute with each other. Being magnetic translation operators, they commute with the free Hamiltonian of a particle of charge $q$. Since our goal is to express $\Delta({\mathbf{r}})$ in the Bloch basis, where the eigenstates of particles with charge $q=e$ are invariant under the square lattice translations ${\mathcal{T}}_{1}={\mathcal{T}}(a\hat{y})$ and ${\mathcal{T}}_{3}={\mathcal{T}}(a\hat{x})$, we also choose the phases of ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$ such that ${\mathcal{T}}_{3}={\mathcal{T}}_{2}^{2}{\mathcal{T}}_{1}^{-1}$ is identically satisfied. Consistent with these conditions, we find that

[TABLE]

Now we impose the magnetic translation symmetry requirements on $\Delta({\mathbf{r}})$. For ${\mathcal{T}}_{1}$, we have

[TABLE]

where we have defined $Q=2\pi/a$. Therefore, we can write

[TABLE]

For ${\mathcal{T}}_{2}$, we have, using $2\pi l_{B}^{2}=a^{2}$ or equivalently $a=Ql_{B}^{2}$,

[TABLE]

Therefore, we have the following form of $\Delta({\mathbf{r}})$:

[TABLE]

We can now find the projection of $\Delta({\mathbf{r}})$ on the single particle Bloch wave-functions. We focus on the lowest Landau level since we are only interested in $C=\pm 1$ bands. We define

[TABLE]

where $\phi_{\pm,{\mathbf{k}}}({\mathbf{r}})$ are the Bloch wave-functions defined in Eq. (1) in the main text. Given the symmetry of $\Delta({\mathbf{r}})$ under magnetic translations ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{3}$, we expect it to be diagonal in Bloch space. Indeed, we find that $\Delta_{00}({\mathbf{k}},{\mathbf{k}}^{\prime})=\Delta_{{\mathbf{k}}}\delta_{{\mathbf{k}},{\mathbf{k}}^{\prime}}$, where

[TABLE]

where $\Delta_{0}=\frac{C_{0}}{\sqrt{2}}$ is a measure of the overall strength of the exciton vortex lattice order parameter.

We can recast $\Delta_{\mathbf{k}}$ in terms of the Jacobi theta function as follows

[TABLE]

The Jacobi theta function has zeros at $z=m+n\tau+1/2+\tau/2$, where $m,n\in\mathbb{Z}$. Restricting to the first BZ, we find that $\Delta_{{\mathbf{k}}}=0$ at ${\mathbf{k}}=\pm{\mathbf{k}}_{0}$, with ${\mathbf{k}}_{0}=(\pi/2,-\pi/2)$. Further, a power series expansion about the zeros shows that $\Delta_{\pm{\mathbf{k}}_{0}+{\mathbf{q}}}=\pm A(q_{x}-iq_{y})+O(q^{2})$ (for some $A\in\mathbb{C}$), indicating that both nodes have the same chirality. The presence of these two nodes in the BZ, which is a topological requirement of $\Delta_{\mathbf{k}}$ arising from hybridization of bands with $C=\pm 1$, is confirmed by plotting the absolute value of $\Delta_{\mathbf{k}}$ in Fig. 2 in the main text.

The nodes in the exciton order parameter are intimately related to the Dirac cones of the $C_{2z}T$-symmetric single-valley Moiré Hamiltonian, and the associated Wannier obstruction Po et al. (2018); Zou et al. (2018); Po et al. (2018); Song et al. (2018); Yuan and Fu (2018); Kang and Vafek (2018); Koshino et al. (2018). To see this, consider a free fermion Hamiltonian with two bands that are isolated from the other bands, such that the momentum space Hamiltonian projected onto the two isolated bands takes the form

[TABLE]

where $|u_{\pm,\textbf{k}}\rangle$ are the periodic parts of the Bloch states in the two-band subspace. The dispersion of the two bands is given by $\pm\sqrt{(\epsilon_{\textbf{k}}-\mu)^{2}+|\Delta_{\textbf{k}}|^{2}}$ (for simplicity, but without loss of generality, we use a particle-hole symmetric projected Hamiltonian). We consider the situation where $|u_{\pm,\textbf{k}}\rangle$ has Chern number $\pm 1$, i.e. $\frac{1}{2\pi}\int_{\textbf{k}}\nabla\times\textbf{A}_{\pm}=\pm 1$, where $\textbf{A}_{\pm}=-i\langle u_{\pm,\textbf{k}}|\nabla|u_{\pm,\textbf{k}}\rangle$. Using the same reasoning as for the exciton order parameter, and the fact that the Hamiltonian has to be periodic over the Brillouin zone, we conclude that $\Delta_{\textbf{k}}$ has two zeros in the Brillouin zone around which its phase winds by $2\pi$. Now imagine tuning $\mu$ from minus infinity to plus infinity. In this process the Chern number of the lowest energy band changes from $+1$ to $-1$, which is only possible if the energy gap between the two bands closes for intermediate values of $\mu$. From the band dispersion, we see that this gap closing will occur precisely at the momenta where the zeros of $\Delta_{\textbf{k}}$ are located. At these points, the nodes in $\Delta_{\textbf{k}}$ give rise to Dirac cones with the same chirality.

V Analytical Hartree-Fock energetics for the lowest Landau level model

In this section, we compute the Hartree-Fock (HF) energy in a variational Slater determinant state for the Hamiltonian $H=H^{i}+H^{p}$, which we recall for completeness.

[TABLE]

In Eq. (38), $u_{0}({\mathbf{q}})$ is the symmetric part of the interaction, while $u_{1}({\mathbf{q}})$ represents the anisotropy. We evaluate $\langle H\rangle$ for a variational Slater determinant state, which can capture the fully valley polarizd Chern insulator, the unpolarized and partially polarized metallic states, and the exciton condensate (both uniform and a vortex lattice) in different limits. Our variational state $\ket{\psi_{MF}}$ may be taken to be the Slater determinant ground state of a mean field Hamiltonian of the form $H_{MF}=\sum_{{\mathbf{k}},\tau,\tau^{\prime}}c^{\dagger}_{{\mathbf{k}},\tau}h_{\tau,\tau^{\prime}}({\mathbf{k}})c_{{\mathbf{k}},\tau^{\prime}}$, where

[TABLE]

Such a state $\ket{\psi_{MF}}$ is characterized by two variational parameters, the polarization $P_{v}$ (determined by $h$) and the strength of the excitonic order parameter $\Delta_{0}$. In the $|h|\gg\text{max}\{|\varepsilon_{{\mathbf{k}}}|\}$ and $\Delta_{{\mathbf{k}}}=0$ limit, our Slater determinant state is thus fully valley polarized with $|P_{v}|=1$, while for $h=\Delta_{{\mathbf{k}}}=0$ we have an unpolarized metal. For $h\neq 0$ and $\Delta_{{\mathbf{k}}}=0$, the state is a partially polarized metal with different chemical potential for the $\tau=\pm$ valleys. For $h=0$ and $\Delta_{\mathbf{k}}\neq 0$, the state is an excitonic condensate or vortex lattice (depending on the precise structure of $\Delta_{\mathbf{k}}$), with $P_{v}=0$ fixed by the discrete $z\rightarrow-z$ symmetry of $H_{MF}$. The most general state will have both $h$ and $\Delta_{\mathbf{k}}$ non-zero.

Evaluating the covariance matrix for $\ket{\psi_{MF}}$ with a chemical potential $\mu$ which fixes the filling of the mean-field bands gives:

[TABLE]

One can indeed check that $P_{\tau,\tau^{\prime}}({\mathbf{k}})$ is a projector matrix, i.e, $P^{2}=P$, as expected for a Slater determinant state. We use these correlators to evaluate $\langle H\rangle$ via Wick’s theorem.

[TABLE]

where we have taken the thermodynamic limit, set the lattice spacing $a=1$ and used the notation $\int_{{\mathbf{k}}}=\int\frac{d^{2}k}{(2\pi)^{2}}$ to denote integration over the first BZ. We now focus on different limits where we can analytically compute the regularized energy density $e^{HF}(P_{v},\Delta_{0})\equiv\lim_{N_{\phi}\rightarrow\infty}\left(\frac{E^{HF}}{N_{\phi}}-\frac{u_{0}(\textbf{0})}{2}\right)$ (where we have subtracted the formally infinite self-energy contribution that is canceled by the positive background) and get physical intuition about the phase diagram and stability of the different phases.

V.1 Competition between metal and valley polarized states

First, we focus on the competition between the metallic state and valley polarized state (setting $\Delta_{0}=0$). The anisotropic part of the interaction $u_{1}({\mathbf{q}})$, while crucial for the excitonic order parameter, does not play a prominent role here other than altering phase boundaries slightly, so we set it to zero for simplicity. In this case, the covariance matrix takes the form

[TABLE]

with a separate Fermi energy $\varepsilon^{\tau}_{F}=\mu+\tau h$ for the two bands ($\tau=\pm 1$). Therefore, the regularized HF energy density is given by

[TABLE]

To intuitively understand the physics, let us consider two extreme limits. For the fully valley polarized state, one of bands is completely full while the other is completely empty. Hence, Eq. (43) evaluates to

[TABLE]

For the unpolarized metal, $f^{+}_{\mathbf{k}}=f^{-}_{\mathbf{k}}\equiv f_{\mathbf{k}}=\Theta(-\varepsilon_{\mathbf{k}})$. Defining $g({\mathbf{q}})=\int_{{\mathbf{k}}}f_{{\mathbf{k}}+{\mathbf{q}}/2}f_{{\mathbf{k}}-{\mathbf{q}}/2}$, Eq. (43) evaluates to

[TABLE]

The function $g({\mathbf{q}})$ is proportional to the overlap of the Fermi surface with itself when shifted by ${\mathbf{q}}$. Hence, $g({\mathbf{q}})$ has a maximum value of 1 at ${\mathbf{q}}=0$ and decreases with ${\mathbf{q}}$ till ${\mathbf{q}}$ is half a reciprocal lattice vector. Since $u_{0}({\mathbf{q}})$ contains the Landau level projection factor $F^{2}({\mathbf{q}})=e^{-q^{2}l_{B}^{2}/2}$, the main contribution to the interaction term comes from $g({\mathbf{q}})$ close to zero, which implies that the unpolarized metal has higher energy than the valley polarized state. In other words, interaction favors valley polarization. On the other hand, the kinetic term from the periodic potential favors the metal, as a full dispersing band costs more energy than two half-filled bands.

V.2 Stability of valley polarized insulator and unpolarized metal to exciton vortex lattice

Following our previous discussion about the metallic phase and the valley polarized insulator, we need to establish that both these phases are stable to an excitonic phase with non-zero $\Delta({\mathbf{r}})$ in presence of anisotropy in the interactions ($u_{1}({\mathbf{q}})\neq 0$). For two $C=+1$ bands, it is well-known that an infinitesimal anisotropy will drive exciton condensation with uniform magnitude. As argued in the main text, our excitonic order parameter $\Delta({\mathbf{r}})$ formed which has electrons from the $C=+1$ band and holes from the $C=-1$ band, will behave like a superconducting order parameter in presence of a uniform magnetic field. Therefore, we can rule out a uniform exciton condensate, but an exciton vortex lattice indeed remains a distinct possibility. Below, we argue that such a phase is also energetically more expensive as long as the anisotropy is small enough.

We start off with the fully valley polarized state, corresponding to a large $h$. We now add a small $\Delta_{\mathbf{k}}$ to see if we gain energy in presence of an arbitrarily weak anisotropy $u_{1}$, while keeping the filling fixed to half. In this limit, the lower ($\beta$) band is still full while the upper ($\alpha$) band is empty, so we can write $\Theta_{{\mathbf{k}},\alpha}=0$ and $\Theta_{{\mathbf{k}},\beta}=1$. We can write the covariance matrix from Eq. (V) as follows:

[TABLE]

Using the form of $P$ from Eq. (46) and writing out the terms in Eq. (LABEL:eq:EHF) in terms of $u_{0}$ and $u_{1}$, we find that (using $\operatorname{Tr}(P({\mathbf{k}}))=1$):

[TABLE]

To check the stability perturbatively, we expand in powers of $|\Delta_{0}|/h$ and consider the difference of energy density $e^{HF}(1,\Delta_{0})$ and $e^{HF}(1,0)$ for the fully polarized state upto quadratic order.

[TABLE]

The first two terms raise the energy, while the last term lowers the energy of our variational state with respect to the fully valley polarized state. If both bands had the same Chern number, a spatially uniform $\Delta$ is allowed so that $\Delta_{\mathbf{k}}=\Delta_{0}~{}\forall~{}{\mathbf{k}}$. In this case, the first term in Eq. (V.2) does not contribute, and we have

[TABLE]

For a uniform exciton condensate with a spatially uniform $\Delta$, we have $\Delta_{\mathbf{k}}=\Delta_{0}~{}\forall~{}{\mathbf{k}}$. Therefore,

[TABLE]

Since the Landau level projection adds a factor of $e^{-{\mathbf{q}}^{2}l_{B}^{2}/2}$ to the bare anisotropy, we have $u_{1}({\mathbf{q}})<u_{1}(\mathbf{0})~{}\forall~{}{\mathbf{q}}\neq\mathbf{0}$. This negative difference in energy density precisely corresponds to the instability of the fully valley polarized phase of the conventional QHFM to uniform intervalley coherence when $u_{1}({\mathbf{q}})>0$. However, for any vortex lattice structure, necessitated by topological constraints of hybridizing opposite Chern bands, we have $\Delta_{\mathbf{k}}$ which is a function of ${\mathbf{k}}$, Therefore, when $u_{1}$ is sufficiently small compared to $u_{0}$ the vortex lattice state has a higher energy than the parent insulator, regardless of the exact nature of the microscopic interactions (as long as both are repulsive). This implies that the fully valley polarized state is robust to the vortex lattice phase.

We now carry out the previous analysis for the unpolarized metal, the second state of our interest. In this case, $P_{v}=0$ (obtained by setting $h=0$), so the covariance matrix is given by

[TABLE]

In this case, the HF energy density evaluates to

[TABLE]

where the labels $\pm$ are shorthand for momenta ${\mathbf{k}}\pm{\mathbf{q}}/2$. We observe that at a total filling of $\nu=1$, the Fermi surface of the $\beta$ band is identical to the Fermi surface of the $\alpha$ band shifted by ${\mathbf{Q}}=(\pi,\pi)$. To prove this, we use $\varepsilon_{{\mathbf{k}}+{\mathbf{Q}}}=-\varepsilon_{\mathbf{k}}$ and $\Delta_{{\mathbf{k}}+{\mathbf{Q}}}=-\Delta_{{\mathbf{k}}}$ (as can be seen from the analytical form of $\Delta_{\mathbf{k}}$ in Eq. (35)).

[TABLE]

This has the very important consequence that the chemical potential $\mu$ is fixed to zero at half filling, and the system behaves like a compensated semi-metal with equally sized electron and hole Fermi surfaces. To analyze the energetics, it is convenient to define two subsets of the BZ. Let $FS_{0}$ be the original diamond shaped Fermi surface of the bands at $\Delta_{\mathbf{k}}=0$, defined by the contours $|k_{x}\pm k_{y}|=\pi$ in the first BZ. Then we define (see Fig. 6)

[TABLE]

Note that ${\mathbf{k}}\in S_{1}$ implies ${\mathbf{k}}+{\mathbf{Q}}\in S_{2}$, so the area of $S_{i}$ ($i=1,2$) in the BZ, which we denote by $\mathcal{A}_{S_{i}}$ are equal. Using these, we compute the kinectic energy term of the vortex lattice phase.

[TABLE]

The second term is positive, and denotes the increase of kinetic energy of our variational state by virtue of distorting the bands. The interaction term can also be analogously split into contributions coming from the original Fermi surface, and those coming from Fermi surface distortions. Further, we need to consider the overlap of Fermi surfaces shifted by a momenta of ${\mathbf{q}}$, but the Landau level projection factors imply that only the overlap at $q\approx 0$ is important when $a\ll l_{B}$. While our lattice has $a/l_{B}=\frac{1}{\sqrt{2\pi}}$ so we are not strictly in this limit, it is nevertheless instructive to look at, as the Fermi surface overlaps can be succinctly expressed in terms of $S_{1}$ and $S_{2}$. Adding all contributions, we finally find that the energy density difference is given by:

[TABLE]

where we have approximated $u_{1}({\mathbf{q}})$ by $u_{1}(l_{B}^{-1})$ to avoid potential singularities at $q=0$. For weakly anisotropic repulsive interactions, $u_{1}>0$ so both terms raise the energy of the vortex lattice variational state with respect to the unpolarized metal. For $u_{1}<0$, the first term raises the energy while the second one lowers it. However, in both cases for a small enough $u_{1}$, the vortex lattice state still has a higher energy and will not be favorable.

VI Numerical Hartree-Fock analysis of magic angle graphene

To confirm that the physical picture discussed in the main text in terms of a Lowest Landau model indeed applies to magic angle graphene, we have numerically solved the Hartree-Fock self-consistency equations in the case where there is a sublattice splitting of $15$ meV on the top graphene layer. We used a dual gate-screened Coulomb potential with a gate distance of $20$ nm, and a dielectric constant of $\epsilon=9.5$. The twist angle was $\theta=1.05^{\circ}$. The simulations were done on a $24\times 24$ momentum grid, keeping six BM bands per spin and valley. The setup of the simulation is exactly the same as that in Ref. Bultinck et al. (2019), and we refer to that paper for more details. We would like to point out that the only assumption that went into the numerics is that the ground state does not break translation symmetry at the Moire scale – every other type of symmetry breaking is allowed to occur and no bias is introduced.

In Fig. 7, we show the self-consistent Hartree-Fock insulating band spectrum at charge neutrality. It has a large bandgap of $\sim 30$ meV. The four lowest conduction bands, which are the bands we focus on in the main text, also remain very flat and do not touch the remote bands. The bandwidth is of the same order as that of the bare BM bands. We also find that at charge neutrality, no symmetries are broken spontaneously. Also, if we initialize the algorithm with the state corresponding to the insulating ground state of the BM model with sublattice splitting, we find that convergence is achieved after only 5 to 10 iterations, which means that this state is already very close to the self-consistent solution. Therefore, the bands do not change significantly in going from the BM ground state to the self-consistent solution.

In Fig. 8, the self-consistent band Hartree-Fock band spectrum is shown at filling $\nu=3$. In 8(a), the total band spectrum is shown, while in 8(b) only the spin up bands are shown, and in 8(c) only the spin down bands are shown. Importantly, we find that system has a clear bandgap, and unbroken valley U$(1)$ symmetry. From Fig. 8 we also clearly see that one spin and valley polarized band of the eight bands around charge neutrality remains unoccupied. This band spectrum therefore agrees with the general physical scenario put forward in the main text.

VII Orbital moment and valley-Zeeman effect

The momentum-dependent orbital moment of electrons in a band labeled by $\alpha$ is given by Xiao et al. (2010); Chang and Niu (1996)

[TABLE]

where $H({\mathbf{k}})$ is the corresponding Bloch Hamiltonian with eigenstates $|u_{\alpha,{\mathbf{k}}}\rangle$ and eigenvalues $\varepsilon_{\alpha,{\mathbf{k}}}$. This orbital moment couples linearly to the out-of-plane component $B^{z}$ of the magnetic field via the orbital Zeeman term

[TABLE]

The average orbital $g$-factor reported in the main text is given by

[TABLE]

where $\mu_{B}$ is the Bohr magneton and $A_{mBZ}$ is the area of the mBZ.