Hybrid Wannier Chern bands in magic angle twisted bilayer graphene and the quantized anomalous Hall effect

We propose a method for studying the strong interaction regimes in twisted bilayer graphene using hybrid Wannier functions, that are Wannier-like in one direction and Bloch-like in the other. We focus on the active bands as given by the continuum model proposed by Bistritzer and MacDonald, and discuss the properties of corresponding hybrid Wannier functions. We then employ the method for a study of the fillings of $\pm 3$ electrons per moir\'e cell using the Hartree-Fock method. We discuss at length different regimes under which a quantized anomalous Hall effect is seen in these two fillings.


I. INTRODUCTION
Heterostructures containing moiré patterns due to incommensurations in multilayers containing graphene and other two dimensional crystals have proven to be very tunable and promising platforms for observing interesting phases that are unprecedented in commensurate graphene systems [1][2][3][4][5][6] .Twisted bilayer graphene (TBG) as the most prominent member has attracted much attention and also has given rise to numerous theoretical studies; however, still many of the different correlation induced phenomena in this system have eluded satisfactory theoretical understanding.
The most important theoretical discovery, probably, was the realization that a low energy theory, a continuum model (CM), could be effectively employed to study the single particle electronic properties of TBG at small twist angles 7 ; in fact, an analysis based on this CM resulted in the prediction of the possibility of strong correlation physics at the magic angle in the first place.Specifically in this CM, the smallness of the twist angle leads to an emergent periodicity in the system -the so called moiré lattice, which has a unit cell length growing like ∼ 1 θ ; such large periodicity in turn leads to formation of Bloch minibands.Interestingly, around the magic angle, the bands closest to the charge neutrality point (CNP) show exceptional flatness and are well separated from other bands.Further including spin and valley degrees of freedom results in eight such bands in total.Since these bands are flat, the correlation between them can play an important role and give rise to interesting correlated phases and thus should be taken into account properly.A possible theoretical approach to this end, is to consider an interacting model consisting of these active bands only, treating the remote bands as inert; we will be taking this route in this work and introduce a basis for the study of strong interactions.
Experimental observations of correlation induced insulating states have been reported in commensurate fillings of these active bands, along with superconducting behav-ior for fillings close to these commensurate values [1][2][3][4][5][6]8,9 . Motvated by these experimental observations, here we pursue a theoretical model consisting of the subspace of the active bands only, in which electronic interactions are also projected onto this subspace; these interactions are local and thus working with local representations of the subspace spanned by active bands is desirable.However, as is well known, a faithful representation preserving manifest symmetries of the active bands using fully localized Wannier functions is difficult 10,11 .Having this in mind, in this work, we work with Hybrid Wannier Functions (HWFs) which are Bloch-like in one direction and localized and Wannier-like in the other.Using this basis is a compromise between locality and symmetry/topology, noting that the wave functions are only localized in one direction, however, as is elaborated later, this ensures that important symmetries like valley and C 2 T (the intravalley symmetry that protects the moiré Dirac points) remain manifest (when not broken at the non-interacting level).Furthermore, one ends up with a quasi-one-dimensional model, with local interactions in one direction, which can be suitable for numerical methods like DMRG 12 .
As we show later, remarkably, full bands of these HWFs when maximally localized automatically exhibit a nonzero Chern number; this means that indeed a suitable collection of full bands of such states can display quantized anomalous Hall effect (QAHE), a phenomenon that has been reported in TBG 5,6 at the filling factor of ν = +3 (we define the filling factor ν to show the number of electrons per moiré cell measured from CNP).This makes the present maximally localized HWFs a natural basis for a corresponding theoretical study.To analyze the effect of the interaction, which is evidently required for stabilizing a full band polarization in the HWF basis, we employ the self-consistent Hartree-Fock (HF) method at the two fillings ν = ±3; these are the fillings where single fully occupied HWF bands of holes or electrons can be candidate many body states respectively.
We perform two separate studies of the effect of arXiv:2007.00134v1[cond-mat.mes-hall]30 Jun 2020 electron-electron interaction; first, we examine how the locality (in one direction only) of the HWFs makes full HWF bands advantageous for the interaction energy penalty when compared with other many body states at the same filling.Specifically, we check if full HWF bands turn out to be solutions of the HF equations when interaction is considered; this ensures that such HWF band polarized states have (at least local) minimal interaction energy compared with other candidate many body states.Second, we study the stability of similar many body states in a model obtained by projection of the full Hamiltonian onto the active bands.We present numerical results on the stability of QAHE in these two settings in a wide range of parameter choices of the models.
There have been other HF studies of the continuum model at various integer filling factors, with the analysis carried out completely using the basis of original Bloch states [13][14][15][16] ; in a subset of these works the remote bands are also kept in the analysis.The present study has the advantage of working directly with a faithful semilocalized representation of the active bands, while providing a continuous description of the QAHE with and without the C 2 T symmetry of TBG.Moreover, in the present analysis, the QAHE appears naturally as polarized bands in the HWF basis and this could provide some more insight into the nature of the Chern bands responsible for this effect.A comparison between these prior HF studies and our results is presented in Appendix E.
The paper is organized as follows: first, in Sec.II, we demonstrate how the maximally localized HWFs are constructed and derive their topological properties.Then, in Sec.III, we present the HF study of the interacting model at the fillings ±3, and the stability of QAHE by varying various parameters is examined.We conclude our results in Sec.IV.

