Maximally-localized Wannier orbitals and the extended Hubbard model for the twisted bilayer graphene

We develop an effective extended Hubbard model to describe the low-energy electronic properties of the twisted bilayer graphene. By using the Bloch states in the effective continuum model and with the aid of the maximally localized algorithm, we construct the Wannier orbitals and obtain an effective tight-binding model on the emergent honeycomb lattice. We found the Wannier state takes a peculiar three-peak form in which the amplitude maxima are located at the triangle corners surrounding the center. We estimate the direct Coulomb interaction and the exchange interaction between the Wannier states. At the filling of two electrons per super cell, in particular, we find an unexpected coincidence in the direct Coulomb energy between a charge-ordered state and a homogeneous state, which would possibly lead to an unconventional many-body state.

It is a challenging problem to theoretically describe the many-body physics in TBG. At θ = 1.05 • , the spatial period of the moiré pattern is more than 10 nm and the number of carbon atoms in a unit cell exceeds 10,000. The electronic property of such a huge and complex system can be calculated efficiently by the effective continuum model which captures the long-wavelength physics associated with the moiré period. 23,28,29,31,[34][35][36] However, the effective continuum energy spectrum still contains a number of energy bands in the low-energy region, and we need one more step to simplify the model so as to exclusively describe the nearly-flat bands at lowest energy.
Actually the nearly-flat bands are separated by the energy gaps from other bands, 1,37,38 making it possible to construct an effective lattice model with well-localized Wannier orbitals purely consisting of the flat band states. Such an effective model was actually predicted by the symmetry analysis, 4 which concludes that the Wannier orbitals should be centered at non-equivalent AB spot and BA spot in the moiré pattern, to form an emergent honeycomb lattice. Arguments and calculations suggesting a honeycomb lattice description have also been put forward in an independent work. 5 To obtain a concrete model with specific parameters, we need to construct the actual Wannier orbitals from the realistic model of TBG.
In this paper, we develop an extended Hubbard model of TBG at the magic-angle (θ = 1.05 • ), based on the effective continuum model. By taking an appropriate linear combination of the Bloch states in the nearly flat bands, we build the Wannier orbitals centered at AB and BA spots, and obtain the effective tight-binding model on the emergent honeycomb lattice. Here we adopted the maximally localized algorithm 39 to minimize the spread of the wave functions. From the obtained Wannier orbitals, we estimate the direct Coulomb energy and the exchange energy between electrons residing at different orbitals. The obtained Wannier state is centered at AB or BA spot, while its maximum amplitude is found to be not at the center, but at three AA spots surrounding the center, as also noticed in Ref. 5. Importantly, the pair of Wannier orbitals that we constructed explicitly has (p x , p y ) on-site symmetry, hence forms a doublet under three-fold rotation around their centers, consistent with the symmetry analysis. 4 Due to this peculiar three-peak form, the electronelectron interaction between the neighboring sites is as important as the on-site interaction term. At the filling of two electrons per super cell, in particular, we find an unexpected coincidence in the direct Coulomb energy between two different many-body states: a homogeneous state where an electron enters every sublattice of the effective honeycomb lattice, and a charge-ordered state where two electrons reside at every two sublattices [ Fig. 7(a) and (b)]. We expect that such competing nature would possibly give rise to a nontrivial many-body ground state. This paper is organized as follows: In Sec. II, we explain the atomic structure of TBG, and in Sec. III, we introduce the effective continuum model and argue the structure of the nearly-flat bands at the magic angle θ = 1.05 • . In Sec. IV, we construct the Wannier orbitals using the maximally localizing method, and obtain the tight-binding model in Sec. V. We consider the electronelectron interaction between the Wannier states in Sec. VI. A brief conclusion is presented in Sec. VII.

