On the nature of the correlated insulator states in twisted bilayer graphene

TL;DR

This study uses self-consistent Hartree-Fock calculations to analyze the correlated insulator states in twisted bilayer graphene, revealing that gaps can open without breaking certain symmetries and that insulating states at specific fillings do not always show a quantized Hall effect.

Contribution

It provides a detailed theoretical analysis of the symmetry properties and conditions for gap opening in twisted bilayer graphene's correlated insulator states.

Findings

Gaps open at neutrality over a wide twist angle range.

Insulating states at certain fillings may not exhibit a quantized Hall effect.

Gapped states can form without breaking the system's  symmetry.

Abstract

We use self-consistent Hartree-Fock calculations performed in the full $\pi$-band Hilbert space to assess the nature of the recently discovered correlated insulator states in magic-angle twisted bilayer graphene (TBG). We find that gaps between the flat conduction and valence bands open at neutrality over a wide range of twist angles, sometimes without breaking the system's valley projected ${\cal C}_{2}{\cal T}$ symmetry. Broken spin/valley flavor symmetries then enable gapped states to form not only at neutrality, but also at total moir\'e band filling $n = \pm p/4$ with integer $p = 1, 2, 3$, when the twist angle is close to the magic value at which the flat bands are most narrow. Because the magic-angle flat band quasiparticles are isolated from remote band quasiparticles only for effective dielectric constants larger than $ \sim 20$, the gapped states do not necessarily break \CT symmetry and as a consequence the insulating states at $n = \pm 1/4$ and $n = \pm 3/4$ need not exhibit a quantized anomalous Hall effect.

Tables

Table 1. Table 1: Total energy per moiré unit cell at neutrality in units of meV calculated in the non-interacting (NI) and self-consistent Hartree-Fock (SCHF) ground states for θ = 1.1 ∘ 𝜃 superscript 1.1 \theta=1.1^{\circ} and ϵ = 5 italic-ϵ 5 \epsilon=5 . Gapped insulating states with and without broken 𝒞 2 ​ 𝒯 subscript 𝒞 2 𝒯 \mathcal{C}_{2}\mathcal{T} compete closely, differing in energy by less than 1 meV per moiré cell. Note that Hartree energies do not play an important role in selecting the broken symmetry state. Both intralayer and interlayer energies grow with the momentum space cut-off, but the difference in energy between non-interacting and interacting states converges.

E(meV)	Intra-layer			Inter-layer			Total
	Hopping	Hartree	Fock	Hopping	Hartree	Fock
NI	1417	0	-44	-2636	0	-278	-1541
SCHF	2372	0	-109	-3452	0	-469	-1658

Equations

$$\frac{A_m}{A}\sum _{\vec{k}}^{mBZ}⟨\mathbf{k}+\mathbf{G},l^{\prime },s^{\prime }|{\Sigma }^F|\mathbf{k}+\mathbf{G},l,s⟩=\sum _{ij}{\Delta }_{\mathbf{G},ij}{\sigma }_{s^{\prime }s}^i{\tau }_{l^{\prime }l}^j.$$

$${\hat{\mathcal{H}}}_0^{\mathrm{K}}=\left(\begin{array}{cc}{h}_{\theta /2}(\mathbf{k})& h_T(\mathbf{r})\\ h_T^{†}(\mathbf{r})& h_{-\theta /2}({\mathbf{k}}^{\prime })\end{array}\right)$$

$$h_{\theta }(\mathbf{k})=-\hslash v_D|\bar{\mathbf{k}}|\left(\begin{array}{cc}0& e^{i({\theta }_{\bar{\mathbf{k}}}-\theta )}\\ e^{-i({\theta }_{\bar{\mathbf{k}}}-\theta )}& 0\end{array}\right),$$

$$h_T(\mathbf{r})=\sum _{j=0}^3{T}_j{e}^{-i{\mathbf{q}}_j\cdot \mathbf{r}}$$

$$T_j={\omega }_0{\sigma }_0+{\omega }_1\cos (j\varphi ){\sigma }_x+{\omega }_1\sin (j\varphi ){\sigma }_y$$

$${\Sigma }_{\alpha ,\mathbf{G};\beta ,{\mathbf{G}}^{\prime }}^H(\mathbf{k})=\frac{1}{A}\sum _{{\alpha }^{\prime }}{V}_{{\alpha }^{\prime }\alpha }({\mathbf{G}}^{\prime }-\mathbf{G})\delta {\rho }_{{\alpha }^{\prime }{\alpha }^{\prime }}(\mathbf{G}-{\mathbf{G}}^{\prime }){\delta }_{\alpha \beta }$$

$${\Sigma }_{\alpha ,\mathbf{G};\beta ,{\mathbf{G}}^{\prime }}^F(\mathbf{k})=-\frac{1}{A}\sum _{{\mathbf{G}}^{\prime \prime },{\mathbf{k}}^{\prime }}$$

$$\times \delta {\rho }_{\alpha ,\mathbf{G}+{\mathbf{G}}^{\prime \prime };\beta ,{\mathbf{G}}^{\prime }+{\mathbf{G}}^{\prime \prime }}({\mathbf{k}}^{\prime })$$

$${\rho }_{\alpha ,\mathbf{G};\beta ,{\mathbf{G}}^{\prime }}(\mathbf{k})=\sum _n{z}_{\beta ,{\mathbf{G}}^{\prime },\mathbf{k}}^{n\ast }{z}_{\alpha ,\mathbf{G},\mathbf{k}}^n$$

