Nematic topological semimetal and insulator in magic angle bilayer graphene at charge neutrality

TL;DR

This paper uses self-consistent Hartree-Fock calculations to explore interaction-driven phases in twisted bilayer graphene at charge neutrality, revealing nematic semimetals and insulators with topological properties influenced by strain.

Contribution

It identifies and characterizes multiple competing interaction-driven phases in twisted bilayer graphene, including topologically nontrivial semimetals and insulators, under realistic conditions.

Findings

Nematic topological semimetals are stabilized by weak strain.

Valley Chern insulators and symmetry-breaking insulators are identified.

Sample variability may explain different experimental observations.

Abstract

We report on a fully self-consistent Hartree-Fock calculation of interaction effects on the Moir\'e flat bands of twisted bilayer graphene, assuming that valley U(1) symmetry is respected. We use realistic band structures and interactions and focus on the charge neutrality point, where experiments have variously reported either insulating or semimetallic behavior. Restricting the search to orders for which the valley U(1) symmetry remains unbroken, we find three types of self-consistent solutions with competitive ground state energy (i) insulators that break $C_2 {\mathcal T}$ symmetry, including valley Chern insulators (ii) spin or valley polarized insulators and (iii) rotation $C_3$ symmetry breaking semimetals whose gaplessness is protected by the topology of the Moir\'e flat bands. We find that the relative stability of these states can be tuned by weak strains that break $C_3$ rotation. The nematic semimetal and also, somewhat unexpectedly, the $C_2 {\mathcal T}$ breaking insulators, are stabilized by weak strain. These ground states may be related to the semi-metallic and insulating behaviors seen at charge neutrality, and the sample variability of their observation. We also compare with the results of STM measurements near charge neutrality.

Equations

$$\begin{gathered} {\mathcal{H}}_{+}=\sum _l\sum _{\mathbf{k}}{f}_l^{†}(\mathbf{k})h_{\mathbf{k}}(l\theta /2)f_l(\mathbf{k}) \\ +\left(\sum _{\mathbf{k}}\sum _{i=1}^3{f}_t^{†}(\mathbf{k}+{\mathbf{q}}_i)T_i{f}_b(\mathbf{k})+h.c.\right). \end{gathered}$$

$$T_j=\left(\begin{array}{cc}{w}_0& w_1{e}^{-(j-1)\frac{2\pi i}{3}}\\ w_1{e}^{(j-1)\frac{2\pi i}{3}}& w_0\end{array}\right),$$

$$P_{\alpha ,\beta }(\mathbf{k})=⟨c_{\alpha }^{†}(\mathbf{k})c_{\beta }(\mathbf{k})⟩,P(\mathbf{k}{)}^2=P(\mathbf{k})=P(\mathbf{k}{)}^{†},$$

$$\begin{gathered} {\mathcal{H}}_{\mathrm{MF}}=\sum _{\mathbf{k}}\{c(\mathbf{k}{)}^{†}[h_0(\mathbf{k})+h_{\mathrm{HF}}(P,\mathbf{k})]c(\mathbf{k}) \\ -\frac{1}{2}\mathrm{tr}{h}_{\mathrm{HF}}(P,\mathbf{k})P^T(\mathbf{k})\}. \end{gathered}$$

$$\begin{gathered} h_{\mathrm{HF}}(P,\mathbf{k})=\frac{1}{A}\sum _{\mathbf{G}}{V}_{\mathbf{G}}{\Lambda }_{\mathbf{G}}(\mathbf{k})\sum _{{\mathbf{k}}^{\prime }}\mathrm{tr}{P}^T({\mathbf{k}}^{\prime }){\Lambda }_{\mathbf{G}}^{†}({\mathbf{k}}^{\prime }) \\ -\frac{1}{A}\sum _{\mathbf{q}}{V}_{\mathbf{q}}{\Lambda }_{\mathbf{q}}(\mathbf{k})P^T(\mathbf{k}+\mathbf{q}){\Lambda }_{\mathbf{q}}^{†}(\mathbf{k}), \end{gathered}$$

$$[{\Lambda }_{\mathbf{q}}(\mathbf{k}){]}_{\alpha ,\beta }=⟨u_{\alpha }(\mathbf{k})|u_{\beta }(\mathbf{k}+\mathbf{q})⟩.$$

$$h_0(\mathbf{k})=h_{\mathrm{BM}}(\mathbf{k})-h_{\mathrm{HF}}(P_0,\mathbf{k})$$

$${\chi }_{C_3}=\frac{1}{N}\sum _{\mathbf{k}}(1-|⟨{\psi }_{O_3\mathbf{k}}|C_3|{\psi }_{\mathbf{k}}⟩{|}^2),$$

$${\Lambda }_{\mathbf{q}}(\mathbf{k})=F_{\mathbf{q}}(\mathbf{k})e^{i{\Phi }_{\mathbf{q}}(\mathbf{k}){\sigma }_z{\tau }_z},$$

$$P(\mathbf{k})=\frac{1}{2}[1+Q(\mathbf{k})],Q(\mathbf{k}{)}^2=1,\mathrm{tr}Q(\mathbf{k})=0$$

$$E_F=-\frac{1}{8A}\sum _{\mathbf{k},\mathbf{q}}{V}_{\mathbf{q}}\mathrm{tr}Q(\mathbf{k}){\Lambda }_{\mathbf{q}}^{†}(\mathbf{k})Q(\mathbf{k}+\mathbf{q}){\Lambda }_{\mathbf{q}}(\mathbf{k}).$$

$$E_F\ge E_{F,\min }=-\frac{1}{A}\sum _{\mathbf{k},\mathbf{q}}{V}_{\mathbf{q}}{F}_{\mathbf{q}}^2(\mathbf{k}).$$

$$Q(\mathbf{k}+\mathbf{q})=e^{i{\Phi }_{\mathbf{q}}(\mathbf{k}){\sigma }_z{\tau }_z}Q(\mathbf{k})e^{-i{\Phi }_{\mathbf{q}}(\mathbf{k}){\sigma }_z{\tau }_z},$$

$$Q_{\mathrm{SM}}(\mathbf{k})={\sigma }_x{e}^{i{\alpha }_{\mathbf{k}}{\sigma }_z{\tau }_z}.$$

$$E_F=-\frac{1}{A}\sum _{\mathbf{k},\mathbf{q}}{V}_{\mathbf{q}}{F}_{\mathbf{q}}^2(\mathbf{k})\cos [{\alpha }_{\mathbf{k}}-{\alpha }_{\mathbf{k}+\mathbf{q}}+2{\Phi }_{\mathbf{q}}(\mathbf{k})].$$

$$E_F=E_{F,\min }+\frac{1}{2A}\sum _{\mathbf{k},\mathbf{q}}{V}_{\mathbf{q}}{\mathbf{q}}^2{F}_{\mathbf{q}}^2(\mathbf{k})({\nabla }_{\mathbf{k}}{\alpha }_{\mathbf{k}}-2{\mathbf{A}}_{\mathbf{k}}{)}^2.$$

$$h_0(\mathbf{k})={\sigma }_{x}f(\mathbf{k})e^{i\varphi (\mathbf{k}){\sigma }_z{\tau }_z}$$

$$E_h[Q]=\frac{1}{2}\sum _{\mathbf{k}}\mathrm{tr}{h}_0(\mathbf{k})Q(\mathbf{k})$$

$$E[Q]=E_h[Q]+E_F[Q]$$

$$[Q(\mathbf{k}),\frac{1}{2}{h}_0(\mathbf{k})-\frac{1}{4A}\sum _{\mathbf{q}}{V}_{\mathbf{q}}(\mathbf{k}){\Lambda }_{\mathbf{q}}(\mathbf{k}{)}^{\ast }Q(\mathbf{k}+\mathbf{q}){\Lambda }_{\mathbf{q}}(\mathbf{k}{)}^T]=0$$

