Model for Metal-Insulator Transition in Graphene Superlattices and Beyond
Noah F. Q. Yuan Liang Fu








Department of Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts 02139, USA


We propose a two-orbital Hubbard model on an emergent honeycomb lattice to describe the low-energy physics of twisted bilayer graphene. Our model provides a theoretical basis for studying metal-insulator transition, Landau level degeneracy lifting and unconventional superconductivity that are recently observed.


Model for Metal-Insulator Transition in Graphene Superlattices and Beyond
    Noah F. Q. Yuan and Liang Fu
    Jan 2018
=========================================================================



The recent discovery of correlated insulator state <cit.> and unconventional superconductivity <cit.> in bilayer graphene with a small twist angle has generated tremendous excitement <cit.>. At small twist angles, the moiré pattern creates a superlattice with a large unit cell comprising more than 10000 atoms, and dramatically modify the low-energy electronic structure. In particular, near certain “magic” twist angles,  four lowest-energy minibands with a total band width on the order of 10meV are separated from excited bands and accommodate a range of carrier densities from charge -4e to 4e per supercell. Due to the strong suppression of kinetic energy in these narrow bands, Coulomb interaction may drive correlated electron phenomena <cit.>. Remarkably, the recent experiments <cit.> on such twisted bilayer graphene (TBG) discovered metal-insulator transition and superconductivity at low temperature by tuning the carrier density, applying the magnetic field or slightly varying the twist angle. These fascinating phenomena show a number of similarities with that of cuprates. Notably, a correlated insulating state occurs below 4K  at the filling of charge ± 2e per supercell. Under electrostatic doping, two superconducting domes appear on both sides of the insulating state, with a maximum transition temperature T_c=1.7K and a record-low carrier density of a few 10^11cm^-2. The mechanism of metal-insulator transition and superconductivity, the nature of the correlated insulating state and the superconducting state are all open questions.

The Hubbard model is the standard model to study metal-insulator transition driven by the competition of kinetic energy and Coulomb interaction <cit.>. It is also believed to capture key features of the cuprate superconductors <cit.>. To study metal-insulator transition in TBG, it is highly desirable to identify the real-space Wannier orbitals for the low-energy miniband and find the corresponding tight-binding and Hubbard model.

This is however a nontrivial task that has not been accomplished so far. Thanks to   extensive studies using various methods <cit.>, the band structure of TBG at small twist angle is known to be rather complex and depend sensitively on microscopic details such as lattice relaxation. Near the so-called magic twist angle, various methods find four nearly-flat minibands at low energy, but differ significantly on important features such as their bandwidth and the gap to excited bands. Therefore, as a first step,
it is important to extract robust and universal features of these narrow minibands from both theoretical calculation and experimental findings on bilayer graphene at small twist angles.



In this work, we demonstrate that the electronic structure of narrow minibands and the effect of Coulomb interaction in twisted bilayer graphene are essentially captured by a two-orbital Hubbard model, constructed from Wannier orbitals that extend over the size of supercells. We deduce the centers and symmetry of Wannier orbitals by a straightforward symmetry analysis, without explicitly computing their wavefunctions.  Importantly, the centers of these Wannier orbital form an emergent honeycomb lattice. The two types of sublattice sites of this honeycomb lattice correspond to AB and BA regions of twisted bilayer graphene respectively, while the hexagon centers correspond to AA regions, as depicted in Fig. <ref>. At every site of the honeycomb lattice, there are two degenerate Wannier orbitals with p_x- and p_y-like symmetries, forming a doublet under on-site three-fold rotation.  We then construct an effective tight-binding model on this honeycomb lattice, which reproduces key features of the miniband structure of TBG <cit.>. By including Coulomb repulsion, our model provides a useful theoretical basis for studying the metal-insulator transition in TBG as a function of twist angle and carrier density, Dirac fermion reconstruction at charge neutrality, as well as other strongly correlated phenomena such as unconventional superconductivity.




The procedure of our analysis is outlined  in Fig. <ref>. First, from general considerations on the lowest minibands of TBG and tight-binding calculation by Nam and Koshino <cit.>, we infer the band symmetry eigenvalues at all high symmetry points Γ ,K,M of mini Brillouin zone (MBZ) of TBG.
Then we examine all possible positions of Wannier centers (which form a lattice in real space) and symmetries of Wannier orbitals (s-, p-wave etc) to search for solutions consistent with the band symmetries in k-space. Luckily, we find that the band symmetries at all high symmetry points can only be reproduced when Wannier orbitals with (p_x,p_y) on-site symmetry are located on a honeycomb lattice. Based on this important result, we construct the simplest tight-binding model that reproduce the key features of band structure, in particular, the warped Fermi surfaces near the miniband edges at Γ point.

This work is organized as follows. In Section I, we study the band structure and energy eigenstates of TBG through the group-theoretical approach. From the obtained energy eigenstates, in Section II we deduce the positions of Wannier orbital centers and the symmetries of Wannier orbitals. Based on the Wannier orbitals, we then construct tight-binding model in Section III. Then, we address the Hubbard model and metal-insulator transition in Section IV, Dirac fermion reconstruction and Landau level degeneracy breaking in Section V, and make connections between theoretical and experimental results. At last in Section VI, we discuss some open questions such as the nature of correlated insulating and unconventional superconducting phases of TBG.






§ BAND STRUCTURE

The first question we ask is: are the lowest bands on the electron and hole sides separated from excited bands? Only when low- and high-energy bands are separated by a band gap, a theoretical description using a small number of low-energy degrees of freedom is possible.

To answer this question, we look into theoretical calculation and experimental evidence on the band structure of TBG. Different numerical methods all show that band edges of the lowest-energy electron and hole minibands near charge neutrality are located at Γ point of the MBZ. However, conflicting results are found on the gap at Γ point between these lowest bands and higher-energy bands when the twist angle is small. Some numerical calculations report relatively large gaps of 10-20meV on both electron and hole sides <cit.>, while others show that the gap only exists on the electron side <cit.> or even no gap exists <cit.>.

In the experiment performed on TBG near the twist angle θ =  1.08^∘, the conductance is found to be zero over a wide range of electron or hole densities near n=± n_s = ±  2.7 × 10^12cm^-2—the density at which the four lowest bands are completely filled or empty. The measured thermal activation gaps are about 40meV, comparable to the single-particle band gaps found by Nam and Koshino in tight-binding calculation with relaxed lattice structure <cit.>. Hence we conclude that at small twist angles, the lowest bands of TBG are well separated from excited bands above and below.

Our next goal is to understand the band structure and Bloch wavefunction of the lowest minibands analytically. Previous analytical calculations have mostly focused on the minibands near the corners of MBZ <cit.>, where the original Dirac points of graphene remain present, but exhibit a much reduced Fermi velocity. However, a full analysis of energy bands at all high symmetry points remains lacking. 