$$|{\psi }_{n,\mathbf{k}}⟩=\sum _{\alpha ,\mathbf{G}}{z}_{\alpha ,\mathbf{G},\mathbf{k}}^n|{\psi }_{\alpha ,\mathbf{G}+\mathbf{k}}⟩,$$

$$\mathcal{H}={\mathcal{H}}_0^{\mathrm{K}}+{\Sigma }^{\mathrm{H}}+{\Sigma }^{\mathrm{F}}.$$

$$\frac{A_m}{A}\sum _{\vec{k}}^{BZ}⟨\mathbf{k}+{\mathbf{G}}^{\prime },l{,}^{\prime }{s}^{\prime }|{\Sigma }^F|\mathbf{k}+\mathbf{G},l,s⟩=\sum _{ij}{\Delta }_{{\mathbf{G}}^{\prime }\mathbf{G},ij}{\sigma }_{s^{\prime }s}^i{\tau }_{l^{\prime }l}^j.$$

$${\Delta }_{{\mathbf{G}}^{\prime }\mathbf{G},ij}={\Delta }_{{\mathbf{G}}^{\prime }\mathbf{G},ij}^R+i{\Delta }_{{\mathbf{G}}^{\prime }\mathbf{G},ij}^I,$$

$$U_X=-\frac{A_m}{A}\sum _{\mathbf{k}}⟨{\psi }_{c,\mathbf{k}}^{\text{NSC}}|{\Sigma }^F(\delta {\rho }_c^{\text{NSC}})|{\psi }_{c,\mathbf{k}}^{\text{NSC}}⟩,$$

$$H_0(\mathbf{m})=-\frac{1}{2}\hslash v_D{k}_{\theta }{\sigma }_y{\tau }_z+{\omega }_0{\sigma }_0{\tau }_x+{\omega }_1{\sigma }_x{\tau }_x.$$

$$E_{\mathbf{m}}=\pm {\omega }_1\pm \sqrt{{\hslash }^2{v}_D^2{k}_{\theta }^2/4+{\omega }_0^2}.$$

$$H_0(\mathbf{\gamma })=\left(\begin{array}{cccccc}{h}_{-\theta /2}({\Gamma }_0^b)& T_0& 0& 0& 0& T_1\\ T_0& h_{\theta /2}({\Gamma }_0^t)& T_{-1}& 0& 0& 0\\ 0& T_{-1}& h_{-\theta /2}({\Gamma }_1^b)& T_1& 0& 0\\ 0& 0& T_1& h_{\theta /2}({\Gamma }_2^t)& T_0& 0\\ 0& 0& 0& T_0& h_{-\theta /2}({\Gamma }_2^b)& T_{-1}\\ T_1& 0& 0& 0& T_{-1}& h_{-\theta /2}({\Gamma }_1^t)\end{array}\right),$$

$$H_0(\mathbf{\gamma })=H_0^{+}(\mathbf{\gamma })+H_0^{-}(\mathbf{\gamma })+H_0^{+-}(\mathbf{\gamma })+H_0^{-+}(\mathbf{\gamma })$$

$$H_0^{\pm }(\mathbf{\gamma })=\hslash v_D|{\mathbf{k}}_{\mathbf{\gamma }}|\cdot \left(\begin{array}{cccccc}\pm 1+2b& 0& 0& 0& 0& 0\\ 0& \pm 1-2b& 0& 0& 0& 0\\ 0& 0& \pm 1+\sqrt{3a^2+b^2}& 0& 0& 0\\ 0& 0& 0& \pm 1+\sqrt{3a^2+b^2}& 0& 0\\ 0& 0& 0& 0& \pm 1-\sqrt{3a^2+b^2}& 0\\ 0& 0& 0& 0& 0& \pm 1-\sqrt{3a^2+b^2}\end{array}\right).$$

$$H_0^{+-}(\mathbf{\gamma })=\frac{i}{\sqrt{3}}\hslash v_D|{\mathbf{k}}_{\mathbf{\gamma }}|a\cdot \left(\begin{array}{cccccc}-2& 0& 0& 0& 0& 0\\ 0& 2& 0& 0& 0& 0\\ 0& 0& -1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& -1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right).$$

Full text

Introduction.— A small relative twist between adjacent graphene layers produces a triangular lattice moiré pattern with a spatial periodicity that is inversely related to twist angle. It was noticed MorellFlatBand ; BMModel some years ago that at a series of magic twist angles $\theta$, the moiré pattern yields very flat low-energy bands that promise strong electronic correlations. This promise has now been realized thanks to recent experimental studies Tutuc2017 ; CaoInsulator ; CaoSuper ; YoungDean ; Efetov ; Gordon ; Andrei ; Pasupathy ; Yazdani of bilayers with accurately controlled twist angles that exhibit interaction-induced insulating ground states at moiré band BMModel filling factors $n=\pm p/4$, where $p=-3,\ldots,3$ is the total charge per moiré unit cell. The insulating states are flanked by superconducting domes CaoSuper ; Efetov . This exciting discovery has inspired a flurry of theoretical work FuModel ; VishwanathMarch ; Scalettar ; Phillips ; SpinLiquid ; WangModel ; MottAF ; Yang ; PALee ; Vafek ; Rademaker ; KoshinoFu ; Bascones ; Ochi ; FuJune ; AFMonHoneyComb ; XuBalents ; FuSuper ; Juricic ; Das ; Super1 ; Super2 ; Super3 ; Super4 directed toward achieving a more complete understanding of the insulating states and their superconducting satellites. Previous work on the insulating states has been based mainly on an indirect approach that starts by identifying effective lattice models for the flat moiré bands, and then combines these with generalized Hubbard models to address interaction phenomena. In this Letter we explore a different approach.

