Wigner crystallization in topological flat bands

We study the Wigner crystallization on partially filled topological flat bands. We identify the Wigner crystals by analyzing the cartesian and angular Fourier transform of the pair correlation density of the many-body ground state obtained using exact diagonalization. The crystallization strength measured by the magnitude of the Fourier peaks, increases with decreasing particle density. The shape of the resulting Wigner crystals is determined by the boundary conditions of the chosen plaquette and to a large extent independent on the underlying lattice, including its topology, and follows the behavior of classical point particles.


I. INTRODUCTION
In recent years, the possibility of realization of the quantum Hall effect (both integer and fractional) without a net magnetic field was intensely studied on topologically nontrivial energy bands of two dimensional (2D) lattice systems [1]. The nontrivial topology of a band is described by a nonzero value of an integer topological invariant named Chern number [2]. When a band with Chern number C = 0 is fully filled, it exhibits Hall conductivity quantized to an integer multiple of e 2 /h, in analogy to a fully filled Landau level in integer quantum Hall effect (IQHE). Such a system is called a Chern insulator. It was proposed that topologically nontrivial bands can arise entirely without a magnetic field in presence of artificial gauge fields acting on cold atom systems [3,4]. This proposition was later achieved experimentally [5][6][7][8]. Another way to realize such bands experimentally is to combine spinorbit interaction with ferromagnetism [9].
Adiabatic continuity between the FCIs and FQHE states was shown for C = 1 bands [24].
For larger Chern numbers, it was found that an adiabatic connection exists between FCIs and multicomponent FQHE states with a special, color entangled, boundary condition [25].
Moreover, the FCIs can be related to the Hofstadter model-the tight-binding model of a lattice in presence of uniform background magnetic field, which can be regarded as a discretized version of the quantum Hall system [26]. There is no fundamental physical difference between a topological flat band and a subband of the Hofstadter model thus the lattice FQHE states in the Hofstadter model can be considered as fractional Chern insulators (see Ref. 27 and the discussion in Ref. 28). Such states were recently observed in bilayer graphene, which can be regarded as the first experimental demonstration of FCIs [29]. There is a number of propositions of experimental realization of FCIs without a magnetic field, including cold atom [30][31][32][33][34][35] and solid state systems [36][37][38].
The subject of the Wigner crystallization in TFBs remains largely untouched in previous works. Several authors investigated the charge ordering induced by short-range interaction at high filling factors [69][70][71][72][73][74][75]. Phase diagrams of various flat-band models were obtained, showing the competition between the FCI and charge-ordered ground state [71][72][73][74]. Moreover, it was found that the charge ordering can coexist with topological ordering [73,75].
However, contrary to the Landau levels in which the Wigner crystallization occurs at arbitrarily low fillings, the short-range nature of interaction considered in  limits this effect to a certain filling factor.
In this work, we demonstrate the Wigner crystallization of spinless particles populating TFBs, interacting via short-and long-range potentials. We follow the exact diagonalization (ED) approach from Ref. [42][43][44]59 and calculate the exact ground states of variety of finite size systems in torus geometry on kagome, honeycomb and checkerboard lattices. A periodic pattern, corresponding to the Wigner crystal, is found in the pair correlation density (PCD). We analyze it using the Cartesian and angular Fourier transform, finding that the strength of the Fourier peaks -corresponding to the strength of the Wigner crystallizationincreases with decreasing filling factor. While there are differences in the shapes of the WC unit cells related to the range of interaction, the results are to a large extent independent of the lattice type, in consistence with a picture of interacting classical point particles in a continuous space. Finally, we compare the results for trivial and nontrivial bands of the Haldane model, showing no significant differences between them.

II. MODEL AND METHODS
Three lattice models with nearly flat bands are considered: kagome [12], honeycomb (Haldane model) [1,13] and checkerboard [11,13], with parameters chosen such that the lowest band of all three models is topologically nontrivial and nearly flat. For each model we have |C| = 1, where C is the Chern number of the lowest band, thus the same set of FCI phases can in principle be realized at each of them. The general form of a single-particle Hamiltonian is where c † i (c i ) is the creation (annihilation) operator at site i, while t ij , φ ij are modeldependent parameters, explained in Appendix A 1. We consider the systems of dimensions L 1 × L 2 = aN 1 × aN 2 in a torus geometry, with N 1 and N 2 being the number of unit cells in the two directions and a a lattice constant. We fill them with N part particles and apply the density-density interaction of the formV = i,j V (r ij )n i n j , where r ij is the shortest distance between the two atoms i and j, with periodic boundary conditions included [62,63,76].
Note also that the other treatment of interactions in strongly correlated systems have been applied, i.e. the Ewald summation, where a sum over all periodic repetitions is taken into account [59]. It is obvious that both approaches give the same results for sufficiently short interaction range, and it was also shown that periodic images give neglecting contribution for a dipolar type of interaction [76]. Our first choice for V (r) is the screened Coulomb , where α is a parameter describing the range of interaction. In the limit α → ∞ the interaction contain only nearest-neighbor terms, while for α → 0 it converges to unscreened 1/r Coulomb interaction. We consider also the logarithmic interaction defined as V Log β (r) = β−ln(r) β for r ≤ exp(β) and V Log β (r) = 0 otherwise, where short range interaction corresponds to small β, while for β → ∞ it converges to V (r) = 1. Both kinds of interactions are normalized, with V (r) = 1 between nearest neighbors.
We determine the ground state using the exact diagonalization method. We consider a projection of the full Hamiltonian of the system to a subspace of the lowest band, similarly to the lowest Landau level projection in FQHE. That is, we first solve the single-particle problem, and then construct the many-particle configuration basis out of the single-particle wavefunctions belonging to the lowest band. Since the wavefunctions are labeled by the momentum k, and the interaction conserves the total momentum of a many-particle state, we divide the basis into corresponding subspaces and diagonalize the Hamiltonian in each of them separately. We apply the flat band approximation, i.e. we neglect the single-particle dispersion by artificially setting the single-particle energies to zero for all k, which is a common procedure in the research on fractional Chern insulators and topological flat bands. In such a way, the only relevant energy scale in the calculation is two-body interaction strength.
However, for the approximations to be meaningful, the interaction energy scale should be larger than the band dispersion and much smaller than the energy gap. The calculations has been performed using highly parallel ED software utilizing adaptive load-balanced on-the-fly matrix-vector multiplication or Hamiltonian storage in compressed sparse blocks format [77], depending on available system resources, paired with ARPACK eigensolver. The configuration basis of the largest system considered in this work a 7 × 10 plaquette with N part = 7 has size ∼ 1.2 × 10 9 (before division into 70 momentum subspaces in this case).