To do this, we need to work out energy eigenstates of TBG, which are superpositions of states on two layers hybridized through interlayer tunneling.  In real space, the interlayer tunneling takes the general form

    ℋ_T=∑_m, n=A,B∑_ x^m, y^nξ^†( x^m) T( x^m, y^n) η( y^n)+h.c.

where ^m, ^n denote the coordinates of carbon atom sites on layer 1 and 2 respectively and m,n denotes the A/B sublattice. ξ and η are the electron annihilation operators on the two layers respectively.
T describes the tunneling amplitudes between two sites  and  on different layers.

We define electron operator at a given momentum k as follows:
ξ_ k=(ξ_ k^A,ξ_ k^B)^T,η_ k=(η_ k^A,η_ k^B)^T where

    ξ^n_ k =   ∑_ x^n e^i  k· x^nξ(^n),     η^n_ k =   ∑_ y^n e^i  k· y^nη(^n)

with n=A, B. By employing the two-center approximation T( x, y)=T( x- y) <cit.>, we can write Eq (<ref>) in momentum space as

    ℋ_T=∑_ q,𝒢_1,𝒢_2ξ_ q+𝒢_1^†T_ q(𝒢_1,𝒢_2)η_ q+𝒢_2+h.c.

where 𝒢_1,2 is the reciprocal vector of layer 1 or 2 respectively, and T_ q(𝒢_1 ,𝒢_2) is the 2 by 2 interlayer scattering matrix with momentum transfer 𝒢_1 -𝒢_2.


To calculate the energy eigenstates explicitly requires the knowledge of the tunneling operator T. However, when the twist angle is commensurate (as we assume throughout this work), T satisfies certain symmetry conditions, hence at high symmetry points of MBZ the low-energy eigenstates can be classified by the irreducible representations of the corresponding symmetry group. Importantly, essential information on these eigenstates and hence Wannier orbitals can be deduced directly from the symmetry representations. In the following three subsections, we will work out the symmetries of the four lowest minibands at all three high symmetry points: Γ,K,M of MBZ respectively.

Before proceeding, we first introduce the coordinate system. As shown in Fig. <ref>a, the graphene superlattice has a point group D_3 with in-plane three-fold rotation and out-of-plane two-fold rotation. We then choose the axis of two-fold rotation as y-axis, and the axis of three-fold rotation as z-axis. Thus two-fold and three-fold rotations can be denoted as C_2y and C_3z. The center of C_3z, denoted as O in Fig. <ref>a, is located at a vertically aligned A site, which we define as the origin of coordinate == 0. In this coordinate system, the primitive vectors of the graphene superlattice are A_1,2=A(√(3)/2,±1/2) where A=a/[2sin(θ/2)] is the superlattice constant, a is the common lattice constant of both layers and θ is the twist angle. Before band folding, the +K point of the top (bottom) graphene layer in its original Brillouin zone is K_1,2=4π/3a[cos(θ/2),±sin(θ/2)]. After band folding, + K_1 and -  K_2 fold onto one corner of MBZ of TBG, while - K_1 and +  K_2 fold onto the other inequivalent corner, as shown in Fig.<ref>b.


In our coordinate system, the scattering matrix in (<ref>) can now be explicitly written as 

    T_ q^mn(𝒢_1 ,𝒢_2)=t^mn( q)exp[-i(𝒢_1·τ_1^m-𝒢_2·τ_2^n)]

where t^mn( q) is the Fourier transform of T( x^m - y^n), and
the vectors τ^A,B_j specify the position of A/B sublattice sites within the unit cell of layer j, given by

    τ^A_1,2   =    0,  
    τ^B_1,2   =   a/√(3)[∓sin(θ/2),cos(θ/2)].

We emphasize that  the phase factor in the interlayer tunneling matrix (<ref>) comes from the intra-unit-cell part of the Bloch wavefunction of each graphene layer, and must be included to obtain the correct band structure of TBG.

Since the interlayer tunneling term ℋ_ T is much smaller than the full bandwidth of graphene layers, the main effect of ℋ_ T is to modify the low-energy part of the band structure near the original Dirac points ± K_1,2.
Therefore, to calculate the low-energy band structure of TBG, it suffices to restrict q+𝒢_1,2 in the interlayer tunneling term (<ref>) to be close to ± K_1,2 points. Then, the full Hilbert space comprising a large number of momentum points is truncated to a reduced Hilbert space comprising only a few momentum points. This  approximation is known as the continuum theory of TBG <cit.>.



 §.§ Γ Point

We first study the lowest minibands of TBG at Γ point of the MBZ. There are in total three pairs of opposite momenta, denoted as ±Λ_a=1,2,3,  which are located closest to the original Dirac points K_1,2 of individual graphene layers and fold onto Γ of MBZ. These three Λ points are all integer multiples of the two superlattice reciprocal vectors G_1,2=G_θ(±1/2,-√(3)/2),G_θ=8π/√(3)asinθ/2, and related by three-fold rotation symmetry, as shown in Fig. <ref>b.

Interlayer tunneling leads to scatterings among ±Λ points on layer 1 and those on layer 2. Among all possible scattering processes, intervalley scattering between +Λ points to -Λ points requires large momentum transfer, hence is negligibly small when the twist angle is small <cit.>. Thus we only consider intravalley scattering within the three +Λ points that are all close to the original Dirac points + K_1,2, and drop the valley index. Since we focus on + K_1,2 valleys, we only need to consider two out of the four lowest energy eigenstates at Γ point, which are denoted by Ψ_e^Γ on electron side and Ψ_e^Γ on hole side. The other two lowest energy eigenstates coming from - K_1,2 valleys will be time-reversal partners of Ψ_e,h^Γ.

According to our discussion above, to the lowest order approximation, energy eigenstates at Γ are in general superpositons of ξ^m_ k,η^n_ k where k=Λ_a   (a=1,2,3). While the energy spectrum depends on the form of tunneling operator T(^m, ^n), we now deduce purely from general symmetry considerations the essential form of these wavefunctions.

Any energy eigenstate at Γ must belong to one of the three irreducible representations of the point group D_3: the one dimensional A_1 representation, the one-dimensional A_2 representation, and the two-dimensional E representation.
The two-fold rotation C_2y around y-axis maps + K_1,2 to - K_2,1 respectively, i.e., simultaneously interchanges the two valleys and the two layers. Therefore, its action cannot be represented within the subspace of states around one valley. Thus it suffices to consider the subgroup C_3 that acts within +Λ states belonging to the + K valley.