II. ATOMIC STRUCTURE
We define the atomic structure of TBG by starting from AA-stacked bilayer graphene (i.e. perfectly overlapping honeycomb lattices) and rotating the layer 1 and 2 around a pair of registered B-sites by −θ/2 and +θ/2, respectively. We define a 1 = a(1, 0) and a 2 = a(1/2, √ 3/2) as the lattice vectors of the initial AA-stacked bilayer before the rotation, where a ≈ 0.246 nm is the lattice constant of graphene. The corresponding reciprocal lattice vectors are a * 1 = (2π/a)(1, −1/ √ 3) and a * 2 = (2π/a)(0, 2/ √ 3). After the rotation, the lattice vectors of layer l are given by a (l) i = R(∓θ/2)a i with ∓ for l = 1, 2, respectively, where R(θ) represents the rotation by θ. Likewise, the reciprocal lattice vectors become a * (l) i = R(∓θ/2)a * i . With respect to the registered Bsites, TBG has point group D 3 generated by a three-fold in-plane rotation C 3z along z-axis and a two-fold rotation C 2y along y-axis.
In a small angle TBG, the slight mismatch of the lattice periods of two layers gives rise to a long-period moiré interference pattern. The reciprocal lattice vectors for the moiré pattern is given by Figure 1(a) illustrates the atomic structure of TBG with θ = 3.89 • . The lattice structure locally resembles the regular stacking such as AA, AB or BA depending on the position, where AA represents the perfect overlapping of hexagons, and AB (BA) is the shifted configuration in which A 1 (B 1 ) sublattice is right above B 2 (A 2 ). In Fig. 1(a), AA spots are located at the crossing points of the grid lines, and AB and BA spots are at the centers of triangles indicated by dots. Figure 1(b) shows the corresponding folding of the Brillouin zone, where two large hexagons represent the first Brillouin zones of layer 1 and 2, and the small hexagon is the moiré Brillouin zone of TBG. The graphene's Dirac points (the band touching points) are located at K (l) is the valley index. We label the symmetric points of the reduced Brillouin zone asΓ,M ,K andK ′ as in Fig. 1(b).
We can construct the TBG in alternative manners, for example, by rotating around the hexagon centers instead of the B site. In that case, we have the different superlattice structure with point group D 6 . For completeness, we leave the discussion of D 6 structure and other superlattice structures in the Supplementary Material. 40

