Magic-angle Twisted Bilayer Systems with Quadratic-Band-Touching: Exactly Flat Bands with High-Chern Number

Studies of twisted moir\'e systems have been mainly focused on two-dimensional (2D) materials such as graphene with Dirac points and transition-metal-dichalcogenide so far. Here we propose a twisted bilayer of 2D systems which feature stable quadratic-band-touching points and find exotic physics different from previously studied twisted moir\'e systems. Specifically, we show that exactly flat bands can emerge at magic angles and, more interestingly, each flat band exhibits a high Chern number ($C=\pm 2$). We further consider the effect of Coulomb interactions in such magic-angle twisted systems and find that the ground state supports the quantum anomalous Hall effect with quantized Hall conductivity $2\frac{e^2}{hc}$ at certain filling. Furthermore, the possible physical realization of such twisted bilayer systems will be briefly discussed.

Introduction: Twisted moiré systems, especially twisted bilayer graphene (TBG), have attracted enormous attention in recent years due to the emergence of topological flat bands and various interesting phases such as correlation insulators and unconventional superconductivity . Since its experimental discovery, extensive studies of such systems have been done on both experimental and theoretical sides. The theoretical prediction of flat bands in TBG was made by Bistritzer and MacDonald (BM) [27], in their paper the BM Hamiltonian and the moiré band theory were developed to study TBG and other twisted moiré systems. Furthermore, a generalization of the BM model was developed [28] and a more complete description and understanding of the flat bands in a twisted bilayer system were obtained through perturbation theory. Based on the moiré band theory, enormous numbers of studies were done to explore the topological features [29][30][31][32] as well as the interaction effects  of TBG systems; non-trivial topology of the flat bands has been shown, and huge progress has been made in understanding the interacting phases.
In this paper, we investigate a twisted bilayer of systems with quadratic band touching, with focus on the twisted bilayer checkerboard (TBCB) model. The checkerboard lattice model in each single layer was proposed by Sun, Yao, Fradkin, and Kivelson (SYFK) [97] to realize a stable quadratic band touching point (QBTP) Note that for AB-bilayer graphene, the putative QBTP is not stable in the presence of trigonal hopping [110].
We found that such twisted systems can host two exactly flat bands per spin in the chiral limit, and more interestingly, each flat band has nontrivial topology with high Chern number C = ±2. Note that, in contrast to TBG with a total of eight flat bands, there are only four flat bands in the TBCB. In the presence of Coulomb interactions, by projecting them onto the topological flat bands of C = ±2 in TBCB systems similar to the analysis employed for TBG [72,73], we showed that the interaction prefers the ground state with minimum Chern number; at charge-neutrality (ν = 0) the ground state is an insulator with Chern number C = 0, while for ν = ±1 the ground state possesses Chern number ±2 and exhibits the quantum anomalous Hall effect [5,7,16,111]. We further propose a possible optical-lattice realization of the TBCB with topological flat bands, providing a promis- ing route to study various correlated phases in TBCB experimentally.
Quadratic-band-touching model: One prototype model hosting stable quadratic-band-touching points (QBTPs) is the checkerboard model proposed by SYFK [97]. As shown in Fig. 1(a), it can be described by the tight-binding Hamiltonian: H = − i,j t ij c † i c j , with the hopping amplitude t ij between sites i and j. Here we consider nearest-neighbor hopping t and next-nearestneighbor hopping ±t , as shown in Fig. 1(a). Note that the lattice consists of two sublattices labeled as A and B. By performing the Fourier transformation where d z (k) = −4t cos (k x /2) cos (k y /2) and d x (k) = 2t (cos k x − cos k y ) [112]. It is straightforward to obtain the dispersion of two bands: k,± = ± d 2 x (k) + d 2 z (k). By expanding the periodic Bloch Hamiltonian around M = (π, π) where two bands cross (namely k → k + M) and keeping only the lowest orders in k, we obtain The dispersion around M is quadratic and M is a called the quadratic-band-touching point (QBTP). To transform the Hamiltonian into a form with explicit chiral symmetry, we can perform a basis transformation ψ → U ψ with U = e i π 2 σ x and obtainH 0 (k) = U † H 0 (k)U asH Note that the QBTP features a Berry phase of 2π, which is twice of that of a Dirac point. Hereafter, unless stated otherwise, we shall assume t = t/2, so that the dispersion E ± (k) = ± t 2 k 2 . Exactly flat bands at magic angles: It is desirable to investigate novel physics in a twisted bilayer of systems with quadratic-band-touching fermions for the exotic property of the QBTP [97]. In this paper, we consider the twisted bilayer of the checkerboard lattice, and explore its physics such as totally flat bands at magic angles and the high Chern number of those flat bands. The lattice structure of the TBCB lattice is shown in Fig. 1(b).
Here we mainly focus on the low-energy physics of the twisted bilayer system with quadratic band touching by employing the continuum model describing the low-energy band structure around the QBTP M. Using the moiré band theory introduced by Bistritzer and Mac-Donald [27], we obtained the inter-layer hopping matrices: , where G 1 , G 2 are the reciprocal vectors of the lattice, α and β labels the sublattice indices A and B, respectively, M θ represents the rotation by angle θ, and τ α(β) represents the relative coordinates of the sublattice α(β) in the unit cell. For the checkerboard lattice shown in Fig.  1(a), we have τ A = (0, 0) and τ B = ( 1 2 , 1 2 ) in the unit of lattice constant a. Inspired by the TBG theory, we only keep the largest four t p terms, i.e, the terms with With these four hoppings we can construct the moiré Brillouin zone (mBZ) as shown in Fig. 2(b) and the hopping matrices take the form: where T 1 (T 2 ) is the hopping matrix of hopping along the green (brown) lines in the mBZ as shown in Fig. 2(b) [113]. Assuming the chiral limits w AA = 0 and w AB = w [114], we numerically computed the moiré bands and observed exactly flat bands for a series of magic angles as shown in Fig. 3. The band structure is controlled by a single parameter α = w t k 2 θ = wa 2 8t π 2 sin 2 (θ/2) which is proportional to 1/ sin 2 (θ/2). Note that the parameter α is qualitatively different from its counterpart w sin(θ/2) for TBG [114]. As a consequence, the magic angle for the TBCB can be much larger compared with TBG. This property makes the twist angle of the TBCB system easier to be tuned experimentally.
Origin of the exactly flat bands: We now provide an analytical understanding of the origin of exactly flat bands at those magic angles. First, we perform the Fourier transformation and obtain the hopping matrices in real space T (r) = 2 n=1 2T n cos(q n · r), where q 1 = k θ √ 2 (1, 1) and a sin θ 2 . Since the system preserves the chiral symmetry when w AA = 0, we choose the basis Φ(r) = (ψ 1,A , ψ 2,A , ψ 1,B , ψ 2,B ) T , where 1 and 2 are the layer indices and A and B are the sublattice indices, such that the Hamiltonian is given by where D(r) is completely antiholomorphic: The QBTP at M with zero energy is protected by symmetry. Explicitly, the zero-energy wave function φ M satisfying Hφ M = 0 is given by If such a holomorphic function f k (z) exists for every k in the mBZ, the wavefunction of the totally flat band with momentum k is obtained. Note that ψ M satisfy the Moire bound- [114]; consequently, f k (z) must have a simple pole, and such a construction fails in general. However, when ψ M (r) has a zero point at special twist angles, nonsingular f k (z) is permitted, and the exactly flat bands from the construction above can exist. Such a special twist angle at which ψ M (r) has zeros is a so-called magic angle. Indeed, our calculation of ψ M with zero energy at the magic angle shows that ψ M (r) is zero when r is at AA stacking point, as shown in Fig. 3