In the subgroup C_3, every representation is labeled by the eigenvalue of angular momentum L_z taking three possible values: 0, ± 1. According to group theory, an eigenstate Ψ with angular momentum L_z formed by Bloch states from three Λ points can be generally written as

    Ψ^Γ   =   ∑_m∑_ x^mα_m∑_a=1^3e^iΛ_a· x^m+2i/3aL_zπξ( x^m)
       +   ∑_n∑_^nβ_n∑_a=1^3e^iΛ_a· y^n+2i/3aL_zπη( y^n)

where α_m,β_n are complex coefficients.
Without loss of generality in the following we shall consider the energy eigenstate of the conduction mini-band, denoted by Ψ_e^Γ. The analysis of the valence mini-band Ψ_h^Γ is similar.

We first focus on the case of L_z = +1. According to the general expression above, Ψ^Γ can be rewritten in the following suggestive way

    Ψ^Γ   =    e^-i 2π/3[ α_A  U^A( R_c)+α_B  U^B (- R_c) ] 
       +    e^i 2π/3[ β_A L^A (- R_c) +β_B  L^B ( R_c) ].

Here R_c=Ax̂/√(3) is the coordinate of an BA spot closest to the AA spot at the origin, as shown in Fig. <ref>. U and L are electron wavefunctions on the upper and lower layers respectively, defined by

    U^m( R)   =   ∑_ x^me^i K_1· x^m f( x^m- R)ξ( x^m),  
    
    L^n( R)   =   ∑_ y^ne^i K_2· y^nf^*( y^n- R)η( y^n).

Both U,L are the product of intra-unit-cell wavefunction e^i K· r which is fast oscillating and the envelope function f( r- R) given by

    f( r)=e^i( K_1- K_2)· r×{1+e^i G_1· r+e^i G_2· r},

which is slowly varying. f( r) is invariant under three-fold rotation around the origin.

It is straightforward to show that the maxima of the envelop function |f( r)| are located at AA spots n_1 A_1+n_2 A_2   (n_1,2∈ℤ), which form a triangular lattice with primitive vectors A_1,2. Then, the maxima of |f( r ± R_c)| are located at the AB/BA spots ∓ R_c+n_1 A_1+n_2 A_2, which also form a triangular lattice but shifted off the origin by ∓ R_c. Therefore, it follows from Eq. (<ref>) that  the component of Bloch wavefunction Ψ^Γ on the A sublattice has its maxima at BA spots on layer 1 and at AB spots on layer 2, while the component on the B sublattice has its maxima on AB spots on layer 1 and BA spots on layer 2.

This feature is robust and can be understood by symmetry considerations. For a state carrying finite angular momentum such as L_z=+1 considered here, its wavefunction necessarily vanishes at rotation centers n_1 A_1+n_2 A_2, where A sublattice sites on two layers are registered. Therefore, the maxima of A-sublattice component of Ψ^Γ must be away from these AA spots. Furthermore, if there is only a single maximum within a supercell (as is the case here), this maximum can only be located at either AB or BA spots, because these positions are invariant under three-fold rotation with respect to AA spots up to superlattice translations. Finally, we note that under the combination of two-fold rotation C_2y and time-reversal symmetry 𝕋, the two layers are interchanged, while the sublattice and angular momentum L_z quantum numbers are unchanged.  This implies that the maxima of A sublattice wavefunction is symmetric under C_2y, hence must be located at BA spots on layer 1 and AB spots on layer 2, forming the C_2y image of each other.
On the other hand, B sublattice sites do not coincide with rotation centers, and the intra-unit-cell phase factors e^i  K_1 · x^B and e^i  K_2 · y^B already carry the angular momentum -1. To make the total angular momentum L_z=+1, the envelope function of B sublattice wavefunction must carry -1 angular momentum and is therefore peaked at AB spots on layer 1 and BA spots on layer 2. Here we have used the fact that under C_3z, angular momentum is defined modulo 3.

We thus conclude that, to respect the point group D_3 of twisted bilayer graphene, the energy eigenstate at Γ point with angular momentum L_z=+1 has the peculiar property that A or B sublattice wavefunction has its maxima on a honeycomb lattice whose two sublattices correspond to AB and BA spots of the two layers. Similar conclusion is found for eigenstates with L_z=-1, which we leave to the Supplementary Material <cit.>.

In the band structure calculations by Nam and Koshino <cit.>, at Γ point the eigenstates form two doublets—one on electron side and one on hole side, and each doublet splits along Γ M line in MBZ. It is our understanding that such band structure is consistent with the scenario that the two members of the doublet come from ± K valleys of graphene and have L_z=± 1 respectively. In this case, the two-fold degeneracy is protected by D_3 point group and time-reversal symmetry and remains intact even when inter-valley scattering is taken into account.
In contrast, in the scenario of L_z=0, there will be four non-degenerate states at Γ, which is not found in band structure calculations. For completeness, we leave the discussion of L_z=0 case to the Supplementary Material <cit.>.

So far our analysis has been completely based on symmetry considerations. To work out the eigenstates Ψ_e,h^Γ at Γ point explicitly requires a microscopic calculation, which we sketch below. Before turning on interlayer tunneling, the Hamiltonian for graphene states at Λ points is given by

    H_0 =
     v_F∑_a=1^3 [ξ^†_a(_a·σ^1)ξ_a +  η^†_a(_a·σ^2)η_a].

Here ξ_a and η_a are respectively the electron operator of layer 1 and 2 at momentum Λ_a. In H_0, _a and _a are the small momentum differences between Λ_a and the Dirac points on layer 1 and 2 respectively, and v_F is the Fermi velocity.
Due to their difference in orientation, the Pauli operators in the Dirac Hamiltonian of individual graphene layers are “rotated” oppositely:

    σ^1,2 = (σ_xcosθ/2∓σ_ysinθ/2, 
     ±σ_xsinθ/2+σ_ycosθ/2 ),

where σ_x,y are Pauli matrices.

According to Eqs.(<ref>) and (<ref>), the interlayer tunneling between Λ points can be written as

    H_T=∑_a,b=1^3ξ_a^†T_abη_b+h.c.

where T_ab is the scattering matrix obtained by integrating out intermediate high-energy states in Eq.(<ref>) other than the low-energy ±Λ_a  states kept in Eq.(<ref>).

The three-fold rotation of Λ points is realized by the cyclic permutation Λ_a→Λ_a+1 (with Λ_4 ≡Λ_1).
Due to the invariance of H_T under three-fold rotation and with a proper choice of relative phases of ξ_a and η_a, we have T_a+1,b+1=T_ab. Hence essentially there are only three independent scattering matrices among all of them, and we denote them as T_0=T_aa,T_±=T_a± 1,a. Physically, interlayer tunneling associated with T_0 preserves momentum while T_± transfers momentum ±( G_1- G_2) from layer 2 to layer 1. 