III. EFFECTIVE CONTINUUM MODEL
When the moiré period is much longer than the atomic scale, the electronic structure can be described by an effective continuum model. 23,28,29,31,[34][35][36] There the intervalley mixing between ξ = ± can be safely neglected and the total Hamiltonian is block-diagonalized into the two independent valleys. The effective Hamiltonian of continuum model for the valley ξ is written in a 4 × 4 matrix for the basis of (A 1 , B 1 , A 2 , B 2 ) as Here H l (l = 1, 2) is the intralayer Hamiltonian of layer l, which is given by the two-dimensional Weyl equation centered at K (l) ξ point, where ± is for l = 1 and 2, respectively. We take v/a = 2.1354 eV. 34 U is the effective interlayer coupling given by [34][35][36] where ω = e 2πi/3 . Here u and u ′ describe the amplitudes of diagonal and off-diagonal terms, respectively, in the sublattice space. The effective models in the previous studies 34-36 assume u = u ′ , which corresponds to a flat TBG in which the interlayer spacing d is constant everywhere. On the other hand, several theoretical studies predicted that the optimized lattice structure of TBG is actually corrugated in the out-of-plane direction, in such a way that d is the widest in AA stacking region and the narrowest AB / BA stacking region. [41][42][43][44] Here we incorporate the corrugation effect as a difference between u = 0.0797eV and u ′ = 0.0975eV in the effective model, of which detailed derivation is presented in the Appendix A. As we show in the following, the difference between u and u ′ introduces energy gaps between the lowest bands and the excited bands, in a qualitative agreement with the experimental observation. 1,2,37 It was found that the energy gaps isolating the lowest nearly-flat bands are also caused by the in-plane distortion. 38 The calculation of the energy bands and the eigenstates is done in the k-space picture. For a single Bloch vector k in the moiré Brillouin zone, the moiré interlayer coupling hybridizes the graphene's eigenstates at q = k + G, where G = m 1 G M 1 + m 2 G M 2 and m 1 and m 2 are integers. Therefore the eigenstate is written as where X = A 1 , B 1 , A 2 , B 2 is the sublattice index, n is the band index and k is the Bloch wave vector in the moiré Brillouin zone. As the low-energy states are expected to be dominated by the individual graphenes' eigenstates near the original Dirac points, we pick up q's inside the cut-off circle |q − q 0 | < q c , where q 0 is taken as the midpoint between K Since the intervalley coupling can be neglected, the calculation is done independently for each of ξ = ± as we discussed previously. We then numerically diagonalize the Hamiltonian within the limited wave space inside the cut-off circle and obtain the eigenenergies and eigenstates. Figure 2(a) shows the energy band and the density of states of TBG at the magic angle θ = 1.05 • , calculated by this approach. Here in the following, the origin of band energy axis is set to the charge neutral point. The lower panel is the enlarged plot of the zero-energy region where the near-flat bands are located. The black solid line and red dashed line represent the energy bands of ξ = ± valleys, respectively. They are the time-reversal partners to each other, and the energy bands of ξ = − are obtained just by inverting k to −k. The flat band cluster consists of two bands per spin and valley, which are denoted as E 1 (k) and E 2 (k) for the hole side and the electron side, respectively. The overall structure is about 7.5 meV wide in energy axis, and separated from the excited bands by the energy gap of about 14 meV in each of the electron side and the hole side. Figure  2(b) shows the contour plots of E 1 (k) and E 2 (k) for the valley ξ = +. E 1 (k) and E 2 (k) are trigonally warped in the opposite directions, so that E 1 (k) = E 1 (−k) and E 2 (k) = E 2 (−k). The particle-hole symmetry is absent and the E 1 band is wider than the E 2 band. The van-Hove singularity is located at E ≈ −0.11 meV and 0.16 meV, which correspond to the carrier density n/n 0 ≈ −0.78 and 0.63, respectively, with spin and valley included.
is the moiré unit area (the band gap is n/n 0 = ±4) and L M is 13.4nm at θ = 1.05 • . The filling of two electrons / holes per super cell (n/n 0 = ±2) corresponds to E ≈ 0.289 meV and −0.286 meV, respectively, which are indicated by dashed contours in Fig. 2