II. HYBRID WANNIER FUNCTIONS
We will be working with the continuum model introduced in Ref. 7. To take into account the two valleys, two parallel copies of the CM are considered; in each copy, we will focus on the two active bands, closest to CNP.Details of the non-interacting Hamiltonian are presented in Appendix A. The CM has two free parameters in it: i) α ∼ 1 θ w AB , which accounts for the collective effect of interlayer hopping w AB and the twist angle θ, and ii) η = wAA wAB , the ratio of the interlayer tunneling strength in AA and AB regions of the moiré lattice, which encodes how much corrugation is present in the system.We will also consider adding a sublattice symmetry breaking term ∆ σ z to the non-interacting Hamiltonian, where the Pauli matrix σ z is used to address sublattice degrees of freedom; this could account for the effect of aligned hexagonal Boron Nitride (hBN) substrates on the two sides of the TBG sample. 17We will also be using an approximation 18,19 which renders a particle-hole symme-try to the CM; this approximation becomes better at small angles, see Appendix A for details.
Equipped with the full non-interacting content of the model, one can find the Bloch states lying in the middle two active bands for each valley.We take the active bands to be well separated from the remote bands, and thus develop an active-bands-only model.Following the notation and methods introduced in Refs.20 and 21 we will Wannier transform properly chosen Bloch states in only one direction to obtain the maximally localized HWF basis as follows: where |k x ; y c , m, ξ stands for a Hybrid Wannier state, with the indices y c , m, ξ denoting the real space position in the localized direction, the band (orbital), and the valley respectively.The states on the right hand side are linear combinations of the Bloch eigenstates of the non-interacting Hamiltonian at each k: The unitary (in the band basis) matrices U are chosen at each k to ensure that maximal localization is achieved in the y direction ultimately and the procedure is detailed below.Here, a rectangular BZ is chosen as shown in Fig. 1(a) so that the k y sum needed for a one dimensional Wannier transform in ( 1) is performed at each k x .The spin index trivially doubles all manipulations here and thus is suppressed.With the above convention, the allowed values of y c form a one dimensional lattice with lattice spacing equal to half a moiré length ( 1 2 a M = a 1,y ), i.e. y c = j a M 2 where j is an integer.Note that we take this lattice to be identical for different values of k x , and so the above y c values are different from but close to the actual locations of Wannier charge centers (WCC) of HWF states (see below for more information).The transformation of the HWFs under moiré lattice translations is depicted in Fig. 1(a).
In order to obtain maximal localization, one needs to choose the matrices U in (2) properly: to this end, we will use the procedure discussed in Ref. 20 to form the parallel transport basis for the Bloch functions, an approach that is suitable for maximal localization of one dimensional Wannier functions, and in the present study should be carried out for each strip with a definite k x separately.We will use a discretization which will lead to a Bloch momenta lattice with lattice spacings b x , b y , and with N x , N y total points along the two directions.According to this prescription, at each k, the overlap matrices M kx,ky,ξ mn = u kx,ky;m,ξ u kx,ky+by;n,ξ , 4 , − 3 2 < ky < 3 2 ; note that this is contrary to the usual hexagonal choice so that the top and bottom of the BZ are identified, note that this is crucial for the usual properties of the one-dimensional Wannier transform in the y direction to hold.The equations governing the translational properties of the HWFs are also presented.(b) WCC positions (solid black lines) and single band Berry phases in the original Bloch bases (dashed red lines) of the two active bands.One of the plots corresponds to the chiral model and the other to the physical value of η.A small sublattice potential is added, ∆ = 0.19meV.The configuration of the dashed lines and the solid lines mean that the two bands carry +1 and −1 Chern numbers in the original Bloch representation and the parallel transport representation respectively.This is a robust feature present in a wide range of parameter choices.Note that kx is rescaled and instead of plotting the interval [−0.5, 0.5), equivalently [0, 1) is drawn.
are calculated, where as usual |u k;n,ξ shows the unit cell periodic part of an original Bloch function; notice that there is a small displacement in the k y direction in the ket state.Next, redefinitions of Bloch functions are made as shown in (2), with U matrices chosen in a way that the updated M matrix for all k x , k y , ξ attains a form as K kx,ky,ξ γ kx,ξ , where K is Hermitian and γ is diagonal, unitary and independent of k y .This, as discussed in Appendix B, ensures maximal localization in the y direction.
The path ordered product of all M matrices along a strip with a given k x defines its Wilson loop, whose eigenvalues are invariant under a k dependent basis change such as the one in (2).One can show that the K matrices as defined above are equal to the identity matrix to first order in b y and thus, the eigenvalues of each γ kx,ξ matrix above are directly related to the Wannier charge center positions, i.e.Wilson loop eigenvalues, in the strip given by k x .Using this fact, WCCs of HWFs as functions of k x could be found with examples drawn in Fig. 1(b).It could be seen by inspection that, regardless of the set of parameters chosen, there is a +1 winding and a −1 winding of the WCCs for the two HWF bands as k x traverses the BZ. 19,22Noting, based on the above observations, that in the parallel transport basis, the single band Berry phases along each strip with a given k x are equal to the WCC values, leads us to an important implica-tion for the parallel transport basis: given how WCCs behave as functions of k x shown in Fig. 1(b), the two Bloch bands in the parallel transport basis have Chern numbers +1 and −1.This, in other words, means that a fully filled band of maximally localized HWFs exhibits a quantized Hall response.As a result, when addressing the maximally localized HWFs, the terms band, orbital and Chern number could be used interchangeably.
In some special cases, the parallel transport basis can be found explicitly.For instance, when ∆ = 0, there is a C 2 T = σ x K symmetry of the Hamiltonian, where K is the complex conjugation operator; as shown in Appendix B, the combinations e ±iφ k,ξ √ 2 (|ψ k;1,ξ ± i |ψ k;2,ξ ), with the phases φ k,ξ appropriately chosen, form the parallel transport Bloch basis at k, where states |ψ k;m,ξ show C 2 T symmetric Bloch eigenstates.In particular, if one now sets η = 0 to obtain the chiral limit, since the two C 2 T symmetric bands are related 23 by |ψ k;1,ξ = iσ z |ψ k;2,ξ , the parallel transport basis consists of sublattice polarized states; remarkably, even with ∆ = 0 while keeping η = 0 this result holds, i.e. the parallel transport basis consists of sublattice polarized states.By numerical inspection, one can show that each of the two bands in the parallel transport basis is more concentrated on one of the sublattices to a high degree in a one-to-one fashion, even away from the chiral limit.It is worthwhile to mention that the U (4) × U (4) symmetry discussed in Ref. 16 (which states that the interaction term of the Hamiltonian is invariant under rotations of the bands with equal Chern numbers into each other) could be seen readily in the above construction of the parallel transport basis.This along with other symmetries of the CM as seen in the HWF basis are discussed in length in Appendix C.
The HWF basis naturally defines the problem in the geometry of a cylinder.The HWFs form ring shaped wires around the cylinder, since these wave functions are localized in one direction and extended in the other.Each wire is identified with a y c , and is composed of states with different values for their k x , band number, valley number and spin (see Eq. ( 1)).At the non-interacting level, hopping occurs between states in separate wires if they have identical k x , valley number and spin (See Appendix C for details).This hopping decays as the distance between wires along the cylinder is increased.Based on this HWF construction, in the next section we will present a HF study of a model consisting of active bands only with a total Hamiltonian of the form: H kin contains the single particle terms in the Hamiltonian induced by the CM, i.e. the hoppings between different wires as mentioned above.The remaining two terms represent effects of interactions: they are both proportional to e 2 / , where e is the electron charge and is the dielectric constant, and thus vanish in the non-interacting limit.H MF,0 , which is quadratic in fermion operators, is responsible for two separate effects: it takes the effect of filled remote bands into account at a mean field level and it also serves to avoid a double counting of HF terms that are already taken into account in H kin 13,16 (see the discussion at beginning of the next section for more details).Turning to the interaction term H int , we have chosen the electron-electron potential to have a screened coulomb form as , which is further projected onto the active bands.The interaction retains its normal-ordered density-density form with respect to spin, sublattice, layer and valley indices (more details are presented in Appendix C).Note that due to the locality of HWFs, the electron-electron interaction between the wires drops as the distance between them is increased, and thus the total Hamiltonian is local in the direction along the cylinder.
We conclude this section by some remarks regarding the parameter values and conventions used: via dividing the energies and lengths by v F k θ and 1 k θ respectively, we have made them dimensionless, where In this notation, we define the dimensionless interaction strength parameter , where A is the area of a moiré unit cell.Numerically g int = 1.01 0 , and thus, a choice of = 7 0 results in g int = 0.14 as an example.
The model introduced above comprises bands (in the parallel transport basis) that carry nonzero ±1 Chern numbers; thus a quantized Hall signal can be observed at integer filling factors if with some interaction induced effect, a suitable valley and band polarization in the system occurs.As a result, it is natural to utilize the present model to study the physics of quantized anomalous Hall effect (QAHE) seen in some samples of twisted bilayer graphene 5,6 .We will do so in the following section for the two fillings ν = ±3.

