The rich physics of the experimentally observed insulating states in magic angle twisted bilayer graphene (TBG) at integer number of electrons per unit cell and the superconducting phase with finite doping above the insulating states has attracted considerable interest . The single-particle picture predicts a gapless metallic state at electron number ±(3, 2, 1), and hence the insulating states have to follow from many-body interactions. The initial observations of the insulating states [2][3][4][5] were then followed by the experimental discovery by both scanning tunneling microscope [20,21] and transport [6,11,[22][23][24][25] that these states might exhibit Chern numbers, even when the TBG substrate is not aligned with hBN, which would indicate a many-body origin of the Chern insulator.
These remarkable experimental advances have been followed by extensive theoretical efforts aimed at their explanation . Using a strong-coupling approach where the interaction is projected into a Wannier basis, Kang and Vafek [71] constructed a special Coulomb Hamiltonian, of an enhanced symmetry, where the ground state (of Chern number 0) at ±2 electrons per unit cell can be exactly obtained (with rather weak assumptions). In Ref. [110], we have showed that the type of Kang-Vafek type Hamiltonians [71] (hereby called positive semidefinite Hamiltonians -PSDH) are actu-* bernevig@princeton.edu † zhidas@princeton.edu ally generic in projected Hamiltonians, and that the presence of extra symmetries [38,63,71] renders some Slater determinant states to be exact eigenstates of PSDH. We found at zero filling, these states are the ground states of PSDH. At nonzero integer filling, these states are the ground states of the PSDH under weak assumptions (first considered by Kang and Vafek [71]). With a unitary particle-hole (PH) symmetry first derived in Ref. [43], the PSDH projected to the active bands has enhanced U(4) (in all the parameter space) and U(4) × U(4) (in a certain, first chiral limit) symmetries first mentioned in Refs. [71][72][73]. We showed [109,110] that these symmetries are valid for PSDHs of TBG irrespective of the number of projected bands. We also found that, for two projected bands in the first chiral limit (a second chiral limit, of U(4) × U(4) defined in Ref. [109] was also found), ground states of different Chern numbers are exactly degenerate [110]. These ground states are all variants of U(4) ferromagnets (FM) in valley/spin. When kinetic energy is added or away from the chiral limit, the lowest/highest Chern number becomes the ground state in low/high magnetic field, which explains/is consistent with experimental findings [6,11,[20][21][22][23][24][25].
In this paper, we show that the Kang-Vafek type of PSDH also allow, remarkably, for an exact expression of the charge ±1 excitation (relevant for transport gaps) energy and eigenstate, neutral excitation (relevant for the Goldstone and thermal transport), and charge ±2 excitation (relevant for possible Cooper pair binding energy). We show that the charge excitation dispersion is fully governed by a generalized "quantum geometric tensor" of the projected bands, convoluted with the Coulomb interaction. The smallest charge 1 excitation gap is at the Ŵ M point. The neutral, and charge ±2 excitation, on top of every FM ground state can also be obtained as a single-particle diagonalization problem, despite the state having a thermodynamic number of particles. The neutral excitation has an exact zero mode, which we identify with the FM U(4)-spin wave, and whose low-momentum dispersion (velocity) can be computed exactly. The charge ±2 excitations allows for a simple check of the Richardson criterion [112][113][114][115] of superconductivity: we check if states appear below the noninteracting two-particle continuum. We find a sufficient criterion for the appearance/lack of Cooper binding energy in these type of PSDH Hamiltonian systems based on the eigenvalues of the generalized "quantum geometric tensor." We analytically show that, generically, the projected Coulomb Hamiltonians cannot exhibit Cooper pairing binding energy. As such, this implies that either phonons or nonzero kinetic energy are needed for superconductivity. Since the Ref. [103] showed that the kinetic energy bounds on the superexchange energy are less 10 -3 in Coulomb units, the phonon mechanism becomes becomes likely. If however, experimentally, the kinetic energy is stronger, a Coulomb mechanism for superconductivity is still possible. Since we proved that flat bands cannot Cooper pair under Coulomb, a prediction of a Coulomb with nonflat bands mechanism for superconductivity would be that the highest superconducting temperature does not happen at the point of highest density of states DOS. This is in agreement with recent experimental data [25].
We generically consider the TBG system with a Coulomb interaction Hamiltonian projected to the active eight lowest bands (2 per spin-valley flavor) obtained by diagonalizing the single-particle Bistritzer-MacDonald (BM) [1] TBG Hamiltonian (see Appendix A 1 for a brief review, and more detail in Refs. [107,108]). The projected single-particle Hamiltonian reads
where we define η = ± for graphene valleys K and K ′ , s = ↑, ↓ for electron spin, and n = ±1 for the lowest conduction/valence bands in each spin-valley flavor. c † k,n,η,s is the electron creation operator of energy band n, with the origin of k chosen at Ŵ point of the moiré Brillouin zone (MBZ).
The density-density Coulomb interaction, when projected into the active bands of Eq. (1), always takes the form of a positive semidefinite Hamiltonian (PSDH) (see proof in Ref. [109], see also brief review in Appendix A 2):
where tot is the sample area, and G runs over all vectors in the (triangular) moiré reciprocal lattice Q 0 . This Hamiltonian is of a same positive semidefinite form as that Kang and Vafek [71] obtained by projecting the Coulomb interaction into the Wannier basis of the active bands. In this work, we will omit the kinetic energy. Due Ref. [110], the energy splitting between the degenerate ground states of Eq. ( 2) is smaller than 0.1meV per electron. As shown in the rest of this work, the characteristic energy of charged and neutral excitations is about 10 meV. Thus it is safe to neglect the kinetic energy for most of the excitations. But some of the U(4) Goldstone modes might be opened a small gap due to the kinetic energy. We leave this effect of kinetic energy to future studies.
The O q,G operator takes the form O q,G = k,m,n,η,s
V (G + q)M (η) m,n (k, q + G) × ρ η k,q,m,n,s -
where V (q) is the Fourier transform of the Coulomb interaction, ρ η k,q,m,n,s = c † k+q,m,η,s c k,n,η,s is the density operator in band basis, and the -1 2 δ q,0 δ m,n factor is a chemical potential added to respect many-body charge conjugation symmetry (see Appendix A 2 a and Ref. [109]). For theoretical derivations, we shall keep V (q) general except that we assume V (q) 0 and only depends on q = |q|; although for numerical calculations we will take V (q) = 2π e 2 ξ tanh(qξ /2)/ǫq for dielectric constant ǫ(∼6) and screening length ξ (∼10 nm) (see Appendix A 2). In particular, O -q,-G = O † q,G , and thus H I in Eq. ( 2) is a PSDH. An important quantity in Eq. ( 3) for our many-body Hamiltonian are the form factors, or the overlap matrices, of a set of bands m, n (Appendix A 2 a)
where u Qα;nη is the Bloch wave function of band n and valley η (here α = A, B denotes the microscopic graphene sublattices, and Q are sites of a honeycomb momentum lattice with definition in Appendix A 1, see also Ref. [107] for details). A nonzero Berry phase of the projected bands renders the spectra of the PSDH Eq. ( 2) not analytically solvable: the O q,G 's at different q, G generically do not commute (unless in the stabilizer code limit discussed in Refs. [109,110]), and hence the PSDH is not solvable. The properties of the PSDH Eq. ( 2) depend on the quantitative and qualitative (symmetries) properties of the form factors in Eq. ( 4), which are detailed in Refs. [107][108][109] and briefly reviewed in Appendix A 2 a. First, in Ref. [107], we showed that M (η) m,n (k, q + G) falls off exponentially with |G|, and can be neglected for |G| > √ 3k θ , where k θ = 2|K| sin(θ /2) is the distance between the K points of two graphene sheets. Furthermore, we showed in Refs. [108,109] that by gauge fixing the C 2z , T , and unitary particle-hole symmetry P [43], the form factors can be rewritten into a matrix form in the n, η basis as [see Eq. (A11)] M (η) mn (k, q + G) = 3 j=0 (M j ) m,η;n,η α j (k, q + G), (5) where M 0 = ζ 0 τ 0 , M 1 = ζ x τ z , M 2 = iζ y τ 0 , and M 3 = ζ z τ z , and α j (k, q + G) are real scalar functions satisfying [Eqs. (A12) and (A13) in Appendix A 3]. A further simplification [72] happens in a region of the parameter space where the AA interlayer coupling w 0 = 0 205415-2 [72], which is called the (first) chiral limit [37] (a similar simplification occurs in a second chiral limit [109]). In this limit there is another chiral symmetry C anticommuting with the single-particle Hamiltonian, which further imposes the constraints α 1 (k, q + G) = α 3 (k, q + G) = 0 (see Ref. [110] and Appendix A 3). The first chiral limit also allows for the presence of a Chern band basis in which bands of Chern number e Y = ±1 are created by the operators
In Ref. [108], we detail the gauge-fixing for this basis. The Chern basis is also discussed in Refs. [72,74,108]. The form factors under the Chern basis take the simple diagonal form
e Y (k, q + G) = α 0 (k, q + G) + ie Y α 2 (k, q + G). (7) The symmetries of the projected Hamiltonian in the nonchiral (w 0 , w 1 = 0) and two chiral w 0 = 0 or w 1 = 0 limits are important. We will use the matrices ζ a , τ a , s a with a = 0, x, y, z as identity and x, y, z Pauli matrices in (particle-hole related) band, valley and spin-space respectively. In Ref. [109] (short review in Appendix A 2 c), we have showed that the PSDH has a U(4) symmetry in the nonchiral limit (with single-particle representations of generators s ab = {ζ y τ y s a , ζ y τ x s a , ζ 0 τ 0 s a , ζ 0 τ z s a } with a, b = 0, x, y, z in the energy band basis c † k,+,η,s ), and a U(4) × U(4) symmetry in the two chiral-flat limits (with single-particle representations of generators s ab ± = (1 ± e Y )τ a s b /2 in the Chern band basis d † k,e Y ,η,s [108], Appendix A 2 b), mirroring the results obtained by Refs. [42,[71][72][73] for projection into the two active bands. We note that in Ref. [109] we showed these symmetries hold for any number of PH symmetric projected bands. In Appendixes A 2 c and A 2 e, we provide a summary of these detailed results. Adding the kinetic term in the first chiral limit breaks the U(4) × U(4) symmetry of the projected interaction to a U(4) subset(with generators s ab = ζ 0 τ a s b in the energy band basis, (a, b = 0, x, y, z)). The symmetries we found in the first chiral and nonchiral limits agrees with that in Ref. [72], and the relation between our U(4) symmetry generators and those of Kang and Vafek [71] are given in Ref. [109]. We will restrict our study within the nonchiral-flat limit and first chiral-flat limit in this paper. Thus, without ambiguity, we will simply call the first chiral limit the "chiral limit." With these symmetries, in the nonchiral-flat limit (where the projected kinetic Hamiltonian H 0 = 0), one can write down exact eigenstates of the PSDH Eq. ( 2), which we have analyzed in full detail in Ref. [110] and review in Appendix A 3 b. In the nonchiral limit, O q,G is diagonal in η and s, and filling both n = ± bands of any valley/spin gives a Chern number 0 eigenstate for all even fillings ν = 0, ±2, ±4 (along with any U(4) rotation) [110]:
where {η j , s j } are distinct valley-spin flavors which are fully occupied. They form the [(2N M ) (ν+4)/2 ] 4 irreducible representation (irrep) of the nonchiral-flat limit U(4) symmetry group, where [λ p ] 4 is short for the Young tableau notation [λ, λ,
with 0 p 4 identical rows of length λ (see Ref. [110] for a brief review). With M (η) m,n (k, q + G) in Eq. (A11), we have that the state
where N M is the total number of moiré unit cells. For ν = 0, the state Eq. ( 8) is always a ground state as it is annihilated by O q,G [110].
In the first chiral-flat limit (where H 0 = 0 and w 0 = 0), the projected Hamiltonian Eq. ( 2) has as eigenstates, the filled band wave functions [110] (see Appendix A 3 a for brief review): (10) where ν + -ν -= ν C is the total Chern number of the state and ν + + ν -= ν + 4 (0 ν ± 4) is the total number of electrons per moiré unit cell in the projected bands, k runs over the entire MBZ and the occupied spin/valley indices {η j 1 , s j 1 } and {η j 2 , s j 2 } can be arbitrarily chosen. Moreover, these eigenstates of Eq. ( 2) are also eigenstates of O q,G in Eq. ( 3
, where A G is still given by Eq. ( 9). They form the ([ 4) (Young tableaux notation, see Ref. [110]). For a fixed integer filling factor ν, we found that the states with different Chern numbers ν C are all degenerate in the chiral-flat limit [110]. In particular, at charge neutrality ν = 0, the U(4) × U(4) multiplet of | ν + ,ν - 0 with Chern number ν C = ν + -ν -= 0, ±2, ±4 are exact degenerate ground states. At nonzero fillings ν, we cannot guarantee that the ν = 0 eigenstates are the ground states.
In Ref. [110], we found that under a weak condition, the eigenstates Eqs. ( 8) and (10) become the ground states of H I for all integer fillings -4 ν 4 (ν even in Eq. 8). If the q = 0 component of the form factor M (η) m,n (k, G) is independent of k for all G's, i.e., flat metric condition:
then all the states in Eqs. ( 8) and ( 10) become ground states of H I by an operator shift [Eq. (A29)] [71,109,110] (see Appendixs A 3 b and A 3 a). We noted in Ref. [107] that this flat metric condition is always true for G = 0, for which M (η) m,n (k, 0) = δ mn from wave-function normalization. In Ref. [107] we have shown that, around the first magic angle, M (η) m,n (k, G) ≈ 0 for |G| > √ 3k θ for i = 1, 2. Hence, the condition Eq. ( 11) is valid for all G with the exception of the 6 smallest nonzero G satisfying |G| = √ 3k θ . Hence, the condition is largely valid, and our numerical analysis [107] confirms its validity for k in a large part of the MBZ. The idea to impose a similar condition as Eq. ( 11) first used by Kang and Vafek [71] to find the ν = ±2 ground state for their PSDH. Due to a slightly different U(4) symmetry, our U(4) FM states are different, but overlap with the Kang and Vafek ones in the chiral limit, as discussed in detail in Refs. [109,110].
We note that for ν = 0, the states in Eqs. ( 8) and (10) still remain the exact ground states if the flat metric condition Eq. ( 11) is not violated too much [110,111]. This is because they correspond to gapped insulator eigenstates [110,111] when condition Eq. ( 11) is satisfied, and the flat metric condition Eq. (11) has to be largely broken to bring down another state into the ground state. From now on, we "call" Eqs. ( 8) and (10) ground states of the system.
Remarkably, as we will show in the rest of our paper below, one can analytically find a large series of excitations above the ground states Eqs. ( 8) and (10).
Our excitations will be build out of acting with the band creation and annihilation operators on the ground states in Eqs. ( 8) and (10). We first need to compute the commutators in the nonchiral Hamiltonian (see Appendix B, in particular, Appendix B 1)
where we have used the property
n,m (kq, q + G) [109]. In the chiral limit, the same operators read in the Chern basis (see Appendix B 2)
From these equations, we can obtain the commutators of O -q,-G O q,G with the band electron creation operators in the nonchiral case as
and in the first chiral limit in Chern basis as
respectively. Similar relations for [O -q,-G O q,G , c k,n,η,s ] and [O -q,-G O q,G , d k,e Y ,η,s ], where M (η) (k, q + G) → M (η) * (k, -q -G), are derived in Appendix B. The matrix factor P is the convolution of the Coulomb potential and the form factor matrices. In the nonchiral case, P is a matrix given by
In the first chiral limit, it is a number independent on e Y :
where α 0 (k, q + G) and α 2 (k, q + G) are the decomposition of the form factors in Eq. (7). The above commutators and the existence of exact eigenstates Eqs. ( 10) and (8), which are ground states with the flat metric condition Eq. ( 11), allow for the computation of part of the low energy excitations with polynomial efficiency. We now show the summary of the computation for the bands of charge +1, +2 and neutral excitations. The charge -1, -2 excitations can be found in Appendixes C 3 and E 4, respectively.
To find the charge one excitations (adding an electron into the system), we sum the commutators in Eq. ( 14) over q, G, and use the fact that the ground states in Eqs. ( 8) and (10) 9) in their corresponding limits. For any state | in Eqs. ( 8) and (10), we find
where N is the electron number operator, and the matrix
We hence see that, if | is one of the | ν + ,ν - ν Eq. (10) or one of the | ν Eq. ( 8) eigenstates of H I , then c † k,m,η,s | can be recombined as eigenstates of H I with eigenvalues obtained by diagonalizing the 2 × 2 matrix R η mn (k). In the nonchiral case, the eigenstates | ν we found in Ref. [110] (and re-written in Eq. ( 8)) have both active bands n = ±1 in each valley η and spin s either fully occupied or fully empty.
In this case, we can consider two charge +1 states c † k,n,η,s | (n = ±) at a fixed k in a fully empty valley η and spin s. These two states then form a closed subspace with a 2 × 2 subspace Hamiltonian R η (k) defined by Eq. (19). Diagonalizing the matrix R η (k) then gives the excitation eigenstates and excitation energies. Furthermore, at ν = 0, the state | ν=0 in Eq. ( 8) is the ground state of the interaction Hamiltonian H I regardless of the flat metric condition Eq. (11), and hence c † k,n,η,s | ν=0 always gives the charge excitation above the ground state.
If we further assume the flat band condition Eq. ( 11) (or its violation is small enough), all eigenstates | ν become exact ground states and the second row of Eq. ( 19) vanishes (see Appendix C 1 a). Since the U(4) irrep of the ground state
Exact charge ±1 excitations given by the simplified excitation matrix [Eq. (20)] for three different w 0 /w 1 at the twist angle θ = 1.05 • . Here we change w 0 while keeping w 1 = 110 meV fixed. Other parameters are given in Appendix A. These excitations are exact at the charge neutrality point (ν = 0) for generic state and are exact at finite integer fillings if the flat metric condition is satisfied. The charge +1 and -1 excitations are degenerate. The exact charge ±1 excitations obtained using the full excitation matrix [Eq. (19)] without assuming the flat metric condition are given in Figs. 5 and6 in Appendix C 4. The charge gap in those cases shrinks considerably.
charge -1 excitations is derived in in Appendix C 3, where we denote the excitation matrix as R.
As explained in Appendix A 3, when the flat metric condition is satisfied, the second term in R η [Eq. (19)] can be canceled by the chemical potential term (the third term), and thus we obtain a simplified expression for R η independent of ν:
It is worth noting that Eq. ( 20) is exact for ν = 0 even without the flat metric condition Eq. ( 11), because the coefficient A G [Eq. ( 9)] in the second term of R η [Eq. ( 19)] and the chemical potential in the third term of R η vanish at ν = 0. The simplified matrix R η for charge -1 excitation with the flat metric condition Eq. ( 11) is the complex conjugation of R η , i.e., R η mn (k) = R η * mn (k). This shows that, the charge -1 excitations are degenerate with the charge +1 excitations if either ν = 0 or the flat metric condition Eq. ( 11) is satisfied. The charge +1 excitation dispersion determined by Eq. ( 20) (which does not depend on ν) is plotted in Fig. 1.
The parameters used in the calculation to obtain the spectrum are given in Appendix A. We find that, with the flat metric condition imposed, the charge ±1 excitation (Fig. 1) is gapped, and the minimum is at the Ŵ point, with a large dispersion velocity. The exact charge ±1 excitations at different fillings obtained using the full R η matrix [Eq. ( 19)] of realistic parameters [which break the flat metric condition Eq. (11)] are given in Figs. 5 and6 in Appendix C 4.
The degeneracy of the excitation spectrum depends on the filling ν of the ground state. In the nonchiral-flat U(4) limit, R η does not depend on spin, and R + , R -have the same eigenvalues because they are related by the symmetry C 2z P, where P is a single-body unitary PH symmetry (App. A 2 a) [43,108,109]. Thus charge +1 excitations in different valleyspin flavors have the same energy. For the state | ν [Eq. ( 8)], the +1 excitations in the empty (4 -ν)/2 spin-valley flavors are degenerate. Correspondingly, -1 excitations in the occupied (4 + ν)/2 spin-valley flavors are also degenerate.
In the (first) chiral-flat limit, and with the flat metric condition Eq. ( 11) [or at ν = 0 without (11)], the expression for the charged excitations in the Chern basis d † k,e Y ,η,s | ν + ,ν - ν becomes diagonal and independent of e Y [see Appendix C 2 for the chiral-flat limit without the flat metric condition Eq. ( 11)]:
provided that the Chern band e Y (= ±1) in valley η and spin s is fully empty and P(k, q + G) given in Eq. (17). We obtain
The spectrum at the magic angle is shown in Fig. 1. The U(4) × U(4) irrep of the charge +1 excited states with e Y = 1 and e Y = -1 are given by ([
, respectively. The charge -1 excitation details can be found in Appendix C 3.
Since in the chiral-flat limit the scattering matrix R 0 (k) is identity in the e Y space, the excitation has degeneracy in addition to the valley-spin degeneracies. For a state in Eq. (10) with filling ν, the charge +1 and -1 excitations have degeneracies 4 -ν and 4 + ν, respectively.
In this section, we will focus on the charge neutrality point (ν = 0), where the second and third terms in R η (k) [Eq. (19)] vanish, and nonzero integer fillings ν = ±1, ±2, ±3 with the flat metric condition Eq. (11) such that the second and third terms in R η (k) cancel each other. In these cases, R η (k) is a positive semidefinite matrix and hence has non-negative eigenvalues. We are able to obtain some analytical bounds for the gap of the ±1 excitation. Detail calculations are given in Appendix C 1 b. Since charge ±1 excitations in this case are degenerate, our conclusion below for charge +1 excitations also apply to charge -1 excitations.
We rewrite the R η (k) matrix as R η mn (k) = (M (η) † (k)V M (η) (k)) mn , where now M (η) (k) with given η and k is a matrix of the dimension 2N M • N G × 2 (with 2 because we are projecting into the two active TBG bands). N M is number of moiré unit cells, N G is the number of plane waves (MBZs) taken into consideration. By separating the {q, G} = 0 contribution, and using Weyl's inequalities we find in Appendix C 1 b that the energies of the excited states are 1 2 V (q = 0)/ tot . The bound 1 2 V (q = 0)/ tot is small but nonzero for large but finite tot . This shows that the states c † k,n,η,s | are not exactly degenerate to the ground state | (note that we did not prove these are the unique ground states).
The excited states of the PSDH appears to give rise to finite gap charge 1 excitations. The largest gap happens in the atomic limit or a material, where u m (k + q)|u n (k) = δ mn , for which R mn = δ mn q,G V (q + G) = δ mn tot V (r = 0). Hence the gap is 1 2 V (r = 0). Away from the atomic limit, the gap is reduced, but will generically remain finite. We now give an argument for this. Since we know that TBG is far away from an atomic limit-the bands being topological, we expect a reduction in this gap. We perform a different decomposition of the matrix R η mn : we separate it into G = 0 and G = 0 sums (see Appendix C 1 c). The G = 0 part, besides being negligible for |G| √ 3k θ [107], is also positive semidefinite, and the eigenvalues of R η mn are bounded by (and close to) the G = 0 part:
where q is summed over the MBZ, and the inequality means that the eigenvalues of the left-hand side are equal to or larger than the eigenvalues of the right-hand side. We then rewrite the right-hand side as q V (q)(δ mn -G mn η (k, q)), where we call the positive semidefinite matrix G mn (k, q) the generalized "quantum geometric," whose trace is the generalized Fubini-Study metric. For small momentum transfer q, we can show that G mn (k, q) = i j q i q j G mn i j (k) + O(q 3 ). where G mn i j (k) is the conventional quantum geometric tensor (and the Fubini-Study metric) [77,116] defined by
in which m, n ∈ B are energy band indices and i, j are spatial direction indices of the orthonormal vectors u m (k) in a N dimensional Hilbert space, with k being the momentum (or other parameter). The G mn (k, q) tensor quantifies the distance between two eigenstates in momentum space. Generically, we expect [77] that the inner product between two functions at k and k + q to fall off as q increases, leaving a finite term in R η mn (k), the electron gap, at every k. In trivial bands in the atomic limit, the positive semidefinite matrix G mn (k, q) reaches its theoretical lower bound 0 and hence the charge 1 gap is maximal. In topological bands, such as TBG, the quantum metric has a lower bound and hence the charge 1 gap is reduced.
To obtain the charge neutral excitations, we choose the natural basis c † k+p,m 2 ,η 2 ,s 2 c k,m 1 ,η 1 ,s 1 | , where | is any of the exact ground states and/or eigenstates in Eqs. (10) and (8) and p is the momentum of the excited state. The scattering matrix of these basis can be solved as easily as a one-body problem, despite the fact that Eqs. (10)
Exact charge neutral excitations with the flat metric condition being imposed for three different w 0 /w 1 at the twist angle θ = 1.05 • . Here we change w 0 while keeping w 1 = 110meV fixed. Other parameters are given in Appendix A. These excitations are exact at the charge neutrality point (ν = 0) for generic states and are exact at finite integer fillings if the flat metric condition is satisfied. The exact charge neutral excitations at different fillings without imposing flat metric condition are given in Figs. 7 and8 in Appendix D 4. Note the softening of further Goldstone modes from finite to zero w 0 , reflecting the symmetry enhancement of the first chiral limit. The continuum above the Goldstone modes is fundamentally made of independent particle-hole excitations
where R η mn (k) (Eq. ( 19)) and R η mn (k) are the charge ±1 excitation matrices. A valley-spin flavor in | ν (Eq. ( 8)) is either fully occupied or fully empty, thus {η 1 , s 1 } belongs to the valley-spin flavors which are fully occupied, while {η 2 , s 2 } belongs to the valley-spin flavors which are not occupied. Equation (26) shows that the neutral excitation scattering matrix is a sum of the two single-particle energies
) plus an interaction term. By translation invariance, the scattering preserves the total momentum p. The spectrum of the charge neutral excitations at each p is a diagonalization problem of a matrix of the dimension 4N M × 4N M , where the left and right indices are (k + q, m, m ′ ) and (k, m 2 , m 2 ), respectively.
The excitation spectrum with the flat metric condition Eq. ( 11) being imposed, i.e.,with the R η [Eq. (19)] being replaced by the simplified Eq. (20), is shown in Fig. 2. As explained in Appendix III A, the simplified charge ±1 matrices R and R do not depend on the filling ν. Thus the obtained charge neutral excitation dispersion also do not depend on ν. Figure 2 is exact for ν = 0 even when the flat metric condition is not satisfied since Eq. ( 20) is exact for ν = 0. The exact charge neutral excitations at different fillings without imposing the flat metric condition Eq. (11) in the (first) chiral-flat U(4) × U(4) limit. Only ν 0 states are tabulated since the symmetry and Goldstone modes of ν > 0 states are same as the ν < 0 states since they are related by the many-body charge-conjugation operator (Appendix A 2 c) [109]. Only states with ν + ν -are tabulated since | ν+,ν- ν and | ν-,ν+ ν have the equivalent little groups upon interchanging of the two U(4)s, and thus have the same number of Goldstone modes.
Number of GMs Ground states
| , the spectrum branch annihilating (creating) electrons in empty (occupied) states does not exist.
Solving Eq. ( 26) provides us with the expression for the neutral excitations at momentum p on top of the TBG ground states, including the Goldstone mode, whose dispersion relation can be obtained in terms of the quantum geometry factors of the TBG. In general, the scattering matrix is not guaranteed to be positive semidefinite, and negative energy would imply instability of the ground states. However, in a large (physical) range of parameters (Appendix A) of TBG at the twist angle θ = 1.05 • , we find that, as shown in Figs. 2,7, and 8, the energies of charge neutral excitations of the exact ground states | ν in Eq. ( 8) in the nonchiral-flat limit and | ν + ,ν - ν in Eq. ( 10) in the chiral-flat limit are non-negative, implying these are indeed stable ground states. As shown in Figs. 5 and7 and discussed in Appendixes C 4 and D 4, strong (first) chiral symmetry breaking may lead to an instability to a metallic phase.
In Tables II andI, we have tabulated the little group (defined as the remaining symmetry subgroup of the state) and the number of Goldstone modes for each ground state in Eqs. ( 8) and (10). As examples, here we only derive the little groups and number of Goldstone modes for | 1,1 -2 (10) and | -2 (8). The little groups and Goldstone modes for other states can be obtained by the same method. First we consider the ground state | 1,1 -2 in the (first) chiral-flat U(4) × U(4) limit, which has vanishing total Chern number. Recall that the U(4
In each
of the e Y = ±1 sectors, only one U(4) spin-valley flavor is occupied. Hence the little group of the state | 1,1 -2 in each e Y sector is U(1) × U(3), where the U(1) is the phase rotation in the occupied flavor and the U(3) is the unitary rotations within the three empty flavors. Thus the total little group of the state | 1,1 -2 is U(1) × U(3) × U(1) × U(3), which has the rank (number of independent generators) 20. Since the Hamiltonian has a symmetry group U(4) × U(4) which has rank 32, we find the number of broken symmetry generators to be 32 -20 = 12. On the other hand, since all the Goldstone modes we derived are quadratic [similar to the SU(2) ferromagnets, see Eq. ( 30)], it is known that [117] the number of Goldstone modes is equal to 1/2 of the number of broken generators, namely, 12/2 = 6. This is because a quadratic Goldstone mode is always a complex boson, which is equivalent to two real boson degrees of freedom corresponding to 2 broken generators.
Next, we consider the ground state | -2 in the nonchiralflat U(4) limit. Since the U( 4 3), where the U( 1) is within the occupied flavor and the U(3) is within the 3 empty flavors. Hence the number of broken generators is 16 -10 = 6, where 16 and 10 are the ranks of U(4) × U(4) and U(1) × U(3), and the number of (quadratic) Goldstone modes is 6/2 = 3.
In the above paragraph, we have shown that state
in the chiral-flat limit has three more Goldstone modes than | -2 in the nonchiral-flat limit, although their wave functions are identical. This is because, if we slightly go away from the (first) chiral-flat limit towards the nonchiral-flat limit, i.e.,take the parameter 0 < w 0 ≪ w 1 , some branches of the Goldstone modes will be gapped by a finite w 0 , as shown in Figs. 2,7, and 8.
The number of Goldstone modes can also be obtained by examining the scattering matrix in Eq. ( 26). Here we take | 1,1 -2 as an example. As discussed in Appendix IV C, in the first chiral limit, the state -2 using this property of scattering matrix.
Suppose the occupied flavors in
take the values in the two e Y -valley-spin flavors {+1, +, ↑}, {-1, +, ↑}, and {η 2 , s 2 } can only take values in the other three valley-spin flavors in each e Y sector. There are in total 6 nonvanishing channels. Since each channel has an zero mode given by the zero of S e Y 1 ,e Y 1 (k ′ , k; 0), there are 6 Goldstone modes, consistent with the group theory analysis in Table I.
, where the valley-spin flavor {η 1 , s 1 } with Chern band basis e Y 1 is fully occupied and the valley-spin flavor {η 2 , s 2 } and Chern band basis e Y 2 is fully empty. The PSDH scatters the basis to
where the scattering matrix S does not depend on η 1 , η 2 , s 1 , s 2 . The simple commutators between O -q,-G O q,G and fermion creation and annihilation operators in the chiral limit [Eq. ( 15)] lead to a simple scattering matrix. We here focus on the e Y 1 = e Y 2 , p = 0 channel. For generic ν = 0 states and the ν = ±1, ±2, ±3 states with flat metric condition [Eq. ( 11)],
we have
The general expression of S e Y 2 ;e Y 1 (k + q, k; p) for all channels without imposing the flat metric condition is given in Appendix D 2. We first show the presence of an exact zero eigenstate of Eq. ( 28) by remarking that the scattering matrix
This guarantees that the rank of the scattering matrix is not maximal, and that there is at least one exact zero energy eigenstate, with equal amplitude on every state in the Hilbert space:
. More details are given in Appendixes D 3 and D 3 a. The U(4)× U(4) multiplet of this state is also at zero energy. Moreover, the scattering matrix S e Y ;e Y (k + q, k; 0) is positive semidefinite. The details of this proof can be found in App. D 3 a.
Since the p = 0 state has zero energy, for small p, by continuity, there will be low-energy states in the neutral continuum. By performing a k • p perturbation in the p = 0 states in Eq. (D11), one can compute the dispersion of the low-lying states. Full details are given in Appendix D 3 b. In the chiral limit, and imposing the flat metric condition Eq. ( 11) we find, by using α a (k, q + G) = α a (-k, -q -G) for a = 0, 2 and as expected for the Goldstone of a FM, the linear term in p vanishes and
to second order in p. We find the Goldstone stiffness 3. The eigenvalues of the mass tensor of the Goldstone mode in the first chiral limit (31). Here w 1 is in units of meV.
Since C 3z symmetry is unbroken in the ground states in Eq. ( 10), an isotropic mass tensor m i j ∝ δ i j is expected. The eigenvalues of m i j with different values of w 1 are plotted in Fig. 3.
We now derive the charge ±2 excitations. We choose a basis for the charge +2 excitations as c † k+p,m 2 ,η 2 ,s 2 c † -k,m 1 ,η 1 ,s 1 | , where | is any of the exact ground states and or eigenstates in Eqs. ( 8) and (10) (for which (O q,G -A G N M δ q,0 )| = 0) and p is the momentum of the excited state. Hence {η 1 , s 1 }, {η 2 , s 2 } belong to the valley-spin flavors which are not occupied. The details of the commutators of the Hamiltonian and the basis are given in Appendix E. We find
where R η mn (k) are the charge +1 excitation matrices in Eq. (19). We see that the charge +2 excitation energy is a sum of the two single-particle energies plus an interaction energy. By the translational invariance, scattering preserves the momentum (p) of the excited state. The spectrum of the excitations at a given p is a diagonalization problem of a matrix of the dimension 4N M × 4N M . The scattering matrix T of the charge -2 excitations is derived in Appendix E 4. It has the same form as T here except that the charge +1 excitation The spectrum of charge ±2 excitations with the flat metric condition Eq. ( 11) imposed is shown in Fig. 4. By imposing the flat metric condition, we can replace the R η matrix [Eq. (19)] in Eq. ( 33) by the simplified Eq. (20). Since Eq. ( 20) does not depend on ν, the obtained charge +2 excitation dispersion also do not depend on ν. Figure 4 is exact for ν = 0 even when the flat metric condition is not satisfied since Eq. ( 20) is exact for ν = 0. Due to the many-body charge-conjugation symmetry at ν = 0 [109], the charge -2 excitations are degenerate with the charge +2 excitations. Exact charge ±2 excitations without imposing the flat metric condition Eq. ( 11) at different fillings are given in Figs. 9 and 10 in Appendix E 5.
The exact expression of the three-particle excitation spectrum [Eq. (33)] allows for the determination of the Cooper pair binding energy (if any). We notice the scattering matrix Eq. ( 33), T (η 2 ,η 1 ) (k + q, k; p), differs by a sign from the neutral charge energy Eq. (26). It is the sum of energies of two charge +1 excitations at momenta k + p, -k plus an interaction matrix, while Eq. ( 26) is the sum of charge +1 and -1 excitations minus an interaction matrix. This allows us to use the Richardson criterion [112][113][114][115] for the existence of Cooper pairing by examining the binding energy as follows: (34) where E (N ) is the energy of the lowest state at N particles. We now assume that the lowest state of the charge +2 excitation continuum obtained by diagonalizing the matrices Eq. ( 33) is the lowest energy state at two particles above the ground state, which is confirmed by numerical calculations for a range of parameters [111]. We note that Eq. ( 19) is the charge +1 excitation. The lowest energy of the noninteracting two-particle spectrum is 2 min(R), where min(R) is the smallest eigenvalue of R η (k) over k ∈ MBZ and valley flavors η = ±.
Hence we can write the binding energy as min(T ) -2 min(R), where min(T ) represents the minimal eigenvalues of T (η 2 ,η 1 ) (k + q, k; p) over momenta p ∈ MBZ and different valley flavors η 2 , η 1 . For later convenience, we denote the sum of the first two terms of Eq. ( 33) as Appendix
and the last term of Eq. ( 33) as
We therefore have T = T ′ + T ′′ in short notations. Here we have used the time-reversal symmetry:
as explained in Appendix E 3. We use Weyl's inequalities to find sufficient conditions for the presence and absence of superconductivity. In particular, for given p, η 2 , η 1 , the smallest eigenvalue of T (η 2 ,η 1 ) (k + q, k; p) is smaller than the smallest eigenvalue of T (η 2 ,η 1 )′ (k + q, k; p) plus the largest eigenvalue of T (η 2 ,η 1 )′′ (k + q, k; p). Hence we have min(T ) min(T ′ ) + max(T ′′ ) = 2 min(R) + max(T ′′ ).
Therefore a sufficient criterion for the presence of Cooper pairing binding energy is that T (η 2 ,η 1 )′′ (k + q, k; p) has all eigenvalues negative:
On the other hand, for given p, η 2 , η 1 , the smallest eigenvalue of T (η 2 ,η 1 ) (k + q, k; p) is larger than the smallest eigenvalue of T (η 2 ,η 1 )′ (k + q, k; p) plus the smallest eigenvalue of T (η 2 ,η 1 )′′ (k + q, k; p). Hence we have min(T ) min(T ′ ) + min(T ′′ ) = 2 min(R) + min(T ′′ ). Therefore, a sufficient criterion for the absence of Cooper pairing binding energy is that T (η 2 ,η 1 )′′ (k + q, k; p) is positive semidefinite:
From the charge +2 excitation spectra in Figs. 4,9, and 10 we can see that the spectrum of T consists of two parts: the two-particle continuum, which is given by the sums of two charge +1 excitations, and a set of charge +2 collective modes above the two-particle continuum. Thus it seems that T ′′ are always non-negative positive.
In Appendix E 3, we proved that, for the projected Coulomb Hamiltonian with the time-reversal symmetry, the matrix T (η,-η)′′ (k + q, k; p), which corresponds to excitations of two particles from different valley, is positive semidefinite. Thus there is no inter-valley pairing superconductivity of the PSDH H I at the integer fillings ν of the ground states in Eqs. ( 8) and (10). We expect this property to hold slightly away from integer fillings. Since TBG shows superconductivity at ν = 2 or slightly away from integer fillings, our results
show that either kinetic energy or phonons are responsible for pairing.
Here we briefly sketch the proof. We consider the expectation value of T (η,-η)′′ (k + q, k; p) on an arbitrary complex function φ m 2 ,m 1 (k):
As detailed in Appendix E 3, substituting the definition of the M matrix (Eq. ( 4) into Eq. ( 36), we can rewrite the expectation value as
where
For simplicity, we use a and b to represent the composite indices (Q, α). Then T ′′ η φ (p) can be written as ab W † ab V W ab , where now W ab is viewed as an N M × 1 vector and V an N M × N M matrix. Since V is positive semidefinite, for each pair of a, b, the summation over k 1 , k 2 , G 1 , G 2 is non-negative. Thus T ′′ is positive semidefinite since T ′′ η φ (p) 0 for arbitrary φ. In Appendix E 3, we also proved that, for η 2 = η 1 = η, T (η,η)′′ (k + q, k; p) is also positive semidefinite due to the symmetry PC 2z T , with P being the unitary single-body PH symmetry of TBG [43,108]. Therefore neither the intervalley pairing nor the intravalley Cooper pair has binding energy in the projected Coulomb Hamiltonian for any integer fillings ν in the chiral-flat limit, and for any even fillings ν = 0, ±2, ±4 in the nonchiral-flat limit.
In this paper, we have calculated the excitation spectra of a series positive semidefinite Hamiltonians (PSDHs) initially introduced by Kang and Vafek [71] that generically appear [109] in projected Coulomb Hamiltonians to bands with nonzero Berry phases and which exhibit ferromagnetic states as ground states, under weak assumptions [109,110]. These assumptions were also used by Kang and Vafek [71] to find the ν = 2 ground states in TBG. In this paper, we show that not only the ground states, but a large number of lowenergy excited states can be obtained in PSDHs. We obtain the general theory for the charge ±1, ±2 and neutral excitations energies and eigenstates and particularize it to the case of TBG insulating states. We find that charge +1 excitations are gapped, with the smallest gap at the Ŵ point. In both the (first) chiral-flat and nonchiral-flat limits, we find the Goldstone stiffness of the ferromagnetic state, as well as the Cooper pairing binding at integer fillings. In particular, we proved by the Richardson criterion [112][113][114][115] that Cooper pairing is not favored at integer fillings (even fillings when nonchiral) in the flat band limit. Since superconductivity has been observed in experiments with screened Coulomb potentials [7][8][9] (such as at ν = 2), we conjecture the origin of superconductivity in TBG is not Coulomb, but is contributed by other mechanisms, e.g.,the electron-phonon interaction [59,60,85], or due to kinetic terms. In particular, our theorem shows that the Luttinger-Kohn mechanism of creating attractive interactions out of repulsive Coulomb forces is ineffective for flat bands. A similar statement can be made for the superexchange interaction. A finite kinetic energy is hence required for these mechanisms.
In future work, the charge excitation energies of these Hamiltonians will be obtained in perturbation theory with the kinetic terms. A further question, of whether there are other further eigenstates of the PSDHs, remains unsolved.
The single-particle Hamiltonian, symmetries, and properties of the wave functions have been discussed at length in Refs. [43,107,108]. For completeness of notation, we give its expression here, for completeness, but we skip all details. The total single-particle Hamiltonian is
where c † k,Q,η,αs is the creation operator at momentum k (in the moiré BZ -MBZ) in valley η (±), sublattice α (1,2), spin s (↑↓), and moiré momentum lattice Q. The Hamiltonians in the two valleys are
where σ = (σ x , σ y ), σ * = (-σ x , σ y ) are Pauli matrices, T j = w 0 σ 0 + w 1 (cos 2π 3 ( j -1)σ x + sin 2π 3 ( j -1)σ y ), and q j = k θ C j-1 3z (0, 1) T ( j = 1, 2, 3) with k θ = 2|K| sin θ 2 being the distance of the Graphene K momenta from the top layer and bottom layer, θ the twist angle, w 0 the interlayer AA hopping, and w 1 the interlayer AB hopping. Q belongs to a hexagonal momentum space lattice, Q ∈ Q ± , where
The eigenstates of Eq. (A1) take the form
where c † k+G,Q,ηαs = c † k,Q-G,ηαs , for any moiré reciprocal wave vector G, and hence
such that c † k+G,nηs = c † knηs . This is the MBZ periodic gauge. In the numerical calculations, we take the parameters θ = 1.05 • , |K| = 1.703 Å-1 , v F = 5.944 eV Å, and w 1 = 110 meV. The projected kinetic Hamiltonian in the flat bands will be denoted by H 0 (without hat), which is given in Eq. (1).
The many-body Hamiltonian, symmetries, and properties of the wave functions, as well as the derivations, have been discussed at length in Refs. [109,110]. For completeness of notation, we give its expression here, for completeness, but we skip all details. The Hamiltonian before projection was derived to be (denoted by a hat) [109,110]:
where
A6) is the total electron density at momentum q + G relative to the charge neutral point. tot is the total area of the moiré lattice, G sums over the moiré reciprocal lattice, and q sums over momenta in MBZ zone.
For the analytic derivations in the current paper, we keep V (r) generic. For the numerical plots of the en-ergy dispersion and other properties, we use twisted bilayer graphene Coulomb interactions screened by the electrons in the two planar conducting gates [71,118]: 2 + n 2 with ξ = 10 nm being the distance between the two gates, U ξ = e 2 /(ǫξ ) = 26 meV (in Gauss units), ǫ ≈ 6 the dielectric constant of boron nitride. The derivation of this interaction was explained at length in Ref. [109,110]. It was also showed that the interaction has nonvanishing Fourier component only for intra-valley scattering to give
For the given values of the parameters, V (q) was plotted in Ref. [108] and is a slowly decreasing function of |q| in the BZ, reaching (around the magic angle) about half of its maximal value as |q| spans the whole MBZ around.
We define for our many-body Hamiltonian the form factors, also called the overlap matrix of a set of bands m, n as
(A8) In terms of which the projected density operator interaction and Hamiltonian to a set of bands denoted by m, n can be written as
For our working convenience, we then define the operator
This allows us to rewrite the projected interaction Hamiltonian H I into the form of Eqs. ( 2) and (3). While most of the projected Hamiltonian properties are valid for any number of projected bands that respect the symmetries of the system (including PH), in TBG at the first magic angle we usually are interested in the projection of the Hamiltonian onto the lowest two flat bands per spin per valley of TBG (8 bands in total). An important step in any calculations-especially numerical-is the gauge-fixing procedure. Different gauges for the wave functions, that make different symmetries of the form factors Eq. (A8) more explicit, can be chosen. This is explained at length in our manuscript Ref. [109], but for completeness we briefly mention them here. We consider only two active bands, the general gauge-fixing mechanism for projection in more than two bands is found in Ref. [109]. To fix the gauge of the Bloch wave functions in Eq. (A3),
, where u Q,α;mη (k) is the solution of the single-particle Hamiltonian h (η) QQ ′ (k), in each η = ±, s = ↑↓ sector, we label the higher energy band by m = + and the lower band by m = -. Due to the spin-SU(2) symmetry, we set the real space wave functions for s = ↑↓ to be identical and omit the index s for the single-particle states. For a symmetry operation g, the sewing matrices B g n ′ η ′ ,nη (k) = ψ gk,n ′ ,η ′ |g|ψ k,n,η , which relate states at momentum k with states at the transformed momentum gk can be consistently chosen to be
for the C 2z T , C 2z , P symmetries of TBG where ζ a , τ a (a = 0, x, y, z) are Pauli-matrices acting on the band and valley indices, respectively [109]. Here P is a unitary single-body PH symmetry that transforms k to -k [43,108]. We leave the other sewing matrices-for C 3z , C 2x -unfixed. With these sewing matrices, once we obtain, by diagonalizing the singleparticle Hamiltonian, the wave functions in the valley η = + for band m = +, then we first fix the C 2z T (in TBG, this is done at the detriment of the wave function being continuous), then we use of the PH to generate the m = -band, while finally, using C 2z symmetry to generate the wave functions in the valley η = -. In Ref. [109], we use the above gauge to determine the most generic form of the M matrix coefficient [Eq. (A8)]. We found that-away from the point w 1 = 0 in parameter space which we call "second chiral limit"-we can decompose the M matrix into four terms [108,109]
where α 0,1,2,3 are real functions which satisfy the following symmetry conditions:
In particular, the combination of Eqs. (A12) and (A13) implies that at q = 0, we have
Note that the same gauge fixing can be found with projection in a larger number of bands [109].
A further simplification occurs in the first chiral limit w 0 = 0 of the single-particle Hamiltonian due to the presence of an extra symmetry [72,108,109]. (A similar simplification takes place in the second chiral limit w 1 = 0, found in Refs. [108,109]). In this limit there is another (chiral) symmetry C of the one-body first-quantized Hamiltonian h QQ ′ (k), which in band space is defined by its sewing matrix:
A gauge choice in which
This allows us to find the wave function of theband at k from the + band at k, in the same valley. In Ref. [109], we prove that the M matrix in the interaction satisfies the chiral symmetry,
Thus with the chiral symmetry, M takes the form
In Ref. [108], we have shown that the two flat bands can be recombined as two Chern bands: The M matrix in the Chern band basis becomes
For later convenience, we have introduced the factors M e Y (k, q + G) and F e Y (k, q + G) to represent the diagonal element and off-diagonal element in the Chern band basis, respectively. In the first chiral limit, where α 1 = α 3 = 0, we have F e Y = 0.
In Ref. [109], we showed that the full projected Hamiltonian H 0 + H I has a many-body charge-conjugation symmetry, P c defined as the single-particle transformation C 2z T P followed by an interchange between electron annihilation operators c and creation operators c † :
:
(A22) The interaction has the charge-conjugation symmetry, and the many-body physics is PH symmetric in this limit.
Using the unitary PH symmetry P introduced in Ref. [43], we demonstrated in Ref. [109] that the projected TBG Hamiltonian has a U(4) symmetry if the kinetic energy is set to zero (flat band limit), for any number of projected bands. This generalizes the U(4) symmetry introduced in Ref. [72] for the two active bands. Using the continuous symmetry operator notation [e iγ ab S ab , H I ] = 0, ∀γ ab ∈ R;
we can generate a full set of U(4) generators given by [109]
(A24) where the Pauli matrices ζ a , τ a , s a with a = 0, x, y, z are identity and x, y, z Pauli matrices in band, valley, and spin space, respectively. For 2N 1 , N 1 ∈ Z projected bands, the generators would be identical, with the ζ representing the + = {1, . . . , N 1 } and -= {N 1 + 1, . . . , 2N 1 } bands.
The kinetic plus the projected interaction term exhibit the Cartan symmetry U(2) × U (2) subgroup of the U(4) symmetry group of the projected interaction, which can be most naturally chosen as the valley spin and charge: Cartan:
e. Enhanced U(4) × U(4) symmetries in the first chiral limit w 0 = 0
In Ref. [109], we demonstrated in detail the presence of two enhanced unitary U(4) × U(4) symmetry in two limits of the single-particle parameter space the first and second chiral limits w 0 = 0 < w 1 and w 1 = 0 < w 0 . For the first chiral limit w 0 = 0, this symmetry was presented in Ref. [72] for the case of two projected bands, but we find that it is maintained, in both chiral limits, for a projection in any number of bands. For the matrix elements in Eq. (A17), the interaction commutes with the following matrices, which form the U(4) × U(4) generators [109]:
The Cartan subalgebra of the chiral U(4
The U(4) implied by C 2z P [Eq. (A24)] is a subgroup of this U(4) × U( 4), but not one of the U(4) factors [109].
In Ref. [108], we found that there exists a further more convenient gauge choice for the wave functions in the chiral limit w 0 = 0, called the Chern basis (an extension to many bands of the Chern basis in Ref. [72] to many bands), in which we choose the single-particle representations of U(4) × U(4) generators as
x, y, z), (A26) which correspond to the first U(4) and second U(4), respectively. Adding the kinetic term in the chiral limit breaks the U(4) × U(4) symmetry of the projected interaction, to a U(4) subset ζ 0 τ a s b , (a, b = 0, x, y, z).
We note that the nonchiral-flat U(4) symmetry and the first chiral-flat U(4) × U(4) symmetry are first identified by Ref. [72]. A similar U(4) symmetry is proposed in Ref. [71], the difference and similarity between which and the symmetries reviewed here is studied in [109].
In Ref. [109], we showed that the eight single-particle basis of the nonchiral-flat U(4) symmetry generators given in Eq. (A24) can be decomposed into two four-dimensional fundamental irreps of the U(4) group, which have ζ y eigenvalues e Y = ±1, respectively, for each momentum k [Eq. (A18)]. The U(4) generators in the Chern band basis are s ab (e Y ) = e Y τ x s a , e Y τ y s a , τ 0 s a , τ z s a , respectively. We also showed [109] that the e Y = +1 irrep and the e Y = -1 irrep are the same-and not conjugate-irrep: the four-dimensional fundamental U(4) irrep represented by a one-box Young tableau labeled by [1] 4 . We presented a detailed review of the U(4) representations related to TBG in Ref. [109], but for the purpose of the current paper, the notation adopted for irreps is the standard Young tableau, conveniently denoted by [λ 1 , λ 2 , • • • ] N , where λ i is the number of boxes in row i (λ i λ i+1 ). The number of boxes in the i-th row is no smaller than that in the (i + 1)th row. The Hook rule then provides the dimensions of each of these irreps. In particular, [1 N ] N is an SU(N) singlet state.
In the first chiral limit w 0 = 0, d † k,e Y ,η,s defined in Eq. (A18) gives the single-particle basis irrep U(4) × U(4) of Eq. (A26). We proved in Ref. [109]
τ a s b , respectively. A similar discussion is provided for the second chiral limit w 1 = 0 in Ref. [109].
In the chiral-nonflat case, the generators U(4) symmetry is given by generators ζ 0 τ a s b and either the original band basis c † k,m,η,s for a fixed band index m(= ±) or the Chern basis d † k,e Y ,η,s with a fixed e Y (= ±1) form a fundamental of U(4).
In Ref. [110], we have examined in detail the exact ground states in the nonchiral-flat U(4) symmetric limit and in the chiral-flat U(4) × U(4) symmetric limit. For completeness we briefly review the results. In Ref. [71], Kang and Vafek first introduced a type of Hamiltonians which have, under some conditions, exact ground states. We have found [110] that any translationally invariant interaction Hamiltonian projected to some active bands can be written in a form of Ref. [71], and that, under some conditions, exact eigenstates and ground states can be found. The key idea of Kang and Vafek for obtaining exact ground states is to rewrite the interacting Hamiltonian into a non-negative form. A state with eigenvalue zero is then ensured to be the ground state.
In Ref. [109], we proved that the projected Coulomb Hamiltonian can be written as
which is non-negative, with
where ρ η k,q,m,n,s = c † k+q,m,η,s c k,n,η,s is the density operator in band basis. Note that O(r) and O(r ′ ) generically do not commute and hence the Hamiltonian is not solvable. The interaction can in general be rewritten as
where N M is the number of moiré unit cells and A G is some arbitrarily chosen G dependent coefficient satisfying A G = A * -G . Note that the first term in Eq. (A29) is nonnegative. In Ref. [110], we found an important condition which can show that eigenstates of H I are in fact ground states of H I . If the q = 0 component of the matrix element in Eq. (A8) M (η) m,n (k, G) is not dependent on k for all G's, i.e., flat metric condition:
then much more information about Eq. (A27) can be obtained. This condition is always true for G = 0, for which M (η) m,n (k, 0) = δ mn from wave-function normalization. In Ref. [107], we have showed that, around the first magic angle, M (η) m,n (k, G) ≈ 0 for, |G| > √ 3k θ for i = 1, 2. Hence, the condition Eq. (A30) is valid for all G with the exception of G for which |G| = √ 3k θ . Hence, the condition is largely valid, and the numerical analysis [111] confirms its validity for a large part of the MBZ. In the below, we will always specify when the condition Eq. (A30) is used. If the Eq. ( A30) is satisfied, one has O 0,G proportional to the total electron number N and the second term in Eq. ( A29) is simply a chemical potential term
where M = tot /N M is the area of moiré unit cell. For a fixed total number of electrons, N = k,m,η,s c † k,m,η,s c k,m,η,s = (ν + 4)N M is a constant, where ν is the filling fraction (number of doped electrons per moiré unit cell) relative to the charge neutrality point, thus the ground state at finite filling is solely determined by the first term which is non-negative.
To build the excitations around a ground state, we review the ground states found in Ref. [110] of the projected Hamiltonian (A27). We proved that in the Chern basis of Eq. (A20) and Eq. (A18), diagonal in the valley index η, spin index s and Chern band index e Y , the projected Hamiltonian (A27) has as eigenstates at integer filling ν the filled band wave functions [without assuming condition Eq. (A30)]:
where ν + -ν -= ν C is the total Chern number of the state, and ν + + ν -= ν + 4 is the total number of electrons per moiré unit cell in the projected bands, with 0 ν ± 4, k running over the entire MBZ. The occupied spin/valley indices {η j 1 , s j 1 } and {η j 2 , s j 2 } can be arbitrarily chosen. These eigenstates of Eq. (A27) are moreover eigenstates of the O q,G operator in Eq. (A28) [110]:
In Ref. [110], we found that the U(4) × U(4) irrep of this multiplet is labeled by
). For a fixed filling factor ν, from Eq. (A33), we found [110] that the states with different Chern number ν C are all degenerate.
At charge neutrality ν = 0, the U(4) × U(4) multiplet of eigenstate state | ν + ,ν - 0 with Chern number ν C = ν + -ν -= 205415-14 0, ±2, ±4 has exactly zero energy and hence are exact degenerate ground states. At nonzero fillings ν, we cannot guarantee that the ν = 0 eigenstates are ground states [without condition Eq. (A30)].
Assuming the flat metric condition Eq. (A30), we showed [110] that we can rewrite the interaction into the form of Eq. (A29), with the coefficient A G = ν √ V (G)ξ (G) in Eq. (A34). By Eq. (A34), we showed that (O q,G -A G N M δ q,0 ) annihilates | ν + ,ν - ν for any ν C = ν + -ν -and thus all the eigenstates | ν + ,ν - ν with any Chern number ν C = ν + -ν - are degenerate ground states at filling ν [110].
Without chiral symmetry, with U(4) symmetry Eq. (A24), O q,G is no longer diagonal in any band basis (such as the Chern basis). Nevertheless, O q,G is still diagonal in η and s and hence filling both m = ± bands, of any valley/spin is still an exact, Chern number 0 eigenstates [110]:
for even fillings ν = 0, ±2, ±4, where {η j , s j } are distinct valley-spin flavors which are fully occupied.
With M (η) m,n (k, q + G) in Eq. (A11), we have the same eigenvalue expression as in Eq. (A34), O q,G | ν = ν √ V (G)δ q,0 k,m,η,s α 0 (k, G)| ν . Along with any U(4) rotation, it is an eigenstate of H I , without using condition Eq. (A30). Moreover, for ν = 0, the state (A35) is always a ground state with or without condition Eq. (A30) [110]. Furthermore if the condition Eq. (A30) is satisfied, by choosing A G = ν √ V (G)ξ (G) we have showed in Ref. [110] that the states in Eq. (A35) are Chern number zero exact ground states. The multiplet of states forms a U(4) irrep [(2N M ) (ν+4)/2 ] 4 [110].
In order to compute the charge 0, ±1, ±2 excitations, a series of commutators are needed. We provide their expressions here.
In the nonchiral case of Eq. (A11), we have
and
where we have used the property [109]
n,m (kq, q + G). From these basic equations, we further find
and
where we define the new matrix element P = V M † M, the convolution of the Coulomb potential and the form factor matrices
These are the commutators needed to obtain the wave functions and energies of the excitations in the nonchiral limit.
In the first chiral limit, we can use the Chern band basis Eq. (A20), where the O q,G is diagonal in the Chern basis. Its form factors do not depend on the valley η and spin s:
In this limit, the commutators between O q,G and the Chern number e Y = ±1 band creation operators become simpler
and
leading to the commutators
where P = V M † M, the convolution of the Coulomb potential and the form factor matrices, takes the chiral limit form
where α 0 (k, q + G), α 2 (k, q + G) are the decomposition of the form factors in Eq. (A20). Notice in the Chern basis, P e Y (k, q + G) does not depend on e Y , so we just denote it as P(k, q + G). These are the commutators needed to obtain the wave functions and energies of the excitations in the first chiral limit.
To look for the charge one excitations (adding an electron into the system), we sum the commutators in Eqs. (B3) and (B10) over q, G and use the shifted Hamiltonian in Eq. (A29). For a generic exact eigenstate | at chemical potential µ satisfying (O q,G -A G N M δ q,0 )| = 0 for some coefficient A G , we find
where N is the electron number operator, and the matrix
is also an eigenstate of H I with eigenvalues obtained by diagonalizing the 2 × 2 matrix R η mn (k). In the case of TBG, this is a 2 × 2 matrix, hence the diagonalization can be done by hand, providing a band of excitations. We note that, the expression c † k,n,η,s | = 0 may vanish and give no charge excitation, for instance, if valley η and spin s is fully occupied. We now delve more into the energies and eigenstates c † k,m,η,s | . Due to the symmetry C 2z P [Eq. (A10)], the M matrix [Eq. (A11)] satisfies M (η) m,n (k, q + G) = mnM (-η) -m,-n (k, q + G). Correspondingly, the R matrix satisfies
Since R + (k) and R -(k) are related by a unitary transformation, they must have the same spectrum.
In the nonchiral limit, the eigenstates | ν we found in Ref. [110] (and re-written in Eq. (A35)) have only fully occupied or fully empty valley η and spin s flavors. For TBG, this means that both active bands m = ± are either full or empty for each valley η and spin s. In this case we can only obtain exact charge +1 excitation at even fillings, i.e.,ν = 0, ±2. We can consider two charge +1 states c † k,n,η,s | (n = ±) at a fixed k in a fully empty valley η and spin s. These two states then form a closed subspace with a 2 × 2 subspace Hamiltonian R η (k) defined by Eq. (C2). Diagonalizing the matrix R η (k) then gives the excitation eigenstates and excitation energies. It is worth noting that, due to Eq. (C3), the spectrum of R η (k) does not depend on η. Since the U(4) irrep of the ground state is | ν is [(2N M ) (ν+4)/2 ] 4 , the U(4) irrep of the charge 1 excited state is given by [(2N M ) (ν+4)/2 , 1] 4 . Furthermore, at ν = 0, the state | ν=0 in Eq. (A35) is the ground state of the interaction Hamiltonian H I and hence c † k,n,η,s | ν=0 is the charge excitation above the ground state. Note that this does not assume the "flat metric condition" (A30) and is hence fully generic.
If we further assume the flat band condition Eq. (A30), the eigenstates | ν become exact ground states, and the chemical potential is given by Eq. (A31). In this case, the 2 × 2 excitation sub-Hamiltonian R η (k) takes a simpler form:
which can be diagonalized to give the band excitation eigenstates and energies above the ground state at each momentum k.
The spectrum at every k is obtained from diagonalizing the matrix R η (k) = G,q V (G + q)M (η) † (k, q + G)M (η) (k, q + G), which depends (up to a convolution with the Coulomb potential), only on the projected band wave functions. This is clearly a sum (over q, G) of positive semidefinite matrices [remember that V (G + q) > 0]. Hence R η (k) is a positive semidefinite matrix, whose eigenvalues are non-negative (expected, since we proved that these are excitations above the ground state). We now find conditions that these excitations are gapped, i.e.,that the matrix R η (k) is positive definite at each k. We now show this by re-writing the R η (k) as
where now M (η) (k) is a matrix of (2N M • N G ) × 2 matrix (with 2 because we are projecting into the two active TBG bands), where N M is number of moiré unit cells, N G is the number of plane waves (MBZs) taken into consideration. In Ref. [107], we have showed that the number of plane waves needed is very small: the matrix elements fall off exponentially with |G| and any contribution above |G| = √ 3k θ is negligible. The matrix elements read
, and its rank is equal to Rank( √ V M η (k)) 2. We show the rank has to be 2 (or in general the number of occupied bands), by the simple argument
The second term is still a positive semidefinite matrix, while the first term is diagonal and has eigenvalues V (0)/2 tot . Hence by Weyl's theorem, the energies of the excited states are V (0)/2 tot . In general, our discussion shows that the states c † k,n,η,s | are not degenerate to the ground state | (note that we did not prove these are the unique ground states). However, we cannot exclude a gapless excitation.
In general, however, it seems that this method gives rise to finite gap charge 1 excitations. The largest gap happens in the atomic limit or a material, where u m k+q |u n k = δ mn , for which R mn = δ mn q,G V (q + G) = δ mn tot V (r = 0). Since we know that TBG is far away from an atomic limit-the bands being topological, we expect a reduction in this gap. However, we argue that this type of charge excitation is always gapped. We perform a different decomposition of the matrix R η mn :
Since the second term is still a semi positive definite matrix by construction, the eigenvalues of R η (k) will be greater or larger than the eigenvalues of the first term. In fact, due to Ref. [107], we know that the eigenvalues of the second term are negligible for |G| 2|b M1 |. Hence, using the math notation of A B for A -B positive semidefinite, we have
We call G mn (k, q) the generalized "quantum geometric tensor," whose trace is the generalized Fubini-Study metric. The property of the generalized quantum geometric tensor/Fubini study metric is that they become the conventional quantum geometric tensor/Fubini study metric for small transfer momentum q. The tensor quantifies the distance between two eigenstates in momentum space. The conventional quantum geometric tensor is defined as
in which m, n are energy band indices and i, j are spatial direction indices of n occupied orthonormal vectors u m (k) in a N dimensional Hilbert space, where k is some parameter. We can show that
Generically, we expect [77] that the overlap between two functions at k and k + q to fall off as q increases, leaving a finite term in R η mn (k), the electron gap, at every k.
For an eigenstate | ν + ,ν - ν defined by Eq. (A32) in the (first) chiral-flat U(4) × U(4) limit, one has the coefficients M(k, q + G) = ζ 0 τ 0 α 0 (k, q + G) + iζ y τ 0 α 2 (k, q + G). Without condition Eq. (A30), the eigenstate | ν + ,ν - ν (which is not necessarily the ground state) satisfies Eq. (A34), which is equivalent to choosing
for Eq. (C2). Using the relation (A12), we can simplify the matrix R η mn (k) defined in Eq. (C2) as In this plot we have used the parameters defined in Appendix A: v F = 5.944 eV Å, |K| = 1.703 Å-1 , w 1 = 110 meV, U ξ = 26 meV, and ξ = 10 nm. Note that the excitation gap is largely reduced from the flat-condition limit. + = -simply means that the chemical potential locates at the middle of conduction and valence bands.
In Figs. 5 and6, the charge ±1 excitations are plotted at different fillings and w 0 /w 1 's for two different screening lengths of the Coulomb interaction [Eq. (A7)], i.e.,ξ = 10 and 20 nm, respectively. The corresponding interaction strengths are U 10 nm = 26 meV and U 20 nm = 13 meV. We have used w 1 = 110 meV in all the calculations and w 0 /w 1 = 0, 0.4, 0.8 for ν = 0, -2 and w 0 /w 1 = 0 for ν = -1, -3.
For ν = 0 and ξ = 10, 20 nm, the charge ±1 gaps are at the Ŵ M momentum and are always larger than 10 meV for different w 0 /w 1 's. For ν = -2 and ξ = 20 nm, the charge ±1 gaps are always larger than 5 meV for different w 0 /w 1 's.
For ν = -2 and ξ = 10 nm, the charge ±1 gaps are finite for w 0 /w 1 = 0, 0.4 but become negative at w 0 /w 1 = 0.8 [Fig. 5(c)]. We find that the gaps close around w 0 /w 1 ≈ 0.75, implying that, with ξ = 10 nm and the other parameters we have used, the ground states in Eq. ( 8) become unstable for w 0 /w 1 ≈ 0.75.
The instability shown in Fig. 5(c) will lead to a metallic phase at ν = -2 in the nonchiral-flat limit with strong chiral symmetry breaking (w 0 /w 1 = 0.8). In the band structure picture, we can understand the charge +1 excitation as the conduction band and the charge -1 excitation as the reverted valence band. The negativities of both imply that the energy of conduction band overlaps with the energy of valence band. w 0 /w 1 =0.4 w 0 /w 1 =0.8 Since we have chosen to fully occupy the valence bands, the ground state energy is not minimized in this case: One can move one particle from the top of valence band to bottom of conduction band to lower the energy. The ground state energy will be minimized by redistributing electrons to occupy the bands below the chemical potential. Due to the overlap between conduction and valence bands, the redistributed band structure will have electron and hole pockets and hence is a metallic state. In the resulted metallic phase at ν = -2, there are two fully empty valley-spin sectors and two partially filled valley-spin sectors. The partially filled sectors contribute to the electron and hole pockets. This state is still invariant under a U(2) subgroup within the fully empty valley-spin sectors. Since U(4) and U(2) have 16 and 4 parameters, there will be 16-4 2 = 6 Goldstone modes.
While the charge ±1 excitations can be obtained by diagonalizing a 2 × 2 matrix, we can obtain the charge neutral, two-body excitations above the ground state. These can be obtained by diagonalizing a 2N M × 2N M matrix, or a one-body problem, despite the state having a thermodynamic number of particles. Due to the fact that we know the exact eigenstates (or ground states) of the system, building excitations of the Hamiltonian on top of these eigenstates (or ground states) becomes a problem of diagonalizing a basis formed only from the excitations. We now obtain the charge neutral excitations, show that they exhibit Goldstone modes with quadratic dispersion-as required by U(4) [or U(4) × U(4)] ferromagnetism, and obtain the stiffness of the Goldstone dispersion in the first chiral limit.
We choose a basis for the neutral excitations
where | is any of the exact ground states and/or eigenstates in Eqs. (A32) and (A35). The scattering matrix of this basis can be solved as easily as a one-body problem, despite the fact that Eqs. (A32) and (A35) hold a thermodynamic number of particles. We first have to compute the commutators:
By rewriting k 2 = k + p and k 1 = k, we can write the scattering equation as
The | are the states | ν in Eq. (A35), and hence η 1 , s 1 belong to the valley-spin flavor/s which are fully occupied, while η 2 , s 2 belong to the valley/spin flavor which are not occupied.
For a generic exact eigenstate | at chemical potential µ satisfying (O q,G -A G N M δ q,0 )| = 0 for some coefficient A G , we find that the scattering matrix reads
where R η mn (k), R η mn (k) are the ±1-excitation matrices in Eqs. (C2) and (C16). We see that the neutral energy is a sum of the two single-particle energies [first row of Eq. (D5)] plus an interaction energy [second row of Eq. (D5)].
The exact expression of the PH excitation spectrum allows for the determination of the Goldstone stiffness. The Goldstone of the U(4) and U(4) × U(4) ferromagnetic ground states is part of the spectrum of the neutral excitation Eq. (D4), and is the state at small momentum p = k 1 -k 2 . We will solve this in the simpler, chiral limit, but it can be obtained in the general, nonchiral limit Eq. (D4).
We now consider the charge neutral excited states reachable by creating one electron-hole pair with total momentum p on the chiral-flat limit eigenstate | ν + ,ν - ν . Assume the valley-spin flavor {η 1 , s 1 } has Chern band basis e Y 1 fully occupied and the valley-spin flavor {η 2 , s 2 } Chern band basis e Y 2 fully empty. We consider the Hilbert space of the following sets of states of momentum quantum number
The O q,G operators in the chiral limit have the simple, diagonal expression of Eq. (B7), which leads to the scattering equation.
The | are the states | ν + ,ν - ν in Eq. (A32), and hence e Y 1 , η 1 , and s 1 belong to the valley-spin flavor/s which are fully occupied, while e Y 2 , η 2 , and s 2 belong to the valley/spin flavor which are not occupied. The scattering matrix in the chiral limit does not depend on η 1 , η 2 :
where M e Y (k, q + G) is given in Eq. (A20) and R η 0 (k) is given in Eq. (C12). If condition Eq. (A30) is satisfied, the eigenstates | ν + ,ν - ν are the ground states, and hence the states Eq. (D6) are the neutral excitations on top of the ground states. Without the condition Eq. (A30), only ν = 0 states are guaranteed to be the ground states, although the others are still eigenstates. With condition Eq. (A30), we have
(D8) Solving Eq. (D6) provides us with the expression for the neutral excitations at momentum p on top of the TBG ground states.
We show that the Goldstone mode of the ferromagnetic ground states are included in the neutral excitations of Eq. (D6) and we obtain their dispersion relation, in terms of the quantum geometry factors of the TBG. We are able to analytically obtain the Goldstone mode if the condition Eq. (A30) holds. We first show the presence of an exact zero eigenstate of Eq. (D6).
We now show that Eq. (D6) has an exact zero energy eigenstate. In order to see this, we remark that the p = 0, e Y 1 = e Y 2 state Eq. (D8) has a scattering matrix
whose elements in every row sums to zero (irrespective of η 1,2 , s 1,2 ):
This guarantees that the rank of the scattering matrix is not maximal, and that there is at least one zero energy eigenstate, with equal amplitude on every
A U(4) × U(4) multiplet of this state is also at zero energy. Moreover, the scattering matrix S e Y ;e Y (k + q, k; 0) is positive semidefinite (as it should, since these eigenvalues are energies of excitations on top of the ground states). For the matrix S e Y ;e Y (k + q, k; 0), we prove that its negative, -S e Y ;e Y (k + q, k; 0), has only nonpositive eigenvalues, and hence S e Y ;e Y (k + q, k; 0) has only non-negative eigenvalues. For -S e Y ;e Y (k + q, k; 0), the diagonal elements -S e Y ;e Y (k, k; 0) are nonpositive, while the off-diagonal elements -S e Y ;e Y (k + q, k; 0) (q = 0) are non-negative. Hence by the Gershgorin circle theorem, all eigenvalues lie in at least one of the Gershgorin disks (which, due to the fact that the matrix is Hermitian, are intervals) centered at -S e Y ;e Y (k, k; 0) and with radiusq =0 S e Y ;e Y (k + q, k; 0). These intervals are
S e Y ;e Y (k + q, k; 0), 0 .
Since q =0 S e Y ;e Y (k + q, k; 0) 0, ∀k, all eigenvalues of -S e Y ;e Y (k + q, k; 0) are nonpositive, and hence all eigenvalues of S e Y ;e Y (k + q, k; 0) are non-negative.
Since the p = 0 state has zero energy, for small p, there will be low-energy states in the neutral continuum. By using the p = 0 states in Eq. (D11), one can compute their dispersion. First, we write the Hamiltonian matrix elements acting on the two particle states
Hence, for small p, the energy of the Goldstone mode is given by the expectation value
As expected for the Goldstone of a ferromagnet, the linear term in p vanishes; by using α a (k, q + G) = α a (-k, -q -G) for a = 0, 2 of Eq. (A13), we find can prove that the linear terms vanish exactly. To second order in p, we find the Goldtone stiffness
In Figs. 7 and8, the charge neutral excitations are plotted at different fillings and w 0 /w 1 's for two different screening lengths of the Coulomb interaction [Eq. (A7)], i.e.,ξ = 10 and 20 nm, respectively. The corresponding interaction strengths are U 10 nm = 26 meV and U 20 nm = 13 meV. We have used w 1 = 110 meV in all the calculations and w 0 /w 1 = 0, 0.4, and 0.8 for ν = 0, -2 and w 0 /w 1 = 0 for ν = -1, -3.
For the parameters we used, the ν = 0 states with w 0 /w 1 = 0, 0.4, 0.8 and ν = -1, -3 states with w 0 = 0 have nonnegative excitations and hence are stable. The ν = -2 states with w 0 /w 1 = 0, 0.4, and 0.8 are also stable for ξ = 20 nm. However, the ν = -2 states for ξ = 10 nm become unstable at w 0 /w 1 = 0.8 [Fig. 7(c)]. The instability can be understood from the instability of the charge ±1 excitations shown in Fig. 5(c). From Fig. 5(c), we can see that the charge +1 excitation is negative at some momenta between Ŵ M and K M and the charge -1 excitation is negative at Ŵ M . Combining a pair of these negative particle and hole one obtains a negative charge neutral excitation. As discussed in the end of Appendix C 4, this instability will lead to a metallic phase.
For the stable ground states, where the spectrum is nonnegative, the charge neutral spectrum consists of a particlehole continuum (the blue area in Figs. 7 and8) and a set of gapped collective modes. In Figs. 7 and8, we only plot the eigenvalues of the scattering matrix. In practice, the existence and degeneracy of an excitation mode also depend on the occupied U(4) flavors [and U(4) × U(4) flavors in the first chiral limit] of the ground state, as discussed in Appendix III A. In particular, the number of Goldstone modes for different ground states are given in Tables II andI. In general, the states | ν + ,ν - ν [Eq. (10)] in the (first) chiral-flat U(4) × U(4) limit with ν = 0, ±2 and ν + = ν -= ν+4 2 (such that it has vanishing Chern number) has more Goldstone modes than the state | ν [Eq. ( 8)] with the same filling ν. From w 0 /w 1 = 0.4 to w 0 /w 1 = 0 for the states at ν = 0, -2, we can clearly see in Figs. 7 and8 that a collective mode is soften and become Goldstone mode, consistent with the theoretical analysis.
The charge ±1 excitations can be obtained by diagonalizing a 2 × 2 matrix; the charge neutral above the ground state can be obtained by diagonalizing a 2N M × 2N M matrix, or a one-body problem, despite the state having a thermodynamic number of particles, due to the fact that we know the exact eigenstates (or ground states) of the system. We now show that the charge +2 excitations can also be obtained by diagonalizing a 2N M × 2N M matrix. The conditions for which Cooper pairing occurs are also obtained.
We choose a basis for the neutral excitations where |ψ is any of the exact ground states/eigenstates Eqs. (A32) and (A35). The scattering matrix of these basis can be solved as easily as a one-body problem. We first have to compute the commutators:
which in detail reads
By rewriting k 2 = k + p and k 1 = -k, we can write the scattering equation as
The | are the states | ν in Eq. (A35), and hence η 1 , s 1 , η 2 , and s 2 belong to the valley/spin flavor which are not occupied. For a generic exact eigenstate | at chemical potential µ satisfying (O q,G -A G N M δ q,0 )| = 0 for some coefficient A G , we find that the T (η 2 ,η 1 ) m 2 ,m;m 1 m ′ (k 1 , k 2 ; q) matrix reads
where R η mn (k), R η mn (k) are the +1 excitation matrices in Eqs. (C2) and (C16). We see that the charge +2 energy is a sum of the two single-particle energies [first row of Eq. (D5)] plus an interaction energy [second row of Eq. (D5)]. The exact expression of the charge +2 excitation spectrum allows for the determination of the Cooper pair binding energy (if any).
We now consider the charge +2 excited states reachable by creating two electron pair with total momentum p on the chiral-flat limit eigenstate | ν + ,ν - ν . Assume the valley-spin flavor {η 1,2 , s 1,2 } has Chern band basis e Y 1 , e Y 2 fully empty. We consider the Hilbert space of the following sets of states of momentum quantum number p
The O q,G operators in the chiral limit have the simple, diagonal expression of Eq. (B7), which leads to the scattering equation
The | are the states | ν + ,ν - ν in Eq. (A32), and hence e Y 1 , η 1 , s 1 , e Y 2 , η 2 , s 2 belong to the valley/spin flavor which are not occupied. The scattering matrix in the first chiral limit does not depend on η 1 , η 2 T e Y 2 ;e Y 1 (k + q, k; p) = δ q,0 (R 0 (k
where M e Y (k, q + G) is given in Eq. (A20) and R η 0 (k) is given in Eq. (C12). If condition Eq. (A30) is satisfied, the eigenstates | ν + ,ν - ν are the ground states, and hence the states Eq. (E7) are the neutral excitations on top of the ground states. Without Eq. (A30), only ν = 0 states are guaranteed to be the ground states, although the others are still eigenstates. With condition Eq. (A30), we have T e Y 2 ;e Y 1 (k + q, k; p) = δ q,0 G,q
where we have used the α 0,2 (k, q + G) = α 0,2 (-k, -q -G). Solving Eq. (E7) provides us with the expression for the charge +2 excitations at momentum p on top of the TBG ground states.
In Sec. V B, we have derived the sufficient conditions for the existence [Eq. (37)] and absence [Eq. (38)] of Cooper pairing binding energies. Now we prove that, in the projected Coulomb Hamiltonian with time-reversal symmetry T and the combined symmetry PC 2z T , where P is the unitary PH symmetry [43,108], T (η 2 ,η 1 )′′ m,m ′ ;m 2 ,m 1 (k + q, k; p) is guaranteed to be positive semidefinite. Thus the condition (38) for the absence of Cooper pairing binding energy is always satisfied. We write T (η 2 ,η 1 )′′ m,m ′ ;m 2 ,m 1 (k + q, k; p), which is defined as the third term in Eq. (E5), as
We first consider the case η 2 = -η 1 = η. Due to the time-reversal symmetry T |u n,η (-k) = |u n,-η (k) in the gauge Eq. (A10) [109], where |u n,η (k) is the 2N Q × 1 vector u Q,α;nη (k) in Eq. (A3), and the definition of the M matrix (Eq. (A8)), we have
Thus we can rewrite T (η,-η)′′ m,m ′ ;m 2 ,m 1 (k + q, k; p) as
We consider the expectation value of T (η,-η)′′ (k + q, k; p) on a complex function φ m 2 ,m 1 (k):
Using Eq. (A4), we have
We then define the matrix
such that T ′′ η φ (p) can be written as
where a and b are the indices of the matrix W (k + G). For each term with given a and b in Eq. (E16), we can view the summation over k 1 , k 2 , G 1 , and G 2 as W † ab V W ab , where now W ab (k + G) is viewed as a vector with the index k + G. Since V (k 2 + G 2k 1 -G 1 ) is a positive semidefinite matrix, W † ab V W ab must be non-negative, and hence T ′′ η φ (p) 0 for arbitrary φ. Therefore T (η,-η)′′ m,m ′ ;m 2 ,m 1 (k + q, k; p) is positive semidefinite at every p.
Then we prove that T (η 2 ,η 1 )′′ m,m ′ ;m 2 ,m 1 (k + q, k; p) with η 2 = η 1 = η is also positive semidefinite. Due to the symmetry PC 2z T |u n,η (k) = n|u -n,η (-k) in the gauge Eq. (A10) [109] and the definition of the M matrix [Eq. (A8)], we have
Thus we can rewrite T (η,η)′′ m,m ′ ;m 2 ,m 1 (k + q, k; p) as
Repeating the calculations starting from Eq. (E13), one can show that T (η,η)′′ m,m ′ ;m 2 ,m 1 (k + q, k; p) must be positive semidefinite. The only difference with the above proof is that the definition of the W matrix becomes
(E19)
Based on the above, the charge -2 excitations are trivial to obtain. We do not give the details, but just the expression for the scattering elements
The | are the states | ν in Eq. (A35), and hence η 1 , s 1 , η 2 , and s 2 belong to the valley/spin flavor which are not occupied. For a generic exact eigenstate | at chemical potential µ satisfying (O q,G -A G N M δ q,0 )| = 0 for some coefficient A G , we find that the T (η 2 ,η 1 ) m 2 ,m;m 1 m ′ (k 1 , k 2 ; q) matrix reads
where R η mn (k) are the -1 excitation matrices in Eq. (C16). We see that the charge -2 energy is a sum of the two single-particle energies [first row of Eq. (E21)] plus an interaction energy [second row of Eq. (E21)]. In particular, for the chiral limit, the scattering matrix elements are identical to those of charge +2, i.e.,Eq. (E9).
In Figs. 9 and 10, the charge ±2 excitations are plotted at different fillings and w 0 /w 1 's for two different screening lengths of the Coulomb interaction (Eq. (A7)), i.e.,ξ = 10 and 20 nm, respectively. The corresponding interaction strengths are U 10 nm = 26 meV and U 20 nm = 13 meV. We have used w 1 = 110 meV in all the calculations and w 0 /w 1 = 0, 0.4, and 0.8 for ν = 0, -2 and w 0 /w 1 = 0 for ν = -1, -3.
The charge +2 (-2) spectrum consists of a two-particle (two-hole) continuum [the blue (red) area in Figs. 9 and 10] and a set of gapped charge +2 (-2) collective modes. The energies in the two-particle (hole) continuum are just sums of two charge +1 (-1) excitation energies. Note that all the charge +2 (-2) collective modes appear above the two-particle (two-hole) continuum, implying the absence of Cooper pairing binding energy, as proved in Sec. V B. lowest states at fillings ν = -3, -1 to the first-order perturbation of the (first) chiral symmetry breaking. In this section, we derive the approximate charge 1 excitations above these perturbative lowest states.
Since these states have fully occupied Chern bands, it is convenient to work in the Chern band basis. We write the O q,G operator as a sum of a (first) chiral preserving term . In the (first) chiral-flat limit, the second term is simply
with E 0 being the ground state energy.
However, the second term also involve excited states when the (first) chiral symmetry is broken. To be specific, we have
The first term on the right-hand side will give the unperturbed ground state energy E 0 , whereas the other three terms yield excited states. Now we approximate it by projecting it into the Hilbert space with a single-particle excitation. Notice that d † k,e Y ,η,s generates a particle in an empty Chern band and the O 1 operator, by definition, generates particles in the empty Chern bands and holes in occupied Chern bands. Thus the terms d † k,e Y ,η,s O 1 -q,-G O 0 q,G and d † k,e Y ,η,s O 0 -q,-G O 1 q,G will at least generate two particles plus one hole. Hence they do not contribute to the projected equation. Now we consider the term d
According to the Wick's theorem, we have
Here : A := A represents the normal ordered form of the operator A with respect to | ν + ,ν - ν , where the operators that annihilate | ν + ,ν - ν is ordered on the right-hand side of the operators that do not. The second term ("• • • ") represent the normal ordered terms with one contraction. The first (second) term either annihilate | ν + ,ν - ν or generate three (two) particles plus two (one) holes. We hence will omit them. The third and fourth terms must vanish since we require the flavor {e Y , η, s} to be empty and hence ν + ,ν - ν |d † k,e Y ,η,s = 0. The last term is the Fock energy correction to the ground state energy. Therefore we conclude
where E 0 is the unperturbed ground state energy and E 0 is the Fock energy correct. The excitation energy is hence given by the spectrum of [H I -µN, d † k,e Y ,η,s ]. Following the calculation in Appendix B, we obtain
Acting the commutator of the interaction Hamiltonian and d † k,e Y ,η,s on the state | ν + ,ν - ν , we have
In the (first) chiral limit, where
, the right-hand side of the above equation only has one particle excitations. However, when the (first) chiral symmetry is broken,
will yield excitations with two 205415-30 particles plus one hole
We now approximate the right-hand side by projecting it into the one particle Hilbert space: We only keep the term satisfying k ′ = k + q, e ′′ Y = e ′ Y , η ′ = η, and s ′ = s:
where n e
Here we have made use of Eqs. (A21) and (A12). With the above approximation, we can write the excitation equation as
where Y ,e Y are limited for the empty Chern band (-1). In this case the R matrix is one-by-one with n -e Y ,η,s = 1. For the excitation in the other (fully empty) valley-spin sectors, the e Y and e ′ Y indices can be either 1 or -1. Hence the R matrix is a two-by-two matrix with n -e Y ,η,s = 0. The calculation for ν = -1 is similar. Due to Ref. [110], the perturbative ground state in the nonchiral-flat limit at ν
There is a fully occupied valley-spin sector and a half-filled valley-spin sector. In the half-filled valley-spin sector, R is one-by-one with n -e Y ,η,s = 1, and in the empty valley-spin sectors, R is two-by-two with n -e Y ,η,s = 0. The approximate charge +1 excitations are shown in Fig. 11.
In the nonchiral-flat limit, the states | ν + ,ν - ν at odd fillings ν = ±1, ±3 are no longer exact eigenstates of the Hamiltonian. However, in Ref. [110] we have shown that the states | 1,0 -3 [or its U(4) rotations] and | 2,1 -1 [or its U(4) rotations] are the lowest states at fillings ν = -3, -1 to the first-order perturbation of the (first) chiral symmetry breaking. Using the same method as in Appendix F 1, in this section, we derive the approximate charge -1 excitations above these perturbative lowest states.
For the same reason in Appendix F 1, the spectrum of excitation is given by the eigenvalues of [H I -µN, d k,e Y ,η,s ]. Following the calculation in Appendix F 1, we have where P η is defined in Eq. (F8). When the (first) chiral symmetry is broken, the term d k+q,e ′ Y ,η,s O q,G | ν + ,ν - ν yields excitations with two holes plus one particle
We now approximate the right-hand side by projecting it into the one hole Hilbert space: We only keep the term satisfying
where n e
Here we have made use of Eqs. (A21) and (A12). With the above approximation, we can write the excitation equation as
where Y ,e Y are limited for the occupied Chern band (+1). In this case the R matrix is one-by-one with n -e Y ,η,s = 0. And there is no hole excitation in the other (fully empty) valley-spin sectors. The calculation for ν = -1 is similar. Due to Ref. [110], the perturbative ground state in the nonchiral-flat limit at ν = -1 is | 2,1 -1 (or its U(4) rotations). There is a fully occupied valley-spin sector and a half-filled valley-spin sector. In the half-filled valley-spin sector R is one-by-one with n -e Y ,η,s = 0 and in the fully occupied valley-spin sector R is two-by-two with n -e Y ,η,s = 1. The approximate charge -1 excitations are shown in Fig. 11.
In the nonchiral-flat limit, the states | ν + ,ν - ν at odd fillings ν = ±1, ±3 are no longer exact eigenstates of the Hamiltonian. However, in Ref. [110], we have shown that the states | 1,0 -3 (or its U(4) rotations) and | 2,1 -1 (or its U(4) rotations] are the lowest states at fillings ν = -3, -1 to the first-order perturbation of the (first) chiral symmetry breaking. In this section, we derive the approximate charge neutral excitations above these perturbative lowest states.
Since these states have fully occupied Chern bands, it is convenient to work in the Chern band basis. Following the calculations in Appendix D 1, we have In the (first) chiral limit, where O -q,-G | ν + ,ν - ν = δ q,0 A -G N M | ν + ,ν - ν , the right-hand side of the above equation only involve excitations with one pair of a particle and a hole. However, when the (first) chiral symmetry is broken, O -q,-G | ν + ,ν - ν will yield additional particle-hole excitations: , where all the creation operators of occupied states and annihilation operators of empty states are on the right-hand side of the creation operators of empty states and annihilation operators of occupied states. The first (normal ordered) term is nonvanishing when acted on | ν + ,ν - ν if the two creation operators are of empty states and the two annihilation operators are of occupied states. However, we will omit this term because it yields two particle-hole pairs. Since we require d k,e Y 1 ,η 1 ,s 1 to be occupied such that the initial state is nonvanishing, there must be ν + ,ν - ν |d k,e Y 1 ,η 1 ,s 1 = 0 and hence the last two terms in the above equation vanish. where the scattering matrix is S η,s,η ′ ,s ′ ;η 2 ,s 2 ,η
The last term is nonzero if both the {η 2 , s 2 } and {η, s} valley-spin flavors are half filled. It couples all the half filled valley-spin flavors, which are independent in the (first) chiral-flat limit, to each other.
According to [110], in the nonchiral-flat limit the state 1,0 -3 (or its U(4) rotations) is still the perturbative ground state at ν = -3. Without loss of generality, we assume the occupied flavor is {+1, ↑, +1}. Then we can divide the neutral excitations into the following sectors.
(1) The half-half sector, where η 2 = η 1 = +1, s 2 = s 1 = ↑. The delta functions in the first term of Eq. (F31) require η = η ′ = +1, s = s ′ = ↑. The delta functions in the second term require η = η ′ and s = s ′ . Since d † η,s,e Y and d η ′ ,s ′ ,e ′ Y must belong to empty and occupied bands, there must also be η = η ′ = +1, s = s ′ = ↑ in the second term. Then it follows that e Y = e Y 2 = -1, e ′ Y = e Y 1 = +1. At given k, p, q, the S matrix is a one-by-one matrix.
(2) The empty-half sector, where {η 2 , s 2 } is an empty valley-spin sector and {η 1 , s 1 } is a half-filled valley-spin sector. The second term of Eq. (F31) vanish due to the delta function δ η 2 η 1 δ s 2 s 1 . The delta functions in the first term require η = η 2 , s = s 2 , η ′ = η 1 , and s ′ = s 1 . It follows that e ′ Y = e Y 1 = +1 and e Y , e Y 2 take values in ±1. At given k, p, and q, the S matrix is a two-by-two matrix.
According to Ref. [110], in the nonchiral-flat limit the state 2,1 -1 [or its U(4) rotations] is still the perturbative ground state at ν = -1. There is a fully occupied valley-spin sector and a half filled valley-spin sector. Without loss of generality, we assume the occupied flavors as {+1, ↑, +1}, {+1, ↑, -1}, and {+1, ↓, +1}. Then we can divide the neutral excitations into the following sectors.
(1) The half-half sector, where η 2 = η 1 = +1, s 2 = s 1 = ↓. The delta functions in the first term of Eq. (F31) require η = η ′ = +1 and s = s ′ = ↓. The delta functions in the second term require η = η ′ and s = s ′ . Since d † η,s,e Y and d η ′ ,s ′ ,e ′ Y must belong to empty and occupied bands, there must also be η = η ′ = +1, s = s ′ = ↓ (the half filled valley-spin sector) in the second term. Then it follows that e Y = e Y 2 = -1, e ′ Y = e Y 1 = +1. At given k, p, q, the S matrix is a one-by-one matrix. (2) The empty-half sector, where {η 2 , s 2 } is an empty valley-spin sector and {η 1 , s 1 } is a half-filled valley-spin sector. The second term of Eq. (F31) vanish due to the delta function δ η 2 η 1 δ s 2 s 1 . The delta functions in the first term require η = η 2 , s = s 2 , η ′ = η 1 , and s ′ = s 1 . It follows that e ′ Y = e Y 1 = +1 and e Y , e Y 2 take values in ±1. At given k, p, and q, the S matrix is a two-by-two matrix.
(3) The half-occupied sector, where {η 2 , s 2 } is the half filled valley-spin sector {+1, ↓} and {η 1 , s 1 } is the fully occupied valley-spin sector {+1, ↑}. The second term of Eq. (F31) vanish due to the delta function δ η 2 η 1 δ s 2 s 1 . The delta functions in the first term require η = η 2 , s = s 2 , η ′ = η 1 , and s ′ = s 1 . It follows that e ′ Y , e Y 1 take values in ±1 and e Y = e Y 2 = -1. At given k, p, and q, the S matrix is a two-by-two matrix.