IV. WANNIER ORBITALS
We construct the localized Wannier orbitals from the Bloch wave functions of the effective model. Since the nearly flat bands are energetically isolated from other bands, we expect that well-localized orbits can be made purely from the flat band states with all other bands neglected. The number of the independent Wannier orbitals in a unit cell coincides with the number of the energy bands taken into account, so we have two Wannier orbitals per spin and valley. According to the symmetry analysis 4 , the two orbitals should be centered at AB and BA spots to form a honeycomb lattice. Our strategy is to first prepare certain initial orbitals centered at AB and BA, and then apply the maximally localized algorithm. 39 The following process is applied to ξ = ± valleys separately, and we omit the valley index ξ hereafter.
The initial wave functions can be prepared as follows. First we fix the global phase factor of the Bloch states in two different ways: In gauge 1, we fix the phase so that ψ B1 nk (r BA ) is real, and in gauge 2, we fix the phase so that ψ A1 nk (r AB ) is real. Here are the positions of BA and AB spots, respectively, measured from the AA spot (0, 0) [ Fig. 1(a)]. We write the Bloch function in the gauge 1 as ψ nk , and that in the gauge 2 as e iφ nk ψ nk , where e iφ nk is the relative phase factor between gauge 1 and 2. We construct the initial Wannier orbitals 1 and 2 by summing the Bloch states of the bands ψ 1k and ψ 2k (corresponding to E 1 (k) and E 2 (k), respectively) as, Here R = n 1 L M 1 + n 2 L M 2 is the moiré lattice vector, and the summation in k is taken over N discrete points in the moiré Brillouin zone. We take N = 18 × 18 in this study. It is straightforward to check the orthonormality, where mn to minimize the spread functional. We put Eq. (5) as the initial value of U (k) mn , and iterate the minimization process until the convergence. In each step, we impose the symmetry constraint to U (k) mn . The optimized Wannier orbitals for the valley ξ = + are illustrated in Fig. 3. Those for the opposite valley ξ = − are given by the complex conjugate. For each of orbital 1 and 2, the top five panels show the contour maps for the squared amplitudes of the total wave function and of the four sublattice components. We actually see that the orbital 1 and 2 are centered at BA and AB positions, respectively, while the maximum of the wave amplitudes are located not at the center, but near three AA spots surrounding the center. This reflects the fact that the Bloch wave functions of the nearly-flat bands are mostly localized AA spot of the moiré pattern. 25,41 The lower panels illustrate the phase of the envelope function F X l (r) (X = A, B and l = 1, 2) on some sample points, where the total wave function is ψ X l (r) = e iK (l) ξ ·r F X l (r). Here the absolute value of F X l (r) is indicated by the radius of a circle, and its phase factor is by the direction of a bar and also by color. Now we see that the envelope functions on different sublattices have different eigenvalues of C ′ 3z , in-plane rotation with respect to its own center. However, noting that the Bloch factor e iK ξ ·r also carries a non-zero eigenvalue of C ′ 3z , the total wave function ψ = (ψ A1 , ψ B1 , ψ A2 , ψ B2 ) is found to be an eigenstate of C ′ 3z with a single eigenvalue. In orbital 1, for example, the C ′ 3z eigenvalue of F X l is (ω, 1, 1, ω * ) for (A 1 , B 1 , A 2 , B 2 ), so that the angular momentum of the envelope function is written as L (env) z = (−1, 0, 0, 1). On the other hand, the C ′ 3z eigenvalue for the Bloch factor e iK ξ ·r can be found by noting that BA spot (the orbital center) coincides with A 1 site and the center of hexagon of layer 2 [ Fig. 1(a)], and then is −1 for all the sublattices. Similarly, we can show L z = −1 also for orbital 2. Since the Wannier functions at the opposite valleys are related by the complex conjugate, we finally conclude that the eigenvalue of C ′ 3z is ω ξ = e ξ2πi/3 for both orbital 1 and 2. Namely orbital 1 and 2 from the same valley ξ have the same nonzero angular momentum L z = −ξ, in accordance with the symmetry analysis. 4 The initial guess of the Wannier orbital in Eq. (5) is closely related to the angular momentum of the envelope function. For the orbital 1, the envelope function of B 1 has zero angular momentum, so that it has a finite amplitude at the orbital center r BA as seen in Fig. 3. It does not contradict with the nonzero total angular momentum L z = −1, because BA spot coincides with the hexagon center of layer 1, but not B 1 site. The finite amplitude at r BA is actually linked to the gauge choice for |R, 1 0 , which requires that ψ B1 nk (r BA ) is real. There all the wave functions add up in the same phase at r BA , so that we have an orbital localized at r BA with finite amplitude. The same is true for the orbital 2, of which envelope angular momentum vanishes at A 1 . The wrong gauge choices (e.g., ψ A1 nk (r BA ) is real) do not make a well localized orbital, because the angular momentum of the Wannier function is forced by the symmetry. Also, the hybridized form of |ψ 1k ± |ψ 2k in Eq. (5) better localizes the wave function than just using |ψ 1k , |ψ 2k . This is similar to monolayer graphene having the same honeycomb lattice structure, where the superposition of the positive and negative energy states is required to have A-site or B-site localized orbital.