With (<ref>) and (<ref>), we can work out the lowest eigenstates Ψ_e^Γ and Ψ_h^Γ at Γ point. An illustrative example of Ψ_e^Γ and Ψ_h^Γ with L_z =+1 is shown in Fig. <ref> where A sublattice component dominates and its maxima are indeed located at the emergent honeycomb lattice formed by AB and BA spots. Details of interlayer coupling Hamiltonian (<ref>) in this example can be found in the Supplementary Material <cit.>.





 §.§ K points

Similar analysis as shown above applies to energy eigenstates at other high symmetry points in the MBZ. Consider the corners of MBZ, denoted by ± K where K≡ K_1- K_2. Since time-reversal symmetry relates bands at ± K by complex conjugation, it suffices to study + K only. Since + K_1 and - K_2 in the large Brillouin zone fold onto + K in the MBZ, low-energy states of TBG at + K come predominantly from the states of layer 1 at Dirac point + K_1 (denoted by ξ^A,B_+) and the states of layer 2 at Dirac point - K_2 (denoted by η^A,B_-). When the intervalley scattering is negligibly small, to leading order approximation we obtain four degenerate zero-energy states at + K, as found in previous studies.

At K point, the point group D_3 is preserved and hence the eigenstates can also be labeled by representations A_1,A_2,E of D_3. Under three-fold rotation, we find that states at A sites of each layer have zero angular momentum, and ξ^A_+ +η^A_- will furnish the trivial representation A_1 while ξ^A_+ -η^A_- will furnish the 1D representation A_2. The remaining two states (ξ^B_+ ,η^B_-) at B sites have finite, opposite angular momenta. Furthermore, under two-fold rotation C_2y, the type of sublattice is unchanged but + K_1 and - K_2 are interchanged. As a result, (ξ^B_+ ,η^B_-) furnish the 2D representation E of D_3.

Previous studies found that the correction to the above energy eigenstates from excited states of individual graphene layers will not lift the degeneracy at Dirac point but only reduces the Fermi velocity. We shall come back to this point later.



 §.§ M points

At M points, only in-plane two-fold rotation C_2y is present. Thus one can classify eigenstates at M points according to the eigenvalues ± 1 of C_2y. Energy eigenstates at M points are from the state with momentum ±1/2( K_1 + K_2) on both layers. Since C_2y interchanges the two layers and the two momenta, we obtain two energy eigenstates with C_2y=+1 and two with C_2y=-1.

The representations of the eigenstates at high symmetry points Γ,K,M are summarized in Table. <ref>.



§ WANNIER ORBITALS AND CENTERS

In this section we deduce the center and symmetry of Wannier orbitals from symmetries of band eigenstates at high symmetry points through consistency check. While the detailed form of Wannier orbitals are to be obtained by Fourier transform of band eigenstates over the MBZ, the center and symmetry of Wannier orbitals are robust features that only depend on symmetries of energy bands at all high symmetry points.

Among all possible configurations of Wannier orbitals, we first consider the case of degenerate (p_x ,p_y)-like Wannier orbitals at the honeycomb lattice formed by AB and BA spots as shown in Fig. <ref>. This choice is motivated by the property of energy eigenstates Ψ^Γ shown earlier. We show below that  the case of (p_x ,p_y) Wannier orbitals on the honeycomb lattice yields the correct band symmetries at all high symmetry points of MBZ shown in Table. <ref>.

To see this explicitly, it is useful to construct the electron Bloch state with momentum k∈MBZ in p_x± i p_y orbital basis  as follows

    c_ k,τ=∑_je^i k· R_jc_jτ

where τ=± denotes the Wannier orbital with L_z= ± 1 under three-fold rotation around its own center, R_j is the coordinate of honeycomb lattice site j, c_jτ annihilates an electron at site j with orbital τ, and the sum is over all sites of the honeycomb lattice.




At Γ point k= 0, the plane-wave phase factor e^i k· R_j=1 does not contribute angular momentum, and the angular momentum of Bloch  wavefunctions c_ 0,τ comes solely from the p_x ± i p_y on-site symmetry of Wannier orbital c_jτ, leading to the desired angular momentum L_z=τ =± 1. The two-fold rotation C_2y interchanges layer indices 1 and 2, and maps p_x+ip_y to -p_x+ip_y and vice versa. As a result, the p_x+ip_y Wannier orbital at AB spot and the p_x-ip_y orbital at BA spot will furnish 2D representation E of the D_3 point group. The same holds for the p_x+ip_y orbital at BA spot and p_x-ip_y orbital at AB spot. This leads to two doublets as energy eigenstates at Γ point, consistent with our results in Table. <ref>.


At K point k= K_1- K_2, 
the plane-wave phase factor e^i k· R_j on the honeycomb lattice generates nonzero angular momentum, which adds to the orbital angular momentum of p_x ± i p_y Wannier orbitals. To see this explicitly, we decompose the Bloch state (<ref>) as the superposition of Bloch waves in different sublattices

    c_ k,τ=c_ k,τ,1+c_ k,τ,2,   c_ k,τ ,n≡∑_j∈ L_ne^i k· R_jc_jτ,

where L_1,2={∓ R_c+n_1 A_1+n_2 A_2|n_1,2∈ℤ} denote the set of AB/BA spots respectively.

We then find under the action of C_3z around the origin,

    c_ K,τ,1   →   ∑_j∈ L_1e^i K· C_3z R_je^-2iτπ/3c_jτ
       =   ∑_j∈ L_1e^i K· R_j+2iπ/3e^-2iτπ/3c_jτ
       =   e^-2i(τ -1)π/3c_ K,τ,1

and hence the Bloch state c_ K,τ,1 has angular momentum τ -1 with an additional -1 angular momentum from Bloch wave phase factor e^i k· R_j. Likewise, the Bloch state c_ K,τ,2 has angular momentum τ +1 under three-fold rotation C_3z. With τ=± 1 we find the four Bloch states c_ K,τ,n fall into two categories: one pair of doublet (c_ K,-,1,c_ K,+,2) and two singlets c_ K,+,1,c_ K,-,2. Notice that the two-fold rotation C_2y maps c_ K,τ,1 to -c_ K,-τ,2 and vice versa, we conclude that the doublet (c_ K,-,1,c_ K,+,2) furnishes the 2D representation E, while the singlets c_ K,+,1∓ c_ K,-,2 furnish the 1D representation A_1,2 respectively. The case with - K point can be obtained through time-reversal operation. 

At M points k=1/2( K_1+ K_2), the little group is C_2. As discussed in section I, the Abelian group C_2 can only host 1D representations denoted as A and B. Since M=1/2( K_1+ K_2) is invariant under C_2y up to a reciprocal vector in MBZ, we find c_ M,+,1- c_ M,-,2 and c_ M,-,1- c_ M,+,2  belong to A representation while c_ M,+,1+ c_ M,-,2 and c_ M,-,1+c_ M,+,2  belong to B representation.

