Symmetry, maximally localized Wannier states, and low energy model for the twisted bilayer graphene narrow bands

We build symmetry adapted maximally localized Wannier states, and construct the low energy tight binding model for the four narrow bands of the twisted bilayer graphene. We do so when the twist angle is commensurate, near the `magic' value, and the narrow bands are separated from the rest of the bands by energy gaps. On each layer and sublattice, every Wannier state has three peaks near the triangular Moire lattice sites. However, each Wannier state is localized and centered around a site of the honeycomb lattice that is dual to the triangular Moire lattice. Space group and the time reversal symmetries are realized locally. The corresponding tight binding model provides a starting point for studying the correlated many-body phases.

We build symmetry adapted maximally localized Wannier states, and construct the low energy tight binding model for the four narrow bands of the twisted bilayer graphene. We do so when the twist angle is commensurate, near the 'magic' value, and the narrow bands are separated from the rest of the bands by energy gaps. On each layer and sublattice, every Wannier state has three peaks near the triangular Moire lattice sites. However, each Wannier state is localized and centered around a site of the honeycomb lattice that is dual to the triangular Moire lattice. Space group and the time reversal symmetries are realized locally. The corresponding tight binding model provides a starting point for studying the correlated many-body phases.
Introduction. The discovery of superconductivity and correlated insulator(s) in the 'magic' angle twisted bilayer graphene [1,2] has resulted in a remarkable flurry of theoretical activity [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. Central theoretical challenges are to understand the nature and the mechanism of the insulator(s) and the superconductor. Should the most prominent insulating states -which onset at two electrons/holes per triangular Moire unit cell, i.e. at quarter filling of the four narrow bands -be thought of as a largely featureless Mott state in which charge motion is arrested by the Coulomb repulsion, or is a spontaneously broken symmetry responsible for the charge gap? Is the superconductivity unconventional in that it breaks some of the lattice symmetries and perhaps originates from the electron-electron repulsion without a major role from electron-phonon interaction, or is it conventional?
In order to address the above questions, it is necessary to first construct a realistic, but simple, model of the electron motion in the narrow bands. As pointed out in Refs. [5,6], this is not an obvious task. When the twist angle is commensurate, the Moire pattern becomes periodic and leads to the triangular super-lattice, see e.g. [19]. At small twist angles, a unit cell contains a large number of carbon atoms and consequently the Moire Brillouin zone (MBz) becomes small. Indeed, the low energy band structure of the twisted bilayer graphene (tBG) differs in important aspects from that of two isolated monolayers due to the sizable interlayer tunneling. The four bands around the charge neutrality point have a strongly reduced bandwidth and Fermi velocity. When the twist angle is fine-tuned to the 'magic' values, the band-width becomes very narrow (but non-zero), the Fermi velocity at Dirac cones vanishes, and the quadratic band touching points appear at the corners of the MBz [1].
Although the local charge density at quarter filing is peaked at the triangular Moire lattice sites [1], as recognized in Refs. [5,6] the salient features of the narrow band structure cannot be recovered unless the Wannier states (WSs) are centered at the dual honeycomb sites. We prove this using different arguments below. In ad-dition, we diagonalize a microscopic tight binding model with the large number of atoms in the unit cell according to the prescription by Moon and Koshino [20]. Based on the layer and the microscopic carbon sublattice structure of the resulting Bloch states at the MBz center, we construct the initial ansatz for the localized WSs which we project onto the Hilbert space spanned by the four narrow bands [21]. By construction, our ansatz realizes the lattice and the time reversal symmetries locally, and forms a non-trivial representation of the site symmetry group. The result is then used as the initial step in the iterative procedure of Marzari and Vanderbild [21] to construct maximally localized, yet symmetry adapted [22], WSs. They are then used to construct the low energy tight binding model. Several theories have been proposed to address the insulating and superconducting phases [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. The closest to ours are Refs. [5,6]. However, there are also important differences. In the theory of Ref. [6], the valley U (1) symmetry -and its spontaneous breaking -plays an important role. Such valley symmetry, together with the product of C 2 and time reversal, is claimed to be an obstruction to building a tight binding model for the four narrow bands [6]. In our microscopic construction we only have the three-fold rotation about the axis formed by the AA stacked carbon atoms (C 3 ), the two-fold rotation about the axis perpendicular to the two atoms (C 2 ), and the time reversal symmetry (see Fig. 1(a)). We find the same group representations of the Bloch states at the high symmetry MBz points as conjectured in Ref. [5]. Although the WSs were not constructed explicitly in Ref. [5], the WS symmetry was insightfully deduced and is in agreement with our findings. The three-peak structure of the WSs which we find explicitly (see Fig. 3(c)) was also recognized in Ref. [6] and dubbed 'fidget spinner'.
Superlattice and band structure of the tBG: For a commensurate twist angle the Morie pattern can be specified by two integers (m, n), see e. g. [19]. The primitive translation vectors L 1 = ma 1 + na 2 and L 2 = −na 1 + (m + n)a 2 , where a 1 and a 2 are the primitive vectors of the single layer graphene lattice. As shown in Fig. 1, the triangular super-lattice sites are the positions of AA stacking. The point group symmetry operations form the D 3 group generated by C 3 and C 2 . This leads to nontrivial symmetry representations of the Bloch states at the high symmetry points in MBz, especially at Γ (k = 0) and K (k = 4π We calculate the band structure based the microscopic model of Ref. [20], which gives the values of the intralayer and the interlayer carbon-carbon tunneling amplitudes. Their tight binding Hamiltonian is written as where c ri and c † ri are the annihilation and creation operators of the electron at the carbon site r i . The detailed parameters are reproduced from Ref. [20] in the supplementary material (SM) for completeness. The MBz contains three high symmetry points Γ, K and K . The time reversal symmetry (TRS) transform K and K into each other and leaves Γ invariant.
As illustrated in Fig.2, this model contains four narrow bands with very small bandwidths near the charge neutrality point where the zero of energy has been defined. Depending on the value of the twist angle, these four bands may or may not be separated by an energy gap from the other bands in the spectrum. When m − n = ±1 mod 3 [19], at the K point, two bands form a Dirac cone and the remaining two bands are split by a tiny gap (< 0.01meV). These four Bloch states at K form a two-dimensional representation (E) and two onedimensional representations (A 1 and A 2 ) of the group D 3 [5], consistent with the degeneracy described above. The Bloch states at the center of the MBz, Γ, are doubly- degenerate; the energy difference between the two pairs defines the (narrow) band width. The doublets are the two-dimensional representations (E) of the group D 3 [5]. Using to represent the phase factor = exp(i2π/3) [23], we choose the two components of each doublet to transform as the eigenstates of C 3 with the eigenvalues of either or * , and label the four Bloch states at Γ as ψ Γ,E±, ±1 . Here E ± refers to the doublet with higher (lower) energy and ±1 refers to the component of the doublet which has the eigenvalue of ( * ) under C 3 . While the two components of each doublet are the eigenstates of C 3 , they transforms into each other under C 2 and the TRS. We wish to stress that there is no simple transformation which relates the two doublets at different energy i.e. ψ Γ,E± . This can be seen in Fig. 3(a) and 3(b) where |ψ Γ | 2 are plotted.
Wannier states: Our next step is to construct the localized WSs by applying the projection method [21]. For this purpose, it is necessary that the four bands are separated by a gap from all others. The experiments of Ref. [1,2] determined that the closest simple commensurate values are m = 30 and n = 31. However, the four bands produced by Eqn. 1 are gapped only near the band maximum, not near the band minimum; this is also seen in Ref. [1] Fig. 1. Such connection with the bands below contradicts the experimental finding that the four bands of interest are separated from either side by insulating states [1]. Therefore, we construct the WSs for the case of m = 25 and n = 26 (with the twist angle θ = 1.30 • ); the four bands are then separated by a gap on both sides. We expect that the values of the hopping parameters of the low energy Hamiltonian at the 'magic' angle to be almost the same, and, importantly, can be fine tuned to it by slight modification. We confirmed that the quadratic band touching at K, which can be taken to be the defining property of the 'magic' angle, can be realized in such a way.
As mentioned, it is crucial to identify the positions of the WSs. One naive choice is to place centers of all four states on the triangular Moire superlattice sites. With this option, WSs transform as where i, j = 1, · · · 4 are the indices of the WSs, R is the position of the triangular super-lattice site, and g is the symmetry operation. The Bloch state ψ i,k is the linear superposition of the WSs. Under the same symmetry operation g, we find It is interesting to study the special case when the momentum is symmetry invariant, i.e. Γ and K in the MBz. We immediately conclude that the Bloch states should transform as U (g), and therefore, the Bloch states should transform in the same way at Γ and K. As we pointed out, the four Bloch states transform as two doublets at Γ, and one doublet and two singlets at K. This proves that the symmetry of the Bloch states cannot be reproduced if all the WSs are placed at the sites of the triangular super lattice.
The argument above suggests that the centers of the four WSs should be placed at non-equivalent sites (Wyckoff positions) to reproduce the symmetry representations at Γ and K. A better choice is to place them at the centers of the equilateral triangles ( Fig. 1(b)), which form the dual honeycomb lattice [5,6]. Note that each triangular superlattice unit cell contains two honeycomb lattice sites. The two WSs, w 1 and w 2 , should be placed at one site, and w 3 and w 4 at another site.
To illustrate the symmetry of the Wannier states, we start by modifying the Eq. (2) for the dual honeycomb lattice [22], where R and R are still the triangular lattice translation vectors, the latter depends only on g and the WS index i. The Eq. 3 now takes the form [22] g|ψ i,k = |ψ j,gk e −igk·R U ji (g) .
Note that the extra phase factor, e −igk·R , now differentiates between Γ (where it is 1) and K (in general nontrivial). For g = C 3 and k = Γ the matrix U must be diagonal i.e. all four WSs must be eigenstates of the C 3 followed by a lattice translation, with the same eigenvalues as those of |ψ i,Γ . We therefore choose the w 1,4 and w 2,3 to have the eigenvalues and * , respectively. Next, because C 2 interchanges the two non-equivalent Wyckoff positions and the C 3 eigenvalues, we can set C 2 w 1 = w 3 and C 2 w 2 = w 4 , see Fig. 1(b). Finally, the time reversal symmetry does not change the position of the WSs, but it does conjugate the eigenvalue of C 3 . Therefore, T w 1 = w 2 and T w 3 = w 4 . These transformation rules together with translation symmetry enforce the symmetry of any low energy model.
As shown in Figs. 3(a)-3(b), we found that the magnitudes of the Bloch states at Γ display a smooth structure in real space when separated out by the layer and the microscopic carbon sublattice. This observation, along with the above considerations, suggests that a good initial ansatz for w 1 can be constructed as follows: first imagine placing a Gaussian-like cutoff centered at the first dual honeycomb site on ψ Γ,E+, , but only on the top layer and sublattice A, and the bottom layer and sublattice B. The amplitudes at the top layer and the sublattice B, and the bottom layer and the sublattice A, are taken from the similarly cut-off ψ Γ,E−, . This guarantees good overlap with the Bloch states. C 2 now generates w 3 , and the TRS generates w 2 ; when the TRS is applied to w 3 it finally gives w 4 .
In the next step, we use the initial ansatz as an input to the Wannier90 program [24] with site symmetry enforced, on the 30 × 30 k-space mesh, and after 200 iterations obtain the four maximally localized WSs. Fig. 3(c) shows the shape of the resulting |w 1 | 2 on different layers and different sublattices. As seen, w 1 is well localized and centered around the dual honeycomb lattice site; it also displays three different peaks, located around the triangular lattice sites. This is consistent with the local density of states obtained by DFT calculations which are also peaked around the triangular lattice sites. We checked that the remaining WSs obtained in this way are related by the mentioned symmetry: w 2 = w * 1 , w 3 = C 2 w 1 and w 4 = w * 3 .
Tight binding model : The tight binding model based on the maximally localized WSs can be readily con-structed. The on-site term must be of the form where f † i,R and f i,R are the creation and annihilation operators of the WSs w i,R . This is because C 3 prohibits mixing between w 1 and w 2 , and w 3 and w 4 ; C 2 and TRS then force a single real parameter µ. In contrast, there are two such parameters in Ref. [6]. The nearest neighbor hopping term is between w 1,2;Ri and the neighboring w 3,4;Rj . The most general term allowed by symmetry is where t 1 is real and t 1 are, in general, a complex number (see SM). The structure of WSs, seen in Fig. 3(c), suggests that the overlap of the next and even the next-next nearest neighboring WSs is sizable, and thus cannot be neglected even in the minimal model. The detailed analysis of the symmetry constraints on the further range hopping is now straightforward. It is presented in the SM, where we also study the quality of the fit as a function of the hopping range. The blue points in Fig. 2 shows the interpolated band structure obtained from the Wan-nier90 program using all the hoppings. The agreement with the the microscopic tight binding model based on DFT calculations (red solid lines) suggests that our low energy model accurately reproduces the main physics. Conclusion: In this paper, we presented a method for constructing symmetry adapted maximally localized Wannier functions and the corresponding low energy model for the four narrow bands of the tBG near the 'magic' angle. The WSs have three peaks around the Moire triangular lattice sites, but are centered at the dual honeycomb lattice sites. They form non-trivial representations of the site symmetry group. Our model provides a firm basis for further study of the many-body effects.
We Supplementary material for "Symmetry, maximally localized Wannier states, and low energy model for the twisted bilayer graphene narrow bands"

I. MICROSCOPIC TIGHT BINDING MODEL
We use the same model as in Ref. [1] to produce the band structure. The Hamiltonian is where we the first and the second term in the Hamiltonian H are for the intralayer and the interlayer tunneling, respectively. We set V 0 ppπ = −2.7eV, V 0 ppσ = 0.48eV. a 0 = 0.142nm is the distance between the two nearest neighbor carbon atoms on the same layer. The decay length for the hopping is δ = 0.319a 0 . The hopping with d > 4a 0 is exponentially small and thus is neglected in the model. The band structure produced by this model is illustrated as red solid lines in Fig. 2.

II. PROJECTION METHOD
In this section, we will explain the projection method we used to produce the localized WSs as the input of the Wannier90 program. We follow the approach in Ref. [2]. As explained in the text, we first construct the trial functions |f i (i = 1, · · · 4), which transform in the same way as the WSs. These trial states are not necessarily orthogonal or normalized. For the Bloch states |ψ i,k , we define the matrix A(k) ij = ψ i,k |f j . The states are smooth in k, because the arbitrary k-dependent phase cancels in the projector. Smoothness in k is required in order for WSs to be localized in real space. However, they are not orthonormal. To construct the orthonormalized k-smooth Bloch-like states, we define the matrix S(k) = A † (k)A(k), and In practice, we apply the singular value decomposition to the matrix A(k) = U (k)D(k)V † (k), where the matrices U (k) and V (k) are unitary, and D(k) is diagonal. It is easy to show that A(k)S −1/2 (k) = U (k)V † (k) is unitary, and thus, |ψ is orthogonal and normalized. With the projection method, the WSs are We use the above as an input to Wannier90 program with site symmetry enforced to obtain the maximally localized symmetry adapted WSs.

A. Symmetry Constraints
In this subsection, we discuss the most general form of the hopping amplitudes allowed by symmetry. As mentioned in the main text, the symmetries are C 3 , C 2 , and the TRS. Since we have already explained the constraints on the on-site hybridization term in the main text, we first study the next term which is the hopping between w 1,2;Ri and the neighboring w 3,4;Rj . Thus, the Hamiltonian is of the form It should be invariant under all the symmetry transformations. First, consider C 3 which brings WSs w R into w R in a different unit cell. In addition, w 1 and w 4 have the eigenvalue of , and w 2 and w 3 have the eigenvalue of * . The C 3 invariance of the Hamiltonian forces t n 14 = t n 14 = t n 14 , t n 23 = t n 23 = t n 23 , t n 13 = t n 13 = * t n 13 , and t n 24 = * t n 24 = t n 24 .
C 2 transforms w 1 ↔ w 3 and w 2 ↔ w 4 , and brings L 1 → L 2 − L 1 , and L 2 → L 2 . Combined with the hermiticity of the Hamiltonian, the C 2 invariance leads to t n 14 = (t n 23 ) * , t n 13 = (t n 13 ) * and t n 24 = (t n 24 ) * . Finally, the TRS enforces t n 13 = (t n 24 ) * . Combining all constraints, we set t n 13 = t 1 and t n 14 = t 1 , where t 1 is real and t 1 is in general, a complex number. Thus, the nearest neighbor hopping term can be written as It seems that t 1 is a complex number. If we apply a gauge transformation w 1,3 → e iθ w 1,3 and w 2,4 → e −iθ w 2,4 , the hopping constant t 1 is invariant but t 1 → e 2iθ t 1 . Thus, the phase of t 1 can be always removed by choosing a particular gauge of the WSs. Therefore, there are only two free parameters for the nearest neighbor hopping [3]. Next, consider the next nearest neighbor hopping H nn : Let us first consider the symmetry constraints on t nn ij when i, j = 1, 2. The hopping constants t nn ij (i, j = 3 or 4) can be obtained by applying C 2 symmetry operation. Therefore, the next-nearest neighbor hopping can be described by two complex numbers t 2 = t nn,1 The constraints are very similar to the one for the nearest neighbor hopping. We found C 3 enforces 0.0012 + 0.0099i Table S1. The hopping constants between w1,2 and w3,4. All the numbers are in the units of meV.  Table S2. The hopping constants between w1,2. All the numbers are in the units of meV.