V. EFFECTIVE TIGHT-BINDING MODEL
From the Wannier orbitals and the energy bands, we can derive the effective tight-binding model to exactly reproduce the dispersion of the nearly-flat bands. In a straightforward calculation, the hopping integral between the Wannier orbitals is written as whereÛ (k) represents the matrix U (k) mn , and we used ψ n ′ k ′ |H|ψ nk = δ nn ′ δ kk ′ E n (k). In Fig. 4(a) and (b), we plot the hopping integrals from orbital 1 and 2, respectively, for the valley ξ = +. Here the honeycomb lattice represents the network of BA spots (orbital 1) and AB spots (orbital 2). The radius of the circle at each lattice point indicates the absolute value of the hopping integral from the origin (green circle at the center) to that point, and the direction of the bar represents the phase in the complex plane. The effective tight-binding model for the valley ξ = − is just given by taking the complex conjugate. The list of the hopping integrals for the valley ξ = + is included in Supplementary Material. 45 To understand this effective tight-binding model, we need to analyze the symmetry properties of Wannier orbitals under point group D 3 , as in Ref. 4. Recall that orbitals 1 and 2 have nonzero angular momentum L z = −ξ at the valley ξ. Furthermore, under two-fold rotation C 2y which interchanges two graphene layers, we find orbital 1 from valley ξ is mapped to orbital 2 from valley −ξ and vice versa as shown in Fig. 3. Hence we can regard orbitals 1 and 2 from valley ξ as the p-wave-like orbitals p ξ ≡ p x + iξp y residing on BA and AB spots respectively. The angular momentum of p ξ orbital is L z = −ξ whether its center is at BA or AB spot, which is consistent with L z of orbital 1 and 2. Under C 2y the two graphene layers and hence BA and AB spots are interchanged, and (p x , p y ) → (−p x , p y ) or p ξ → −p −ξ . In other words, C 2y interchanges p ξ orbital at BA spot and p −ξ orbital at AB spot, which reproduces the symmetry transformation of orbital 1 and 2 under C 2y .
Once we identify the symmetries of orbital 1 and 2, the tight-binding model then describes hopping among (p x , p y ) orbitals on the honeycomb lattice formed by BA and AB spots, which reads where c iξ annihilates a p ξ -orbital electron at site i, r ij is the vector from site i to j, and t(r), φ(r) are as shown in Fig. 4 (a) and (b). The symmetry group of the tight-binding model of Eq.
(8) is G = D 3 ×U(1)×SU(2)×T , where D 3 is the point group of TBG, which acts jointly on lattice sites and (p x , p y ) orbitals, U(1) acts in orbital space, SU(2) acts in spin space and T is the time-reversal symmetry. As discussed in Ref. 4, the microscopic origin of this orbital U(1) symmetry is that at small twist angles the intervalley coupling is strongly suppressed, leading to this approximate valley conservation that exists independent of crystal symmetries. The hopping integral t(r) roughly decays with increasing r = |r|. To include dominant contributions, we consider the nearest five hopping integrals t 1 to t 5 shown in Fig. 4 (a) and (b), which are within the range r √ 3L M . Notice that the subscript are not labeled according to r. In the present model, we have t 1 ≈ 0.331 meV, t 2 ≈ (−0.010 ± 0.097i) meV, t 3 ≈ 0.016 meV, t 4 ≈ 0.036 meV, and t 5 ≈ 0.119 meV. Figure 5 presents the band structure in the effective tight-binding models with (a) t 1 and t 2 , (b) t 1 , t 2 and t 5 and (c) all the hopping parameters within the distance r < 9L M . Dashed line indicates the original energy band of the effective continuum model.
With hopping terms t 1 and t 2 only, tight-binding model (8) where c i = (c i,x , c i,y ) T with c i,x(y) annihilating an electron with p x(y) -orbital at site i, c jξ = (c jx + iξc jy )/ √ 2. µ is the on-site chemical potential,t 2 = Re(t 2 ), t ′ 2 = Im(t 2 ), and the sum over ij ′ includes bonds with length √ 3L M along three directionsx, C 3zx and C 2 3zx . The minimum tight-binding model (9) gives rise to a spectrum with Dirac nodes atK,K ′ points. Notice that t 1 denotes hopping between two sublattices and we can always make t 1 real by properly choosing the relative phase between sublattices. The t ′ 2 term describes the hexagonal warping effect in orbital space, which is responsible for band splittings alongΓM lines as shown in Fig. 5.
The symmetry group G allows finite gaps atK,K ′ points. 4 However, due to the approximate sublattice symmetry at small twist angles, 31 we can introduce an additional Z 2 symmetry g : c Rξ → c −R,−ξ , which combines twofold rotation in real space and chirality flip in orbital space. In the presence of g and original symmetry group G, the gapless Dirac nodes atK,K ′ points are guaranteed. The minimum model (9) satisfies both G and g.
With this additional Z 2 symmetry g, we then consider additional hopping terms t 3 , t 4 and t 5 . As finite Im(t 3 ) obeys G while violates g, we find Im(t 3 ) =0 from our numerical calculation of hopping integrals. The nonzerõ t 3 ≡ Re(t 3 ), t 4 , t 5 terms preserve both G and g, and quantitatively modify the band structure of the minimum model (9). In fact including t 1 to t 5 we have m −1 e,h = We can also incorporate the effect of intervalley coupling in the effective tight-binding model by introducing U(1)-breaking hopping terms such as those in Ref. 4, which may explain Landau level degeneracy lifting in experiments. A detailed analysis of Dirac nodes and mass generation will be presented in a forthcoming work.