$$Q(\mathbf{k})=Q_{C_2\mathcal{T}\mathrm{I}}(\mathbf{k})\cos \theta (\mathbf{k})+Q_{\mathrm{SM}}(\mathbf{k})\sin \theta (\mathbf{k})$$

$$|\sin \theta |=\frac{E_h[Q_{\mathrm{SM}}]}{2(E_F[Q_{\mathrm{SM}}]-E_{F,\min })}$$

$$E[Q]=E_{F,\min }-\frac{E_h[Q_{\mathrm{SM}}{]}^2}{4(E_F[Q_{\mathrm{SM}}]-E_{F,\min })}$$

$$\nu (\mathbf{r})=8\left(\frac{{\rho }_{\mathrm{LB}}(\mathbf{r})}{{\rho }_{\mathrm{LB}}(\mathbf{r})+{\rho }_{\mathrm{UB}}(\mathbf{r})}-\frac{1}{2}\right)\in [-4,4],$$

$${\mathcal{H}}_{+}=\sum _l\sum _{\mathbf{k}}{\tilde{f}}_l^{†}(\mathbf{k})h_{\mathbf{k}}(l\theta /2){\tilde{f}}_l(\mathbf{k})+\left(\sum _{\mathbf{k}}\sum _{i=1}^3{\tilde{f}}_t^{†}(\mathbf{k}+{\mathbf{q}}_i)T_i{\tilde{f}}_b(\mathbf{k})+h.c.\right).$$

$$h_{\mathbf{k}}(\theta )=\hslash v_F\left(\begin{array}{cc}0& (k_x-i{k}_y)e^{i\theta }\\ (k_x+i{k}_y)e^{-i\theta }& 0\end{array}\right),$$

$$T_1=(1+\beta )\left(\begin{array}{cc}{w}_0& w_1\\ w_1& w_0\end{array}\right),T_2=\left(\begin{array}{cc}{w}_0& w_1{e}^{-2\pi i/3}\\ w_1{e}^{2\pi i/3}& w_0\end{array}\right),T_3=\left(\begin{array}{cc}{w}_0& w_1{e}^{2\pi i/3}\\ w_1{e}^{-2\pi i/3}& w_0\end{array}\right),$$

$$C_3{\tilde{f}}_{\mathbf{k}}{C}_3^{-1}=e^{-\frac{2\pi }{3}i{\tau }_z{\sigma }_z}{\tilde{f}}_{C_3\mathbf{k}},(C_2\mathcal{T}){\tilde{f}}_{\mathbf{k}}(C_2\mathcal{T}{)}^{-1}={\sigma }_x{\tilde{f}}_{\mathbf{k}},{\mathcal{M}}_y{\tilde{f}}_{\mathbf{k}}{\mathcal{M}}_y^{-1}={\sigma }_x{\mu }_x{\tilde{f}}_{M_y\mathbf{k}},$$

$$\mathcal{T}{\tilde{f}}_{\mathbf{k}}{\mathcal{T}}^{-1}={\tau }_x{\tilde{f}}_{-\mathbf{k}}.$$

$$c_{\mu ,n}^{†}(\mathbf{k})=\sum _{\mathbf{G},I}{u}_{\tau ,n;\mathbf{G},I}(\mathbf{k})f_{\mu ,I}^{†}(\mathbf{k}+\mathbf{G}),$$

Full text

I Introduction

The discovery of interaction-driven insulating and superconducting behavior in twisted bilayer graphene (TBG) Cao et al. (2018a, b) has inspired intensive efforts to understand this behavior Xu and Balents (2018); Po et al. (2018a); Isobe et al. (2018); Thomson et al. (2018); You and Vishwanath (2018); Kang and Vafek (2018); Xie and MacDonald (2018); Lin and Nandkishore (2019); Dodaro et al. (2018); Padhi et al. (2018); Wu et al. (2018); Lian et al. (2018a); Zou et al. (2018); Roy and Juričić (2019) and to find related systems which exhibit similar phenomenology Zhang et al. (2019a); Khalaf et al. (2019); Mora et al. (2019); Cea et al. (2019); Bi et al. (2019); Lee et al. (2019). This work has started to bear fruit with several groups announcing similar observations in TBG samples Yankowitz et al. (2019); Sharpe et al. (2019); Lu et al. (2019); Codecido et al. (2019) as well as other Moiré materials Chen et al. (2019a, b); Liu et al. (2019); Shen et al. (2019); Cao et al. (2019). The basic mechanism underlying the enhancement of correlation in these materials is understood to originate from the long-wavelength Moiré pattern leading to quenching of the electron kinetic energy manifested in flat energy bands Bistritzer and MacDonald (2011); Lopes dos Santos et al. (2012). Nevertheless, the nature of the observed correlated insulating states remains under debate Po et al. (2018a); Thomson et al. (2018); Isobe et al. (2018); Kang and Vafek (2018); Xie and MacDonald (2018).

Early experiments found clear signature of a correlated insulating state at half-filling 111We follow the standard convention of measuring the filling relative to charge neutrality where the Dirac points of the original graphene sheet reside. Complete filling corresponds to $\nu_{T}=4$ while completely empty $\nu_{T}=-4$ electrons per Moiré unit cell, accounting for both spin and valley degeneracy. Half filling corresponds to $\nu_{T}=\pm 2$ Cao et al. (2018a); Yankowitz et al. (2019). Subsequently, insulating ferromagnetic states were also observed at quarter and three-quarter fillings Yankowitz et al. (2019); Sharpe et al. (2019). On the other hand, in all these experimentsCao et al. (2018a, b); Yankowitz et al. (2019) insulating behavior was absent at charge neutrality (CN) where signatures of semimetallic behavior were observed instead. In contrast, a recent experiment surprisingly found an insulator at CN with a transport gap exceeding those at 1/2, 1/4 and 3/4 fillings Lu et al. (2019).

On the theory side, it was realized early on Po et al. (2018a, b) that a simple Mott picture for the insulating phase is complicated by the band topology which prohibits the construction of localized orbitals describing the flat bands while preserving all the symmetries. Various orders have been proposed to account for the insulating states. At charge neutrality, a $C_{2}\mathcal{T}$ symmetry breaking insulator, with Chern number $\pm 1$ for each spin and valley flavour was proposed in Xie and MacDonald (2018), along with a $C_{2}\mathcal{T}$ symmetry preserving insulator that is believed to require mixing with remote bands Xie and MacDonald (2018). An intervalley coherent order Po et al. (2018a) was proposed as a candidate for the insulating state at half-filling while nematic orders were discussed in Dodaro et al. (2018); Thomson et al. (2018); Sboychakov et al. (2019) and ferromagnetic ordering was proposed in Kang and Vafek (2018). In the presence of explicit $C_{2}$ symmetry breaking induced by a substrate, valley or spin polarized insulator with valley resolved Chern numbers have been discussed Zhang et al. (2019a); Bultinck et al. (2019a).

In this letter, we perform a self-consistent Hartree-Fock mean field analysis for the screened Coulomb interaction projected onto the flat bands.