III. QUANTIZED ANOMALOUS HALL EFFECT IN TWISTED BILAYER GRAPHENE
In this section, we present two separate HF studies in which different choices of H MF,0 are utilized.We focus on the filling ν = ±3 and explore the stability of QAHE phase in these two different schemes.Before we go into the detail, we first discuss how the HF procedure is carried out in general.
The HF procedure is implemented as follows: we fix the filling and seek a Slater determinant many body state, composed of single particle states |ψ l = α ψ l,α |α , that minimizes the expectation value of the Hamiltonian (4), where α, β, . . .denote the HWF basis indices k x , y, ξ, m, s (the states |α will be normalized in this section).One seeks |ψ l states by transforming the Hamiltonian (4) written in the form into a single particle HF Hamiltonian, wherein the interaction term is transformed into In the above, a, b, . . .(contrary to α, β, . ..) show the HWF indices except k x .Notice that we have dropped the x subscript from k x and will do so from now on; the k dependent matrices P have the form P (k) aa = l ψ * l,ka ψ l,ka .It has, furthermore, been assumed that the translational symmetry around the cylinder is not broken.
The above HF Hamiltonian depends on its own eigenstates and thus we aim to obtain them iteratively: starting from a well chosen initial many body state, at each iteration step, P matrices are updated using the eigenstates found in the previous step; a ν dependent number of these eigenstates with lowest eigenvalues participate in forming the P matrices.The resulting HF Hamiltonian is then diagonalized to yield the updated set of eigenvalues and eigenstates.This procedure is continued until convergence is achieved.We obtain the sought HF many body state as a slater determinant of the converged eigenstates with lowest HF eigenvalues.Moreover, the nearby eigenvalues above and below the "Fermi energy" could be used to give estimates of the actual energies needed for adding or removing an electron at this filling (Koopmans' theorem 24 ) 25 .
The two approaches mentioned at the beginning of this section are taken into account by two different choices for H MF,0 in the Hamiltonian (4).In the first study, Sec.III A, we examine the motivation with which the HWF basis was introduced: the interaction energy of different many body states are compared with H MF,0 = 0.In particular, the energy of the state that is described as a full band of electrons (ν = +3) or holes (ν = −3) in the HWF basis is compared with other HF many body states.Note that this choice of H MF,0 = 0 results in a competition between the interaction energies and the band structure energies as given by the CM; the latter, which could also be viewed as the hopping term in the HWF basis, is kept in the analysis so that one attains a measure for defining strong and weak interaction regimes.
In the second study, in Sec.III B, on the other hand, we take 26 where the α, α summation is done over all states in the active bands, but the partial summation over β, β (indicated by the prime on the sum) ranges only over those states in the active bands that are below the CNP of the CM.Note that the latter states when written in terms of the HWF basis will not be band diagonal.By taking H MF,0 to have the form in Eq. ( 7), we are taking two separate effects into account: first, a mean field potential induced by the filled remote bands.The second effect, instead, has to do with the fact that within HF, the electron/hole dispersion will only agree (at best) at one filling with the dispersion given by the term H kin .We take that point to be the CNP of the CM bands in the second study, i.e. we assume that the CNP dispersion given by the CM, describing single electron or single hole excitation energies on top of the CNP, is unaltered by HF (see Appendix C for discussion).In order for this to be true, a HF effect of all filled bands (including remote and active bands) at the CNP is subtracted.The combination of these two effects results in a cancellation of the mean field effect of the filled remote bands and thus one ends up with the form in (7) with only the mean field effect of active filled bands subtracted.
In the next two subsections, we present our numerical results corresponding to these two studies.

A. First study
In this subsection, we consider a model in which H MF,0 = 0, wherein a competition between electronelectron interactions and the noninteracting hopping in the HWF basis enables us to tune the model into and out of the strong coupling regime.Previous studies, working on generic models similar to the one used in this subsection, have shown analytically that in strong coupling limits, valley polarization in these two filling factors is expected 27,28 .Here, we present a more thorough HF study of the Hamiltonian, trying to identify different regimes in which QAHE could be achieved.
In a given setting, we say that the QAHE is stabilized through HF if two requirements are met: i) if we start from a fully spin-valley-band-polarized state, and the HF iterative steps lead to a final HF state that is only perturbatively different from the spin-valley-band-polarized state, and ii) the final HF solution properties in large enough systems do not change considerably as the system size is varied.
This means that the QAHE state is at least a minimum of energy; we have also tried perturbing the final HF state in different ways to examine the stability of the HF solutions and we have observed that the final many body states generally show a high level of stability.In each setting we start with strong interactions first and see if the QAHE state is stabilized, and then continue to lower the interaction strength.We will, furthermore, use periodic boundary conditions along the cylinder.For the numerical results presented in this paper, the system is chosen to have N x = N y = 20.