At small twist angles the electronic structure of twisted bilayer graphene can be accurately described using a continuum model BMModel ; SantosEarly in which single-particle electronic states with a four-level spin/valley internal flavor degree-of-freedom are approximated by four-component envelope function spinors that specify $\pi$-orbital amplitudes on the bilayer’s four sublattices. The simplest version of the continuum model BMModel adds a spatially periodic interlayer hopping term to isolated layer $\pi$-orbital Dirac models. This moiré band Hamiltonian is spin-independent, and its projections onto graphene’s two valleys are related by time-reversal symmetry. Up to an overall energy scale, its spectrum depends on a single twist-angle dependent parameter $\alpha=w/\hbar vk_{\theta}$ where $w\approx 110$ meV is an inter-layer tunneling amplitude, $v\approx 10^{6}$ m/s is the Dirac velocity, $k_{\theta}=2K\sin(\theta/2)$ is the momentum separation between the Brillouin-zone (BZ) corners in different layers, and $K$ is the single-layer BZ corner momentum magnitude.

The valley-projected moiré flat bands in TBG occur not singly, but in valence/conduction pairs connected by two symmetry-protected linear Dirac band crossings. Importantly the two Dirac points of weakly coupled bilayers carry the same chirality. This property implies VishwanathWannierObstructions ; VishwanathMarch ; Bernevig that the moiré flat bands can be described only by tight binding models with at the very least four orbitals per spin/valley flavor per moiré unit cell. Recent work VishwanathAugust suggests that faithful descriptions of interaction physics using generalized Hubbard models may require the inclusion of at least eight bands per flavor, limiting the motivation for, approximate lattice models.

In this Letter we report on a study of the correlation induced insulator states in magic angle TBG that starts directly from the moiré band continuum model and accounts for the long-range of the Coulomb interaction between electrons. Our principle results are summarized in Fig. 1. In this figure the largest values of interaction strength parameter $\epsilon^{-1}$ correspond to screening by a surrounding hexagonal boron nitride dielectric only. In practice interactions are sample dependent and always weakened by nearby gates. We find that gapped states occur at neutrality when the effective fine structure constant $\alpha^{*}=e^{2}/\epsilon\hbar v_{D}^{*}$ exceeds $\sim 0.4$ for twist angle $\theta$ above the magic value, and almost always for twists that are smaller. Here $v_{D}^{*}$ is the reduced Dirac velocity of the twisted bilayer which vanishes as the magic angle is approached BMModel . We attribute the smaller value of the critical fine structure constant in twisted bilayer graphene than in the corresponding single-layer graphene ($\alpha^{*}\sim 1$ Min ) calculations to the non-uniform spatial distribution of flat band orbitals, which enhances interaction effects. As in single layer grapheneAllanMonolayer , these gapped states break $\mathcal{C}_{2}\mathcal{T}$ symmetry and have non-zero Berry curvatures, whereas the the gapless states preserve $\mathcal{C}_{2}\mathcal{T}$ symmetry. For $n=\pm p/4\neq 0$, gapped states are enabled by broken spin/valley flavor symmetries and occur over a much narrower range of twist angles. When the interaction strength is sufficiently strong, gapped states can be opened without breaking the $\mathcal{C}_{2}\mathcal{T}$ symmetry that protects band crossings when the Hamiltonian is projected onto the strongly correlated flat bands. This property is significant because it has implications for the occurrence of quantized anomalous Hall effects at band filings $n=\pm 1/4$ and $n=\pm 3/4$.

Mean field theory.—Our theoretical approach is guided by the experimental CaoInsulator discovery of insulating states in magic angle TBG that are naturally explained by broken symmetries that lift the four-fold spin/valley degeneracy of the band Hamiltonian, and do not require broken translational symmetry. In most cases the ground states of insulators can be described using Hartree-Fock mean-field theory. The difficulty in the TBG case compared to the familiar case of atomic-scale insulators, is that the Hamiltonian does not contain strong attractive potential terms within each unit cell that select particular high-weight atomic or ionic configurations. To understand the nature of the insulating states, we must perform unbiased Hartree-Fock calculations in the full footnote $\pi$-orbital Hilbert space, placing no restrictions on the flavor or position dependence of the model’s four-component envelope function spinors.

A typical self-consistent Hartree-Fock calculation result, for the point $\epsilon=5$ and $\theta=1.1^{\circ}$ marked by a red circle in Fig. 1(a), is summarized in Fig. 2, which illustrates quasiparticle dispersion, topology, and $\mathcal{C}_{2}\mathcal{T}$ breaking order parameters, and in Table 1 in which the ground state energies of the non-interacting and interacting cases are partitioned into intralayer and interlayer tunneling and interaction contributions. The technical details of these calculations are described in the Supplemental Material. All energies are expressed relative to the energy of the non-interacting state in the absence of inter-layer tunneling, and a neutralizing background charge density is assumed. As shown in Fig. 2, we find separate self-consistent gapped solutions with and without $\mathcal{C}_{2}\mathcal{T}$ symmetry breaking. The $\mathcal{C}_{2}\mathcal{T}$ symmetry broken solution features moiré bands with well-defined non-zero Berry curvatures and sublattice polarizations $\rho_{B}-\rho_{A}$, and is lowest in energy in most regions of the phase diagram as shown in Fig. 1(a).