We focus on the CN point because of its pivotal role in determining the entire phase diagram. We will discuss other fillings in subsequent work. For simplicity, we restrict our attention to orders that conserve the U(1) valley charge as well as translation invariance at the scale of the Moiré unit cell. Our results include the expected spin-polarized and valley-polarized insulators, which break no other symmetries. In addition we observe a strong tendency to breaking spatial rotation symmetries. We find a $C_{2}\mathcal{T}$-breaking insulator and two distinct $C_{2}\mathcal{T}$-symmetric semimetallic phases which break $C_{3}$-symmetry . The $C_{2}\mathcal{T}$-breaking phase has a Chern number of $C=\pm 1$, per flavor (i.e. valley and spin). Different spin/valley orderings then lead to various ground states ranging from quantum anomalous Hall (QAH) to quantum valley Hall (QVH) or quantum spin Hall (QSH) insulators with very similar Hartree-Fock energies. Of these, the QVH insulator breaks $C_{2}$ symmetry, allowing for a direct coupling to the $C_{2}$-breaking hBN substrate, potentially favoring it in situations where the sample and substrate are aligned. The gapless $C_{3}$-breaking phases are obtained by bringing the two Dirac cones from the Moiré $K$ and $K^{\prime}$ very close to the $\Gamma$ point. Instead of merging and opening a gap as one might normally expect, the Dirac points remain gapless since they carry the same chirality Po et al. (2018a), a consequence of descending from the Dirac points of graphene from the same valley for the two layers. This topological protection prevents them from annihilating, resulting in a gapless semimetallic state. Thus, the metallic nature at CN in this scenario is intimately tied to the topological properties of the magic angle flat bands.

We investigate the effect of small explicit $C_{3}$ symmetry breaking which can arise in real samples due to strain Bi et al. (2019), and show that it strongly influences the competition between different symmetry broken phases, favoring one of the $C_{3}$-breaking semimetallic phases. Surprisingly, the $C_{2}\mathcal{T}$-breaking insulators also exhibits a strong susceptibility to $C_{3}$ symmetry breaking. The energies of the insulating $C_{2}\mathcal{T}$-breaking and the semimetallic $C_{2}\mathcal{T}$-preserving states are lowered compared to the spin/valley ferromagnets, and approach each other quickly as the value of the $C_{3}$-breaking parameter is increased. Our results suggest that these two states are candidate ground states in the presence of very small explicit $C_{3}$ symmetry breaking which is likely to exist in experiments. Further competition between these two phases is likely to be settled by small sample-dependent details, potentially explaining the realization of an insulator in some samples and a semimetal in other samples.

II Problem setup

The single-particle physics is described by the Bistritzer-MacDonald (BM) model Lopes dos Santos et al. (2012); Bistritzer and MacDonald (2011), which employs a continuum approximation close to $K$ and $K^{\prime}$ for a pair of graphene sheets rotated relative to each other by an angle $\theta$. The Hamiltonian for the $K$ valley is given by:

[TABLE]

Here, $l=t/b\simeq\pm 1$ is the layer index, and $f_{l}({\bm{k}})$ is the $K$-valley electron originated from layer $l$. $h_{{\bm{k}}}(\theta)$ is the monolayer graphene $K$-valley Hamiltonian with twist angle $\theta$ (see Appendix for details) and ${\bm{q}}_{1}$ is defined as $K_{b}-K_{t}$ with $K_{l}$ denoting the $K$-vector of layer $l$. ${\bm{q}}_{2}=O_{3}{\bm{q}}_{1}$ is the counterclockwise $2\pi/3$ rotation of ${\bm{q}}_{1}$, and ${\bm{q}}_{3}=O_{3}{\bm{q}}_{2}$. Finally, the interlayer coupling matrices are given by

[TABLE]

with $w_{0}$ and $w_{1}$ denoting intrasublattice and intersublattice hopping, respectively. Due to lattice relaxation effects, which shrink the AA regions relative to the AB regions, the value of $w_{0}$ at the magic angle is about 75% of $w_{1}$ Nam and Koshino (2017); Koshino et al. (2018). Throughout this work, we will use the values $w_{1}=110$ meV and $w_{0}=82.5$ meV. Explicit $C_{3}$ symmetry breaking is implemented via the substitution $T_{1}\rightarrow(1+\beta)T_{1}$ Zhang et al. (2019b).

To study possible correlated insulating states, we employ a momentum-space self-consistent Hartree-Fock mean field theory. The momentum-space description allows us to focus on the pair of flat bands at CN, thereby evading the difficulties associated with the real space Wannier obstruction rooted in the fragile topology of these bands Po et al. (2018c, a, b).Restricting the analysis to the flat bands is only justified in the limit when both the bandwidth and interaction strength are much smaller than the bandgap. For realistic interactions, the interaction strength is of the same order as the bandgap which means that symmetry-broken states involving remote bands cannot be ruled out Xie and MacDonald (2018). However, given the strong angle dependence and appearance of correlated states only near the magic angle Cao et al. (2018a, b); Yankowitz et al. (2019), it is likely that the relevant symmetry-broken phases experimentally originate mostly from the two flat bands which justifies our approximation. We note, however, that our approach cannot capture possible changes in the topology of the flat bands arising from mixing with remote bands which were argued to be relevant in Ref. Xie and MacDonald (2018).

The Hartree-Fock (HF) mean field theory is defined in terms of the projector

[TABLE]

where $c_{\alpha}({\bm{k}})$ is the annihilation operator for an electron at momentum ${\bm{k}}$ and $\alpha=(n,\tau,s)$ is a combined index for band, valley and spin, respectively. For a gapped or semimetallic phase at CN, $P$ satisfies $\mathop{\mathrm{tr}}P({\bm{k}})=4$ for all ${\bm{k}}$ points222For semimetals, this is true everywhere except for the gapless points where the projector is not defined. The HF mean field Hamiltonian has the form

[TABLE]

Here, $c({\bm{k}})$ is a column vector in the index $\alpha$, $h_{0}({\bm{k}})$ denotes the single particle Hamiltonian and $h_{\rm HF}({\bm{k}})$ is given by

[TABLE]

where $A$ is the total area. The ${\bm{k}}^{\prime}$, ${\bm{G}}$, and ${\bm{q}}$ summations range over the first Brillouin zone, reciprocal lattice vectors, and all momenta, respectively. $V_{\bm{q}}$ is the interaction potential which we take to be a single-gate-screened Coulomb interaction $V_{\bm{q}}=\frac{e^{2}}{2\epsilon\epsilon_{0}q}(1-e^{-2qd_{s}})$ with dielectric constant $\epsilon=7$ and screening length equal to the gate distance $d_{s}\approx 40$ nm.

The first term in (5) is the Hartree term while the second is the Fock term. The matrix $\Lambda_{\bm{q}}({\bm{k}})$ contains the form factors for the single-particle

[TABLE]

where $\alpha$ and $\beta$ denote a combined index for spin, valley and band. The HF self-consistent analysis starts by proposing an ansatz for the projector $P({\bm{k}})$, then substituting in the Hamiltonian (4) which is then used to compute the new projector. This procedure is iterated until convergence is achieved.

One important subtlety in the HF approach is that the band structure depends on the filling even without symmetry breaking. This follows from the fact that the form factor matrix $\Lambda_{\bm{q}}({\bm{k}})$ is not diagonal in the band index for ${\bm{q}}\neq 0$ since the Bloch wavefunctions $u_{\alpha,{\bm{k}}}$ from different bands are not orthogonal at different momenta. In addition, the Hartree term contains a trace over all filled bands which also affects the dispersion of the empty bands. This means that the band structure obtained from the BM model is only valid at a specific filling which determines the references point for out analysis. At this point, it will be assumed that interaction effects are already included in the parameters of the effective model which can be obtained by fitting to ab initio calculations or comparing to STM data away from the magic angles Choi et al. (2019); Jiang et al. (2019); Kerelsky et al. (2019). A natural choice of the reference point, which we adopt throughout this letter, is the CN point. This means that the single particle Hamiltonian $h_{0}({\bm{k}})$ is given by

[TABLE]