1.
η = 0, short range interaction We now start to present our numerical HF results.We construct the basis of HWFs by forming the parallel transport basis on a finite lattice in k-space as discussed in Appendix B, and then make a Wannier transform along y for each k x .
The chiral model, i.e. when η = 0, is first considered, in which absolutely flat bands are achieved at the magic angle.We start with the small ξ limit so that the electronelectron interaction V int (r) is very short ranged.More realistic longer range interaction is considered later.Moreover, we also take the twist angle different from but close to the magic value so that the bands exhibit a nonzero small width.
In the first setting outlined above, or concretely with the choices η = 0 and ξ = 0.17a M , numerical analysis shows that the QAHE is generically stabilized at ν = −3 at large interaction strength, see Fig. 2(a,b), where the density of different flavors along with HF eigenvalues (energies) are shown for an instance where the interaction plays the dominant role.Note that because of the nature of the HWF basis, this is a C 2 T broken many body state, despite the fact that this symmetry is present at 0.58, 1.9 0.58, 19 0.58, 0 0.59, 1.9 0.57, 1.9 0.56, 1.9 0.55, 1.9 (d) FIG. 2. The filling ν = −3, η = 0, and ξ = 0.17a M plots within the first study, the angles are chosen around the magic value (α 0 = 0.586 corresponding to θ = 1.05 • ) (a) Density of electrons (number of electron per unit length of the cylinder) vs kx (momentum across the cylinder) for different flavors where s, ξ, m stand for spin, valley, and band.α = 0.58, ∆ = 0 and furthermore g int = 0.05 have been chosen here, this value for the latter corresponds to being close to strong interaction limit since the band width is very small.Almost full polarization is seen here.(b) The HF eigenvalues (energies), with the same parameters as described in (a).The Fermi surface is shown with the dashed red line, one can observe a HF gap which will be used as a criterion for determining whether QAHE has been stabilized under HF iterations or not.(c) The same plot as in (b) with ∆ = 1.9meV.One can observe a second gap that is formed above the interacting one, due to the relatively large ∆ chosen; it separates the states belonging to the opposite sublattices.(d) The HF gap divided by g int as a function of g int for several parameter choices.It can be inferred that approaching the magic angle and making ∆ larger makes the QAHE more HF stable.
the noninteracting level.We define a HF gap as the lowest unoccupied HF eigenvalue minus the highest occupied eigenvalue, this quantity when divided by the interaction strength g int serves as a good qualitative measure of whether and to what extent the polarized state is stabilized under HF.A plot of such gaps as functions of interaction strength for several parameter choices is shown in Fig. 2(d).Note that the polarized state continues to exhibit HF stability as the interaction is lowered but becomes unstable when the interaction energy per particle becomes roughly comparable to the band width.Moreover, we consider a range of ∆ from small to large values (always smaller than the noninteracting gap to remote bands); as shown in Fig. 2(d), regardless of the value of ∆, large interaction strength stabilizes the QAHE, while in the range of small interaction strength, larger ∆ results in a more stable polarization.In addition, at intermediate interaction strength, a second gap between HF eigenvalues, apart from the one induced by the interaction, is visible due to the large sublattice potential and scales with ∆ (see Fig. 2(c)); obviously, we will keep track of the former to study stability of QAHE.Also, as is expected and also shown in Fig. 2(d), tuning the twist angle away from the magic value results in weaker stabilization of QAHE and generally larger interaction is needed to stabilize the QAHE.
At the filling of +3, on the other hand, starting from large values of interaction strength, with the present settings, the QAHE state is not stabilized.However, upon decreasing the interaction strength, interestingly, when the interaction energy per particle becomes comparable to the band width, a narrow interval of interaction strength allows for the QAHE to be stabilized although it gets unstable again for smaller interactions (see Fig. 3(a)).This observation holds true irrespective of the value of ∆.
The above discrepancy between the two filling factors indicates that there is a particle-hole asymmetry in the system with the current choice of the Hamiltonian, although the non-interacting Hamiltonian is chiral and thus particle-hole symmetric with and without ∆.This asymmetry could be understood by noting the following fact within the active-bands-only model we have chosen to work with here, i.e. the choice of H MF,0 = 0: starting from the extreme cases, there is a difference between a single electron at ν = −4 and a single hole at ν = +4, in that, the hole senses an additional potential due to the presence of eight full bands of electrons.In the same fashion, a single hole senses an additional k-dependent potential at ν = +3 when compared with an electron at ν = −3, and thus some k values in the hole bands could be preferred over others; more details can be found in Appendix D. This single hole potential is interaction induced and thus becomes stronger as the interaction is raised.One can argue that destabilization of QAHE in the strong interaction limit of the filling +3 presented above happens exactly due to this potential; holes prefer to occupy some momenta more than others.As we will see below, using a longer range interaction could weaken this asymmetry.
Let us mention two important points regarding the particle hole transformation of the many body state here before moving on: our choice of H MF,0 = 0 here means that single electron excitations on top of ν = −4, receive no HF correction in their dispersion.Had we chosen another form for H MF,0 , so that the holes at ν = +4 experienced no change in dispersion from the CM, we would have gotten the same theory but with electrons replaced with holes; in other words, using this prescription for H MF,0 will yield the particle hole transformed version of the present model with H MF,0 = 0. Additionally here we only discussed the model at η = 0, where there is a chi- FIG. 3. Normalized gaps as functions of g int with longer range interaction and nonzero η in the first study.α = 0.58, and ∆ = 1.9meV have been chosen.(a) Here both of the fillings are considered still in the chiral limit (η = 0).The QAHE is more stable at larger screening length; interestingly, it is stabilized even for ν = +3 for large interaction, if large enough screening length is chosen.Note that this feature is lost for larger η values and in particular η phys = 0.8, as discussed in the main text.(b) Gaps are drawn for ν = −3 with ξ = 1.2aM here, as η is increased.Increasing η makes QAHE less stable until it is not stabilized at all even at large g int at η ≈ 0.9 − 0.95.
ral symmetry in the model, while for generic η, there is an approximate particle hole symmetry in the CM which plays a similar role.With this particle hole symmetry one can repeat the above considerations for nonzero η as well, i.e. show that the symmetry between holes and electron at the two fillings ±3 is broken within the present model and also that in a particle hole transformed version of the model, holes will play the role of electrons (see Appendix D for details).
2. η = 0, longer range interaction Upon using longer range interactions, which are more realistic, some of the results presented above are altered: a longer range interaction does not change the picture at ν = −3 much, i.e. with the use of longer range interactions, the QAHE is still stabilized at large interaction and stability is lost at small enough interaction strength (see Fig. 3(a)).However at the filling +3 the effect is more remarkable: the narrow range of the polarized states is made wider.As can be seen in Fig. 3(a), above some intermediate screening length, even at large interaction the QAHE state is stabilized.Generally, we have observed that increasing ξ makes QAHE more stable.
3. η = 0, away from the chiral limit We take another step toward making the model more realistic by choosing η to be nonzero and increasing it to the physical value η phys ≈ 0.8 29 ; the physics at ν = −3 stays similar to a high extent even up to η phys .However, at larger η, i.e. η ≈ 0.9 − 0.95, one starts to observe that the HF iterations do not stabilize the QAHE at this filling even with largest interactions.This means that the spinvalley polarized state ceases to be a local minimum in energy (among Slater determinant states) even when the interaction plays the dominant role.As can be observed in Fig. 3(b) the HF gap becomes smaller as η is increased.For the filling of +3, on the other hand, we observed that although large interaction of long enough range stabilizes QAHE at small η, for larger η and in particular for the physical value this ceases to be true no matter how long range the interaction is chosen.
Next, we turn our attention to a second study with a projected Hamiltonian.

