An O(N) Ab~initio Calculation Scheme for Large-Scale Moiré Structures
Tan Zhang, Nicolas Regnault, B. Andrei Bernevig, Xi Dai, Hongming Weng
AAn O ( N ) Ab initio
Calculation Scheme for Large-Scale Moir´e Structures
Tan Zhang, Nicolas Regnault,
2, 3, ∗ B. Andrei Bernevig, Xi Dai, and Hongming Weng
1, 5, 6, † Beijing National Laboratory for Condensed Matter Physics,and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China Department of Physics, Princeton University, Princeton, NJ 08544, USA Laboratoire de Physique de l’Ecole normale sup´erieure,ENS, Universit´e PSL, CNRS, Sorbonne Universit´e,Universit´e Paris-Diderot, Sorbonne Paris Cit´e, Paris, France Department of Physics, Hong Kong University of Science and Technology, Kowloon, Hong Kong School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China
We present a two-step method specifically tailored for band structure calculation of the small-angle moir´e-pattern materials which contain tens of thousands of atoms in a unit cell. In the firststep, the self-consistent field calculation for ground state is performed with O ( N ) Krylov subspacemethod implemented in OpenMX. Secondly, the crystal momentum dependent Bloch Hamiltonianand overlap matrix are constructed from the results obtained in the first step and only a smallnumber of eigenvalues near the Fermi energy are solved with shift-invert and Lanczos techniques.By systematically tuning two key parameters, the cutoff radius for electron hopping interaction andthe dimension of Krylov subspace, we obtained the band structures for both rigid and corrugatedtwisted bilayer graphene structures at the first magic angle ( θ = 1 . ◦ ) and other three largerones with satisfied accuracy on affordable costs. The band structures are in good agreement withthose from tight binding models, continuum models, plane-wave pseudo-potential based ab initio calculations, and the experimental observations. This efficient two-step method is to play a crucialrole in other twisted two-dimensional materials, where the band structures are much more complexthan graphene and the effective model is hard to be constructed. Introduction. —Twisted bilayer graphene (TBG) hasrecently been intensely researched in both theory and ex-periment [1–7]. When the twisted angle ( θ ) is decreasedfrom large value down to the so-called first “magic angle”( θ = 1 . ◦ ), the system starts exhibiting correlated insu-lating phases [8–13] and unusual superconductivity [12–15] at special fillings in experiments. It is widely believedthat these novel transport properties originate from thenearly flat bands (FBs) that develop near the charge neu-trality point (Fermi energy of TBG) at the magic angle[2]. The measured bandwidth of the FBs is in the orderof 10-20 meV from tunneling spectroscopy experiments[11, 16–18]. Angle-resolved photoemission spectroscopymeasurements [19] find that the FBs are separated bya gap of 30-50 meV from both higher and lower energy“passive” bands. Although theoretical effective modelshave been shown to capture the essential properties ofthe FBs in TBG, the parameters of effective models areestimated according to experimental results which ofteninclude the correlation effects of Coulomb interaction. Itis not easy to get the accurate electronic structures purelyfrom ab initio calculation without involving any approx-imation for small angle cases. Untangling these correla-tions from the underlying single-particle band structureupon which they act becomes difficult and ab initio stud-ies of this problem are desirable.To solve this problem, researchers attempted to per-form ab initio calculations within density functional the- ∗ [email protected] † [email protected] ory (DFT) for TBG. Uchida et al . [20] carried out large-scale DFT calculations for structural relaxation and bandstructure within a real-space scheme. They found theFermi velocity of the Dirac electron dramatically reducedtoward zero when θ < ◦ , and it vanished at θ = 1 . ◦ and then slightly increased for the smaller θ . Song et al .[21] performed single-particle DFT calculations for rigidTBG using the Vienna Ab initio Simulation Package(VASP) based on plane-wave pseudo-potential method.They calculated the eigenenergies only at Γ and K mo-menta in the moir´e Brillouin zone (MBZ) when θ (cid:54) . ◦ ,while the whole band structure along high symmetricallines for θ = 1 . ◦ and larger ones. Recently, Cantele et al . [22, 23] obtained the relaxed TBG structures andthe bands at the magic angle and other larger ones byVASP. Their results show that the FBs at the first magicangle are connected with the atomic displacements origi-nating from the interlayer van der Waals interaction, andthe main contribution to the experimental gaps at theΓ point is the out-of-plane displacements. The compu-tational cost of all the aforementioned approaches is inorder of N , namely O ( N ), where N is the number ofatoms per unit cell, and N is more than tens of thousandwhen θ is around the first magic angle. The further studyof TBG with smaller angles within such DFT calculationis unaffordable.In this paper, we develop a two-step DFT calculationworkflow which can dramatically reduce the computa-tional cost while achieve satisfied accuracy. We performthis two-step calculation for both rigid and corrugatedTBG down to the first magic angle. The computationalcost is of O ( N ), which is much reduced comparing with a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b Check wherther the SCF calculations are converged
Yes NoObtain the real spaceHamiltonian and overlap matrixLocate the Fermi energy by countingthe occupation number at K pointCalculate 80 eigenenergies near theFermi energy by Lanczos techniquesConstructe the HamiltonianH( k ) and overlap matrix S( k ) by FFTObtain the band structurePerform O(N)
DFT SCF calculations by OpenMX
FIG. 1. (Color online) Calculation workflow contains the SCFcalculation by O ( N ) Krylov subspace method implemented inOpenMX and band structure calculation with shift-invert andLanczos techniques. that of VASP. The band structure, bandwidth of FBsand the band gap at K point have been shown anddiscussed. They are in good agreement with those ob-tained from tight-binding (TB) models, continuum mod-els, VASP calculations, as well as experimental observa-tions. This method can be applied on TBG with evensmaller angles and other twisted two-dimensional mate-rials, whose effective models are hard to be constructed,such as twisted trilayer graphene, twisted double bi-layer graphene, twisted bilayer transition-metal dichalco-genides. Workflow of the method. —The workflow is shown inFig. 1. The first step is to perform self-consistent field(SCF) calculation within DFT using the O ( N ) methodimplemented in OpenMX [24], which is an ab initio soft-ware package based on localized pseudo-atomic orbitals(PAOs) and pseudo-potential method. The details aboutthe crystal structure of rigid and corrugated TBG andthe calculation method are introduced in the Supplemen-tal Material [25]. The convergence criterion is set as thedifference in total energies of the last two SCF loops be-ing smaller than 10 − Hartree. If the SCF calculationis not converged within 200 loops, the parameters con-trolling the mixing of charge densities in previous loopswill be modified and restart the SCF calculation. If theSCF calculation is converged, we obtain the real spaceHamiltonian and the overlap matrix both in PAO ba-sis set. Then, the momentum k dependent Hamiltonian H ( k ) and overlap matrix S ( k ) can be constructed usingfast Fourier transformation (FFT), which is very similarto the interpolation method based on Wannier functionsconstructed from DFT band calculation. The dimensionof H ( k ) and S ( k ) is quite big in small angle case, and tosolve its eigenvalues is very challenge. Due to the local- ba site i site jr h O K (c)(b) B a nd g a p a t K po i n t ( m e V ) B a nd w i d t h ( m e V ) r h = 6 Å r h = 8 Å r h = 10 Å O K r h = 6 Å r h = 8 Å r h = 10 Å (a) FIG. 2. (Color online) The top view of crystal structure ofTBG with i = 16, θ = 2 . ◦ . It contains 3268 atoms andeach carbon atom site has thirteen localized PAOs as basisset. Two red circles represent truncated clusters centered onC atoms at sites i and j with r h = 8 ˚A. Bandwidth of FBs (b)and band gap at K point (c) obtained with different param-eters of hopping range r h and dimension of Krylov subspace O K for rigid TBG at θ = 1 . ◦ . ized PAO basis set and the cutoff of real space hoppinginteraction among electrons, H ( k ) and S ( k ) are sparsematrices. Lanczos technique combined with the shift-invert spectral transformation is suitable for solving onlya small number of eigenvalues (out of the large number ofeigenvalues at small angles) near a selected value of en-ergy, which is much faster than direct full diagonalizationfor all eigenvalues. The details about this techniques areintroduced in the Supplemental Material [25]. The Fermienergy is located by counting the occupation number ofthe bands at K by direct diagonalization. Finally, eightyeigenenergies near the Fermi level are calculated at a se-ries of momenta in the MBZ by shift-invert and Lanczostechniques and thus the band structures can been ob-tained. We have also implemented ScaLAPACK librarysupport for calculation, allowing us to obtain the fullspectrum down θ = 1 . ◦ providing the other consistencycheck. O ( N ) method. —The conventional DFT methods arenot efficient for large systems like TBG when θ is closeto magic angle. We employ the Krylov subspace O ( N )method implemented in OpenMX to perform the SCFcalculation. This method has been well explained inRef. 26, and we just briefly introduce the parts closelyrelated to our calculation. For a large system, the atomsinside a sphere centered on each site are physically trun-cated into small clusters, as shown in Fig. 2 (a). Theradius of the sphere is so-called hopping range r h . If thedistance between the neighboring sites and the centralatom is less than r h , the hopping of electrons betweenthem is considered. For each truncated cluster, the vec-tor space U C defined by basis functions in the truncatedcluster is mapped to a Krylov subspace U K . In general,the dimension of U K ( O K ) is smaller than that of U C ,and the dimension of U C is much smaller than that of to-tal vector space U F spanned by all of the basis functionsin the whole unit cell. Thus, the computational cost issignificantly reduced compared to the conventional DFTmethods.The key input parameters are r h and O K , which de-termine the balance between accuracy and computationalcost. To set their values properly, we test them for rigidTBG at θ = 1 . ◦ using this two-step method. The band-width of FBs and the band gap at K obtained from differ-ent r h and O K have been shown in Figs. 2 (b) and 2 (c).For r h = 6, 8 and 10 ˚A, the truncated cluster contains84 atoms, 168 atoms and 222 atoms, whose dimensionof Hamiltonian matrix is 1092, 2184, 2886, respectively,and the results are converged when O K = 1000, 2000 and2600, respectively, which are about 8.9% smaller than thedimension of the Hamiltonian matrix ( U C ) on avarage.However, it is noted that the cutoff radius of PAOs forbasis function of C atom is 3.18 ˚A, and the maximum dis-tance between two C sites with their PAOs overlapped is6.36 ˚A. The r h = 6 ˚A is too small for the truncated clus-ters. On the other hand, r h = 10 ˚A is too large, whichincludes too many useless non-orthogonal basis functionswithout any hopping interaction with the central C atomin the truncated clusters. This leads to the over com-pleteness of the basis set and gives out quite strange re-sults. Therefore, r h = 8 ˚A is quite reasonable choice. Infact, the band structures with these three r h are somehowdifferent from each other. We check these bands againstthose from TB models, continuum models and VASP cal-culations. The most reliable one is from calculations with r h = 8 ˚A. Because the cutoff radius of PAOs and the to-tal number of atoms in a truncated cluster with the same r h hardly change at different θ , the appropriate r h and O K should be independent of θ . This is verified by ourfurther tests for both rigid and corrugated structures atdifferent θ . Hence, we use r h = 8 ˚A and O K = 2000 incalculations for all TBG structures. Calculation efficiency. —We show the computationalcost of rigid TBGs at four different θ values in Fig. 3.These calculations are performed with Intel Xeon CPUE3-1240 v6 family at base frequency of 3.7 GHz. Byusing O ( N ) method, the computational cost, in unit ofcore · hour, of one SCF loop calculation almost linearly in-creases with the number of atoms per unit cell. BecauseTBG contains tens of thousands of atoms in a unit cellat small angles, the O ( N ) method is considerably fasterthan the conventional DFT methods of O ( N ). The costfor obtaining 80 eigenenergies near Fermi energy at onemomentum is also nearly O ( N ), because the dimensionsof Hamiltonian matrix and overlap matrix are of O ( N ).In Fig. 3, one SCF loop takes about 4.3 core · h, 6.6 core · h,8.7 core · h, 14.6 core · h, and one band calculation at one C o s t o f on e S C F l oop ( c o r e ∙ h ) Atoms number per unit cell
One SCF loopEigenenergies C o s t o f e i g e n e n e r g i e s ( c o r e ∙ h ) FIG. 3. (Color online) The computational cost of one SCFloop (red dots) and selected eighty eigenenergies solved atone momentum (blue squares) for rigid TBG of four twistedangles. The unit cell contains 3268, 5044, 6628 and 11164 car-bon atoms at θ = 2 . ◦ , 1 . ◦ , 1 . ◦ and 1 . ◦ , respectively. momentum point needs about 21 core · h, 49 core · h, 67core · h, 133 core · h at θ = 2 . ◦ , θ = 1 . ◦ , θ = 1 . ◦ , θ = 1 . ◦ , respectively. We also get the similar resultsfor corrugated structures. As a reference, Ref. 21 statedthat the cost of about 576 core · h, 1728 core · h and 17280core · h is needed to finish one band calculation at onemomentum point by VASP at θ = 2 . ◦ , θ = 1 . ◦ and θ = 1 . ◦ , respectively. Our band calculation is about130 times faster than that of VASP at θ = 1 . ◦ . It isalso noted that the total computational cost by VASPrequires 2 . × core · h to 4 . × core · h including onerelaxation calculation, one SCF calculation and band cal-culations at seven momentum points of TBG at θ = 1 . ◦ [22, 23]. Results of the method. —The band structures for rigidand corrugated TBGs at four twist angles are shown inFig. 4. The red points represent the bands calculatedby the two-step method, and the red lines are plotted asguide to the eye. The blue lines are bands obtained byTB models [27]. They are in quite good agreement, aswell as with those from continuum models [2, 28]. Wefirstly focus on the bands of rigid structures. The band-width of the four bands (two per valley in the continuummodels) near the Fermi energy monotonically decreaseswith θ as shown in Fig. 5 (a). When θ = 1 . ◦ , thebandwidth decreases to 21.0 meV and the FBs appearat the Fermi energy (see Fig. 4(d)). They are separatedfrom the lower energy bands by a gap of 8.4 meV. Thisis a little different from those calculated from both TBmodel and continuum model where the FBs entangledwith the lower energy bands at Γ. The other differenceis that there exists a gap among these four bands at K point from this method and VASP [21], because the inter-valley coupling in TBG is not ignored in DFT calculation,while it is four-fold degenerate at K point in model cal-culations where the inter-valley coupling is ignored. Thisgap gradually decreases with θ and becomes very smallat the first magic angle (0.12 meV) of θ = 1 . ◦ as shownin Fig. 5 (b). (a) M K Г E n e r gy ( e V ) M i = 16, θ = 2.00° M K Г -0.2-0.1 0 0.1 0.2 M G K M E n e r g y ( e V ) MRigid structure Corrugated structure -0.2-0.1 0 0.1 0.2 M G K M E n e r g y ( e V ) -0.1-0.05 0 0.05 0.1 M G K M E n e r g y ( e V ) (b) M K Г E n e r gy ( e V ) M i = 20, θ = 1.61° -0.1-0.05 0 0.05 0.1 M G K M E n e r g y ( e V ) M K Г MRigid structure Corrugated structure -0.1-0.05 0 0.05 0.1 M G K M E n e r g y ( e V ) (c) M K Г E n e r gy ( e V ) M i = 23, θ = 1.41° -0.1-0.05 0 0.05 0.1 M G K M E n e r g y ( e V ) M K Г MRigid structure Corrugated structure -0.1-0.05 0 0.05 0.1 M G K M E n e r g y ( e V ) (d) M K Г E n e r gy ( e V ) M i = 30, θ = 1.08° -0.1-0.05 0 0.05 0.1 M G K M E n e r g y ( e V ) M K Г MRigid structure Corrugated structure
FIG. 4. The band structures obtained from our two-step method (red dots) and TB model (blue lines) for rigid (left) andcorrugated (right) TBGs at twisted angles: (a) θ = 2 . ◦ , (b) θ = 1 . ◦ , (c) θ = 1 . ◦ and (d) θ = 1 . ◦ . B a nd g a p a t K po i n t ( m e V ) (b) (a) B a nd w i d t h ( m e V ) RigidCorrugatedRigidCorrugated
FIG. 5. (Color online) Bandwidth of the four bands nearFermi level (a) and band gap at K point (b) for rigid (reddots) and corrugated (blue squares) TBGs at four twistedangles: θ = 2 . ◦ , 1 . ◦ , 1 . ◦ , and 1 . ◦ . Next, we discuss the results for corrugated structures.The corrugation enlarges both the bandwidth and thegap at K point when θ > . ◦ , but it reduces them atthe magic angle θ = 1 . ◦ (see Fig. 5). We focus on theband structure at θ = 1 . ◦ . The four bands becomeFBs with bandwidth of about 13.2 meV, and the gap at K point decreases to 0.08 meV. The two gaps separat-ing FBs from the higher and lower energy bands increasedue to the corrugations [27, 29], being 15.3 meV and11.3 meV, respectively. These results are consistent withthose from TB models [27, 29–34], continuum models[27, 35], VASP calculations with fully relaxed structures[22, 23] and experimental observations [8, 11, 14, 16–19]. This implies that atomic corrugations in TBG play animportant role in the electronic structures, and that themain effects of corrugations in real materials can be sim-ulated by modeled corrugations. Conclusions. —We design a two-step workflow to carryout DFT calculation of band structures for large sys-tems like TBGs with tens of thousands of atoms. In thefirst step, a SCF calculation using O ( N ) Krylov subspacemethod implemented in OpenMX is performed to obtainthe Hamiltonian and overlap matrix in real space. Inthe second step, the momentum dependent Bloch Hamil-tonian and overlap matrix are constructed in the senseof interpolation, and only a small number of eigenener-gies near the Fermi energy are solved by shift-invert andLanczos techniques. The key input parameters in O ( N )SCF calculation are r h and O K , which can be systemat-ically tuned to balance the accuracy and computationalcost. We find that r h = 8 ˚A and O K = 2000 are suitablefor TBG and that these values do not depend much on θ and corrugation. The computational cost linearly in-creases with the number of atoms per unit cell for bothrigid TBGs and those with corrugation. For rigid TBG at θ = 1 . ◦ , one SCF loop requires about 14.6 core · h andone band calculation of 80 eigenvalues at one momen-tum point requires about 133 core · h on our computers.This is affordable for even larger system with smallerangles. The band structures obtained at four twistedangles by this method are consistent with those of TBmodels, continuum models, VASP calculations and theexperimental observations. This method can be appliedto other twisted materials and could become an efficienttool when no effective model is available.We acknowledge the supports from the National Nat-ural Science Foundation (Grant No. 11925408), theMinistry of Science and Technology of China (GrantsNo. 2016YFA0300600 and 2018YFA0305700), the Chi-nese Academy of Sciences (Grant No. XDB33000000),the K. C. Wong Education Foundation (GJTD-2018-01),the Beijing Natural Science Foundation (Z180008), andthe Beijing Municipal Science and Technology Commis-sion (Z191100007219013). B. A. B and N. R. were alsosupported by the DOE Grant No. DE-SC0016239, the Schmidt Fund for Innovative Research, Simons Investi-gator Grant No. 404513, the Packard Foundation, theGordon and Betty Moore Foundation through Grant No.GBMF8685 towards the Princeton theory program, anda Guggenheim Fellowship from the John Simon Guggen-heim Memorial Foundation. Further support was pro-vided by the NSF-EAGER No. DMR 1643312, NSF-MRSEC No. DMR-1420541 and DMR-2011750, ONRNo. N00014-20-1-2303, Gordon and Betty Moore Foun-dation through Grant GBMF8685 towards the Princetontheory program, BSF Israel US foundation No. 2018226,and the Princeton Global Network Funds. [1] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H.Castro Neto, Phys. Rev. Lett. , 256802 (2007).[2] R. Bistritzer and A. H. MacDonald, Proceedings of theNational Academy of Sciences , 12233 (2011).[3] G. Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Cas-tro Neto, A. Reina, J. Kong, and E. Y. Andrei, NaturePhysics , 109 (2010).[4] D. S. Lee, C. Riedl, T. Beringer, A. H. Castro Neto,K. von Klitzing, U. Starke, and J. H. Smet, Phys. Rev.Lett. , 216602 (2011).[5] P. San-Jose, J. Gonz´alez, and F. Guinea, Phys. Rev.Lett. , 216802 (2012).[6] P. Maher, C. R. Dean, A. F. Young, T. Taniguchi,K. Watanabe, K. L. Shepard, J. Hone, and P. Kim,Nature Physics , 154 (2013).[7] B. M. Hunt, J. I. A. Li, A. A. Zibrov, L. Wang,T. Taniguchi, K. Watanabe, J. Hone, C. R. Dean, M. Za-letel, R. C. Ashoori, and A. F. Young, Nature Commu-nications , 948 (2017).[8] 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).[9] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney,K. Watanabe, T. Taniguchi, M. A. Kastner, andD. Goldhaber-Gordon, Science , 605 (2019).[10] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora,R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. vonOppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge,Nature Physics , 1174 (2019).[11] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian,M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi,J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Na-ture , 95 (2019).[12] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir,I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang,A. Bachtold, A. H. MacDonald, and D. K. Efetov, Na-ture , 653 (2019).[13] E. Codecido, Q. Wang, R. Koester, S. Che, H. Tian,R. Lv, S. Tran, K. Watanabe, T. Taniguchi, F. Zhang,M. Bockrath, and C. N. Lau, Science Advances (2019),10.1126/sciadv.aaw9770.[14] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,E. Kaxiras, and P. Jarillo-Herrero, Nature , 43(2018).[15] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watan- abe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean,Science , 1059 (2019).[16] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora,R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. vonOppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge,Nature Physics , 1174 (2019).[17] Y. Xie, B. Lian, B. J¨ack, X. Liu, C.-L. Chiu, K. Watan-abe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Na-ture , 101 (2019).[18] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule,J. Mao, and E. Y. Andrei, Nature , 91 (2019).[19] S. Lisi, X. Lu, T. Benschop, T. A. de Jong, P. Stepanov,J. R. Duran, F. Margot, I. Cucchi, E. Cappelli,A. Hunter, A. Tamai, V. Kandyba, A. Giampietri,A. Barinov, J. Jobst, V. Stalman, M. Leeuwenhoek,K. Watanabe, T. Taniguchi, L. Rademaker, S. J. van derMolen, M. P. Allan, D. K. Efetov, and F. Baumberger,Nature Physics (2020), 10.1038/s41567-020-01041-x.[20] K. Uchida, S. Furuya, J.-I. Iwata, and A. Oshiyama,Phys. Rev. B , 155451 (2014).[21] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A.Bernevig, Phys. Rev. Lett. , 036401 (2019).[22] P. Lucignano, D. Alf`e, V. Cataudella, D. Ninno, andG. Cantele, Phys. Rev. B , 195419 (2019).[23] G. Cantele, D. Alf`e, F. Conte, V. Cataudella, D. Ninno,and P. Lucignano, Phys. Rev. Research , 043127 (2020).[24] T. Ozaki, Phys. Rev. B , 155108 (2003).[25] See Supplemental Material at for detailed informationabout the crystal structure of rigid and corrugated TBG,the calculation method, and shift-invert and Lanczostechniques..[26] T. Ozaki, Phys. Rev. B , 245101 (2006).[27] J. Liu, J. Liu, and X. Dai, Phys. Rev. B , 155415(2019).[28] P. Moon and M. Koshino, Phys. Rev. B , 205404(2013).[29] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi,K. Kuroki, and L. Fu, Phys. Rev. X , 031087 (2018).[30] N. N. T. Nam and M. Koshino, Phys. Rev. B , 075311(2017).[31] S. Carr, S. Fang, P. Jarillo-Herrero, and E. Kaxiras,Phys. Rev. B , 085144 (2018).[32] N. Leconte, S. Javvaji, J. An, and J. Jung, arXive-prints , arXiv:1910.12805 (2019), arXiv:1910.12805[cond-mat.mes-hall].[33] M. Angeli, D. Mandelli, A. Valli, A. Amaricci, M. Capone, E. Tosatti, and M. Fabrizio, Phys. Rev.B , 235137 (2018).[34] A. O. Sboychakov, A. L. Rakhmanov, A. V. Rozhkov,and F. Nori, Phys. Rev. B , 075402 (2015). [35] S. Carr, S. Fang, Z. Zhu, and E. Kaxiras, Phys. Rev.Research1