where $h_{\rm BM}({\bm{k}})$ is the BM Hamiltonian and $P_{0}$ is the projector corresponding to symmetry unbroken state obtained by filling the lower band of the BM Hamiltonian at CN. We note that in the different approach of Ref. Xie and MacDonald (2018) where multiple bands are included, the reference point was taken in the limit of decoupled layers. In our complementary approach of projection onto the flat bands it is not possible to implement such a choice.

Using CN as our reference point implies that the bands at empty or full filling should include some HF corrections leading to a modified band structure shown in Fig. 1 together with the resulting density of states (DOS). We notice that the separation between the two peaks in the DOS is about 10-15 meV in agreement with the measured DOS in STM experiments Choi et al. (2019, 2019); Jiang et al. (2019); Xie et al. (2019). Thus, our approach provides an explanation for the discrepancy between the experimentally measured peak separation and the expectation based on the BM model whose bandwidth close to the magic angle is much smaller (1-3 meV).

III Symmetry-broken phases

The interacting TBG Hamiltonian is characterized by the following symmetries Po et al. (2018a): spinless time-reversal $\mathcal{T}$ mapping the two valleys, SU(2) spin rotation in each valley, U(1) valley charge conservation and $C_{6}$ symmetry as well as a mirror symmetry which switches layers and sublattices but acts within each valley. Of these, only spin rotation, $C_{3}$, mirror and $C_{2}\mathcal{T}$ act within each valley. At integer filling, different correlated insulating phases can emerge by breaking some of these symmetries. Time-reversal symmetry is broken by valley polarized (VP) states, where the filling of the two valleys is different. Spin rotation symmetry is broken by spin polarization (SP) leading to ferromagnetic order. U(1) valley charge conservation is broken in the presence of intervalley coherent (IVC) superposition of states from the two valleys. $C_{2}$ symmetry is broken by sublattice polarization which gaps out the Dirac points at the Moiré $K$ and $K^{\prime}$ points.

Breaking $C_{3}$ symmetry alone does not generally lead to a gapped phase since it only moves the Dirac points away from the Moiré $K$ and $K^{\prime}$ without gapping them out. Even strong $C_{3}$ breaking does not result in the merging and gapping of the Dirac nodes, in contrast with other familiar band structures such as single layer graphene. This follows from the fact that the chirality of the Dirac nodes, which is well defined in the presence of $C_{2}{\mathcal{T}}$ symmetry, is the equal Po et al. (2018a); Zou et al. (2018), rather than opposite. Alternately, one can phrase the argument as follows. The topology of the two flat bands is captured by the second Stiefel-Whitney invariant $w_{2}$ Ahn et al. (2019); Po et al. (2018c); Lian et al. (2018b). This invariant is protected by $C_{2}\mathcal{T}$ and only depends on the flat band eigenstates which are unaffected by any symmetry breaking that does not involve other bands (which is the main assumption in this work). Indeed, the $w_{2}$ invariant must be trivial for a single isolated band, implying we cannot separate the pair of connected flat bands 333In general, the total second Stiefel-Whitney number of two isolated bands is not equal to the sum of the second Stiefel-Whitney numbers separately for the two bands; there is an additional contribution from the first Stiefel-Whitney classes. In the special case we considered here, the $C_{3}$ symmetry of the two-band subspace guarantees this additional term vanishes.. Thus the non-trivial $w_{2}$ invariant implies that the two Dirac points cannot be removed without breaking $C_{2}\mathcal{T}$.

Before presenting the numerical results, let us make the following observations. First, it is relatively easy to show that a state with uniform full spin or valley polarization is always a self-consistent solution to the HF equations at CN for sufficiently narrow bands. These two states have the same energy in the absence of intervalley Hund’s coupling Zhang et al. (2019a); Lee et al. (2019). The IVC state is known to have higher energy than SP/VP states for isolated bands with non-vanishing valley Chern number Zhang et al. (2019a); Bultinck et al. (2019a). These arguments do not generalize to TBG where the extra symmetries of the problem complicate the discussion (see Appendix and Ref. Bultinck et al. (2019b)). Nevertheless, we will exclude IVC orders from our numerical analysis (which is equivalent to assuming unbroken U(1) valley charge conservation) since it leads to significant simplification by allowing us to focus on a single flavor (spin and valley). Different diagonal spin-valley orders can then be generated from the single-flavor solution by applying different symmetries. The energy competition with U(1) valley symmetry broken phases is considered in Bultinck et al. (2019b).

IV Results

The results for the self-consistent HF analysis are provided in Fig. 2 showing the energies of the different solutions as a function of the $C_{3}$ symmetry breaking parameter $\beta$. There are two types of gapped solutions corresponding to either flavor-polarized (spin/valley) states or $C_{2}\mathcal{T}$-breaking ($C_{2}\mathcal{T}$I) insulators. There are three gapped solutions corresponding to spin-polarized (SP), valley-polarized (VP) and $C_{2}\mathcal{T}$-breaking ($C_{2}\mathcal{T}$I) insulators. The latter does not break $C_{3}$ when $\beta=0$ but develops a large $C_{3}$ breaking component for $\beta\neq 0$. The extent of $C_{3}$ symmetry breaking can be quantified by defining

[TABLE]

which vanishes for any $C_{3}$-symmetric state. Here, $|\psi_{{\bm{k}}}\rangle$ are the occupied single-particle eigenstates of the HF Hamiltonian with momentum ${\bm{k}}$ (for a given flavor). The value of $\chi_{C_{3}}$ for the $C_{2}\mathcal{T}$I state is shown in Fig. 2 as a function of $\chi_{C_{3}}$ for the corresponding non-interacting states arising from explicit $C_{3}$-breaking parameter $\beta\neq 0$ in the BM model. We can clearly see from the figure that a relatively small $\chi_{C_{3}}(BM)\sim 10^{-3}$ in the non-interacting states induces induces a much larger $C_{3}$ symmetry breaking of almost two orders of magnitude in the $C_{2}\mathcal{T}$I state. This serves to show that the $C_{2}\mathcal{T}$I has very large susceptibility to $C_{3}$ symmetry breaking. In the following, we will refer to this state for positive and negative $\beta$ as $C_{3}C_{2}\mathcal{T}$I$y$ and $C_{3}C_{2}\mathcal{T}$I$x$, respectively. In addition to this insulating state, there are two distinct $C_{2}$-preserving semimetallic phases which spontaneously break $C_{3}$ even for $\beta=0$ which we denote by $C_{3}$S$x$ and $C_{3}$S$y$, since they have Dirac points along $k_{x}$ and $k_{y}$, respectively. We notice that these ground states are similar to the ones obtained within a 10-band model Ref. Choi et al. (2019) which used a site-local ansatz for the interactions. Examples of the band structures for the $C_{3}$S$y$ and $C_{3}C_{2}\mathcal{T}$I$y$ states are shown in Fig. 3.

The $C_{2}\mathcal{T}$I phase obtained here is characterized by a Chern number of $\pm 1$ for a given spin and valley. The nature of the resulting phase depends on the precise symmetries which are broken as follows: (i) if $\mathcal{T}$ is broken but $C_{2}$ and spin rotation are preserved, we obtain a quantum anomalous Hall (QAH) insulator with total Chern number $C=\pm 4$, (ii) if spin rotation is broken we obtain a quantum-spin Hall (QSH) insulator with opposite Chern number $C_{s}=\pm 2$ for opposite spins, (iii) if $C_{2}$ is broken but $\mathcal{T}$ and spin rotation symmetry are preserved, then a quantum valley Hall (QVH) insulator obtains with opposite Chern number $C_{v}=\pm 2$ for opposite valleys, or (iv) we can additionally break spin rotation symmetry to obtain a state where opposite spins within the same valley also have opposite Chern numbers resulting in a quantum spin-valley Hall (QSVH) insulator where flipping either spin or valley flips the Chern number. There can also be states with total Chern number $\pm 2$, but these do not lead to any new physics or a lower energy, therefore we restrict to the above four possibilities for simplicity. The QH, QSH and QVSH preserve either $C_{2}$ or a combination of $C_{2}$ and some internal symmetries, thus the only state expected to couple to the $C_{2}$-breaking hBN substrate is the QVH which is likely to be energetically favored in aligned samples444We thank A. Thomson for correspondence on this point. The energies of the four possible $C_{2}\mathcal{T}$ states are very similar differing only by very small $\sim 10^{-3}$ meV Hartree corrections. Similarly, the flavor-polarized state can be spin-polarized, breaking SU(2) spin rotation, valley-polarized breaking $\mathcal{T}$ or a spin-valley-locked states which breaks both but preserves a combination of $\mathcal{T}$ and $\pi$-spin rotation. These states are degenerate on the mean field level and are only distinguished by intervalley Hund’s coupling Zhang et al. (2019a); Lee et al. (2019) neglected in this study.