VI. ELECTRON-ELECTRON INTERACTION
We can calculate the electron-electron interaction parameters between the Wannier orbitals directly from the wave functions obtained above. The direct Coulomb interaction V and the exchange interaction J between x y AB (orbital 2) BA (orbital 1) |R, m and |R ′ , m ′ are defined by where ǫ is the dielectric constant induced by the electrons in other bands and by the external environment (e.g., the substrate). The direct term is the classical Coulomb interaction and it works for any combinations of spin and valley. On the other hand, the exchange interaction works only for the same spin and the same valley. Rigorously speaking, the exchange term between different valleys (and the same spin) is not exactly zero, but there the integral of e i(K+−K−)·(r−r ′ ) /|r − r ′ | in Eq. (11) becomes much smaller than that for the same valley, so we neglect it. We label the direct interaction terms at different distances as V 0 , V 1 , V 2 · · · as in Fig. 6(a), where V 0 is the on-site interaction, V 1 is the nearest neighbor interaction TABLE I. Direct interaction Vn and the exchange interaction Jn for the Wannier orbitals in units of e 2 /(ǫLM). The definition of V0, V1 · · · is presented in Fig. 6(a). V (approx) n is the direct interaction terms estimated by the point-charge approximation (see the text). and so forth. Similarly the exchange terms can be labeled as J 1 , J 2 · · · , where J 0 does not exist due to the Pauli principle. The calculated interaction parameters are listed in Table I. Here we notice that the on-site interaction V 0 is not much greater than others, but it is in a similar magnitude to the nearest-neighbor interaction V 1 . The further interactions V 2 and V 3 are more than half of V 0 . This is quite different from usual Hubbard-type models where V 0 dominates the interaction effect. The peculiar distance dependence of Coulomb interaction in this model is closely related to the three-peak structure of the Wannier orbital. For a single electron, each of three peaks accommodates the electric charge of −e/3, and thus the interacting potential between two electrons can be written as a summation over the nine combinations of those fractional charges. The direct Coulomb potential between two fractional charges located at the same peak (i.e., the "on-site interaction" for the fractional charges) is u 0 ≈ (e/3) 2 /ǫ/(0.28L M ), while the potential between different peaks is well approximated by that for the point charges, i.e., (e/3) 2 /ǫ/r, where r is the distance between the peak centers. The direct interaction terms estimated by this approximation are presented as V (approx) n in Table  I, where the error is found to be 1% or less.
Obviously the dominant contribution to V n comes from the on-site part u 0 . As shown in Fig. 6(b), the electrons located at the same orbital (V 0 ) share all the three peaks, the nearest neighbor configuration (V 1 ) the two peaks, and the next nearest ones (V 2 , V 3 ) a single peak. Therefore, the on-site interaction of the fractional charges included in V 0 , V 1 , V 2 , V 3 is 3u 0 , 2u 0 , u 0 , u 0 , respectively, and this explains the dominant part of the relative amplitudes of V n 's. These relatively long-range electronelectron interactions can potentially modify the hopping parameters and hence renormalize the low-energy band structure.
At the filling of two electrons per super cell, in particular, the triangular charge distribution results in an unexpected coincidence in the direct Coulomb energy between two different many-body states shown in Fig. 7, with (a) a homogeneous state where an electron resides at every sublattice of the honeycomb lattice, and (b) a charge-ordered state where two electrons enter every two sublattice. It may seem that the direct Coulomb energy in (b) is greater than in (a) because of the double occu-pancy. However, since an electron at the honeycomb site is actually composed of three 1/3 charges at the triangle corners, the states (a) and (b) have nearly identical charge distribution as shown in the right panels, where the charge of −2e is registered to every AA spots. Considering that the direct Coulomb interaction is very well approximated by the simple point-charge model as argued above, the total direct energies of (a) and (b) must be nearly equal. The competing nature of the two completely different states may suggest a nontrivial manybody ground state. For further consideration, we need to include the exchange interaction and also the kinetic energy. Lastly, Fig. 7(c) illustrates an excitation from the state (a), where an electron is transferred from a single honeycomb site to another. This actually corresponds to a pair creation of the fractional charges (±1/3)e as shown in the right panel. This is another intriguing property at this filling factor.