In Table 1 we note that the ratio of the cost in intra-layer tunneling energy, to the energy gain from inter-layer tunneling in the non-interacting ground state is $1:2$, the ratio that is obtained when inter-layer tunneling is treated as a weak perturbation. This observation is consistent with the property BMModel that the first magic angle in twisted bilayer graphene is accurately predicted by perturbation theory. Secondly we observe that at neutrality both non-interacting and interacting ground states have almost uniform charge density, even though the flat band wavefunctions are spatially peaked near AA positions in the moiré pattern. The absence of a Hartree energy at neutrality is related to the moiré band Hamiltonian’s approximate particle-hole symmetry, and is quite distinct from what would be obtained if the Hilbert space were truncated to include only the lowest energy flat bands. Finally we note that the condensation energy of the gapped state, which we define as the difference between ground state energy and the expectation value of the Hamiltonian in the non-interacting ground state, is $\sim 117$ meV per moiré period and originates mainly from enhanced interlayer exchange energies. At this value of $\epsilon^{-1}$ and $\theta$, total energy minimization including interactions adjusts the ground state so as to enhance interlayer tunneling and interlayer exchange energies at a cost in the intra-layer hoping energy, and a substantial part of the ground state rearrangement occurs in higher energy (remote) valence bands. A description of interaction physics in terms of the single-particle flat bands alone is sufficient only for $\epsilon^{-1}\lesssim 0.04$. Even in this case, however, it is necessary to include the Hartree GuineaInteractions and exchange self-energies from the frozen negative energy sea which lower the energies of states near the moiré BZ (mBZ) $\kappa$ and $\kappa^{\prime}$ points (see Fig. 2) relative to those near $\gamma$ and therefore contributes to quasiparticle band dispersion.

Band topology, insulating states, and broken $\mathcal{C}_{2}\mathcal{T}$ symmetry.—The moiré band model BMModel captures the microscopic tight-binding model’s $\mathcal{D}_{6}$, time reversal $\mathcal{T}$, and $U(1)$ valley symmetries. When the Hamiltonian is projected to a single valley, both $\mathcal{T}$ and $\mathcal{C}_{2}=(\mathcal{C}_{6})^{3}$a are lost because they map states between valleys. We are left only with the combined symmetry $\mathcal{C}_{2}\mathcal{T}$, the three-fold rotational symmetry $\mathcal{C}_{3}$ and a two-fold rotation with respect to the $x$-axis $\mathcal{M}_{x}$. As in monolayer graphene, the $\mathcal{C}_{3}$ symmetry guarantees Dirac points at both $\bm{\kappa}$ and $\bm{\kappa^{\prime}}$ in the mBZ. We find that near magic angle $\mathcal{C}_{3}$ is already broken at the weakest interactions we consider and that because of the flatness of the magic-angle bands, the Dirac point positions rapidly move close to $\gamma$ where the bands are most dispersive. The band topology evolution with $\theta$ and $\epsilon^{-1}$ is sensitive not only to interactions, but also to the details of the non-interacting band model (see Supplemental Material and XieMacDonald ). Unless otherwise specified, we have taken $T_{AA}/T_{AB}=0.8$, where $T_{AA}$ and $T_{AB}$ are the continuum model’s intra-sublattice and inter-sublattice hopping parameters, to account for corrugation and strain effects Kaxiras . Even though this choice gives a gap between flat bands and remote bands in the non-interacting limit, near the magic angle a small interaction strength is sufficient to pull the highest valence (lowest conduction) band down (up) in energy to touch the remote bands. The band topology at $\epsilon=5$ is illustrated in Fig. 2(b) by plotting Berry curvature in the $\mathcal{C}_{2}\mathcal{T}$ broken state as a function of moiré band momentum. Because of the involvement of remote bands, $\mathcal{C}_{2}\mathcal{T}$ symmetry no longer guarantees degeneracies between the first conduction and valence bands VishwanathWannierObstructions ; VishwanathMarch . Breaking $\mathcal{C}_{2}\mathcal{T}$ symmetry does however lift degeneracies between flat and remote bands, and generates corresponding Berry curvature hot spots that are visible in Fig. 2(b). The difference in condensation energy between $\mathcal{C}_{2}\mathcal{T}$ breaking and preserving states is extremely small. It follows that near magic twist angles $\mathcal{C}_{2}\mathcal{T}$ symmetry breaking is not essential for gap formation at moderate and stronger interaction strengths.

We characterize states that do break $\mathcal{C}_{2}\mathcal{T}$ symmetry by performing a Pauli matrix expansion of mBZ averages of the non-local Fock exchange self-energy, defining

[TABLE]