B. Second study
In this subsection we work with a Hamiltonian that is obtained by projecting an interacting Hamiltonian onto the subspace of active bands, and the zero point of the HF approach is chosen to be at the CNP of the moiré bands, i.e. we will use Eq.(7).Unlike the previous case, this choice results in a particle hole symmetry between the many body states at the two filling factors +ν and −ν (see Appendix D for details), and therefore we will focus on ν = −3 only in this study.Note that this particle hole symmetry is present regardless of the value of η and is reflected in the HF spectrum.As an illustration, we present two sets of converged HF energies shown in Fig. 4(a) with η = 0 and η = 0.8.The HF energies at ν = +3 and ν = −3 are related by the particle hole transformation which in particular involves a k x → −k x transformation.Notice that for η = 0, these two sets of HF energies are also related by the chiral symmetry.
Fig. 4(b) shows the HF gap as a function of interaction strength for several parameter choices.We observe that for large interaction strength at small η, QAHE is stabilized under HF iterations but this does not happen for larger η.In particular, the QAHE phase is HF stable at η phys = 0.8 for a small windows of parameter choices and is absent when η = 0.5.We generically see a bump in the rescaled gap as shown in Fig. 4(b) when interaction strength is comparable to the noninteracting energies.This occurs due to a partial cancellation between quadratic terms of the Hamiltonian that arise from H kin and H MF,0 ; it is indeed this same effect that gives rise to the narrow window exhibiting QAHE at η phys .Furthermore, the destabilization of the QAHE at larger η values for large g int is attributed to the fact that H MF,0 has quadratic terms that scale with g int , and these terms prefer states with particular k x values over others.
The above results have focused on the ∆ = 0 limit.One can also consider finite ∆ or even the large ∆ limit.The latter limit is defined by requiring that the noninteracting gap due to ∆ is not closed by the interaction induced effects, hence schematically g int ∆.But this limit also results in similar behavior of the HF stability to the present model and we will not present these  In all these four subplots, we take α = 0.58, ∆ = 0, ξ = 0.5a M .The interaction strength is chosen as g int = 0.05 and 0.002 for the left and right columns respectively.Plots in each column have the same parameter values except for the filling factor: the first row corresponds to ν = −3 and the second row to ν = +3.It can be seen that the HF energies at the two fillings ν = ±3 can be transformed into each other by the particle hole symmetry of the CM, which in particular needs kx → −kx.It can be seen that for η = 0, the two sets of HF energies are also related by the chiral symmetry of the CM.(b) The gaps as functions of g int .α = 0.58, ∆ = 0 are chosen for this plot and η and ξ are varied.
numerical results here.