VII. CONCLUSION
An extended Hubbard model is obtained for the nearly flat band in the low-angle TBG by starting from the Bloch states in a realistic continuum model. The obtained Wannier localized state is centered at AB or BA spot to form a honeycomb lattice. The wave function of the Wannier orbital takes a triangular form which peaks at three AA spots surrounding the center, and it leads to a competition between the on-site interaction and the neighboring interaction. At the filling of two electrons per super cell, in particular, we have an unusual degeneracy of the a charge-ordered state and a homogeneous state, which implies an nontrivial nature of the ground state. The detailed studies for the many-body ground states in this model will be left for future works. Note added: After completion of the present study, we have come to notice that a recent preprint which also reports the maximally-localized Wannier states for TBG. 46 We derive the effective interlayer interaction of a corrugated TBG in Eq. (3) following the method in Ref. 34. We start from the single-orbital tight-binding model for p z orbital of carbon atoms. We assume that the transfer integral between any two orbitals is written in terms of the Slater-Koster form as, Here e z is the unit vector perpendicular to the graphene plane, a 0 = a/ √ 3 ≈ 0.142 nm is the distance of neighboring A and B sites on graphene, and d 0 ≈ 0.335 nm is the interlayer spacing of graphite. The parameter V 0 ppπ is the transfer integral between the nearest-neighbor atoms on graphene, and V 0 ppσ is the transfer integral between vertically located atoms on the neighboring layers of graphite. We take V 0 ppπ ≈ −2.7 eV, V 0 ppσ ≈ 0.48 eV, to fit the dispersions of monolayer graphene. 34 Here r 0 is the decay length of the transfer integral, and is chosen as 0.184a so that the next nearest intralayer coupling becomes 0.1V 0 ppπ . To construct the Hamiltonian matrix, we define the Bloch wave bases as where the position R A l (R B l ) runs over all A(B) sites on the layer l(= 1, 2), N is the number of monolayer's unit cell in the whole system, and k is two-dimensional Bloch wave vector defined in the first Brillouin zone of monolayer on the layer l.
When we rotate one graphene layer to another by a small twist angle θ, the local lattice structure in the moiré pattern is approximately viewed as a non-rotated bilayer graphene, where the displacement δ slowly depends on the position r in accordance with 38 δ(r) = [R(θ/2) − R(−θ/2)]r.
The interlayer matrix element for valley ξ is then approximately written by U X ′ X [K ξ , δ(r)]. 34 Using Eqs. (A8) and (A10), we obtain where we used the relationship a * i · δ(r) = G M i · r. Now we see that Eq. (A11) is periodic in r with the moiré reciprocal vectors G M i . SinceŨ X ′ X (q) rapidly decays in q, we only need a few Fourier components in Eq. (A11).