At $\beta=0$, the energies of the insulating VP/SP and $C_{2}\mathcal{T}$I states are very close and smaller by about 1 meV per particle than the energies of the two semimetallic $C_{3}$ breaking states. An explanation for this fact is provided in the next section in the limit where the intrasublattice hopping is switched off Tarnopolsky et al. (2019). In this limit, we can establish a rigorous bound for the Fock energy which plays the dominant role in the energy competition. We will show that the SP/VP and the $C_{2}\mathcal{T}$I states saturate the bound in this limit and should thus be approximately degenerate. The energies of the $C_{3}$S states are higher than these two. This analysis explains why these states are expected to be close in energy, but it is generally insufficient to capture the details of the energy competition which depends sensitively on the intrasublattice hopping $w_{0}$ and can only be determined numerically. A more detailed analysis of the effects of dispersion and finite sublattice symmetry breaking $w_{0}$ is given in Ref. Bultinck et al. (2019b).

Once $\beta$ becomes non-zero, the energy of the $C_{2}\mathcal{T}$ breaking state is reduced relative to the VP/SP state. Furthermore, one of the two $C_{3}$-breaking states (depending on the sign of $\beta$) goes down in energy becoming more energetically favorable to the VP/SP state around $\beta=\pm 3\times 10^{-4}$. For larger values of $\beta\gtrsim 4\times 10^{-4}$, the energies of the insulating $C_{2}\mathcal{T}$-breaking phase and the semimetallic $C_{2}\mathcal{T}$-preserving phase approach the same value. This indicates that even very small explicit $C_{3}$-breaking picks out these two states as the main candidates for the ground state at CN. Assuming that such small explicit $C_{3}$ symmetry breaking exists in real samples due to strain, our analysis leads to the conclusion that the ground state of TBG at CN is either a $C_{2}\mathcal{T}$-breaking insulator or a $C_{2}\mathcal{T}$-symmetric semimetal. Both states strongly break $C_{3}$ symmetry and are very close in energy. It is worth noting that, in the presence of disorder, massless Dirac fermions may also emerge from the spatial domain walls of locally $C_{2}\mathcal{T}$ breaking insulating regions Thomson and Alicea (2019).

V Energy competition in chiral limit

The existence of several states close in energy can be explained by considering the simplified chiral model of Ref. Tarnopolsky et al. (2019) where the Moiré intrasublattice hopping is switched off. In such limit, the band becomes strictly flat at charge neutrality in the absence of interaction and the model has an extra sublattice symmetry i.e. the single-particle Hamiltonian anticommutes with the sublattice operator $\sigma_{z}$. As a result, the flatband wavefunctions can be chosen to be eigenfunctions of $\sigma_{z}$ such that the wavefunction $\sigma_{z}u_{A/B,\tau,{\bm{k}}}=\pm u_{A/B,\tau,{\bm{k}}}$. In this limit, we can show (see Appendix for details) that the form factor matrix has the simple form

[TABLE]

where $F_{\bm{q}}({\bm{k}})$ and $\Phi_{\bm{q}}({\bm{k}})$ denoting the magnitude and phase of the form factor for the flatband wavefunction in sublattice A and valley $K$.

Let us now compare the mean field energy for different phases. First, we notice that, even in the chiral limit where the non-interacting band dispersion vanishes, a sizable dispersion would be generated by the interaction Bultinck et al. (2019b). It will be instructive, however, to start by ignoring the interaction generated band dispersion and focusing on the Hartree-Fock energy. At a fixed integer filling, the Hartree term does not play a role in the energy competition between phases so we can focus on the Fock term. We will find it convenient to write the projector $P({\bm{k}})$ as

[TABLE]

In terms of $Q$, the Fock energy is given by

[TABLE]

up to a constant. We note that $\langle A,B\rangle=\mathop{\mathrm{tr}}AB$ defines a positive definite inner product on the space of hermitian matrices. Using Cauchy-Schwarz inequality, we get

[TABLE]

This inequality is satisfied if and only if $Q({\bm{k}}+{\bm{q}})$ is parallel to $\Lambda_{{\bm{q}}}({\bm{k}})Q({\bm{k}})\Lambda^{\dagger}_{{\bm{q}}}({\bm{k}})$ for every ${\bm{k}}$ and ${\bm{q}}$ which implies

[TABLE]

A ${\bm{k}}$-independent flavor (spin,valley) polarized state obviously satistfies this constraint since $\sigma_{z}\tau_{z}$ is diagonal in flavor, thereby saturating the Fock bound. This state can either be spin-polarized $Q=s_{z}$, valley polarized $Q=\tau_{z}$ or spin-valley locked $Q=\tau_{z}s_{z}$.

Next let us discuss $C_{2}\mathcal{T}$ breaking insulating solutions. A $C_{2}\mathcal{T}$ breaking solution to Eq. 13 is obtained by taking $Q_{\bm{k}}=\sigma_{z}$ (since $C_{2}\mathcal{T}$ is off diagonal in the sublattice index). The resulting state saturates the Fock bound and is characterized by a Chern number $\pm 1$ in the lower band. Its energy competition with the flavor-polarized states is settled by the effects of the single-particle term $h_{\bm{k}}$ and $w_{0}$ which are neglected here. In our numerics, we found that these states are very close in energy for the realistic model with their energy difference only sensitive to the strain parameter $\beta$. Depending on its structure in flavor space, the $C_{2}\mathcal{T}$-breaking insulator may correspond to one of four states: (i) $Q=s_{0}\tau_{0}\sigma_{z}$ breaks $C_{2}$ but not $\mathcal{T}$ and corresponds to a valley Hall state, (ii) $Q=s_{0}\tau_{z}\sigma_{z}$ breaks $\mathcal{T}$ but not $C_{2}$ and corresponds to a quantum Hall state with Chern number $\pm 4$, (iii) $Q=s_{z}\tau_{0}\sigma_{z}$ corresponds to a valley-spin-Hall state where the Chern number is invariant under flipping spin and valley, and (iv) $Q=s_{z}\tau_{z}\sigma_{z}$ corresponds to a spin-Hall state.

Finally, let us consider semimetallic states which break neither flavor nor $C_{2}\mathcal{T}$ symmetry. These can be generally described (within each flavor) by the order parameter

[TABLE]

Substituting in the Fock energy (11) yields

[TABLE]

To simplify further, we make the reasonable assumption that $F_{\bm{q}}({\bm{k}})$ decays quickly with the relative momentum ${\bm{q}}$ which enables us to expand $\Phi_{\bm{q}}({\bm{k}})$ and $\alpha_{{\bm{k}}+{\bm{q}}}$ in ${\bm{q}}$. Noting that $\Phi_{\bm{q}}({\bm{k}})={\bm{q}}\cdot{\bm{A}}_{\bm{k}}+O({\bm{q}}^{2})$ with ${\bm{A}}_{\bm{k}}$ the Berry connection $-i\langle u_{A,K,{\bm{k}}}|\nabla_{\bm{k}}|u_{A,K,{\bm{k}}}\rangle$, we get