where $A$ is the system area, $A_{m}$ is the moiré unit cell area, $(i,j)=0,\ldots,3$ are Pauli matrix labels, $(s^{\prime}s)=A,B$ are sublattice labels and $(l^{\prime}l)=t,b$ (top, bottom) are layer labels. (The role of self-energy terms that are off-diagonal in reciprocal lattice vector as $\mathcal{C}_{2}\mathcal{T}$-breaking order parameters is discussed in the Supplemental Material.) The slow fall-off of the order parameter’s reciprocal lattice vector expansion reflects the spatial scale of quasiparticle wavefunction variation within the moiré unit cell. The largest symmetry-breaking self-energies are proportional to $\sigma^{z}_{s^{\prime}s}$ and $\tau^{0}_{l^{\prime}l}$, i.e. they are layer independent mass terms that favor one sublattice over the other. (We have sought self-consistent solutions with large $\sigma^{z}_{s^{\prime}s}\tau^{z}_{l^{\prime}l}$, self-energies but find that they are not stable.) The Berry curvatures plotted in Fig.2(b) are large near $\bm{\gamma}$, not near $\bm{\kappa},\bm{\kappa^{\prime}}$, as they would be if the same self-energy were added to a weakly-coupled-layer band Hamiltonian at a larger twist angle. We find that the flat bands sometimes have non-zero Chern numbers, implying that quantized anomalous Hall effects can occur ZhangChern when band occupations are valley-dependent.

Flavor symmetry breaking.— Because interactions normally induce gaps between conduction and valence bands, which then remain relatively flat, states with spin/valley flavor dependent band occupancies can be insulating. For example at $n=+p/4$, a state with the first conduction band occupied for $p$ flavors and empty for the remaining flavors is metastable if the exchange-energy shift of the conduction band upon occupation $U_{X}$ exceeds the conduction band width. A rough estimate based on non-self-consistent Hartree-Fock calculations (see Supplemental Material) yields $U_{X}\simeq 250{\rm meV}/\epsilon$, increasing slowly with twist angle. The flat band width, on the other hand, increases rapidly when the magic angle is exceeded, so that insulating states are restricted to the immediate vicinity of the magic angle. Fig. 3 illustrates the quasiparticle bands that emerge from a typical fully self-consistent mean-field calculation for a broken flavor symmetry insulator. In mean-field theory, coupling between flavors occurs only through the Hartree potential which is absent at neutrality and attractive at $AA$ sites in the moiré pattern at negative band filling factors. It follows that the energies and wave functions of the Hartree-Fock quasiparticle states of one flavor depend only weakly on the band-fillings of the other flavors.

Discussion.— Guided by earlier work Jung2014 , we anticipate that there is generally an energetic preference (not captured in continuum models) for states in which opposite valleys are occupied equally. insulating states at $n=\pm 1/4$ and $\pm 3/4$, must however break valley symmetry and are likely to be maximally valley polarized, which implies quantization of the anomalous Hall effect. It therefore follows from our calculations that quantized anomalous Hall effects in graphene bilayers can occur at quarter band fillings, but that it may not either because $\mathcal{C}_{2}\mathcal{T}$ symmetry is not broken or because the Chern number happens to equal zero. Indeed there is evidence Efetov experimentally that some insulating states at the quarters are Chern insulators and some are not. Further experimental work that maps out how this behavior depends on twist angles and distances to gates will be necessary to make a detailed comparison with mean-filed theory. The important role of thermal and quantum fluctuations of the collective fields present in the insulating states will be discussed elsewhere.

Note added.—Three independent related papersAshvin19 ; ValleyWave ; Dai19 which appeared after submission of our work report related results from mean-field theory and are complimentary to this paper.

Appendix A Supplemental Material

A.1 Self-consistent Hartree-Fock theory

for twisted bilayer graphene

Our mean field calculation is based on the continuum model proposed in BMModel . Below we first explain the formalism for the flavorless case, i.e. the valley projected and spinless case, and then extend it to include both spin and valley degrees of freedom. We choose the convention that the top and the bottom layers are respectively rotated counter-clockwise and clockwise around the $\hat{z}$ axis by a small angle $\theta/2$. The relative displacement between layers prior to rotation does not affect the spectrum and is ignored in our tunneling Hamiltonian. The quasiparticle wave functions in this theory are four component spinors, $\psi_{\alpha,\bm{k}}(\bm{r})$, $\alpha=\{A1,\ B1,\ A2,\ B2\}$, where $A(B)$ in the first index specifies sublattice and $1(2)$ in the second index specifies layer. We solve the mean-field equations by expanding each component in a plane-wave basis consistent with the continuum model’s spatial periodicity, which matches that of the moiré pattern.

The spinless single particle Hamiltonian projected onto valley K takes the form:

[TABLE]

where $\hat{h}_{\pm\theta/2}$ are the Dirac Hamiltonians for isolated rotated graphene layers,

[TABLE]

$\theta_{\bar{\bm{k}}}$ is the orientation angle of momentum measured from the Dirac point $\bar{\bm{k}}=\bm{k}-\bm{K}_{\theta}$. ($\bm{K}_{\pm\theta/2}$ is the Dirac momentum of top(bottom) layer.) Due to the moiré periodic modulation, interlayer tunneling is accompanied by momentum boosts $\bm{q}_{j}=\bm{0}$, $\bm{b}_{1}$ or $\bm{b}_{2}$ for $j=0,1,2$ respectively. $\bm{b}_{1,2}=(\pm 1/2,\sqrt{3}/2)4\pi/(\sqrt{3}a_{M})$ are the basis vectors of moiré reciprocal lattice, where $a_{M}=a/(2\sin(\theta/2))$ is the lattice constant of moire pattern and $a$ the lattice constant of monolayer graphene. The spatial modulation of the tunneling Hamiltonian takes the form:

[TABLE]

where

