Exact Diagonalization for Magic-Angle Twisted Bilayer Graphene
EExact Diagonalization for Magic-Angle Twisted Bilayer Graphene
Pawel Potasz
Department of Physics, The University of Texas at Austin, Austin, TX 78712, USA andDepartment of Physics, Wroclaw University of Science and Technology, 50-370 Wroclaw, Poland
Ming Xie and A. H. MacDonald
Department of Physics, The University of Texas at Austin, Austin, TX 78712, USA (Dated: February 5, 2021)We report on finite-size exact-diagonalization calculations in a Hilbert space defined by thecontinuum-model flat moir´e bands of magic angle twisted bilayer graphene (MATBG). For moir´eband filling 3 > | ν | >
2, where superconductivity is strongest, we obtain evidence that the groundstate is a spin ferromagnet. Near | ν | = 3, we find Chern insulator ground states that have sponta-neous spin, valley, and sublattice polarization, and demonstrate that the anisotropy energy in thisorder-parameter space is strongly band-filling-factor dependent. We emphasize that inclusion of theremote band self-energy is necessary for a reliable description of MATBG flat band correlations, andcomment on possible mechanisms for superconductivity. Introduction: — Near a magic twist angle, the widthof bilayer graphene’s low energy moir´e bands shrinks[1, 2] by an order of magnitude or more, allowing in-teractions to play a prominent role in shaping electronicproperties. The flat bands form an octet that is the di-rect product of two-fold spin, valley, and band or sub-lattice degrees of freedom. The recent discovery of su-perconductivity and interaction-induced Chern and triv-ial insulator states [3–21] in magic-angle twisted bilayergraphene (MATBG) has motivated a rich body of the-oretical work [22–71]. It is already clear that althoughMATBG states share properties with doped and undopedMott insulators in conventional crystals, they also havea relationship to integer and fractional quantum Hall(FQH) states that flows from flat-band topological prop-erties [41, 57, 58]. The MATBG octet is closely analogousto the spin/valley/layer octet of Bernal bilayer graphene[72–75], but its orbital degree of freedom is more com-plex.In this Letter we report on exact-diagonalization cal-culations performed in a Hilbert space defined by con-tinuum model [2] moir´e bands. Our calculations confirm[15–17, 20] the spin, valley, and sublattice polarization ofthe | ν | = 3 insulating ground state, demonstrate that theanisotropy energy associated with this order is stronglyfilling factor dependent, and provide evidence for spin-polarized ground states for | ν | ∈ (2 , | ν | ∈ (2 , | ν | = 3. Panel(b) shows that the ground at | ν | = 3 is a spin and valleypolarized doublet formed by states with opposite sensesof spontaneous sublattice polarization. These states areknown to be Chern insulators and are accurately ap-proximated by Hartree-Fock theory. The ground state ofthe system with one charge added to (or removed from)the | ν = 3 | ground state (panel (c)) is still fully spin-polarized, but completely loses its valley polarization. As shown in panel (d) these states nevertheless have pre-cisely integer occupation numbers for all momenta, butonly when summed over valleys. The states with addedand removed charge have easy-plane valley order; we at-tribute the sudden change in anisotropy to the strongband/sublattice dependence of the single-particle Hamil-tonian at momenta near the γ -point in the moir´e Bril-louin zone. The sublattice polarization-properties (panel(e)) of the ground states near | ν | = 3, discussed furtherbelow, are revealed by the responses to sublattice andvalley dependent potentials illustrated in panels (b,c). Flat Band Projected Exact Diagonalization: —Becauseof large Dirac velocities, the electronic density-of-statesof an isolated neutral graphene sheet has a minimum atneutrality and is small over a broad energy range, al-lowing interaction effects to be described perturbatively.When magic-angle moir´e bands [2] are formed, strongelectronic correlations emerge and perturbative analy-ses are less reliable. Exact-diagonalization (ED) of theHamiltonian is a powerful non-perturbative method tostudy strong correlations, but, because the many-bodyHilbert space grows exponentially with system size, ispractical only when the single-particle Hamiltonian canbe truncated to a reasonably small dimension, typicallywith at most several tens of single-particle states. InMATBG the spectral isolation of the eight flat bands ofinterest (flat conduction and valence bands for each offour spin/valley flavors) motivates projection to an oc-cupation number subspace in which all remote valencebands in graphene’s negative-energy sea are fully occu-pied, all remote conduction bands are empty, and occu-pation numbers are allowed to fluctuate only within theflat bands. This strategy leads to a low-energy effectiveHamiltonian that acts entirely in the flat-band Hilbert- a r X i v : . [ c ond - m a t . s t r- e l ] F e b -0.2 -0.1 0.0 0.1 0.2-172.2-172.0-171.8-171.6-171.4 E [ m e V ] m [meV] -0.2 -0.1 0.0 0.1 0.2-159.5-159.0-158.5-158.0 E [ m e V ] m [meV] N el =9N el =10 (1) (1) (2) 𝑚 = −𝑚 ′ 𝑚 = 𝑚 ′ (2) (a)(b)(c)(e) (d) K K’ B A K K’ B A BA BA n E [ m e V ] N el D E D E (1)(2) K(1)(2) BAK BA FIG. 1: Spin, valley, and sublattice order vs. electron num-ber N el in finite-size MATBG with M = 9 moir´e unit cells:(a) ∆ E = E min ( S max ) − E min ( S < S max ) (meV/unit cell),where S = S K + S K is total spin and S max = N el / E = E min ( P v = 0) − E min ( P v = 1) where P v = | N K − N K | / ( N K + N K ). (b) Response of the groundstate energy of the valley polarized state at | ν | = 3 to anexternal field that couples to sublattice polarization. (c) Re-sponse of low energy states to a valley-odd sublattice field at N el = 10. The right panels in (b) and (c) schematically il-lustrate the sublattice polarizations induced in states (1) and(2) by the corresponding fields. (d) Ground state momentumspace occupation numbers projected to valley K (left), andtraced over valley (right) at N = 10. (e) Schematic illustra-tion of N el = 10 easy plane ferromagnet states. The groundstate (1)(left) has the same sense of sublattice polarizationin the two valleys, and at this system size, strong hybridiza-tion between the two sublattice polarization states (AA andBB). The first excited state (2) (right) has opposite senseof sublattice polarization in opposite valleys (AB) and (BA),and much weaker hybridization between the two degeneratestates at this system size. The spectra show evidence of reso-nant tunneling between ground and exited states confined tothe AB and BA valleys. space: H eff = X i ,i [ (cid:15) i δ i ,i + Σ i ,i ] c † i c i + 12 X i ,i,j ,j h i , j | V | i, j i c † i c † j c j c i , (1)where h i , j | V | i, j i is a two-body Coulomb interactionmatrix element, i , i, j , j label flat band states, (cid:15) i is aneigenvalue of the single-particle twisted-bilayer graphene Hamiltonian [2] including the interlayer tunneling contri-bution that is responsible for flat band formation, andΣ i ,i = X v [ h i , v | V | i, v i − h i , v | V | v, i i ] − X ¯ v [ h i , ¯ v | V | i, ¯ v i − h i , ¯ v | V | ¯ v, i i ] , (2)which we refer to the remote band self energy, accountsfor Hartree and exchange interactions with states v inthe frozen negative energy sea. In Eq. 2 the sum over¯ v in the regularization term is over the frozen valencebands of a neutral bilayer with no-interlayer tunneling[39]. As we shall emphasize, the remote band self-energyplays an essential role in MATBG physics and cannot beneglected.Because the bilayer has two honeycomb sublatticesper layer, the flat-band single-particle states are four-component sublattice spinors labelled by spin, valley, andmomentum in the moir´e Brillouin-zone. The two-particlematrix elements [39] are spin-independent, conserve val-ley at each vertex, and conserve crystal momentum. Incombination with the band dispersion (cid:15) i they define amany-body Hamiltonian that acts in a flat-band Hilbertspace with 2 × × × M single-particle states. ( M is thenumber of moir´e unit cells contained in the finite areaover which periodic boundary conditions are applied.)The three factors of 2 account for the four spin/valleyflavors and for the property that the flat bands occur asa conduction/valence pair, and together limit the largestvalue of M that can be considered.Fig. 2 plots the empty-system valence and conductionbands obtained by diagonalizing (cid:15) i δ i ,i + Σ i ,i for finitesize systems containing (a) M = 9 and (b) M = 36 moir´eunit cells at twist angle θ = 1 .
1, interaction strengthparameter ε − = 0 .
05, and inter-layer hopping parame-ters ω = 110 meV amd ω = 0 . ω for inter-sublatticeand intrasublattice tunneling. (Each band has a two-foldspin-degeneracy and two valley partners ( K K ) relatedby time-reversal symmetry. See [76] and Ref. [39] for amore detailed discussion of band and screening ( (cid:15) ) pa-rameters.) We see in Fig. 2 that the main features ofthe flat bands are captured already at M = 9, used formost of our ED calculations. The valence band is ex-tremely flat with a width less than 5 meV and a localmaximum at the moir´e Brillouin-zone γ -point, The con-duction band is much broader with a width larger than25 meV and a sharp peak at γ . Importantly the energyseparation between conduction and valence bands at γ iscomparable to the largest energy scales in the problem.The conduction band width greatly exceeds the few meVwidth [33] obtained when the remote band self-energy isneglected, but does not strongly influence the band sepa-ration. The bands have a three-fold rotational symmetrywhich is more clearly revealed for M = 36.The remote band self-energy reshapes the bands prin-cipally by shifting energies near γ upward, relative to -15.01.518.0-13.03.019.0 valence -16.0-14.0-12.0 (a) -15.0-13.5-12.0 (b) γ κκ' E [meV] E [meV]E [meV] E [meV]
FIG. 2: Quasiparticle energies of empty flat bands for (a) M = 9 ( N x = 3, N y = 3) and (b) M = 36 ( N x = 6, N y = 6)moir´e unit cell systems at twist angle θ = 1 .
1, interactionstrength parameter (cid:15) − = 0 .
05, and inter-layer intrasublat-tice and intersublattice hopping parameters ω = 110 meVand ω = 0 . ω . The valence band is on the left and theconduction band on the right. The color maps are interpola-tions from the discrete finite-size momenta indicated by theblack squares. Empty squares are used to mark equivalentmomenta. those near κ , κ . The relative shifts occur primarily be-cause the Hartree potential from the remote bands is at-tractive near the AA positions where states near γ haveless weight [38, 39]. The sharp contrast between the con-duction and valence band widths in these empty-banddispersions does not imply strong particle-hole asym-metry. Indeed the model we will study is very nearlyparticle-hole symmetric, and the relative widths of thebands is reversed when we describe flat band states interms of interacting holes instead of interacting electrons[76]. Instead, the upward shift at γ works in concertwith weaker electron-electron repulsion matrix elementsfor states near γ [76] that reduce their Coulomb energypenalty as the flat bands are filled. The ED results inthis MS were calculated for the M = 9 system, which issufficiently large to capture the important distinction be-tween states near γ and those in the rest of the Brillouin-zone.The many-body Hamiltonian separates into decoupledblocks labelled by the number of electrons in each valley N K and N K , valley-dependent total ( S K and S K ) andazimuthal spin ( S zK and S zK ) quantum numbers, and to-tal crystal momentum ( K x , K y ). The separate spin quan-tum numbers for the two valleys valley apply because themodel is invariant under independent valley-dependent spin-rotations.Unlike the model we study, experimental samples ex-hibit considerable particle-hole asymmetry. For exam-ple, the Chern insulator states we discuss below tend tobe more prominent at positive than at negative fillingfactors. The assymmetry is thought [40] to be due tonon-local corrections to the interlayer tunneling modelwe employ. The relationship of our findings to exper-iment is addressed more fully in the discussion sectionbelow. Numerical Results: — Our first important result is re-lated to the regime in which | ν | ∈ (2 , | ν | < | ν | ∈ (2 ,
3) interval, which correspondsto the N el ∈ [10 ,
18] in our flat-band projected ED cal-culation. For N el = 10 , ,
12 full Hilbert space calcu-lations confirm that the ground state is maximally spin-polarized, as illustrated in Fig. 3. For larger N el we canshow that the fully spin-polarized state is lower in en-ergy than the corresponding fully valley-polarized state.Some of these conclusions rest on extrapolations fromcalculations performed in a selected subspace of the fullED Hilbert space, as explained in the supplementary ma-terial [76]. The conclusion that the ground state is fullyspin-polarized helps constrain and simplify potential the-ories of superconductivity.For N el = 9 ( | ν | = 3), we are able to fully explorenearly all subspaces, including all with particles dis-tributed over three flavors, and subspaces with particlesdistributed over four flavors provided one of the flavors isfilled by at least five particles. We find that the groundstate is fully spin and valley polarized, and well approx-imated by a single Slater determinant. For example, wefind that the maximum deviation from unit momentum-state occupation across the Brillouin-zone is 0.04. Theground state appears as a quartet with nearly degener-ate doublets for each sense of valley polarization. Bystudying the response of this doublet to a sublattice de-pendent potential m σ z , where σ z acts on the sublatticedegree-of-freedom in both layers, we see that for a givenvalley polarization the doublet is formed by states withopposite sublattice polarizations and that there is observ-able hybridization between these states. It is known fromHartree-Fock theory that these states are Chern insula-tors with Chern number magnitudes | C | = 1 and signsdetermined by the sign of the product of the valley andsublattice polarization. In Refs. [43, 44] the two stateswith the same Chern number are described by a σ -model electrons S K S K' S K , S K ' N el P v N el FIG. 3: Ground state spin and valley quantum numbers as afunction of electron number N el . Top: Total spin S K and S K in each valley. Bottom: Valley polarization P v . Integer bandfilling ν = − N el = 9 highlighted by a dashed line. S = S K + S K is total spin. in which only the orientation of the corresponding pseu-dospin is retained as a relevant degree-of-freedom. TheChern insulator at | ν | = 3 [41, 42, 60, 61] can be viewedas a simple ferromagnet formed from these pseudospins.From Fig. 1(b) we conclude that h σ z i ∼ .
25, implyingthat P sub = h σ z i /N el ∼ .
25, in agreements with previ-ous Hartree-Fock results [42], and that the Hamiltonianmatrix element for collective tunneling between stateswith opposite senses of spontaneous sublattice polariza-tion (which is expected to fall exponentially with systemsize) is ∼ .
058 meV for N el = M = 9 and (based on aseparate calculation) ∼ . N el = M = 16calculation. Easy-Plane Valley Anisotropy: — In Fig. 1(c) andFig. 3 we see that valley polarization is completely lostwhen we add or remove one electron from the N el = 9valley and spin polarized ferromagnet. We attribute thisbehavior to the strong band splitting at γ , which has anout-sized influence on valley anisotropy by suppressingthe band-mixing degree of freedom. An important el-ement of our interpretation is the observation that oursystem has only U (1) and not SU (2) valley symmetry.In the language of magnetism our system has uniaxialvalley anisotropy, which allows easy axis or easy-planevalley magnetism. Our conclusion that the state at γ plays a crucial role is supported by the property that the excitation spectra at N = 8, where the γ -state is empty,and at N = 10, where the γ -state is doubly occupied, arenearly identical. Our calculations confirm that easy axissublattice/band order is present for both easy-axis andeasy-plane valley anisotropy, with four degenerate clas-sical states distinguished by the sublattice polarizationof K and K valley components in the easy plane case.In Fig. 1(c) the ground state responds most strongly tosublattice potentials that are identical in the two valleys,demonstrating AA or BB sublattice polarizations. Thesetwo classical states should have identical energies, andwe conclude from the ED spectra that the tunneling be-tween them is large at this system size. We associatethe excited state doublets in Fig. 1(c) with AB and BAsublattice polarizaiton for valleys KK’. This interpreta-tion is supported by strong response to valley-odd sub-lattice potentials. Our ED results demonstrate that themany-particle tunneling matrix element between thesesublattice states is greatly reduced compared to tunnel-ing between degenerate AA and BB states. In this casethe ED spectra exhibit resonant tunneling non-only be-tween ground states, but also between excited states ofthe isolated AB and BA sectors. Discussion: — Our calculations show that MATBGground state energies are generally speaking well approxi-mated by unrestricted Hartree-Fock approximations thatallow spin, valley, and sublattice symmetries to be bro-ken. In the top panel of Fig. 4 we show the dependenceof the correlation energy on electron number in the sub-space with full spin and valley flavor polarization over thefull range of available filling factor for that flavor between ν f = − ν f = 1. The correlation energy, defined asthe difference between the ED ground state energy andthe minimum energy single-Slater determinant, vanisheswhen the orbital doublet is empty ( ν f = −
1) and full( ν f = 1), and also reaches an extremely small value inthe insulating state at ν f = 0. These results suggest in-sulating states at all integer filling factors with a bandfilling per flavor equal to -1, 0, or 1 are accurately ren-dered by Hartree-Fock, even when symmetry is broken bychoosing different band filling factors for different flavors.For a given total integer filling factor a variety of differentstates, characterized by different flavor-dependent fillingfactors and senses of sublattice polarization, are expectedto compete closly in energy. The states have different to-tal Chern numbers with | C | = 1 for | ν | = 3, 0 or 2 for | ν | = 2, 1 or 3 for | ν | = 1, and 0, 2 or 4 for ν = 0. (The | ν | = 1 and | ν | = 0 cases are outside of the reach of ED.)In Fig. 4 we see that the correlation energy is larger awayfrom integer filling factors. We expect that this trendwill be stronger in sectors with less flavor polarization,and that Hartree-Fock therefore overestimates the ten-dency to break flavor symmetries. Insulating states at | ν = 1 | , will therefore compete with metallic states withno broken flavor symmetries that have much larger cor-relation energies. The difference in energy between the E c o rr [ m e V ] q [] E corr P s ub P EDsub P HFsub
Correlation Energy
Filling factor N el =9 E [ m e V ] N el E Polcorr =E Poltot -E PolHF D E tot =E tot - E Poltot E c o rr [ m e V ] , P s ub q [] E corr FIG. 4: The correlation energy per moir´e period as a functionof filling factor (top) and as a function of the twist angle for N el = 9 together with corresponding sublattice polarizations(bottom). Within this range of angles, both exact diagonal-ization and Hartree-Fock calculations predict full flavor po-larization. Top: The green triangles correspond to one flavorfull-spin-polarization calculations and E corr = E Poltot − E PolHF ,where E Poltot is the ground state energy from exact diagonal-ization calculations and E PolHF is the mean-field energy fromself-consistent Hartree-Fock calculations. The black squaresshow to the energy difference between the ground state en-ergy and the lowest energy state in the full flavor polarizationsector. The black dashed line marks the filling factor ν = − N el = 9 with greentriangles for the correlation energy and black squares red cir-cles for ED and HF sublattice polarizations. true ground state and the ground state in the fully po-larized sector increases quickly for N el >
9, showing thatthe energy cost of valley polarization quickly increases.The appearance of insulating states at integer fillingfactors depends on screening environment, twist angle,and band-structure details that we have not fully ex-plored here. For example in the bottom panel of Fig. 4we illustrate how the correlation energy of the | ν | = − θ ,the screening parameter (cid:15) used in our calculations in- fluences quantitative conclusions. Interactions in theflat bands of MATBG are screened by the surround-ing hexagonal Boron Nitride (hBN) dielectric, by thenearby electrical gates, and by transitions between theflat and remote bands. Strictly speaking, the latter twoeffects yield wavevector-dependent contributions to thedielectric constant with gates dominating in importanceat small wavevectors and screening within the bilayerdominating a larger wavevectors. Quantiative propertyprediction for particular samples will require accuratelyknown twist angles and geometries, and band structuremodel that account for both non-local interlayer tunnel-ing and, when appropriate, the influence of the adjacentboron nitride. Acknowledgment:
This work was supported by DOEBES under Award DE-FG02-02ER45958. PP acknowl-edges financial support by the Polish National Agencyfor Academic Exchange (NAWA). The authors acknowl-edge helpful interactions with A. Bernevig, N. Regnault,T. Senthil, A. Vishwanath, and F. Wu. ED calculationswere performed using the Wroc law Center for Networkingand Supercomputing and the Texas Advanced Comput-ing Center (TACC). [1] E. Suarez Morell, J. D. Correa, P. Vargas, M. Pachecoand Z. Barticevic, Phys. Rev. B , 121407 (2010).[2] R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad.Sci. U.S.A. , 12233 (2011).[3] K. Kim, A. DaSilva, S. Huang, B. Fallahazad, S. Lar-entis, T. Taniguchi, K. Watanabe, B. J. LeRoy, A. H.MacDonald, and E. Tutuc, Proc. Natl. Acad. Sci. U.S.A. , 3364 (2017).[4] 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, Nature , 80 (2018).[5] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,E. Kaxiras, and P. Jarillo-Herrero, Nature , 43(2018).[6] M. Yankowitz, S. Chen, H. Polshyn, K. Watanabe, T.Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science , 11059 (2019).[7] X. Lu, P. Stepanov, W. Yang, M. Xie, M. Ali Aamir, I.Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang,A. Bachtold, A. H. MacDonald and D. K. Efetov, Nature , 656 (2019).[8] A. Uri, S. Grover, Y. Cao, J. A. Crosse, K. Bagani,D. Rodan-Legrain, Y. Myasoedov, K. Watanabe, T.Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, andE. Zeldov, Nature , 47 (2020).[9] L. Balents, C. R. Dean, D. K. Efetov and A. F. Young,Nature Physics , 725–733 (2020).[10] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S.Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig andA. Yazdani, Nature , 198–202 (2020).[11] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R.Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, Ady Stern, E. Berg, P. Jarillo-Herrero and S.Ilani Nature , 203–208 (2020).[12] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney,K. Watanabe, T. Taniguchi, M. A. Kastner, and D.Goldhaber-Gordon, Science , 605 (2019).[13] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J.Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F.Young, Science , 900 (2020).[14] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe,T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov,and D. K. Efetov, Nature 583, 375 (2020).[15] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y.Andrei, arXiv:2007.03735[16] I. Das, X. Lu, J. Herzog-Arbeitman, Z.-D. Song, K.Watanabe, T. Taniguchi, B. A. Bernevig, and D. K. Efe-tov, arXiv:2007.13390[17] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, P.Jarillo-Herrero arXiv:2008.12296v1[18] H. S. Arora, R. Polski, Y. Zhang, A. Thomson, Y.Choi, H. Kim, Z. Lin, I. Z.Wilson, X. Xu, J.-H. Chu,K.Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge,Nature , 379 (2020).[19] Y. Cao, D. Rodan-Legrain, J. Min Park, F. N. Yuan,K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, P.Jarillo-Herrero, arXiv:2004.04148[20] P. Stepanov, M. Xie, T. Taniguchi, K. Watanabe, X. Lu,A. H. MacDonald, B. Andrei Bernevig, and D. K. Efetov,arXiv:2012.15126[21] A. Rozen, J. M. Park, U. Zondiner, Y. Cao, D.Rodan-Legrain, T. Taniguchi, K. Watanabe, Y. Oreg,Ady Stern, E. Berg, P. Jarillo-Herrero and S. Ilaniarxiv:2009.01836v2[22] N. F. Q. Yuan and L. Fu, Phys. Rev. B , 045103 (2018)[23] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Phys.Rev. X , 031089 (2018).[24] H. Guo, X. Zhu, S. Feng, and R. T. Scalettar, Phys. Rev.B , 235453 (2018).[25] B. Padhi, C. Setty, and P. W. Phillips, Nano Lett. ,6175 (2018).[26] V. Y. Irkhin and Y. N. Skryabin, JETP Lett. , 651(2018).[27] J. F. Dodaro, S. A. Kivelson, Y. Schattner, X.-Q. Sun,and C. Wang, Phys. Rev. B , 075154 (2018).[28] T. Huang, L. Zhang, and T. Ma, arXiv:1804.06096.[29] C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Phys.Rev. Lett. , 217001 (2018)[30] X. Y. Xu, K. T. Law, and P. A. Lee, Phys. Rev. B ,121406(R) (2018).[31] J. Kang and O. Vafek, Phys. Rev. X , 031088 (2018).[32] L. Rademaker and P. Mellado, arXiv:1805.05294.[33] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K.Kuroki, and L. Fu, Phys. Rev. X , 031087 (2018).[34] J. M. Pizarro, M. J. Calderon, and E. Bascones,arXiv:1805.07303.[35] M. Ochi, M. Koshino, and K. Kuroki, Phys. Rev. B ,081102 (2018).[36] H. Isobe, N. F. Q. Yuan, and L. Fu, Phys. Rev. X ,041041 (2018).[37] A. Thomson, S. Chatterjee, S. Sachdev, and M. S.Scheurer, Phys. Rev. B , 075109 (2018).[38] F. Guinea and N. R. Walet, Proc. Natl. Acad. Sci. U.S.A. , 13174 (2018).[39] Ming Xie, Allan H. MacDonald Phys. Rev. Lett. ,097601 (2020). [40] Ming Xie, Allan H. MacDonald arXiv:2010.07928v1[41] N. Bultinck, S. Chatterjee, and M. P. Zaletel, Phys. Rev.Lett. , 166601 (2020)[42] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vish-wanath,and M. P. Zaletel, Phys. Rev. X , 031034(2020)[43] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, andA. Vishwanath, arXiv:2004.00638v2[44] E. Khalaf, N. Bultinck, A. Vishwanath, M. P. ZaletelarXiv:2009.14827[45] S. Chatterjee, N. Bultinck, and M. P. Zaletel, Phys. Rev.B , 165141 (2020)[46] S. Chatterjee, M. Ippoliti, and M. P. Zaletel,arXiv:2010.01144[47] D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel,and Nick Bultinck, arXiv:2012.09885v1[48] C. Xu and L. Balents, Phys. Rev. Lett. , 087001(2018).[49] F. Wu, A. H. MacDonald, and I. Martin,arXiv:1805.08735.[50] B. Roy and V. Juricic, arXiv:1803.11190.[51] S. Ray and T. Das, arXiv:1804.09674.[52] T. J. Peltonen, R. Ojajarvi, and T. T. Heikkila,arXiv:1805.01039.[53] Y.-Z. You and A. Vishwanath, arXiv:1805.06867.[54] X.-C. Wu, , K. A. Pawlak, C.-M. Jian, and C. Xu,arXiv:1805.06906.[55] M. Fidrysiak, M. Zegrodnik, and J. Spalek, Phys. Rev.B , 085436 (2018).[56] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H.Castro Neto, Phys. Rev. Lett. , 256802 (2007).[57] H. C. Po, H. Watanabe, and A. Vishwanath, Phys. Rev.Lett. , 126402 (2018).[58] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Phys.Rev. B , 195455 (2019).[59] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath,Phys. Rev. Lett. , 106405 (2019).[60] J. Liu, J. Liu, and X. Dai, Phys. Rev. B , 155415(2019).[61] Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T.Senthil, Phys. Rev. B , 075127 (2019).[62] G. W. Semenoff, Phys. Scr. , 014016 (2011).[63] J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Phys.Rev. B , 205414 (2014).[64] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F.Young, arXiv:1911.13302v1[65] A. Kumar, M. Xie, and A. H. MacDonaldarXiv:2010.05946v1[66] B. Andrei Bernevig, Z. Song, N. Regnault, and B. Lian,arXiv:2009.11301v1[67] Z.-D. Song, B. Lian, N. Regnault, and A. B. Bernevig,arXiv:2009.11872v1[68] B. Andrei Bernevig, Z. Song, N. Regnault, and B. Lian,arXiv:2009.12376v1[69] B. Lian, Z. Song, N. Regnault, D. K. Efetov, A. Yazdani,and B. A. Bernevig, arXiv:2009.13530v1[70] B. Andrei Bernevig, B. Lian, A. Cowsik, F. Xie, N. Reg-nault, and Z. Song, arXiv:2009.14200v1[71] F. Xie, A. Cowsik, Z. Song, B. Lian, B. Andrei Bernevig,and N. Regnault, arXiv:2010.00588v2[72] Maxim Kharitonov, Phys. Rev. B , 195435 (2012).[73] F. Zhang, J. Jung, G. A. Fiete, Q. Niu, and A. H. Mac-Donald, Phys. Rev. Lett. , 156801 (2011).[74] A. H. MacDonald1, J. Jung, and F. Zhang, Phys. Scr. , 014012 (2012).[75] Y. Barlas, R. Cˆot´e, K. Nomura, A. H. MacDonald, Phys.Rev. Lett. , 097601 (2008). [76] Supplementary materials upplementary MaterialsExact Diagonalization for Magic Angle Twisted Bilayer Graphene Pawel Potasz,
1, 2
Ming Xie, and A. H. MacDonald Department of Physics, The University of Texas at Austin, Austin, TX 78712, USA Department of Physics, Wroclaw University of Science and Technology, 50-370 Wroclaw, Poland
COULOMB MATRIX ELEMENTS
Exact diagonalization (ED) calculations account for alldirect and exchange terms, and also for the many val-ley and momentum conserving processes in which pairsof electrons are scattered between different occupationnumber states. In Fig. S1 direct (top) and exchange(bottom) Coulomb matrix elements are shown as a func-tion of their single-particle state momentum and orbitallabels. We label the single-particle states by i = k y N x + k x for the valence band and i = N x ∗ N y + k y N x + k x for theconduction band, where k = k x N x b + k y N y b , b , are recip-rocal lattice vectors, and N x and N y are the integers thatspecify the numbers of unit cells in the finite size systemin the primitive lattice vector directions ( M = N x × N y is a total number of unit cells) in the finite-size system.Here we use the terms valence and conduction bands torefer respectively to the lower and higher energy statesobtained by diagonalizaing the single-particle Hamilto-nian including the remote band self-energy contribution.The state indices i = 0 and i = 9 correspond to the γ point. Clearly direct interactions with these states areweaker than all other Coulomb matrix elements. This re-flects the fact that wavefuncions around the γ point aremore spread in a real space comparing to k -points fromother regions of the Brillouin zone. In the map of ex-change interactions one can distinguish stronger interac-tion between states with momentum k and − k indicatedby a blue color with values above 2 . SUBLATTICE POLARIZATIONS AND BANDOCCUPATIONS
A bilayer has four sublattices α = { A , B , A , B } ,where A, B distinguish the two sublattices of of onegraphene honeycomb layer and 1 , n α ( r ) = h Ψ + α ( r )Ψ α ( r ) i , (1)where the field operators areΨ + α ( r ) = 1 √ M X k ,n φ ∗ k ,n,α ( r ) a + k ,n,α , (2) s t a t e i nd e x state index s t a t e i nd e x state index exchange FIG. S1: A map of direct Coulomb (top) and exchange (bot-tom) matrix elements within the flat bands. State indices h , i correspond to the valence band with a state index i , i = k y N x + k x , and h , i to the conduction band with astate index i , i = N x ∗ N y + k y N x + k x . The γ points cor-respond to the state index i = 0 for the valence band and i = 9 for the conduction band. Direct Coulomb matrix ele-ments from i = 0 and i = 9 to all the states have clearly lowervalues in comparison to all other matrix elements. and a + k ,n,α are creation operators of particles with mo-mentum k and band index n , and φ ∗ k ,n ( r ) are Blochwavefunctions expressed in a plane wave basis of moir´ereciprocal lattice vectors G : φ k ,n,α ( r ) = X G z nα, G ( k ) e i ( k + G ) r , (3)where z nα, G ( k ) are the expansion coefficients of the con-tinuum model spinors.The many-body state is defined as | ψ MBl i = P s A ls | s i a r X i v : . [ c ond - m a t . s t r- e l ] F e b with coefficients A ls obtained from diagonalization Hamil-tonian matrix, where | s i = Q i s a + k is | i are occupationnumber configurations corresponding to Slater determi-nant (SD)s wavefunctions. Using the definitions above,the density per moir´e cell on sublattice α is n MBα = 1 M X n s ,n p X k , G z n s G ,α ( k ) z n p G ,α ( k ) A ∗ ls A lp ( − r (4)where z n s are coefficients of an occupied band n s of s -thbasis configuration, and a factor ( − r comes from anti-commutation of operators and r is a number of exchangesfor a given set ( n s , n p ). The sublattice polarization is de-fined as P sub = | n MBA + n MBA − n MBB − n MBB | P α n MBα . (5)We can also define a total valence band occupation n V B ,by restricting indices n s to ones from the valence bandonly in Eq. 4 and summing over all sublattices α . Withthese definitions, P sub and n V B are quantities per parti-cle, thus normalized to one in the case of full sublatticepolarization or full valence band occupation.All sublattices are equivalent and an external field isnecessary to break the symmetry between them. In or-der to induce a finite sublattice polarization, a sublatticemass term is added, described by the HamiltonianH m = X k ,n,α,µ m µ σ αα z a + k ,n,α a k ,n,α , (6)where momentum k is restricted to a single valley µ = K, K , σ z is a Pauli matrix in the z direction in spaceof sublattices and we consider two cases, one with samemass m = m K = m K and one with opposite mass terms m = m K = − m K in the two valleys.Fig. S2 shows the total valence band occupation n V B (black squares) and the sublattice polarization of theground state as a function of a particle number N el .(Here the valence band is defined as the lower band of flatband states when the remote band self-energy is includedin the single-particle Hamiltonian.) We evaluate sub-lattice polarization P sub for the ED many-body groundstate (red circles) and for the case in which only one fla-vor is occupied (green triangles), both calculated for afinite sublattice mass m = 1 meV. For a small filling,the mixing between the valence and conduction bands isweak. For N el < sub , whichis close to zero, confirming approximately equal filling ofall aublattices, and four flavors. Moving closer to integerfilling ν = − N el = 9, the valence band occupa-tion fraction drops below n V B = 0 . N el = 7,with a tiny increase of sublattice polarization, and spin VB P sub P polsub n V B , P s ub N el Band mixing, valence band occupation
FIG. S2: The valence band occupation n V B and sublatticepolarization P sub of the many-body ground state obtain fromexact diagonalization as a function of the filling factor for asublattice mass m = 1 . P polsub for a flavor polarized lowest energy state isalso shown. A vertical dash line indicates an integer fillingfactor ν = − N el = 9 electrons. polarization (see MS). Interactions start playing the cen-tral role as the integer filling factor is approached. At N el = 9 the ground state has clear spontaneous sublat-tice polarization as emphasized in the main MS. This isseen in these relatively large m calculations as a peakin P sub . Similar behavior is observed for the single fla-vor case (green triangles), except at N el = 8, where thepeak in the sublattice polarization is the highest. Thisis related to a formation of a state with a uniform fill-ing of all momentum points except the one at γ , whichremains nearly empty. This behavior is a reminder thatthe interactions of particles filling in γ with other elec-trons are significantly weaker compared to interactionsbetween particles filling other momenta. Further studieswill be needed to determine how universal these detailsare when the twist angle and the band structure modelare changed, however we anticipate that the details ofsublattice and valley order are sensitive to changes inband filling factor that change occupation numbers near γ , where both the wavefunctions and the dispersion areanomalous. Above integer filling ν = −
3, sublattice po-larizations deceases, and the fractional valence band oc-cupation increases.
FILLING ν = − The insulating state for N el = 9 is a spontaneous sub-lattice symmetry broken state. In Fig. S3 we showthe total flat band occupation (left) and the occupa-tion projected onto valence band states (right). Thedeviation from equal occupation of each momentum isless than 0.04. However the valence band occupationchanges strongly with momentum, indicating that theground state consists of k-dependent nontrivial mixturesof the original valence and conduction band states. In thevalence band occupation numbers, one can also notice asymmetry with respect to the k y = 0 axis.We analyze the response of the four lowest many-bodyenergy states to a sublattice potential mass term in thetop panel of Fig. S4. The ground state is nearly doubledegenerate, with a splitting related to collective tunnel-ing between states with opposite senses of spontaneoussublattice polarization that should vanish in a thermody-namic limit. These two lowest energy states reveal lineardependence of energy on mass, while two upper statesremain relatively unaffected by the sublattice field, asalso shown in the main MS. We show also the corre-sponding HF energy dependence on m , obtained for thesame M = 9 system size. The comparison allows us todetermine the correlation energy E corr , indicated by avertical red arrow, which remains almost constant with m . The magnitude of the sublattice polarization can bedetermined by P sub = ( dEdm ) /N el (middle panel). Thecurve corresponding to the lowest state (black squares)quickly saturates at P EDsub ∼ .
25, while the curve forthe second state (green triangles) reaches P sub ∼ . m = 0 .
02 and drops down with further increases in m due to level-repulsion by the third state. The HF resultshave slightly larger P HFsub ∼ . P sub = h σ z i (bottom panel). Sublattice polarizationsfrom ED smoothly increases with m , approaching theHF result for sufficiently large m . The blue diamondsillustrate the single particle (SP) sublattice polarization.In summary, from both definitions of sublattice polariza-tion, and both methods, ED and HF, sublattice polar-ization is P sub ∼ .
25, and close to a value in a ther-modynamic limit obtained from HF (inset in the bottompanel).Fig. S5 shows the low energy many-body spectrumfor N el = 9 as a function of valley polarization P v . Theground state is maximally flavor polarized, and energygradually increases with a decrease of valley polariza-tion. The energy gap seen between the four lowest energystates and higher energy states is an artifact of plottingonly the four lowest states from every total momentumsubspace. The energy gap for P v = 1 between the lowestenergetic doublet, below E = − . E = − . P v = 0 .
33 do K↓ Full occupation VB occupation
FIG. S3: Total flat band occupation (left) and valence bandoccupation (right) for the N el = 9 ground state. The filledsquares correspond to the discrete k -points of the M = 9finite-size system and empty squares to the correspondingpoints in neighboring Brillouin zones. not exhibit a linear response to mass fields, while the en-ergy states around E = 158 . N K = 6, N K = 3). ONE PARTICLE CHARGE EXCITATION FROMFILLING ν = − We investigate charge excitation from the ν = − N el = 8, or adding, N el = 10, one par-ticle. In Fig. S6 many-body spectra as a function ofvalley polarization are shown. Both types of charge ex-citation have fully spin polarized and fully valley unpo-larized ground states. For N el = 8, the energy of low-est excitations gradually increases with valley polariza-tion, and a quite similar spectrum is seen for N el = 10,with the largest energetic cost of the lowest excitationfor full valley polarization. Similarity between ± N el = 10, the two extra particles added to the sys-tem with N el = 8, fill γ points from the two valleys,and occupations of other momenta remain almost un-changed. This is consistent with the Coulomb matrixelements, where the particles filling the γ point interactvery weakly with particles at other momenta. Total oc-cupation summed over two valleys is close to one for ± ∼ .
001 for N el = 8 and ∼ .
011 for N el = 10. The occupation indicated by the purple colorfor N el = 10 is close to two, ∼ .
94. In Fig. S8 lowenergy sublattice mass dependence for (a) N el = 8 and(b) N el = 10 for the same (left) and opposite (right)sublattice mass fields in two valleys are shown. The low-est state slightly response to the even mass field and iscompletely not affected by the odd mass field. This isan evidence of the character of this state, with equal oc- -0.2 -0.1 0.0 0.1 0.2-0.3-0.2-0.10.00.10.20.3 P s ub m [meV] ED state (1) ED state (2) HF-0.2 -0.1 0.0 0.1 0.2-159.5-159.0-158.5-158.0 E [ m e V ] m [meV] N el =9 (1)(2) ED vs HF HF P s ub m [meV] P EDsub P HFsub P SPsub
Derivative E corr P s ub FIG. S4: Top panel: Total energy as a function of a thesublattice mass term for four low-energy many-body statesobtained from exact diagonalization (ED) (black squares) andHartree-Fock (HF) energy (red circles). The vertical arrowshows the correlation energy E corr . Middle panel: Sublatticepolarization calculated using the relation P sub = ( dEdm ) /N el .Bottom panel: Sublattice polarization P sub as a function ofthe sublattice mass term for ED, HF, and single particle (SP)ground states calculated using the definition P sub = h σ z i .The inset show extrapolation of HF sublattice polarization toa thermodynamic limit. cupation of the opposite sublattice in two valleys, andslight imbalance between two sublattice within a valley.Three doublets of higher energy states show a clear pat-tern in the odd mass field, a linear response, and energyseparation between them that increases approximatelyquadratically, and the splitting within a doublet decreas-ing, with a doublet index. This suggest that these aresets of bound states separated by a finite size tunnellingbarrier. We obtain sublattice polarization of the lowestdoublet from P sub = ( dEdm ) /N el , for the opposite mass E [ m e V ] P v E [ m e V ] K x *N y +K y P v N el =10 E [ m e V ] P v N el =9 FIG. S5: Low energy many-body spectrum as a function ofvalley polarization P v for N el = 9. E [ m e V ] P v E [ m e V ] P v E [ m e V ] K x *N y +K y P v N el =10N el =8 FIG. S6: Low energy many-body spectrum as a function ofvalley polarization P v for N el = 8 (left) and N el = 10 (right). field P sub ( N el = 8) ∼ .
077 and P sub ( N el = 10) ∼ . P sub ( N el = 9) ∼ . FILLING FACTORS BETWEEN ν = +3 AND ν = +4 We perform particle-hole transformation and investi-gate electronic properties in a vicinity of the full flatbands. The quasiparticle band dispersions for a singlehole are shown in Fig. S9 for the same set of parametersas in the case of empty bands considered in this work.The band dispersion is very similar to the one for elec-trons confirming particle-hole symmetry in the system.We depopulate these bands, or equivalently create holes,and analyze properties of the ground state between fill-ing factors ν = +3 and ν = +4. The integer filling factor ν = +3 is for a system with N h = 9 holes. The holesphase diagram is shown in Fig. S10. Similarities withthe electron phase diagram are clearly visible. For low K K’ K+K’K K’ K+K’ (a) (b)
FIG. S7: Total flat band occupation of the ground state for(a) N el = 8 and (b) N el = 10, for valley K (left), valley K (middle), and both valleys K + K (right). An occupationindicated by a purple color for N el = 10 is close to two, ∼ . M = 9 k -points from a discretizedfirst Brillouin zone and empty squares to corresponding pointsfrom neighboring Brillouin zones. -0.2 -0.1 0.0 0.1 0.2-172.2-172.0-171.8-171.6-171.4-171.2 E [ m e V ] m [meV] -0.2 -0.1 0.0 0.1 0.2-172.2-172.0-171.8-171.6-171.4-171.2 E [ m e V ] m [meV] -0.2 -0.1 0.0 0.1 0.2-145.0-144.5-144.0 E [ m e V ] m [meV] -0.2 -0.1 0.0 0.1 0.2-145.0-144.5-144.0 E [ m e V ] m [meV] Same mass Opposite mass (a)(b) N el =8 N el =8N el =10 N el =10 FIG. S8: A sublattice mass dependence of the (a) N el = 8and (b) N el = 10 four lowest energy states with the samemass field m = m K = m K (left) and opposite mass field m = m K = − m K (right). densities, no clear flavor polarization is observed witha transition to full spin polarization for N h >
6. Onecan notice no full valley depolarization for N h = 10, incontrast to the electron case. Summarizing, for a modelconsidered in this work, the system has the electron-holesymmetry up to a very good approximation. -18.0-1.515.0 -13.03.019.0 valence Momentum maps holes -19.5-17.5-15.5 (a) -15.0-13.5-12.0 (b) γ κκ' FIG. S9: Quasiparticle energies of full flat bands for M = 9( N x = 3, N y = 3) moir´e unit cell systems at twist angle θ = 1 .
1, interaction strength parameter (cid:15) − = 0 .
05, and inter-layer hopping parameters ω = 110 meV for inter sublatticetunneling and ω = 0 . ω for intrasublattice tunneling. Thevalence band is on the left and the conduction band on theright. The color maps are interpolations from discrete meshesindicated by the black squares; empty squares are used tomark momentum points that are equivalent to those includedin the discrete mesh.
10 9 8 7 6 5 4 3 2 1 0012345 S K S K' S K , S K ' N h
10 9 8 7 6 5 4 3 2 1 00.00.51.0 P v N h
10 9 8 7 6 5 4 3 2 1 001234510 9 8 7 6 5 4 3 2 1 00.00.10.20.310 9 8 7 6 5 4 3 2 1 0 S K S K' S K , S K ' N h E G a p [ m e V ] N h P v N h electrons holes FIG. S10: The ground state properties as a function of a fillingfactor for holes. Top: Total spin S K and S K in each valley.Bottom: Valley polarization P v . The integer filling ν = +3 isindicated by a dash line and corresponds to N h = 9. FLAT BANDS AT CHARGE NEUTRALITY
Fig. S11 shows the quasi-particle bands at chargeneutrality. In contrast to the empty flat band quasi-particles, the bottom of the valence band is at γ , as inthe non-interacting model, the bands have approximateparticle-hole symmetry, however the energy difference oftwo bands at γ is significantly larger, around 30 meV,in contrast to few meV without the self-energies effect ofremote bands. -1.47.215.7-13.03.019.0 valence Momentum maps CN -16.2-7.41.3 (a) -15.0-13.5-12.0 (b) γ κκ' FIG. S11: Quasiparticle energies at charge neutrality for M =9 ( N x = 3, N y = 3) moir´e unit cell systems at twist angle θ =1 .
1, interaction strength parameter eps − = 0 .
05, and inter-layer hopping parameters ω = 110 meV for inter sublatticetunneling and ω = 0 . ω for intrasublattice tunneling. Thevalence band is on the left and the conduction band on theright. The color maps are interpolations from discrete meshesindicated by the black squares; empty squares are used tomark momentum points that are equivalent to those includedin the discrete mesh. MANY-BODY EXACT DIAGONALIZATIONCALCULATIONSGround state quantum numbers estimation
The total Hilbert space of four possible flavors and two,valence and conduction bands, in a translationally invari-ant system is divided into smaller subspaces, with a givennumber of electrons in each valley N K and N K , and la-belled by quantum numbers: total spin S K and S K , az-imuthal spin S zK and S zK , and total momentum ( K x , K y )in two directions determined by modulo of a number ofunit cells N x and N y along two Moire lattice vectors. Thebasis is constructed in an occupation number configura-tions, distributing particles with a given spin z azimuthalquantum numbers on Hartree-Fock quasi-particle ( k x , k y )momentum states. The many-body Hamiltonian is diag-onalized in S zK and S zK subspaces, and we do not rotatea Hamiltonian matrix to a S K and S K basis as this isan additional computational cost, but rather determinethe ground state total spins based on S zK and S zK energyeigenvalues degeneracy. All of our exact diagonalization(ED) calculations are limited only by a maximal matrixsize of a given subspace. The largest subspace corre-sponds to the lowest possible pair of S zK , S zK , becausethey contain states with quantum numbers for all totalspins S K , S K . Using a notation ( N K ↓ , N K ↑ , N K ↓ , N K ↑ )we can distinguish possible particle distributions overfour flavors: the first two particle numbers correspondto two spins in valley K and the last two particle num-bers to two possible spins in valley K .We present a method of determining the ground statequantum numbers based on a system with N el = 9 parti-cles, ν = −
3. We have five possible N K and N K particledistributions with the lowest spin z quantum numbers: -159-158-157-156 (4,3,2,0) (4,2,3,0) (3,2,4,0) E [ m e V ] N conf FULL ED
FIG. S12: The lowest energy eigenvalues as a function of anumber of configurations N conf taken into truncated Hilbertspace calculations for distributions of N el = 9. Top: (4 , , , P v = 5 / , , ,
0) with P v = 3 / , , ,
0) with P v = 1 / (4 , , ,
0) with N K = 9 and N K = 0, (4 , , ,
0) with N K = 8 and N K = 1, (4 , , ,
1) with N K = 7 and N K = 2, (3 , , ,
1) with N K = 6 and N K = 3, (3 , , , N K = 5 and N K = 4, taking into account thefact that due to time-reversal and SU(2) × SU(2) symme-tries in the system, exchanging two valleys or two spinsgives the subspace with identical energy spectra, e.g.(4 , , , , , , , , , , , ,
4) have the samespectra, but energy spectra for (4 , , ,
0) and (4 , , , , , ,
0) and(4 , , ,
0) subspaces, but not for (4 , , , , , , , , , z subspaces full EDcalculations are accessible, e.g. for a particle distributionwith N K = 7 and N K = 2 we can do full ED calcula-tions for (7 , , , , , , , , , , , , S K and S K for the lowest eigenstate. This is estimated based on S zK and S zK degeneracy forming a given S K and S K multi-plet. In this case, one can not be completely sure whetherthe total spins of the ground state are not contained inthe subspace N K = 7 and N K = 2 and correspond to S K = 0 . S K = 0, due to lack of possibility of cal-culations for (4 , , , S K and S K gradually changes with their magnitude. Indeed, weobserve such a situation in most of the cases, the ener-gies of the largest possible S K and S K are usually thelowest, and when S K and S K decrease, their energy in-creases up to the highest energy with the lowest S K and S K , suggesting that full spin polarization is quite robustfor different fillings. Truncated Hilbert subspace calculations
The occupation number configurations can be sortedby their Hamiltonian expectation values and a truncatedHilbert space can be obtained by neglecting higher en-ergetic ones. One can determine approximately valleypolarization P v of the ground state by comparing thelowest energies of two or three different particle distribu-tions for calculations in a truncated Hilbert space, ana-lyzing results as a function of a truncated Hilbert spacesize. The eigenvalues within a given subspace convergemonotonically with a number of configurations taken,giving always an upper bound of the energy of a givenstate. Tempo of convergence usually depends on a ma-trix size, but considering different particle distributionsamong valleys but with the same or similar number ofparticles within a given flavor, gives matrices of the sameor similar size and a reliable comparison. While it maybe difficult to reach a satisfactory approximation to theenergy of the ground state by taking only a few per-centage from a total number of configurations, but basedon the energy difference between the lowest states as afunction of a truncated Hilbert space size, we can es-timate which of particle distributions is favorable, andhas a tendency to have the lower energy, as this differ-ence usually remains constant or even increases. We canperform similar comparisons for pairs or triads of parti-cle distributions getting approximate evaluation of val-ley polarization P v . In Fig. S12 we show energies as afunction of truncated Hilbert space size for particle dis-tributions (4 , , ,
0) ( P v = 5 / , , ,
0) ( P v = 3 / , , ,
0) ( P v = 1 / P v = 5 / ν = −
3, the lowestenergy in a given sector of valley polarization is alwayslower when compared with a smaller valley polarizationsector. After determining valley polarization P v , one cancompare the lowest eigenvalues from different large-spin z subspaces obtained from full ED and follow an estima-tion procedure described in a previous subsection in orderto approximately determine total spins S K and S K . Filling ν = − Using truncated Hilbert space calculations method de-scribed in the previous section, we analyze the nature of ν = −
2. Fig. S13 shows the low energy spectrum as afunction of Hilbert space size given by a number of config-urations taken N conf comparing three different particledistributions with the highest azimuthal quantum num-bers S zK and S zK . Two lowest states correspond to fullyvalley unpolarized state (9 , , ,
0) and valley unpolar-ized state (10 , , ,
0) with intervalley excitation from thestate (9 , , , , , , -278-276-274-272 (9,0,9,0) (10,0,8,0) (9,9,0,0) E [ m e V ] N conf FULL ED n=-2
Niu=-2 E g a p [ m e V ] N conf FULL ED
FIG. S13: Low energy spectrum as a function of a numberof configurations N conf taken into truncated Hilbert spacecalculations for three different particle distributions of N el =18, (9 , , ,
0) with P v = 0 (black squares), and (9 , , , P v = 1 (red circles), (10 , , ,
0) with P v = 1 //