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

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.

Early experiments found clear signature of a correlated insulating state at half-filling [34] [1,23]. Subsequently, insulating ferromagnetic states were also observed at quarter and three-quarter fillings [23,24]. On the other hand, in all these experiments [1,2,23] 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 [25].
On the theory side, it was realized early on [4,35] 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 T symmetry breaking insulator, with Chern number ±1 for each spin and valley flavour was proposed in [9], along with a C 2 T symmetry preserving insulator that is believed to require mixing with remote bands [9]. An intervalley coherent order [4] was proposed as a candidate for the insulating state at half-filling while nematic orders were discussed in [6,11,36] and ferromagnetic ordering was proposed in [8]. 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 [17,37].
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 valleypolarized insulators, which break no other symmetries. In addition we observe a strong tendency to breaking spatial rotation symmetries. We find a C 2 T -breaking insulator and two distinct C 2 T -symmetric semimetallic phases which break C 3 -symmetry . The C 2 T -breaking phase has a Chern number of C = ±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 2breaking hBN substrate, potentially favoring it in situations where the sample and substrate are aligned. The gapless C 3breaking phases are obtained by bringing the two Dirac cones from the Moiré K and K very close to the Γ point. Instead of merging and opening a gap as one might normally expect, the Dirac points remain gapless since they carry the same chirality [4], 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 [21], 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 Tbreaking insulators also exhibits a strong susceptibility to C 3 symmetry breaking. The energies of the insulating C 2 Tbreaking and the semimetallic C 2 T -preserving states are low-ered 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 [32,33], which employs a continuum approximation close to K and K for a pair of graphene sheets rotated relative to each other by an angle θ. The Hamiltonian for the K valley is given by: Here, l = t/b ±1 is the layer index, and f l (k) is the Kvalley electron originated from layer l. h k (θ) is the monolayer graphene K-valley Hamiltonian with twist angle θ (see Appendix for details) and q 1 is defined as K b − K t with K l denoting the K-vector of layer l. q 2 = O 3 q 1 is the counterclockwise 2π/3 rotation of q 1 , and q 3 = O 3 q 2 . Finally, the interlayer coupling matrices are given by 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 [38,39]. 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 → (1 + β)T 1 [40].
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 [4,35,41].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 [9]. However, given the strong angle dependence and appearance of correlated states only near the magic angle [1,2,23], 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. [9].
The Hartree-Fock (HF) mean field theory is defined in terms of the projector where c α (k) is the annihilation operator for an electron at momentum k and α = (n, τ, s) is a combined index for band, valley and spin, respectively. For a gapped or semimetallic phase at CN, P satisfies tr P (k) = 4 for all k points [42]. The HF mean field Hamiltonian has the form Here, c(k) is a column vector in the index α, h 0 (k) denotes the single particle Hamiltonian and h HF (k) is given by where A is the total area. The k , G, and q summations range over the first Brillouin zone, reciprocal lattice vectors, and all momenta, respectively. V q is the interaction potential which we take to be a single-gate-screened Coulomb interaction V q = e 2 2 0 q (1 − e −2qds ) with dielectric constant = 7 and screening length equal to the gate distance d s ≈ 40 nm.
The first term in (5) is the Hartree term while the second is the Fock term. The matrix Λ q (k) contains the form factors for the single-particle where α and β denote a combined index for spin, valley and band. The HF self-consistent analysis starts by proposing an ansatz for the projector P (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 Λ q (k) is not diagonal in the band index for q = 0 since the Bloch wavefunctions u α,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 [43][44][45]. 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 (k) is given by where h 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. [9] 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 [43,44,46]. 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 [4]: spinless time-reversal 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 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 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 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 T symmetry, is the equal [4,15], 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 [41,47,48]. This invariant is protected by C 2 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 [49]. Thus the non-trivial w 2 invariant implies that the two Dirac points cannot be removed without breaking C 2 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 [17,22]. The IVC state is known to have higher energy than SP/VP states for isolated bands with non-vanishing valley Chern number [17,37]. These arguments do not generalize to TBG where the extra symmetries of the problem complicate the discussion (see Appendix and Ref. [50]). 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 spinvalley 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 [50].

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 β. There are two types of gapped solutions corresponding to either flavorpolarized (spin/valley) states or C 2 T -breaking (C 2 T I) insulators. There are three gapped solutions corresponding to spin-polarized (SP), valley-polarized (VP) and C 2 T -breaking (C 2 T I) insulators. The latter does not break C 3 when β = 0 but develops a large C 3 breaking component for β = 0. The extent of C 3 symmetry breaking can be quantified by defining which vanishes for any C 3 -symmetric state. Here, |ψ k are the occupied single-particle eigenstates of the HF Hamiltonian with momentum k (for a given flavor). The value of χ C3 for the C 2 T I state is shown in Fig. 2 as a function of χ C3 for the corresponding non-interacting states arising from explicit C 3breaking parameter β = 0 in the BM model. We can clearly see from the figure that a relatively small χ C3 (BM ) ∼ 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 T I state. This serves to show that the C 2 T I has very large susceptibility to C 3 symmetry breaking. In the following, we will refer to this state for positive and negative β as C 3 C 2 T Iy and C 3 C 2 T Ix, respectively. In addition to this insulating state, there are two distinct C 2 -preserving semimetallic phases which spontaneously break C 3 even for β = 0 which we denote by C 3 Sx and C 3 Sy, 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. [43] which used a site-local ansatz for the interactions. Examples of the band structures for the C 3 Sy and C 3 C 2 T Iy states are shown in Fig. 3. The C 2 T I phase obtained here is characterized by a Chern number of ±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 T is broken but C 2 and spin rotation are preserved, we obtain a quantum anomalous Hall (QAH) insulator with total Chern number C = ±4, (ii) if spin rotation is broken we obtain a quantum-spin Hall (QSH) insulator with opposite Chern number C s = ±2 for opposite spins, (iii) if C 2 is broken but T and spin rotation symmetry are preserved, then a quantum valley Hall (QVH) insulator obtains with opposite Chern number C v = ±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 ±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 samples [51]. The energies of the four possible C 2 T states are very similar differing only by very small ∼ 10 −3 meV Hartree corrections. Similarly, the flavor-polarized state can be spin-polarized, breaking SU(2) spin rotation, valley-polarized breaking T or a spin-valleylocked states which breaks both but preserves a combination of T and π-spin rotation. These states are degenerate on the mean field level and are only distinguished by intervalley Hund's coupling [17,22] neglected in this study.
At β = 0, the energies of the insulating VP/SP and C 2 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 [52]. 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 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. [50].
Once β becomes non-zero, the energy of the C 2 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 β) goes down in energy becoming more energetically favorable to the VP/SP state around β = ±3 × 10 −4 . For larger values of β 4 × 10 −4 , the energies of the insulating C 2 T -breaking phase and the semimetallic C 2 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 T -breaking insulator or a C 2 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 T breaking insulating regions [53].

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. [52] 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 σ z . As a result, the flatband wavefunctions can be chosen to be eigenfunctions of σ z such that the wavefunction σ z u A/B,τ,k = ±u A/B,τ,k . In this limit, we can show (see Appendix for details) that the form factor matrix has the simple form where F q (k) and Φ q (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 [50]. 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 (k) as In terms of Q, the Fock energy is given by up to a constant. We note that A, B = tr AB defines a positive definite inner product on the space of hermitian matrices. Using Cauchy-Schwarz inequality, we get This inequality is satisfied if and only if Q(k + q) is parallel to Λ q (k)Q(k)Λ † q (k) for every k and q which implies A k-independent flavor (spin,valley) polarized state obviously satistfies this constraint since σ z τ z is diagonal in flavor, thereby saturating the Fock bound. This state can either be spin-polarized Q = s z , valley polarized Q = τ z or spinvalley locked Q = τ z s z .
Next let us discuss C 2 T breaking insulating solutions. A C 2 T breaking solution to Eq. 13 is obtained by taking Q k = σ z (since C 2 T is off diagonal in the sublattice index). The resulting state saturates the Fock bound and is characterized by a Chern number ±1 in the lower band. Its energy competition with the flavor-polarized states is settled by the effects of the single-particle term h 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 β. Depending on its structure in flavor space, the C 2 T -breaking insulator may correspond to one of four states: (i) Q = s 0 τ 0 σ z breaks C 2 but not T and corresponds to a valley Hall state, (ii) Q = s 0 τ z σ z breaks T but not C 2 and corresponds to a quantum Hall state with Chern number ±4, (iii) Q = s z τ 0 σ z corresponds to a valleyspin-Hall state where the Chern number is invariant under flipping spin and valley, and (iv) Q = s z τ z σ z corresponds to a spin-Hall state.
Finally, let us consider semimetallic states which break neither flavor nor C 2 T symmetry. These can be generally described (within each flavor) by the order parameter Substituting in the Fock energy (11) yields To simplify further, we make the reasonable assumption that F q (k) decays quickly with the relative momentum q which enables us to expand Φ q (k) and α k+q in q. Noting that Φ q (k) = q · A k + O(q 2 ) with A k the Berry connection −i u A,K,k |∇ k |u A,K,k , we get 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 [50,52], 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 A k will be smooth but the phase α k will necessarily wind by 4π 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 [54]. 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 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 σ z [50]) becomes strongly peaked at the Γ point and for the realistic parameter w 0 /w 1 ≈ 0.75, the vortices prefer to sit close to the Γ 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 T , valley and particle-hole symmetries, the form of the dispersion can be reduced to [50] h 0 (k) = σ x f (k)e iφ(k)σzτz (17) where f (k) is a positive real function which vanishes at the two Dirac points and φ(k) winds by 2π around each of the Dirac points. These points lie at the Moiré K M and K M for β = 0 but move away for finite β. The energy contribution of h 0 (k) for a state described by the matrix Q is 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) 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 which satisfies the condition Q 2 = 1 due to the anticommutation of Q C2T I and Q SM . The k dependence of θ can be determined from the self-consistency condition (20). In the following, we will assume for simplicity that θ is k independent and choose it to minimize the the total energy (19) leading to where we used the fact that Q C2T 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 θ 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 Thus, the C 2 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 SM ]|). It also explains their high C 3breaking susceptibility which mainly arises from their Q SM component.