[TABLE]

The second term in the energy is non-negative and implies that the semimetal does not generally satisfy the bound. In fact, this term has the form of the energy of a superconductor in a magnetic field in momentum space. Due to the non-trivial Chern number of the sublattice-polarized bands Tarnopolsky et al. (2019); Bultinck et al. (2019b), the Bloch wavefunctions cannot be chosen to be simultanuously smooth and periodic over the Brillouin zone. If we choose them to be smooth, the Berry connection ${\bm{A}}_{\bm{k}}$ will be smooth but the phase $\alpha_{\bm{k}}$ will necessarily wind by $4\pi$ around the Brillouin zone, leading to at least two vortices in the Brillouin zone where the projector is undefined due to the vanishing of the gap. These will yield a finite contribution to the second term which implies that the semimetal never satisfies the Fock bound in the ideal limit 555Notice that alternatively, if we choose a periodic gauge, $\alpha_{\bm{k}}$ can be made continuous but the Berry connection will be singular leading also to finite energy contribution. The location of the vortices which minimize this term is determined by the competition between the weak (logarithmic) repulsion between the vortices in 2D and the tendency of the vortices to sit wherever the magnetic field corresponding to ${\bm{A}}$, i.e. the Berry curvature, is maximal. As shown in Fig. 5, in the chiral limit $w_{0}=0$, the Berry curvature is relatively uniform so that vortices prefer to sit far apart at the K and K*′* point, thereby preserving $C_{3}$ symmetry. However, for finite $w_{0}/w_{1}$, the Berry curvature (in the basis which diagonalizes the sublattice operator $\sigma_{z}$ Bultinck et al. (2019b)) becomes strongly peaked at the $\Gamma$ point and for the realistic parameter $w_{0}/w_{1}\approx 0.75$, the vortices prefer to sit close to the $\Gamma$ point. This explains the energetic advantage of $C_{3}$ symmetry breaking in the realistic parameter regime considered in the numerics.

Finally, let us consider the effect of the dispersion (which is mostly generated by the interaction). Using $C_{2}\mathcal{T}$, valley and particle-hole symmetries, the form of the dispersion can be reduced to Bultinck et al. (2019b)

[TABLE]

where $f({\bm{k}})$ is a positive real function which vanishes at the two Dirac points and $\phi({\bm{k}})$ winds by $2\pi$ around each of the Dirac points. These points lie at the Moiré $K_{M}$ and $K^{\prime}_{M}$ for $\beta=0$ but move away for finite $\beta$. The energy contribution of $h_{0}({\bm{k}})$ for a state described by the matrix $Q$ is

[TABLE]

We notice that this term favors the semimetal order parameter with $Q$ given by (14) but vanishes for all the order parameters which satisfy the Fock bound. This might suggest none of these states benefit energetically from dispersion but this turns out not to be true as we will show below.

We start by writing the total energy functional in terms of $Q$ (ignoring the Hartree term which does not play a role in the competition between phases)

[TABLE]

with $E_{h}[Q]$ given by (18) and $E_{F}[Q]$ given by (11). The Hartree-Fock self consistency equation is taken by setting the variation of $E[Q]$ relative to $Q$ to zero (subject to the constraint $Q^{2}=1$) which leads to the equation

[TABLE]

From this equation, we see that the flavor polarized states which are characterized by $[Q,h_{0}]=0$ remain solution of the self-consistency equation for non-zero $h_{0}$ whereas the $C_{2}\mathcal{T}$ breaking insulators which are characterized by $\{Q,h_{0}\}=0$ are no longer solution to the self-consistency equation. To write a $C_{2}\mathcal{T}$ breaking solution, we write

[TABLE]

which satisfies the condition $Q^{2}=1$ due to the anticommutation of $Q_{C_{2}\mathcal{T}\rm I}$ and $Q_{\rm SM}$. The ${\bm{k}}$ dependence of $\theta$ can be determined from the self-consistency condition (20). In the following, we will assume for simplicity that $\theta$ is ${\bm{k}}$ independent and choose it to minimize the the total energy (19) leading to

[TABLE]

where we used the fact that $Q_{C_{2}\mathcal{T}\rm I}$ saturates the Fock bound. We notice that the denominator is the same as the second term on the right hand side of (16) which is guaranteed to be finite and roughly of the order of the interaction scale. This implies that $\theta$ is small whenever the interaction dominates the dispersion. In this case, the order parameter described by (21) is very close to the pure insulator order parameter. The corresponding energy is given by

[TABLE]

Thus, the $C_{2}\mathcal{T}$-breaking insulators actually benefit from the dispersion in quadratic order by developing a small component proportional to the semimetal order parameter. This explains why their energy is decreased in response to explicit $C_{3}$-breaking which energetically favors the semimetal (thus increasing $|E_{h}[Q_{\rm SM}]|$). It also explains their high $C_{3}$-breaking susceptibility which mainly arises from their $Q_{\rm SM}$ component.

VI Consequences for Experiment

The possibility of $C_{3}$ breaking at CN is consistent with several recent reports Choi et al. (2019); Kerelsky et al. (2019); Jiang et al. (2019); Xie et al. (2019) which observed direct evidence of $C_{3}$ breaking in STM measurements. To check the compatibilty of these measurements with our mean-field solutions, we compute the DOS for the four possible $C_{3}$-breaking states (arising for positive or negative values of $\beta$) in Fig. 4b. We see that in all cases the global DOS consists of two broad peaks (which are sometimes further split into two) separated by about 60-80 meV, in agreement with the STM measurements Kerelsky et al. (2019); Choi et al. (2019). The DOS alone however, is insufficient to distinguish the insulating and semi-metallic state, since both have very low DOS close to zero energy. One way to distinguish the two is to compute the local filling fraction defined as Jiang et al. (2019)

[TABLE]

where $\rho_{\rm LB/UB}({\bm{r}})$ denote the integrated local DOS from the upper layer for the lower/upper band. $\nu({\bm{r}})$, shown in Fig. 4, exhibits clear $C_{3}$-symmetry breaking pattern for all the four phases. The patterns for the $C_{3}$S and $C_{3}C_{2}\mathcal{T}$I QVH phases can be distinguished by $C_{2}$ symmetry which is visibly present in the former but absent in the latter (the pattern of the other $C_{3}C_{2}\mathcal{T}$I is obtained by symmetrizing the QVH result relative to $C_{2}$ with the result being very similar to the one for the $C_{3}$S breaking states up to an overall factor). The $C_{3}$S pattern is qualitatively similar to that measured in Ref. Jiang et al. (2019), with the density vanishing at the Kagomé lattice site lying on the mirror plane and changing sign twice around it, while having non-vanishing magnitude with opposite signs on the other two Kagomé lattice sites. Such structure is generic for any $C_{3}$-breaking phase which preserves mirror and particle-hole symmetries 666The BM Hamiltonian has an approximate particle-hole symmetry at small angles and a combination of $C_{2}$ and some local symmetry. Hejazi et al. (2019)

VII Conclusion

In conclusion, we have performed a momentum-space self-consistent Hartree-Fock analysis to uncover the nature of the symmetry-broken phase in twisted bilayer graphene at charge neutrality. In addition to insulating states corresponding to spin, valley or sublattice polarization, we found two $C_{3}$-breaking $C_{2}\mathcal{T}$-symmetric semimetallic solutions. Our main finding is that the existence of very small explicit $C_{3}$-breaking energetically favors one of these $C_{2}\mathcal{T}$-symmetric metallic state together with $C_{2}\mathcal{T}$-breaking insulating states. Both sets of states have similar energy within HF, strongly break $C_{3}$ symmetry and are consistent with the density of states measured in STM experiments. They can be experimentally distinguished in transport measurements or by comparing space-resolved local filling fractions in STM. We propose these two states as candidates for the insulating and conducting states observed in different experiments at the CNP and suggest that the competition between the two is settled by small details that are likely sample-dependent.