[TABLE]

Note that we have chosen a different convention for $T_{j}$ than in Ref. BMModel, which sets the origin of real space at the center of the AA region. The magnitude of interlayer tunneling is taken to be $\omega_{1}\equiv T_{AB}=110meV$ and $\omega_{0}\equiv T_{AA}=0.8\omega_{1}$. Below $\bm{k}$ is understood to be restricted to the first moiré Brillouin zone (mBZ) ($\bm{k}\in$ mBZ). For each $\bm{k}$ we employ the plane-wave expansion basis $|\psi_{\alpha,\bm{G},\bm{k}}\rangle$ where $\bm{G}=m\bm{b}_{1}+n\bm{b}_{2}$ and $m,n$ are integers. The single-particle Hamiltonian $\mathcal{\hat{H}}^{\rm{K}}_{0}$ has both terms that are diagonal in reciprocal lattice vector and terms that are off-diagonal in reciprocal lattice vector.

We take electron-electron interactions into account using the self-consistent Hartree-Fock method. In a plane wave basis, the Hartree and Fock self-energies are

[TABLE]

and

[TABLE]

where $\delta\rho=\rho-\rho_{\rm{iso}}$ is the density matrix defined relative to that of isolated rotated graphene layers each filled up to the charge neutrality point. The density-matrix $\rho$ is defined in a plane wave basis as:

[TABLE]

where the summation is over filled bands. We have also used $\delta\rho_{\alpha\beta}(\bm{G})\equiv\sum_{\bm{k},\bm{G}_{1}}\delta\rho_{\alpha,\bm{G}_{1}+\bm{G};\beta,\bm{G}_{1}}(\bm{k})$ in Eq. S5. $z_{\alpha,\bm{G},\bm{k}}^{n}$ is a numerical eigenvector in the plane-wave expansion,

[TABLE]

obtained by diagonalizing the total Hamiltonian

[TABLE]

Although interactions alter the lowest conduction and valence bands most strongly, we find it is generally essential to include a relatively large number of remote conduction and valence bands. The numerical results presented in the MS used a reciprocal space cut-off that allowed a total of 148 bands. This reciprocal space cut-off is sufficient to achieve converged insulating state condensation energies for the range of $\epsilon$ values we have considered. Spin and valley flavor degree of freedom are easily restored simply by observing that the Hartree self-energy has contributions from all occupied quasiparticle states of all flavors whereas the exchange self-energy has contributions from occupied quasiparticles of the same flavor only. We note that the Fock self-energy is strongly non-local and that this non-locality plays an essential role in how interactions select the character of the insulating states. No local exchange-correlation approximation can describe the interaction physics of twisted bilayers reliably.

A.2 Decomposition of the Fock self-energy and symmetry analysis

The Fock self-energy dominates the interaction physics in the SCHF Hamiltonian because it is responsible for the broken symmetries of TBG. In the following we discuss the expansion of the Fock self-energy in terms of layer and sublattice Pauli matrices. For each flavor we can classify terms in this expansion according to whether or not they violate the $\mathcal{C}_{2}\mathcal{T}$ symmetry. $\mathcal{C}_{2}\mathcal{T}$ symmetry is a property of the non-interacting Hamiltonian and must be broken in order to generate quasiparticle bands with non-zero momentum-space Berry curvatures and Chern numbers.

We can quite generally write:

[TABLE]

For $\bm{G}^{\prime}=\bm{G}$, $\Delta_{\bm{G},ij}\equiv\Delta_{\bm{G}\bm{G},ij}$ is always real as discussed in the main text. For $\bm{G}^{\prime}\neq\bm{G}$, $\Delta_{\bm{G}^{\prime}\bm{G},ij}$ is complex and, since $\Sigma^{F}$ is Hermitian, must satisfy $\Delta_{\bm{G}^{\prime}\bm{G},ij}=\Delta^{*}_{\bm{G}\bm{G}^{\prime},ij}$. If we write

[TABLE]

the $\Delta^{R}_{\bm{G}^{\prime}\bm{G},ij}$ terms have the same symmetry classification under $\mathcal{C}_{2}\mathcal{T}$ as the $\Delta_{\bm{G}\bm{G},ij}$ terms, i.e., states with finite $\Delta^{R}_{\bm{G}^{\prime}\bm{G},ij}$ where $(ij)\in\mathcal{S}\equiv\{(02),(12),(22),(30),(31),(33)\}$ break $\mathcal{C}_{2}\mathcal{T}$ symmetry. On the other hand, the $\Delta^{I}_{\bm{G}^{\prime}\bm{G},ij}$ terms have the opposite classification because of the complex conjugation operation in the time reversal operator such that terms $\Delta^{I}_{\bm{G}^{\prime}\bm{G},ij}$ with $(ij)$ outside $\mathcal{S}$ break $\mathcal{C}_{2}\mathcal{T}$ symmetry.

A.3 Band topology evolution with twist angle and interaction strength

Fig. S1a shows the non-interacting gap between the flat band and remote bands which depends on both twist angle $\theta$ and on the ratio $T_{AA}/T_{AB}$. At magic angle $\theta=1.1^{\circ}$ (marked by the vertical black dashed line), flat bands are in touch with remote bands in the original BM continuum model $T_{AA}/T_{AB}=1$. To approximate the effect of strain in interacting electron calculations we choose $T_{AA}/T_{AB}=0.8$ unless otherwise stated. In Fig. S1b we shows the energy separation between the highest interaction-modified quasiparticle valence band and the remote valence bands as a function of twist angle $\theta$ and the interaction strength $\epsilon^{-1}$.