VI. CONSEQUENCES FOR EXPERIMENT
The possibility of C 3 breaking at CN is consistent with several recent reports [43][44][45][46] 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 β) 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 [43,45]. 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 [44] where ρ LB/UB (r) denote the integrated local DOS from the upper layer for the lower/upper band. ν(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 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 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. [44], 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 [55] and a combination of C 2 and some local symmetry. [56] VII. CONCLUSION In conclusion, we have performed a momentum-space selfconsistent 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 3breaking C 2 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 T -symmetric metallic state together with C 2 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 spaceresolved 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.
(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: Here, l = t/b ±1 is the layer index, andf l (k) is the K-valley electron originated from layer l. The sublattice index σ is suppressed, thus eachf l (k) operator is in fact a two-component column vector. h k (θ) is the monolayer graphene K-valley Hamiltonian with twist angle θ: where v F = 9.1 × 10 5 m/s is the Fermi velocity. Let K l be the K-vector of layer l, then q 1 is defined as is the counterclockwise 120 • rotation of q 1 , and q 3 = O 3 q 2 . Finally, the three matrices T i are given by where we introduced an explicit C 3 -breaking parameter β. 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 off l (k) is measured from K l , i.e. the monolayer UV momentum associated tof l (k) is in fact K l + k. It is convenient to also choose a common momentum reference point for the two layers. For example, we can define f t (k) =f t (k − q 3 ) and f b (k) =f b (k + q 2 ), such that the arguments for both f t and f b are measured from K t − q 3 , indicated as γ 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 q contributes a diagonal block h k−q (±θ/2) to the Hamiltonian for the MBZ momentum k, where the sign is determined by the sublattice that 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 θ = 1.05 • , 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 H ± is invariant under the following symmetries In addition, the two valleys are related by time-reversal symmetry given by Here, σ, τ and µ denote the Pauli matrices in sublattice, valley and layer spaces, respectively. 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 † α (k) be the creation operator for the energy eigenstate in the band structure with internal flavor µ and band index n, where µ = (τ, s) is a collective index including both valley τ = ± and spin s =↑ / ↓, and n = 1, 2 represents the lower and upper bands, respectively. Also let f † µ,I (q) be the "elementary" continuous fermion with monolayer momentum q, flavor µ = (τ, s) and I = (l, σ) representing layer and sublattice, then c † and f † are related to each other by the k-space wave functions as follows: where 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 τ,n;G,I (k) for all k in some MBZ, c † (k) are defined in terms of the f † (q) for those k, and whenever necessary, we define c † (k + G) = c † (k) for any reciprical lattice vector G, which is equivalent to defining u τ,n;G,I (k + G 0 ) = u τ,n;G+G0,I (k). Note that the momentum argument for f † is unconstrained since we are using the continuum theory for monolayers of graphene. We choose the normalization {f µ,I (q), f † µ ,I (q )} = δ µµ δ II δ qq (suppose the system size is finite), and u τ,n (k)|u τ ,n (k) := G,I u * τ,n;G,I (k)u τ ,n ;G,I (k) = δ τ τ δ nn , which imply {c µ,n (k), c † µ ,n (k )} = δ µµ δ nn δ kk when k, k 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: where q is not restricted to the first Brillouin zone. The form factors satisfy just from the definition, and also has the property λ mn,τ ;q (k) = λ mn,−τ ;−q (−k) * (B4) due to the time-reversal symmetry.
The interaction Hamiltonian is given by where A is the total area of the system and V (q) is the momentum space interaction potential, related to the real-space one by V (q) := d 2 rV (r)e −iq·r . Depending on the number of gates, V (q) takes the following form in the SI units: 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 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 (A3) [52]. 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 Γ = σ z . In addition, it is invariant under the particle-hole symmetry Pf k P −1 = iσ x µ yf † −k (modulo a small basis rotation gauging away the θ dependence). This means that we can choose the wavefunctions for the flat bands to be eigenfunctions of the sublattice operator σ z i.e. completely sublattice polarized [50]. The wavefunctions can then be labelled by their sublattice index σ =A/B. This means that the sublattice off-diagonal components of the form factor vanishes, i.e. λ AB;τ,G (k, k ) = λ BA;τ,G (k, k ) = 0. Furthermore, the action of C 2 T is given by which implies that λ AA;τ,q (k) = λ BB;τ,q (k) * . Finally, we can use PT symmetry to restrict the form factors further by noting that it exchanges positive and negative energy eigenstates in opposite valleys and sublattices PT u A,τ,k = e iη A,τ,k u B,−τ,k , PT u B,τ,k = e −iη B,τ,k u A,−τ,k which implies that λ A/B,τ,q (k) = λ B/A,−τ,q (k). Summarizing these conditions, we can write where F q (k) = |λ AA,+,q (k)| and φ q (k) = arg λ AA,+,q (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 where α, β range over spin, valley and band indices. For an insulator or a semimetal, the number of filled states is k independent and equal 4. This means that the order parameter P k is a projector satisfying P (k) 2 = P (k) = P † (k), tr P (k) = 4 (C5) The Hartree-Fock energy can then be written as (using properties of the chiral limit) where we defined the form factor matrix Λ q (k) in terms of the combined index α = (s, τ, n) for spin, valley, and band as The form factors decay in the separation 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 G = 0 term. For the Hartree term, this equals 4V (0)N and is independent of the order parameter P (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 (k) can generally be written as Since the Fock term is the largest contribution to the mean-field energy, let us now neglect the Hartree term as well as the singleparticle term h 0 . We note that A, B = tr AB defines a positive definite inner product on the space of hermitian matrices. Using Cauchy-Schwarz inequality, we get This inequality is satisfied if and only if P (k + q) is parallel to Λ q (k)P (k)Λ † q (k) for every k, q. This means P (k + q) = e iΦq(k)σzτz P (k)e −iΦq(k)σzτz ⇒ Q(k + q) = e iΦq(k)σzτz Q(k)e −iΦq(k)σzτz (C12)

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 [37] and other Moiré materials which lack C 2 symmetry [17,22]. Here, we will show that a similar argument fails [57] in C 2 -symmetric twisted bilayer graphene due to the extra particlehole symmetry P which is exact in the chiral limit. A more detailed discussion of such phases is provided in Ref. [50]. To clarify the importance of 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 k are not related, i.e. λ A/B,+,q (k) = λ B/A,−,q (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 In terms of valley resolved blocks of P (k), i.e. P = P + P 12 the second condition can be written as P 2 + + P 12 P 21 = P + , P 2 − + P 21 P 12 = P − .
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 k,q V q [|λ A,+,q (k)| 2 tr P + (k) 2 tr P + (k + q) 2 + |λ A,−,q (k)| 2 tr P − (k) 2 tr P − (k + q) 2 + |λ A,+,q (k)||λ A,−,q (k)| tr P o (k) 2 tr P o (k + q) 2 ] ≥ − 1 2A k,q V q [|λ A,+,q (k)| 2 tr P + (k) 2 + |λ A,−,q (k)| 2 tr P − (k) 2 + |λ A,+,q (k)||λ A,−,q (k)| tr P o (k) 2 ] = − 1 2A k,q V q [|λ A,+,q (k)| 2 tr(P + (k) − P 12 (k)P 21 (k)) + |λ A,+,q (k)| 2 tr(P − (k) − P 21 (k)P 12 (k)) + 2|λ A,+,q (k)||λ A,−,q (k)| tr(P 12 (k)P 21 (k))] Here, we used the Cauchy-Shwarz inequality to go from the first to the second line. We then used the geometric-arithmetic mean inequality √ xy ≤ 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 (C16) 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 P symmetry forces the second term to vanish which makes the bound (C16), while correct, inconclusive to rule out IVC states.