Acknowledgements.

We thank Eva Andrei, Allan MacDonald, Adrian Po, A. Thomson and M. Xie for helpful discussions. S. L., E. K., J. Y. L. and A. V. were supported by a Simons Investigator Fellowship and by NSF-DMR 1411343. E. K. was supported by the German National Academy of Sciences Leopoldina through grant LPDS 2018-02 Leopoldina fellowship.

Appendix A Single-particle physics: Bistritzer-Macdonald model

Our starting point is the Bistritzer-MacDonald (BM) model of the TBG band structure Bistritzer and MacDonald (2011), which we now briefly review. We begin with two layers of perfectly aligned (AA stacking) graphene sheets extended along the $xy$ plane, and we choose the frame orientation such that the $y$-axis is parallel to some of the honeycomb lattice bonds. Now we choose an arbitrary atomic site and twist the top and bottom layers around that site by the counterclockwise angles $\theta/2$ and $-\theta/2$ (say $\theta>0$), respectively. When $\theta$ is very small, the lattice form a Moiré pattern with very large translation vectors; correspondingly, the Moiré Brillouin zone (MBZ) is very small compared to the monolayer graphene Brillouin zone (BZ), as illustrated in Fig. 6. In this case, coupling between the two valleys can be neglected. If we focus on one of the two valleys, say $K$, then the effective Hamiltonian is given by:

[TABLE]

Here, $l=t/b\simeq\pm 1$ is the layer index, and $\tilde{f}_{l}({\bm{k}})$ is the $K$-valley electron originated from layer $l$. The sublattice index $\sigma$ is suppressed, thus each $\tilde{f}_{l}({\bm{k}})$ operator is in fact a two-component column vector. $h_{{\bm{k}}}(\theta)$ is the monolayer graphene $K$-valley Hamiltonian with twist angle $\theta$:

[TABLE]

where $v_{F}=9.1\times 10^{5}$ m/s is the Fermi velocity. Let $K_{l}$ be the $K$-vector of layer $l$, then ${\bm{q}}_{1}$ is defined as $K_{b}-K_{t}$. ${\bm{q}}_{2}=O_{3}{\bm{q}}_{1}$ is the counterclockwise $120^{\circ}$ rotation of ${\bm{q}}_{1}$, and ${\bm{q}}_{3}=O_{3}{\bm{q}}_{2}$. Finally, the three matrices $T_{i}$ are given by

[TABLE]

where we introduced an explicit $C_{3}$-breaking parameter $\beta$. We take $w_{1}=110$ meV and $w_{0}=82.5$ meV. The difference between $w_{0}$ and $w_{1}$ reflects the effect of lattice relaxation. Note that the argument of $\tilde{f}_{l}({\bm{k}})$ is measured from $K_{l}$, i.e. the monolayer UV momentum associated to $\tilde{f}_{l}({\bm{k}})$ is in fact $K_{l}+{\bm{k}}$. It is convenient to also choose a common momentum reference point for the two layers. For example, we can define $f_{t}({\bm{k}})=\tilde{f}_{t}({\bm{k}}-{\bm{q}}_{3})$ and $f_{b}({\bm{k}})=\tilde{f}_{b}({\bm{k}}+{\bm{q}}_{2})$, such that the arguments for both $f_{t}$ and $f_{b}$ are measured from $K_{t}-{\bm{q}}_{3}$, indicated as $\gamma$ in the right panel of Fig. 6.

One intuitive way of thinking about this effective Hamiltonian is to imagine a honeycomb lattice of Dirac points in the momentum space, as shown in Fig. 7b, where the two sublattices correspond to the two layers. A Dirac point at momentum ${\bm{q}}$ contributes a diagonal block $h_{{\bm{k}}-{\bm{q}}}(\pm\theta/2)$ to the Hamiltonian for the MBZ momentum ${\bm{k}}$, where the sign is determined by the sublattice that ${\bm{q}}$ belongs to. The off-diagonal blocks $T_{i}$ are nothing but the nearest-neighbor couplings of these Dirac points.

When the twist angle is near the magic angle $\theta=1.05^{\circ}$, two isolated flat bands per spin and valley appear near the charge neutrality (CN) Fermi energy, shown in Fig. 7a. These two bands are the focus of the current work.

The single particle Hamiltonian within each valley $\mathcal{H}_{\pm}$ is invariant under the following symmetries

[TABLE]

In addition, the two valleys are related by time-reversal symmetry given by

[TABLE]

Here, $\bm{\sigma},\bm{\tau}$ and $\bm{\mu}$ denote the Pauli matrices in sublattice, valley and layer spaces, respectively.

Appendix B Projecting the interaction onto the flat bands

In the following, we derive the form of the interaction when projecting onto the two flat bands. Since these two bands have a Wannier obstruction, we can only write such projected interaction in $k$-space. Let $c^{\dagger}_{\alpha}({\bm{k}})$ be the creation operator for the energy eigenstate in the band structure with internal flavor $\mu$ and band index $n$, where $\mu=(\tau,s)$ is a collective index including both valley $\tau=\pm$ and spin $s=\uparrow/\downarrow$, and $n=1,2$ represents the lower and upper bands, respectively. Also let $f^{\dagger}_{\mu,I}({\bm{q}})$ be the “elementary” continuous fermion with monolayer momentum ${\bm{q}}$, flavor $\mu=(\tau,s)$ and $I=(l,\sigma)$ representing layer and sublattice, then $c^{\dagger}$ and $f^{\dagger}$ are related to each other by the $k$-space wave functions as follows:

[TABLE]

where ${\bm{G}}$ is a Moiré reciprocal lattice vector. In the above expression, we are already using the fact that the wave functions are spin-independent. Once we choose a gauge of $u_{\tau,n;{\bm{G}},I}({\bm{k}})$ for all ${\bm{k}}$ in some MBZ, $c^{\dagger}({\bm{k}})$ are defined in terms of the $f^{\dagger}({\bm{q}})$ for those ${\bm{k}}$, and whenever necessary, we define $c^{\dagger}({\bm{k}}+{\bm{G}})=c^{\dagger}({\bm{k}})$ for any reciprical lattice vector ${\bm{G}}$, which is equivalent to defining $u_{\tau,n;{\bm{G}},I}({\bm{k}}+{\bm{G}}_{0})=u_{\tau,n;{\bm{G}}+{\bm{G}}_{0},I}({\bm{k}})$. Note that the momentum argument for $f^{\dagger}$ is unconstrained since we are using the continuum theory for monolayers of graphene. We choose the normalization $\{f_{\mu,I}({\bm{q}}),f^{\dagger}_{\mu^{\prime},I^{\prime}}({\bm{q}}^{\prime})\}=\delta_{\mu\mu^{\prime}}\delta_{II^{\prime}}\delta_{{\bm{q}}{\bm{q}}^{\prime}}$ (suppose the system size is finite), and $\langle u_{\tau,n}({\bm{k}})|u_{\tau^{\prime},n^{\prime}}({\bm{k}})\rangle:=\sum_{{\bm{G}},I}u^{*}_{\tau,n;{\bm{G}},I}({\bm{k}})u_{\tau^{\prime},n^{\prime};{\bm{G}},I}({\bm{k}})=\delta_{\tau\tau^{\prime}}\delta_{nn^{\prime}}$, which imply $\{c_{\mu,n}({\bm{k}}),c^{\dagger}_{\mu^{\prime},n^{\prime}}({\bm{k}}^{\prime})\}=\delta_{\mu\mu^{\prime}}\delta_{nn^{\prime}}\delta_{{\bm{k}}{\bm{k}}^{\prime}}$ when ${\bm{k}},{\bm{k}}^{\prime}$ are confined in the MBZ. For the purpose of projecting the interaction into these two bands, it is convenient to introduce the form factor notation:

[TABLE]

where ${\bm{q}}$ is not restricted to the first Brillouin zone. The form factors satisfy

[TABLE]

just from the definition, and also has the property

[TABLE]

due to the time-reversal symmetry.

The interaction Hamiltonian is given by

[TABLE]

where $A$ is the total area of the system and $V({\bm{q}})$ is the momentum space interaction potential, related to the real-space one by $V({\bm{q}}):=\int d^{2}{\bm{r}}V({\bm{r}})e^{-i{\bm{q}}\cdot{\bm{r}}}$. Depending on the number of gates, $V({\bm{q}})$ takes the following form in the SI units:

[TABLE]

where the screening length $d_{s}$ is nothing but the distance from the graphene plane to the gate(s). Projecting onto the two narrow bands, this Hamiltonian has the form

[TABLE]

Appendix C Hartree-Fock analysis in the chiral limit

The chiral limit of the BM model is obtained by switching off the $w_{0}$ term in (27) Tarnopolsky et al. (2019). In this limit, the bands become exactly flat at the magic angle and the eigenstates of the Hamiltonian have a simple form similar to the Landau levels on a torus In the chiral limit, the single particle Hamiltonian anticommutes with the chiral (sublattice) symmetry operator given by $\Gamma=\sigma_{z}$. In addition, it is invariant under the particle-hole symmetry $\mathcal{P}\tilde{f}_{\bm{k}}\mathcal{P}^{-1}=i\sigma_{x}\mu_{y}\tilde{f}^{\dagger}_{-{\bm{k}}}$ (modulo a small basis rotation gauging away the $\theta$ dependence). This means that we can choose the wavefunctions for the flat bands to be eigenfunctions of the sublattice operator $\sigma_{z}$ i.e. completely sublattice polarized Bultinck et al. (2019b). The wavefunctions can then be labelled by their sublattice index $\sigma=$A/B. This means that the sublattice off-diagonal components of the form factor vanishes, i.e. $\lambda_{AB;\tau,{\bm{G}}}({\bm{k}},{\bm{k}}^{\prime})=\lambda_{BA;\tau,{\bm{G}}}({\bm{k}},{\bm{k}}^{\prime})=0$. Furthermore, the action of $C_{2}\mathcal{T}$ is given by

[TABLE]

which implies that $\lambda_{AA;\tau,{\bm{q}}}({\bm{k}})=\lambda_{BB;\tau,{\bm{q}}}({\bm{k}})^{*}$. Finally, we can use $\mathcal{P}\mathcal{T}$ symmetry to restrict the form factors further by noting that it exchanges positive and negative energy eigenstates in opposite valleys and sublattices

[TABLE]

which implies that $\lambda_{A/B,\tau,{\bm{q}}}({\bm{k}})=\lambda_{B/A,-\tau,{\bm{q}}}({\bm{k}})$. Summarizing these conditions, we can write

[TABLE]

where $F_{\bm{q}}({\bm{k}})=|\lambda_{AA,+,{\bm{q}}}({\bm{k}})|$ and $\phi_{\bm{q}}({\bm{k}})=\arg\lambda_{AA,+,{\bm{q}}}({\bm{k}})$

We now investigate the Hartree-Fock solutions by looking for the minima of the Hartree-Fock energy. The Hartree-Fock energy is defined in terms of the order parameter

[TABLE]

where $\alpha,\beta$ range over spin, valley and band indices. For an insulator or a semimetal, the number of filled states is ${\bm{k}}$ independent and equal 4. This means that the order parameter $P_{\bm{k}}$ is a projector satisfying

[TABLE]

The Hartree-Fock energy can then be written as (using properties of the chiral limit)

[TABLE]

where we defined the form factor matrix $\Lambda_{\bm{q}}({\bm{k}})$ in terms of the combined index $\alpha=(s,\tau,n)$ for spin, valley, and band as

[TABLE]

The form factors decay in the separation ${\bm{q}}$ with a characteristic scale which is typically smaller than the size of the Brillouin zone. This means that the sums in the Hartree and Fock terms are dominated by the ${\bm{G}}=0$ term. For the Hartree term, this equals $4V(0)N$ and is independent of the order parameter $P({\bm{k}})$. Thus, the Hartree term has little effect on the competition between different phases and and will be neglected in the following.

We now consider possible types of symmetry-broken orders at charge neutrality. The order parameter $P({\bm{k}})$ can generally be written as

[TABLE]

Since the Fock term is the largest contribution to the mean-field energy, let us now neglect the Hartree term as well as the single-particle term $h_{0}$. We note that $\langle A,B\rangle=\mathop{\mathrm{tr}}AB$ defines a positive definite inner product on the space of hermitian matrices. Using Cauchy-Schwarz inequality, we get

[TABLE]

This inequality is satisfied if and only if $P({\bm{k}}+{\bm{q}})$ is parallel to $\Lambda_{{\bm{q}}}({\bm{k}})P({\bm{k}})\Lambda^{\dagger}_{{\bm{q}}}({\bm{k}})$ for every ${\bm{k}}$, ${\bm{q}}$. This means

[TABLE]

C.1 Intervalley coherent states

Let us now briefly discuss states which break U(1) valley symmetry. These states were argued to be energetically unfavorable in the context of twisted bilayer graphene with an aligned hBN substrate Bultinck et al. (2019a) and other Moiré materials which lack $C_{2}$ symmetry Zhang et al. (2019a); Lee et al. (2019). Here, we will show that a similar argument fails 777We thank Mike Zaletel for discussions on this point in $C_{2}$-symmetric twisted bilayer graphene due to the extra particle-hole symmetry $\mathcal{P}$ which is exact in the chiral limit. A more detailed discussion of such phases is provided in Ref. Bultinck et al. (2019b). To clarify the importance of $\mathcal{P}$, we will assume that the symmetry is not present in which case, the form factors in the two opposite valleys at the same momentum ${\bm{k}}$ are not related, i.e. $\lambda_{A/B,+,{\bm{q}}}({\bm{k}})\neq\lambda_{B/A,-,{\bm{q}}}({\bm{k}})$. Furthermore, we will assume unbroken time-reversal symmetry. The case of time-reversal symmetry breaking can be addressed similarly. The projector for an IVC state can be split into a diagonal and off-diagonal component in valley space

[TABLE]

In terms of valley resolved blocks of $P({\bm{k}})$, i.e.

[TABLE]

the second condition can be written as

[TABLE]

Since the form factors are diagonal in valley space, the Fock energy can be written as a sum of a term with only diagonal part and one with only off-diagonal parts as

[TABLE]

Here, we used the Cauchy-Shwarz inequality to go from the first to the second line. We then used the geometric-arithmetic mean inequality $\sqrt{xy}\leq\frac{x+y}{2}$ to go from the second to the third. We now see that whenever the second term in the last line is non-zero, we would conclude that the Fock energy for the IVC is larger than the energy bound for the valley polarized or unpolarized phases by an amount which is proportional to the valley off-diagonal part of the order parameter. This energy difference in (52) is generally not expected to be small unless there is a symmetry relating the same mometa in the two valleys (or equivalently a symmetry relating opposite momenta in the same valley). This is precisely the reason why the existence of particle-hole $\mathcal{P}$ symmetry forces the second term to vanish which makes the bound (52), while correct, inconclusive to rule out IVC states.