In Fig. S2 we plot the momentum-space Berry curvature of the non-interacting continuum model bands with a small $\sigma_{3}\tau_{0}$ mass term added to break $\mathcal{C}_{2}\mathcal{T}$ symmetry and reveal the Dirac points. Near the magic angle the band topology is altered relative to the band topology of weakly coupled layers at larger twist angles by transferal of Dirac points from remote bands to the lowest conduction and valence bands. The band topology is further altered by electron-electron interactions as illustrated for the self-consistent mean field ground state discussed in the main text (Fig. 2(b)). For a twist angle just above the magic angle Fig. S2 shows that the weak-coupling $\pm\pi$ Berry phases near $\kappa$ and $\kappa^{\prime}$ are surrounded by three new Dirac points with the opposite chirality bands. The change in topology occurs at a slightly higher twist angle than illustrated here and is accompanied by band touchings involving remote bands that occur at momenta near $\gamma$. A detailed discussion of the dependence of non-interacting band topology on twist angle and on details of the interlayer hopping model will be presented elsewhere XieMacDonald .

A.4 Exchange-energy shift and nature of quasiparticle bandgap

Exchange interaction shifts the quasiparticle bands of a particular flavor downward in energy as the electron filling of this flavor increases. The quasiparticle bandgap of the flavor symmetry breaking states thus forms between the same active bands of different flavors, provided that the renormalized bandwidth is less than the amount of exchange-energy shift. Here we estimate the exchange-energy shift in a non-self-consistent (NSC) manner. As an example, we examine the case $p=1$, i.e., with the lowest conduction band of one particular flavor filled and the rest of the flavors filled to neutrality.

Define the amount of energy lowered by exchange self-energy from the filled flat conduction band as

[TABLE]

where $\psi^{\textrm{NSC}}_{c,\bm{k}}$ is the lowest conduction band wavefunction considering the self-energy of the frozen remote bands and $\delta\rho^{\textrm{NSC}}_{c}$ the corresponding density matrix. $\Sigma^{F}$ is the Fock self-energy defined in Eq. S6. Fig. S3 compares $U_{X}$ and the renormalized bandwidth $W_{c}$ of the filled conduction band as twist angle is varied. It captures the overall trend of the self-consistent value for $U_{X}-W_{c}$ that peaks near to but above the magic angle.

The exact value of the insulating gap depends on the details of the quasiparticle dispersion and must be determined by a converged self-consistent calculation. Fig. S4 shows quasiparticle band dispersions near band extrema for different flavor symmetry breaking states. For $n=-1/4$ and $-1/2$, the bandgap occurs between flat valence bands of filled and empty flavors and is indirect in momentum. However, for $n=-3/4$, because the remote valence band is slightly higher than the filled flat valence band, the bandgap becomes direct.

A.5 Quasiparticle dispersion of $n=-1/4$ and valley polarized $n=-1/2$ states

In this section, we present quasiparticle dispersions for SCHF ground states at filling factors $n=-1/4$ and $n=-1/2$ with valley polarization. Fig. S5 shows the overall quasiparticle dispersion for the two cases where the same set of parameters are used as in Fig. 3 (a) and (b). Both states are insulating with energy gap $E_{g}=3.16$ meV in the valley polarized $n=-1/2$ state and $E_{g}=0.06$ meV in the $n=-1/4$ state.

At $n=-1/2$, mean field states with full filling of either two of the four flavors, e.g., the valley polarized and the spin polarized states, are metastable and have degenerate total energies. We expect that including valley-anisotropic terms, such as short range Coulomb interactions and electron-phonon interactions, will lift the degeneracy and lead to ground states with inter-valley coherence. A detailed study of these interaction terms will be included in our future work.

A.6 Total energy difference $\mathcal{C}_{2}\mathcal{T}$ preserving and breaking states at neutrality

The self-consistent Hartree-Fock calculation yields solutions either preserve or break $\mathcal{C}_{2}\mathcal{T}$ symmetry. In Fig. S6, we present the total energy difference between the lowest energy $\mathcal{C}_{2}\mathcal{T}$ preserving solution and the lowest energy solution as a function of interaction strength for the magic angle case $\theta=1.1^{\circ}$, and for the large twist angle case $\theta=1.6^{\circ}$. In the range of interaction strength where the energy difference vanishes, we only found $\mathcal{C}_{2}\mathcal{T}$ preserving solutions. Near the magic angle, the energy gained tends to be small even when $\mathcal{C}_{2}\mathcal{T}$ symmetry is broken. For $\theta=1.6^{\circ}$. starting from small interaction strength limit, there is a phase transition from a gapless state with $\mathcal{C}_{2}\mathcal{T}$ symmetry to a gapped $\mathcal{C}_{2}\mathcal{T}$ breaking state which occurs at about $\epsilon^{-1}\approx 0.06$. For magic angle case, this transition occurs first at a very small interaction strength which is not captured by our discrete plot, but does not produce a large ground state energy reduction. At twist angle $\theta=1.6^{\circ}$, the cusp of the curve after the transition is believed to be associated with the complicated band crossing occurring between flat and remote bands.

A.7 Flat band projected Hartree and Fock self-energy