We further show in the Supplementary Material <cit.> that from an exhaustive case-by-case study, any other type of lattice or orbital symmetry is incompatible with the band symmetry at all high symmetry points. Thus, we conclude that the Wannier orbitals of TBG have p_x,y symmetries and form a honeycomb lattice whose sites correspond to AB and BA spots.





The exact form of Wannier orbitals of TBG are to be constructed from energy eigenstates at all momenta in the MBZ, which we leave for a future work. We expect that these Wannier orbitals will extend over multiple supercells of TBG.





§ TIGHT-BINDING MODEL


In this section we construct a tight-binding model on the honeycomb lattice with (p_x,p_y)-like orbitals as an effective model for the four lowest bands of TBG. Our construction is guided by the D_3 point group symmetry of TBG and the  valley U(1) symmetry that emerges in the absence of inter-valley scattering.

As is well-known, p-orbitals are allowed to have two types of hopping, σ (head-to-tail) and π hopping (shoulder-to-shoulder). When the σ and π hopping amplitudes are the same, we arrive at the simplest tight-binding model of p-orbitals on the honeycomb lattice

    H_0=-∑_iμ^†_i·_i+∑_⟨ ij⟩    t_1 [^†_i·_j + h.c.]
    
    +∑_⟨ ij⟩' t_2[^†_i·_j+ h.c.]

where _i= (c_i,x, c_i,y)^T with c_i,x(y) annihilating an electron with p_x(y)-orbital at site i. t_1 is the real hopping amplitude between nearest neighbor sites of different sublattices. t_2 denotes real hopping amplitudes between fifth nearest neighbors with bond length √(3)A, or equivalently, second nearest neighbors within the same sublattice. With only the t_1 term, the band structure is particle-hole symmetric, while the t_2 term breaks this symmetry.
μ is the on-site chemical potential. Further neighbor hoppings can be included as well.

The Hamiltonian (<ref>) has emergent SU(4) symmetry and hence every band is four-fold degenerate including orbital and spin. This result however does not match with numerical band structure by Nam and Koshino <cit.> where the orbital degeneracy is found broken along Γ M lines. To describe such effect, we include the following Hamiltonian

    H_1   =   ∑_⟨ ij⟩'t_2'[ (^†_i×_j)_z+ h.c.]
       =   -i∑_⟨ ij⟩'t_2' (c^†_i+c_j+-c^†_i-c_j-)+ h.c.

In terms of the chiral basis c_±=(c_x± ic_y)/√(2) associated with p_x± ip_y orbitals, the hopping terms in (<ref>) become imaninery values. As illustrated in Fig. <ref>b, orbitals of different chirality experience finite and opposite magnetic fluxes, and the model as a whole preserves time-reversal symmetry. Microscopically the chiral Wannier orbitals c_± originate from states near different valleys ± K of graphene. When intervalley coupling is negligible, the orbital U(1) symmetry c_τ→ e^i τϕ c_τ is respected.
H_1 breaks SU(4) symmetry down to U(1)×SU(2), where U(1) refers to orbital chirality conservation and SU(2) refers to spin rotation symmetry.



In principle, we can further include the Hamiltonian breaking orbital U(1) symmetry as follows

    H_2 =∑_⟨ ij⟩ t'_1 [ ^†_i· e^∥_ij  e^∥_ij·_j -^†_i· e^⊥_ij  e^⊥_ij·_j+ h.c.]

where e^∥,⊥_ij denote in-plane unit vectors in the direction parallel and perpendicular to the nearest neighbor bond ⟨ ij⟩ respectively. This symmetry-breaking term arises when intervalley scattering is included, and leads to Dirac mass generation which we will discuss later in Section V.

Symmetry-breaking terms (<ref>) and (<ref>) have important effects on the band structure. Denoting τ as Pauli matrices acting in chiral orbital index and σ as Pauli matrices acting in the sublattice index of the honeycomb lattice, then in the chiral basis the tight-binding Hamiltonian H_tb leads to the k· p Hamiltonian near Γ point in MBZ

    H( k)   =   H_1( k)+H_2( k),
    
    H_1( k)   =   t_1(3-A^2/4| k|^2)σ_x
    +v k·τσ_y-μ,
    
    H_2( k)   =   t_2(6-9A^2/2| k|^2)+λ(k_+^3 + k_-^3) τ_z,

where v=√(3)/2At_1',  λ =-3√(3)/8A^3 t'_2,   k_±=k_x± ik_y.

Interestingly, the t_2' term (k_+^3+k_-^3)τ_z  is the orbital analog of hexagonal warping in spin-helical surface states of topological insulators <cit.>.
In the full MBZ, the hexagonal warping term vanishes along Γ K and MK lines and becomes finite along Γ M lines. Hence orbital degeneracy is preserved by hexagonal warping term along Γ K and MK lines, consistent with numerical band structures <cit.>. The t_1' term further beaks the orbital degeneracy at any point of the MBZ except Γ, K, which corresponds to the reduced symmetry group D_3×SU(2) where SU(2) is spin rotation group and D_3 is point group acting jointly on lattice sites and (p_x,p_y) orbitals. These effects can be seen from band structures in Fig. <ref>c and d where the band structure in Fig. <ref>c with t_1'=0 agrees with calculations by Nam and Koshino <cit.>, and the band structure in Fig. <ref>d with finite t_1' shows no band degeneracy except Γ,K points where D_3 is respected and 2D representation is realized.





§ HUBBARD MODEL AND METAL-INSULATOR TRANSITION

From numerical calculations on band structure of TBG, we expect the lowest band width and hence hopping parameters t_1,t_2,t'_1,t'_2 are of order of meV. With such small kinetic energy, the effect of interaction have to be considered. In terms of Wannier orbitals, we write down the general form of symmetry-allowed two-orbital on-site Hubbard interaction Hamiltonian

    H_UJ=U∑_i,τn_iτ↑n_iτ↓+U'∑_in_ixn_iy
       +   J∑_i,s,s'c^†_ixsc^†_iys'c_ixs'c_iys
    +J'∑_i,τ≠τ'c^†_iτ↑c^†_iτ↓c_iτ'↓c_iτ'↑

where U,U',J,J' denote the intraorbital Coulomb, interorbital Coulomb, exchange, and pair-hopping interactions, respectively. τ=x,y is the orbital index and s=↑,↓ is the spin index. Here n_iτ s=c_iτ s^†c_iτ s is the number of electrons at site i with orbital τ and spin s, and n_iτ=n_iτ↑+n_iτ↓. 
Since we choose real p_x,y orbitals, we further find J=J' denotes the Hund's coupling.


In general, the density-density interactions should be comparable with each other but much larger than Hund's coupling, U∼ U'≫ J,J'. At the simplest level of approximation, we set U=U' and write down the two-orbital Hubbard model as H=H_0+H_1+H_2+H_U where the on-site interaction is

    H_U=1/2 U∑_i(n_i-1)^2   (n_i=n_ix+n_iy),

