# Fileset

[著者最終稿.pdf](https://mdr.nims.go.jp/filesets/da7c33ad-5945-4992-9931-f8426743874f/download)

## Creator

[Toshikaze Kariyado](https://orcid.org/0000-0002-3746-6803), Ashvin Vishwanath, Zhu-Xi Luo

## Rights

©2025 American Physical Society[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[Single-band square-lattice Hubbard model from twisted bilayer <math>  <msub>    <mi>C</mi>    <mn>568</mn>  </msub></math>](https://mdr.nims.go.jp/datasets/9a52295b-3c50-408f-a4c2-1b605fcd1904)

## Fulltext

Single-band square lattice Hubbard model from twisted bilayer C568Toshikaze Kariyado,1 Ashvin Vishwanath,2 and Zhu-Xi Luo2, 31International Center for Materials Nanoarchitectonics (WPI-MANA),National Institute for Materials Science, Tsukuba 305-0044, Japan2Department of Physics, Harvard University, Cambridge, MA 02138, USA3School of Physics, Georgia Institute of Technology, Atlanta, Georgia 30332, USA(Dated: December 24, 2025)We propose twisted homobilayer of a carbon allotrope, C568, to be a promising platform to realizecontrollable square lattice single-band extended Hubbard model. This setup has the advantageof a widely tunable t′/t ratio without adding external fields, and the intermediate temperaturet ≪ T ≪ U regime can be easily achieved. We first analyze the continuum model obtained fromsymmetry analysis and first-principle calculations, and calculate the band structures. Subsequently,we derive the corresponding tight-binding models and fit the hopping parameters as well as theCoulomb interactions. When displacement field is applied, anisotropic nearest neighbor hoppingscan further be achieved. If successfully fabricated, the device could be an important stepping stonetowards understanding high-temperature superconductivity.CONTENTSI. Introduction 1II. Modeling from first-principles 2III. Twisted bilayer model 4A. Displacement field 4IV. Tight binding model 5A. The hopping parameters 5B. Hubbard U parameters 6V. Discussion 6Acknowledgements 7A. Density functional theory 8B. Lower harmonics expansion and the least square fit8C. Interlayer distance 8D. Emergent symmetry L 9References 9I. INTRODUCTIONOne of the most highly prized goals in condensed mat-ter physics is to realize controllable high-temperature su-perconductivity. The cuprate materials comprise a fa-mous family of high-temperature superconductors [1],which are fascinating but also challenging to understand.Recently, the rise of van der Waals heterostructures ashighly controllable quantum simulators of various quan-tum systems [2], especially Hubbard models [3, 4] natu-rally lead us to wonder: can van der Waals materials sim-ulate cuprate physics and particularly high-temperaturesuperconductivity?It is widely believed that the simplest effective modelthat captures essential cuprate physics is that of thesingle-band square-lattice Hubbard model , see for ex-ample [5]. To simulate the latter, a Moiré superlatticewith square unit cells is necessary. This can be achievedthrough various setups, including twisted homobilayersof square lattice materials, heterobilayers of square lat-tice materials with mismatched lattice constants, two-dimensional materials from other crystal groups twistedat special angles, or multilayer heterostructures. For ex-ample, recently the possibility of using 90◦-twisted ho-mobilayer of the GeSe family to simulate cuprate physicshave been proposed [6]. While the monolayer only hasC2v symmetry, the 90◦-twisted stacking generates an ef-fective moiré square lattice model. The tunability of thispromising setup, however, is limited as the angle is fixedand the moiré scale is not significantly larger than thatof the parent lattice scale. Across various possible setupsmentioned above, the twisted homobilayer stands out asthe most minimalist and versatile choice.To the end of simulating a single-orbital square latticeHubbard model, the monolayer material should have itsband extremum sitting at the Brillouin zone center (Γpoint) or corner (M point). We will focus on the lattercase as it is believed to have stronger interference effects[7], leading to deeper moiré potentials. Ref. [8] used aminimal theoretical model to examine this possibility andfound a widely tunable t − t′ − U Hubbard model whenthe displacement field is present. While some monolayersquare lattice materials have been fabricated in the lab(see for example [9–12]), simple non-degenerate M -pointmaterials are sparse [13] and have not been experimen-tally realized to our best knowledge. One promising can-didate is the monolayer carbon allotrope C568, shown tobe stable via first-principle calculations by independentgroups [14–17]. Its geometry contains pentagon, hexagonand octagon plaquettes, therefore the name. The con-duction band minimum is predicted to be at the Γ-point,while the valence band maximum is at the M -point. Al-though hosting a square lattice structure, the material is2non-planar and buckled such that the C4z symmetry isabsent. Instead, there is a S4 symmetry which is the com-position of C4z and reflection with respect to the mono-layer plane. This distinguishes C568 from the previoustheoretical analysis in ref. [8].In this work, we examine twisted homobilayers of C568.Through symmetry and first-principle analyses, we writedown the bilayer continuum Hamiltonian with realisticparameters. Upon studying the band structures, bothwith and without external displacement fields, we derivethe corresponding effective moiré Hubbard models. Wefind the effective model without field to be that of thesingle-band Hubbard model on the square lattice, unlikethe case considered in Ref. [8]. Upon tuning the twistingangles, the nearest and next-nearest neighbor hoppingstrengths cross at around 8.5◦. A widely tunable t′/t ra-tios plays an important role in realizing superconductiv-ity in square lattice t−t′−U Hubbard model [18–21] andt−t′−J models [22–24], and also leads to the potential offrustrated magnetism [25–28]. When displacement fieldis further added, an anisotropic nearest neighbor hoppingtx ̸= ty can be reached, leading to imbalanced band flat-tening in the two directions, also found in moiré systemswith generic symmetry, see for example [7, 29, 30]. Weestimate the magnitude of both the onsite and the longerranged the Coulomb interactions. At smaller twisting an-gles and fields, the system is in the strongly interactingregime. For example, when θ = 3◦ and D = 0, we havethe ratio U/tx ∼ 105, and tx+y/tx = −0.7639.We conclude the introduction by pointing out some ofthe other works on moiré square lattice systems whichserve different purposes, including: twisted bilayers ofuniform [31] or staggered flux states [32], states withquadratic band touching [33], twisted bilayers of cuprates[34–38] and FeSe [39].II. MODELING FROM FIRST-PRINCIPLESThe lattice structure of monolayer C568 is shown inFigs. 1(a) and 1(b). Interestingly, the network consistsof mixture of sp3 and sp2 bondings, where the four sp3bonds locate at the center of unit cell in Fig. 1(b). Due tothe sp3 nature, the monolayer is buckled as clearly seenin Fig. 1(a). A key monolayer symmetry is roto-inversionS4 associated with the sp3 bonds, which can equivalentlybe expressed as C4zMxy, where Mxy is a reflection withrespect to the plane of monolayer. Note that S4 automat-ically gives C2z with respect to the S4 center. (Due to thelattice translation symmetry, there appears the other S4center at the corner of a unit cell.) The monolayer breaksinversion symmetry. The other important symmetry isreflection with respect to the plane perpendicular to thelayer and running along horizontal or vertical directionsin Fig. 1(b). The plane passes through either the center(M1) or the corner (M2) of a unit cell. There is also aC2 axis along the diagonal direction in Fig. 1(b). Thisin-plane C2 operation can be decomposed into Mxy anda reflection with respect to the diagonal lines in Fig. 1(b),which we name Md.Monolayer C568 has been predicted [14] to be a semi-conductor with its valence top at theM -point in the Bril-louin zone. These crystalline and electronic structuresare reproduced in our density functional theory (DFT)calculations, see Appendix A for computational details.In this paper, we focus on electronic structures near thevalence top.The moiré pattern in twisted bilayer stems from theposition dependence of relative shift τ between two lay-ers. For small twist angles, a moiré system is locally wellapproximated by untwisted bilayer with relative lateralshift τ , and investigating τ -dependence of band energynear the valence top in the untwisted bilayer gives infor-mation for constructing a moiré twisted bilayer model.Since monolayer C568 breaks spatial inversion symmetry,there are two types of stacking patterns: one with twoidentical copies of monolayers and one between two in-version images of each other. The former still breaks in-version symmetry after stacking while the latter restoresinversion symmetry after stacking. In this paper, we fo-cus on the former with the aim of realizing the squarelattice Hubbard model in twisted bilayers.For different relative shifts τ , the interlayer distancesdz(τ ) can be different. Therefore, such corrugation effectin the moié bilayers must be taken into account whencomparing the bilayer band structures for different τ ’s.For a fixed stacking pattern τ , we calculate the corre-sponding total energy using DFT for various interlayerdistances, and fix the optimal dz(τ ) to be the one thatminimizes the energy. (See Appendix C for details.) Be-cause of layer doubling, there are two bands near the va-lence top, and the energies of these two bands are used inextracting parameters for an effective model. Figure 1(c)shows bilayer band structures for several selected τ s. Fo-cusing on the valence top at M-point, sizable τ depen-dence is observed, indicative of strong band reconstruc-tion in twisted bilayer cases.The effective continuum model for the untwisted bi-layer expanded near the Brillouin zone corner kM can beexpanded in the layer basis asH = H0(kM )σ0 + Vpot(τ ), Vpot(τ ) =3∑a=0Va,kM(τ )σa,(1)where σa labels the Pauli matrices acting in the layerspace and σ0 is identity. The functions H0 and Vpot bothdepend on the base momentum kM and are subject tosymmetry constraints of the bilayer system. Vpot furtherdepends on the relative displacement between the twolayers τ .The associated energies areE±(kM , τ ) = H0(kM )+V0,kM(τ )±(3∑a=1Va,kM(τ )2)1/2.(2)3−5.4−5.2−5.0−4.8−4.6M← X Γ→Energy[eV]−7−6−5−4−3Y Γ X M Γτ = (0, 0)(0, .5)(.5, 0)(.5, .5)Energy[eV](a)(b)S4S4C2M2 M1(c)(d)FIG. 1. Side view (a) and the top view (b) of the lattice structure of monolayer C568. In the top view, a unit cell is alsoindicated. The yellow and cyan dotted lines are for the reflections M1 and M2, respectively. The dash-dotted lines are forC2 = MdMxy. The red cross marks are for the S4 centers. (c) Band structures of untwisted bilayers with several selectedrelative shifts τ . In the picture, τ is indicated in the unit of a0. (d) The close up view of the boxed area in (c).We choose the zero-energy level H0(kM ) = 0 for conve-nience and will omit the subscript kM in Va from now on.The V -functions can be determined by comparing theseenergies with the DFT results.The symmetry of Vpot is essential for constructing aneffective model. In the monolayer case, roto-inversionsymmetry S4 = C4zMxy. For bilayers, Mxy inter-changes two layers, therefore the S4 constraint readsVpot(C4zτ ) = σ1Vpot(τ )σ1. (3)More specifically V0,1(C4zτ ) = V0,1(τ ) and V2,3(C4zτ ) =−V2,3(τ ). The monolayer in-plane C2 = MdMxy, andleads toVpot(Mdτ ) = σ1Vpot(τ )σ1. (4)Specifically this gives constraints V0,1(Mdτ ) = V0,1(τ )and V2,3(Mdτ ) = −V2,3(τ ). Lastly, since the M-point isa high symmetry point in the Brillouin zone, quantum in-terference can eliminate interlayer tunneling for specific τpoints [7, 40]. In this case, the reflection symmetries M1and M2 in fig. 1 in monolayer lead to V1,2(τ ) = 0 whenτ is at the boundary of a unit cell. In addition, on thediagonal lines in the unit cell which satisfy Mdτ = τ , V3vanishes. Also, since M-point is time reversal invariantmomentum, there exists a basis such that the Hamilto-nian becomes real. In that basis, V2 is zero and there re-main three functions to be determined. We remark thatif the inversion symmetry was present, then the V3(τ )term would not have been allowed, which was the case in[8].However, easily obtained inputs from DFT are onlytwo functions, E±(kM , τ ). To overcome this issue, weadopt lower harmonics expansion of Vi respecting crys-talline symmetries, and fix required coefficients by theleast square fitting to DFT data. (See Appendix B fordetails.) The prepared lower harmonics expansions areV0(τ ) = γ0 + γ1(cos τ̃x + cos τ̃y) + 2γ2 cos(τ̃x) cos(τ̃y)+ γ3(cos 2τ̃x + cos 2τ̃y)+ 2γ4(cos(2τ̃x) cos(τ̃y) + cos(τ̃x) cos(2τ̃y)),V1(τ ) =α1 cosτ̃x2cosτ̃y2+ α2(cos3τ̃x2cosτ̃y2+ cosτ̃x2cos3τ̃y2),V3(τ ) =β1(cos τ̃x − cos τ̃y) + β2(cos 2τ̃x − cos 2τ̃y)+ β3(cos 3τ̃x − cos 3τ̃y),(5)where τ̃µ = 2πτµ/a0. These expansions satisfy all thesymmetry constraints described above.When fitting with the DFT data, coefficients in V0 canbe determined using a data set (E++E−)/2. The remain-ing data set (E+ − E−)/2 can be used for determiningtwo functions V1 and V3: at the boundary of a unit cell,V1 = 0 and (E+ − E−)/2 is solely determined from V3;on the other hand, on the diagonal lines in a unit cellthat satisfies Mdτ = τ , the symmetry constraints leadsto V3 = 0, and (E+ − E−)/2 is solely determined fromV1. Using these facts, we obtainα1 = 0.194, α2 = −0.052, γ0 = −2.292,γ1 = −0.011, γ2 = 0.021, γ3 = −0.007, γ4 = −0.001,β1 = −0.024, β2 = −0.008, β3 = −0.011,(6)all in the unit of eV. The lower harmonics expansionand the DFT results are compared in Fig. 2. Consid-ering the simplicity of the approximation and limitationin converting E± data set to the coefficients, the matchbetween the data and the lower harmonics expansion isremarkable.40.00.51.0 0.00.51.0−2.5−2.4−2.3−2.2−2.1E−E+τx/a0τy/a0Energy[eV]FIG. 2. E+(kM , τ ) and E−(kM , τ ) obtained with the lowerharmonics expansion (lines) and DFT (points).III. TWISTED BILAYER MODELIn this section, we specify to the rigid twist case withdispacement τ (r) = 2 sin θ2 ẑ × r. The moiré lattice con-stant is thus aM ≈ a0/θ at small twisting angles. Wewill denote the moiré lattice vectors by R1,R2 and themoiré reciprocal lattice vectors by q1, q2. The continuumHamiltonian readsH =(H−(−i∇) + V0 + V3 eiδk·rV1e−iδk·rV1 H+(−i∇) + V0 − V3),(7)where H±(−i∇) = ℏ22m∗(−i∇ − kM ± δk/2)2 and δk =2 sin θ2 (ẑ × kM ) = (q1 − q2)/2, which is depicted in Fig.3.When V0 = V3 = 0 and in (7) and keeping only theleading term in V1, i.e., V1 = α1 cos(πa0τ1) cos(πa0τ2) =α2 cos(πaMy) cos( πaMx), the model reduces to that studiedin [8]: there is an emergent symmetry at small twistingangles L0 = σxeiσzδk·r. This symmetry forces the bandmaxima to be degenerate EΓ = EM leads to the van-ishing of nearest neighbor hopping in the effective tightbinding model defined on the moiré lattice. Adding V0as well as the sub-leading terms in V1 don’t change thephysics.In C568, inversion symmetry is broken, such that a fi-nite V3 term is allowed which breaks L0. However, adifferent emergent symmetry, which involves the four-fold rotation, becomes present at small twisting an-gles. To understand the effect of the emergent symme-try on the moiré momentum space, we apply the fol-lowing convenient unitary transformation U = eik0·rσ0where k0 = kM + (q1 + q2)/4 , such that H ′ = U†HUhas similar form as in (7) but with H± replaced byH ′± = ℏ22m∗(−i∇− kM + k0 ± δk/2)2. We can now labelany momentum using its shift from k0 : q = k−k0, wherek0 is the Γm point in the shifted Brillouin zone plottedin the right panel of figure 3. Then we conveniently haveH ′+ =ℏ22m∗(−i∇+ q + q1/2)2,H ′− =ℏ22m∗(−i∇+ q + q2/2)2.(8)It is easy to check (details are shown in appendix D) thatthe following operation leaves the Hamiltonian invariant:L = e−i2q2·r(σ3+σ0)/2σ1C4z, H ′(q) = LH ′(Rπ/2q)L−1.(9)Mbot.MtopΓδkΓmMmXmδkYmq1q2FIG. 3. Left: Brillouin zones of the top and bottom layers.Dashed blue square is the moiré Brillouin zone, and dashedred square is the shifted moiré Brillouin zone under the trans-formation U . Right: High symmetry points in the shiftedmoiré Brillouin zone. Γm point has momentum k0.Consequently, as long as the band is non-degenerate,the energy of the Hamiltonian H ′ must satisfyE(q) = E(Rπ/2q), (10)which leads to E(Xm) = E(Ym) in the shifted moiréBrillouin zone shown in fig. 3.Using the fitted parameters found in (6) and keepingonly those whose magnitudes are larger than 0.01 eV ,we show the band structure E(q) along high symmetrylines in fig. 4: the band top lives at the Mm point, andthe band bottom is at the Γm point in the shifted Bril-louin zone. The shoulders at Xm and Ym points have thesame energies consistent with (10). The band structure4 clearly resembles that of the tight-binding Hamiltonianon square lattice with nearest neighbor hopping, whichwe will explain in more details in section IVA. This is incontrast to the case in ref. [8] where an external field isnecessary for reaching the same dispersion. At θ = 3.6◦,the flat moiré band is separated from the neighboringband by a large gap of ∼ 19.5 meV, which is about 28times the bandwidth of the top flat band.A. Displacement fieldUpon adding a displacement field Dσz, the emergentL symmetry which involves layer exchange is absent. In5FIG. 4. Left: The top moiré band of (7) at θ = 3.6◦ (blue).This clearly resembles the dispersion of tight binding modelon square lattice (red dashed). Right: the next band is well-separated.FIG. 5. Band structures at different values of displacementfields. Left: the top two moiré bands are well-separated atD = −120 meV. Right: the top two moiré bands cross at theΓ−X high symmetry line when D = −170 meV, but are stillwell separated from other bands.the shifted moiré Brillouin zone plotted in the right panelof fig. 3, this translates into the asymmetry EX ̸= EY .Taking into account all the terms in the Hamiltonian, asmall and negative D already shifts the band extrema:EY > EM > EΓ > EX . Further increasing the mag-nitude, when D = −0.17 eV, band crossings betweenthe top two moiré bands are induced, and beyond whichEM < EΓ. The band crossings always happen alongthe Brillouin zone boundary which respects the reflec-tion symmetry. See Fig. 5 for typical band structures.For the a positive D, the evolution of band structures issimilar but with X ↔ Y exchanged.IV. TIGHT BINDING MODELIn this section we analyze the effective tight-bindingmodel for the twisted bilayer. At small twisting angles,it is sufficient to include only the nearest neighbor tx, tyand the next nearest neighbor hoppings tx+y, tx−y. Timereversal symmetry requires them to be real. Reflectionsymmetries with respect to the x- and y- axes requirethe next nearest neighbor hoppings, tx+y = tx−y. Thecorresponding dispersion is:E(k) = E0 − tx cos kx − ty cos ky − 2tx+y cos kx cos ky.(11)0.01.02.03.04.00 3 6 9 12 15tx−tx+y−t2xt2x+yHopping[meV]Twist angle [degree]−0.20.00.20.40 20 40 60 80 100θ = 3◦txty−tx+yD [meV](a) (b)FIG. 6. Hopping parameters derived using full parametersand Wannier functions. (a) Twist angle dependence. tx, tx+y,t2x, and t2x+y correspond to the first to the fourth nearestneighbor hoppings in the square lattice. (b) Displacementfield dependence for θ = 3◦. Insets: total charge (sum ofthe charges on the upper and the lower layer) for a Wannierfunction at ϕ = 3◦ with (a) D = 0 and (b) D = 50 meV.The emergent L symmetry, when present, imposes thatthe nearest neighbor hopping tx = ty.A. The hopping parametersThe four parameters {E0, tx, ty, tx+y} can be eas-ily solved by requiring that at high symmetry pointsΓm, Xm, Ym,Mm defined in fig. 3, the energies computedfrom the twisted bilayer model and that from the tightbinding model (11) should match with each other. Thismethod works well at small twisting angles where thelonger ranged hoppings can be ignored. For example atD = 0 and θ = 3.6◦, we have tx = ty = 0.171 meV,t1 = 0.014 meV. We plot the dispersions calculated bothfrom the tight-binding model and the lowest-harmonicstwisted bilayer model in fig. 4(a) and they match nicely.While one can also plot t as a function of twisting angleand displacement field using this simple method, belowwe will perform a more careful numerical examination.The elaboration is done in the following aspects: (i)taking all γi, αi, and βi into account, and (ii) deriving theWannier orbitals for the highest energy band. In the pa-rameter range showing no crossing between the top andsecond bands, the highest energy band is well isolatedfrom the other bands, and it is justified to have reason-ably localized Wannier orbitals. In practice, we prepareBloch wave functions for the highest energy band in thecontinuum Hamiltonian Eq. (7) (with full account of γi,αi, and βi) on 13× 13 regular grid in the moire Brillouinzone, and project them to an ansatz Gaussian functionto fix the phases required in constructing Wannier func-tions. (This procedure is the same as the projection oninitial guess functions in Ref. [41].)The obtained hopping parameters are plotted in Fig. 6.6The angle dependence without displacement field is sum-marized in Fig. 6(a). A characteristic crossover from theregime |tx| < |tx+y| to the regime |tx| > |tx+y| takes placeat around 8.5◦ as the twist angle is decreased. Within theinvestigated range of the twist angle, |t2x| and |t2x+y| aresignificantly smaller than |tx| and |tx+y|, indicating thevalidity of the model Eq. (11). The displacement field Ddependence at θ = 3◦ is shown in Fig. 6(b). As we haveseen in the simplified model, tx = ty for D = 0 due tothe emergent L symmetry, and finite D lifts this degener-acy. Interestingly, there is a point where ty = 0, i.e., al-most purely one-dimensional state. (There remains slighttwo-dimensional nature due to small but finite |tx+y| andlonger range hoppings.) In the large D side, |tx| and |ty|grows much faster than |tx+y|.The insets in Fig. 6 show the typical Wannier functionswhere the color corresponds to the sum of the wave func-tion weights for the upper and the lower layers, whichwe write ρ(r). The S4 symmetry involves the layer ex-change, but if we focus on the sum over the two layers,ρ(r) should look C4 symmetric forD = 0, which is consis-tent with the inset of Fig. 6(a). On the other hand, finiteD breaks the S4 symmetry, ρ(r) should look distorted,which is consistent with the inset of Fig. 6(b).B. Hubbard U parametersHaving the Wannier functions, the matrix elements forCoulomb interaction, in other words the Hubbard UR pa-rameters can be estimated, where R represents relativecoordinate between two Wannier functions. See for ex-ample reference [42] for the detailed method. We haveadapted the formulaUR =∫d2rρ(r +R)V (r − r′)ρ(r′)=∫d2q(2π)2ρqρ−qVqe−iq·R,(12)where ρq is a Fourier component of ρ(r), andVq =e22ϵϵ0|q|tanh(|q|ξ). (13)Here, we assume that the Coulomb interaction is screenedby top and bottom metallic gates placed with the dis-tance ξ from the sample [43–45]. From Eq. (13), we cansee that the screening effect is significant when the lengthscale 1/|q| is much larger than ξ. Note that this is an es-timation of the upper limit of U because the screeningeffect from higher energy bands in addition to the onefrom the substrate has been ignored.The twist angle dependence of U ≡ U0 without dis-placement field is plotted in Fig. 7(a). Here, we employϵ = 6, a typical value for a hBN substrate, to take into ac-count the screening effect from the substrate. The Hub-bard U is typically in the range of a few hundreds of meV,01002003000 3 6 9 12 15(a)010 1 2 3 4 5θ = 5◦10◦UR/U0|R|/aMξ = 5nmξ = 3nmU[meV]θ◦01002003000 20 40 60 80 100(b)10◦5◦2◦D [meV]FIG. 7. Hubbard U estimated using the Wannier functions.(a) Twist angle dependence for ξ = 3nm and ξ = 5nm. Insetshows |R| dependence at ξ = 3nm. (b) Displacement fielddependence for selected twisted angles with ξ = 5nm.which is significantly larger than the hopping parame-ters. As the twist angle gets smaller, U is also decreasingsince the moiré length scale becomes larger. Figure 7(b)is for the displacement field dependence, showing thatU only weakly or moderately depends on D. The insetof Fig. 7(a) is for |R| dependence of UR with ξ = 3nm,showing rapid decay of UR. The decay rate is determinedby the screening effect from the metallic gates, and de-pends on the twist angle, since the length scale is differentfor different twist angles.In order to have an intuition on the size of U , we es-timate the Coulomb interaction via U(ϵ) = e24πϵϵ0l, wherel is the size of the Wannier orbital. The latter can beobtained from the harmonic oscillator approximation ofthe leading moiré potential α1 :ℏ22m∗1l2=12α1π2 l2a2M. (14)Inserting realistic numbers m∗ = 0.64m0, a0 = 5.725Å,we can solve l = 1.6515 nm at θ = 3◦ and D = 0. Ifthe relative permittivity is taken to be ϵ = 6, it leadsto U = 0.1453 eV, which is significantly larger than thebandwidth of 0.2984 meV for the same parameter choice.If instead one plugs in the l found from Wannier func-tion calculations described above, at θ = 3◦ we wouldhave l = 2.9562 nm and U = 0.0812 eV. These estimatesmatch with the orders of magnitude in first-principle cal-culation, but the detailed discrepancies could be due tothe oversimplified treatment where we only take into ac-count α1 and the spherical approximation of the poten-tial.V. DISCUSSIONWe proposed twisted homobilayer C568 to be a promis-ing device in realizing the (extended) single-band square-7lattice Hubbard model. If successfully fabricated, the de-vice could shed important light on the mechanism of hightemperature superconductivity realized in other squarelattice materials. At negligible displacement field, thesystem is naturally in the strong coupling limit, exhibit-ing a remarkably high U/t ratio. For instance, whenthe twisting angle is at 3◦ and 9◦, the U/t ratio is ofthe order 105 and 102, respectively. One can then probethe intermediate temperature regime, of t ≪ T ≪ U ,which is hard to achieve in electronic materials. Tak-ing into account additional screening effects from higherenergy bands, we expect these U/t values to decrease,potentially reaching a minimum value of U/t ∼ 10, aregime closely resembling cuprate physics. Additionally,the bilayer offers great flexibility: by tuning the twist-ing angle, a wide range of ratios between next-nearest-neighbor and nearest-neighbor hoppings can be achieved.Furthermore, the introduction of a displacement field in-duces anisotropy between hoppings in different directionstx, ty. The orientation of this anisotropy can be rotatedby 90◦ simply by reversing the sign of the vertical electricfield. Notably, within an experimentally accessible range,the hopping amplitudes can change sign upon increasing|D|, with ty (tx) crossing from positive to negative whenD > 0 (D < 0).In addition to reproducing the target single-bandsquare lattice Hubbard model, twisted homobilayer C568can host interesting new phenomena. One possibility istunable metal-insulator transition. Fig. 8 plots the band-width of the top moiré band. Away from band crossing,the bandwidth reflects the magnitudes of the hoppingamplitudes in the system. Comparing with the onsiteCoulomb energy shown in Fig. 7, we see that at largertwisting angles and larger displacement fields (such thatthe system is highly anisotropic), the Coulomb interac-tion can have the same order of magnitude as that ofthe bandwidth. Consequently, there can potentially bea metal-insulator transition tuned by displacement field.However, we have also shown in Fig. 5 that the toptwo moiré bands will cross at larger displacement fields|D| ∼ 170 meV. Therefore, a full description of the poten-tial transition requires a two-band tight-binding model.Since the large parameter space is beyond the scope ofthe current paper, we leave a careful examination for fu-ture work.Another interesting variation of the model is the fol-lowing. In this paper, we focus on the stacking of twoidentical copies of a monolayer. In the case of the othertype of stacking, i.e., stacking between inversion images,both of Eqs. (3) and (4) reduce toVpot(τ ) = σ1Vpot(τ )σ1, (15)which kills V2,3 and induces no constraint on V0,1. Whenthere is no constraint, there is no reason for low energyeffective model to have the symmetry of the square lat-tice. Namely, we expect |tx| ≠ |ty| even in the absenceof a displacement fields. Thus, the successful fabricationof this family of structures would both shed importantFIG. 8. Left: Band width W as a function of twisting anglewhen displacement field D = 0. Right: Band width as afunction of D for different twisting angles.light on old problems and stimulate a host of new inves-tigations.Finally, we briefly discuss the potential for experi-mental realization of this proposed twisted homobilayersystem. Theoretical predictions [14] indicate that C568not only exhibits stability at room temperature but alsopossesses a formation energy of 423 meV/atom, whichis considerably lower than that of other prominent 2Dcarbon allotropes such as graphdyine (823 meV/atom),graphyne (675 meV/atom), and pentagraphene (968meV/atom). Notably, even graphdiyne with high forma-tion energy has already been experimentally synthesizedon a copper surface through a cross-coupling reaction in-volving hexaethynylbenzene [46]. This precedent makesit highly reasonable to anticipate the successful synthe-sis of C568. Moreover, in light of the promising potentialapplications in high-performance nanodevices and as an-odes in lithium-ion batteries (see, for instance, [17]), theauthors maintain an optimistic outlook regarding the ex-perimental fabrication of C568.ACKNOWLEDGEMENTSZ.-X. L. thanks Trithep Devakul, Eslam Khalaf, PhilipKim and Jedediah Pixley for enlightening conversationsand Pavel Volkov and Paul Eugenio for collaborations ona related project. This research is partially supportedby the Simons Collaboration on Ultra Quantum Matterwhich is a grant from the Simons Foundation (651440, Z.-X. L.), and JSPS KAKENHI Grant Number JP24K06968(T.K.). This work was performed in part at the AspenCenter for Physics, which is supported by National Sci-ence Foundation grant PHY-2210452. The part of cal-culations in this study have been done using the Numer-ical Materials Simulator at NIMS, and the facilities ofthe Supercomputer Center, the Institute for Solid StatePhysics, the University of Tokyo.8Appendix A: Density functional theoryIn this study, the density functional theory calcula-tions have been done using Quantum Espresso package[47, 48]. For the exchange-correlation functional, we haveadopted rev-vdW-DF2 functional [49] to take van derWaals forces into account. Pseudopotentials are fromPSlibrary [50, 51], and we use projector augmented wavemethod to include core state information in pseudopo-tentials. In order to handle our target two-dimensionalmaterial, a 30Å thick supercell in z-direction is used. Theparameters for our continuum model are derived witha naive supercell method, but we have confirmed thatdipole field effects are small and our conclusions are notaffected.Appendix B: Lower harmonics expansion and theleast square fitLet us assume τ dependence of some physical quantityQ(τ ) can be written asQ(τ ) =N∑l=1clgl(τ ) (B1)where each gl(τ ) function is some combination of har-monic functions that respects symmetry and periodicityof the system. Here, N is the number of coefficients, orin other words, the number of degrees of freedom. Onthe other hand, let us assume that we have data (fromDFT) for that quantity D(τ ) on some finite grid in τspace. Basically, cl is simply fixed by minimizingΞ =∑τ∈grid|Q(τ )−D(τ )|2. (B2)However, we sometimes want to apply constraints. Forinstance, we may take a special care about the values atτ with high symmetry, e.g., τ at zone corners or zoneboundaries. Then, Ξ is minimized under the conditionthat Q(τ ) values at selected τ points are pinned down toD(τ ). Let Nc be the number of constraints. Then, Ncshould not exceed N . For convenience, Nc constrainedgrids are labeled by i. We first define Nc ×Nc-matrix Â,Nc × (N −Nc)-matrix B̂, and Nc-dimensional vector DasAi,l = gl(τi), (l ≤ Nc) (B3)Bi,l−Nc= gl(τi), (l > Nc) (B4)Di = D(τi). (B5)Using these entities, Nc-dimensional vector C(0) and(N −Nc)×Nc-matrix C(1) are defined asC(0) = Â−1D, (B6)Ĉ(1) = t(Â−1B̂). (B7)Now, new functions g̃l(τ ) for l > Nc are introduced asg̃l(τ ) = gl(τ )−Nc∑l′=1C(1)l−Nc,l′gl′(τ ), (B8)and the constrained function Q(τ ) is written asQ(τ ) =Nc∑l=1C(0)l gl(τ ) +N∑l>Ncclg̃l(τ ). (B9)Note that we have g̃l(τi) = 0, which is confirmed asg̃l(τi) = Bi,l−Nc −Nc∑l′=1(Â−1B̂)l′,l−NcAi,l′= Bi,l−Nc − (ÂÂ−1B̂)i,l−Nc = 0,(B10)and we have Q(τi) = D(τi), which is confirmed asQ(τi) =Nc∑l=1(Â−1D)lAi,l= (ÂÂ−1D)i = D(τi).(B11)To minimize Ξ, we take derivative of Ξ with respect tocl, leading to∂Ξ∂cl= 2∑l′Ml,l′cl′ − 2Xl, (B12)whereMl,l′ =∑τ∈gridg̃l(τ )g̃l′(τ ), (B13)Xl =∑τ∈gridg̃l(τ )(D(τ )−∑l′C(0)l′ gl′(τ )). (B14)Then, the optimized cl can be obtained fromc = M̂−1X. (B15)Appendix C: Interlayer distanceTo obtain E±(kM , τ ) within DFT, we need to estimatethe interlayer distance dz for each τ . We first derive dzdependence of the total energy within DFT, and choosedz that realizes minimum total energy. While scanningover dz, the intralayer coordinates are fixed. Note that ifin-plane relaxation is also allowed, τ can flow into the onewith lower total energy, which prevents us from getting τdependence. In practice, dz-scan is done over finite set ofdz values. Then, for each τ , the obtained data are usedto fit an approximated expressionEtot(dz) = αe−β(dz−d0) − γ(d0dz)δ. (C1)90.00.51.0 0.00.51.03.603.804.004.204.40fittingDFTτx/a0τy/a0dz[Å]FIG. 9. Interlayer distance obtained by the lower harmonicsexpansion (lines) and DFT (points).Here, d0 is chosen to be 4Å, and {α, β, γ, δ} are the ad-justable parameters in the fitting. The fitting is donewith FindFit function in Mathematica.The τ dependence is approximated bydz(τ ) = c0 + c1(cos τ̃x + cos τ̃y)+ c2(cos(τ̃x + τ̃y)− cos(τ̃x + τ̃y))+ c3(cos 2τ̃x + cos 2τ̃y)+ c4(cos(2τ̃x + τ̃y) + cos(τ̃x + 2τ̃y)+ cos(2τ̃x − τ̃y) + cos(τ̃x − 2τ̃y)).(C2)The coefficients are fixed by the constrained least squarealgorithm above. Here, DFT calculations are done on8×8 regular grids in a unit cell, and the constraint is suchthat the fitting faithfully reproduce the DFT results atτ = (0, 0), (a0/2, 0), and (a0/2, a0/2), which correspondto the unit cell center, the center of the unit cell edge,and the unit cell corner, respectively. Using the obtainedcoefficients,c0 = 3.896, c1 = −0.093, c2 = 0.007 (C3)c3 = 0.021, c4 = −0.028, (C4)which are in the unit of Å, the fitting and the DFT resultscompare as in Fig. 9, showing satisfactory match betweenthem.Appendix D: Emergent symmetry LIn the convenient basis described near equation (8),the Hamiltonian isH ′(q) =(H ′−(q) + V0 + V3 eiδk·rV1(r)e−iδk·rV1(r) H ′+(q) + V0 − V3). (D1)Under C4z rotation, different terms in the Hamiltoniantransforms asC4zH′+(q)C−14z=ℏ22m∗[−i∇+R−π/2q +R−π/2q1/2]2=ℏ22m∗[−i∇+R−π/2q − q2/2]2=eiq2·rH ′−(R−π/2q)e−iq2·r.(D2)Similarly,C4zH′−(q)C−14z=ℏ22m∗[−i∇+R−π/2q +R−π/2q2/2]2=ℏ22m∗[−i∇+R−π/2q + q1/2]2= H ′+(R−π/2q).(D3)As for the moiré potentials,V0,1(Rπ/2r) = V0,1(r), V3(Rπ/2r) = −V3(r). (D4)Putting everything together, the whole Hamiltoniantransforms asC4zH̃(q)C−14z = V0σ0 − V3σ3+(H ′+(R−π/2q) eiδk·Rπ/2rV1e−iδk·Rπ/2rV1 eiq2·rH ′−(R−π/2q)e−iq2·r))(D5)The second line can be further simplified to(H ′+(R−π/2q) eiR−π/2δk·rV1e−iR−π/2δk·rV1 eiq2·rH ′−(R−π/2q)e−iq2·r))=(H ′+(R−π/2q) e−iδk·re−iq2·rV1e+iδk·reiq2·rV1 eiq2·rH ′−(R−π/2q)e−iq2·r)) (D6)where we have usedR−π/2δk = (−q1 − q2)/2 = −δk − q2. (D7)Further, applying e−i2q2·r(σ3+σ0)σ1, we obtainLH̃(q)L−1 =e−i2q2·r(σ3+σ0)σ1C4zH̃(q)C−14z σ1ei2q2·r(σ3+σ0)=H̃(R−π/2q).(D8)[1] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, andJ. Zaanen, From quantum matter to high-temperaturesuperconductivity in copper oxides, Nature 518, 179(2015).[2] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J.Millis, J. Hone, C. R. Dean, D. Basov, A. N. Pa-supathy, and A. Rubio, Moiré heterostructures as acondensed-matter quantum simulator, Nature Physics17, 155 (2021).10[3] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Hub-bard model physics in transition metal dichalcogenidemoiré bands, Phys. Rev. Lett. 121, 026402 (2018).[4] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak,K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan,et al., Simulation of hubbard model physics in wse2/ws2moiré superlattices, Nature 579, 353 (2020).[5] C. Proust and L. Taillefer, The remarkable underlyingground states of cuprate superconductors, Annual Re-view of Condensed Matter Physics 10, 409 (2019).[6] Q. Xu, N. Tancogne-Dejean, E. V. Boström, D. M.Kennes, M. Claassen, A. Rubio, and L. Xian, Engineering2d square lattice hubbard models in 90◦ twisted ge/snx(x=s, se) moiré supperlattices (2024), arXiv:2406.05626[cond-mat.str-el].[7] T. Kariyado and A. Vishwanath, Flat band in twistedbilayer Bravais lattices, Phys. Rev. Research 1, 033076(2019).[8] P. M. Eugenio, Z.-X. Luo, A. Vishwanath, and P. A.Volkov, Tunable t − t′ − U hubbard models in twistedsquare homobilayers (2024), arXiv:2406.02448 [cond-mat.str-el].[9] D. Ji, S. Cai, T. R. Paudel, H. Sun, C. Zhang, L. Han,Y. Wei, Y. Zang, M. Gu, Y. Zhang, et al., Freestandingcrystalline oxide perovskites down to the monolayer limit,Nature 570, 87 (2019).[10] K. Shigekawa, K. Nakayama, M. Kuno, G. N. Phan,K. Owada, K. Sugawara, T. Takahashi, and T. Sato, Di-chotomy of superconductivity between monolayer fes andfese, Proceedings of the National Academy of Sciences116, 24470 (2019).[11] Y. Yu, L. Ma, P. Cai, R. Zhong, C. Ye, J. Shen, G. D.Gu, X. H. Chen, and Y. Zhang, High-temperature su-perconductivity in monolayer bi2sr2cacu2o8+ δ, Nature575, 156 (2019).[12] Q.-Y. Wang, Z. Li, W.-H. Zhang, Z.-C. Zhang, J.-S.Zhang, W. Li, H. Ding, Y.-B. Ou, P. Deng, K. Chang,et al., Interface-induced high-temperature superconduc-tivity in single unit-cell fese films on SrTiO3, ChinesePhysics Letters 29, 037402 (2012).[13] N. Mounet, M. Gibertini, P. Schwaller, D. Campi,A. Merkys, A. Marrazzo, T. Sohier, I. E. Castelli, A. Ce-pellotti, G. Pizzi, et al., Two-dimensional materials fromhigh-throughput computational exfoliation of experimen-tally known compounds, Nature nanotechnology 13, 246(2018).[14] B. Ram and H. Mizuseki, C568: A new two-dimensionalsp2 − sp3 hybridized allotrope of carbon, Carbon 158,827 (2020).[15] Q. Gao, H. Sahin, and J. Kang, Strain tunable bandstructure of a new 2d carbon allotrope C568, Journal ofSemiconductors 41, 082005 (2020).[16] K. Lu, T. Wang, X. Li, L. Dai, J. Yin, and Y. Ding,A new 2d carbon allotrope C568 as a high-capacityelectrode material for lithium-ion batteries, Fullerenes,Nanotubes and Carbon Nanostructures 30, 385 (2022),https://doi.org/10.1080/1536383X.2021.1944119.[17] W.-H. Zhao, F.-Y. Li, H.-X. Zhang, R. I. Eglitis, J. Wang,and R. Jia, Doping at sp3-site in me-graphene (c568)for new anodes in rechargeable li-ion battery, AppliedSurface Science 607, 154895 (2023).[18] H.-C. Jiang and T. P. Devereaux, Superconductiv-ity in the doped hubbard model and its interplaywith next-nearest hopping t′, Science 365, 1424 (2019),https://www.science.org/doi/pdf/10.1126/science.aal5304.[19] C.-M. Chung, M. Qin, S. Zhang, U. Schollwöck, and S. R.White (The Simons Collaboration on the Many-ElectronProblem), Plaquette versus ordinary d-wave pairing inthe t′-hubbard model on a width-4 cylinder, Phys. Rev.B 102, 041106 (2020).[20] Y.-F. Jiang, J. Zaanen, T. P. Devereaux, and H.-C.Jiang, Ground state phase diagram of the doped hub-bard model on the four-leg cylinder, Phys. Rev. Res. 2,033073 (2020).[21] M. Danilov, E. G. van Loon, S. Brener, S. Iskakov, M. I.Katsnelson, and A. I. Lichtenstein, Degenerate plaque-tte physics as key ingredient of high-temperature super-conductivity in cuprates, npj Quantum Materials 7, 50(2022).[22] X. Lu, F. Chen, W. Zhu, D. N. Sheng, and S.-S. Gong,Emergent superconductivity and competing charge or-ders in hole-doped square-lattice t−J model, Phys. Rev.Lett. 132, 066002 (2024).[23] S. Jiang, D. J. Scalapino, and S. R. White, Ground-statephase diagram of the t− t′−J model, Proceedings of theNational Academy of Sciences 118, e2109978118 (2021),https://www.pnas.org/doi/pdf/10.1073/pnas.2109978118.[24] F. Chen, F. D. M. Haldane, and D. N. Sheng, Globalphase diagram of d-wave superconductivity in thesquare-lattice t − J model, Proceedings of the Na-tional Academy of Sciences 122, e2420963122 (2025),https://www.pnas.org/doi/pdf/10.1073/pnas.2420963122.[25] H.-C. Jiang, H. Yao, and L. Balents, Spin liquid groundstate of the spin- 12square J1-J2 heisenberg model, Phys.Rev. B 86, 024424 (2012).[26] Y. Nomura and M. Imada, Dirac-type nodal spin liq-uid revealed by refined quantum many-body solver usingneural-network wave function, correlation ratio, and levelspectroscopy, Phys. Rev. X 11, 031034 (2021).[27] X. Qian and M. Qin, Absence of spin liquid phase in theJ1 − J2 heisenberg model on the square lattice, Phys.Rev. B 109, L161103 (2024).[28] A. Rückriegel, D. Tarasevych, and P. Kopietz, Phase di-agram of the J1−J2 quantum heisenberg model for arbi-trary spin, Phys. Rev. B 109, 184410 (2024).[29] T. Kariyado, Twisted bilayer bc3: Valley interlockedanisotropic flat bands, Phys. Rev. B 107, 085127 (2023).[30] P. A. Volkov, J. H. Wilson, K. P. Lucht, and J. H. Pixley,Magic angles and correlations in twisted nodal supercon-ductors, Phys. Rev. B 107, 174506 (2023).[31] Y. Soeda, K. Asaga, and T. Fukui, Moiré landau levelsof a C4-symmetric twisted bilayer system in the absenceof a magnetic field, Phys. Rev. B 105, 165422 (2022).[32] Z.-X. Luo, C. Xu, and C.-M. Jian, Magic continuum in atwisted bilayer square lattice with staggered flux, Phys.Rev. B 104, 035136 (2021).[33] M.-R. Li, A.-L. He, and H. Yao, Magic-angle twisted bi-layer systems with quadratic band touching: Exactly flatbands with high chern number, Phys. Rev. Res. 4, 043151(2022).[34] O. Can, T. Tummuru, R. P. Day, I. Elfimov, A. Damas-celli, and M. Franz, High-temperature topological super-conductivity in twisted double-layer copper oxides, Na-ture Physics 17, 519 (2021).[35] S. Y. F. Zhao, X. Cui, P. A. Volkov, H. Yoo, S. Lee, J. A.Gardener, A. J. Akey, R. Engelke, Y. Ronen, R. Zhong,G. Gu, S. Plugge, T. Tummuru, M. Kim, M. Franz,https://doi.org/10.1103/PhysRevLett.121.026402https://doi.org/https://doi.org/10.1146/annurev-conmatphys-031218-013210https://doi.org/https://doi.org/10.1146/annurev-conmatphys-031218-013210https://arxiv.org/abs/2406.05626https://arxiv.org/abs/2406.05626https://arxiv.org/abs/2406.05626https://arxiv.org/abs/2406.05626https://arxiv.org/abs/2406.05626https://doi.org/10.1103/PhysRevResearch.1.033076https://doi.org/10.1103/PhysRevResearch.1.033076https://arxiv.org/abs/2406.02448https://arxiv.org/abs/2406.02448https://arxiv.org/abs/2406.02448https://arxiv.org/abs/2406.02448https://doi.org/https://doi.org/10.1016/j.carbon.2019.11.062https://doi.org/https://doi.org/10.1016/j.carbon.2019.11.062https://doi.org/10.1088/1674-4926/41/8/082005https://doi.org/10.1088/1674-4926/41/8/082005https://doi.org/10.1080/1536383X.2021.1944119https://doi.org/10.1080/1536383X.2021.1944119https://arxiv.org/abs/https://doi.org/10.1080/1536383X.2021.1944119https://doi.org/https://doi.org/10.1016/j.apsusc.2022.154895https://doi.org/https://doi.org/10.1016/j.apsusc.2022.154895https://doi.org/10.1126/science.aal5304https://arxiv.org/abs/https://www.science.org/doi/pdf/10.1126/science.aal5304https://doi.org/10.1103/PhysRevB.102.041106https://doi.org/10.1103/PhysRevB.102.041106https://doi.org/10.1103/PhysRevResearch.2.033073https://doi.org/10.1103/PhysRevResearch.2.033073https://doi.org/10.1103/PhysRevLett.132.066002https://doi.org/10.1103/PhysRevLett.132.066002https://doi.org/10.1073/pnas.2109978118https://doi.org/10.1073/pnas.2109978118https://arxiv.org/abs/https://www.pnas.org/doi/pdf/10.1073/pnas.2109978118https://doi.org/10.1073/pnas.2420963122https://doi.org/10.1073/pnas.2420963122https://arxiv.org/abs/https://www.pnas.org/doi/pdf/10.1073/pnas.2420963122https://doi.org/10.1103/PhysRevB.86.024424https://doi.org/10.1103/PhysRevB.86.024424https://doi.org/10.1103/PhysRevX.11.031034https://doi.org/10.1103/PhysRevB.109.L161103https://doi.org/10.1103/PhysRevB.109.L161103https://doi.org/10.1103/PhysRevB.109.184410https://doi.org/10.1103/PhysRevB.107.085127https://doi.org/10.1103/PhysRevB.107.174506https://doi.org/10.1103/PhysRevB.105.165422https://doi.org/10.1103/PhysRevB.104.035136https://doi.org/10.1103/PhysRevB.104.035136https://doi.org/10.1103/PhysRevResearch.4.043151https://doi.org/10.1103/PhysRevResearch.4.04315111J. H. Pixley, N. Poccia, and P. Kim, Time-reversalsymmetry breaking superconductivity between twistedcuprate superconductors, Science 382, 1422 (2023),https://www.science.org/doi/pdf/10.1126/science.abl8371.[36] X.-Y. Song, Y.-H. Zhang, and A. Vishwanath, Doping amoiré mott insulator: A t − J model study of twistedcuprates, Phys. Rev. B 105, L201102 (2022).[37] R. Haenel, T. Tummuru, and M. Franz, Incoherenttunneling and topological superconductivity in twistedcuprate bilayers, Phys. Rev. B 106, 104505 (2022).[38] P. A. Volkov, J. H. Wilson, K. P. Lucht, and J. H. Pix-ley, Current- and field-induced topology in twisted nodalsuperconductors, Phys. Rev. Lett. 130, 186001 (2023).[39] P. M. Eugenio and O. Vafek, Twisted-bilayer FeSe andthe Fe-based superlattices, SciPost Phys. 15, 081 (2023).[40] R. Akashi, Y. Iida, K. Yamamoto, and K. Yoshizawa,Interference of the Bloch phase in layered materials withstacking shifts, Phys. Rev. B 95, 245401 (2017).[41] N. Marzari and D. Vanderbilt, Maximally localized gen-eralized Wannier functions for composite energy bands,Phys. Rev. B 56, 12847 (1997).[42] V. V. Mazurenko, S. L. Skornyakov, A. V. Kozhevnikov,F. Mila, and V. I. Anisimov, Wannier functions and ex-change integrals: The example of Licu2o2, Phys. Rev. B75, 224408 (2007).[43] T. Cea, N. R. Walet, and F. Guinea, Electronic bandstructure and pinning of Fermi energy to Van Hove sin-gularities in twisted bilayer graphene: A self-consistentapproach, Phys. Rev. B 100, 205113 (2019).[44] J. Kang and O. Vafek, Non-Abelian Dirac node braidingand near-degeneracy of correlated phases at odd inte-ger filling in magic-angle twisted bilayer graphene, Phys.Rev. B 102, 035161 (2020).[45] T. Cea and F. Guinea, Band structure and insulatingstates driven by Coulomb interaction in twisted bilayergraphene, Phys. Rev. B 102, 045107 (2020).[46] G. Li, Y. Li, H. Liu, Y. Guo, Y. Li, and D. Zhu, Archi-tecture of graphdiyne nanoscale films, Chem. Commun.46, 3256 (2010).[47] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ-cioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris,G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis,A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello,L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentz-covitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of ma-terials, J. Phys. Condens. Matter 21, 395502 (2009).[48] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau,M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni,D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo,A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio,A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer,U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawa-mura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri,M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé,D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitso-nen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari,N. Vast, X. Wu, and S. Baroni, Advanced capabilities formaterials modelling with Quantum ESPRESSO, J. Phys.Condens. Matter 29, 465901 (2017).[49] I. Hamada, van der Waals density functional made accu-rate, Phys. Rev. B 89, 121103(R) (2014).[50] A. Dal Corso, Pseudopotentials periodic table: From Hto Pu, Comput. Mater. Sci. 95, 337 (2014).[51] https://dalcorso.github.io/pslibrary/.https://doi.org/10.1126/science.abl8371https://arxiv.org/abs/https://www.science.org/doi/pdf/10.1126/science.abl8371https://doi.org/10.1103/PhysRevB.105.L201102https://doi.org/10.1103/PhysRevB.106.104505https://doi.org/10.1103/PhysRevLett.130.186001https://doi.org/10.21468/SciPostPhys.15.3.081https://doi.org/10.1103/PhysRevB.95.245401https://doi.org/10.1103/PhysRevB.56.12847https://doi.org/10.1103/PhysRevB.75.224408https://doi.org/10.1103/PhysRevB.75.224408https://doi.org/10.1103/PhysRevB.100.205113https://doi.org/10.1103/PhysRevB.102.035161https://doi.org/10.1103/PhysRevB.102.035161https://doi.org/10.1103/PhysRevB.102.045107https://doi.org/10.1039/B922733Dhttps://doi.org/10.1039/B922733Dhttps://doi.org/10.1088/0953-8984/21/39/395502https://doi.org/10.1088/1361-648x/aa8f79https://doi.org/10.1088/1361-648x/aa8f79https://doi.org/10.1103/PhysRevB.89.121103https://doi.org/https://doi.org/10.1016/j.commatsci.2014.07.043 Single-band square lattice Hubbard model from twisted bilayer C568 Abstract Contents Introduction Modeling from first-principles Twisted bilayer model Tight binding model Discussion Acknowledgements Density functional theory Lower harmonics expansion and the least square fit Interlayer distance Emergent symmetry L References