(f).
Derivation of the first magic angle: Based on the requirement that at the magic angle the wave function φ M (r) = 0 for r at the AA stacking point, we can analytically derive the parameter α corresponding to the magic angle. Solving the equation Dψ M (r) = 0 perturbatively in the parameter α < 1, we obtain the spinor wave function where · · · represents higher-order terms in α and u j (r) carries momentum mq 1 + nq 2 with |m| + |n| = j.
Up to the second order u 2 , one can get the solution u 1 (r) = ∓2α [cos(q 1 ·r)+cos(q 2 ·r)] and u 2 (r) =  Another (probably more intuitive) way to derive the magic angle is to require the vanishing of the inverse effective mass of the fermions at the QBTP M which is qualitatively different from the TBG system, which only requires the vanishing of the Fermi velocity at magic angles. Requiring the vanishing of the inverse effective mass of the fermions at the QBTP M, we obtain the first magic-angle parameter α = 1 √ 12 ≈ 0.289, which is also quite close to the value obtained numerically in Fig. 3(a). Details of computing the inverse effective mass are shown in Appendix B.
High-Chern number of exactly flat bands: To analyze the topology of these flat bands, we first calculate the Wilson loop's winding number of the TBCB lattice. For the Bloch states in the moiré Brillouin zone, to define the Wilson loop, we need to first restore the periodicity of the Bloch states and thus need to introduce the extra embedding matrix V G Q,Q = δ Q−G,Q . Consequently, the Wilson loop for the TBCB lattice is calculated by where U k1,k2 = U k = (|u 1k , ..., |u N k ) with k = (k 1 , k 2 ) [29,32,115]. We assume that b 1 and b 2 are the reciprocal vectors of the moiré lattice, k = k 1 b 1 + k 2 b 2 . We keep k 1 unchanged and vary k 2 to obtain the flow of the Wilson loop spectrum along k 1 . We find that the Wilson loop winds from ±1 to ∓1 twice while k 1 goes from −π to π, so the winding of the Wilson loop at the first magic angle is ±2, as shown in the Fig. 4. This is also true for other magic angles that we identified. When adding more bands into the Wilson loop calculation such as the middle six bands, the winding is still preserved which suggests that the topology is stable. In fact, the antiunitary particle-hole symmetry P (see the Appendix C) protects the degeneracy of the Wilson bands at k 1 = 0/π in the same way as in TBG [32]. In the non-interacting limit, the degeneracy of the two flat bands for the TBCB is protected by the timereversal symmetry T . To give a simple illustration of the topology with the existence of interactions, we apply a weak0time-reversal-symmetry-breaking term which preserves all other symmetries except the mirror symmetry: Aσ y sin( kx 2 ) sin( ky 2 ). The degeneracy is then lifted, while the flatness of bands is still well preserved. Assuming A/w AB = 0.5, we have calculated the Berry curvature of the two flat bands and the corresponding Chern number: C n = 1 2π mBZ B n d kx d ky . Our result shows that the lower band hosts a Chern number of ∓2, while the upper band host a Chern number of ±2 where the sign of the Chern number depends on the sign of A. For a spinless TBCB system, this is the only possible way for the two bands to split; the topology of the system is highly nontrivial with high Chern number C = ±2. It is noteworthy that this high Chern number of ±2 is realized in the flat bands of a twisted bilayer systems [116] with stable QBTPs, and in our system, the QBTPs have only one valley, which will lead to different physics.
Correlation effect: Interactions can play an essential role in the physics of twisted bilayer systems [33-40, 42-47, 49-79]. Here, we consider the Coulomb interactions: where δρ q = r e iq·r (ρ r − 1 2 δ q,0 δ G,0 ), G represents the moiré reciprocal lattice vectors in the Brillouin zone (BZ) of the original lattice, and V (q) = πd 2 U d tanh(d|q|/2) d|q|/2 is the screened Coulomb potential with U d = e 2 /( d) [72]; d is the distance between the TBG and the top or bottom gate, and ρ r is the charge density at r. To solve its low-energy physics, one can project the Hamiltonian onto the subspace of the two flat bands. To do so, we employ fermion operators in the mBZ energy band basis c † n,s (k) = Qα u Q,α,n (k)f † α,s (k + Q), where Q ∈ Q ± and Q ± is the collection of the sites of layer l = ± of the mBZ as plotted in Fig. 2. Here, n is the moiré band index, and n = ±1 represents the two flat bands. Due to its nontrivial topology, we cannot define a symmetric smooth and periodic wave function u Q,α,n (k) [68]. Here we adopt a periodic gauge that satisfies u Q,α,n (k + b i ) = u Q−G,α,n (k). Similar to the treatment of interactions in TBG [72,73], after the projection into the flat band subspace the interacting Hamiltonian for the TBCB system is written as where Ω tot is the total area of the TBCB system, and M mn (k, q+G) = α,Q∈Q± u * Q−G,α,m (k + q)u Q,α,n (k). (11) As discussed in the Appendix C, considering the C 2z , T , and P symmetries, the M matrix is constrained as follows: M (k, q + G) = ζ 0 α 0 (k, q + G) + iζ y α 2 (k, q + G), (12) where ζ represents the Pauli matrix for the two-flat-band subspace, α 0 (k, q+G) and α 2 (k, q+G) are real numbers with the constraints α 0 (k, q+G) = α 0 (k+q, −q−G) and α 2 (k, q + G) = −α 2 (k + q, −q − G) and α a (k, q + G) = α a (−k, −q − G) for a = 0, 2 (see details in the Appendix C).
In the Chern band basis [32,72] which is the eigenstate of the flat bands with Chern number C = 2e: , where e = ±1, we can rewrite the operator O q,G in a diagonal way where M e (k, q + G) = α 0 (k, q + G) + ieα 2 (k, q + G). ν + 2 = ν + + ν − is the total filling factor of the system with ν ± being the integer filling of the Chern bands with Chern number C = ±2. It is clear that the state Ψ ν+,ν− ν carries a Chern number of 2(ν + −ν − ) and different states with the same ν are degenerate.
In real TBCB systems, it is difficult to tune the intrasublattice hopping w AA to be strictly zero. When w AA = 0, neither the particle-hole symmetry nor the chiral symmetry holds (see Appendix C for details). Thus the M matrix takes a general form In the Chern basis, the operator O q, where F e (k, q + G) = ieα 1 (k, q+G)+α 3 (k, q+G). When the system has an even filling factor ν and assuming the flat band approximation, the ground state of the operator , which has a zero Chern number. For odd filling factors ν, after taking O 1 q,G as a perturbation, the degeneracy of the different Chern states will be lifted, and O 1 q,G prefers the ground state with minimum Chern number. For instance, when ν = ±1, the system possesses a ground state with broken time-reversal symmetry and high Chern number C = ±2.
Discussions and concluding remarks: In this paper, we proposed a twisted bilayer system of fermions with C 4 -symmetry-protected quadratic band touching, which can exhibit exactly flat bands with high Chern numbers C = ±2. Our system's symmetry and the stable QBTPs are noteworthy aspects of this study of the twisted graphene system. The origin of the exactly flat band is related to the anti-holomorphic property of the Hamiltonian in the chiral limit with t = 2t . At the first magic angle, the flatness of the topological bands is rather robust against deviation from the exactly flat conditions; that is, the topological bands in the middle exhibiting a high Chern number of ±2 are nearly flat for a wide range of parameters. See details in the Appendix A.
Such a TBCB lattice may be realized by loading cold atoms into a specially-designed optical-lattice system [117,118]. It has been proposed that the twisted square lattice can be realized by introducing four states (labeled by spin ±1/2 and two layers) and constraining each "layer" by a set of square optical lattices that differ by polarization and a small twisting angle [119]. If 2π fluxes are added to the square plaquettes such that the hopping amplitude along its diagonal is −t , the TBCB lattice maybe experimentally realized, and the quantum anomalous Hall effect associated with the flat band with high Chern number may be observed. Furthermore, away from integer filling, it is also possible to realize interesting phases such as unconventional superconductivity and fractional Chern insulators, which are left for future studies.
Acknowledgments As discussed in the main text, the exactly flat band criteria of the TBCB require t = 2t and w AA = 0. In this appendix, we show how the flatness of the two flat bands is affected by these two parameters near the first magic angle. We calculated the bandwidth of the middle two flat bands while varying t or w AA . The results are shown in Fig. S1. Notice that if w AA = 0, the particlehole symmetry is broken, and thus the flat bands deviate from zero energy. As shown in Fig. S1(a), as t varies from 500 meV to 1500 meV, the bandwidth of the flat bands varies from 0 to 2.5 meV which is relatively small. As w AA varies from 0 to 2.0 meV [ Fig. S1(c)], the bandwidth of the flat bands varies from 0 to 3 meV. The flat bands are quite robust against the deviation of the parameters.