and H_0,1,2 are given in (<ref>), (<ref>) and (<ref>).


As mentioned earlier in Section II, the Wannier orbitals are expected to extend over the supercell. Hence electron-electron interaction between nearby sites

    H_V=∑_ijV_ijn_in_j

is also present and may be important. Finally we note that the strength of Coulomb repulsion is reduced by screening from excited bands of TBG which span a wide range of energies from 10 meV to 10 eV. Furthermore, in typical graphene devices, the long-range Coulomb interaction is screened by nearby metallic gates at a distance comparable to the size of the supercell. We thus expect a model with a few short-range interactions is sufficient.


We now address the metal-insulator transition in the two-orbital honeycomb Hubbard model and connect our findings to experimental results. With two sites per unit cell and two orbitals at every site, our honeycomb lattice model can accommodate up to 8 electrons per unit cell, which corresponds to the complete filling of the miniband in TBG. The charge neutrality point of TBG corresponds to on average n=4 electrons per unit cell or equivalently q=2 electrons per site of honeycomb lattice in our model. The Mott insulator state found in experiments occurs at n=2 electrons/holes per unit cell, or equivalently q=1 electron/hole per site.

Generally speaking, in Hubbard models with an integer average number q of electrons per site, a transition from metal to Mott insulator is induced by increasing the ratio of Coulomb repulsion and bandwidth U/t. When the twist angle approaches the magic value, the miniband bandwidth decreases very rapidly <cit.>. Therefore varying the twist angle slightly induces a bandwidth-controlled metal-insulator transition.
Importantly, in two-orbital Hubbard models, the critical value of U depends on q. In the SU(4) symmetric limit, the critical coupling U_c for q=2 is known to be larger than that for q=1, because orbital fluctuations are largest in the former case <cit.>. This is consistent with the experimental finding of the insulating state at q=1, but not at q=2 (charge neutrality).



Transport measurement reveals that the insulating state in TBG has a thermal activation gap of 0.3meV, and is destroyed by the application of an in-plane or perpendicular magnetic field of 8T corresponding to a Zeeman energy of 0.5meV <cit.>. This implies that the charge gap of the correlated insulator is much smaller than the bandwidth. Therefore, the insulating state is a weak Mott insulator close to a metal-insulator transition, ruling out the possibility of U≫ t. Moreover, before it becomes highly insulating at low temperature,  above 4K resistivity increases linearly with temperature. This behavior is characteristic of high-temperature bad-metal regime of a Mott insulator with or without doping, seen in many oxides and in numerical study of Hubbard model <cit.>.











§ DIRAC FERMION RECONSTRUCTION AT CHARGE NEUTRALITY


We now highlight a interesting finding in the experimental paper <cit.>, whose significance may have not been fully recognized.
Over a wide range of densities around charge neutrality, Landau level degeneracy is  found to be 4-fold instead of 8-fold, as expected from the 8 degenerate massless Dirac fermions coming from graphene valley, electron spin, and ± K points of MBZ.


Our work offers a possible solution to this problem.
As we showed previously in Section I and III,
after band folding the Dirac doublet at + K_1 in layer 1 and the one at - K_2 in layer 2 will fall on top of each other at + K of MBZ. From symmetry considerations, the resulting four states at + K consist of one doublet (E)  and two singlets (A_1 and A_2) representation of the D_3 point group. The same result applies to the - K point. Therefore half of the Dirac fermions in the doublet representation are protected by symmetry, while the other half in the singlet representation are unstable and can become massive if perturbations are considered.

First, the interlayer tunneling Hamiltonian (<ref>) contains terms that scatter states at different valleys on different layers, which generates single-particle mass for the half unstable Dirac fermions. However, as we mentioned earlier,  the intervalley scattering requires large momentum transfer and is very weak when the twist angle is small.

On the other hand, when the twist angle approaches the magic value, the miniband bandwidth decreases rapidly <cit.>, lattice relaxation and Coulomb interaction become important and can strongly renormalize the single-particle band dispersion. In particular, we envision that these correlation effects can significantly enhance the single-particle intervalley scattering  and hence the corresponding Dirac mass term, without breaking any lattice symmetry. If this scenario is correct, we expect the correlation-enhanced Dirac mass in TBG near charge neutrality will be controlled by the twist angle, becoming very small away from the magic angle. This prediction can be tested by systematically studying Landau level degeneracy in TBG with different twist angles.

In terms of our tight-binding model, the mass generation term corresponds to t_1' term in (<ref>) which break U(1) symmetry in orbital space. As we discussed in Section. III and shown in Fig. <ref>a and b, for t_1' ≠ 0, the unstable Dirac fermions are gapped out at low-energy and only half of Dirac fermions remain robust. As a result, the  Landau level degeneracy becomes 4-fold, coming from ± K points of MBZ and electron spin.



§ OUTLOOK


The two-orbital honeycomb Hubbard model proposed in this work provides a theoretical  framework for studying correlated electron phenomena in graphene superlattices. Many important questions remain to be addressed, among which we highlight a few. First, the nature of the correlated insulator ground state at q=1. It may exhibit long-range order, such as spin, orbital or valence bond solid order. Alternatively, because of the proximity to metal-insulator transition and/or the presence of orbital fluctuations <cit.>, the correlated insulator may be a quantum spin liquid. In this regard, it is worth noting that in the special limit of our Hubbard model with SU(4) symmetry, the effective Hamiltonian for the correlated insulator at strong coupling is the SU(4) generalization of Heisenberg model (Kugel-Khomskii model) on the honeycomb lattice. Analytical and numerical studies of this model indicate the lack of any long-range order <cit.>, and suggests a possible quantum spin liquid with gapless neutral fermionic excitations <cit.>. Thermal conductivity measurements can tell the existence or absence of gapless neutral excitations in the correlated insulator state of TBG.

Our model also provides a starting point for studying superconductivity in TBG. The strong electron repulsion disfavors on-site pairing, and opens the possibility of unconventional pairing symmetry. We leave the study of superconductivity in the two-orbital honeycomb Hubbard model to a future work.




§ ACKNOWLEDGMENT


We thank Pablo Jarillo-Herrero, Yuan Cao and Mikito Koshino for invaluable discussions. This work is supported by DOE Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award DE-SC0010526. LF is partly supported by the David and Lucile Packard Foundation.

Note added: Recently the prediction of this work based on symmetry analysis that  Wannier orbitals of the lowest mini-bands of twisted bilayer graphene have  (p_x, p_y) symmetry and form a honeycomb lattice has been verified by explicit numerical calculations <cit.>. These works also found tight-binding models of the same form as ours. In contrast, a different honeycomb model was proposed in Ref. <cit.>, where orbitals have two different on-site energies, instead of being degenerate as dictated by the D_3 lattice symmetry of twisted bilayer graphene.