Discussion
We have considered two different models with different H MF,0 , in a manner that the HF zero point is chosen at ν = −4 and ν = 0 (CNP) respectively.The former case, as is discussed in Appendix D, is actually also related to a model with HF zero point choice of ν = +4, by a particle hole transformation.More specifically, the many body states at the filling factor ν that are obtained using the HF zero point of ν = −4, with an appropriate replacement of electrons with holes, are equivalent to states at filling factor −ν, if the HF zero point is moved to ν = +4.On the contrary, the choice of ν = 0 as the zero point, i.e. the case in the second study, is always particle-hole symmetric.
We would like to emphasize that the particle hole symmetry discussed above is expected to be broken even at the noninteracting level when actual physical effects like lattice relaxation are taken into account.However, we should bear in mind that if the particle hole symmetry is not broken at the noninteracting level, the interactions will also keep it intact.On the other hand, the current experiments exhibiting QAHE 5,6 only observe the effect at ν = +3, and not at ν = −3, which is an indication of particle hole symmetry breaking.On a phenomenological level, this makes us speculate that among the three different cases discussed above, the HF zero point choice of ν = +4 is probably most relevant to the physics seen in the samples exhibiting QAHE.Let us mention that ultimately within the framework of this paper, we cannot argue in favor of any of the above three choices.However we note that a definitive answer to this issue needs further study of several other effects that are neglected here; in particular, consideration of the particle hole symmetry breaking effects (such as lattice relaxation as mentioned above) and also more careful treatment of the effects of the filled remote bands could play decisive roles in determining which of the above choices (if any) could serve as a consistent physical model describing the relevant physics.

