{'fraction_non_alphanumeric': 0.049433633459192564, 'fraction_numerical': 0.01750411462871527, 'mean_word_length': 3.9271131463461173, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Haldane Statistics for Fractional Chern Insulators with an Arbitrary Chern number

In this paper we provide analytical counting rules for the ground states and the quasiholes of fractional Chern insulators with an arbitrary Chern number. We first construct pseudopotential Hamiltonians for fractional Chern insulators. We achieve this by mapping the lattice problem to the lowest Landau level of a multicomponent continuum quantum Hall system with specially engineered boundary conditions. We then analyze the thin-torus limit of the pseudopotential Hamiltonians, and extract counting rules (generalized Pauli principles, or Haldane statistics) for the degeneracy of its zero modes in each Bloch momentum sector.


I. INTRODUCTION
As the canonical example of topological order, the fractional quantum Hall (FQH) effect was originally discovered in two-dimensional electron gas subject to a strong perpendicular magnetic field. 1,2 Recently, several groups demonstrated numerically that these stronglycorrelated phases also exist in a topological flat band characterized by a non-zero Chern number C, even in the absence of a magnetic field. [3][4][5] This discovery of the so-called fractional Chern insulators (FCI) generated enormous interest. 6,7 Subsequent numerical studies [8][9][10][11][12][13][14][15][16][17] quickly confirmed the presence of more intricate single-component FQH states in lattice models [18][19][20][21] , such as the Read-Rezayi series [10][11][12]22,23 and the compositefermion states 14,15,24 . Powerful techniques from the study of FQH, including density algebra, 10,25-29 entanglement spectrum, 5,11,30,31 parton construction, [32][33][34] and the Hamiltonian theory of composite fermions, 35,36 were introduced to understand the topological ground state of FCI and the nature of its excitations. [37][38][39][40][41][42][43][44] Possible experimental realizations have also been proposed. 45,46 Most of the above progress dealt with a topological band with Chern number C = 1, which is essentially the same 18 as the continuum FQH in a periodic potential [47][48][49] . The strongly-correlated physics in a C > 1 Chern band [50][51][52][53][54] turned out much richer than the conventional FQH, due to the interplay between topological order and lattice structure. 33,[55][56][57][58] Barkeshli and Qi 55 mapped a C > 1 Chern band to a C-component lowest Landau level (LLL) using hybrid Wannier states 59 , and suggested the possibility to realize multicomponent FQH states in a single Chern band. Numerical studies 53,56,57,60 indeed found clear signature of such states, including the color SU(C) version of the Halperin 61 and the non-Abelian spin-singlet states 62 (NASS), but also identified qualitative deviations from these states, 56,57 which implies a more complex structure than proposed in Ref. 55. In a previous paper, 63 we proposed to understand these new features as the consequences of a special set of boundary conditions associated with the LLL mapping. In the simplest case, this alternative bound-ary condition can be understood as a color-dependent magnetic flux insertion. We demonstrated that the multicomponent LLL in a new Bloch basis can be seen as a single manifold with constant Berry curvature and Chern number C. Using pseudopotential Hamiltonians, we constructed model states for FCI with an arbitrary Chern number, and found high overlaps with the exact ground states. Crucially, our model states correctly capture the anomalous features in the particle entanglement spectrum of the C > 1 FCI that make our states distinct from the conventional multicomponent FQH states.
In this paper we provide details of the mapping between a Chern band and a multicomponent LLL, and demonstrate the distinctive features of our pseudopotential Hamiltonian due to the new boundary conditions. We construct, in a C-component LLL, a momentumspace Bloch basis and a hybrid Wannier basis that mimic the lattice counterparts. Both bases entangle the real space and the internal color space. Using the explicit one-body wave functions for the bases, we derive the representation of the projected density operators in both bases. We define model states as the exact zero modes of the pseudopotential Hamiltonian built from the projected density operators. As we demonstrated in our previous paper 63 , the Bloch basis is useful for numerical studies as it preserves the full lattice symmetry. The hybrid Wannier basis, on the other hand, facilitates the analysis of the pseudopotential Hamiltonian.
We give a detailed analysis of the simplest bosonic pseudopotential Hamiltonian for the Halperin colorentangled states. We show that the pseudopotential Hamiltonian reduces to almost classical electrostatics in the hybrid Wannier basis, when we take the so-called thin-torus limit [64][65][66][67][68][69][70][71][72][73][74][75][76][77] and carry out truncations motivated by previous numerical results. 56,57 This enables us to write down the form of its zero modes in this limit. However, in contrast to most well-known FQH states such as Laughlin and Read-Rezayi, a purely classical thin-torus description is not possible. We pinpoint the key difference from the conventional multicomponent FQH due to a subtle twist in the hybrid Wannier states, and detail the procedure to compute the total Bloch momentum of each zero mode. The resulting algorithm correctly pre-arXiv:1310.6354v1 [cond-mat.str-el] 23 Oct 2013 dicts the degeneracy of the FCI quasiholes in each lattice momentum sector, without resorting to numerical diagonalization, and can be seen as the extension of the generalized Pauli principle 78,79 to the color-entangled states.

II. ONE-BODY STATES IN A MULTICOMPONENT LOWEST LANDAU LEVEL
In this section we construct one-body bases in a multicomponent LLL that mimic the Bloch and the Wannier bases in a Chern band with an arbitrary Chern number C. 86 We consider a C-component (generalized spin) electron moving on a torus with a perpendicular uniform magnetic field. The major difference between our approach and the usual treatment of the multicomponent LLL problem is the adoption of a new set of boundary conditions. This alternative choice entangles together the C components and enables us to construct a single manifold of Bloch states with Chern number C. In contrast to the usual picture of multicomponent LLL as C separate manifolds (one for each of the C components) each with unity Chern number, our bases provide a natural foundation for the mapping to a single Chern band with an arbitrary Chern number C. The central result of this Section is Eq. (27), the expansion of the electron density operator in the Bloch basis.

A. Translations Operators
We consider electrons with C internal (color) degrees of freedom For simplicity, we work on a rectangular torus spanned by L x = L xx and L y = L yŷ , where L x and L y are the two fundamental cycles of the torus, andx andŷ are orthonormal. The torus is pierced by a magnetic field in the −ê z direction, B = ∇ × A = Bê z with B < 0. We denote by e < 0 the charge of the electron. The magnetic length is l B = /(eB). We define the total number of fluxes N φ penetrating the torus by Here we do not assume N φ to be an integer as in the original treatment of the Landau level on a toroidal geometry 80 . As we will see soon, the alternative set of boundary conditions we pick only requires This integer is equal to the dimension of the one-body Hilbert space in the lowest Landau level. We define the magnetic translation operator where is the guiding center momentum. The translation T (a) commutes with the one-body Landau Hamiltonian H = (−i ∇ − eA) 2 /(2m) but not with the translation T (b) at a different displacement, As argued in the introduction, we need to make contact between the multicomponent Landau level states and the Bloch states in a Chern band. For the latter, we consider a single Bloch band with Chern number C in a tightbinding model on a lattice with N x × N y unit cells. The band has a total of N x N y one-body states, one at each lattice momentum in the N x × N y Brillouin zone (BZ). To make contact with this lattice system, we first look in the Landau level for a pair of commuting translation operators that also resolve an N x × N y BZ. To this end, we tune the magnetic field to match the number of onebody states, and we consider the magnetic translations over a fictitious N x × N y unit cell structure of the continuous torus, namely, The operator T x (resp. T y ) has N x (resp. N y ) different eigenvalues. As opposed to the C = 1 case, however, for generic C they do not commute due to the N φ /(N x N y ) = 1/C flux over each fictitious plaquette, To compensate for this, we define the 'clock and shift' operators Q and P over the internal (color) Hilbert space by P |σ = |σ + 1 (mod C) , Q|σ = e i2πσ/C |σ . (10) Both operators are unitary, and they satisfy This leads to a pair of commuting composite operators We will refer to this pair as the 'color-entangled' magnetic translation operators. For the (color-neutral) Landau Hamiltonian, both operators are good symmetries, and they resolve an N x × N y Brillouin zone once we specify the boundary conditions. Notice that in general [T (L x ), T y ] = 0, [T (L y ), T x ] = 0. This means that we have to abandon the usual boundaries 80 T (L α ) = 1, α = x, y. Instead, we adopt the color-entangled generalization T Nα α = 1, namely, This alternative set of boundary conditions make it possible to construct two sets of basis states in the one-body Hilbert space with desirable properties spelled below.

B. Bloch and Wannier Bases
We define the Bloch states |k as the simultaneous eigenstates of T x and T y within the LLL, with k = (k x , k y ) ∈ Z 2 . The N x N y states within the first Brillouin zone have distinct eigenvalues under T α , and they constitute the Bloch basis in the N x N y = CN φ -dimensional Hilbert space of the C-component LLL. We now look for the explicit wave function x, y, σ|k for these basis states. We specialize to the Landau gauge A = Bxŷ. Consider the states |X, k y with X, k y ∈ Z defined by the real-and internal-space wave function 87 Here X, k y are state labels taking integer values, while x, y are real space coordinates taking continuous values, and σ ∈ Z C is a discrete coordinate in the internal color space. It is not hard to see that |X, k y belongs to the lowest Landau level, as the above wave function can be recast in the form f (x + iy, σ) e −x 2 /(2l 2 B ) . Moreover, we find that |X, k y is periodic in X, but with a twist in k y : These relations are reminiscent of the flow of hybrid Wannier states in a Chern insulator 37 . Moreover, as we prove in Appendix A, the color-entangled magnetic translations [Eq. (12)] have a representation on |X, k y similar to the representation of the lattice translations on the hybrid Wannier states, namely, We thus refer to these states as the hybrid Wannier states in the C-component LLL. It is easy to see the states with X ∈ [0 .. N x ) and k y ∈ [0 .. N y ) are linearly independent. We emphasize that unless N x is divisible by C, these states are not color eigenstates, in contrast to the states studied in Ref. 55. We want to define the Bloch states in the LLL as a Fourier sum of the hybrid Wannier states, 88 From Eqs. (17) and (18), we find that the simultaneous eigenvalue equation in (14) indeed holds. These states are periodic in k x , but only quasi-periodic in k y , The latter non-periodicity signals the topological obstruction to a periodic smooth gauge due to the non-zero Chern number of a Landau level. 89

C. Projected Density Operator
The density operator projected to the lowest Landau level plays a central role in the FQH physics, as it is used to define the inter-particle interaction. As we now show, this operator takes a particularly nice form in our Bloch basis.
By definition, the density operator of color σ at r = (x, y) projected to the LLL is given by where φ k (r, σ) = r, σ|k is the wave function of the Bloch state |k defined in Eq. (19), and k 1 , k 2 are each summed over a full BZ. 90 Since ρ(r, σ) must have torus periodicity, we can express it as a Fourier sum, Here, the wave vector q lives on the reciprocal lattice The projected density operator in momentum space for a single color component σ is thus given by where dr is over the torus [0, L x ) × [0, L y ). We define the full projected density operator ρ q by This operator is the building block of a color-neutral interacting Hamiltonian. In Appendix B, we finish the integral in Eq. (25) with the help of the sum over color σ, and prove the main result of this section, It should be noted that when N x is divisible by C, the integral in Eq. (25) can be finished for each σ individually, without the color sum. The above formula can be recast [using Eqs. (2) and (7)] as Note that the dependence on C enters only through the exponent shared by all ρ q and all terms in k .

D. Geometric Phase Structure
The above result suggests that the torus formed by the Bloch states |k is endowed with a rich geometric structure. As usual, the Berry connection between the BZ points k and k + q is defined as (the phase of) the inner product between the periodic part of the Bloch states |k and |k + q . This amounts to the matrix element of the operator e −iq·r between the two states, wherer is the position operator. Notice that this exponentiated position operator, when projected to the lowest Landau level, is nothing but the full density operator ρ q in Eq. (26). Therefore, we can interpret Eq. (27) as the parallel transport in the momentum space implemented by the projected density ρ q .
Define the primitive vectors on the reciprocal lattice g x = (2π/L x , 0) and g y = (0, 2π/L y ), and the shorthand notations ρ α = ρ q=gα and Phase[z] = z/|z| for z ∈ C. At momentum transfer q = g α , the (unitary) exponentiated Berry connection resolves the band geometry, while the norm is the quantum distance between k and k + g α . Notice that the quantum distance does not depend on k. 91 The gauge-invariant Berry phases can be extracted from parallel transport around closed loops of |k states over the BZ torus.
Given that we are interested in the Abelian Berry connection, each contractible loop can be decomposed into a product of loops around single plaquettes. Such plaquette Wilson loops take a particularly nice form for the Bloch states we constructed. Around the plaquette at k, is independent from k. Further, we can define the Berry curvature over a single plaquette 63 where takes the imaginary part in the principal branch log z ∈ (−π, π]. We find that the BZ torus for the multicomponent Landau level has constant Berry curvature and its Chern number is equal to the number of components In addition to the contractible loops, there are two independent non-contractible Wilson loops around the two fundamental cycles of the torus, related to charge polarization. We define The geometric phases over the BZ torus are fully specified by the following quantities For example, W x (1) can be obtained from W x (0) times the product of W (k) around each of the N x plaquettes between k y = 0 and k y = 1 in the first BZ. We can easily add a twist to the color-entangled boundary conditions in Eq. (13), The twist angles γ = (γ x , γ y ) ∈ R 2 implement colorindependent magnetic flux insertions. We incorporate this change by keeping (k x , k y ) ∈ Z 2 , but applying to every equation so far.

E. Twisted Torus
The above results can be directly generalized to a twisted torus. Instead of the rectangular torus spanned by L x = L xx and L y = L yŷ , we consider a torus with twist angle θ, spanned by The number of fluxes N φ is now defined by The reciprocal lattice primitive vectors g α are now defined by and we have the wave vector k = k x g x + k y g y , (k x , k y ) ∈ Z 2 . Once we change the wave functions of the hybrid Wannier states in Landau gauge A = Bxŷ to all of the earlier results still hold with no essential modifications. In particular, the proof in Appendix B can be adapted straight-forwardly (albeit with even more tedious algebra), and in Eq. (27) the density operator requires no formal change except for q = q x g x + q y g y . For the rest of the paper, we return to the rectangular torus.
The results can be similarly generalized to the twisted torus by simple substitutions.

III. PSEUDOPOTENTIAL HAMILTONIAN
With the one-body Bloch and hybrid Wannier bases at hand, we move to the many-body interacting problem. Our ultimate purpose is to build pseudopotential Hamiltonians for FCI with arbitrary Chern number C. As demonstrated in the last section, the multicomponent LLL resembles the Chern band once we impose appropriate boundary conditions that join together the C components. This link enables us to take advantage of the well-developed pseudopotential formalism in the LLL. We construct pseudopotential Hamiltonians (in the same way as those of single-component LLL 81,82 ) in the LLL from the projected density operator ρ q , and obtain its zero modes through numerical diagonalization. Following the usual practice in the FQH literature, 92 we define these zero modes at the FCI model wave functions.
Then, through the mapping between the Bloch states in the LLL and on the lattice, we transcribe these LLL wave functions to the lattice. The resulting trial wave functions can be directly compared with the FCI ground states obtained numerically for lattice Hamiltonians. As demonstrated in our earlier paper 63 , this approach yields model Hamiltonians adiabatically connected to the microscopic lattice Hamiltonian, and leads to trial wave functions with the correct total momentum on lattice and very high overlaps with the actual FCI ground states. Our trial wave functions also reproduce the anomalous particle entanglement spectrum as observed in Ref. 57.
The question remains, however, how to predict the total lattice momentum for the trial wave functions (including quasiholes) without numerical diagonalization, similar to the methods developed for the FQH 78,79 . For C = 1, this problem was solved by two of us 10 by combining the generalized Pauli principle 78,79 for singlecomponent FQH states (including quasiholes) with lattice folding. For C > 1, we now have the LLL-to-lattice mapping. What we still lack is a multicomponent version of the generalized Pauli principle. Refs. 83,84 studied this problem for the usual boundary conditions. Due to our modifications to the boundary conditions, their results do not directly apply here.
Fortunately, we can also extract the generalized Pauli principle from the Hamiltonian in the thin-torus limit [64][65][66] . In this limit, the hybrid Wannier orbitals in the LLL become isolated from each other. Specifically, we find from Eq. (16) that the ratio between the width of the hybrid Wannier orbital and the spacing between them scales as Therefore, when the aspect ratio L x /L y satisfies the hybrid Wannier orbitals are so separated that the projected density operator becomes approximately diagonal in the hybrid Wannier basis. As a result, the pseudopotential Hamiltonian built from projected density operators also becomes approximately diagonal in the hybrid Wannier basis. (This is not true for certain non-unitary states 85 .) By analyzing the classical electrostatics of the leading terms in the Hamiltonian, we can obtain the quantum numbers of the Hamiltonian zero modes. (For FQH with the usual boundary conditions, this was done in Refs. 65,67,69.) After the Bloch mapping between FCI and FQH, this will give us a counting rule for the degeneracy of the FCI quasiholes in each lattice momentum sector.
In the rest of this Section, we expand the new pseudopotential Hamiltonian proposed earlier 63 in the Wannier basis, and perform the necessary resummation to make it amenable to proper truncation in the thin-torus limit. The actual truncation and the analysis of the zero modes of the truncated Hamiltonian is left for the next Section.

A. Projected Density in the Hybrid Wannier Basis
We obtain the projected density operator in the hybrid Wannier basis by plugging the Fourier transform Eq. (19) into Eq. (27), Notice that the phase factor depends on the summation variables X, k y only through the linear combination XN y + k y C, which is proportional to the center position of the hybrid Wannier orbital |X, k y [Eq. (16)], wherex is the position operator in the x direction. This motivates us to index these orbitals by their center position. In the following, we introduce an alternative labeling |j, s for the Wannier states. The j index gives the center position of the Wannier state while the s index plays a role similar (but not identical) to the color index σ. As we will see in the next Section, the projected interaction decays exponentially when the difference in the j indices between two particles increases. As seen from Eq. (16), the hybrid Wannier state |X, k y depends on (X, k y ) only through and in the exponential and the Kronecker-δ in Eq. (16), respectively. For integers X, k y , the linear combination XN y + k y C must be an integer multiple of the greatest common divisor (GCD) Therefore, we introduce two integer labels For future convenience, we also define integers We emphasize that j and s are not independent. This can be seen by examining the solutions to the first equation in Eq. (49). For a given j, if (X, k y ) is a solution, then all the solutions can be parametrized as (50)], corresponding to n = 1, . . . , C in X +nC/ C. For a given j, we denote this set of C allowed values of s by A useful property is which follows from the fact that j → j+d can be achieved by k y → k y + 1 without touching X. Plugging Eq. (49) k y x mod L x into Eq. (16), we find that indeed we can relabel the hybrid Wannier states modulo the identification An example is given in Fig. 1. It is not hard to see that this mapping is bijective, although we cannot easily write down an explicit formula for the solution (X, k y ) to Eq. (49) at a given (j, s). We denote the k y solution formally as Then, the representation of the color-entangled magnetic translations T α in the |j, s basis can be constructed indirectly from Eq. (18), T y |j, s = e −i2πky(j,s)/Ny |j, s .
The wave functions for |j, s can be obtained from Eq. (16), In parallel to Eq. (17), |j, s is periodic in s but quasiperiodic in j, As we will see soon, this twist in s when shifting j is the main issue that sets the current problem apart from the usual multicomponent FQH. 83 We now want to expand the projected density operator in the relabeled hybrid Wannier basis. On the one hand, notice that due to the quasi-periodicity of |X, k y [Eq. (17)], the double sum of (X, (44) can be shifted to any set of N x N y points in the Z 2 plane, as long as the corresponding hybrid states are independent from each other. On the other hand, notice that label a set of N x N y hybrid Wannier states that are independent from each other for any given j 0 ∈ Z. Therefore, we can rewrite the double sum in Eq. (44) as a sum over the above set. Since increasing k y by q y while keeping X constant sends (j, s) to (j + q y d, s), we have where the primed sum is over j ∈ [j 0 .. j 0 + M ) for an arbitrary j 0 ∈ Z, with M = N x N y / C [Eq. (50)]. The appearance of j +q y d, s| requires special attention: when we shift j + q y d back to [j 0 .. j 0 + M ) using Eq. (58), the s index must be changed accordingly, by N x (mod C). This boundary effect dictates that ρ q is not diagonal in s unless N x is divisible by C, which discourages a seemingly plausible interpretation of s as an effective spin index in general.

B. Interacting Hamiltonian
We consider only interactions between a pair of colorneutral densities ρ q . The relevance of such interactions to the Chern insulators was justified numerically in our previous paper 63 . Such interactions can be specified in terms of the Haldane pseudopotentials. Higherbody pseudopotentials 82 can be implemented in the same spirit. We consider only the first two pseudopotentials (V 0 , V 1 ) being non-negative, with all V m>1 = 0. The interaction strength at momentum transfer q then reads and the Hamiltonian is given by Here q is summed over the infinite reciprocal lattice. As shown in our previous paper 63 , the color-entangled generalizations of the bosonic/fermionic Halperin singlet states and the corresponding quasihole states are defined as the exact zero modes of the above Hamiltonian (using V 1 = 0 for the bosonic case). These states are distinct from the usual Halperin states due to the color-entangled boundary conditions inherent in ρ q . Through numerical diagonalization, we can obtain these zero modes, and then transcribe them to the lattice system of an arbitrary Chern insulator using the one-body mapping between the LLL Bloch states and the lattice Bloch states. We now attempt to achieve an analytic understanding of this Hamiltonian, by exploiting its assumed adiabatic connectivity 63 to the thin-torus limit.
We first plug Eq. (60) into Eq. (62) and write H in the relabeled hybrid Wannier basis, where M and d are defined in Eq. (50), and for q = (2πq x /L x , 2πq y /L y ), we have We want to massage the above expansion of H to a form amenable to justified truncation in the thin-torus limit. The main obstacle is obviously the oscillatory factor e −i2πqx(j1−j2+qyd)/M in the coefficient. This can be removed in exchange for a Gaussian factor by performing a Poisson resummation over q x , which does not appear in the index of the creation/annihilation operators. After some straightforward but tedious algebra in Appendix C, we find Terms in the expansion of the pseudopotential Hamiltonian. Here we illustrate the example of d = 3. Each dot represents a term (qy, ∆) in Eq. (69). The weight of each term decays exponentially in its distance from the origin. The dashed circle marks the empirical threshold for truncation (qyd) 2 + (∆ − qyd) 2 = d 2 . The solid black dots inside are the density-density terms in Eq. (72), while the four solid gray dots contain the pair hopping and the density-density terms in Eqs. (75) and (85).
where ∆ is summed over an interval of length M centered around q y d, and we have defined the shorthand

IV. THIN-TORUS ANALYSIS
In Eq. (65), the Hamiltonian has been organized into groups of density-density or pair hopping terms. The strengths of the terms decay exponentially in the limit This is exactly the thin-torus limit in Eq. (43). In the following, we perform a proper truncation of the Hamiltonian in this limit and analyze the degeneracy and quantum numbers of its zero modes. The thin-torus analysis is a well-known, powerful technique to tackle the strongly-correlated physics in singlecomponent FQH effect. 65,67,69,77 In the thin-torus limit, the pair hopping terms die off quickly, and the Hamiltonian becomes classical, dominated by density-density terms and thus solvable. (This is not true for certain non-unitary states 85 .) One can obtain the correct degeneracy of the ground states and extract their total momenta simply by minimizing the classical electrostatic energy and completely ignoring the pair hoppings. By assumed adiabatic connectivity, 63 the results must also apply to the isotropic limit. The thin-torus analysis thus provides an intuitive picture for the 'root partitions' and the underlying generalized Pauli principle of Refs. 78,79. Our multicomponent pseudopotential Hamiltonian with color-entangled boundaries (65) turns out to be considerably more complicated, due to the essential role played by the pair hopping terms. As we will see soon, the largest pair hopping terms have strengths comparable to the subleading density-density terms. Keeping only the leading density-density terms results in too many zero modes compared with the numerical studies 56,57,63 . The correct ground-state degeneracy is recovered only after we put back the largest pair hoppings, which turn out to be of similar strength as some of the density-density terms. This indicates that the thin-torus limit of our multicomponent pseudopotential Hamiltonian cannot be described by classical electrostatics alone. The useful result of this Section is a set of rules [Sec. IV D] that correctly predict the degeneracy and total lattice momenta of FCI ground states (including quasiholes). This is illustrated by explicit examples in Secs. IV E and IV F.

A. Truncation of Bosonic Hamiltonian
Numerical studies in Refs. 56,57 found gapped FCI phases of bosons at filling ν = 1/(C + 1) with (C + 1)fold degenerate ground states, stabilized by on-site interactions projected to a topological flat band with Chern number C. In the following we specialize to the simplest case of bosons and try to understand the ground states of the pseudopotential Hamiltonian at filling ν = 1/(C + 1) and with quasiholes. Setting V 0 = N φ L y /L x > 0 and where the primed sum of j is over for an arbitrary j 0 ∈ Z [Eq. (59)], while ∆ is summed over the interval of length M given in Eq. (66).
In the β 1 limit, we can safely truncate the sum over n to a single term at n = 0, if we assume that M/d = N φ 1. Further, only the terms with q y ∼ 0 and ∆ − q y d ∼ 0 have a significant contribution, since the coefficients decay exponentially with respect to the (squared) Euclidean distance from q y d = ∆ − q y d = 0, as illustrated in Fig. 2. The 4-boson ψ † ψ † ψψ operator can be either density-density interaction or pair hopping. We find that the terms with q y = 0 are all density-density interactions, while the strongest pair hopping terms may appear at |q y | = 1, ∆ = q y d, with Euclidean distance In light of the previous studies 65,67,69,77 , we first examine the effect of the terms with R 2 (q y , ∆) < d 2 . They can be collected into where the number operator n j is defined by Recall from Eq. (51) that S j is the set of all allowed values of s for ψ j,s at a given j, and this set contains C different values. Also, recall from Eq. (50) that d C = C. By solving the simple electrostatics, we find that the zero modes of H <d 2 with highest density appear at filling ν = 1/C. This leads to much more than (C + 1) zero modes at filling 1/(C + 1), inconsistent with the findings from numerical diagonalization of actual FCI Hamiltonians. 56,57 . This is a clear signal that we should include more terms in the truncated Hamiltonian.
In the following we analyze the effect of the next strongest terms in Eq. (69), with Euclidean distance R 2 (q y , ∆) = d 2 . They are located at (|∆ − q y d|, |q y d|) = (0, d) and (d, 0), represented by the four solid gray dots in Fig. 2. In the next section we provide detailed analysis of the simplest case with d = 1. The results for general d will be presented afterwards.

B. Effect of Pair Hopping Terms: d = 1
In this subsection we specialize to the simplest case d = 1, illustrated in Fig. (3). In this case N y is divisible by C [Eqs. (48) and (50) collect together the terms nicely, recall from Eq. (52) that S j = S j+1 at d = 1, and note that we can take advantage of the freedom in Eq. (70) to shift the range of the primed sum over j. We then find The four terms in the above brackets are labeled by (q y , ∆), and explicitly they are given by Further, notice that we can combine the above four terms into a single product, b † j,s,s b j,s,s , where the pair annihilation operator is given by b j,s,s = ψ j,s ψ j+1,s + ψ j+1,s ψ j,s .
This combination is the key to the enumeration of zero modes as we detail below. Together with the densitydensity terms in Eq. (72), the bosonic pseudopotential Hamiltonian takes the truncated form The residual terms are exponentially small for β 1. When C = 1, the s index can take only a single value s = 0, reducing b j,s,s to 2ψ j,0 ψ j+1,0 . This includes the case of Chern number C = 1. The truncated Hamiltonian becomes very simple: Its zeros modes have no more than one boson in two consecutive orbitals. We thus recover the familiar result 65,78 for the bosonic Laughlin state at half filling.
We now come back to the case with generic C. We look for the constraints on the zero modes of the above truncated Hamiltonian in Eq. (79). Due to the two-body nature of the interaction, we only need to consider a pair of bosons at a time, with j indices being j 1 , j 2 . In Eq. (79), each term in the summation is positive-semidefinite by itself. This means that to find the zero modes of Eq. (79), we only need to identify the null space of each term individually, and then take their intersection. From the density-density terms, we find that in a zero mode we must have This amounts to a minimal distance between adjacent bosons along the j axis, with no discrimination of the s indices. The pair hopping terms b † b in Eq. (79) kick in only when the equality sign is taken in Eq. (81), as is evident from Eq. (78). Specifically, b † b enforces in a zero mode the antisymmetrization of the s indices between bosons with |j 1 − j 2 | = 1, We emphasize that the ψ † 's are bosonic operators. It is easy to verify that the above antisymmetrized form is indeed annihilated by b † b, whereas the symmetrized form acquires a positive energy 2e −β . To find the zero modes for a system of N bosons, we need to perform the above procedure on each pair of bosons. This is explained in more details in Sec. IV D, and illustrated by an example in Sec. IV F.
One last subtlety comes from the quasi-periodicity of the j index [Eq. (58)]. The orbitals at j +M are identified with those at j, but there is a possible mismatch between the s indices, For the density terms, this does not make much trouble since n j = n j+M after the summation of the s index over S j [Eq. (73)]; we just need to enforce the minimal distance condition [Eq. (81)] across the quasi-periodic boundary j = 0 mod M . For the pair hopping terms, however, we have to be more careful about the s index mismatch. We have to first shift their j indices (by integer multiples of M ) such that |j 1 − j 2 | = 1 before we can apply the antisymmetrization in Eq. (82). More explicitly, if |j 1 − j 2 + M | = 1 for example, then the correct antisymmetrization can be either of the following two equivalents, but not Eq. (82) anymore. This is the only reason why we were not able to consistently implement 63 the exclusion principle for conventional multicomponent FQH model states 83,84 for the color-entangled system.

C. Effect of Pair Hopping Terms: General d
The analysis for general d is not much different from d = 1. Here we just state the essential results. The pair hopping and density-density terms at (q y d) 2 + (∆ − q y d) 2 = d 2 can be merged together, where the two-body annihilation operator is given by b j,s,s = ψ j,s ψ j+d,s + ψ j+d,s ψ j,s .
Combined with the density-density terms in Eq. (72), the leading terms in the bosonic pseudopotential Hamiltonian in the limit of β 1 are The zero modes of the truncated Hamiltonian satisfy the following pairwise constraints. First, for a pair of bosons with j indices being j 1 and j 2 , we must have Here the difference in j is understood with the quasiperiodic identification j ∼ j + M . When the equality in Eq. (88) holds, the two bosons are further subject to an antisymmetrization in the s indices. For the simplest case |j 1 − j 2 | = d, we need Eq. (82), whereas for |j 1 − j 2 + M | = d, we need either of the two equivalents in Eq. (84). When C = 1, as s can take only one value, this antisymmetrization consistently reduces to an electrostatic repulsion at distance |j 1 − j 2 | = d (and also |j 1 − j 2 + M | = d).

D. Counting Rule for Degeneracy and Momenta
Following the above constraints, we can enumerate all the zero modes of the truncated Hamiltonian for a given system size and a given number of particles, in the form where A antisymmetrizes the s indices as follows. As noted earlier, for any pair of particles a and b in a zero mode, we must have |j a − j b | ≥ d, and when the equality holds, we need to carry out antisymmetrization over the s indices (s a , s b ). Obviously, if we have j 1 − j 2 = d and j 2 − j 3 = d, then we need to antisymmetrize over (s 1 , s 2 , s 3 ). More generally, if we have a cluster of m consecutive particles satisfying j a − j a+1 = d, we need a full antisymmetrization over all the s indices of these m particles.
The last remaining step is to group these zero modes according by the total Bloch momentum and count the degeneracy per momentum sector. The resulting degeneracy is linked by the Bloch mapping 63 to the degeneracy of FCI ground states per lattice momentum sector. This largely follows the same procedure as detailed in Ref. 10. We represent by lowercase k α the Bloch momenta of individual particles in the α = x, y direction, and by uppercase the total Bloch momentum of the many-body system (the summation is over particles). We denote by T cm α the center-of-mass color-entangled magnetic translations, i.e. applying T α simultaneously on all the particles. Then, the total Bloch momentum K α can be read off from the eigenvalue of T cm α , The action of T cm α on the zero modes in Eq. (89) is spelled out in Eq. (56).
There are four points to make here. First, the zero modes in the form of Eq. (89) are automatically eigenstates of T cm y . Evidently each term in the antisymmetrization A individually is an eigenstate of T cm y . Moreover, the eigenvalues have to be the same for all those terms. This follows from the linearity of Eq. (49): to find the total k y of all particles, we only need to know the total j and s; the actual association of between j and s does not matter. Second, under the action of T cm x , the zero modes in Eq. (89) form closed orbits. This follows from the fact that T cm x commutes with the (truncated) pseudopotential Hamiltonian, and thus preserves its null space. More directly, one can easily verify that the constraints on the zero modes described in Secs. IV B and IV C are invariant under the action of T cm x (namely X → X + 1, or |{j, s} → |{j + N y / C, s + 1} ), and that the action of T cm x always brings one zero mode in the form of Eq. (89) to another zero mode in the same form. Third, each action of T cm x along the orbit is associated with a sign, since a term in the antisymmetrization A in Eq. (89) may be brought to a term with the opposite sign. 93 Fourth, all the zero modes in an orbit under T cm x share the same eigenvalue under T cm y . This is a direct consequence of [ T cm x , T cm y ] = 0.
For each zero mode in the form of Eq. (89), we can directly compute the total K y momentum by just looking at a single term in the antisymmetrization A. We can group together the zero modes by the value of K y = k y mod N y . Then, within each group, we successively apply T cm x on each zero mode and further break them into disjoint orbits. Consider an orbit consisting of n zero modes |0 , |1 , . . . , |n − 1 of the form in Eq. (89). They are linked together by T cm x |r = g r |r + 1 mod n , r ∈ [0 .. n), with g r = ±1 determined from the action of T cm x on the antisymmetrization in Eq. (89). The n eigenstates of T cm x are linear recombinations of these n states in the form of Fourier sums. Without actually writing down the linear recombinations, we can directly obtain the eigenvalues. By successively applying the above equation, we find with g = n r g r . This fixes the n eigenvalues of T cm x to be the n distinct n-th roots of g. If g = 1, the total K x momenta of the zero modes are whereas if g = −1, they are The numbers on the right hand side of the above equation are guaranteed to be integers: Since [ T cm x ] Nx is the identity operator per the color-entangled boundary condition [Eq. (13)], we must have N x /n ∈ Z, and also g Nx/n = 1.
Our end goal is an analytic algorithm to obtain the degeneracy of the zero modes in each Bloch momentum sector. This request is more modest than to find the actual expression of the zero modes in each sector, and the above procedure can be further simplified. For example, we do not need to actually write down the zero modes as in Eq. (89). We only need to keep track of the structure of clusters of consecutive particles with j a − j a+1 = d, as noted below Eq. (89), and the set of s indices in each cluster. An open-source reference implementation can be found at http://fractionalized.github.io. We have tested our algorithm extensively against the total Bloch momenta of the actual ground states obtained from numerical diagonalization for various system sizes, and found perfect agreement across all cases.
All of them satisfy either |j 1 − j 2 | = d (first two) or |j 1 − j 2 + M | = d (last one), and are thus subject to full antisymmetrization of the s indices (s 1 , s 2 ). Since there are only two allowed values of s, we can already see that there are only 3 zero modes in the form of Eq. (89). We now go through them one by one. First, consider (j 1 , j 2 ) = (0, 1). Using Eq. (82), we find that the only possible (s 1 , s 2 ) antisymmetrization is Here the double bracket |·, · distinguishes the manybody zero mode from the one-body basis state |j, s , and the subscript of the creation operator ψ † denotes (j, s).
Similarly, for (j 1 , j 2 ) = (1, 2), we find The case of (j 1 , j 2 ) = (0, 2) satisfies |j 1 − j 2 + M | = d rather than |j 1 − j 2 | = d. So we use Eq. (84) rather than Eq. (82), and find Notice that after we bring the j indices back to [0 .. M ) using Eq. (83), the s indices on the second line are not in an explicit antisymmetrized form. This manifests the core difference of our problem from the usual FQH: When the lattice size is incommensurate with the Chern number, we cannot consistently distinguish the C families of Wannier states, since the flow of Wannier centers are entangled on the quasi-perioidic boundary.
Using the lookup table in Eq. (98), we find that the total K y momenta of the three zero modes are all equal to 0 mod 2, consistent with Eq. (96). To compute the K x momentum, we need to find out the action of the center-of-mass translation T cm x on these states. For our example, Equation (56) reduces to T x |j, s = |j + 1, s + 1 .
We thus find the representation of T cm x on the zero modes: Notice that we can evaluate T cm x |0, 2 using either line in Eq. (102); the results are guaranteed to be the same by the consistency between Eqs. (56) and (58).
The three zero modes thus form a single orbit under the successive action of T cm x . They can be recombined to form eigenstates of T cm x . To find the total K x momenta of the recombined states, we can either follow the procedure detailed in the last subsection, or we can brute-force diagonalize T cm x . The representation matrix of T cm x over the three zero modes reads From its eigenvalues {1, e i2π/3 , e −i2π/3 }, we find the total K x momenta of the three recombined zero modes to be 0, 1, 2 mod 3. In summary, we reproduce the correct total Bloch momenta in Eq. (96).

F. An Example with Quasiholes
Next, we consider a slightly more complicated example with quasiholes. For a system of 3 bosons with C = 2, the densest zero modes of our pseudopotential Hamiltonian occur at filling ν = 1/(C + 1) = 1/3, i.e. 3 bosons in 9/2 fluxes. The fractional flux is possible thanks to the color-entangled boundary conditions in Eq. (13). We now add 1 2 flux to each color component and consider N x × N y = 5 × 2 and N φ = 9/2 + 1/2 = N x N y /C = 5. This leads to a set of 10-fold degenerate quasihole states at zero energy, with one mode in each momentum sector . This can be seen in the numerical diagonalization results in Fig. 5.
We now show how to obtain this counting using our algorithm. The basic procedure is the same as the previous example. We first compute C = GCD(C, N y ) = 2, d = C/ C = 1, and M = N x N y / C = 5. Equation (49) again reduces to Eq. (97), and we have two allowed values of s (0 and 1), denoted by ↓ and ↑. Then we can build the (X, k y ) ↔ (j, s) lookup table, Again, the last line in the table has a flipped s index due to the quasi-periodic boundary condition in j [Eq. (58)]. Compared with the previous example, the enumeration of the zero modes in the form of Eq. (89) has an extra complication. Let us first apply the |j a − j b | ≥ 1 rule between each pair of bosons. We find two groups of allowed (j 1 , j 2 , j 3 ) configurations, (0, 1, 2), (1, 2, 3), (2,3,4), (3,4,0), (4, 0, 1), We can similarly write down the other 8 zero modes. This gives the correct 10-fold degeneracy. Using the lookup table in Eq. (106), we can compute the K y lattice momentum for each zero mode and construct the representation matrix of the color-entangled center-of-mass translation operator T cm x in exactly the same manner as in the previous example. This reproduces the correct degeneracy in each momentum sector. We leave details of this last step for the interested readers.

V. CONCLUSION
In this paper we have studied the pseudopotential model Hamiltonian for FCI with an arbitrary Chern number. We establish a one-body mapping between a Chern band with Chern number C, and a C-component LLL with specially engineered boundary conditions. The new boundary conditions lead to an alternative set of pseudopotential Hamiltonians, and the corresponding zero modes define new model wave functions. By taking the thin-torus limit and keeping only the leading densitydensity and pair hopping terms, we are able to analytically solve the pseudopotential Hamiltonian and obtain its zero modes. By analyzing the representation of the center-of-mass translation operators, we derive an algorithm to directly compute the total Bloch momenta of the degenerate zero modes. As we showed in our last paper 63 , our pseudopotential Hamiltonian is adiabatically connected to the lattice FCI Hamiltonian, and the its zero modes serve as good trial wave functions for the FCI ground states. In particular, there is a 1-to-1 correspondence between the trial wave function and the FCI ground state in each momentum sector. Therefore, our counting algorithm can be used to obtain the total lattice momenta of the FCI ground states (including quasiholes) without diagonalizing the FCI Hamiltonian, for Abelian FCI states at filling ν = 1/(C + 1).
In this Appendix we prove Eq. (18), the representation of T x and T y in the hybrid Wannier basis |X, k y .
In Landau gauge A = Bxŷ, the magnetic translation operators T x and T y defined in Eq. (8) have the real-space representation Acting on a trial state |ψ , they transform the real-space wave function x, y|ψ by x, y|T x |ψ = e i2π N φ Nx y Ly x − L x /N x , y|ψ , x, y|T y |ψ = x, y − L y /N y |ψ . (A2) Plugging these into the Landau-gauge definition of |X, k y in Eq. (16) and using Eq. (2), we find x, y, σ|T x |X, k y = x, y, σ + 1|X + 1, k y , x, y, σ|T y |X, k y = e −i2π(ky/Ny+σ/C) x, y, σ|X, k y . Since the clock-and-shift operators Q, P defined in Eq. (10) are unitary, we have Putting Eqs. (A3) and (A4) together, we find the action of T x = T x P and T y = T y Q to be particularly simple, x, y, σ| T x |X, k y = x, y, σ|X + 1, k y , x, y, σ| T y |X, k y = e −i2πky/Ny x, y, σ|X, k y .

Appendix B: Projected Density in Bloch Basis
In this Appendix we prove Eq. (27), the expansion of the projected density operator in the Bloch basis, proof which, due to lack of space, was not included in Ref. 63.
We first derive a simpler form for the Bloch wave function φ k (r, σ) = r, σ|k . When we plug Eq. (16) into Eq. (19), we have a double sum over X, m. However, notice that (X, m) in the double sum can always be combined into X + mN x , thanks to XN y /C + mN φ = (X + mN x )N y /C and e i2πXkx/Nx = e i2π(X+mNx)kx/Nx (recall that N φ = N x N y /C). This enables us to merge the double sum into a single sum of X +mN x over Z. The Kronecker-δ enforcing σ = X + mN x mod C suggests we split X + mN x → nC + σ with n summed over Z. This leads to the final form of the Bloch wave function, This wave function indeed depends only on σ mod C (by a re-shift in the dummy variable n), and it has the quasi-periodicity in k y as in Eq. (20). We now plug this into ρ q,σ defined in Eq. (25).
Notice that the summations of k 1 and k 2 over BZ in the above equation are independent. To accommodate the Kronecker-δ in the above equation, we set the summation of k 1y over [0 .. N y ), and the summation of k 2y over [q y .. N y + q y ). Then, the Kronecker-δ above can be decomposed into two separate Kronecker-δ's, enforcing And we have It is easy to check that ρ q,σ is indeed invariant under a shift of the dummy variable k y → k y + N y . We now tackle the dx integral in the bracket. We can collect terms and complete the square in the exponential. After some trivial but tedious algebra, the integrand becomes e −q 2 l 2 B /4 e −i2πqx(ky+qy/2+σNy/C+nNy)/N φ exp − x l B − 2πl B L y k y + σ C N y + nN y + 1 2 Here we have used 2πl 2 B N φ = L x L y [Eq. (2)], and The projected density can thus be written as Ny ky |k 1x , k y k 2x , k y + q y | e −i2πqx(ky+qy/2)/N φ Z n e i2π(nC+σ)(k2x−k1x−qx)/Nx Notice that we can shift the integration interval to This moves the dependence on n from the integrand to the integration limits (and also the exponential prefactor e i2π(nC+σ)(k2x−k1x−qx)/Nx ). We want to sew together the integrals for all n so that we can finish the Gaussian integral, but the integration intervals for different n are overlapping and cannot be joined head to tail in general, unless N x is divisible by C. However, recall that (to have symmetries P, Q) we restrict the interacting Hamiltonian to be color-neutral, so we are interested only in ρ q = C σ ρ q,σ . The color sum saves us. Notice that the dependence on (σ, n) is all through the combination nC + σ. We can merge the two sums over σ and n into a single sum over integers, nC + σ → m: Ny ky |k 1x , k y k 2x , k y + q y | e −i2πqx(ky+qy/2)/N φ Z n e i2πn(k2x−k1x−qx)/Nx Nx L x ) is covered by the integral for N x times, and during the N x times, the exponential prefactor runs through all the N x values of e i2πt(k2x−k1x−qx) for t ∈ [0 .. N x ). In formula, we have Z n · · · dx · · · = Nx n e i2πn(k2x−k1x−qx)/Nx The "mod N x " does not lead to any problem, since |k is periodic in k x . Finally, we arrive at Eq. (27): BZ k e −i2πqx(ky+qy/2)/N φ |k k + q|.

Appendix C: Pseudopotential Hamiltonian Reorganized
In this Appendix we prove Eq. (65), the reorganized expression for the pseudopotential Hamiltonian in Eq. (63) suitable for truncation. ψ † j1,s1 ψ † j2,s2 ψ j2−qyd,s2 ψ j1+qyd,s1 , (C1) where G V (j 1 − j 2 , q y ) is defined by Through a Poisson resummation, we can easily prove the general formula Z qx e −πAq 2 Setting A = L y /(N φ L x ), ξ = (k + q y d)/M , and defining To handle G V (k, q y ) in Eq. (C2), we need to be able to insert powers of q 2 x into the q x sum. This can be achieved by taking partial derivative with respect to πLy N φ Lx = βd 2 on Eq. (C4). For the simple case of V q = 4πl 2 B [V 0 + V 1 · (1 − q 2 l 2 B )] as in Eq. (61), we find Plugging this back to Eq. (C1), we get ψ † j1,s1 ψ † j2,s2 ψ j2−qyd,s2 ψ j1+qyd,s1 , (C6) where C = L x /(N φ L y ) is an inconsequential overall factor. At last, recall from Eq. (59) that the range of summations over j 1 and j 2 each contain an arbitrary shift. We can keep the outer sum over j 1 general, while make a convenient choice for the inner sum over j 2 . We define ∆ = j 2 − j 1 and rewrite the above equation as We make this special choice for the ∆ sum to facilitate later truncations in the thin-torus limit β 1. This proves Eq. (65).