99
Cao1 Yuan Cao, Valla Fatemi, Ahmet Demir, Shiang Fang, Spencer L. Tomarken, Jason Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori and P. Jarillo-Herrero, Nature, doi:10.1038/nature26154 (2018).

Cao2 Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras and Pablo Jarillo-Herrero, Nature, doi:10.1038/nature26160 (2018).

Mele E. J. Mele, Nature News and Views, doi:10.1038/d41586-018-02660-4.

Xu Cenke Xu and Leon Balents, arXiv:1803.08057.

Volovik G.E. Volovik, arXiv: 1803.08799 .

MacDonald Rafi Bistritzer and Allan H. MacDonald, PNAS 108, 12233 (2011).

MIT Masatoshi Imada, Atsushi Fujimori, and Yoshinori Tokura, Rev. Mod. Phys. 70, 1039 (1998).

Simons J.P.F. LeBlanc, Andrey E. Antipov, Federico Becca, Ireneusz W. Bulik, Garnet Kin-Lic Chan, Chia-Min Chung, Youjin Deng, Michel Ferrero, Thomas M. Henderson, Carlos A. Jiménez-Hoyos, E. Kozik, Xuan-Wen Liu, Andrew J. Millis, N.V. Prokof’ev, Mingpu Qin, Gustavo E. Scuseria, Hao Shi, B. V. Svistunov, Luca F. Tocchio, I. S. Tupitsyn, Steven R. White, Shiwei Zhang, Bo-Xiao Zheng, Zhenyue Zhu, and Emanuel Gull (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 5, 041041 (2015).

Stanford S. Raghu, S. A. Kivelson, and D. J. Scalapino, Phys. Rev. B 81, 224505 (2010).

Stripe Bo-Xiao Zheng, Chia-Min Chung, Philippe Corboz, Georg Ehlers, Ming-Pu Qin, Reinhard M. Noack, Hao Shi, Steven R. White, Shiwei Zhang, Garnet Kin-Lic Chan, Science 358, 1155 (2017).

key Nguyen N. T. Nam and Mikito Koshino, Phys. Rev. B 96, 075311 (2017).

TB1 Pilkyung Moon and Mikito Koshino, Phys. Rev. B 85, 195458 (2012).

TB2 E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B 82, 121407(R) (2010).

TB3 G. Trambly de Laissardière, D. Mayou and L. Magaud, Nano Lett. 10, 804 (2010).

TB4 G. Trambly de Laissardière, D. Mayou, and L. Magaud, Phys. Rev. B 86, 125413 (2012).

TB5 A. O. Sboychakov, A. L. Rakhmanov, A. V. Rozhkov, and Franco Nori, Phys. Rev. B 92, 075402 (2015).

DFT1 Shiang Fang and Efthimios Kaxiras, Phys. Rev. B 93, 235153 (2016).

DFT2 Kazuyuki Uchida, Shinnosuke Furuya, Jun-Ichi Iwata, and Atsushi Oshiyama, Phys. Rev. B 90, 155451 (2014).

DFT3 Sylvain Latil, Vincent Meunier, and Luc Henrard, Phys. Rev. B 76, 201402(R) (2007).

Mele2 E. J. Mele, Phys. Rev. B 81, 161405(R) (2010).

Neto J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).

Neto2 J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 86, 155449 (2012).

Koshino0 Pilkyung Moon and Mikito Koshino, Phys. Rev. B 87, 205404 (2013).

SM See Supplementary Material for eigenstates at Γ point and consistency check for other types of tight-binding models.

LF Liang Fu, Phys. Rev. Lett. 103, 266801 (2009).

Cao3 Y. Cao, J. Y. Luo, V. Fatemi, S. Fang, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Phys. Rev. Lett. 117, 116804 (2016).

George Serge Florens and Antoine Georges, Phys. Rev. B 70, 035114 (2004).

Kotliar S. Florens, A. Georges, G. Kotliar, and O. Parcollet, Phys. Rev. B 66, 205102 (2002).

badmetal Edward Perepelitsky, Andrew Galatas, Jernej Mravlje, Rok Žitko, Ehsan Khatami, B. Sriram Shastry, and Antoine Georges, Phys. Rev. B 94, 235115 (2016).

Zaanen Louis Felix Feiner, Andrzej M. Oleś, and Jan Zaanen, Phys. Rev. Lett. 78, 2799 (1997).

Zhang Y.Q. Li, M. Ma, D.N. Shi, and F.C. Zhang, Phys. Rev. Lett. 81, 3527 (1998).

spinliquid P. Corboz, M. Lajkó, A. M. Lauchli, K. Penc, F. Mila, Phys. Rev. X 2, 041013 (2012).


Vafek
Jian Kang and Oskar Vafek, arXiv:1805.04918.

Koshino
Mikito Koshino, Noah F. Q. Yuan, Takashi Koretsune, Masayuki Ochi, Kazuhiko Kuroki, Liang Fu, arXiv:1805.06819.

Senthil
Hoi Chun Po, Liujun Zou, Ashvin Vishwanath, T. Senthil, arXiv:1803.09742.










§ SUPPLEMENTARY MATERIAL FOR “MODEL FOR METAL-INSULATOR TRANSITION IN GRAPHENE SUPERLATTICES AND BEYOND"




§ POSSIBLE TIGHT-BINDING MODELS

In this section we list all possible tight-binding models and show that only honeycomb lattice with degenerate p_x,y orbitals can reproduce the band symmetry at all high symmetry points.

First, when the centers of Wannier orbitals are not at high symmetry points of the superlattice, the number of lowest bands should be multiples of 6 due to the point group D_3, while we only have four lowest bands from numerical calculations. Thus the Wannier centers should be at high symmetry points under D_3={C_2y^mC_3z^n|m,n∈ℤ}. When Wannier centers are along C_2y rotation axis but not at C_3z rotation invariant points (the so-called saddle points, SP for short), they form a kagome lattice, which contains 3 superlattice sites in each supercell. The number of lowest bands in kagome lattice will be multiples of 3, which is also not consistent with four lowest bands from numerical calculations.

When Wannier centers are not along C_2y rotation axis but at C_3z rotation invariant points (namely AB and BA spots), they will form a honeycomb lattice, whose number of lowest bands is multiples of 2. Notice that AB and BA spots are invariant under C_3z followed by a superlattice translation. When Wannier centers are both along C_2y rotation axis and at C_3z rotation center (namely AA spots), the triangular lattice is formed and the number of lowest bands can be any positive integer. The kagome, honeycomb and triangular lattices are shown in Fig. <ref>.