A. Identification of the Wigner crystal
interaction on a N 1 × N 2 = 6 × 9 kagome plaquette corresponding to ν = Npart N 1 N 2 = 1/9 filling factor. The PCD is made continuous by replacing each site by a Gaussian (see the Appendix A 2). Because our system is a torus, we repeat the plaquette to make the pattern in the PCD more visible. The red triangles mark the position of the fixed electron and its periodic images. Each maximum of the PCD corresponds to one particle forming the WC. There are N part = 6 particles at each plaquette giving five maxima and one fixed particle. They are arranged in a hexagonal crystalline lattice with lattice vectorsã 1 = [6, 0],ã 2 = [3, 3 √ 3] and its Wigner-Seitz unit cell is marked by a white solid hexagon. As a comparison, the unit cell of the underlying kagome lattice defined by the lattice vectors a 1 = [2, 0], a 2 = [1, √ 3] is shown by a yellow solid hexagon, which is nine times smaller, three times in each vector direction.
The crystallization can be confirmed by looking at the plot of Cartesian Fourier transform G c and angular Fourier transform G a , Fig. 1(b) and Fig. 1(c) respectively. In Fig. 1 there is a strong peak at zero frequency, which is the average value of the PCD. Around, there is a number of peaks arranged in a hexagonal lattice, whose lattice vectors areb , reciprocal lattice vectors toã 1 andã 2 , in agreement with the pattern shown in Fig. 1(a). The peaks further away from the origin are weaker because the particles are not perfectly localized (see the Appendix A 3 for a detailed explanation). The shape of the Wigner crystal is also probed using the angular Fourier transform in Fig. 1(c). The k θ = 0 component is related to the value of the PCD averaged over the full angle. It is zero at r = 0, then it increases and reaches a maximum at r = L 1 /2 corresponding to the distance between the fixed electron and six nearest particles. Moreover, at this radius we also see a clear component at k θ = 6 as a result of a six-fold rotational symmetry of the Wigner crystal.
The range of the plot in the radial direction is r ∈ [0, r 0 ], where r 0 = 0.6 max(L 1 , L 2 ), marked with a white dashed circle in Figure 1, to avoid the artifacts arising from the periodic images of the fixed particle. We note that the angular Fourier transform does not always look as clear as in this case. Usually the WC will be neither a perfect hexagon nor a square, hence we would obtain several peaks at frequencies k θ = 2, 4, 6 or higher, possibly at different r values (see the Appendix A 4). Nevertheless, the highest Fourier peak will correspond to the closest symmetry.

B. Wigner crystals on kagome lattice
We move to investigate plaquettes of different size and shape. Fig. 2 (a) compares the shape of the Wigner crystal unit cells on different plaquettes of kagome lattice with screened Coulomb interaction with α = 0.5 (relatively short range interaction). We call this kind of plot a phase diagram. It contains data from a number of plaquettes with sizes from N 1 × N 2 = 4 × 5 to N 1 × N 2 = 7 × 9, each populated with N part = 6 particles.
Their positions on the plot denote their filling factor ν = Npart N 1 N 2 (horizontal axis) and aspect ratio A = N 2 N 1 (vertical axis). The blue shapes are the Wigner-Seitz cells of the Wigner crystal. The N 1 × N 2 = 6 × 9 plaquette described in the previous Subsection is situated at ν ≈ 0.11, A = 1.5. It can be recognized by a perfectly hexagonal unit cell, although   The cross denotes a liquid phase with S being too small to be visible.
In Fig. 2(a) it can be seen that the strength of crystallization increases with decreasing filling factor. On the smallest plaquette, N 1 × N 2 = 4 × 5 (ν = 0.3 and A ≈ 1.25), we observe a state with nearly uniform PCD, which we interpret as a liquid. On the largest plaquette considered in this phase diagram, N 1 × N 2 = 7 × 9 (ν ≈ 0.095 and A ≈ 1.28), the Wigner crystal is the strongest. We do not observe clear liquid-crystal threshold filling factor but this can be related to finite size effects that will be discussed later. The strength of the crystallization depends on the aspect ratio. The WC for N 1 ×N 2 = 6×9 plaquette (ν ≈ 0.111, A = 1.5) is stronger than the one on N 1 × N 2 = 7 × 8 plaquette (ν = 0.107, A = 1.14) although the filling factor of these two is similar. A possible origin of such a dependence is the preference for the hexagonal WC. The perfectly hexagonal unit cell is allowed by the boundary conditions on plaquettes with A = 1.5, for example the N 1 × N 2 = 6 × 9 one.
Indeed, this plaquette has a second strongest WC, hence we can interpret the plot as if this aspect ratio was optimal, i.e. yielding the highest S for fixed ν. Although the N 1 ×N 2 = 7×9 plaquette with A = 1.28 yields a stronger WC, this may be attributed to the general trend of S increasing with the decrease of filling factor. Fig. 2(b) shows the predictions of the WC shape from minimization of the classical energies of point like particles with short range interaction V SC 0.5 by comparing all the Wigner crystals allowed by the boundary conditions. The details of the procedure are described in the Appendix B. There is a good agreement between the resulting WC shapes and the ones obtained from ED, shown in 2(a). We note that in the case of L 1 = L 2 the ground state of the classical model is degenerate. If the degeneracy exists also on the ED level, the Wigner crystallization would not be detected using the product of Fourier peaks. Hence, we decided to exclude the L 1 = L 2 plaquettes from the phase diagram and analyze them separately in the Appendix C.
When we increase the range of the interaction, the strongest WCs deviate from the hexagonal shape. Similar effect is seen also for the logarithmic interaction. For both shortand longer-range V Log we get a good match between classical and ED results. However, for V SC the agreement deteriorates when the screening is decreased. Nevertheless, the shape of the strongest WCs is still the same as predicted classically (see the Appendix D 1).

