Non-Abelian Dirac node braiding and near-degeneracy of correlated phases at odd integer filling in magic angle twisted bilayer graphene
NNon-Abelian Dirac node braiding and near-degeneracy of correlated phases at oddinteger filling in magic angle twisted bilayer graphene
Jian Kang
1, 2, ∗ and Oskar Vafek
2, 3, † School of Physical Science and Technology & Institute for Advanced Study, Soochow University, Suzhou, 215006, China National High Magnetic Field Laboratory, Tallahassee, Florida, 32310, USA Department of Physics, Florida State University, Tallahassee, Florida 32306, USA
We use the density matrix renormalization group (DMRG) to study the correlated electron statesfavored by the Coulomb interaction projected onto the narrow bands of twisted bilayer graphenewithin a spinless one-valley model. The Hilbert space of the narrow bands is constructed from apair of hybrid Wannier states with opposite Chern numbers, maximally localized in one directionand Bloch extended in another direction. Depending on the parameters in the Bistritzer-Macdonaldmodel, the DMRG in this basis determines the ground state at one particle per unit cell to be eitherthe quantum anomalous Hall (QAH) state or a state with zero Hall conductivity which is nearly aproduct state. Based on this form, we then apply the variational method to study their competition,thus identifying three states: the QAH, a gapless C T symmetric nematic, and a gapped C T symmetric stripe. In the chiral limit, the energies of the two C T symmetric states are found to besignificantly above the energy of the QAH. However, all three states are nearly degenerate at therealistic parameters of the Bistritzer-Macdonald model. The single particle spectrum of the nematiccontains either a quadratic node or two close Dirac nodes near Γ. Motivated by the Landau leveldegeneracy found in this state, we propose it to be the state observed at the charge neutrality pointonce spin and valley degeneracies are restored. The optimal period for the C T stripe state is foundto be 2 unit cells. In addition, using the fact that the topological charge of the nodes in the C T nematic phase is no longer described simply by their winding numbers once the translation symmetryis broken, but rather by certain elements of a non-Abelian group that was recently pointed out, weidentify the mechanism of the gap opening within the C T stripe state. Although the nodes at theFermi energy are locally stable, they can be annihilated after braiding with other nodes connectingthem to adjacent (folded) bands. Therefore, if the translation symmetry is broken, the gap at oneparticle per unit cell can open even if the system preserves the C T and valley U (1) symmetries,and the gap to remote bands remains open. I. INTRODUCTION
Since the discovery of correlated insulating phasesand superconductivity (SC) in magic angle twisted bi-layer graphene (TBG) and other moire systems ,tremendous theoretical effort has been devoted to-wards understanding the properties and the mecha-nisms of these correlated electron phenomena . Sig-nificant progress has been achieved in understandingthe topological band properties of this material andother moire systems . Furthermore, several ap-proaches have revealed the similarity betweenthe quantum hall ferromagnetism and the insulatingstates observed at the even integer fillings. However, twoentirely different insulating phases have been observed atthe filling of ν = 3 . While the (quantum) anomalousHall (QAH) state has been readily identified when oneof the layers of the TBG is aligned with the hexagonalboron nitride (hBN) substrate , the observed gappedinsulating state at ν = 3 without the hBN alignment –and without anomalous Hall conductance – is much lessunderstood.The experiments, as well as the band calculations withlattice corrugation , have shown that the TBG near themagic angle contains four spin degenerate narrow bandsseparated from other remote bands by a finite band gap.The four copies of Dirac nodes at K and K (cid:48) per each of the two spin projections of the TBG are protected by C T – two fold rotation about the axis perpendicular tothe graphene plane followed by time reversal – and theconservation of the number of fermions within each valleyi.e. U v (1) . We wish to stress that this means that C T and U v (1) symmetries are necessary for stable nodesto exist, but they are not sufficient . As we discuss below,the spectrum may be gapped despite the presence of the C T and U v (1) symmetries, and despite maintaining thegap to remote bands, when the moire lattice translationsymmetry is broken.A more familiar example of Dirac nodes protected bysymmetries is the monolayer graphene, with two masslessDirac fermions per spin projection and without spin-orbitcoupling. In this case, the nodes are said to be protectedby the time reversal ( T ) and inversion ( I ) symmetries.Nevertheless, strong breaking of the rotational symme-try can in principle result in an insulating state , de-spite preserving T and I throughout the process of gapopening. This happens when the Dirac nodes with oppo-site chirality move across the Brillouin zone, meet, andannihilate.Unlike in the monolayer graphene example, however,the two Dirac nodes in magic angle twisted bilayergraphene have the same chirality and thus cannot be an-nihilated by simply meeting together. Strong breakingof rotational symmetry alone will therefore not producean insulator. Thus, in the simplest scenario with polar- a r X i v : . [ c ond - m a t . s t r- e l ] J u l ized spin and valley, it seems that the gap at the Diracnodes can be opened only by breaking the C T sym-metry. The mentioned QAH, observed at the filling of ν = 3 with hBN alignment, is an example of such C T symmetry breaking. However, as mentioned, without thehBN alignment experiments demonstrated that the sys-tem at ν = 3 is in a gapped state without anomalous Halleffect .One of the goals of this paper is to explore the mech-anism of gap opening in a C T and valley U (1) sym-metric system and how the energy of the resulting statecompetes with QAH at odd integer filling. Such a statewould be insulating and not display the anomalous Halleffect, and thus be consistent with the experiments at ν = 3 without hBN alignment; the connection to thespinless one-valley model is to simply spin and valley po-larize one hole per moire unit cell. We find that in thechiral limit , the density matrix renormalization group(DMRG) identifies the QAH as the ground state. In themore realistic case, however, DMRG always produces anon-QAH state, even when the initial state for the al-gorithm is set to be the QAH state. This result is alsoconfirmed by minimizing the energy of the trial wave-function inspired by studying the correlations in the non-QAH state obtained in DMRG. Our variational analysisdiscovers three competing states: the QAH and two C T symmetric states with dramatically different fermion ex-citation spectra. Furthermore, applying the insights ofrecent work by Wu, Soluyanov and Bzduˇsek , we canidentify the mechanism of the transition between thesetwo C T symmetric states via assignment of the non-Abelian topological charges to Dirac nodes once moirelattice translation symmetry is broken. This naturallyexplains why the gap can be opened while preserving C T and valley U (1) symmetries. Interestingly, these C T symmetric states, with polarized spin and valley de-grees of freedom, are variationally nearly degenerate withthe QAH state even though they are not connected to theQAH by U (4) symmetry (or U (4) × U (4) symmetryin the chiral limit ). As a consequence, the manifoldof the low energy states in the realistic TBG appears tobe larger than QAH-related states. We should also men-tion in passing that our earlier approach based on maxi-mally localized Wannier states in all directions –relationto which we discuss in the section below– did identify aperiod 2 stripe state as an insulating candidate for theodd integer filling .Although not gapped, the single particle excitationspectrum of the C T symmetric nematic state obtainedvariationally is also interesting in that it displays either aquadratic node or two close Dirac nodes near Γ point,i.e. the center of the moire mini-Brillouin zone, whenthe electron-electron interactions dominate the kineticenergy of the narrow band states. This is in sharp con-trast to the single particle spectrum obtained when thekinetic energy of the narrow bands dominates, in whichcase the two Dirac cones sit at the corners of the moiremini-Brillouin zone. In the latter case, the sequence of the Landau levels, restoring the spin and valley degen-eracy, would be ν = ± , ± , ± , . . . , inconsistent withthe experimentally observed sequence near the magicangle ( ∼ . ◦ ) ν = ± , ± , ± , ± . . . . In the formercase, however, the quadratic node at the moire Brillouinzone center would indeed produce the experimentally ob-served sequence because the Landau levels are doubly de-generate at zero energy, and non-degenerate at all otherenergy levels (not including the spin and valley degen-eracy). Two close Dirac nodes would also produce theexperimentally observed sequence , except for a verysmall magnetic field below which the sequence would re-vert to the ν = ± , ± , ± , . . . . In practice, no Lan-dau quantization is seen at very small magnetic field,so two close nodes are also consistent with the data atthe charge neutrality point (CNP). Interestingly, becausethis explanation relies on the electron-electron Coulombinteraction dominating the kinetic energy of the narrowbands, it would suggest that a useful probe of their rela-tive strength at different twist angles is the Landau levelsequence. Indeed, at the higher twist angle ( ∼ . ◦ )the observed sequence reverts to ν = ± , ± , ± , . . . ,suggesting that at this higher angle the kinetic energydominates.Note that our goal is not to identify strictly a single state that has the lowest energy for our Hamiltonian.Rather it is to identify a group of competing low energystates if they lie close in energy . This is because smallterms in the Hamiltonian beyond currently accepted the-oretical models, and beyond control of the experimen-talists, can tip the balance and select different groundstate from this near degenerate group. There are ex-perimental indications that this is indeed happening, inparticular because nominally same fabrication protocolsresult in different phase diagrams, for example amongthe Columbia/UCSB and the Barcelona groups . Ourstrategy is therefore to identify the leading candidatesfor the ground state based on comparing the competingstates’ robust phenomenological properties with existingexperiments.We reach the above conclusions by starting with the(energy eigen-) Bloch states for the narrow band obtainedfrom the Bistritzer-Macdonald (BM) model . This con-tinuum model has two parameters, w and w , relatedto interlayer AA and AB couplings respectively. Due tothe lattice relaxation, w is generally smaller than w ,and w /w ∼ .
83 as obtained by STM . Assuming boththe spin and valley are polarized (i.e. spinless one valleymodel), we consider how the ground state at one-particleper unit cell could depend on this ratio. For each differ-ent value of the ratio, we solve the BM model to obtainthe Bloch states, construct the hybrid WSs, and projectthe Coulomb interactions onto the hybrid WSs. By ne-glecting the impact of the remote bands, the basis ofthe hybrid WSs allows us to run DMRG with projectedinteractions only. In addition, we propose a trial wave-function for the ground state based on the outcome ofDMRG. Starting from this trial state, we minimize theenergy to study the ground states and fermion excitationswith both interactions and kinetic terms.The rest of the paper is organized as follows: in thenext section we describe the continuum model withinwhich we compute the hybrid Wannier states, discusstheir relation to the exponentially localized states in alldirections , and express the kinetic energy and theelectron-electron Coulomb interaction in the hybrid Wan-nier basis. In Section III we describe the results of ourDMRG calculation. In Section IV, we analyze the trialstate inspired by the results from DMRG and compute itssingle fermion excitation spectrum. In section V we anal-yse improved trial states which further lower the energy.We also study their excitation spectrum and its evolu-tion from gapless C T nematic to gapped C T stripeusing the topological methods discussed above. Finally,Section VI is reserved for discussion. Various technicaldetails of our calculations are presented in the Appendix. II. CONTINUUM LIMIT HAMILTONIAN ANDTHE NARROW BAND HYBRID WANNIERSTATES
The starting point of our analysis is the continuumHamiltonian H BM = (cid:126) v f (cid:88) l = t,b (cid:88) p ψ † l, p σ l,θ · p ψ l, p + (cid:88) p (cid:88) j =1 (cid:16) ψ † b, p + q j T j ψ t, p + h.c. (cid:17) , (1)where ψ l, p is the fermion operator that annihilates thestate with the momentum of p on layer l . It contains twocomponents corresponding the two sublattices on eachlayer. σ t/b,θ = e − i θ σ z ( σ x , σ y ) e i θ σ z , with θ being the twist angle. Suppose K t and K b are thetwo Dirac points on two de-coupled layers. The secondterm in Eqn. 1 is the inter-layer coupling with q = K b − K t = k θ (0 , −
1) and k θ = | K t − K b | = 2 | K t | sin θ/
2, and q , = k θ ( ± √ , ). In addition, T j = w σ + w (cid:18) cos 2 π ( j − σ + sin 2 π ( j − σ (cid:19) . The Bloch state is labeled by its crystal momentum, k = p l + K l , where p l is the momentum at the layer l .Although each Bloch state contains multiple p l ’s in theBM model, these momenta with the same layer index l differ from each other only by reciprocal lattice vectors,and thus k is uniquely defined if it is restricted in theBrillouin zone (BZ).As discussed in Refs. , the parameter w is ameasure of the tunneling within the AA regions while w within AB/BA regions. In our calculation, w is fixed to (a) (b) - - - - - - - - - - - - - - - - (c) FIG. 1. (a) The schematic plot for the moire unit cell, withthe lattice vectors of L and L . At the ‘magic’ angle, theirlength is about 13nm. The colored circles refer to the AA,AB, and BA regions. Due to the lattice relaxation effects, theAB and BA regions become larger than AA region. (b) Theschematic plot of the moire BZ. g and g are two reciprocallattice vectors, and the region enclosed by dashed line is thefirst BZ. The band dispersions along the arrowed path withdifferent w /w are illustrated in (c). Experimentally relevantvalue is w /w ∼ .
83. The energies are normalized by v f k θ ,where k θ = 2 | K | sin θ/ | K | = 4 π/ (3 a ), a = 0 . nm , and θ is the twist angle. be 0 . v f k θ and w is allowed to vary . The relativearea of AA to AB/BA is ∼ .
83 as measured in STM .The spectrum of this Hamiltonian is shown in the Fig-ure. 1(c) for a range of parameters w /w starting fromthe chiral limit where w = 0 and the narrow band isexactly flat .We assume that the Coulomb interaction ∼ meV acts mainly in the subspace of the narrow bands whereits effects may be non-perturbative . We also assumethat its mixing of the remote bands can be treated per-turbatively due to the presence of the gap between thenarrow band and the remote bands; its value is at least ∼ − . Therefore, in order to study the effects ofthe electron-electron interactions, we previously con-structed a complete and orthogonormal basis for the nar-row band which is exponentially localized in all direc-tions . To do this, we used a microscopic tight-bindingmodel at a commensurate twist angle and D symme-try . We found that the Coulomb interaction projectedonto such basis leads to a homogeneous SU(4) ferromag-netic state at ν = ±
2. We also proposed a ferromagneticperiod 2 stripe state at ν = 3 as a good candidate for theinsulating state observed at this filling without the hBNalignment.The narrow bands obtained within the continuumHamiltonian (1) carry non-trivial (fragile) topology .This is due to discarding the (expectedly small) mixingbetween the valleys, thus making the particle numberwithin each valley separately conserved. H BM is indeedinvariant under U v (1), and due to the invariance of H BM under the C T symmetry transformation, the non-trivialtopology of the narrow bands is intimately linked to thecombined symmetry protecting two Dirac cones with thesame winding number . However, unlike in the caseof a Chern insulator, the non-trivial topology here does not obstruct the construction of exponentially localizedWannier states (WS) in both directions , but it does ob-struct such states from transforming in a simple way un-der both U v (1) and C T . For example, if the WSs trans-form simply under U v (1) by acquiring an overall phase,then they cannot simply acquire a phase under C T . Be-cause the transformations which relate the Bloch statesof H BM and such WSs are perfectly unitary, no infor-mation is lost, and the C T transformed WS can stillbe expressed exactly as a linear superposition of the ex-ponentially localized WSs in both directions. The roleof the non-trivial topology is to prevent this linear su-perposition to be confined to a single site. Instead, the C T transformed WS is reconstructed from a linear su-perposition of WSs whose centers lie within the regionsurrounding the transformed WS. The size of such re-gion is determined by the exponential decay length ofthe WS, and the convergence towards full symmetry isachieved exponentially fast with increasing such region .In this respect it is perhaps helpful to reiterate that if theproblem is solved on a microscopic tight-binding latticewith ∼ carbon sites within the unit cell insteadof in the continuum approximation, in, say, the D con-figuration, then U v (1) are C T are emerging, but theyare not exact; the exponentially localized WSs in bothdirections obtained in Ref. thus transform simply un-der all exact symmetries of the starting model. This ap-proach based on exponentially localized WSs allowed usto obtain an explicit understanding of the form of thereal space interaction and, importantly, to identify thegeneralized spin-valley ferromagnetism as the dominantordering tendency in the strong coupling limit. We alsolinked this tendency to the nontrivial topological bandproperties . In order to gain a clearer understanding of the effectsof the Coulomb interaction on the U v (1) and C T sym-metries of the low energy states, in this paper we choseto work in a Wannier basis which is localized only in onedirection. In the other direction, our hybrid WSs behaveas extended Bloch waves (see Fig.2). Additional advan-tages of this basis are that the topology of the narrowbands of H BM is more transparent , and that stateswith broken translational symmetry in the localized di-rection can be readily described. Moreover, in the basis ofthe hybrid Wannier orbitals the QAH state is completelyunentangled. Because at w /w = 0 the QAH state canbe analytically shown to be the exact ground state ofprojected interactions , and because it is gapped inthis model, it is stable at small but finite w /w . We cantherefore study within DMRG whether it ‘melts away’as w /w increases beyond a critical value by initializingthe DMRG with QAH. Since, as we will see, it does, weknow with certainty that there is a quantum phase tran-sition into a state different from QAH even for finite bonddimension which limits every numerical calculation, be-cause low bond dimension would favor QAH. The disad-vantage is the complicated form the Coulomb interactiontakes in the hybrid Wannier basis making its effect lesstransparent. A. Hybrid WSs for the narrow bands
We follow the approach outlined in Ref. to constructthe hybrid WS, which are maximally localized along the L -direction and extended Bloch waves along the L -direction (see Fig. 1a).Such hybrid WSs are the eigenstates of the projectedposition operator satisfying periodic boundary condi-tions ˆ O = ˆ P e − i N g · r ˆ P , (2)where ˆ P is the projection operator onto the narrowbands. g is the primitive vector of the reciprocal lat-tice, and N is the number of unit cells along the direc-tion of L in the entire lattice with periodic boundaryconditions. We thus haveˆ O | w ± ( n, k g ) (cid:105) = e − πi N ( n + (cid:104) x ± (cid:105) k / | L | ) | w ± ( n, k g ) (cid:105) . (3)The hybrid WSs | w α ( n, k g ) (cid:105) are labeled by their mo-mentum k along g which is conserved by ˆ O and theindex n of the unit cell along L (see Appendix for de-tails of the derivation ); α = ± k . When k is close to k = 0 or 1,the hybrid WSs contain two peaks centered around AAin the localized direction of L so that (cid:104) x (cid:105) ± ≈ . | L | ,whereas the hybrid WSs with k close to 0 . L . (a) (b) FIG. 2. The hybrid WSs in real space with w /w = 0 .
85 and the Chern index of −
1. “Top/Bot” refers to the top and bottomlayer respectively. “A/B” refers to the sublattice A and B respectively. (a) With k = 0, the state contains two neighboringpeaks in AA regions along the localized direction L (See Fig. 1(a)). (b) With k = 0 .
5, the state contains only one peak in thelocalized direction L . In both cases, the state mainly occupies sublattice B, showing that the Chern index is related to thesublattice polarization. The (cid:104) x ± (cid:105) k physically represents the average of the po-sition operator within each 1D unit cell whose depen-dence on the conserved momentum k is shown in theFig. 3. Such shapes were previously obtained in Ref. .The two curves display the winding numbers of ± k increases from 0 to g i.e. the averageposition of one set of states slides to the right and theother set of states to the left under the increase of thewavenumber k , similar to Landau gauge Landau levelstates in opposite magnetic field . This makes thenontrivial topology of the system explicit: within eachvalley, the two narrow bands of H BM can be decomposedinto one Chern +1 band and one Chern − w /w are topologically the same, the shapes of (cid:104) x ± (cid:105) k clearly differ for different w /w . As seen in the Fig. 3,the slope of (cid:104) x (cid:105) near k = 0 . w /w . In the chiral limit where w /w = 0, theslope is almost the same as for Landau states, whilenear the more realistic value w /w = 0 .
8, the slopeat k = 0 .
25 nearly vanishes and the curve is very flatand thus insulating-like throughout most of the BZ. Itis only close to the BZ boundary that the winding num-bers are established. As shown below, the nature of themany-body ground state in the strong coupling limit issensitive to the shape of the (cid:104) x ± (cid:105) k curves, not just theirtopology. We carefully choose the phases of the hybrid WSs sothat the states | w ± ( n, k g ) (cid:105) are continuous functions ofthe momentum k and also satisfy the following proper-ties: | w ± ( n, ( k + 1) g (cid:105) = | w ± ( n ± , k g ) (cid:105) (4) C T | w ± ( n, k g (cid:105) = | w ∓ ( − n, k g ) (cid:105) (5) C (cid:48)(cid:48) | w ± ( n, k g (cid:105) = e − πink | w ∓ ( n, (1 − k ) g ) (cid:105) (6)ˆ T L | w ± ( n, k g ) (cid:105) = | w ± ( n + 1 , k g ) (cid:105) (7)ˆ T L | w ± ( n, k g ) (cid:105) = e − πik | w ± ( n, k g ) (cid:105) , (8)where C (cid:48)(cid:48) is the two fold rotation around the in-plane x axis as shown in Fig. 1(a) and ˆ T L , are translationoperators by L , . Note that the phase on the right side ofEqn. 6 cannot be removed, because C (cid:48)(cid:48) does not commutewith ˆ T L and, as a consequence, the extra phase becomesnecessary as long as the unit cell index n is non-zero. B. Kinetic energy and the Chern Bloch states
The kinetic energy can be written in the hybrid Wan-nier basis as H kin = (cid:88) nn (cid:48) k (cid:88) αα (cid:48) = ± t αα (cid:48) ( n − n (cid:48) , k ) d † α,n,k d α (cid:48) ,n (cid:48) ,k , (9) - - - - - - - - - - - - FIG. 3. The phase of the eigenvalue of the Wilson loop oper-ator of the two valley-polarized narrow bands with differentvalues of w /w . Note that (cid:104) x (cid:105) has an odd winding numberof ± k changes from 0 to 1, illustratingthe nontrivial topological properties of the bands. Blue andred curves are (cid:104) x (cid:105) / | L | with Chern index of +1 and − where the 1D hopping matrix elements are t αα (cid:48) ( n − n (cid:48) , k ) = (cid:104) w α ( n, k g ) | H BM | w α (cid:48) ( n (cid:48) , k g ) (cid:105) , (10)and where d † α,n,k creates the hybrid WS with the Chernindex α , the unit cell n , and the momentum k g . Thekinetic energy operator is diagonal in k , but not in n . Dueto C T symmetry whose action on our basis follows (5),and the fact that H BM is Hermitian, it is straightforwardto show that t ++ ( δn, k ) = t ∗−− ( − δn, k ) = t −− ( δn, k ) . (11)There are no additional constraints on the hopping con-stants t + − imposed by C T . Also, because of C (cid:48)(cid:48) sym-metry, we have t + − ( δn, k ) = e πikδn t − + ( δn, − k ) . (12)The expression (11) guarantees that the 2 × t αα (cid:48) ( n − n (cid:48) , k ) does not contain the Pauli matrix σ = (cid:18) − (cid:19) . This in turn allows us to study the windingnumber of the two Dirac points in H kin by defining Blochstates via the Fourier transform of the hybrid WSs: | φ ± ( q, k ) (cid:105) = 1 √ N (cid:88) n e πiqn | w ± ( n, k g ) (cid:105) , (13)and expressing the kinetic energy operator in this Blochbasis. It is important to emphasize that the states (13) are not kinetic energy eigenstates, but they do satisfyBloch condition as can be seen by acting with the trans-lation operators:ˆ T L | φ ± ( q, k ) (cid:105) = 1 √ N (cid:88) n e πiqn | w ± ( n + 1 , k g ) (cid:105) = e − πiq | φ ± ( q, k ) (cid:105) (14)and ˆ T L | φ ± ( q, k ) (cid:105) = e − πik | φ ± ( q, k ) (cid:105) . (15)As is seen from Eq.(13), the states φ ± ( q, k ) are smoothand periodic functions of q with the period 1. Moreover,because the hybrid WSs were constructed to be continu-ous functions of k and satisfy (4), we also have | φ ± ( q, k + 1) (cid:105) = e ∓ πiq | φ ± ( q, k ) (cid:105) . (16)This means that | φ ± ( q, k ) (cid:105) are Bloch states and carryChern numbers ± b α,q,k = 1 √ N (cid:88) n e − i πqn d α,n,k , (17)we can now express the kinetic energy as H kin = (cid:88) αα (cid:48) = ± (cid:88) q,k t αα (cid:48) ( q, k ) b † α,q,k b α (cid:48) ,q,k (18)= (cid:88) k,q (cid:18) b + ,q,k b − ,q,k (cid:19) † (cid:32) (cid:88) µ =0 n µ ( q, k ) σ µ (cid:33) (cid:18) b + ,q,k b − ,q,k (cid:19) , (19)where t αα (cid:48) ( q, k ) = (cid:80) δn t αα (cid:48) ( δn, k ) e πiqδn . As pointedout above n = 0. From (12) we also find n ( q, k ) = n ( q + k, − k ) ,n ( q, k ) = − n ( q + k, − k ) . (20)Fig. 4 shows the sign of n (red) and n (blue) as afunction of the momentum q and k in the BZ. We seethat the two Dirac points have the same chirality inthat going from, say, ++ to + − , we encircle either one ofthe Dirac nodes clockwise. Naively, this seems to violatethe fermion doubling theorem based on which we expectopposite chirality of the Dirac nodes . However, thistheorem assumes not only that the Hamiltonian H kin ( k )in Eqn. 19 is smooth, but also that it is periodic in themomentum space. Periodicity in q is guaranteed by (13),but not in k as shown by Eqn. 16. Indeed, n would suffera sign change at a step discontinuity if we were to identify k = 0 with k = 1, as can be seen in Fig. 4. The samechirality of two Dirac points characterizes the nontrivialtopology of the narrow bands. It prevents construction ofexponentially localized WSs in both directions if we alsoinsist that each originates within a single valley (i.e. novalley mixing) and with a simple transformation under C T , because in such case the kinetic energy would besmooth, periodic in the BZ, and the σ matrix would beabsent . However, as discussed above, such states canbe constructed if we relax the mentioned requirements. FIG. 4. Chiralities of two Dirac points. The bold red andblue curves show where n and n , defined in Eqn. 19, vanishrespectively, and the colored “ ± ” shows the sign of corre-sponding n or n in the region separated by the bold curves.The two Dirac points K and K (cid:48) are at the intersections of thetwo colored curves and have the same chiralities, as can bedetermined by how n and n change their signs going aroundthe points.
1. Symmetry and Fermion Spectrum
Before proceeding to the detailed calculations, we firstsummarize the impact of various symmetry breaking onthe fermion spectrum and thus provide a qualitative un-derstanding of our results. The kinetic Hamiltonian H kin in Eqn. 19 produces two C T symmetry protected Diracnodes with the same chirality at the corner of the BZ .Without breaking the C T or the translation symme-try, the system is metallic even in the strong couplingregime. Because the x -direction electrical current density j x and the perpendicular electric field E y have oppositeparities under C T transformation, the Hall conductiv-ity, defined by the formula j x = σ xy E y , always vanishesin a C T symmetric system. Breaking C T symmetrycan open a gap in our single flavor model at half fill-ing, corresponding to ν = 3 if the fermions with differentspins or different valleys are assumed to be filled. Thisgapped phase could be either QAH if the masses ( n inEqn. 18) at two nodes are the same, or a topologicallytrivial phase if these two masses are opposite, consistentwith the flipped Haldane model picture of Refs. . Asshown later in the text, our numerical calculation canonly find the QAH phase when C T symmetry is spon-taneously broken, suggesting that the phase with oppo-site masses is not energetically favored by the interac-tions. Furthermore, a C T symmetric stripe phase isalso found to be energetically favored by the interactionsand gapped; as mentioned it must have vanishing Hallconductivity. C. Coulomb interaction energy in the hybridWannier basis
We start from the gate-screened Coloumb interaction,with two metallic gates placed distance ξ above and be- low the TBG,ˆ V = 12 (cid:88) r r (cid:88) µν (cid:32)(cid:88) l V intra ( r − r ) : ˆ ρ lµ ( r )ˆ ρ lν ( r ) : + (cid:88) l (cid:54) = l (cid:48) V inter ( r − r ) : ˆ ρ lµ ( r )ˆ ρ l (cid:48) ν ( r ) : , (21)where V intra ( r ) ( V inter ( r )) is the gate-screened Coulombpotential for two point charges separated by in-plane dis-tance r , and located in the same (different) graphene lay-ers. The graphene layers in the TBG are assumed to beseparated by a small distance d ⊥ of the order of a coupleof carbon lattice spacings . l is the layer index, and µ isthe index combining spin and sublattice degrees freedom(as mentioned, we ultimately study a spinless model, sothis is just for generality). ˆ ρ ( r ) is the charge density at r and : ˆ ρ ˆ ρ : is the normal ordered operator ˆ ρ ˆ ρ . TheFourier transform of such gate-screened Coulomb inter-actions is V intra ( q ) ≈ V inter ( q ) ≈ e π(cid:15) πq tanh qξ qd ⊥ (cid:28)
1. At large momentum qd ⊥ (cid:38)
1, the chargedensity of our single valley model, ρ ( q ), becomes negligi-bly small, and thus, the Coulomb interaction with largemomentum transfer can be safely neglected. (In the twovalley case, it also peaks at the momentum difference be-tween the valleys and there the decrease of V intra/inter ( q )with increasing q makes such terms smaller, see e.g.Ref ). Eqn. 22 is used in all the following analysis andnumerical calculations with ξ set to 10nm.To obtain the projected Coulomb interaction, we firstproject the bare fermion creation and annhilation oper-ator to the constructed hybrid WSs: c † µ ( r ) −→ (cid:88) β = ± (cid:88) k,n w ∗ β,n,k ( µ, r ) d † β,n,k = (cid:88) β = ± (cid:88) k,n w ∗ β, ,k ( µ, r − n L ) d † β,n,k , (23)where d † β,n,k creates the hybrid WS | w β ( n, k g ) (cid:105) withthe wavefunction of (cid:104) r | w β ( n, k g ) (cid:105) = w β,n,k ( r ). Notethat we now absorb the layer index µ and the sublatticeindex apparent in H BM , Eq. (1), into the four component‘spinor’ w β,n,k ( r ). The projected interaction becomesˆ V int = (cid:88) ββ (cid:48) γγ (cid:48) (cid:88) n n n n (cid:88) k k (cid:48) p p (cid:48) J γγ (cid:48) ,n n ,p p (cid:48) ββ (cid:48) ,n n ,k k (cid:48) d † βn k d † γn p d γ (cid:48) n p (cid:48) d β (cid:48) n k (cid:48) (cid:88) i ∈ Z δ k + p ,p (cid:48) + k (cid:48) + i . (24)The term in the last line ensures the momentum conser-vation along g . It is worth emphasizing that the ob-tained interaction in Eqn. 24 has been numerically foundto be sizable even if the difference of unit cell indices | n i − n j | ≤ i, j ). Different from thewavefunction of the lowest Landau level (LLL), the con-structed hybrid Wannier state, shown in Fig. 2, containstwo peaks along L , leading to significant overlap be-tween two hybrid Wannier states with consecutive unitcell indices. Correspondingly, the projected Coulomb in-teractions decays exponentially only when | n i − n j | ≥ III. DMRG
We consider a system with the size of N L × N L and choose the open boundary condition along L andanti-periodic boundary condition along L . Therefore,the momentum indices of the hybrid WSs | w ± ( n, k g ) (cid:105) take the values: k = i + N with i = 0 , , · · · , N − n = − N , − N , · · · , N . Since we study the quasi-1D system with DMRG, the hybrid WSs are arrangedin a one-dimensional chain with each site indexed as i + nN . Also, each site contains two hybrid WSs, la-beled by β = ±
1. In the DMRG calculation, N = 6,and N = 30, the bond dimension set to be 2000 andthe truncation error no more than 10 − . Calculationswere performed using the ITensor Library to study theground states at the half filling of the spinless one-valleymodel, i.e. at the average occupation of one particle perunit cell. The Hamiltonian studied in DMRG includesonly the electron-electron interactions with the kineticterms H kin neglected. Because of the complicated formof the projected interaction, ITensor produces a ratherlarge matrix product operator (MPO), 3 − w and w in the BM model (1).When w = 0, i.e. the interlayer intra-sublattice hop-ping vanishes, the system is in the chiral limit withtwo Chern bands located on different sublattices. Conse-quently, the ground state has been shown to be the QAHstate . As w /w increases, the two Chern bands startto spatially overlap, leading to the scattering among themand frustrating the QAH. Nevertheless, because QAH isa gapped phase, it is stable with respect to a small in-crease of w /w from the chiral limit. Whether or notit collapses and a different state is favored as w /w reaches ∼ .
83 is the purpose of our DMRG calculation.We should note, however, that the increased propensitytowards a many body insulating state with increasing w /w could also be intuited from the shapes of thephases of Wilson loop eigenvalues. In Fig.(3) we see thatthey progressively flatten, suggesting that a good corre-lation hole can be built when the two hybrid WSs are FIG. 5. The order parameter (cid:104) n (cid:105) given by Eqn. 26 of theground state as obtained by DMRG. |(cid:104) n (cid:105)| = 1 when theground state is QAH for w /w (cid:46) .
7. However, (cid:104) n (cid:105) ≈ w /w (cid:38) .
8, suggesting the vanishing Hall conductivityin the system. coherently (and equally) distributed among the Chern+1 and Chern − C T and possibly insulate.This intuitive picture turns out to be consistent withour DMRG calculation. With N up to 6, DMRG findsQAH as the ground state when w /w ≤ .
7, i.e. themany-body ground state turns out to be a product stateof the hybrid WSs | w β,n,k (cid:105) with the same Chern index β : | Ψ (cid:105) GS = (cid:89) n,k d † + ,n,k |∅(cid:105) or (cid:89) n,k d †− ,n,k |∅(cid:105) . (25)With w /w = 0 . w /w = 0 .
85, DMRG alwaysproduces a non-QAH state with translation symmetrybreaking even if the initial state is set to be the QAHstate. In this parameter regime, however, the DMRGcalculation does not result in a fully converged groundstate, in that the details of the final state are sensitiveto the choice of the bond dimension; this is despite theentanglement entropy through the middle bond never go-ing above 0 .
82. However, several interesting propertiesare found to be common among all the obtained stateswhich we now discuss.
Vanishing (cid:104) n (cid:105) : To illustrate the difference betweenthe QAH states and the state obtained when w /w ≥ .
8, we define the order parameter (cid:104) n (cid:105) : (cid:104) n (cid:105) = 1 N (cid:88) n,k (cid:104) d † + ,n,k d + ,n,k − d †− ,n,k d − ,n,k (cid:105) GS , (26)where N = N × N is the total number of particles inthe system. We found that (cid:104) n (cid:105) changes dramaticallywhen w /w is between 0 . .
8. When w /w ≤ . (cid:104) n (cid:105) ≈ ±
1, consistent with the QAH states described inEqn. 25. With w /w ≥ . (cid:104) n (cid:105) quickly drops to 0,suggesting that DMRG gives a topologically trivial statewith vanishing Hall conductivity. Product State : On each site labeled by the two indices n and k , we also find that the fermion occupation numberin the ground state is almost 1 when 0 ≤ w /w ≤ . - - - - - (a) - - - - - (b) - - - - - (c) - - - - - (d) FIG. 6. The fermion correlation between different sites with w /w = 0 .
85 in the state obtained by DMRG. We used 6 k-pointsand 31 n-points with the bond dimension of 2000 and the truncation error of 10 − . one, and two particles on several typical sites ( n, k ) when w /w = 0 .
85, where the particle number operator ˆ N n,k on site ( n, k ) is d † + ,n,k d + ,n,k + d †− ,n,k d − ,n,k . The proba-bility of having 0, 1, and 2 particles on the site ( n, k ) iscalculated with the following formula: P ( ˆ N n,k = 0) = 12 (cid:104) (1 − ˆ N n,k )(2 − ˆ N n,k ) (cid:105) GS (27) P ( ˆ N n,k = 1) = (cid:104) ˆ N n,k (2 − ˆ N n,k ) (cid:105) GS (28) P ( ˆ N n,k = 2) = 12 (cid:104) ˆ N n,k ( ˆ N n,k − (cid:105) GS . (29)As shown in the Table. I, the probability of having 0 or2 particles on each site are negligible, suggesting that theDMRG produced state can be well approximated by theproduct of one-particle state on each site. To further jus-tify this statement, we calculated the equal-time fermioncorrelation between different sites, i.e. (cid:104) d † α,n,k d β,n (cid:48) ,k (cid:48) (cid:105) GS .This correlation is shown in Fig. 6 when w /w = 0 . n, k ) = (0 , ), and listthe absolute value of the fermion correlation with various( n (cid:48) , k (cid:48) ). When ( n, k ) is fixed to be (0 , ), the correla-tion is also listed in Fig. 6(b) and 6(d). Although this ( n, k ) (0 ,
112 ) (0 ,
14 ) (0 ,
512 ) (0 ,
712 ) (0 ,
34 ) (0 , P ( ˆ N n,k = 0) 0 .
038 0 .
017 0 .
018 0 .
018 0 .
018 0 . P ( ˆ N n,k = 1) 0 .
922 0 .
966 0 .
965 0 .
964 0 .
965 0 . P ( ˆ N n,k = 2) 0 .
040 0 .
017 0 .
017 0 .
018 0 .
017 0 . P of having zero, one, or two par-ticles on several typical sites labeled by ( n, k ) when w /w =0 .
85. The small probability of having zero or two particlesshows the negligible charge fluctuations away from the one-particle occupancy. correlation increases with k close to 0 or 1, the off-sitecorrelation is generally found to be tiny. In addition, wefound that that the correlation exponentially decays asa function of | n − n (cid:48) | , suggesting that this is a gappedphase.Overall, the dominant one-particle occupancy and tinyoff-site correlation of the DMRG produced state suggestthat the DMRG produced wavefunction can be well ap-0proximated by the following formula: | Ψ GS (cid:105) ≈ (cid:89) n,k (cid:16) u n,k d † + ,n,k + v n,k d †− ,n,k (cid:17) |∅(cid:105) (30)with | u n,k | + | v n,k | = 1 for all the sites labeled by ( n, k ). Phase of (cid:104) d † + ,n,k d − ,n,k (cid:105) : We also found that the phaseof (cid:104) d † + ,n,k d − ,n,k (cid:105) in the non-QAH phase can be describedby a functionarg (cid:16) (cid:104) d † + ,n,k d − ,n,k (cid:105) (cid:17) ≈ arg (cid:18) v ( n, k ) u ( n, k ) (cid:19) ≈ f ( n ) π , where f ( n ) = ±
1. Obviously, the phase of (cid:104) d † + ,n,k d − ,n,k (cid:105) depends on the choice of the phase of the constructedhybrid WSs, and thus is not U (1) gauge invariant. Oncethe U (1) phase of the hybrid WSs is fixed , f ( n ) is foundto depend only on n but does not show any regular pat-tern in DMRG produced final state. In the next section,we will see that the magnitude of this phase, π , is re-produced by minimizing the (cid:104) H (cid:105) with the trial state inEqn. 30.As stated previously, the DMRG produced state doesnot converge and is very sensitive to the bond dimension.Additionally, we found that the C T local order param-eter (cid:104) d † + ,n,k d + ,n,k − d †− , − n,k d − , − n,k (cid:105) strongly depends on n and does not show any spatially periodic pattern. Thismay come from the strong competition between variouslow energy states, such as QAH and other C T symmet-ric states. As suggested by the drop of (cid:104) n (cid:105) in Fig. 5 andanalysis in the following sections, the system undergoesa first order phase transition from QAH to a non-QAHstate when w /w reaches approximately 0 .
8. Since theglobal order parameter (cid:104) n (cid:105) vanishes in the non-QAHstate, we suspect that C T symmetry may be still locallyconserved in this state, leading to zero Hall conductivity.When w /w becomes slightly larger than 0 .
8, the systemis in the non-QAH regime but close to the phase transi-tion point. As a consequence, these states are still almostdegenerate, leading to strong competition among them.Furthermore, the QAH state, as discovered in DMRG, isfavored at the open ends of the system. Therefore, thesmall advantange of the non-QAH states in the bulk andthe superiority of the QAH state at the boundary maydrive the system into the intermediate phase without reg-ular spatial patterns in the DMRG produced state.
IV. ANALYSIS OF THE TRIAL STATE INEQN.(30)A. Ground State
Because the trial ground state (Eqn. (30)) suggested byDMRG is a product state at each site ( n, k ), it is straight-forward to analyze its energy variationally by implement-ing the Wick’s theorem. This allows us to increase N ,the number of k -points, which is limited to N ≤ E N = (cid:104) Ψ GS | λ ˆ H kin + ˆ V int | Ψ GS (cid:105) , (31)with the constraint | u ( n, k ) | + | v ( n, k ) | = 1 for every n and k , and allowing λ to increase continuously form 0.Furthermore, we seek a solution periodic in the unit cellindex n , i.e. u ( n, k ) = u ( n + n p , k ) , v ( n, k ) = v ( n + n p , k ) , where n p is the period. The C T symmetric state satis-fies the constraint that u ( n, k ) = v ∗ ( − n, k ) e iθ ( n,k ) v ( n, k ) = u ∗ ( − n, k ) e iθ ( n,k ) (32)and the C (cid:48)(cid:48) symmetric state should satisfy u ( n, k ) = v ( n, − k ) e iθ (cid:48) ( n,k ) v ( n, k ) = u ( n, − k ) e iθ (cid:48) ( n,k ) . (33)We then optimize the energy allowing both C T and C (cid:48)(cid:48) symmetries to be broken and allowing the period n p tobe as large as 8. Numerically, we found three types ofthe solutions with the symmetry listed in Table. II. Solution Translation C T C (cid:48)(cid:48) C T broken Conserved Broken Broken C T Conserved Conserved Broken if λ (cid:46) . λ (cid:38) . C T Broken Conserved BrokenStripe ( n p = 2) ( T L C (cid:48)(cid:48) Conserved)TABLE II. The symmetry of three obtained solutions with thetrial states given by Eqn. 30, where λ , defined in Eqn. 31 is thescaling factor of the kinetic terms. Whether the one fermionspectrum of the states is gapless or gapped is specified later inthe text and figures; we only focus on the symmetry breakinghere. For N = 6, the C T broken state is found to be theground state for w /w (cid:46) . C T symmetric state with brokentranslation symmetry is found to be the ground statefor w /w ≈ .
85, corresponding to the C T -symmetricperiod-2 stripe. Although this stripe state breaks C (cid:48)(cid:48) symmetry, the combination of translation along L and C (cid:48)(cid:48) , i.e. T L C (cid:48)(cid:48) is still conserved. This result is consistentwith the one obtained using DMRG for the same value of N , in that the (cid:104) n (cid:105) vanishes at roughly the same valuesof w /w . Moreover, up to a π , the k - dependence of thearg( v ( n, k ) /u ( n, k )) obtained by the variational methodis the same as the one by DMRG. The k and n depen-dence of this phase will be more thoroughly discussed inthe next subsection.1 FIG. 7. The energy of various correlated states with the trial function in Eqn. 30. The energies are normalized by U = e / (4 π(cid:15)L m ), where L m = | L | . Left: the energies of four different states: C T broken state, C T symmetric nematic state, C T symmetric period-2 stripe state, and the semi-metal; the semi-metal is defined as the (non-interacting) state obtained bydiagonalizing the kinetic energy only. The results are plotted vs (cid:104) E K (cid:105) SM ≡ λ (cid:104) ˆ H kin (cid:105) SM as λ increases from 0 and (cid:104) . . . (cid:105) SM isthe expectation value in the semi-metal state. Right: the energies of three nearly degenerate states, and the transition betweenthem are marked with colored stars. As seen here, at w /w = 0 .
85, these three states are nearly degenerate. In contrast, at w /w = 0 .
3, the near degeneracy is lifted in favor of QAH whose energy is 0 . U ≈ C T symmetricstates (see the corresponding Fig. S2 in the Appendix ). We can also obtain another variational state by im-posing the translational symmetry ( n p = 1) and C T invariance. We refer to this state as C T -nematic. Thebest variational energies of these states are compared inFig. 7 as we change λ , plotting the result per particleas a function of the kinetic energy in the half-filled non-interacting semi-metallic state, (cid:104) E K (cid:105) SM ≡ λ (cid:104) ˆ H kin (cid:105) SM ,in units of U = e N/ (4 π(cid:15)L m ). The interaction constant U ≈ L m is the moire superlattice con-stant and (cid:15) is the dielectric constant of hBN. Becausethe simple non-interacting semi-metal state diagonalizesthe kinetic energy –and therefore optimizes it– we in-clude in the plot the expectation value of λ ˆ H kin + ˆ V int inthis state. Although, the non-interacting semi-metal isclearly not competitive in the range of the parameters ofinterest to us, its expected energy does provide us with ameasure of near degeneracy among the competing states.For N = 16 and w /w = 0 .
85 the Fig. 7 thus com-pares the energies of three different variational states: C T broken state, C T symmetric period-2 stripe phase,and C T symmetric nematic state. As seen, all thesethree states are nearly degenerate. As a measure ofhow close the energies of the three states are, we dividethe energy difference between the competitive states bythe energy difference between the ground state and thenon-interacting semi-metal. Without the kinetic termsin the Hamiltonian, we find that the normalized en-ergy difference between QAH and C T nematic state isonly 0 . U / (0 . U ) ≈ .
02, and the normalized en-ergy difference between QAH and C T stripe state is0 . / . ≈ .
06, in favor of QAH. As seen in Fig.7, theenergies of the competitive states are even closer whenthe kinetic energy terms are included.We find that the ground state is always translation-ally invariant. Additionally, the C T symmetry is bro-ken when the kinetic energy (cid:104) E K /N (cid:105) SM < . U , and fully gapped when (cid:104) E K /N (cid:105) SM < . U , suggestingthat the state we found is QAH for small kinetic en-ergy and turns into an anomalous Hall metal when0 . < |(cid:104) E K (cid:105) SM | /N U < .
4. It eventually evolves intoa normal metal with vanishing Hall conductivity when |(cid:104) E K (cid:105) SM | /N U > .
4. Nevertheless, the energies of thetwo C T symmetric states in Table. II are very close tothe energy of the QAH state in all the parameter regimeswe have calculated. As we will discuss later, the energyof the C T symmetric states can be further lowered byimproving the form of the variational states. This neardegeneracy necessitates the inclusion of all three differ-ent states as the candidates for the ground state at oddinteger filling.Because the anomalous Hall state seems to have beenruled out in experiments on magic angle TBG at ν = 3without the alignment with hBN, and because as we willsee below C T period-2 stripe state can be fully gapped,we consider it as a candidate for the Chern-0 insulat-ing state experimentally observed at ν = 3. QAH state,on the other hand, can be favored by breaking C sym-metry, and thus is the state discovered at the same fill-ing but aligning the system with the hBN substrate. Inaddition, the C T nematic state, being gapless, simul-taneously breaks C rotation symmetry and possessesthe interesting pattern of the Landau level degeneracy.Therefore, after including the spin and valley degree offreedom, we propose the C T nematic state as a candi-date for the gapless state at the charge neutrality point(CNP).2 B. C (cid:48)(cid:48) Symmetry
The QAH states can be approximated as | Ψ QAH (cid:105) ≈ (cid:89) n,k d † + ,n,k |∅(cid:105) or (cid:89) n,k d †− ,n,k |∅(cid:105) . (34)Since the hybrid states transform as Eqn. 6 under C (cid:48)(cid:48) ,this state obviously breaks C (cid:48)(cid:48) symmetry (in addition to,of course, C T ).If the C T symmetric state is translationally invariantor has the period of 2 unit cells, the u ( n, k ) and v ( n, k )in the trial wavefunction Eq. 30 can be written as u ( n, k ) = 1 √ e iφ ( n,k ) v ( n, k ) = 1 √ e − iφ ( n,k ) . (35)If the state is further C (cid:48)(cid:48) symmetric, e iφ ( n,k ) = ± e − iφ ( n, − k ) = ⇒ φ ( n, k ) + φ ( n, − k ) = 0 or π . (36) (a) (b) FIG. 8. (a) φ ( n, k ), defined in Eqn. 35, in the C T nematicstate with the trial function in Eqn. 30 and N = 16. λ is thescaling factor of the kinetic energy, defined in Eqn. 31. (b) φ ( n, k ) in the C T stripe state when n is even. If n is odd, φ ( n, k ) can be obtained from the relation in Eqn. 39. Fig. 8 illustrates the k dependence of the phase φ ( n, k )in the two C T symmetric states. In the C T nematicstate, due to the translation symmetry, φ ( n, k ) is inde-pendent of n . With the interactions only, λ , defined in Eqn. 31, vanishes, and as the Fig. 8(a) shows, the statecan be approximated as | Ψ N (cid:105) ≈ (cid:89) n,k √ (cid:16) e i π d † + ,n,k + e − i π d †− ,n,k (cid:17) |∅(cid:105) . (37)Obviously, this state breaks the C (cid:48)(cid:48) symmetry. With in-creasing λ , the phase φ ( k ) becomes smaller, and eventu-ally vanishes when λ (cid:38) .
8, and thus the C (cid:48)(cid:48) symmetryis recovered. In particular, for the BM model includingboth the interaction and kinetic terms without any scal-ing, λ = 1, and thus the state can be approximated as | Ψ N (cid:105) ≈ (cid:89) n,k √ (cid:16) d † + ,n,k + d †− ,n,k (cid:17) |∅(cid:105) . (38)Our calculation shows that the C T stripe state is al-ways invariant under T L C (cid:48)(cid:48) transformation. This leadsto the relation φ ( n, k ) = − φ ( n + 1 , − k ) (39)relating the phase φ ( n, k ) with even n and the phase withodd n . When λ vanishes, Fig. 8(b) shows that the wave-function of the stripe state can be approximated as | Φ sN (cid:105) ≈ (cid:89) m,k √ (cid:16) e i π d † + , m,k + e − i π d †− , m,k (cid:17) × √ (cid:16) e − i π d † + , m +1 ,k + e i π d †− , m +1 ,k (cid:17) |∅(cid:105) . (40)Obviously, φ ( n, k ) (cid:54) = − φ ( n, − k ) and thus the C (cid:48)(cid:48) sym-metry is broken in this state.Similar to the nematic phase, the magnitude of φ ( n, k )decreases with increasing λ . When λ (cid:38)
24, the stripestate satisfies the relation φ ( n, k ) = − φ ( n, − k ) andthus becomes C (cid:48)(cid:48) symmetric. C. Excitation spectrum
If the trial state (Eqn. 30) does not break the transla-tion symmetry, it can be written as | Ψ N (cid:105) = (cid:89) q,k (cid:16) u ( k ) b † + ,q,k + v ( k ) b †− ,q,k (cid:17) |∅(cid:105) , (41)where b ± ,q,k is defined in Eqn. 17 as the Fourier transformof the fermion operator d ± ,n,k .To construct the one particle and hole excited states,we delocalize a linear combination of d † + ,n,k and d †− ,n,k : | Ψ N +1 ( q, k ) (cid:105) = (cid:16) v ∗ ( k ) b † + ,q,k − u ∗ ( k ) b †− ,q,k (cid:17) | Ψ N (cid:105) (42) | Ψ N − ( q, k ) (cid:105) = ( u ∗ ( k ) b + ,q,k + v ∗ ( k ) b − ,q,k ) | Ψ N (cid:105) . (43)The variational energies of these excited states are E N ± ( q, k ) = (cid:104) Ψ N ± ( q, k ) | H | Ψ N ± ( q, k ) (cid:105) , (44)3 (a) (b) (c) FIG. 9. The direct gap ∆ d ( q, k ) of three states with the trialstate in Eqn. 30 with w /w = 0 .
85. The kinetic energyis set to be (cid:104) E K (cid:105) SM = 0. The energies are normalized by U = e / (4 π(cid:15)L m ). and the gap is given by∆ = min( E N +1 ( q, k ) − E N ) − max( E N − E N − ( q (cid:48) , k (cid:48) )) . (45)Fig. 9 illustrates the direct gap ∆ d , defined as∆ d ( q, k ) = ( E N +1 ( q, k ) − E N ) − ( E N − E N − ( q, k )) , in the BZ when w /w = 0 .
85 and the kinetic terms are set to be 0. Interestingly, we found the C T sym-metric nematic phase is gapless with nodes around Γ .The robustness of the nodes has been discussed in Ref. and additional properties will be presented in the nextsection. The gap can be opened by breaking C T sym-metry. A typical example of this case is the QAH statehaving a gap of order U as shown in Fig. 9(a). In thenext subsection, we will show that a gapped C T sym-metric state can be obtained by breaking the translationsymmetry. FIG. 10. The single fermion excitation gap of three nearlydegenerate states. The energies are normalized by U = e / (4 π(cid:15)L m ). As mentioned before, we obtained three different typesof solutions. While the former two are translationallyinvariant, the last one is C T symmetric but breaks thetranslation symmetry with the period of 2. Consequently,the coefficients u ’s and v ’s in the trial state satisfy u (2 n, k ) = u (0 , k ) v (2 n, k ) = v (0 , k ) (46) u (2 n + 1 , k ) = u (1 , k ) v (2 n + 1 , k ) = v (1 , k ) . (47)The corresponding stripe state can be described by thefollowing wavefunction: | Ψ sN (cid:105) = (cid:89) n,k (cid:89) m =0 (cid:16) u ( m, k ) d † + , n + m,k + v ( m, k ) d †− , n + m,k (cid:17) |∅(cid:105) . (48)Similar to the translationally invariant state, we can alsodescribe this stripe state in the basis of Chern Blochstates. For notational convenience, we introduce anotherset of fermion operators: f ± ,q,k, = 1 √ (cid:16) b ± ,q,k + b ± ,q + ,k (cid:17) (49) f ± ,q,k, = 1 √ (cid:16) b ± ,q,k − b ± ,q + ,k (cid:17) , (50)with 0 ≤ q < . Up to an overall phase, the stripe statein Eqn. 48 can be written as | Ψ sN (cid:105) = (cid:89) k ∈ [0 , q ∈ [0 , /
2) 1 (cid:89) m =0 (cid:16) u ( m, k ) f † + ,q,k,m + v ( m, k ) f †− ,q,k,m (cid:17) |∅(cid:105) . (51)4Similar to the case of the translationally invariant groundstate, the one particle and hole excited states are builtwith delocalized linear combination of fermion operators d †± ,n,k : | Ψ sN +1 ,m ( q, k ) (cid:105) = (cid:16) v ∗ ( m, k ) f † + ,q,k,m − u ∗ ( m, k ) f †− ,q,k,m (cid:17) | Ψ sN (cid:105) (52) | Ψ sN − ,m ( q, k ) (cid:105) = ( u ∗ ( m, k ) f + ,q,k,m + v ∗ ( m, k ) f − ,q,k,m ) | Ψ sN (cid:105) . (53)To obtain the fermion spectrum, we construct the two2 × (cid:0) H sN ± ( q, k ) (cid:1) m m = (cid:104) Ψ sN ± ,m ( q, k ) | H | Ψ sN ± ,m ( q, k ) (cid:105) . (54)The energy and the wavefunction of the electron and holeexcites states are obtained by diagonalizing these ma-trices H sN ± ( q, k ). At each momentum, we obtain twoeigenvalues E sN ± , ( q, k ) and E sN ± , ( q, k ). The gap ofthe period-2 stripe ground state is given by∆ = (cid:0) min( E sN +1 , ( q, k ) , E sN +1 , ( q, k )) − E sN (cid:1) − (cid:0) E sN − min( E sN − , ( q (cid:48) , k (cid:48) ) , E sN − , ( q (cid:48) , k (cid:48) )) (cid:1) . (55)Fig. 10 shows the magnitude of the gap in the C T sym-metric stripe state. This gap is found to be (cid:38) U atvanishing kinetic energy, but decreases with increasingkinetic energy. It vanishes before evolving into the C T nematic phase. The opening and closing of this gap willbe discussed in much more detail in the next section,where we expose the non-Abelian topological aspectsof this process.Fig. 9(c) plots the direct gap ∆ d , defined as∆ d ( q, k ) = (cid:0) min( E N +1 , ( q, k ) , E N +1 , ( q, k )) − E N (cid:1) − (cid:0) E N − min( E N − , ( q, k ) , E N − , ( q, k )) (cid:1) , (56)making it is obvious that ∆ d ( q, k ) = ∆ d ( q + , k ) due tothe breaking of the translation symmetry with the periodof 2. V. GENERALIZED TRIAL STATESA. Translationally Invariant State
The trial function in Eqn. 41 is not the most generalform for the translationally invariant state, as the coef-ficients u ’s and v ’s are independent of the momentumcomponent q . This comes from the complete absence ofthe correlations between hybrid WSs on different sites inour trial state ( Eqn. 30 ). To improve the trial state, weconsider the following wavefunction | Ψ N (cid:105) = (cid:89) q,k (cid:16) u ( q, k ) b † + ,q,k + v ( q, k ) b †− ,q,k (cid:17) |∅(cid:105) , (57) FIG. 11. The energies of three nearly degenerate states at w /w = 0 .
85 with the more general trial wavefunction givenby Eqn. 57, except for the stripe state which is given by thewavefunction (Eqn. D1) in the Appendix . where u and v depend on both q and k and satisfy | u ( q, k ) | + | v ( q, k ) | = 1. If this state is C T sym-metric, u ’s and v ’s should have the same magnitude,i.e. | u ( q, k ) | = | v ( q, k ) | = 1 / √ E N ( u, v ) = (cid:104) Ψ N | H | Ψ N (cid:105) with respect to u ’s and v ’s at various momenta. Themomentum mesh in the BZ is chosen to be q = i + 1 / N and k = j + 1 / N with i = 0 , , · · · , N − j = 0 , , · · · , N − N and N are taken to be 16. Our calculations still findtwo different solutions, C T broken and C T nematicsolution. Compared with the trial state in Eqn. 30, bothsolutions have lower energies. Interestingly, the C T ne-matic state now has a lower energy than the C T brokenstate, although the difference between these two solutionsare again tiny. The C T broken solution is obtained bysearching a local minimum near the product state with u ( q, k ) = 1 and v ( q, k ) = 0.The fermion spectrum is also calculated in the sameway as shown by Eqn. 42 – 45. Fig. 12 shows the directgap of the C T broken state inside the BZ. With smallkinetic energy, this state is fully gapped over the wholeBZ. The global C T breaking order parameter, definedas (cid:104) n (cid:105) = 1 N (cid:88) q,k (cid:104) Ψ N | b † + ,q,k b + ,q,k − b †− ,q,k b − ,q,k | Ψ N (cid:105) , (58)is found to be almost 1 in this state. Therefore, it is iden-tified as the QAH state. Interestingly, the gap minimumis always located at the Γ point. This is dramatically dif-ferent from the non-interacting state, in which the gap islargest at Γ and closes at K and K (cid:48) .In contrast, the C T nematic state is always gapless,and thus, can never be the experimentally observed in-sulating phase at ν = 3. The plot of the direct gap inthe BZ in Fig. 13 has shown a node (or nodes) either at,5 - - - - - FIG. 12. The direct gap ∆ d ( q, k ) of the QAH state. Thisstate is given by Eqn. 57. or very close to, Γ point if the kinetic terms are set to be0. Additionally, the direct gap quickly increases to ∼ U once the momentum is away from Γ. To understand theproperties of nodes in the C T nematic phase and thegap opening in the C T broken phase, we consider theself-consistent equations obtained as δ (cid:104) Ψ N | H | Ψ N (cid:105) δu ∗ ( q, k ) = E ( q, k ) u ( q, k ) (59) δ (cid:104) Ψ N | H | Ψ N (cid:105) δv ∗ ( q, k ) = E ( q, k ) v ( q, k ) . (60)This is equivalent to the minimization of (cid:104) Ψ N | H | Ψ N (cid:105) .The E ( q, k ) is the Lagrange multiplier, needed because ofthe constraints | u ( q, k ) | + | v ( q, k ) | = 1 for each ( q, k ).This equation can be written in the matrix form, H eff ( q, k ) (cid:32) u ( q, k ) v ( q, k ) (cid:33) = E ( q, k ) (cid:32) u ( q, k ) v ( q, k ) (cid:33) , (61)where H eff ( q, k ) is a Hermitian 2 × u ’s and v ’s. E ( q, k ) is an eigenvalue of this ma-trix, and, by Koopman’s theorem, the direct gap ∆ d ( q, k )is calculated as the difference between the two eigenval-ues of the matrix H eff ( q, k ).The Hermitian matrix H eff is obtained as( H eff ( q, k )) ββ (cid:48) = (cid:104) φ β ( q, k ) |F| φ β (cid:48) ( q, k ) (cid:105) , (62)where | φ β ( q, k ) (cid:105) is the Bloch state defined in Eqn. 13,and F is a Hermitian operator independent of the mo-mentum . Note that the Bloch state | φ ± ( q, k ) (cid:105) has thewinding number of ± F has no winding number, the matrix element (cid:0) H eff ( q, k ) (cid:1) + − has the winding number of − . As a consequence, it contains, at least, either aquadratic node or two Dirac nodes inside the BZ . For a C T symmetric state, ( H eff ( q, k )) ++ = ( H eff ( q, k )) −− , and thus a node appears as long as the off-diagonal ma-trix element ( H eff ( q, k )) + − = ( H eff ( q, k )) ∗− + vanishes.This state, therefore, must be a gapless state. On theother hand, the QAH state breaks the C T symme-try, and the two diagonal elements ( H eff ( q, k )) ++ and( H eff ( q, k )) −− become unequal, leading to the openingof a gap. - - - - - (a) - - - - - (b) FIG. 13. The direct gap ∆ d ( q, k ) for C T nematic statewith w /w = 0 .
85 at (a) (cid:104) E K (cid:105) SM = 0 and (b) (cid:104) E K (cid:105) SM = − . U . It is interesting to investigate phenomenological conse-quences of the nodes in this state. The arguments abovesuggest the fermion spectrum contains either a quadraticband touching point or two close Dirac nodes with lineardispersion. As shown in Fig. 13(a), our numerical calcu-lation suggests the existence of a quadratic node whenthe weight of kinetic terms (cid:104) E K (cid:105) SM vanishes. In thiscase, the Landau levels are doubly degenerate at zeroenergy and non-degenerate at all other energy levels .It is worth noting that such Landau level degeneracy,plus the degeneracy brought by valley and spin degreesof freedom, produces the filling patterns of Landau fan6observed in the experiments ν = ± , ± , ± , . . . . Weshould also point out that the possibility of two very closeDirac nodes in this system cannot be ruled out due to theinsufficient resolution of the momentum mesh. But thisdoes not affect such Landau level filling pattern as long asthe two Dirac nodes are close enough and the magneticfield is not too small. Moreover, Fig. 13(a) illustratesthat the density of states monotonically increases as afunction of energy as the filling changes away from theneutrality point, a feature also qualitatively consistentwith experiments .Besides the Landau fan pattern, the Fig. 13(a) alsoshows that the fermion spectrum in this state breaks C symmetry and therefore we refer to it as the C T nematicphase. To understand why C symmetry is broken, weexpand the effective Hamiltonian around Γ. It is moreconvenient to express the momenta as the complex num-bers and introduce z = k x − ik y . Up to the quadraticterms, the effective Hamiltonian can be approximated as H eff ( k ) = (cid:32) λ ( z − z )( z − z ) λ ∗ ( z − z ) ∗ ( z − z ) ∗ (cid:33) , (63)where λ is a complex constant. It is obvious that the theHamiltonian contains two nodes at z and z with thesame chirality. Under C rotation, the two Bloch statesat Γ transform trivially , but the momentum z = k x − ik y obtains a phase of e − i π/ . Therefore, H eff ( k )must break the C rotation, and the resulting gaplessphase is nematic.Another interesting feature is the location and robust-ness of the node close to Γ. Fig. 14 illustrates the di-rect gap of the C T nematic state when w /w = 0 . ∼ ), this C T nematic state contains a quadraticnode at or very close to Γ when the weight of the ki-netic term vanishes, but evolves into two well-separatedDirac nodes with a small weight of kinetic term. At thelarger ratio, w /w = 0 .
85, and the same weight of ki-netic terms (cid:104) E K (cid:105) SM = 0 . U , no splitting of nodes orthe movement of nodes can be identified in Fig. 13(b).It seems that the nodes are trapped in a deep potentialwell at Γ with large w /w . This may be related withthe steep slope of the Wilson loop eigenvalue at k closeto 0 with w /w (cid:38) .
8, or more explicitly, the high peakof the Berry curvature of ± B. C T symmetric period- stripe state In this subsection, we investigate the properties of the C T symmetric period-2 stripe phase. Here, we only con-sider a subset of the such states, | Ψ s (cid:105) , that can be writtenin the form of a product states, so that the Wick’s theo-rem applies. The most general form of such states can bewritten in a relatively simple expression with four free pa-rameters specifying two points on an abstract unit sphere - - - - - (a) - - - - - (b) FIG. 14. The direct gap ∆ d ( q, k ) of the C T nematic statewith w /w = 0 . (cid:104) E K (cid:105) SM = 0 and (b) (cid:104) E K (cid:105) SM = − . U . Note the splitting of the node near Γ. at each momentum . By minimizing E s = (cid:104) Ψ s | ˆ H | Ψ s (cid:105) ,we found a local minimum of E s with the state | Ψ s (cid:105) breaking the translation symmetry. As shown in Fig. 11,this state has the energy E s slightly higher than the C T nematic phase, and still lower than the QAH state.But the energy differences between this stripe state andother states are found to be very small, no more than0 . U ≈ . ν = 3. With some degree of valley mixing itmay also be possible to further lower the energy of sucha state.Fig. 15 shows the momentum dependent direct gap,defined in Eqn. 56 but calculated with general trial wave-function Eqn. D1. It is obvious that the gap is periodic∆ d ( q, k ) = ∆ d ( q + , k ). To have a deeper understand-7 - - - - - FIG. 15. The direct gap ∆ d ( q, k ) of the C T stripe phase at E K = 0. It is clear that this state does not contain any nodes. ing of the fermion spectrum in this phase, we considerthe self-consistent equations, which can be written as H seff ( q, k ) u ( q, k ) u ( q + , k ) v ( q, k ) v ( q + , k ) = E ( q, k ) u ( q, k ) u ( q + , k ) v ( q, k ) v ( q + , k ) , (64)with H seff being a 4 × | η ( q, k ) (cid:105) = (cid:18) | φ + ( q, k ) (cid:105) , | φ + ( q + 12 , k ) (cid:105) , | φ − ( q, k ) (cid:105) , | φ − ( q + 12 , k ) (cid:105) (cid:19) . (65)The matrix element of the effective Hamiltonian of thestripe state can then be written as (cid:0) H seff ( q, k ) (cid:1) ij = (cid:104) η i ( q, k ) |F s | η j ( q, k ) (cid:105) , (66)where F s is an operator independent of the momentum .In addition, C T symmetry leads to the form H seff ( q, k ) = (cid:15) ( q, k ) δ ( q, k ) ∆ ( q, k ) ∆ ( q, k ) δ ∗ ( q, k ) (cid:15) ( q + , k ) ∆ ( q, k ) ∆ (cid:48) ( q, k )∆ ∗ ( q, k ) ∆ ∗ ( q, k ) (cid:15) ( q, k ) δ ∗ ( q, k )∆ ∗ ( q, k ) (∆ (cid:48) ( q, k )) ∗ δ ( q, k ) (cid:15) ( q + , k ) (67)where ∆ (cid:48) ( q, k ) = ∆ ( q + , k ). The ground state is ob-tained by solving the eigenvalue problem of the Hermi-tian matrix H seff . Being a 4 × E ( q, k ) ≤ E ( q, k ) ≤ E ( q, k ) ≤ E ( q, k ). The direct gap ∆ d ( q, k ) can be calculated as E ( q, k ) − E ( q, k ). With the gauge chosen in Eqn. 16, we find the matrixelements satisfy the following conditions :∆ ( q + 12 , k ) = ∆ ( q, k )∆ ( q + 12 , k ) = ∆ (cid:48) ( q, k )∆ (cid:48) ( q + 12 , k ) = ∆ ( q, k )∆ ( q, k + 1) = e i πq ∆ ( q, k )∆ (cid:48) ( q, k + 1) = e i πq ∆ (cid:48) ( q, k )∆ ( q, k + 1) = − e i πq ∆ ( q, k ) (cid:15) ( q + 1 , k ) = (cid:15) ( q, k ) . (68)Based on the boundary conditions above, it is easy to seethat the matrix element ∆ ( q, k ) has the winding numberof 1 around the stripe BZ (0 ≤ q < and 0 ≤ k <
1. Gapped Spectrum
To understand the gapped fermion spectrum of thisstripe phase, we first consider a special case in which∆ = δ = 0. Then, H seff ( q, k ) = (cid:15) ( q, k ) 0 0 ∆ ( q, k )0 (cid:15) ( q + , k ) ∆ ( q, k ) 00 ∆ ∗ ( q, k ) (cid:15) ( q, k ) 0∆ ∗ ( q, k ) 0 0 (cid:15) ( q + , k ) . (69)It is obvious that this effective Hamiltonian matrix canbe decomposed into two 2 × (cid:15) (cid:48) ( q, k ) = ( (cid:15) ( q, k ) − (cid:15) ( q + , k )). The four energies E ( q, k ) , . . . E ( q, k ) are E ± , = (cid:15) ( q, k ) + (cid:15) ( q + , k )2 ± (cid:113) | ∆ ( q, k ) | + ( (cid:15) (cid:48) ( q, k )) , (70)where the subscript 1(2) is the index of the degeneratebands. Therefore, the direct gap can be calculated as∆ d ( q, k ) = 2 (cid:113) | ∆ ( q, k ) | + ( (cid:15) (cid:48) ( q, k )) . Since ∆ ( q, k )has the winding number of 1 around the stripe BZ, itmust contain a zero point at a momentum ( q , k ). Be-cause, in general, (cid:15) ( q , k ) (cid:54) = (cid:15) ( q + , k ), this state isfully gapped. The addition of small δ ( q, k ) and ∆ ( q, k )will not close this gap.It is also interesting to study how the double degen-eracy between the two low (high) energy bands is liftedby δ ( q, k ) and ∆ ( q, k ). For this purpose, we first write8down the eigenstates of the Hamiltonian in Eqn. 69: | ψ +1 (cid:105) = (cid:18) cos θ , , , sin θ e − iφ (cid:19) T | ψ − (cid:105) = (cid:18) − sin θ e iφ , , , cos θ (cid:19) T | ψ +2 (cid:105) = (cid:18) , sin θ e iφ , cos θ , (cid:19) T | ψ − (cid:105) = (cid:18) , cos θ , − sin θ e − iφ , (cid:19) T (71)with e iφ = ∆ / | ∆ | and cos θ = (cid:15) (cid:48) / (cid:112) | ∆ | + (cid:15) (cid:48) . Ap-plying the first order perturbation theory with degener-ate states, we obtain H + ij = (cid:104) ψ + i | H | ψ + j (cid:105) (72) H − ij = (cid:104) ψ − i | H | ψ − j (cid:105) (73) H = δ ∆ δ ∗ (cid:48) ∆ ∗ δ ∗ (cid:48)∗ δ . (74)For simplicity, consider the effects of ∆ and ∆ (cid:48) only.We obtain H +12 = ∆ cos θ (cid:48)∗ (cid:18) ∆ | ∆ | (cid:19) sin θ H − = ∆ (cid:48)∗ cos θ (cid:18) ∆ ∗ | ∆ | (cid:19) sin θ H +11 = H +22 and H − = H − because of the C T symmetry. Applying the boundary conditions listed inEqn. 68, we obtain θ ( q + 12 , k ) = π − θ ( q, k ) (77) H +12 ( q + 12 , k ) = (cid:0) H +12 ( q, k ) (cid:1) ∗ e iφ ( q,k ) (78) H +12 ( q, k + 1) = e πiq H +12 ( q, k ) . (79)Although the boundary conditions cannot determine theexact winding number of H +12 around the stripe BZ, theyrestrict the parity of winding number to be even . Asa consequence, the two bands above the CNP can havewinding numbers of 0 , ± , ± , · · · . This conclusion is stillvalid with the inclusion of δ terms .Finally, we should point out that the set of the eigen-states in Eqn. 71 is ill-defined if, at a particular mo-mentum ( q (cid:48) , k (cid:48) ) in the stripe BZ, ∆ ( q (cid:48) , k (cid:48) ) = 0 and (cid:15) (cid:48) ( q (cid:48) , k (cid:48) ) < has the windingnumber of 1 in the stripe BZ, and thus must vanishat a momentum ( q , k ) in the stripe BZ. If this is theonly momentum at which it vanishes, and (cid:15) (cid:48) ( q , k ) < ≤ q < < k < ( q + , k ) = 0 and (cid:15) (cid:48) ( q + , k ) = − (cid:15) (cid:48) ( q , k ) > ( q, k ) acciden-tally vanishes at multiple momenta, the conclusions arestill valid .Our C T stripe state obtained variationally is found tobe close to this limiting case, in that ∆ ( q, k ) dominatesover other matrix elements in most of the BZ. Further-more, when ∆ vanishes at the momentum ( q , k ), thedirect gap ∆ d comes from (cid:15) (cid:48) ( q , k ), and δ and ∆ arenegligible close to ( q , k ).
2. Non-Abelian Topological Charge of Dirac nodes
It is helpful to study another limiting case in which δ ( q, k ) = 0 and (cid:15) (cid:48) ( q, k ) = 0. The effective Hamiltonian H seff ( q, k ) matrix can then be written as H seff ( q, k ) = (cid:32) q, k )∆ ∗ ( q, k ) 0 (cid:33) with ∆( q, k ) = (cid:32) ∆ ( q, k ) ∆ ( q, k )∆ ( q, k ) ∆ (cid:48) ( q, k ) (cid:33) . (80)Since this Hamiltonian anti-commutes with (cid:32) I × − I × (cid:33) , the spectrum is particle-hole sym-metric, and thus the energy must be 0 at the half-filling(CNP). With the boundary conditions given in Eqn. 68,we can readily show that det(∆( q, k )) has the windingnumber of 2 around the stripe BZ, implying thatdet(∆( q, k )) vanishes at two momenta. Therefore,the spectrum must contain two nodes at zero energy.Consider the two zero modes | φ ( q, k ) (cid:105) and | φ ( q, k ) (cid:105) of a single node. The U (1) gauge of these two modescan be chosen such that they are invariant under C T symmetry, i.e. C T | φ (cid:105) = | φ (cid:105) and C T | φ (cid:105) = | φ (cid:105) . Dueto C T symmetry, the effective Hamiltonian H seff withthe matrix element of (cid:0) H seff (cid:1) ij = (cid:104) φ i | ˆ H seff | φ j (cid:105) contains only σ and σ terms. This is still true even withthe addition of (cid:15) ( q, k ) and δ ( q, k ) because the Hamilto-nian should still be C T symmetric. Therefore, introduc-ing a small (cid:15) and δ terms can only shift the position ofnodes without opening a gap. Furthermore, the non-zerowinding number of det(∆) seems to suggest that the chi-rality of two zero energy nodes is the same, naively imply-ing that they cannot be annihilated by meeting together.This contradicts the result of the previous subsection, inwhich the gapped phase is clearly robust.Before resolving this apparent contradiction, it is help-ful to investigate any possible crossings between the up-per two bands. We can focus on the upper two bands,9 (a) (b) FIG. 16. Schematic plot showing how the two orange nodesaround the charge neutrality point (CNP) of a spinless onevalley model with the same topological charge can meet andannihilate each other in the presence of the red (blue) nodesformed by the two lower (upper) bands, respectively. Thetwo nodes still carry the same topological charge if they movealong the dashed path and can’t annihilate, but they havethe opposite charge and can meet, annihilate, and open agap if they move along the solid path. (a) The red and bluenodes coincide because of the particle-hole symmetry. In thiscase, the two orange nodes always carry the same topologicalcharge no matter how they move, and so can’t annihilate (b)By breaking the particle-hole symmetry, the two orange nodescan move along the solid path and carry the opposite topo-logical charge when they meet together. As a consequence,they can annihilate each other and open a gap. because whenever there is a Dirac point between the twoupper bands at the momentum ( q (cid:48) , k (cid:48) ) with the energyof E , the lower two bands also cross at the same mo-mentum with the energy of − E due to the particle-holesymmetry. As a consequence, the matrix ∆ † ∆ = E I × at ( q (cid:48) , k (cid:48) ). Therefore, we obtain two constraints for ma-trix elements at ( q (cid:48) , k (cid:48) ): | ∆ | = | ∆ (cid:48) | and ∆ ∗ ∆ + ∆ ∗ ∆ (cid:48) = 0 . (81)Consider the term f ( q, k ) = ∆ ∗ ∆ + ∆ ∗ ∆ (cid:48) . Applyingthe boundary conditions listed in Eqn. 68, we obtain f ( q + 12 , k ) = f ∗ ( q, k ) , f ( q, k + 1) = − f ( q, k ) . (82)Similar to the previous subsection, this boundary condi-tion does not determine the exact winding number, butrestricts the parity of the winding number to be odd.Therefore, f ( q, k ) contains an odd number of zero points.In addition, notice that f ( q, k ) = 0 when ∆ vanishes atthe odd number of momenta because it has the windingnumber of 1 inside the stripe BZ. As a consequence, thenumber of momentum points at which f ( q, k ) = 0 and∆ ( q, k ) (cid:54) = 0 must be even, as is the number of possiblecrossings between the upper (lower) two bands becausethen the first equation in (81) is automatically satisfied.To resolve the contradiction involving nodes with equalchirality connecting the two middle bands ( E ( q, k ) and E ( q, k )) and the possibility of a gap between the twomiddle bands, we follow Ref. who described the topo-logical properties of nodes in a multiple band system with PT symmetry in 3D as well as C T symmetry in 2D. Itis well-known that the topological charge associated witha node in a two band system with C T symmetry can bedescribed by an integer winding number, or an elementin Z group. However, as Ref. insightfully points out,this description must be modified in a system with morebands. For example, in a three-band system with C T symmetry, the topological charge of nodes should not bethought of as an integer but as a quaternion. With N bands, it is an element in ¯ P N , Salingaros vee group of realClifford algebra C(cid:96) ,N − Ref. . For our C T symmetricperiod 2 stripe, N = 4. The nodes thus anti-commutewith each other if they are from the consecutive bands,and commute otherwise. The two nodes annihilate witheach other if they meet and carry opposite charges, buteven if their charges start out opposite, after braidingone of them with an anti-commuting node, the chargecan change sign and the resulting pair can consequentlyannihilate.Fig. 16 illustrates how the two nodes with the sametopological charge can meet and annihilate with eachother in such a four band system. For notational conve-nience, the bands are (still) labeled by positive integerscounted from the lowest energy to the highest one. Thetwo orange points are the nodes formed by band 2 and 3.We also assume that the system contains the two Diracnodes connecting bands 1 and 2, labeled by red color inFig. 16. If the system is particle-hole symmetric, it alsocontains the two Dirac points at the same momentumbut connecting bands 3 and 4, labeled by blue color. Inthis case, there is no path along which the orange nodecan change its topological charge since any closed loopcontains an even number of nodes formed by neighbor-ing bands (consistent with the winding number 2 of thedeterminant found above). If the particle-hole symmetryis broken, the red and blue nodes move relative to eachother. As shown in Fig. 16(b), there exists a loop en-closing an odd number of nodes from neighboring bands.Therefore, the topological charge of the nodes connect-ing the two middle bands becomes opposite if they meetalong the solid path. As a consequence, the system willbe in the gapped phase without breaking C T symmetry.In order to convincingly show how the gapped state canbe obtained from the gapless state containing two Diracpoints with the same topological charge, we construct afour-band toy model which should accurately capture theregion of the momentum space near the zero(s) of variousterms, but not the periodicity in the full BZ. The matrixelements of the effective Hamiltonian (67) are thus set to∆ = λe iβ ( z − z )( z − z ) (83)∆ = (1 − λ ) z (84)∆ (cid:48) = λe iβ (85) (cid:15) ( q, k ) = − (cid:15) ( q + 12 , k ) = (cid:15) (86) δ ( q, k ) = 0 , (87)where z = p x + ip y is a complex variable, and β , λ , z ,0 FIG. 17. The annihilation of two (orange) nodes in the toymodel given by Eqns. 67, 83 – 87. The convention to label thenodes and bands is the same as the one in Fig. 16. ∆ , givesthe direct gap between band 2 and 3, and thus vanishes atthe orange points. Starting with particle-hole symmetry whenthe two orange nodes have the same topological charge, theircharges become opposite if they move along the solid curveand meet together. If they move along the dashed curves,their topological charges are still the same. The two nodes canannihilate each other only when they have opposite charges. z and (cid:15) are the parameters of this model. It is obvi-ous that the matrix element ∆ in this model has thewinding number of 1, and the determinant of the ∆ ma-trix, det(∆) = ∆ ∆ (cid:48) − ∆ has the winding number of2. It is worth emphasizing that ∆ ∗ ∆ + ∆ ∗ ∆ (cid:48) has theodd winding number of 1, consistent with the analysisabove. This toy model does not satisfy the boundaryconditions listed in Eqn. 68 and should not be extendedto the whole stripe BZ. But, as mentioned, it should ac-curately describe the effective Hamiltonian in a regionof the stripe BZ enclosing the zero points of the matrixelements.In the calculations, we arbitrarily set β = π/ λ = 0 . z = 2 + i , and z = 2 + 0 . i , and vary (cid:15) to study theannihilation of the two nodes at the CNP. The conven-tion for labeling the nodes and bands in Fig. 17 is thesame as the one in Fig. 16. At (cid:15) = 0, when the systemis particle-hole symmetric, the nodes (red points) con-necting the bands 1 and 2, coincide with the nodes (bluepoints) connecting the bands 3 and 4. As a consequence,the two nodes (orange points) connecting the bands 2and 3 with the same topological charge cannot annihi-late each other. As shown in Fig. 17, increasing (cid:15) leadsto the particle-hole symmetry breaking. Thus the nodesconnecting the bands 1 and 2 start to move relative tothe nodes connecting the bands 3 and 4, and thereforeallowing the two orange nodes to annihilate if they meet FIG. 18. The “worldlines” of the Dirac nodes as the param-eter (cid:15) , defined in Eqs.(67,83-87), varies from −∞ to ∞ . Thecurves are obtained using the same C T symmetric modelas in the Fig.(17) with matching colors for the nodes (weomit the label for the p x - p y plane but it should be under-stood). The arrows of the colored curves follow the prescrip-tion for the orientation reversal discussed in a similar contextin Ref. and represent the topological charges of the cor-responding nodes. Starting from (cid:15) = −∞ , there are two bluelines representing the two nodes connecting the bands 3 and4 (ordered in energy), with the same topological charge. Theorange loops corresponds to the nodes connecting bands 2 and3 which, as we prove in the text, must be parallel at (cid:15) = 0.The red lines correspond to the nodes connecting bands 1 and2. Orange anti-commutes with the blue and the red, the lat-ter two commute. Note the orientation reversal each time aline passes under an anti-commuting line; it forces the parallelorientation of blue and the red at large | (cid:15) | . along the dashed path. The toy model shows that thetwo orange nodes annihilate each other when (cid:15) ≈ . providesa general rule to identify the change of the topologicalcharge when the nodes move in the momentum space asour parameter (cid:15) varies. From a fixed vantage point, suchchange happens when a node worldline passes under ananti-commuting worldline. This orientation reversal is il-lustrated in Fig. 18, with the arrow representing the topo-logical charge of the nodes. If the two same-color arrowshave the same orientation at fixed (cid:15) then the topologicalcharges of these two nodes are the same. Otherwise, thetwo charges are opposite. The two nodes with the samecolor can annihilate each other only when they carry theopposite topological charges as they meet.Thus, starting from (cid:15) = −∞ , the system contains two1nodes connecting bands 3 and 4 with the same topologi-cal charge as indicated by the two blue parallel arrows inFig. 18. Later, a pair of red and orange nodes are gener-ated with opposite topological charges. Notice that boththe red curve and the blue curve (which mutually com-mute) pass through the orange loop only once. There-fore the arrows on the blue and the red curves changetheir orientation only once, but the arrows on the orangecurves change twice. As a consequence, the topologicalcharge of the two orange nodes are eventually still op-posite and thus can meet at (cid:15) ∼ ± . VI. DISCUSSION
In this work, we studied the possible ground statesof the TBG near the magic angle at odd integer filling,inspired by the experiments near ν = 3. Using the hy-brid WSs constructed within the spin and valley polar-ized BM model, the projected Coloumb interactions areincluded in the Hamiltonian to study the possible phasesin the strong coupling limit. DMRG identifies differentground states as we vary the ratio of two interlayer hop-ping parameters w and w in the Bistrizer-MacDonaldmodel. When this ratio is small ( ≤ . w /w ≥ .
8, DMRG identifies a state different fromQAH as having lower energy. Surprisingly, this state canbe well approximated as a product state of the hybridWSs, and thus motivates our study of the competitionbetween various phases by minimizing the energy of theDMRG inspired trial wavefunction. Our variational cal-culations discover three nearly degenerate states, QAH, C T nematic state, and C T period-2 stripe state. Thetiny energy difference among them leads to strong compe-tition between these states. Consequently, the manifoldof the low energy states should include all of them even ina spin and valley polarized model, suggesting rich physicsbeyond the U (4) × U (4) manifold .To further obtain the properties of these variousphases, we also calculate their fermion spectrum. TheQAH state is a gapped state, with the minimal gap atΓ point, which is almost certainly further favored by thehBN alignment. On the other hand, the C T nematicstate is a gapless state, with a quadratic node or twovery close Dirac nodes near Γ point, and thus must breakthe C symmetry. This state, as we discussed, has theLandau level degeneracy of 2 , , , · · · . We propose this C T nematic state, with spin and valley degeneracy re- stored, as the candidate for the gapless state observedat the charge neutrality point (CNP), because the fillingfactors of the Landau fan in this state are consistent withthe pattern experimentally observed at the CNP .While the nodes in the C T nematic phase have beenassumed to be generally protected by C T and valley U (1) symmetries, we found that these nodes can be liftedby only breaking the moire translation symmetry, with-out breaking the C T and valley U (1) symmetries, andwithout closing the gap to the remote bands. Our cal-culation shows that a gapped C T period-2 stripe stateis nearly degenerate with the QAH and C T nematicstates, and thus is a candidate state for the ground stateat the filling of ν = 3 without the hBN alignment. Tounderstand how the gap is opened in C T stripe phase,we present an analysis of the topological properties of theDirac nodes in the C T nematic state. During the tran-sition from the C T nematic state to the C T symmetricperiod 2 stripe state, remarkably, the topological chargeassociated with these nodes should not be described bytheir (Abelian) winding number, but by elements of (non-Abelian) Salingaros vee group of real Clifford algebra C(cid:96) , . Since it is a non-Abelian group, the topologicalcharge of these nodes depends on how they are braidedwith other nodes away from the CNP. Therefore, a gap atCNP can be opened even without breaking the C T andvalley U (1) symmetry. We expect that this mechanismis general and applies when spin and valley degrees offreedom are fully restored, in which case a gap at odd in-teger filling may not necessitate translational symmetrybreaking.Finally, the mechanism discussed makes it apparentthat the gap opening in a C T symmetric, but moiretranslation symmetry broken, state relies on the non-Abelian topological charges of the Dirac nodes, whichis effective only if the particle-hole symmetry is broken.Otherwise, the node lines providing the non-trivial braid-ing are glued together and the equal chirality nodes atthe neutrality cannot annihilate. Therefore, if the par-ticle hole symmetry is a good symmetry, the C T sym-metric state must remain gapless even when translationsymmetry is broken. This means that it is in principlepossible to be in the strong coupling limit, have weakparticle-hole symmetry breaking, and end up in a statewhich has a gap parametrically smaller than U , the scaleset by the Coulomb repulsion. ACKNOWLEDGMENTS
We thank B. Andrei Bernevig, Leon Balents, KasraHejazi, Nicolas Regnault, Tomo Soejima, and MichaelZaletel for discussions. We are especially grateful toHitesh Changlani and P. Myles Eugenio for help with theDMRG calculations. J. K. is supported by Priority Aca-demic Program Development (PAPD) of Jiangsu HigherEducation Institutions and was partially supported bythe National High Magnetic Field Laboratory through2NSF Grant No. DMR-1157490 and the State of Florida.O. V. was supported by NSF DMR-1916958. Part ofthis work was performed while the authors visited theAspen Center for Physics which is supported by the Na-tional Science Foundation grant PHY-1607611, and the Kavli Institute of Theoretical Physics which is supportedin part by the National Science Foundation under GrantNo. NSF PHY-1748958. J. K. also thanks the Kavli In-stitute for Theoretical Sciences for hospitality during thecompletion of this work. ∗ [email protected] † [email protected] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken,J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T.Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, “Correlated insulator behaviour at half-fillingin magic-angle graphene superlattices,” Nature , 43(2018). Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,E. Kaxiras, and P. Jarillo-Herrero, “Unconventional super-conductivity in magic-angle graphene superlattices,” Na-ture , 80 (2018). M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watan-abe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean,“Tuning superconductivity in twisted bilayer graphene,”,Science , 1059 (2019). A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney,K. Watanabe, T. Taniguchi, M. A. Kastner, andD. Goldhaber-Gordon, “Emergent ferromagnetism nearthree-quarters filling in twisted bilayer graphene,” Science , 605 (2019). M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu,K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young,“Intrinsic quantized anomalous hall effect in a moire het-erostructure,” Science (2019), 10.1126/science.aay5533. A. Kerelsky, L. J McGilly, D. M. Kennes, L. Xian, M.Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone,C. Dean, et al., “Maximized electron interactions at themagic angle in twisted bilayer graphene”, Nature , 95(2019). S. L. Tomarken, Y. Cao, A. Demir, K. Watanabe, T.Taniguchi, P. Jarillo-Herrero, and R. C. Ashoori, “Elec-tronic compressibility of magic-angle graphene superlat-tices,” Phys. Rev. Lett. , 046601 (2019). Y. Xie, B. Lian, B. Jack, X. Liu, C.-L. Chiu, K. Watanabe,T. Taniguchi, B. A. Bernevig, and A. Yazdani, “Spectro-scopic signatures of many body correlations in magic-angletwisted bilayer graphene,” Nature , 101 (2019). Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J.Mao, and E. Y. Andrei, ”Charge order and broken rota-tional symmetry in magic-angle twisted bilayer graphene,”Nature , 91 (2019). X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das,C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bach-told, A. H. MacDonald, and D. K. Efetov, “Superconduc-tors, orbital magnets and correlated states in magic-anglebilayer graphene,” Nature , 653 (2019). D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S.Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A.Yazdani, “Cascade of transitions between the correlatedelectronic states of magic-angle twisted bilayer graphene,”arXiv:1912.06145. U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R.Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Op- pen, A. Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani,“Cascade of Phase Transitions and Dirac Revivals in MagicAngle Graphene,” arXiv:1912.06150 Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F.Young, “Decoupling superconductivity and correlated in-sulators in twisted bilayer graphene,” arXiv:1911.13302 P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe,T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levi-tov, D. K. Efetov, “The interplay of insulating and su-perconducting orders in magic-angle graphene bilayers,”,arXiv:1911.09198. Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora,R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F.von Oppen, K. Watanabe, T. Taniguchi, S. Nadj-Perge,“Electronic correlations in twisted bilayer graphene nearthe magic angle”, Nat. Phys. , 1174 (2019) C. Shen, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian,Y. Chu, K. Watanabe, T. Taniguchi, R. Yang, Z.Y. Meng,D. Shi, G. Zhang, “Observation of superconductivity withtc onset at 12K in electrically tunable twisted double bi-layer graphene,” arXiv:1903.06952. X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, K. Watanabe, T.Taniguchi, A. Vishwanath, and P. Kim, “Spin-polarizedcorrelated insulator and superconductor in twisted doublebilayer graphene,” arXiv:1910.04654. Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M.Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero,“Electric field tunable correlated states and magneticphase transitions in twisted bilayer-bilayer graphene,”arXiv:1903.08596. G. Chen, A. L. Sharpe, E. J. Fox, Y.-H. Zhang, S. Wang,L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi,Z. Shi, T. Senthil, D. Goldhaber-Gordon, Y. Zhang, F.Wang, “Tunable correlated chern insulator and ferromag-netism in trilayer graphene/boron nitride moire superlat-tice,” arXiv:1905.06535. G. Chen, L. Jiang, S. Wu, B. Lyu, H. Li, B. L. Chittari, K.Watanabe, T. Taniguchi, Z. Shi, J. Jung, Y. Zhang, andF. Wang, “Evidence of a gate-tunable Mott insulator in atrilayer graphene moire superlattice,” Nature Physics ,237 (2019). G. Chen, A. L. Sharpe, P. Gallagher, I. T. Rosen, E. J.Fox, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi,J. Jung, Z. Shi, D. Goldhaber-Gordon, Y. Zhang, andF. Wang, “Signatures of tunable superconductivity in atrilayer graphene moire superlattice,” Nature , 215(2019). R. Bistritzer and A. H. MacDonald, “Moire bands intwisted double-layer graphene,” Proc. Natl. Acad. Sci.U.S.A. , 12233 (2011). C. Xu and L. Balents, “Topological superconductivityin twisted multilayer graphene,” Phys. Rev. Lett. ,087001 (2018). J. Kang and O. Vafek, “Symmetry, maximally localized Wannier states, and a low-energy model for twisted bilayergraphene narrow bands,” Phys. Rev. X , 031088 (2018). M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi,K. Kuroki, and L. Fu, “Maximally localized wannier or-bitals and the extended hubbard model for twisted bilayergraphene,” Phys. Rev. X , 031087 (2018). H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, “Origin ofmott insulating behavior and superconductivity in twistedbilayer graphene,” Phys. Rev. X , 031089 (2018). C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, “Chi-ral spin density wave and d + id superconductivity in themagic-angle-twisted bilayer graphene”, Phys. Rev. Lett. , 217001 (2018). M. Ochi, M. Koshino, K. Kuroki, “Possible correlated insu-lating states in magic-angle twisted bilayer graphene understrongly competing interactions”, Phys. Rev. B , 081102(2018). J. F. Dodaro, S. A. Kivelson, Y. Schattner, X. Q. Sun, andC. Wang, “Phases of a phenomenological model of twistedbilayer graphene,” Phys. Rev. B , 075154 (2018). H. Guo, X. Zhu, S. Feng, and R. T. Scalettar, “Pair-ing symmetry of interacting fermions on a twisted bilayergraphene superlattice”, Phys. Rev. B , 235453 (2018). Hiroki Isobe, Noah F. Q. Yuan, and Liang Fu, “Uncon-ventional superconductivity and density waves in twistedbilayer graphene,” Phys. Rev. X , 041041 (2018). Louk Rademaker and Paula Mellado, “Charge-transfer in-sulation in twisted bilayer graphene,” Phys. Rev. B ,235158 (2018). J. W. F. Venderbos and R. M. Fernandes, “Correlationsand electronic order in a two-orbital honeycomb latticemodel for twisted bilayer graphene,” Phys. Rev. B ,245103 (2018). F. Guinea and N. R Walet, “Electrostatic effects, banddistortions, and superconductivity in twisted graphene bi-layers,” Proc. Natl. Acad. Sci. U.S.A. , 13174 (2018). Q.-K. Tang, L. Yang, D. Wang, F.-C. Zhang, and Q.-H. Wang, “Spin-triplet f-wave pairing in twisted bilayergraphene near -filling,” Phys. Rev. B , 094521 (2019). J. Gonzalez and T. Stauber, “Kohn-luttinger superconduc-tivity in twisted bilayer graphene,” Phys. Rev. Lett. ,026801 (2019). J. Kang and O. Vafek, “Strong coupling phases of partiallyfilled twisted bilayer graphene narrow bands,” Phys. Rev.Lett. , 246401 (2019). K. Seo, V. N. Kotov, and B. Uchoa, “Ferromagnetic Mottstate in twisted graphene bilayers at the magic angle,”Phys. Rev. Lett. , 246402 (2019). Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T.Senthil, “Nearly flat chern bands in moire superlattices,”Phys. Rev. B , 075127 (2019). J. Y. Lee, E. Khalaf, S. Liu, X. Liu, Z. Hao, P. Kim, andA. Vishwanath, “Theory of correlated insulating behav-ior and spin-triplet superconductivity in twisted doublebilayer graphene,” Nat. Commun. , 1 (2019). X.-C. Wu, A. Keselman, C.-M. Jian, K. A. Pawlak, and C.Xu, “Ferromagnetism and spin-valley liquid states in moirecorrelated insulators,” Phys. Rev. B , 024421 (2019). M. Xie, and A. H. MacDonald, “On the nature of thecorrelated insulator states in twisted bilayer graphene,”arXiv:1812.04213. N. Bultinck, S. Chatterjee, and M. P. Zaletel, “Anoma-lous hall ferromagnetism in twisted bilayer graphene,” arXiv:1901.08110. S. Liu, E. Khalaf, J. Y. Lee, and A. Vishwanath, “Nematictopological semimetal and insulator in magic angle bilayergraphene at charge neutrality,” arXiv:1905.07409. Y. Alavirad and J. D. Sau, “Ferromagnetism and its sta-bility from the one-magnon spectrum in twisted bilayergraphene,” arXiv:1907.13633. T. Huang, L. Zhang, and T. Ma, “Antiferromagneticallyordered Mott insulator and d + id superconductivity intwisted bilayer graphene: A quantum Monte Carlo study,”Sci. Bull. , 310 (2019). S. Chatterjee, N. Bultinck, and M. P. Zaletel, “Symme-try breaking and skyrmionic transport in twisted bilayergraphene,” arXiv:1908.00986. D. V. Chichinadze, L. Classen, and A. V. Chubukov,“Nematic superconductivity in twisted bilayer graphene,”arXiv:1910.07379. N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vish-wanath, and M. P. Zaletel, “Ground State and HiddenSymmetry of Magic Angle Graphene at Even Integer Fill-ing,” arXiv:1911.02045. J. Liu and X. Dai, “Correlated insulating states and thequantum anomalous Hall phenomena at all integer fillingsin twisted bilayer graphene”, arXiv: 1911.03760. R. M Fernandes and J. W. F. Venderbos, ”Nematicity witha twist: rotational symmetry breaking in a moir´e superlat-tice,” arXiv:1911.11367. Y. Zhang, K. Jiang, Z. Wang, and F. C. Zhang, “Cor-related insulating phases of twisted bilayer graphene atcommensurate filling fractions: a Hartree-Fock study,”arXiv:2001.02476. F. Wu and S. Das Sarma, “Collective Excitations of Quan-tum Anomalous Hall Ferromagnets in Twisted BilayerGraphene,” Phys. Rev. Lett. , 046403 (2020). C. Repellin, Z. Dong, Y.-H. Zhang, T. Senthil, “Fer-romagnetism in narrow bands of moire superlattices,”arXiv:1907.11723. L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, “Bandstructure of twisted bilayer graphene: Emergent symme-tries, commensurate approximants, and Wannier obstruc-tions,” Phys. Rev. B , 085435 (2018). G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, “Ori-gin of magic angles in twisted bilayer graphene,” Phys.Rev. Lett. , 106405 (2019). J. Ahn, S. Park, and B.-J. Yang, “Failure of nielsen-ninomiya theorem and fragile topology in two dimensionalsystems with space-time inversion symmetry: Applicationto twisted bilayer graphene at magic angle,” Phys. Rev. X , 021013 (2019). Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. AndreiBernevig, “All magic angles in twisted bilayer graphene aretopological,” Phys. Rev. Lett. , 036401 (2019). K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Ba-lents, “Multiple topological transitions in twisted bilayergraphene near the first magic angle,” Phys. Rev. B ,035111 (2019). J. Liu, J. Liu, and X. Dai, “The pseudo-Landau-level rep-resentation of twisted bilayer graphene: band topology andthe implications on the correlated insulating phase,” Phys.Rev. B , 155415 (2019). T. Neupert, L. Santos, S. Ryu, C. Chamon, and C. Mudry,“Topological Hubbard Model and Its High-TemperatureQuantum Hall Effect”, Phys. Rev. Lett. , 046806(2012). V. M. Pereira, A. H. Castro Neto, and N. M. R. Peres,“Tight-binding approach to uniaxial strain in graphene,”Phys. Rev. B , 045401 (2009). Q. Wu, A. A. Soluyanov, and T. Bzduˇsek, “Non-Abelianband topology in noninteracting metals,” Science ,1273 (2019). E. McCann and V. I. Falko, “Landau-Level Degeneracyand Quantum Hall Effect in a Graphite Bilayer,” Phys.Rev. Lett. 96, 086805 (2006). Y. H. Zhang, H. C. Po, and T. Senthil, “Landau level de-generacy in twisted bilayer graphene: Role of symmetrybreaking” Phys. Rev. B , 125104 (2019). Y. Cao, J. Y. Luo, V. Fatemi, S. Fang, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, andP. Jarillo-Herrero “Superlattice-Induced Insulating Statesand Valley-Protected Orbits in Twisted Bilayer Graphene”Phys. Rev. Lett. , 116804 (2016). L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young,“Superconductivity and strong correlations in moir´e flatbands”, Nat. Phys. , 725 (2020). N. N. T. Nam and M. Koshino, “Lattice relaxation and en-ergy band modulation in twisted bilayer graphene,” Phys.Rev. B , 75311 (2017). P. Lucignano, D. Alf´e, V. Cataudella, D. Ninno, and G.Cantele, “Crucial role of atomic corrugation on the flatbands and energy gaps of twisted bilayer graphene at the magic angle θ ∼ . ◦ ”, Phys. Rev. B , 195419 (2019). N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D.Vanderbilt, “Maximally localized Wannier functions: The-ory and applications,” Rev. Mod. Phys. X. Wang and O. Vafek, “A diagnosis of explicit symmetrybreaking in the tight-binding constructions for symmetry-protected topological systems,” arXiv:2002.02057. R. Yu, X. L. Qi, A. Bernevig, Z. Fang, and X. Dai, “Equiv-alent expression of Z topological invariant for band insula-tors using the non-Abelian Berry connection,” Phys. Rev.B , 075119 (2011). R. Resta, “Quantum-Mechanical Position Operator in Ex-tended Systems,” Phys. Rev. Lett. , 1800 (1998). See Appendix for more detailed derivation of the hybridWannier states, gate-screened Coulomb interaction, ener-getics of various states at w /w = 0 .
3, parameterizationof the C T stripe state, and the calculation of windingnumbers of the matrix elements of the effective Hamilto-nian. O. Vafek, and A. Vishwanath, “Dirac Fermions in Solids:From High-Tc Cuprates and Graphene to Topological In-sulators and Weyl Semimetals,” Ann. Rev. of Cond. Mat.Phys., , 83 (2014). ITensor Library (version 2.0.11) http://itensor.org Q. Wu, A. A. Soluyanov, and T. Bzduˇsek, SupplementaryMaterials for “Non-Abelian band topology in noninteract-ing metals,” Science , 1273 (2019).
Appendix for “Non-Abelian Dirac node braiding and near-degeneracy of correlatedphases at odd integer filling in magic angle twisted bilayer graphene”
Appendix A: Hybrid Wannier states
In this section, we illustrate our approach in detail to construct the hybrid Wannier states(WS). They are theeigenstates of the position operator generalized to the periodic boundary conditions, and projected onto the narrowbands, ˆ O = ˆ P e − i N g · r ˆ P (A1)The obtained states are maximally (exponentially) localized in one direction and extended Bloch states in the otherdirection . The projection operator, ˆ P , onto the narrow band composite can be written in terms of the (energyeigenstate) Bloch states | k , m (cid:105) as (cid:80) k ,m | k , m (cid:105)(cid:104) k , m | , where k resides in the first Brillouin zone (BZ), m is the bandindex and ˆ g is a primitive vector of the reciprocal lattice. The Bloch state overlaps can be expressed in terms ofoverlaps of the periodic part of the Bloch function, u k ,m ( r ), as (cid:104) k , m | e − i N g · r | k (cid:48) , m (cid:48) (cid:105) = (cid:88) G δ k (cid:48) , k + N ˆ g + G (cid:88) r ∈ uc e i G · r u ∗ k ,m ( r ) u k (cid:48) ,m (cid:48) ( r ) , (A2)where G is a reciprocal lattice vector. The operator ˆ O thus mixes only states in the first BZ which lie along the lineparallel to ˆ g . If we focus on one such line along g , and choose N = N to be the number of unit cells along thedirection of L , the operator which we wish to diagonalize isˆ O η = (cid:88) mm (cid:48) N − (cid:88) j =0 | η g + jN g , m (cid:105) Λ mm (cid:48) ( η, j ) (cid:104) η g + j + 1 N g , m (cid:48) | + | η g + N − N g , m (cid:105) Λ mm (cid:48) ( η, N − (cid:104) η g , m (cid:48) | (cid:19) , (A3)5where 0 ≤ η < mm (cid:48) ( η, j ) = (cid:88) r ∈ uc u ∗ η g + jN g ,m ( r ) u η g + j +1 N g ,m (cid:48) ( r ); j (cid:54) = N − mm (cid:48) ( η, N −
1) = (cid:88) r ∈ uc e − i g · r u ∗ η g + N − N g ,m ( r ) u η g ,m (cid:48) ( r ) . (A5)We seek a solution of the form N − (cid:88) j =0 (cid:88) m α η,j,m | η g + jN g , m (cid:105) , so that (A6)ˆ O η N − (cid:88) j =0 (cid:88) m α η,j,m | η g + jN g , m (cid:105) = (cid:15) η N − (cid:88) j =0 (cid:88) m α η,j,m | η g + jN g , m (cid:105) . (A7)Clearly, this leads to the equations that (cid:88) m (cid:48) Λ mm (cid:48) ( η, j ) α η,j +1 ,m (cid:48) = (cid:15) η α η,j,m with j = 0 , · · · , N − (cid:88) m (cid:48) Λ mm (cid:48) ( η, N − α η, ,m (cid:48) = (cid:15) η α η,N − ,m . (A8)This gives us an eigenvalue problem (cid:88) m (cid:48) W mm (cid:48) ( η ) α η, ,m (cid:48) = ( (cid:15) η ) N α η, ,m , (A9)where the matrix W ( η ) = Λ( η, η, · · · Λ( η, N − N → ∞ , W ( η ) is a unitary matrix. Inpractice, for finite N , we perform an singular value decomposition (SVD) of each Λ( η, j ) = U ( η, j )Σ( η, j ) V † ( η, j )where Σ( η, j ) is diagonal (and very close to 1) and U and V are unitary, and replace Σ with a unit matrix. Once wehave the eigenvectors α η, and the eigenvalues (cid:15) N η , we can construct the remaining α ’s as α η,j = ( (cid:15) η ) j Λ − ( η, j − . . . Λ − ( η, − ( η, α η, = ( (cid:15) η ) j (Λ( η, η, . . . Λ( η, j − − α η, . (A10)Note that if (cid:15) N η is an eigenvalue, then replacement (cid:15) η → e − πi nN (cid:15) η also gives an eigenvalue for an arbitrary integer n . Thus, we set the phase of (cid:15) η to be in the range of − π/N ≤ arg( (cid:15) η ) < π/N . For the two bands of interest, thereare two eigenvectors α ± . Our hybrid Wannier states are therefore | w ± ( n, η g ) (cid:105) = ce iχ ± η (cid:18) α ± ,m | η g , m (cid:105) + . . . + e − πi jnN (cid:15) jη (Λ( η, η, . . . Λ( η, j − − mm (cid:48) α ± ,m (cid:48) | η g + jN g , m (cid:105) + . . . + e − πi ( N − nN (cid:15) N − η (Λ( η, η, . . . Λ( η, N − − mm (cid:48) α ± ,m (cid:48) | η g + N − N g , m (cid:105) (cid:19) , (A11)with repeated m, m (cid:48) indices summed. c is a positive number added for normalization. Note that under translationby L , ˆ T L | w ± ( n, η g ) (cid:105) = | w ± ( n + 1 , η g ) (cid:105) , which holds if χ ± η is n -independent. Under the translation by L ,ˆ T L | w ± ( n, η g ) (cid:105) = e − πiη | w ± ( n, η g ) (cid:105) . C T Symmetry of the hybrid Wannier states
We first fix the phase of our Bloch states by choosing them to be eigenstates of C T with a unit eigenvalue. Withthis convention, it is clear that (cid:104) r | k , m (cid:105) = (cid:104)− r | k , m (cid:105) ∗ = ⇒ u k ,m ( r ) = u ∗ k ,m ( − r ) . (A12)This guarantees that Λ’s, and therefore W , are real, and because for each valley and spin we have only two bands, itis a 2 × W = e iθσ where σ = (cid:32) − ii (cid:33) . The eigenstates of σ can bechosen to be √ (1 , ± i ) T with W eigenvalues e ± iθ . We choose α +0 in such a way that its eigenvalue winds by 2 π as η α − in such a way that its eigenvaluewinds by − π as η goes from 0 to 1, corresponding to the Chern − η we may thus haveeither √ (1 , i ) T or √ (1 , − i ) T as α +0 , for either choice, complex conjugation interchanges α +0 and α − , and their O η eigenvalues (cid:15) ± η . The action of C T on the Chern ± | w ± ( n, η g ) (cid:105) is thereforeˆ C T | w ± ( n, η g ) (cid:105) = ˆ C T ce iχ ± η (cid:18) α ± ,m | η g , m (cid:105) + . . . + e − πi jnN (cid:0) (cid:15) ± η (cid:1) j (Λ( η, η, . . . Λ( η, j − − mm (cid:48) α ± ,m (cid:48) | η g + jN g , m (cid:105) + . . . + e − πi ( N − nN (cid:0) (cid:15) ± η (cid:1) N − (Λ( η, η, . . . Λ( η, N − − mm (cid:48) α ± ,m (cid:48) | η g + N − N g , m (cid:105) (cid:19) = ce − iχ ± η (cid:18) α ∓ ,m | η g , m (cid:105) + . . . + e πi jnN (cid:0) (cid:15) ∓ η (cid:1) j (Λ( η, η, . . . Λ( η, j − − mm (cid:48) α ∓ ,m (cid:48) | η g + jN g , m (cid:105) + . . . + e πi ( N − nN (cid:0) (cid:15) ∓ η (cid:1) N − (Λ( η, η, . . . Λ( η, N − − mm (cid:48) α ∓ ,m (cid:48) | η g + N − N g , m (cid:105) (cid:19) . (A13)Therefore, our hybrid WSs satisfy the constraintˆ C T | w ± ( n, η g ) (cid:105) = | w ∓ ( − n, η g ) (cid:105) . (A14)as long as χ + η = − χ − η . The C T symmetry is therefore implemented “on-site” in the hybrid Wannier state basis,where n and k ≡ η g are the generalized “sites”.
2. Continuity of the hybrid WSs
Before proceeding to the discussion of other symmetries, we should address the continuity of he hybrid WSs. Werequire that | w ± ( n, η g ) (cid:105) should be continuous in terms of η and | w ± ( n, (1 + η ) g ) (cid:105) = | w ± ( n ± , η g ) (cid:105) . For thispurpose, we consider how to fix the phase factor χ + η = − χ − η .First, notice that the phase of the Bloch states are almost fixed by C T except their signs: C T | k , m (cid:105) = | k , m (cid:105) ⇐⇒ C T ( −| k , m (cid:105) ) = −| k , m (cid:105) , To further remove this sign freedom, we apply the constraints that (cid:60) (cid:32) (cid:88) r ∈ uc u ∗ η g ,m ( r ) u ( η + δη ) g ,m ( r ) (cid:33) > . (A15)With this convention set for Bloch states, it is obvious that the hybrid WSs defined in Eqn. A11 are continuous interms of η as long as χ ± η are smooth functions of η .It is also interesting to investigate | w ± ( n, (1 + η ) g (cid:105) . In the most general form, the Bloch states | (1 + η ) g + jN g , m (cid:105) = e iθ η,j,m | η g , jN g , m (cid:105) . (A16)Applying C T on both sides, we found that e iθ η,j,m = e − iθ η,j,m , ie. θ η,j,m = 0 or π . Since the phase of Bloch states with j = 0 are fixed by Eqn. A15, θ η,j =0 ,m should be a smooth function of η . As a consequence, θ η,j =0 ,m is independentof η . Numerically, we found that θ η,j =0 ,m always vanishes for both bands.Now, consider | w ± ( n, η g ) (cid:105) defined in Eqn. A11. As the hybrid WSs carry the Chern indices ± , the correspondingeigenvalues of the projected position operator O have the properties of (cid:15) ± η = (cid:15) ± η e ∓ i π/N . Furthermore, since | (1 + η ) g + jN g , m (cid:105) = | η g + jN g , m (cid:105) = ⇒ u (1+ η ) g + jN g ,m ( r ) = e − i g · r u η g + jN g ,m ( r ) and Λ mm (cid:48) (1 + η, j ) = Λ mm (cid:48) ( η, j ) | w ± ( n, (1 + η ) g ) (cid:105) = e iχ ± η (cid:0) α ± ,m | η g , m (cid:105) + . . . + e − πi j ( n ± N (cid:0) (cid:15) ± η (cid:1) j (Λ( η, η, . . . Λ( η, j − − mm (cid:48) α ± ,m (cid:48) | η g + jN g , m (cid:105) + . . . + e − πi ( N − n ± N (cid:0) (cid:15) ± η (cid:1) N − (Λ( η, η, . . . Λ( η, N − − mm (cid:48) α ± ,m (cid:48) | η g + N − N g , m (cid:105) (cid:19) = ⇒ | w ± ( n, (1 + η ) g ) (cid:105) = e i ( χ ± η − χ ± η ) | w ± ( n ± , η g ) (cid:105) (A17)7Now, it is easy to see that we can set χ ± η = 0 so that the hybrid WSs | w ± ( n, η g ) (cid:105) is a smooth function of η , andsatisfies | w ± ( n, (1 + η ) g ) (cid:105) = | w ± ( n ± , η g ) (cid:105) and C T | w ± ( n, η g ) (cid:105) = | w ∓ ( − n, η g ) (cid:105) . Assuming the Bloch states are properly normalized, (cid:104) g , m | g (cid:48) , m (cid:48) (cid:105) = δ gg (cid:48) δ mm (cid:48) , we obtained (cid:104) w α ( n, η g ) | w α (cid:48) ( n (cid:48) η (cid:48) g ) (cid:105) ∝ δ αα (cid:48) δ nn (cid:48) δ ηη (cid:48) and (cid:104) w α ( n, η g ) | w α ( n, η g ) (cid:105) = N c (A18)We choose that c = ( N ) − / for normalization. C (cid:48)(cid:48) Symmetry of the Hybrid WSs
We also wish to find how the C (cid:48)(cid:48) transformation acts on the hybrid Wannier states. To this end we note that forthe Bloch states, ˆ C (cid:48)(cid:48) | η g + jN g , m (cid:105) ∝ | (1 − η ) g + ( jN − η ) g , m (cid:105) . (A19)Let η = hN , where h is an integer. Then, because Bloch states at wavevectors related by a reciprocal lattice vectorare identical, and C (cid:48)(cid:48) ( x L + y L ) C (cid:48)(cid:48) = x ( L − L ) + y L = ⇒ C (cid:48)(cid:48) e − i N g · r C (cid:48)(cid:48) = e − i N g · r (A20)= ⇒ ˆ C (cid:48)(cid:48) ˆ O η ˆ C (cid:48)(cid:48) = ˆ O − η (cid:0) C (cid:48)(cid:48) (cid:1) =1 ====== ⇒ ˆ O − η ˆ C (cid:48)(cid:48) | w ± ( n, η g ) (cid:105) = (cid:15) ± η e − i πn/N ˆ C (cid:48)(cid:48) | w ± ( n, η g ) (cid:105) . (A21)Note that the state at 1 − η , obtained from the state at η has the same eigenvalue (cid:15) η . But, because the phase of theeigenvalues of W winds by ± π as η changes from 0 to 1, we must have (cid:0) (cid:15) ± η (cid:1) N = (cid:0) (cid:15) ∓ − η (cid:1) N . Therefore, up to a phasefactor, the ± C (cid:48)(cid:48) : e iφ ± η ( n ) | w ∓ ( n, (1 − η ) g ) (cid:105) = ˆ C (cid:48)(cid:48) | w ± ( n, η g ) (cid:105) . Consider the case when n = 0, we found C (cid:48)(cid:48) | w + (0 , η g ) (cid:105) = e iφ + η (0) | w − (0 , (1 − η ) g ) (cid:105) = ⇒ C (cid:48)(cid:48) | w − (0 , η g ) (cid:105) = e − iφ + η (0) | w + (0 , (1 − η ) g ) (cid:105) (A22)The second formula is derived by applying C T to both sides of the first formula. Thus, φ + η (0) = − φ − η (0). By applying C (cid:48)(cid:48) to both sides of the first formula, we obtain | w + (0 , η g ) (cid:105) = e iφ + η (0) C (cid:48)(cid:48) | w − (0 , (1 − η ) g ) (cid:105) = ⇒ φ − − η (0) = − φ + η (0) = ⇒ φ ± η (0) = − φ ∓ η (0) = φ ± − η (0) . Since the hybrid WSs are smooth with respect to η , φ ± η (0) are also smooth functions of η . We redefining the phasefactor χ ± η → χ ± η − φ ± η (0) /
2. Therefore, C (cid:48)(cid:48) | w ± (0 , η g ) (cid:105) = | w ∓ (0 , (1 − η ) g ) (cid:105) , C T | w ± ( n, η g ) (cid:105) = | w ∓ ( − n, η g ) (cid:105)| w ± ( n, (1 + η ) g ) (cid:105) = | w ± ( n ± , η g ) (cid:105) . (A23)Because ˆ C (cid:48)(cid:48) ˆ T n L = ˆ T n ( L − L ) ˆ C (cid:48)(cid:48) , we haveˆ C (cid:48)(cid:48) | w ± ( n, η g ) (cid:105) = ˆ C (cid:48)(cid:48) ˆ T n L | w ± (0 , η g ) (cid:105) = ˆ T − n L ˆ T n L ˆ C (cid:48)(cid:48) | w ± (0 , η g ) (cid:105) = ˆ T − n L ˆ T n L | w ∓ (0 , (1 − η ) g ) (cid:105) = ˆ T − n L | w ∓ ( n, (1 − η ) g ) (cid:105) = e − πinη | w ∓ ( n, (1 − η ) g ) (cid:105) . (A24)This means that C (cid:48)(cid:48) is also implemented “on-site”.8 - - - - - - - (a) - - - (b) FIG. S1. The kinetic energy vector fields ( Eqn. A28 ) obtained from (locally) smooth gauge using hybrid Wannier states (left) n ( q, k ) and (right) n ( q, k ). The bold green lines correspond to the points where n , ( q, k ) = 0. The intersection of the greenlines results in the Dirac nodes. As shown in Fig. 4, the chirality of the two Dirac points is the same, seemingly contradictingthe fermion doubling theorem. Note, however, that one of the assumptions of the theorem does not hold, namely, the vectorfields n , ( q, k ) need to be smooth and periodic. As this figure shows, n ( q, k ) is periodic but is not globally smooth, instead,it has a branch-cut discontinuity at k = 0 and k = 1.
4. Kinetic energy H kin = (cid:88) n,n (cid:48) (cid:88) k (cid:88) α,α (cid:48) = ± (cid:104) w α ( n, k g ) | H | w α (cid:48) ( n (cid:48) , k g ) (cid:105) d † α,n,k d α (cid:48) ,n (cid:48) ,k = t αα (cid:48) ( n − n (cid:48) , k ) d † α,n,k d α (cid:48) ,n (cid:48) ,k . (A25)It is obvious that t αα (cid:48) ( n, k ) = t ∗ α (cid:48) α ( − n, k ). Due to C T , we have t ++ ( n − n (cid:48) , k ) = (cid:104) w + ( n, k g ) | H | w + ( n (cid:48) , k g ) (cid:105) = (cid:104) w − ( − n, k g ) | H | w − ( − n (cid:48) , k g ) (cid:105) ∗ = (cid:104) w − ( − n (cid:48) , k g ) | H | w − ( − n, k g ) (cid:105) = t −− ( n − n (cid:48) , k ) (A26)Moreover, due to C (cid:48)(cid:48) t + − ( n − n (cid:48) , k ) = (cid:104) w + ( n, k g ) | ˆ C (cid:48)(cid:48) H ˆ C (cid:48)(cid:48) | w − ( n (cid:48) , k g ) (cid:105) = e πik ( n − n (cid:48) ) (cid:104) w − ( n, (1 − k ) g ) | H | w + ( n (cid:48) , (1 − k ) g ) (cid:105) = e πik ( n − n (cid:48) ) t − + ( n − n (cid:48) , − k ) . (A27)Fourier transforming gives the hybridized band operators, d ± ,n,η = N − / (cid:80) N − q =0 e πinq b ± ,q,k , in terms of which thekinetic energy has the form H kin = (cid:88) k,q (cid:32) b + ,q,k b − ,q,k (cid:33) † (cid:32) n ( q, k ) + n ( q, k ) n ( q, k ) − in ( q, k ) n ( q, k ) + in ( q, k ) n ( q, k ) − n ( q, k ) (cid:33) (cid:32) b + ,q,k b − ,q,k (cid:33) = (cid:88) αα (cid:48) = ± (cid:88) q,k t αα (cid:48) ( q, k ) b † α,q,k b α (cid:48) ,q,k , (A28)where t αα (cid:48) ( q, k ) = (cid:88) δn t αα (cid:48) ( δn, k ) e πiqδn (A29)With Eqn. A26, it is obvious that t ++ ( q, k ) = t −− ( q, k ), ie. n ( q, k ) = 0. Furthermore, with Eqn. A27 by C (cid:48)(cid:48) symmetry( ), t + − ( q, k ) = (cid:88) δn t + − ( δn, k ) e πiqδn = (cid:88) δn t − + ( δn, − k ) e πikδn e πiqδn = t − + ( q + k, − k )= ⇒ n ( q, k ) = n ( q + k, − k ) , n ( q, k ) = − n ( q + k, − k ) . (A30)This in turn implies that, if the two Dirac nodes are present, then the winding numbers of at the two Dirac nodes are thesame. The hybrid Wannier states thus provide means to construct a locally smooth gauge with on-site representationof the C T and C (cid:48)(cid:48) symmetries. However, the gauge is not globally smooth, in that there are branch-cuts at k = 0and k = 1. Otherwise, there would be an additional pair of Dirac nodes at the location of the branch-cuts, both nodeswith the same chirality, cancelling the overall chirality in the Brillouin zone.9 Appendix B: Gate-Screened Coulomb Interaction
In this section, we calculate the metallic gate-screened Coulomb interaction in the bilayer system, with two graphenelayers separated by the distance d ⊥ , located in the middle of two gates. Assuming the distance between two gates is ξ , due to the gate screening effect, the screened Coulomb interaction is V intra ( r ) = e π(cid:15) hBN ∞ (cid:88) n = −∞ ( − ) n (cid:112) r + ( nξ + (( − ) n − d ⊥ / = e π(cid:15) hBN ∞ (cid:88) n = −∞ (cid:32) (cid:112) r + (2 nξ ) − (cid:112) r + ((2 n + 1) ξ − d ⊥ ) (cid:33) (B1) V inter ( r ) = e π(cid:15) hBN ∞ (cid:88) n = −∞ ( − ) n (cid:112) r + ( nξ + (( − ) n + 1) d ⊥ / = e π(cid:15) BN ∞ (cid:88) n = −∞ (cid:32) (cid:112) r + (2 nξ + d ⊥ ) − (cid:112) r + ((2 n + 1) ξ ) (cid:33) (B2)To calculate the Fourier transform of V ( r ), notice that ˆ d r e i k · r (cid:112) r + r = 1 π ˆ d r d z e i k · r r + z + r = 1 π ˆ d r e i k (cid:107) · r r + r = 2 ˆ d r d cos θ r e ikr cos θ r + r = 2 ik ˆ d r rr + r (cid:0) e ikr − e − ikr (cid:1) = 2 ik ˆ ∞−∞ d r re ikr r + r = 2 πk e − kr = ⇒ V intra ( q ) = e π(cid:15) πq ( e qd ⊥ − e − qξ )( e − qd ⊥ − e − qξ )1 − e − qξ (B3)= ⇒ V inter ( q ) = e π(cid:15) πq e qd ⊥ ( e − qd ⊥ − e − qξ ) − e − qξ (B4)When qd ⊥ (cid:28)
1, it is clear that V intra ( q ) ≈ V inter ( q ) ≈ e π(cid:15) πq tanh (cid:18) qξ (cid:19) Appendix C: Energetics at w /w = 0 . FIG. S2. The energy of various states with the trial function in Eqn. 30 at w /w = 0 .
3. The energies are normalized by U = e / (4 π(cid:15)L m ). The figure includes the energies of four different states: C T broken state, C T nematic state, C T period-2stripe state, and the semi-metal state obtained by minimizing the kinetic energy only. The QAH state is clearly the groundstate, with 0 . U ≈ C T symmetric states, suggesting that the system with w /w = 0 . Appendix D: Parameterization of C T Stripe Phase
In this section, we discuss how to parameterize the C T symmetric period-2 stripe phase. As mentioned in themain text, we consider the states that can be written in the product form so that the Wick’s theorem can be applied,ie. | Ψ s (cid:105) = (cid:89) k ∈ [0 , q ∈ [0 , / χ † ( q, k ) χ † ( q, k ) |∅(cid:105) (D1)with χ † i ( q, k ) = u i ( q, k ) b † + ,q,k + u i ( q + 12 , k ) b † + ,q + ,k + v i ( q, k ) b †− ,q,k + v i ( q + 12 , k ) b †− ,q + ,k i = 1 , . (D2)Since the many-body state is C T symmetric, the effective Hamiltonian H seff in Eqn. 64 for the stripe phase is also C T symmetric. The two vectors ϕ i ( q, k ) = (cid:18) u i ( q, k ) , u i ( q + 12 , k ) , v i ( q, k ) , v i ( q + 12 , k ) (cid:19) T (D3)for i = 1 and 2 can be chosen to be the eigenstates of H seff and thus also be C T symmetric. This leads to theconstraints that v i ( q, k ) = u ∗ i ( q, k ) and v i ( q + , k ) = u ∗ i ( q + , k ). Therefore, we can write the vectors as ϕ i ( q, k ) = (cid:32) ψ i ( q, k ) ψ ∗ i ( q, k ) (cid:33) with ψ i ( q, k ) = (cid:32) u i ( q, k ) u i ( q + , k ) (cid:33) . The normalization gives ψ † i ( q, k ) ψ i ( q, k ) = 1 /
2. Furthermore, the two vectors ϕ ( q, k ) and ϕ ( q, k ) are orthogonalto each other, leading to another constraint that ψ † ( q, k ) ψ ( q, k ) + c.c. = 0. Therefore, ψ † ( q, k ) ψ ( q, k ) is purelyimaginary. We write ψ = √ | ˆ n, ↑(cid:105) , meaning it is the eigenstate of the operator ˆ n · σ with the eigenvalue of 1 / √ ψ is also important, because the state ϕ is changed with an additional phase addedto ψ . We also write ψ ( q, k ) = √ ( iα | ˆ n ↑(cid:105) + β | ˆ n ↓(cid:105) ) with | α | + | β | = 1. It is obvious that α must be real to ensurethe orthogonality condition. Thus, ψ = e i π ˆ n (cid:48) · σ ψ = i ˆ n (cid:48) · σ ψ with ˆ n (cid:48) being an arbitrary three dimensional unitvector.However, this configuration for ϕ and ϕ still contains redundancy, since any real orthogonal transformation thatmixes these two vectors leaves the two dimensional subspace unchanged. This transformation can be written as ψ (cid:48) = cos ω ψ − sin ω ψ = cos ω ψ − i sin ω n (cid:48) · σ ψ = e − i ω ˆ n (cid:48) · σ ψ ψ (cid:48) = i ˆ n (cid:48) · σ ψ (cid:48) . Thus we can rotate the spinor ψ and ψ around the ˆ n (cid:48) unit vector, and obtain an equivalent many-body state. Since ψ is acquired by rotating ψ by π around the same unit vector n (cid:48) , we can always choose the spinor ψ = √ | ˆ n ↑(cid:105) sothat ˆ n has the same azimuthal angle as ˆ n (cid:48) . Therefore, we can parameterize these two unit directions asˆ n = (sin( θ + θ ) cos φ , sin( θ + θ ) sin φ , cos( θ + θ )) and ˆ n (cid:48) = (sin θ cos φ , sin θ sin φ , cos θ ) . This gives the spinors ψ and ψ as ψ = 1 √ e iφ (cid:32) cos θ + θ sin θ + θ e iφ (cid:33) ψ = i ˆ n (cid:48) · σψ = i √ e iφ (cid:32) cos θ − θ sin θ − θ e iφ (cid:33) . (D4)Therefore, the state in Eqn. D1 can be described by 4 parameters living on S × S at every momentum. The projectorcan be written as ϕ ϕ † + ϕ ϕ † = 12 θ cos θ e − iφ sin θ cos θ − e iφ sin θ sin θ e i ( φ +2 φ ) cos θ sin θ e iφ cos θ sin θ − cos θ cos θ e i ( φ +2 φ ) cos θ sin θ e i ( φ + φ ) sin θ sin θ − e − iφ sin θ sin θ e − i ( φ φ ) cos θ sin θ θ cos θ e iφ sin θ cos θ e − i ( φ +2 φ ) cos θ sin θ e − i ( φ + φ ) sin θ sin θ e − iφ cos θ sin θ − cos θ cos θ . (D5)1 Appendix E: Self-Consistent Equation
In this section, we present the detailed derivation of the self-consistent equations 61 and 67, as well as the expressionof the operator F in Eqns. 62 and 66. First, we consider the state with translation symmetry. For a state given byEqn. 57, the fermion correlation function is (cid:104) c † ( r ) c ( r (cid:48) ) (cid:105) = (cid:88) q,k (cid:88) α,β g ∗ α,q,k ( r ) g β,q,k ( r (cid:48) ) M αβ ( q, k ) with M ( q, k ) = (cid:32) | u ( q, k ) | u ∗ ( q, k ) v ( q, k ) v ∗ ( q, k ) u ( q, k ) | v ( q, k ) | (cid:33) , (E1)where g α,q,k ( r ) = (cid:104) r | φ α ( q, k ) (cid:105) is the wavefunction of the Chern Bloch state at the momentum of ( q, k ) with the Chernnumber α . Since the trial state is a product state, by Wick’s theorem (cid:104) ˆ V int (cid:105) = 12 ˆ d r d r (cid:48) V ( r , r (cid:48) ) (cid:104) c † ( r ) c † ( r (cid:48) ) c ( r (cid:48) ) c ( r ) (cid:105) = 12 ˆ d r d r (cid:48) V ( r , r (cid:48) ) (cid:2) (cid:104) c † ( r ) c ( r ) (cid:105)(cid:104) c † ( r (cid:48) ) c ( r (cid:48) ) (cid:105) − (cid:104) c † ( r ) c ( r (cid:48) ) (cid:105)(cid:104) c † ( r (cid:48) ) c ( r ) (cid:105) (cid:3) (E2) (cid:104) ˆ H kin (cid:105) = ˆ d r d r (cid:48) t ( r , r (cid:48) ) (cid:104) c † ( r ) c ( r (cid:48) ) (cid:105) (E3) E = (cid:104) ˆ H kin (cid:105) + (cid:104) ˆ V int (cid:105) (E4)Notice that u ( q, k ) and v ( q, k ) show only in the matrix M αβ ( q, k ). After some calculations, we obtain (cid:32) δEδu ∗ ( q,k ) δEδv ∗ ( q,k ) (cid:33) = H eff ( q, k ) (cid:32) u ( q, k ) v ( q, k ) (cid:33) (E5)( H eff ( q, k )) αβ = ˆ d r d r (cid:48) V ( r , r (cid:48) ) (cid:2) (cid:104) c † ( r ) c ( r ) (cid:105) g ∗ α,q,k ( r (cid:48) ) g β,q,k ( r (cid:48) ) − (cid:104) c † ( r (cid:48) ) c ( r ) (cid:105) g ∗ α,q,k ( r ) g β,q,k ( r (cid:48) ) (cid:3) + ˆ d r d r (cid:48) t ( r , r (cid:48) ) g ∗ α,q,k ( r ) g β,q,k ( r (cid:48) ) (E6)Note that we have used the relation V ( r , r (cid:48) ) = V ( r (cid:48) , r ).Therefore, ( H eff ( q, k )) αβ can be written as (cid:104) φ α ( q, k ) |F| φ β ( q, k ) (cid:105) with F ( r , r ) = (cid:104) r |F| r (cid:105) = δ ( r − r ) ˆ d r V ( r , r ) (cid:104) c † ( r ) c ( r ) (cid:105) − V ( r , r ) (cid:104) c † ( r ) c ( r ) (cid:105) + t ( r , r ) , (E7)where the fermion correlation is given by Eqn. E1. It is clear that the operator F is independent of the momentum( q, k ), and thus has the winding number of 0. Since V ( r , r ) = V ( r , r ) and t ( r , r ) = t ( r , r ), F is a hermitianoperator.In addition, if the state is C T symmetric, (cid:104) c † ( r ) c ( r ) (cid:105) = (cid:104) c † ( − r ) c ( − r ) (cid:105) and (cid:104) c † ( r ) c ( r ) (cid:105) = (cid:104) c † ( − r ) c ( − r ) (cid:105) ∗ = ⇒ F ( r , r ) = F ∗ ( − r , − r ) (E8)Combined with hermiticity of the operator F , we conclude F ( r , r ) = F ( − r , − r ).For the period-2 stripe state described by Eqn. D1, we can follow the same approach to find the expression of F .For this purpose, define the wavefunction to be g q,k ( r ) = (cid:18) (cid:104) r | φ + ( q, k ) (cid:105) , (cid:104) r | φ + ( q + 12 , k ) (cid:105) , (cid:104) r | φ − ( q, k ) (cid:105) , (cid:104) r | φ − ( q + 12 , k ) (cid:105) (cid:19) As a consequence, the fermion correlation function can be expressed as (cid:104) c † ( r ) c ( r (cid:48) ) (cid:105) = (cid:88) k ∈ [0 , q ∈ [0 , /
2) 2 (cid:88) i =1 4 (cid:88) α,β =1 g ∗ α,q,k ( r ) g β,q,k ( r (cid:48) ) ( M i ( q, k )) αβ with M i ( q, k ) = (cid:0) ϕ i ( q, k ) (cid:1) ∗ (cid:0) ϕ i ( q, k ) (cid:1) T , (E9)2where ϕ i ( q, k ) is defined in Eqn. D3. The energy can be calculated by applying the Wick’s theorem, and we obtainthe same expression as Eqns. E2 – E4. As a consequence, (cid:16) δEδu ∗ ( q,k ) δEδu ∗ ( q + ,k ) δEδv ∗ ( q,k ) δEδv ∗ ( q + ,k ) (cid:17) T = H seff ( q, k ) (cid:16) u ( q, k ) u ( q + , k ) v ( q, k ) v ( q + , k ) (cid:17) T (E10) (cid:0) H seff ( q, k ) (cid:1) αβ = ˆ d r d r (cid:48) V ( r , r (cid:48) ) (cid:2) (cid:104) c † ( r ) c ( r ) (cid:105) g ∗ α,q,k ( r (cid:48) ) g β,q,k ( r (cid:48) ) − (cid:104) c † ( r (cid:48) ) c ( r ) (cid:105) g ∗ α,q,k ( r ) g β,q,k ( r (cid:48) ) (cid:3) + ˆ d r d r (cid:48) t ( r − r (cid:48) ) g ∗ α,q,k ( r ) g β,q,k ( r (cid:48) ) (E11)Thus, (cid:16) H seff ( q, k ) (cid:17) αβ = (cid:104) η α ( q, k ) |F s | η β ( q, k ) (cid:105) , and the operator F s for the stripe state has the same expression asEqn. E7, but with the fermion correlation given in Eqn. E9. Notice that F s is also hermitian. Similar to the nematicphase, if the stripe state is also C T symmetric, F s ( r , r ) = F s ( − r , − r ). In the next section, this property willbe used to derive a series of periodic properties of the effective Hamiltonian. Appendix F: Properties of the Effective Hamiltonian of the C T Stripe State
In this section, we discuss the properties of the matrix elements of H seff given in Eqn. 67 and E11. As illustrated inthe previous section, the matrix element can be expressed as Eqn. 66. With the gauge of Chern Bloch states chosenin Eqn. 16, (cid:15) ( q, k + 1) = (cid:15) ( q, k ) (cid:15) ( q + 1 , k ) = (cid:15) ( q, k ) δ ( q + 12 , k ) = δ ∗ ( q, k ) δ ( q, k + 1) = − δ ( q, k )∆ ( q, k + 1) = e πiq ∆ ( q, k ) ∆ ( q + 1 , k ) = ∆ ( q, k ) ∆ ( q + 12 , k ) = ∆ ( q, k ) ∆ ( q, k + 1) = − e πiq ∆ ( q, k )(F1)Notice that the definition of the operator F s is needed to prove ∆ ( q + , k ) = ∆ ( q, k ):∆ ( q + 12 , k ) = (cid:104) φ + ( q + 12 , k ) |F s | φ − ( q, k ) (cid:105) = ˆ d r d r g ∗ + ,q + ,k ( r ) F s ( r , r ) g − ,q,k ( r )= ˆ d r d r g − ,q + ,k ( − r ) F s ( − r , − r ) g ∗ + ,q,k ( − r )= ∆ ( q, k ) (F2)It is obvious that ∆ has the winding number of 1 around the stripe BZ.Next, we consider how the double degeneracy between the two low (high) energy bands is lifted by both δ ( q, k ) and∆ ( q, k ). Based on Eqns. 71 – 74, the first order perturbation gives H +12 ( q, k ) = ∆ cos θ δ sin θ ∆ | ∆ | + ∆ (cid:48)∗ sin θ (cid:18) ∆ | ∆ | (cid:19) (F3) H − ( q, k ) = ∆ (cid:48)∗ cos θ − δ sin θ ∆ ∗ | ∆ | + ∆ sin θ (cid:18) ∆ ∗ | ∆ | (cid:19) (F4)Applying the boundary conditions listed in Eqn. F1, we obtain H +12 ( q, k + 1) = e i πq H +12 ( q, k ) H +12 ( q + 12 , k ) = (cid:0) H +12 ( q, k ) (cid:1) ∗ (cid:18) ∆ | ∆ | (cid:19) H − ( q, k + 1) = e − i πq H − ( q, k ) H − ( q + 12 , k ) = (cid:0) H − ( q, k ) (cid:1) ∗ (cid:18) ∆ ∗ | ∆ | (cid:19) (F5)Now, we prove that the winding number of H ± around the strip BZ must be even. Consider the closed contour C = (0 , − ( , − ( , − (0 , − (0 , δθ = 1 i ˆ d q ∂ q H +12 ( q, H +12 ( q, , δθ = 1 i ˆ d k ∂ k H +12 ( , k ) H +12 ( , k ) , δθ = 1 i ˆ d q ∂ q H +12 ( q, H +12 ( q, , δθ = 1 i ˆ d k ∂ k H +12 (0 , k ) H +12 (0 , k ) . H +12 ( q, k + 1) = e πiq H +12 ( q, k ), it is easy to prove that δθ + δθ = − π . For notationalconvenience, introduce φ as the phase of ∆ . Applying the boundary condition H +12 ( q + , k ) = (cid:0) H +12 ( q, k ) (cid:1) ∗ e iφ ( q,k ) ,we obtain δθ = 1 i ˆ d k ∂ k H +12 ( , k ) H +12 ( , k ) = 1 i ˆ d k ∂ k (cid:16)(cid:0) H +12 (0 , k ) (cid:1) ∗ e iφ (0 ,k ) (cid:17)(cid:0) H +12 (0 , k ) (cid:1) ∗ e iφ (0 ,k ) = 1 i ˆ d k (cid:32) ∂ k (cid:0) H +12 (0 , k ) (cid:1) ∗ (cid:0) H +12 (0 , k ) (cid:1) ∗ + 2 i∂ k φ (0 , k ) (cid:33) = − i ˆ d k ∂ k H +12 (0 , k ) H +12 (0 , k ) + 2 ˆ d k ∂ k φ (0 , k ) = δθ + 2 ˆ d k ∂ k φ (0 , k ) (F6)Since ∆ (0 ,
1) = − ∆ (0 , ´ d k ∂ k φ (0 , k ) = (2 n + 1) π . In addition, the boundary condition H +12 (0 ,
1) = H +12 (0 , δθ = 2 mπ . Therefore, the total change of the phase is δθ + δθ + δθ + δθ = − π + 2 δθ + 2 ˆ d k ∂ k φ (0 , k ) = 2 π ( − m + (2 n + 1)) = 4 π ( n − m ) . (F7)Therefore, the winding number must be even. Similar derivation can be done for H − , and we obtain the sameconclusion.As mentioned in the text, the formula in Eqns. F3 and F4 are ill-defined at a particular momentum ( q (cid:48) , k (cid:48) ) at which∆ ( q (cid:48) , k (cid:48) ) = 0 and (cid:15) (cid:48) ( q (cid:48) , k (cid:48) ) <
0. Obviously, Eqns. F3 and F4 are well defined in the region enclosed by the closedcontours C and γ q (cid:48) ,k (cid:48) . The winding numbers on the contour C are still given by Eqn. F7. On the small contour γ q (cid:48) ,k (cid:48) ,Eqns. F3 and F4 are dominated by their last terms. Since θ ( q (cid:48) , k (cid:48) ) = π , the winding numbers of H ± on the contour γ q (cid:48) ,k (cid:48) are ± C and γ q (cid:48) ,k (cid:48) are still even. As a consequence, the numbers of Dirac points of H ±±