We summarize the analysis of honeycomb and triangular lattices in Table. <ref>. Since in D_3 point group the angular momentum is defined modulo 3, under D_3 higher order orbitals will transform the same as (p_x ,p_y) orbitals, s or p_z orbital, depending on the details of the orbitals. For example, (d_x^2-y^2,d_xy) orbitals are the same as (p_x ,p_y) orbitals under D_3. In other words, the analysis of (p_x ,p_y) orbitals, s and p_z orbital covers all types of orbitals. For honeycomb lattice with (p_x ,p_y) orbitals the analysis is given in the maintext and the band symmetry at all high symmetry points can be correctly reproduced. Any other orbital or lattice will not reproduce the band symmetry at all high symmetry points as shown in Table. <ref>. Some combinations cannot reproduce Γ point while some fail at K point.

For honeycomb lattice with s or p_z orbital, denote a_i as the annihilation operator at site i, the Bloch state at momentum k can be constructed as

    a_ k=∑_ie^i k· R_ia_i,

where R_i is the coordinate of honeycomb lattice site i.
We find under C_3z, the Bloch state at Γ point transforms as

    a_Γ≡∑_ia_i→∑_i'a_i'
    =a_Γ,   R_i'=C_3z R_i.

Hence the eigenstate at Γ point will be a singlet, which is not consistent with the doublet representation found in numerical band structure.

For triangular lattice with (p_x ,p_y) orbitals, the Bloch state is

    c_ k,τ=∑_ie^i k· R_ic_iτ

where τ=± denotes the angular momentum and hence the type of Wannier orbital, R_i is the coordinate of triangular lattice site i, c_iτ annihilates an electron at site i with orbital τ. We find under C_3z, the Bloch state at + K point transforms as

    c_ K,τ→∑_ie^i K· C_3z R_ie^-2iτπ/3c_iτ=∑_ie^i K· R_ie^-2iτπ/3c_iτ=e^-2iτπ/3c_ K,τ

Hence the eigenstates at K point will be two doublets {E,E}, which is not consistent with the representation {A_1,A_2,E} from microscopic considerations. Notice that at + K point, the Bloch wave function in triangular lattice is invariant under C_3z, while the Bloch wave function in honeycomb lattice has finite angular momentum L_z=± 1 under C_3z, as shown in Fig. <ref>.



For triangular lattice with s or p_z orbital, denote a_i as the annihilation operator at site i, the Bloch state at momentum k can be constructed as

    a_ k=∑_ie^i k· R_ia_i,

where R_i is the coordinate of triangular lattice site i.
We find under C_3z, the Bloch state at Γ point transforms as

    a_Γ≡∑_ia_i→∑_i'a_i'
    =a_Γ,   R_i'=C_3z R_i.

Hence the eigenstate at Γ point will be a singlet, which is not consistent with the doublet representation found in numerical band structure.





§ ILLUSTRATIVE EXAMPLES OF EIGENSTATES AT Γ POINT

Here we show illustrative examples of energy eigenstates Ψ_e,h^Γ at Γ point with L_z=-1,+1,0. In all these examples there is accidental particle-hole symmetry and states at A sublattice always dominate.


We find that the lowest two eigenstates Ψ_e,h^Γ both have angular momentum L_z=+1 if we choose

    T_0=v_F/5√(3)a[    1    1;    1 -0.6 ],
    T_±=v_F/5√(3)a[      1    ω^±;    ω^± 0.6ω^∓ ],

where ω =e^2iπ/3. The eigenenergies and real-space probability distribution of Ψ_e,h^Γ are shown in Fig. <ref>. It can be found that

    Ψ^Γ_e   =   α U^A( R_c)+β U^B(- R_c)+α^*L^A(- R_c)+β^*L^B( R_c),
    Ψ^Γ_h   =   α U^A( R_c)+β U^B(- R_c)-α^*L^A(- R_c)-β^*L^B( R_c)

where (|α|,|β|)=(0.2770,0.1643), and (argα,argβ)=(0.5250,0.0214). The eigenenergies and real-space probability distribution of Ψ_e,h^Γ of this case are shown in Fig. 4 of the maintext.

When we choose

    T_0=v_F/5√(3)a[   1   1;   1 0.6 ],
    T_±=v_F/5√(3)a[   1   1;   1 0.1 ],

we find that the lowest two eigenstates Ψ_e,h^Γ both have angular momentum L_z=-1 and

    Ψ^Γ_e   =   α U^A(- R_c)+β U^B( 0)+α^*L^A( R_c)+β^*L^B( 0),
    Ψ^Γ_h   =   α U^A(- R_c)+β U^B( 0)-α^*L^A( R_c)-β^*L^B( 0)

where (|α|,|β|)=(0.2943,0.2065), and (argα,argβ)=(0.1566,0.2446). 

We note that according to Eq.(<ref>), the component of Bloch wavefunction Ψ^Γ on the A sublattice has its maxima at AB spots on layer 1 and at BA spots on layer 2, while the component on the B sublattice has its maxima around AA spots on both layers. Following the same symmetry argument in the maintext we can understand that maxima of A sublattice component should be at AB and BA spots.
In contrast, B sublattice sites are not at rotation centers, and the intra-unit-cell phase factors e^i  K_1 · x^B and e^i  K_2 · y^B already carry the angular momentum -1. As a result, the envelope function of B sublattice wavefunction carries zero angular momentum and is therefore allowed to be peaked at AA spots.

The maxima of envelope part of Ψ_e^Γ and Ψ^Γ_h form a honeycomb superlattice consisting of two layers of corrugated triangular superlattices in both cases of L_z=± 1. Around each maximal spot, the state is dominated by states from central A sites and states from surrounding B sites, resulting in a three-fold structure. When L_z=+1, the AB bonds will be pointing away from three-fold rotation centers. When L_z=-1, the AB bonds connecting probability centers between A-site and B-site states will be pointing towards three-fold rotation centers.



The case with L_z=0 can be realized with

    T_0=v_F/5√(3)a[   1   1;   1 0.6 ],  
    T_±=v_F/5√(3)a[ 0.1   1;   1   1 ],

and is shown in Fig. <ref>. And the wavefunctions are found as following

    Ψ^Γ_e    =   α U^A( 0)+β U^B( R_c)+α^*L^A( 0)+β^*L^B(- R_c)
    Ψ^Γ_h    =   α U^A( 0)+β U^B( R_c)-α^*L^A( 0)-β^*L^B(- R_c).

where (|α|,|β|)=(0.2450,0.3949), and (argα,argβ)=(0.2329,0.2094).

In the case of L_z=0, we find a triangular lattice is formed by envelope peaks at AA spots. Thus it is tempting to construct a two-orbital tight-binding model on the triangle lattice. However, it turns out that the band symmetries at high symmetry points cannot be correctly obtained in this case since we will obtain s and p_z orbitals at each site of the triangular lattice, and obstructions may be encountered in constructing tight-binding models.