To have a better understanding of the different roles played by Hartree and Fock interaction self-energies on the low energy flat bands, we plot the projection of the self-energies to the Hilbert space of the non-interacting flat bands. We note this is done in a non-self-consistent fashion in which the self-energy are constructed from non-interacting density matrices. Fig. S7 plots the different components of the projected self-energies at neutrality and in the case in which the flat bands are all empty. At neutrality, the Hartree self-energy is vanishingly small and almost does not affect the flat band dispersion. However, the Fock self-energy plays an important role in changing the flat band states by lowering the energy of the filled valence band more strongly relative to the energy of the conduction band. In contrast, at empty filling of the flat valence bands, Hartree energy becomes finite. Because the Hartree self-energy is common to all the flavors, it does not play a role in splitting the energy between filled and empty flavors. The Fock self-energy is diagonal in flavor and differs between filled and empty flavors which result in exchange-energy shift responsible for flavor symmetry breaking as discussed in previous section. It does however alter the quasiparticle dispersions by lowering the energies of states that are more strongly localized near AA-stacking positions in the moiré pattern relative to states that are less strongly peaked.

A.8 Non-interacting eigenstates at the high symmetry momenta $\bm{\gamma}$ and $\bm{m}$

A minimal non-interacting Hamiltonian at the $\bm{m}$ point can be constructed by keeping the lowest order terms in Eq. S1:

[TABLE]

This Hamiltonian acts on the four-component spinor $(\psi_{A1,\bm{k}_{\bm{m}}},\psi_{B1,\bm{k}_{\bm{m}}},\psi_{A2,\bm{k}_{\bm{m}}},\psi_{B2,\bm{k}_{\bm{m}}})^{T}$ where $\bm{k}_{\bm{m}}$ is the momentum of $\bm{m}$. $\mathcal{H}_{0}^{K}(\bm{m})$ can be rewritten in a block diagonalized form by simultaneous rotation around the $y$ axis in both sublattice and layer pseudo-spin spaces so that $\sigma_{x}\rightarrow\sigma_{z}$, $\sigma_{z}\rightarrow-\sigma_{x}$ ,$\tau_{x}\rightarrow\tau_{z}$, and $\tau_{z}\rightarrow-\tau_{x}$. We then obtain the eigenstate energies,

[TABLE]

The two states close in energy to the charge neutrality point have energies $\pm(\omega_{1}-\sqrt{\hbar^{2}v_{D}^{2}k_{\theta}^{2}/4+\omega_{0}^{2}})$ which suggests that a band inversion can happen at $\bm{m}$ as twist angle $\theta$ is decreased when $\omega_{1}<\omega_{0}$.

At the high symmetry momentum $\bm{\gamma}$, the lowest order Hamiltonian couples degenerate Dirac states of the isolated layers with momentum $\Gamma^{t}_{i}=\{\bm{\gamma},\bm{\gamma}+\bm{b}_{1},\bm{\gamma}+\bm{b}_{1}+\bm{b}_{2}\}$ from the top layer and $\Gamma^{b}_{i}=\{\bm{\gamma},\bm{\gamma}+\bm{b}_{2},\bm{\gamma}+\bm{b}_{1}+\bm{b}_{2}\}$ from the bottom layer, where $i=0,1,2$. States at these momenta are coupled in a cyclic order and the resulting Hamiltonian takes the form,

[TABLE]

which acts on a twelve-component spinor (two sublattices at each of the six momenta.) We find that $H_{0}(\bm{\gamma})$ can be diagonalized by first rotating from the sublattice basis to the basis which diagonalize the Dirac Hamiltonian $h_{\pm\theta/2}$ at each of the six momenta. One can group the positive ($+$) and negative($-$) energy Dirac states and further rotate within each group to their eigen-basis. The final Hamiltonian can be written as,

[TABLE]

where

[TABLE]

The off-diagonal block takes the form

[TABLE]

Here $|\bm{k}_{\bm{\gamma}}|=k_{\theta}$ is the distance from either of the two Dirac points to $\bm{\gamma}$. We have defined $a=\frac{\sqrt{3}\omega_{0}}{2\hbar v_{D}|\bm{k}_{\bm{\gamma}}|}$ and $b=\frac{\omega_{1}}{\hbar v_{D}|\bm{k}_{\bm{\gamma}}|}$. The Hamiltonian can then be diagonalized and has eigenenergies that are either singly degenerate with $E=\pm\hbar v_{D}|\bm{k}_{\bm{\gamma}}|\cdot(\sqrt{1+4a^{2}/3}\pm 2b)$ or doubly degenerate with $E=\pm\hbar v_{D}|\bm{k}_{\bm{\gamma}}|\cdot(\sqrt{1+a^{2}/3}\pm\sqrt{3a^{2}+b^{2}})$. The degeneracy of these energies agrees with dimension of the irreducible representations of the little group at $\bm{\gamma}$.

The energies of the flat bands at the $\bm{\gamma}$ point for weakly coupled bilayers are singly degenerate and correspond to $\pm\hbar v_{D}|\bm{k}_{\bm{\gamma}}|\cdot(\sqrt{1+4a^{2}/3}-2b)$. For strong interlayer coupling or small twist angles, the doubly degenerate states become the states that are lowest in energy at $\bm{\gamma}$ and the flat bands will be in contact with the remote bands at $\gamma$ (see Fig. S1(a)). This delicate evolution of band topology with twist angle and hopping parameter ration is partially illustrated in Fig. S2.