Appendix B. Another Derivation of the First Magic Angle
We now provide another derivation of the first magic angle by requiring the vanishing of the inverse effective mass of the fermions at the QBTP. (Note that this is in contrast to the magic-angle definition of vanishing Fermi velocity in TBG.) Here we only consider the nearest four M points of the bottom layer (red) to the center point of the top layer (blue) in the mBZ, as shown in Fig. 2 of the main text. We write the ten-band Hamiltonian in the momentum space for these five points in the mBZ: where 4 ) satisfies the Schrödinger equation: From Eq. (S2b) we obtain that ψ b,i = (E − H 0 (q i )) −1 T i ψ t , from which we can get the effective Schrödinger equation for ψ t : (S3) Neglecting the E 2 and k n for n > 2 terms as small and noticing that q 2 i = 1, we get Substituting the T i we have obtained before, one can get the effective Hamiltonian with t eff = 1−12α 2 1+4α 2 t . When α = 1 The Hamiltonian of the twisted bilayer checkerboard lattice is where h lθ/2 (k) is the kinetic term of the checkerboard lattice with a twist angle lθ/2 from the x axis [l = +1 (l = −1) for the upper (lower) layer] and has the form of Eq.(1) in the main text. Let σ and τ represent the Pauli matrix for the sublattice degrees of freedom and the layer degrees of freedom respectively. The Hamiltonian with the moiré BZ as shown in Fig. 2(b) respects the following spatial symmetries.
C 2x/y symmetry The Hamiltonian also processes the particle-hole symmetry and time-reversal symmetry Particle-hole symmetry (S10) Notice that here the particle-hole symmetry is a rigorous one but will be broken when w AA = 0 which is different from TBG. The system also preserves chiral symmetry if w AA = 0 with the operator: σ y T . With these symmetries, we can fix the gauge of the wave function. We introduce the sewing matrix B g (k) for the C 2z , T and P symmetries.
[D (P )] u n (k) = B P (k) m,n u m (−k). (S13) The sewing matrix can be simplified as B C2z (k) m,n = δ m,n e iϕ C 2z (k) n , B T (k) m,n = δ m,n e iϕ T n (k) , B P (k) m,n = δ −m,n e iϕ P n (k) . (S14) These three symmetry operators can be combined to obtain two independent symmetry operations C 2z P and C 2z T which keep k unchanged. The corresponding sewing matrices are defined by the following equations (S15) The symmetry operations C 2z P and C 2z T satisfy the properties Thus we can adopt the following k-independent sewing matrices B C2zT (k) m,n = δ m,n , B C2zP (k) m,n = − sgn(n)δ −m,n .
(S17) These sewing matrices can also be expressed by the Pauli matrix for the two flat bands where ζ represents the Pauli matrix for the two-flat-band subspace. We have chosen a similar form to the sewing matrix of the TBG systems adopted in Ref. [72], and the difference is that TBCB systems do not have two valleys. The wave function and thus the M matrix introduced in the main text and Appendix D, M mn (k, q+G) = α,Q∈Q± u * Q−G,α,m (k + q)u Q,α,n (k), (S19) are also constrained by the two symmetries, C 2z T and C 2z P , with the sewing matrices we obtained in Eq. (S18), Thus the M matrix takes the form where α 0 (k, q + G) and α 2 (k, q + G) are real numbers. Besides, from the definition of the M matrix in the main text, the M matrix also satisfies the Hermiticity condition which means that α i (k, q + G) satisfy α 0 (k, q + G) = α 0 (k + q, −q − G) We can also fix the relative gauge between the wave functions with momentum k and those with momentum −k by C 2z symmetry. Notice that in the TBCB system, C 2z , T , and P symmetries commute with each other, and we can choose the sewing matrix for these three symmetries Thus, the M matrix also has the following constraint implied between the momentum k and the momentum −k: which implies that When the hopping w AA = 0, the particle and chiral symmetries are broken, and Eq. (S20) no longer holds. Constrained by the real condition, the M matrix takes a more general form, M (k, q + G) = ζ 0 α 0 (k, q + G) + ζ x α 1 (k, q + G) + iζ y α 2 (k, q + G) + ζ z α 3 (k, q + G) = M 0 (k, q + G) + M 1 (k, q + G), (S27) where α i (k, q + G) (i = 0, 1, 2, 3) are real numbers.
Similar to the chiral case introduced above, now α i (k, q + G) are also constrained by the Hermiticity condition and the C 2z symmetry: α a (k, q + G)α a (k + q, −q − G) for a = 0, 1, 3, where M m,n (k, q + G) = α Q∈Q± u * Q−G,α,m (k + q)u Q,α,n (k). (S32) The interacting Hamiltonian now is written in a semi-positive definite form where O q,G = ks m,n=±1 V (q + G)M m,n (k, q + G) c † m,s (k + q)c n,s (k) − 1 2 δ q,0 δ m,n . Notice that the number of the electron N is conserved; thus we are able to introduce a Lagrange multiplier A G When the flat metric condition [72,73] M m,n (k, G) = ξ(G)δ m,n is satisfied or filling factor ν = 0, the last two terms in Eq. (S34) are constant which depend on N . In this way, one can easily conclude that the ground state of the interacting Hamiltonian satisfies the equation To solve the ground state, one only needs to solve Eq. (S35). In general, the Flat Metric Condition is not strictly satisfied except for G = 0. Fortunately, when the flat metric condition is not largely violated, the ground states which satisfy Eq. (S35) persist as long as the gap between the ground states and exciting states is not closed. Since the wave function decreaseS exponentially as G increases for the moiré Hamiltonian [28], one can assume that the flat metric condition is not largely violated and the ground state derived above is the real ground state of the system. Future studies can adopt the real-space projection method [50,122] to confirm our conclusion for ground states.