C. Wigner crystal on other lattices
In Fig. 3 we analyze the liquid -crystal transition on all three lattices: (a) kagome, (b) honeycomb, (c) checkerboard. The crystallization strength is now measured by the angular transform by computing the Fourier components at k θ = 2, 4, 6 and choosing the value of the strongest one. This value is normalized by dividing it by the maximum value of k θ = 0 Fourier component within the range r ∈ (0, r 0 ), with r 0 = 0.6 min(L 1 , L 2 ) as defined previously. Clearly, k θ = 4 and k θ = 6 corresponds to square and hexagonal WCs, k θ = 2 describes WCs elongated in one of directions. Since for some plaquettes we obtain a stripe ordering, which is not rotationally invariant and hence has nonzero angular Fourier components, we marked the plaquettes with empty symbols corresponding to no clear Wigner crystal. We consider interaction V SC 0.3 , which has slightly larger range in comparison to previous results with α = 0.5, because on a checkerboard lattice shorter-range interactions lead to appearance of PCD patterns other than WC at low filling factors (see the Appendix   6x9  5x9  5x8  5x7  4x8  4x7  4x6  4x5   7x8  6x8  6x7  5x6   7x9  6x9  5x9  5x8  5x7  4x8  4x7  4x6  4x5   7x8  6x8  6x7  5x6   7x9  6x9  5x9  5x8  5x7  4x8  4x7  4x6  4x5   7x8  6x8  6x7  The results are also comparable for short-range and long-range logarithmic interactions. We observe significant differences between the WCs on both lattices only if we consider the Coulomb interaction with small screening. A more detailed description of the results for different interaction parameters is presented in Appendix D 1. We note that Wigner crystallization in a presence of kagome or honeycomb lattice (pinning arrays) was considered for vortices in a superconductor [78,79]. These vortices behave like classical particles and significant differences in a crystallization pattern are observed between the kagome and the honeycomb lattices. However, the setup considered in Refs.
78, 79 allows also the particles to locate at interstitial positions, which is not possible in our models. Additionally, they considered filling factors much larger than in our work, leading to much larger Wigner lattice constant. As the Wigner lattice constant grows, the influence of the lattice decreases, because the particle positions become less discretized. We note that this may be the reason why we do not observe significant lattice effects. However, it is important to emphasize that we investigate small system sizes, much smaller than in Refs. 78,79, and also small number of particles, thus we do not rule out the possibility of the existence of larger differences between the lattices for larger systems.
The WCs on the checkerboard lattice ( Fig. 3(c)) differ from the ones on two other lattices.
This stems from the fact that its Bravais lattice is square rather than hexagonal, hence the shape of the plaquettes is different. This results in a different set of WCs allowed by the boundary conditions. At low filling factors, hexagonal WCs are the strongest, but elongated hexagonal WCs appear also, as a nearly-regular hexagon can not be fitted in some plaquettes

D. Finite-size effects
To investigate the dependence of the Wigner crystallization on particle number, we consider systems with N part different than 6. In Appendix E, plaquettes with N part = 4, N part = 5 particles are investigated. We find a good agreement between the classical model and ED results even for long range Coulomb interaction. In general, these results are consistent with the ones for N part = 6 particles. It is important to note that the Wigner crystals allowed by the boundary conditions are different for every value of N part . This means that our results depend strongly on the geometric factors. For example, the optimal aspect ratio to fit a hexagonal WC with N part = 4 is 1, not 1.5 as in case of N part = 6. Now, we want to analyze liquid-WC transition regardless of the shape of WC. To find out how the Wigner crystallization is affected by the finite size effects, we compare the results for N part = 4, 5, 6 described above and complement them also with results for N part = 7. In  Figure 4 shows the crystallization strength S on kagome lattice for the Coulomb interaction V SC 0.5 . It can be seen that the curves corresponding to different particle numbers have a similar behavior, increasing with lowering a filling factor. The rapid increase of the crystallization strength S with decreasing filling factors ν starts to occur at ν ≈ 0.15, i.e. close to ν = 1/7, although the curves for N part = 6, 7 are shifted towards lower filling factor with respect to curves for N part = 4, 5. We note that the plot for N part = 7 ends at plaquette 6 × 9, with relatively high filling factor ν ≈ 0.13. This is because on the plaquettes 7 × 10 and 7 × 11, which are closest to 6 × 9 in terms of aspect ratio from all the N 1 = 7 plaquettes, we do not observe the Wigner crystallization. We interpret this result as a signature of the sensitivity of the Wigner crystal made of 7 particles to the aspect ratio of the plaquette. This may be connected with the fact that one cannot realize a nondegenerate hexagonal Wigner crystal with 7 particles.
The analysis of finite size effects for other lattices and for the long-range potential V SC 0.0 is presented in Appendix F. The behavior of the S vs. ν curves is similar to what is shown in Figure 4. We note that neither in Fig. 4 nor in results in Appendix F we do not observe the liquid-crystal transition becoming more abrupt as the number of particles increases.
However, this does not necessarily mean that in the thermodynamic limit the transition will be continuous. We note that the numbers of particles investigated by us are rather small. Moreover, the behavior of the Wigner crystal depends strongly on the geometry of the sample. Thus, the reliability of the extrapolation to the infinite system is limited. Our results do not allow to determine whether the continuous nature of the transition persists in the thermodynamic limit, or is just a consequence of the small size of investigated system.
We note that the finite size effects can influence not only the profile of S vs. ν curves, but also the shape of the Wigner crystals. We analyze this effect in Appendix F. Also, we do not rule out the possibility that there are effects which are not captured by our calculation due to the small size of plaquettes. For example, it might occur that structural changes in the Wigner crystal can happen for larger systems and that the phase diagrams of larger systems are richer than the ones we obtained.

E. Band topology
To check how the band topology influences our results, we compared the Wigner crystallization of N part = 6 on trivial and nontrivial Haldane model. We have found no significant differences between these two cases (see the Appendix G). This can be contrasted with earlier results for ν = 1/3 and ν = 2/3, where the topology is important in the description of the system, as the phase diagram contains both charge ordered and topologically ordered phases [71,73,74], however we consider lower filling factors, where FCI phases are less stable. We think that the WC-to-FCI transition can be triggered by modifying the interaction, in analogy to varying the pseudopotential parameters in FQHE.

IV. SUMMARY AND CONCLUSIONS
In summary, we have shown that the Wigner crystallization occurs in topological flat bands for low particle densities in all three considered lattice models and with a variety of interaction parameters determining the interaction range. The Wigner crystallization strength increases smoothly with decreasing filling factor. In our finite-size calculation, the WC shape depends strongly on the size and shape of the plaquette and the number of particles, which determine WCs allowed by the boundary conditions. The WC shapes were to a large extent independent on the details of the lattice type and followed the predictions made by comparing the classical energies of crystals of point-like particles in a continuous space.
The underlying lattice is important only for certain aspects of the Wigner crystallization, such as the phase diagram for unscreened Coulomb interaction and the WC deformations on checkerboard lattice.
We do not observe a sharp threshold below which the crystallization starts, but this can be related to finite size effects, which can not be eliminated from calculations presented in this work. However, we can summarize that in all our systems with various lattice models, particle numbers and interaction types, the strong Wigner crystals always occur at the lowest filling factors. The rapid increase of crystallization strength with decreasing filling factor starts at filling ν = 1/7 or higher. Also, we note that the agreement between the classical model and ED results exists despite the finite size effects. If it persists in the thermodynamic limit, the resulting Wigner lattice for an infinite system with an interaction V SC will be hexagonal [80,81].
We have found no significant influence of band topology on the formation of the Wigner crystals. This is in contrast to earlier results obtained for ν = 1/3 and ν = 2/3 with shortrange interaction and is consistent with the observation that the long-range interaction usually destroys the FCIs. The single-particle Hamiltonian of the kagome model [12] reads where c † i (c i ) is a creation (annihilation) operator at site i, , , denote the first and the second neighbors, respectively, t 1 and t 2 are the real parts of first and second neighbour hoppings, λ 1 , λ 2 are their imaginary parts, and ν ij = ±1 depending on the direction of hopping (see Fig. 5(a)).

The Hamiltonian of the Haldane model [1] is
where t 1 and t 2 are magnitudes of the first and second neighbor hoppings, respectively, φ ij = ±φ is a complex phase with a sign depending on the direction of hoppings, shown in Fig. 5(b). ǫ i ± ǫ is the staggered onsite potential, +ǫ on red sublattice and −ǫ on the blue one.
The checkerboard model [11,12] is described by the Hamiltonian where t ′ ij = ±t 2 depends on the sublattice and the direction of the hopping, as indicated in Fig. 5(c), t α , with α = 1, 2, 3 denoting the absolute values of αth-neighbor hopping. The nearest-neighbor hopping contains a complex term with a phase φ ij = ±φ, where the sign corresponds to clockwise or counterclockwise direction of the hopping.
In all three models, the parameters can be tuned so that the lowest band is topologically nontrivial with |C| = 1 and nearly flat [11][12][13]. In the course of this work, we use for kagome model t 1 = −1, t 2 = 0.3, λ 1 = 0.6, λ 2 = 0, t 1 = 1, for honeycomb model t 2 =  We consider finite-size systems in torus geometry, i.e. we investigate finite plaquettes of N 1 × N 2 unit cells with periodic boundary conditions. The lattice is defined by lattice vectors a 1 , a 2 , so the dimensions of the plaquette are L 1,2 = |a 1,2 |N 1,2 . For all the lattices we consider, we have |a 1 | = |a 2 | = a. The scale of |a 1,2 | is determined by the distance d between the nearest neighbor sites, which we fix to be d = 1.

Pair correlation density
Having obtained the ground state |ψ , we calculate the pair correlation density (PCD) defined in the discrete basis of sites, describing probability of finding a particle at site j assuming that there is a fixed particle at site i. We make it continuous by replacing every site by a Gaussian, where r is the vector connecting atom i and a given point in space, i.e. we take the site i as the origin of our coordinate system, and σ is the width of the Gaussian, which we choose to be σ = 0.5. The choice of starting site i does not affect the results significantly, as the exact-diagonalization eigenstates are translationally invariant. To find the Wigner crystal, we discretize this function on a Cartesian or polar grid and perform the Fourier transform using the Fast Fourier Transform algorithm.

The Cartesian Fourier transform
If we choose the Cartesian grid, we perform the Fourier transform in both directions and obtain the Fourier coefficients where P denotes the area of the plaquette, and k is the wave vector. Because the system is periodic, the k vectors can have only discrete values k = p N 1 b 1 + q N 2 b 2 , with b 1,2 being the reciprocal lattice vectors corresponding to the real-space lattice defined by a 1,2 , and p, q being arbitrary integers.
The Wigner crystal is defined by lattice vectorsã 1,2 . Because our system is a finite-size torus, only a subset ofã 1,2 vectors is allowed by the boundary conditions. Moreover, since we fix the number of particles N part , the number of PCD maxima within the plaquette should be equal to N part − 1. Otherwise, the state is not a Wigner crystal but another charge ordering.
Ideally, the Wigner crystal would consist of point particles arranged in a lattice, with PCD where with m, n being arbitrary integers and δ(r) being the Dirac delta. The delta at the origin is subtracted because the fixed particle is not included in the pair correlation function.
The Fourier transform of G 0 would be an infinite sum of periodically arranged Dirac deltas, whereb 1,2 are the reciprocal lattice vectors of WC, each of them given by a pair of two Not every choice ofp i ,q i is permitted, as they should yield a correct number of PCD maxima.
The Fourier transforms we obtain in ED calculations are not as ideal as G c 0 (k) for two reasons. First, the particles have finite spatial dimensions. This can be seen on a simple example of particles described by Gaussians of width σ W . Then, the PCD will be a convolution The Fourier transform of G G0 is a multiplication of G 0 (k) and a Gaussian in a momentum Therefore, the spatial delocalization makes the Fourier peaks decay with increasing distance from the origin -an effect which is visible in Fig. 1(b) of the main text.
Another source of distortion from the ideal periodic pattern is the fact that the fixed particle is not included in the pair correlation density. The subtracted delta in Eq. A6 and Gaussian in Eq. A7 will give rise to additional Fourier components at k vectors not belonging to the reciprocal lattice of the Wigner crystal. Similar efect is observed in our numerical results. The spurious Fourier components are visible as the bright "cloud" around the origin in Fig. 1(b) of the main text.
We use the magnitude of the Fourier peaks as the measure of the strength of the Wigner crystallization. The WC has to be periodic in two directions, hence we should observe at least two nonzero peaks. Therefore we choose our measure to be a product of two peaks 1,2 are the two reciprocal lattice vectors defining the WC of a given type indexed by i. If the PCD is non-periodic in at least one direction, this product will vanish. We do not know which WC will be present on which plaquette. Therefore, we first list all the possible N W Wigner crystals and their lattice vectors. For example for N part = 6 particles on kagome lattice N W = 8. Since the dimensions L 1,2 of the plaquettes differ, these vectors will be different at each of them. Nevertheless, they will be defined by the same (p,q) pairs. To determine which WC is present on the plaquette, we check which pair of reciprocal lattice vectors gives the highest product S i of the Fourier components. This product is then taken as the crystallization strength S.
Several comments need to be made here. First, to compare the results for different plaquettes, the Fourier spectrum has to be normalized, which is done by dividing it by the k = [0, 0] component. Secondly, the "holes" in the PCD corresponding to the fixed electron may introduce nonzero Fourier components atb (i) 1,2 defining the Wigner crystals even if there is in fact no WC. Indeed, some of the small unit cells in Fig. 2(a) of the main text do not correspond to WCs. However, if strong WC is present, the peaks due to WC will dominate over the spurious Fourier components, as can be seen in Fig. 1(b)

The angular Fourier transform
Another choice of discretization of G(r) is the polar grid. Then, the Fourier transform is taken only along the angular direction, and the Fourier components are given by where k θ is the angular frequency. The k θ = 0 component is related to the average PCD at radius r, while all the others allow to distinguish the lattice symmetry. In the case of a nearly-hexagonal or nearly-square WC, the Fourier transform will contain a strong component at k θ = 6 or k θ = 4, respectively. As noted in the main text, it would occur at the radius equal to the distance between the first particle and the six or four nearest particles. At this radius, the zeroth component would exhibit its first maximum.
The transform is not meaningful at large r. The "holes" in PCD due to the presence of periodic images of fixed electron introduce at least 2-fold rotational symmetry and therefore nonzero Fourier component even for perfectly isotropic liquid state. Therefore, we have to introduce a cutoff r 0 . Strictly speaking, the influence of the periodic images of fixed electron starts at half the distance to the closest of them, i.e. r = 0.5 max(L 1 , L 2 ). However, we note that often a particle is located at this distance or even further, therefore the cutoff has to be slightly larger. We choose r 0 = 0.6 max(L 1 , L 2 ).
Also, we note that the angular Fourier transform does not always look as clear as in Fig. 1(c) of the main text (see Fig. 6). If the Wigner lattice is not close to neither hexagonal nor square symmetry, we would obtain several strong Fourier components at even frequencies (the odd components will vanish at least approximately because all the possible Wigner lattices have a 2-fold rotational symmetry). Moreover, if |ã 1 | = |ã 2 | the maxima of different Fourier components may occur at different radii, hence the PCD may exhibit different symmetries at different r (see Fig. 6). To determine which rotational symmetry (2, 4-or 6-fold) is the closest one, we compute the maximal value of Fourier components with k θ = 2, 4, 6 in the range [0, r 0 ]. The k θ at which the value is the highest indicates the symmetry of WC. We use this value as an alternative measure of crystallization strengthS.
However, since the magnitude of Fourier components depends on the mean particle density, we normalize it by dividing by the maximal value of k θ = 0 component in the range [0, r 0 ].
As we noted in the main text,S can be nonzero even if the system is not a WC (for example a stripe phase would also have 2-fold rotational symmetry). Therefore we have to select the WCs first, can be done visually by looking at the PCD plot, or comparing with the results for Cartesian Fourier transform.

Appendix B: Classical model
We compare the shapes of WCs obtained from the exact diagonalization calculation to predictions made using a simple classical model. The classical energy of a set of point particles is given by where the indices i and j run over all the particles, and r ij is the shortest interparticle distance on the torus. The classical prediction of the WC shape is found by calculating this energy for every Wigner lattice allowed by the boundary conditions, and choosing the one in which E is minimal. We do not take the underlying lattice into account, i.e. the particle position is not restricted to lattice sites, and is determined only by the Wigner lattice.
Such a model allows also for introduction of patterns other than the perfect crystal. We will consider several such shapes, parameterized by a single number δ (e.g. the displacement of some particle from ideal crystal positions). For each pattern like this, the energy is minimized with respect to δ and then compared with the energies of other patterns and WCs.
We note that for the logarithmic interaction, the particles may not interact classically if β is too small. Then the classical model may have several zero-energy ground states.
However, the interaction may still exist at the quantum level, possibly because the particles are not perfectly localized, and their positions are restricted to lattice sites. For example, for N part = 6 particles on kagome lattice we have a degenerate classical ground state at β < 1.82, although the ED calculations yield a nondegenerate WC even when β ∼ 1.4. Because of this effect, the exact diagonalization results cannot be compared to classical predictions for certain values of β. Such a problem is not present in screened Coulomb interaction, whose exponential tail always lifts the degeneracy.

Appendix C: Degeneracy
The plaquettes with aspect ratio A = N 2 N 1 = 1 were omitted in our analyses of N part = 6 case. This is because the ground state will always be degenerate. For example, for the plaquettes with hexagonal Bravais lattice, the N W = 8 possible Wigner crystals can be divided into two sets of WCs with the same classical energy, one consisting of six WCs, the other of two.
Indeed, the results for L 1 = L 2 honeycomb plaquettes obtained with certain interactions can be interpreted in such a way. There are six degenerate ground states, none of which yields a clear Wigner crystal in the pair correlation density. Instead, pairs of these states have similar, stripe-like PCD. This does not mean that the Wigner crystallization does not occur.
The ground state obtained in the exact diagonalization procedure may be a superposition of degenerate groundstates. We interpret each of the stripe-like patterns as two Wigner lattices superimposed (see Fig. 7). For the kagome 7 × 7 plaquettes the ground state is also 6-fold degenerate, and the sum of their PCDs has some similarities with a superposition of all six Wigner lattices. Moreover, at smaller plaquettes we obtain a similar PCD pattern, but the ground state is 3-or 1-fold degenerate. Even if these states are indeed a superposition of Wigner crystals, we cannot measure the crystallization strength, as we would have to take into account a combination of six reciprocal-space lattices. Therefore we decided to exclude the L 1 = L 2 plaquettes from our considerations.  for kagome and honeycomb lattices, respectively. One can clearly see that there are more differences between these two than between (a) and (b) subfigures. In general, the similarity betwen WCs on kagome and honeycomb lattices lowers with decreasing α. However, even if α = 0 (Fig. 8(c) and (d)), there is a considerable similarity if one limits the comparison to strong WCs only. The N 1 × N 2 = 7 × 9, 6 × 9 and 5 × 9 plaquettes (i.e. the ones with strongest WCs in (c)) yield the same shape of WC on both lattices. Decreasing α leads also to deterioration of the accuracy of the classical predictions. Nevertheless, the WC shapes on the three plaquettes mentioned above are in agreement with classical results. Also, the classical model correctly predicts that increasing the range of interaction makes the WCs at lower aspect ratios deviate from the hexagonal shape, even if the exact shape of WC unit cell does not agree with ED results.
For logarithmic interaction, such a deterioration does not happen. We investigated the logarithmic interaction on kagome lattice with β between 1.4 and 3.0 and found that at small β the WCs seem to prefer the hexagonal shape, while for higher β the WCs at small aspect ratios are closer to rectangular shape. This behavior is also well captured by the classical model, as long as β is large enough that the particles interact classically. Although the details of the transition differ in classical and ED approaches, their results agree well or even perfectly at its "end points" at high and low β. Also, we found that the WC shapes for kagome lattice are similar to the ones for honeycomb lattice for both short (β = 1.3 honeycomb, β = 1.4 kagome) and longer range interaction (β = 3.0 on both lattices).

The checkerboard lattice
The checkerboard lattice is more difficult to analyze, as, in addition to Wigner crystal, liquids and stripe patterns one observes also another type of charge ordering. We call it a "Wigner pattern" (WP) to emphasize that it consist of well-localized particles, but exhibit no periodicity other than the periodicity of the torus. In general, many Wigner patterns are possible, but in our calculations we encounter only one. We call it "half-elongated", since it resembles the half-elongated triangular tiling of the plane. It consists of rows of triangles and squares, with two rows of triangles per one row of squares, with particles located in their corners (see Fig. 9(a)). Obviously, the aspect ratio of the plaquette usually does not allow the triangles and squares to be regular polygons, so the pattern is always squeezed or stretched. Also, we observe WCs in which the particles deviate from their ideal positions in the crystal lattice, but the displacement is small enough for the Wigner lattice to be identified (see Fig. 9(b), (c), (d)). We will call these "deformed WCs".
The existence of these effects makes it more difficult or even impossible to measure the crystallization strength. The half-elongated WP cannot be described by two Fourier peaks, so we can only check visually whether it exists or not. The deformed WCs, if they are close enough to the perfect lattice, will have nonzero Fourier components corresponding to that crystal, so they may be visible using the procedure described in the main text. We have investigated the checkerboard lattice with screened Coulomb interaction with α = 0, 0.1, ..., 1 and logarithmic with β = 1.2, 1.4, ..., 3.0. For sufficiently long-range interaction the Wigner crystals are common. On three plaquettes, N 1 × N 2 = 4 × 7, N 1 × N 2 = 5 × 6 and N 1 ×N 2 = 6×7, we encounter deformations, but they are small enough for the crystallization to be seen from Fourier peaks. The shapes of WCs (including the deformed ones) are the same for both interaction types on all the plaquettes. The maximum of crystallization strength occurs again at N 1 × N 2 = 7 × 9 plaquette. When the range of the interaction is decreased, more and more WPs and/or deformations start to appear, starting from low fillings and low aspect ratios. Also, for a small number of plaquettes with low fillings, we observe a charge ordering which is neither WC nor WP, as it does not correspond to six well-localized particles. Figure 9(a) shows a comparison of the classical prediction of particle positions with the exact-diagonalization PCD for a N 1 × N 2 = 7 × 9 checkerboard plaquette with V Log 1. 6 . A good agreement between those two results is seen. In general, the classical model correctly describes the emergence of the half-elongated WP at the qualitative level. For longer-range interaction it predicts no WPs. They emerge, starting with high fillings and low aspect ratios, when the interaction range is decreased. On the quantitative level, the model does relatively well for the screened Coulomb interaction V SC . For example, for α = 0.9 and α = 1.0 the classical model predicts half-extended WP on five plaquettes (N 1 × N 2 = 7 × 8, 7 × 9, 6 × 7, 6 × 8, 5 × 6), in four of which it exists also in quantum results (all the above except N 1 × N 2 = 5 × 6). For logarithmic interaction its performance is worse. For example, for V Log 1.6 it predicts half-elongated WPs at four plaquettes (N 1 × N 2 = 7 × 9, 6 × 8, 6 × 9, 5 × 9), while in exact diagonalization it exist on three (N 1 × N 2 = 7 × 8, 7 × 9 and 6 × 9), and only two are guessed correctly.
On the other hand, the classical model fails to describe the deformed WCs. This can be seen in Fig. 9(b) and (c). In both subfigures, the classical model predicts no deformation, although they exist on the ED level. Similar behaviour is observed in the case of longerrange V Log , and V SC regardless of α. For short-range V Log , the model predicts too many deformations. Although in several cases it correctly predicts their shape ( Fig. 9(d)) usually the prediction is wrong. This suggests that the deformations arise rather due to the presence of the lattice. Also, we note that the deformation of the type shown in Fig. 9(b) exists only when N 1 is odd (N 1 × N 2 = 5 × 6, 5 × 8, 5 × 9, 7 × 9 plaquettes) while the one in Fig. 9(c) only for even N 1 and odd N 2 (N 1 × N 2 = 4 × 7, 6 × 7). This suggests a commensuration effect, although the number of plaquettes is too small to determine it.
We have investigated the same plaquettes as described above filled with 4 or 5 particles.
When the number of particles is changed, different Wigner crystals are allowed by the boundary conditions. However, they still follow, to large extent, the behavior of classical particles. Figure 10 shows the phase diagram for kagome lattice with N part = 4 and V SC 0.5 along with the classical predictions. Note that the L 1 = L 2 plaquettes are now included, because they do not yield degenerate Wigner crystals. The Wigner-Seitz cells of the WCs tend to be close to hexagonal for low aspect ratio (with a perfect hexagon for aspect ratio 1), while for higher aspect ratio they deviate from this shape. The agreement between classical and ED results is good. We have investigated N part = 4 on kagome and honeycomb lattice provided that the interaction is sufficiently long-range so that it does not yield degenerate ground states. It is perfect or nearly perfect (at most one plaquette predicted wrong) for logarithmic interaction, and slightly worse for the screened Coulomb potential, where typically there are two or three plaquettes where the predicted shape was different from the one in ED.
On the checkerboard lattice, we do not encounter any Wigner patterns, but the deformations of WCs are present. Again, we try to parameterize them using a single parameter and include in the classical model. However, the predictions obtained in such a way do not reproduce the ED results. Moreover, we again note that there are two types of deformations which tend to occur mostly when N 1 is even and N 2 is odd, and vice versa. This strengthens our suggestion that this is a commensuration effect, and at least some deformations are due crystal is the same regardless of interaction parameters in the whole range we investigated (α ∈ [0, 1], β ∈ [1. 4,3], changing by 0.2) and is predicted by the classical model with 100% accuracy. What is interesting is also the disappearance of Wigner crystals at higher aspect ratios A for V Log β≥2.0 . The WCs are not replaced by Wigner patterns, but rather by stripe-like PCD patterns. 5 particles on checkerboard lattice are much more difficult to analyze, as every possible Wigner crystal is two-fold degenerate due to reflection symmetry.
Indeed, for some plaquettes and some interaction parameters we observe a PCD which can be interpreted as two such WCs superimposed. Also, we find PCDs which may be a superposition of degenerate WPs or deformed WCs. Due to the degeneracies, we decide to exclude the 5-particle checkerboard cases from our analysis.

Appendix F: Finite-size effects -details
To gain some insight on the finite size effects, we compare the results for different particle numbers. Figure 11 shows comparison of the crystallization strength vs filling factor plots for  four cases: (a) kagome lattice with short-range interaction V SC 0.5 , (b) the honeycomb lattice with long-range interaction V SC 0.0 , (b) the honeycomb lattice with long-range interaction V SC 0.0 , (d) the checkerboard lattice with long-range interaction V SC 0.0 . The results in all the subfigures of this Figure involve the results for N part = 4, 5, 6, described in the previous Appendices and in the main text. Additionally, for kagome lattice we performed the calculation with N part = 7, whose results are included in Fig 11 (a). Also, we note that in Fig. 11 (d) we plot only two curves, as the N part = 5 case leads to degeneracy on the checkerboard lattice, and that for this lattice we study only the long-range interaction, as the short-range one leads to the presence of Wigner patterns at N part = 6.
The results shown in all four subfigures of Fig. 11 subfigures show an agreement between the crystallization strengths obtained for different particle numbers. This agreement is better for checkerboard ( Fig. 11 (d)) and honeycomb ( Fig. 11 (b) and (c)) lattices than for kagome lattice (Fig. 11). We do not observe the transition getting more sharp as the system size increases. However, as noted in the main text the extrapolation to the thermodynamic limit cannot be performed reliably, especially when the result depend strongly on sample geometry.
The finite-size effects influence also the shape of the Wigner crystal. It is difficult to investigate this effect systematically, as the boundary conditions rarely allow the formation of the Wigner crystals with the same shape and with different N part . We have such a possibility only on three pairs of plaquettes: (i) 7 × 6 with N part = 4 and 7 × 9 with N part = 6, (ii) 6 × 6 with N part = 4 and 6 × 9 with N part = 6, (iii) 5 × 6 with N part = 4 and 5 × 9 with N part = 6. Figure 12 shows the results for pair (i) for kagome lattice and V SC interaction. In Fig 12 (a) we plot the pair correlation density for a 6 × 7 plaquette with N part = 4 with V SC 0.5 interaction. The white shape is the Wigner-Seitz cell of the Wigner crystal. This result can be compared with Fig 12 (b), which shows the pair correlation density for the 7 × 9 plaquette with N part = 6. The unit cell of the Wigner crystal is the same as in Fig. 12 (a), suggesting that the finite-size effects do not influence the shape of the Wigner crystal. The situation becomes different when we consider the unscreened Coulomb interaction. In such a case, for the 6 × 7 plaquette with N part = 4 we obtain a PCD pattern indistinguishable from the one in Fig 12 (a). However, for the 7×9 plaquette with N part = 6, we obtain the PCD shown in Fig 12 (c), different from the one in Fig 12 (b). Thus, for the long-range interaction the finite-size effects are stronger and influence the shape of the WC. A similar behavior is seen for pair (ii): the same unit cell of WC is obtained on both plaquettes if the interaction range is short, but for the long range they become different. On the other hand, for pair (iii) we get different WCs on the two plaquettes for both short and long range interaction. That is, even for the short-range interaction the finite-size effects can influence the shape of the Wigner crystal. Similar results are obtained also for honeycomb lattice. Thus, we conclude that on kagome and honeycomb lattices we observe strong finitesize effects for long range interaction and moderate finite-size effects for short-range one.
For checkerboard lattice, we obtain strong finite-size effecs for both short-and longrange interaction. For short-range interaction, the Wigner crystals is the same on the two plaquettes only in pair (iii). On the two other pairs, the larger plaquette usually contains the Wigner pattern, which cannot be present on the smaller one. For long-range interaction, in all three pairs we obtain a rectangular WC on the smaller plaquettes and a more hexagonal one at the larger plaquettes. The only exception is the logarithmic interaction V Log 1.8 , for which we get a WC with the same unit cell on each pair of plaquettes.
Thus, we conclude that our calculation is prone to finite-size effects. We cannot perform a reliable extrapolation to thermodynamic limit neither for the crystallization strength nor for the crystal shape. We note that the finite-size effects related to WC shape are present also on the classical level -for example, while in the N part = 6 case the classical model predicts a non-hexagonal WC for unscreened Coulomb interaction, the infinite-plane classical Wigner crystal is hexagonal [80]. On the other hand, as the similarity between the classical and quantum results exists for all the system sizes we investigated regardless of the number of particles and interaction types, it is possible that it will hold also in the thermodynamic limit. This does not have to be the case, as it may occur that, for example, the lattice effects will be more visible as N part increases. However, if it is, and if Wigner crystal exists in the thermodynamic limit, we can expect that it will be hexagonal for both screened and unscreened Coulomb interaction basing for the infinite-plane results [80,81].
Appendix G: Comparison with trivial system There are two reasons to suspect that the topological properties of the flat bands may affect the Wigner crystallization. The first is the possible occurrence of FCIs on these lattices. In the course of our analysis, several plaquettes allowed for the occurence of the Laughlin fillings ν = 1/5 or ν = 1/7. At some of them, for certain interaction parameters, the lowest energy states obey the FCI counting rules [82]. Nevertheless, for most of them the pair correlation density is not uniform, it is either WC, WP, a stripe pattern or a different, but non-uniform charge ordering. The only cases in which we are not able to disprove the presence of an FCI by looking at the pair correlation density are kagome 5 × 5 plaquettes with 5 particles (ν = 1/5) for some values of interaction parameters. However, this plaquette allows for a degenerate WC and hence is excluded from our analyses. Moreover, even if this state is an FCI, it is not a stable one, as we do not observe it for similar interaction on other ν = 1/5 plaquettes. Therefore, we can neglect the presence of FCIs in our analysis.
The second reason are the constraints on particle localization forced by nontrivial topology. It is impossible to localize the Wannier function in both dimensions if the Chern number is nonzero [83]. Therefore, it may mean that the Wigner crystallization in the trivial lattice would be stronger. To check this hypothesis, we have performed the calculation for trivial honeycomb system with nonzero staggered potential ǫ = 0.15. We have chosen four data points representing a shorter-and longer-range version of both interactions: V SC 0 , V SC 0.5 , V Log 1.3 , V Log 3.0 . Comparing the phase diagrams with the ones of nontrivial honeycomb lattice discussed previously, we discover that the shapes of Wigner crystals are exactly the same, and there were only minor changes in the crystallization strength. Therefore, we conclude that the Chern number of the flat band has no significant effect on the Wigner crystallization. We note that this is not the effect of the band mixing due to strong interaction, as all the results are obtained using the band-projected exact diagonalization.