IV. CONCLUSION
To summarize, we introduce the hybrid Wannier basis in the continuum model of TBG and study the strong interaction effect by using the self consistent Hartree Fock approximation.We focus on the filling factors ±3 and investigate the stability of QAHE phases at these two fillings.Interestingly, we observe that stability of the QAHE depends crucially on the zero point choice of the HF dispersion.In the range of physically relevant parameter choices we see that the QAHE is most robustly stabilized at large interaction strength under HF for the zero point choices of ±4, and the corresponding filling factors of ν = ±3.We note that the QAHE is observed in experiments on magic angle TBG at ν = +3.Moreover, we numerically observe that reducing the sublattice potential, reducing the screening length and tuning away from the magic angle, generically make the QAHE less stable.In particular, the weakened stability by reducing the sublattice potential is consistent with experiments 5,6 , which have observed the effect in TBG samples with aligned hexagonal Boron Nitride substrates, which is believed to induce a sublattice potential.
Further development of the present method can be envisioned.The HWF basis we introduced in this paper might be used to find possible fractional phases at noninteger fillings and other interesting phases at integer filling factors.Another possible application, for which HWFs are particularly well suited, is to address situations containing spatially varying configurations such as domain walls between different symmetry broken states or states with one-dimensional "stripey" translational symmetry breaking.We leave these directions for the future study.from the CNSI, MRL: an NSF MRSEC (DMR-1720256) and NSF CNS-1725797.The research of KH and LB was supported by the NSF CMMT program under Grants No. DMR-1818533, and we benefited from the facilities of the KITP, NSF grant PHY-1748958.XC acknowledges support from the DARPA DRINQS program.constant k x line, so that all (except for the last one completing the 1D loop) M matrices become Hermitian; this is done by the Gauge transformation u kx,ky;1,ξ , u kx,ky;2,ξ → u kx,ky;1,ξ , u kx,ky;2,ξ • W V † ky−by . . .W V † ky,0 , (B1) where • denotes a matrix multiplication in the space of bands, and the k x , ξ indices on W and V matrices are suppressed.In traversing the BZ in the y direction once, one is able to define the accumulated matrix note that this matrix naively gives the prescription for a change of basis at k y,0 , the point one started with.However, we would like to end up with the state we started with so that a smooth Bloch basis is defined throughout the 1D Brillouin zone.One can achieve this if a series of actions are taken: at all k y points, a unique basis change is made with the unitary matrix that diagonalizes Λ kx,ξ , i.e. the matrix V λ , where V † λ ΛV λ = λ with λ a diagonal matrix.This last basis change updates all of the M matrices (except the last one at −k y,0 − b y , more on this below) to have the form of a Hermitian matrix.The Hermitian matrices are proportional to unity to first order in b y , and this ensures that Ω OD shown above vanishes to first order in lattice spacings.However, there remains band-diagonal total Berry phases in this new basis which are invariant under single band gauge transformations; these are the inverses of the eigenvalues of the Λ matrix defined above and are at this stage accumulated in the last M matrix, i.e. at −k y,0 − b y .One should make a band-diagonal gauge transformation (a phase redefinition) to ensure that this Berry phase is distributed evenly along the one-dimensional Brillouin Zone to make Ω D vanish as well.This final (band diagonal) gauge transformation results in the final form Kγ for the M matrices, with a Hermitian K and a diagonal unitary γ for the M matrices.
Note that an evenly distributed Berry phase means that γ kx,ξ is independent of k y and in fact equal to λ − 1 Ny .Furthermore, bear in mind that the matrix Λ † is equal to the path ordered product of M matrices to first order in b y for each k x and thus is equal to the Wilson loop at k x to this order.Noting that eigenvalues of the Wilson loop operators are related to the WCCs of the final bands means that eigenvalues of γ kx,ξ take the form e 2πi Ny y kx ,n , where y kx,n denotes the Wannier charge centers at k x in units of 1 2 a M .The U matrices defined in Eq. ( 2), can be explicitly derived as: where all right hand side matrices are evaluated at k x , ξ.
Finally, we discuss further the special cases mentioned in the main text: • In the case where ∆ = 0, due to the C 2 T symmetry of the Hamiltonian, one can work with Bloch eigenstates of Hamiltonian that are C 2 T symmetric.Any inner product of two C 2 T eigenstates is real; this means that the M matrices have the form exp [iµ y m k,ξ b y ] + O(b 2 y ), where µ y acts in the two dimensional band space.Thus every SVD operator V W † could be taken to be equal to y ) and furthermore V λ could be taken to be the matrix that diagonalizes µ y .Additionally, the integrals ± dk y m k,ξ yield the single band total Berry phases of the two bands in the parallel transport basis which should be distributed evenly along the strip with k x .All this means that the states e ±iφ k,ξ √ 2 , form the parallel transport basis, if the phases are chosen properly to distribute the single band Berry phases evenly along the y direction.
• In the case of η = 0, regardless of the value of ∆, the sublattice polarized states form the parallel transport basis.One can argue for this as follows: starting from ∆ = 0, we note that states with opposite sublattice polarizations automatically have zero contribution to Ω OD .Suitable single band gauge transformations are furthermore needed to minimize Ω D as well.On the other hand, we know that the two bands in the chiral limit are related by 23 : |ψ k,ξ,1 = i σ z |ψ k,ξ,2 .This means that adding the term σ z ∆ to the Hamiltonian does not change the subspace of active bands.And thus previously found sublattice polarized states still represent the active bands subspace, and with suitable single band phase redefinitions will form the parallel transport basis.
It is important to note that addition of ∆ does not change Wilson loop matrices for each k x , and thus the phases chosen for ∆ = 0 in the parallel transport basis remain valid choices for nonzero ∆ as well.
Second, we discuss the terms shown in the main text by H MF,0 .Although remote bands are not treated as dynamical, a proper projection of the interacting Hamiltonian onto the active bands needs inclusion of an induced mean field potential due to the filled remote bands on the electrons in the active bands.This contribution will have the form: where α, α run over active bands and β, β run over remote bands below CNP.
In addition to that, as discussed in the main text, we have taken the zero point of the HF to be at the CNP of the moiré bands.This means that at CNP, the single electron/hole dispersions as given by the CM should be unaltered under HF.In order for this to be true, we subtract the HF effect of the moiré CNP noninteracting state from the Hamiltonian.The addition of these two effects will result in the form given in Eq. ( 7) in the main text for H MF,0 .
There is a subtlety in the projection approach outlined above; with the above projected model at hand, we have considered changing the interaction strength in our study presented in the main text, this alters the coefficients of both H int and H MF,0 (a change in the dielectric constant, for example, could result in this effect).However, such a change will result in a different single electron/hole potential according to (C7); in particular, the single layer Fermi velocity v F and the interlayer tunneling parameters w AA , w AB will be renormalized, and other single particle terms will be induced or altered, these can include for example nonlinearities in the single layer dispersion in general, etc. .A change in v F , w AA , w AB parameters will tune the system away from the magic angle range.In this work, we have assumed that such change could be compensated by a change in the twist angle so that the system is tuned back to the new magic value for the twist angle as the interaction strength is altered.We have furthermore assumed that other induced effect (such as the monolayer nonlinear dispersion) could also be corrected by some means or are negligible and do not result in an appreciable effect.These assumptions allows us to also change the interaction strength in the terms correcting the zero point of our HF, and we will be left with the form in Eq. ( 7) with the interaction strength altered.
• Symmetries: -The C 2 T , when present, acts on the parallel transport basis as stated in the main text, transforms one band to the other with k unchanged: Note that y is an integer times a M 2 and we have used the translational properties shown in Fig. 1(a).-The particle hole symmetry exchanges the two bands of the HWFs basis, taking k x to −k x .In the parallel transport basis the states can be related by this transformation as follows: r, στ ψk,m,ξ = (−1) m i e −i(2ξq h )•r (−1) 1+τ (−x, y), στ ψ(−kx,ky), m,ξ . (C10) The factor (−1) m i in the above equation can be derived in the C 2 T symmetric case explicitly; it furthermore can be maintained in the ∆ = 0 case as well by appropriate phase redefinitions.The above property, furthermore, results in the symmetry of WCC positions under k x ↔ −k x , as seen in Fig. 1(b).
-The time reversal symmetry also relates the HWF states in the valleys in the following fashion: the HWF basis; it is straightforward to work out the sublattice potential form as well, since HWF basis is sublattice polarized: We also note that φ k,ky,2,ξ = −φ k,ky,1,ξ , regardless of the value of ∆.Now, it is straightforward to check that the terms in H kin , including the ∆ term, have the same form in terms of d operators as that in terms of c operators, where they are defined as in the following particle hole transformations: This means that H kin is particle hole symmetric with the above prescription.Spin indices could be trivially added to the above terms.
One should furthermore consider the interaction term; the interaction in general can be written as follows: The interaction thus takes the following form in terms of the d operators: where for simplicity a change of notation has been made 1 ≡ (k 1 , y 1 , ξ, s, m), and so forth.The first term above is identical in form to the original interaction Hamiltonian in terms of c operators.However, there are terms quadratic in d, the single hole terms, that were not present in the original Hamiltonian.Note that the terms on the third row are constant.The terms on the second row, on the other hand, turn out to be k x -dependent and thus impose a single hole potential; this is the origin of the particle hole asymmetry between the fillings ±3, as was discussed in the main text.Note that unlike this situation, in the usual Hubbard model with a single band, nearest neighbor hopping and constant on-site interaction for example, the analogue of this term is just a redefinition of the chemical potential.
Interestingly, the U (4) × U (4) symmetry of the chiral model (η = 0) discussed in this work can also be seen in the HWF basis (and also the parallel transport basis) as discussed in the present paper; for general η, when C 2 T is present, an interaction-only model consisting of active bands only displays a U (4) × U (4) symmetry (see Appendix C).Finally, we consider Ref. 14, where a HF study taking all bands into account has been presented.The authors consider several filling factors, and in particular, they are able to see a QAHE at ν = ±3; the presence of a significant sublattice potential is crucial for the QAHE to materialize.This is in contrast to the present work, where the presence of a sublattice potential can make the QAHE stronger, but it is not necessary for the occurrence of the required flavor polarization.Within our study, we observed that a larger interaction strength could compensate for the absence of the sublattice potential.

FIG. 1 . 3 4
FIG. 1.(a) The moiré lattice in real space and the corresponding BZ.A rectangular BZ is chosen − √ 3 4 < kx < √ 3 t e x i t s h a 1 _ b a s e 6 4 = " i Q x O g e a M 7 N 0

FIG. 4 .
FIG. 4.The HF results obtained within the projected model, i.e. the model used in our second study.(a) The converged HF energies normalized by g int vs kx.In all these four subplots, we take α = 0.58, ∆ = 0, ξ = 0.5a M .The interaction strength is chosen as g int = 0.05 and 0.002 for the left and right columns respectively.Plots in each column have the same parameter values except for the filling factor: the first row corresponds to ν = −3 and the second row to ν = +3.It can be seen that the HF energies at the two fillings ν = ±3 can be transformed into each other by the particle hole symmetry of the CM, which in particular needs kx → −kx.It can be seen that for η = 0, the two sets of HF energies are also related by the chiral symmetry of the CM.(b) The gaps as functions of g int .α = 0.58, ∆ = 0 are chosen for this plot and η and ξ are varied.