Signatures of many-body localization of quasiparticles in a flat band superconductor

We construct a class of exact eigenstates of the Hamiltonian obtained by projecting the Hubbard interaction term onto the flat band subspace of a generic lattice model. These exact eigenstates are many body states in which an arbitrary number of localized fermionic particles coexist with a sea of mobile Cooper pairs with zero momentum. By considering the dice lattice as an example, we provide evidence that these exact eigenstates are in fact manifestation of local integrals of motions of the projected Hamiltonian. In particular the spin and particle densities retain memory of the initial state for a very long time, if localized unpaired particles are present at the beginning of the time evolution. This shows that many-body localization of quasiparticles and superfluidity can coexist even in generic two-dimensional lattice models with flat bands, for which it is not known how to construct local conserved quantities. Our results open new perspectives on the old condensed matter problem of the interplay between superconductivity and localization.


I. INTRODUCTION
Strongly correlated quantum many-body systems remain a foremost challenge in physics and hold the key to understand fascinating phenomena such as hightemperature (high-T c ) superconductivity.[1] The microscopic origin of high-temperature superconductivity remains a topic of active research because it is still very difficult to accurately simulate the model Hamiltonians that are believed to describe the relevant low-energy degrees of freedom of high-T c superconducting materials.One of these model Hamiltonians, the Fermi-Hubbard model [2] which describes copper-based high-T c superconductors (so-called cuprates), has become the favorite test subject in the field of strongly correlated systems.[3,4] Numerical exact diagonalization is particularly challenging in the case of the Hubbard model with repulsive interactions, due to many competing orders.[5][6][7] The type of order that manifests in the ground state is sensitive to the size of the finite cluster used in exact diagonalization, making the extrapolation to the thermodynamic limit problematic.On the other hand, quantum Monte Carlo suffers from the sign problem in the case of repulsive interactions.[3] To overcome these difficulties, new approaches have been developed, such as the quantum simulation of the Hubbard model using ultracold gases in optical lattices [8,9], in addition to new advanced numerical methods based on tensor network states, [10,11] for instance.Recent numerical results indicate that superconducting long-range order does not seem to occur in the repulsive Hubbard model on a simple square lattice with only nearest-neighbor hoppings.[12] Thus, to properly account for the superconductive properties of copper-based superconductors it becomes important to consider more realistic model Hamiltonians, for instance, by augmenting the pure Hubbard model with next-nearest-neighbor * sebastiano.peotta@aalto.fiSchematic of the dice lattice.Six-fold coordinated sites (hub sites, orbital index α = 1, 2) are denoted by hexagons, while three-fold coordinated sites (rim sites) are denoted by triangles.The rim sites can be further divided into two triangular sublattices denoted by up-(α = 3, 5) and down-pointing triangles (α = 4, 6), respectively.The bonds denote nearest-neighbour hoppings which are all real and have the same absolute value t = 1.The color of the bond denotes the sign of the hopping amplitude, grey for +t and red for −t.The rectangular box is the magnetic unit cell, given by the fundamental vectors aj=1,2 (29).Two Wannier functions of the two lowest flat bands of the dice lattice Hamiltonian are also shown.The Wannier function w nl , with n = 1, 2, is centered on the hub site α = n in unit cell l and is nonzero only on the same hub site and the adjacent rim sites, therefore it is compactly localized.These "flower states"form the orthonormal basis used in the expansion of the projected Hamiltonian (7).
hoppings, [13,14] or by considering a lattice model with more than one band.In the case of cuprates, a threeband model describing both the copper d-orbitals and the oxygen p-orbitals in the copper-oxide planes of cuprates seems the most appropriate.[15][16][17][18][19][20] Multiband lattice models contain more than one or-arXiv:2302.06250v3[cond-mat.supr-con]22 Dec 2023 bital per unit cell, and are therefore more difficult to simulate numerically.On the other hand, they can harbor new qualitative effects.For instance, in a multiband lattice model, the strongly correlated regime can be achieved by reducing the bandwidth of a partially filled band to zero, obtaining a so-called flat band.[21] Crucially, due to interfering hopping paths that lead to particle localization, the vanishing of the bandwidth is not necessarily accompanied by the vanishing of the hopping amplitudes in the lattice model, as it would in the case of a single band/orbital lattice model.This means that the Wannier functions that span the flat band subspace can have a large overlap and by projecting the interaction term on this subspace one obtains a nontrivial purely quartic Hamiltonian, that is a linear combination of products of four fermionic field operators.[22,23] The range of the terms in this Hamiltonian is determined by how fast the Wannier functions decay with distance.
The Hubbard interaction term projected on the flat band subspace of a lattice model, called the projected Hamiltonian from here on, is the subject of the present work.This class of quartic Hamiltonians essentially poses a many-body problem and are in general, difficult to study.Nevertheless, it is possible to derive several interesting exact results.One example is provided in this work, in which we derive exact eigenstates of generic projected Hamiltonians in the spin imbalanced case, that is for different number of spin up and down particles N ↑ ̸ = N ↓ .This new result builds on previous related work.[23][24][25] Recently, there has been a lot of interest on the topic of superconductivity and superfluidity occurring in flat or quasi-flat bands.The original motivation is that the critical temperature of the superconducting transition is enhanced by the high density of states of a quasi-flat band [26,27] (the density diverges in the strict flat band limit, nevertheless, the critical temperature remains finite).Moreover, lattice models with flat bands have been realized with ultracold gases in optical lattices, which are a flexible platform for exploring superfluid phenomena.Two notable examples are the Lieb lattice, [28,29] which also describes the atomic structure of the copperoxide planes of cuprates, and the kagome lattice.[30,31] Proposals to realize other lattice models with flat bands, such the dice lattice, have been put forward.[32] Finally, the discovery of superconductivity in twisted bilayer graphene [33] has stimulated a lot of theoretical and experimental efforts in the new field of moiré materials, [34] in which flat bands and band structures with exotic properties are engineered through the precise twisting of layers of two-dimensional materials, such as graphene, with respect to each other.The same idea has also been implemented in ultracold gas experiments.[35] Flat bands are rather peculiar, since, in the absence of interactions they are necessarily insulating regardless of the filling.Thus, one way to characterize flat bands is as translationally invariant or disorder-free Anderson insulators.On the other hand, a partially filled flat band generally becomes a superconductor with a high critical temperature as soon as an attractive Hubbard interaction is switched on.[36][37][38][39] The reason being that two-body bound states become delocalized and their inverse effective mass is proportional to the interaction strength, and to the overlap between distinct Wannier functions of the flat band.[23,[40][41][42] Recent works [24,25] have pointed out that, even if the ground state of a partially filled flat band in the presence of interactions is superconducting/superfluid, the quasiparticle excitations can be localized, namely they have infinite effective mass.Quasiparticle excitations carry a nonzero spin and have fermionic statistics.These localized quasiparticles coexist with mobile two-body bound states, which are bosons.In other words, the dispersion of quasiparticles is flat as in the noninteracting case.Evidence of localization of quasiparticles in flat band superconductors mainly comes from the analytical results obtained in Refs.24 and 25 and from the quasiparticle dispersion obtained from the Bogoliubov-de Gennes Hamiltonian of mean-field theory, for instance, in the case of the Lieb lattice.[36] The phenomenon of localization of quasiparticles in flat band superconductors draws an interesting connection with many-body localization, a subject which has attracted considerable interest, [43][44][45][46][47][48][49][50] thanks to its observation in ultracold gas experiments.[51][52][53] Loosely speaking, many-body localization is simply Anderson localization in the presence of interactions.More precisely, it is now established that the defining property of manybody localization is the presence of an extensive number of local integrals of motion that completely block particle transport and prevent thermalization in a quantum many-body system.[47,[54][55][56][57] These local integrals of motion generally appear for large enough disorder, however, examples of disorder-free systems in which manybody localization occurs are also known.[58][59][60][61][62] Methods to design lattice models possessing local integrals of motion from the bottom up have been proposed as well.[63,64] Interestingly, for some specific one-dimensional lattice models with flat bands, it is possible to analytically construct local integrals of motion, also called conserved quantities, thereby establishing a rigorous connection between the localization of quasiparticles in flat band superconductors and many-body localization.[24] The main difference with standard many-body localization is that the mobile two-body bound states present ergodic behavior and thermalize, contrary to the localized quasiparticles.Analogous analytical results regarding local conserved quantities in the case of lattice models in two or more dimensions is lacking at present.This important open question is addressed here with an exact diagonalization study of a specific lattice model with flat bands, the dice lattice, [65][66][67][68][69][70][71][72] see Fig. 1.Through the analysis of the eigenstate spectrum and the time evolution of three-and four-particle states, we find evidence for the presence of local conserved quantities in this model.
The present work is organized as follows.In Section II, we provide new analytic results for the projected Hamiltonian of generic lattice models with flat bands.In order the keep our work self-contained, in Section II A we introduce the concept of projected Hamiltonian for the attractive Hubbard interaction and derive its exact ground state in the spin balanced case N ↑ = N ↓ (N σ is the number of particles with spin σ).The method used to obtain this result is slightly different from the original work [23,25] and is more suited to derive the new exact eigenstates in the spin imbalanced case |N ↑ − N ↓ | ≥ 1, which are presented in Section II B. In these exact eigenstates, localized quasiparticles are placed at arbitrary positions on top of a sea of Cooper pairs with zero momentum.The special case |N ↑ − N ↓ | = 1 has been already discussed in Ref. 25.The assumptions under which these results hold are presented in detail in Section II A, particularly, the so-called uniform pairing condition introduced in Ref. 23.
In Section III, a specific lattice model with flat bands, the two-dimensional dice lattice with a perpendicular magnetic field, [67,68] is introduced.For a specific value of the magnetic flux per unit cell, the two lowest bands of the dice lattice are perfectly flat and degenerate.A special property of the dice lattice is that the flat band subspace is spanned by Wannier functions that can be taken to be compactly localized, that is, nonvanishing only on a finite number of lattice sites.This property makes the dice lattice particularly amenable to analytical and numerical studies.Indeed, in Section III A, we find that the general analytical results presented in Section II remain valid in the case of the dice lattice even if the uniform paring condition is partially lifted.This is a straightforward consequence of the compact nature of the Wannier functions.In Section III B, we provide a convenient representation of the projected Hamiltonian of the dice lattice.The projected Hamiltonian can be written in terms of operators that belong to three different classes: i) localized spins, ii) on-site singlets and iii) bond singlets.The on-site singlets and the bond singlets are two-body bound states moving on triangular and kagome lattices, respectively.The projected Hamiltonian also contains a term that converts on-site singlets to bond singlets and vice versa.In preparation for the numerical results presented in the following sections, the two-body problem for the projected Hamiltonian of the dice lattice is solved in Appendix A.
Section IV presents numerical results obtained by exact diagonalization of the projected Hamiltonian of the dice lattice for particle number N ↑ + N ↓ > 2 and N ↑ ̸ = N ↓ .First, in Section IV A the spectrum of the Hamiltonian is investigated.It is found that, while the ground state is perfectly degenerate, in agreement with the general analytical results of Section II, the excited states form multiplets that are only quasidegenerate.We interpret this as evidence of local conserved quantities, with the lifting of the degeneracy possibly due to finite size effects.In Section IV B, the focus is on the nonequilibrium dynamics.The main result is that the spin and particle densities do not thermalize, but rather retain memory of the initial positions of the single unpaired particles.This is another piece of evidence for the presence of local conserved quantities enforcing quasiparticle localization.
Finally, in Section V, we summarize and discuss the main results of this work and point out interesting directions for future studies.

II. EXACT MANY-BODY EIGENSTATES IN LATTICE MODELS WITH FLAT BANDS
A. Flat band projected Hamiltonian In this section, we introduce the technique of projecting the Hubbard interaction term of a lattice Hamiltonian onto the subspace corresponding to a flat band or a group of degenerate flat bands, which are obtained by diagonalizing the noninteracting term of the same Hamiltonian.This method is just the first step in a systematic expansion known as the Schrieffer-Wolff transformation.[23] The small parameter in this perturbative expansion is U/E gap ≪ 1, where U is the coupling constant of the interaction, in our case the Hubbard interaction, and E gap is the energy gap, that is, the minimal energy interval separating the group of degenerate flat bands from all other bands.The projected Hamiltonian provides an accurate description of the low-energy degrees of freedom of the many-body system if the parameter U/E gap is small enough.This is called the isolated band limit.
Using the algebra of projected field operators, it is shown how to construct exact eigenstates of the projected Hamiltonian with an arbitrary number of Cooper pairs with zero momentum.The derivation of this result, presented originally in Ref. 23 in the case where the particle number is not fixed (grand canonical ensemble) and in Ref. 25 for fixed particle number (canonical ensemble), is repeated here with some variations since it serves as a starting point for obtaining a new class of exact eigenstates in Section II B. The notation used here is very close to that of Refs.23 and 24.
Consider a generic translationally invariant lattice model with a Hamiltonian of the form Ĥ = Ĥfree + Ĥint , where Ĥfree is the noninteracting term and Ĥint an Hubbard interaction term of the form with ĉiασ fermionic field operators.Here, we focus on the case of attractive interactions U α > 0. The vector of integers i = (i 1 , i 2 ) T labels the unit cells and α = 1, . . ., N orb labels the different orbitals inside the unit cell.For instance, the dice lattice, which is studied in detail in Sections III and IV, contains N orb = 6 orbitals in its unit cell, as shown in Fig. 1.For spin-1/2 fermions, the spin index takes two possible values σ =↑, ↓.Without loss of generality, we restrict our presentation to twodimensional systems for definiteness.
The noninteracting term Ĥfree is quadratic in the field operators and can be diagonalized by solving the corresponding single-particle problem.If we also assume that the component σ of the spin is a conserved quantity, then the noninteracting term can be written in a diagonalized form as where fnkσ are new fermionic operators associated to the eigenstates ψ nkσ (i, α) of the noninteracting (singleparticle) Hamiltonian.Due to translational invariance, these eigenstates are Bloch plane waves labeled by the band index n and the quasimomentum k = (k x , k y ) T .We assume that the bands in a given set F are all degenerate and flat, that is, We require the Hamiltonian to be time-reversal symmetric, therefore all the flat bands are spin degenerate.The group of degenerate flat bands is also isolated, that is, separated from all other bands by a finite band gap E gap .
Since we are interested in the low energy properties of the system when the flat bands in the set F are partially filled and U α /E gap ≪ 1, it is useful to introduce projected field operators [23] The projected field operators are obtained by expanding the field operators ĉiασ in the basis of Bloch plane waves and retaining only those terms corresponding to the flat bands in the set F. An alternative expansion is in terms of the Wannier functions w nσ (i − l, α) and the associated field operators dnlσ .The Wannier functions are obtained from the Bloch plane waves as with N c the number of unit cells in the lattice.[73] If, and only if, the bands are flat, will the Wannier functions be eigenstates of the noninteracting Hamiltonian.
The Wannier functions provide a convenient basis of local wave functions on which the projected Hamiltonian can be expanded, as shown below in the case of the dice lattice.An important property of Wannier functions is that a complete basis for the subspace of a given band n is obtained by taking the translates of just a single Wannier function, that is A good description of the low energy properties for small U α /E gap is given by the projected Hamiltonian, in our case is the Hubbard interaction term projected on the subspace of the degenerate flat bands.The projected Hamiltonian expanded in the Wannier function basis reads Here, the sum over the band indices is restricted to the group of degenerate flat bands (n 1 , . . ., n 4 ∈ F).
It is important to note that the projected field operators defined in (4) satisfy modified anticommutation relations where P σ is the single-particle operator projecting on the F subspace with spin σ.As a consequence of the above comutation relations, we have the identity n2 iασ = P σ (0, α, α)n iασ , with niασ = c † iασ ciασ , (11) which should be compared with n2 j = nj , valid for the usual fermionic operators ĉj , ĉ † j and n = ĉ † j ĉj .In order to construct exact eigenstates of the projected Hamiltonian, we introduce the creation operator of a Cooper pair in a zero momentum state The equivalence of the three different expansions of the operator b † shown above is a consequence of time-reversal symmetry, which for spin-1/2 particles amounts to the following relations between wave functions with opposite spin The following commutation relations are essential for what follows and are only valid in the case of timereversal symmetry, [ To prove our results, we need an additional assumption, namely that the following condition, introduced for the first time in Ref. 23, holds Here S is the set of orbitals on which the wave functions of the group of degenerate flat bands F are nonzero.Note that this condition, called the "uniform pairing condition", [23,25] can always be met by adjusting the orbitaldependent Hubbard couplings U α .Assuming time-reversal symmetry, conservation of the spin component σ and the uniform pairing condition, it is shown that b , b † are ladder operators for the projected Here we have used, in order, the commutation relations (16) and then ( 8) and ( 9) and finally the uniform pairing condition (19).It follows that the states with N p Cooper pairs in the same zero momentum state, are exact eigenstates of the projected Hamiltonian with energy −E p N p .We denote by |∅⟩ the vacuum, the state with no particles, while the positive Cooper pair binding energy E p > 0 has been introduced in (19).A simple argument, [23,25] not repeated here, shows also that these are eigenstates with minimal energy for any fixed number of pairs N p .

B. Exact eigenstates with localized quasiparticles
In this section we present a new result, that is a new class of exact eigenstates of the projected Hamiltonian, characterized by the presence of unpaired localized particles on top of a sea of Cooper pairs.The proof is very straightforward and uses the fact that the Cooper pair creation operator b † ( 12) is a ladder operator, as shown in (21).The exact eigenstates in (22) are generated by repeatedly applying the Cooper pair creation operator to the vacuum |∅⟩, which is an eigenstate of H int .Instead of the vacuum, one can start from the following eigenstates of the projected Hamiltonian where I σ is an arbitrary set of Wannier functions.The second equation above follows from the fact that the Hubbard interaction term is nonzero only if particle with opposite spins are present.The state |I σ ⟩ corresponds to an arbitrary arrangement of particles all with the same spin σ, therefore it is a trivial eigenstate of the projected Hamiltonian with eigenvalue zero.We can then apply (21) to show that the states of the form are also eigenstates of H int , that is Note that the energy eigenvalue does not depend on the set I σ , only on the number of pairs N p , exactly in the same way as for the eigenstates in (22).These states are all orthogonal and their normalization constant can be easily calculated.[25] The case in which I σ consists of a single element, that is a single unpaired particle, has already been derived in Ref. 25 using a different method.Despite the simplicity of the proof, the result in ( 26) is rather nontrivial since the exact eigenstates (25) describe mobile Cooper pairs coexisting with localized particles all with the same spin, which are arbitrarily arranged in space.These localized particles are a source of disorder for the Cooper pairs.Therefore, this result is a promising starting point for investigations on the interplay of superconductivity and disorder from a genuine many-body perspective.However, one has to first answer the question on whether unpaired particles remain localized if the Cooper pairs have a finite momentum, for instance if a supercurrent is present in the system.One way to address this question is by considering specific lattice models with flat bands, such as the dice lattice shown if Fig. 1, which is the subject of the following sections.

III. DICE LATTICE
The goal of this section is to study the Hamiltonian obtained by projecting the Hubbard interaction term on the two lowest degenerate flat bands of the dice lattice.The dice lattice with its hopping amplitudes shown in Fig. 1 has a band structure composed of six flat bands, which are two by two degenerate.[24] A convenient property of the dice lattice is that the Wannier functions of the two lowest flat bands can be taken to be compactly localized and not just exponentially decaying.Compactly localized means that a Wannier function of a lattice model is nonzero only on a finite number of sites as shown in Fig. 1.Thus, the projected Hamiltonian ( 7) contains only a finite number of nearest-neighbor terms.Due to its special properties, in particular, the existence of compactly localized Wannier functions, the dice lattice with the hopping amplitudes shown in Fig. 1 has been studied in a number of previous works, [24,[65][66][67][68][69][70][71][72] particularly within the context of superconducting wire networks and Josephson junction arrays.[74][75][76][77][78] If the sign of the the hopping matrix elements shown in Fig. 1 is ignored, the underlying Bravais lattice is triangular with fundamental vectors Here, a is the length of the edge of an elementary rhombus in the dice lattice, that is, the distance between two nearest-neighbor lattice sites.Any pair of vectors chosen from the above three is sufficient to generate the triangular lattice.However, for future calculations, it is convenient to introduce a third vector v 3 , which is a linear combination of the other two.The magnetic field is encoded in a lattice model by means of Peierls phases in the hopping amplitudes.In the case of the dice lattice shown in Fig. 1, the Peierls phases are all real, that is, they are simply a sign ±1, corresponding to the bond colors in the figure.Therefore, the magnetic flux through each elementary rhombus of the lattice is a half-flux quantum.The uniform magnetic field partially breaks the translational symmetry of the dice lattice and the magnetic unit cell has an area twice as large as the unit cell of the triangular lattice (27).The underlying Bravais lattice is rectangular and a possible choice for the fundamental vectors is For our calculations, we only need the the Bloch functions of the two lowest degenerate flat bands [24] We denote with |g nk ⟩ a vector with components g nk (α), that is, |g nk ⟩ = g nk (1), g nk (2), . . ., g nk (N orb ) T and we In the above equations, the normalization factor of the Bloch functions is given by c = 1 where is the energy of the two lowest flat bands of the dice lattice.The parameter ε h is the on-site energy of the hub sites, while the rim sites have zero on-site energy (see Fig. 1 for the definition of hub and rim sites).The unit of energy is the positive hopping between the nearestneighbors in the dice lattice.Due to the fact that the Bloch functions in (31) and (32) are polynomials in the coefficients e ±iki , the Wannier functions, obtained from (5), are compactly localized.The Wannier function generated by |g n,k ⟩ is centered on hub site α = n = 1, 2 and has nonvanishing weight only on this hub site and the adjacent rim sites, as shown in Fig. 1.

A. Uniform pairing condition and exact eigenstates in the dice lattice
The uniform pairing condition ( 19)-( 20) is required to show that the Cooper pair operator b † is a ladder operator for the projected Hamiltonian of a generic lattice (21).As a consequence of the compact localization of the dice lattice Wannier functions, it is shown here that the uniform pairing condition can be partially relaxed.Indeed, the commutation relation in (21) is valid in the case of the dice lattice if the orbital-dependent Hubbard couplings take the following form Here, U r and U h are the coupling constants of the Hubbard interaction on the rim and hub sites, respectively, and can take arbitrary values.For the dice lattice, the diagonal elements of the projector P are given by Therefore, under the conditions ( 35)-( 36), the uniform pairing condition is satisfied for the sublattices formed by the hub and rim sites separately, but is not generally satisfied if α = 1, 2 and β = 3, 4, 5, 6 in (19).
The key result that allows us to modify the derivation in (21) and make it work for arbitrary values of the coupling constants U r and U h is the following Note how the sum over the orbitals is restricted to the hub sites only.The above equation is a consequence of the fact that distinct Wannier functions have zero overlap on the hub sites.From (38), it is easy to show that the commutation relation [ H int , b † ] = −E p b † holds with a modified expression for the Cooper pair binding energy E p for arbitrary values of U r and U h .
Another consequence of the compact nature of the Wannier functions of the dice lattice is that it is possible to construct an even wider class of eigenstates similar to (25) with the difference being that the localized particles can have both up or down spins.This is the case since particles with opposite spin localized on nonoverlapping Wannier functions do not interact.Thus, the following states are eigenstates with eigenvalue zero of the projected Hamiltonian under the conditions that all the Wannier functions in the set I ↑ have zero overlap with the Wannier functions in the set I ↓ .However, this is not true in general for arbitrary lattice models, since the Wannier functions have usually exponential tails and particle with opposite spins do interact slightly even if far apart.We can then apply the ladder operator b † to the above states and obtain the following exact eigenstates of the dice lattice projected Hamiltonian The energy eigenvalue of this state is given by −E p N p as in (26) and is independent from the number, position and spin of the localized particles.

B. Projected Hamiltonian of the dice lattice
Here, we provide the projected Hamiltonian of the dice lattice in a form which is rather compact and convenient 1 2 FIG. 2. Graphical representation of the term Ĥtri.+ Ĥkag. in the projected Hamiltonian Hint (48).The Wannier functions w n=1,l and w n=2,l are represented by the green and red hexagons, respectively.One unit cell (shown as the blue rectangle) contains one green site and one red site, corresponding to the two nonequivalent Wannier functions.The dashed lines connecting green and red sites with each other represent the nearest-neighbor interaction and hopping terms of on-site pairs and on-site spins in Ĥtri.(49).In particular the on-site pairs hop between nearest-neighbors in the triangular lattice composed by green and red sites.The black sites sit on the middle of the bonds connecting the green and red sites and form a kagome lattice.The operator B+ ⟨n 1 l 1 ,n 2 l 2 ⟩ creates a singlet on the bond ⟨n1l1, n2l2⟩, thus the bond singlets live on the kagome lattice formed by the black dots.The black and red bonds represent terms of the form B+ in the bond singlet hopping Hamiltonian Ĥkag.(50).The sign of the hopping amplitude of the bonds connecting the black sites is given by −s(n1l1|n2l2, n3l3) (see Eq. ( 51)) and is equal to −1 for the black bonds and to +1 for the red bonds.The terms of the form B+ for subsequent considerations and computations.To this end, we introduce three different sets of operators that are linear combinations of products of two Wannier function operators dnlσ and d † nlσ .The first is the set of on-site spin operators, that is, the spin operators relative to a single Wannier function.The on-site spin operators of the Wannier function labeled by nl are defined as In the first equation above, we have introduced the spinresolved occupation number operator ρnlσ = d † nlσ dnlσ of a Wannier function.For later use, it is convenient to also introduce the occupation number operator summed over the spin components as ρnl = ρnl↑ + ρnl↓ .
The second set is composed of operators that create and annihilate a pair of particles with opposite spins on the same Wannier function.These operators have been introduced, for instance, in Ref. 24 and are defined as We call these operators the on-site pair operators (or on-site singlets) and it is easy to check that they obey the same SU(2) algebra as the on-site spin operators.Moreover, the on-site spin operators and the on-site pair operators commute with each other: [ Ŝα nl , Bβ nl ] = 0 for α, β = x, y, z.
The third and final set is composed of operators of the form The Wannier functions are centered on the hub sites of the dice lattice and thus form a triangular lattice.This triangular lattice is shown in Fig. 2 as green and red sites, corresponding to the w n=1,l and w n=2,l Wannier functions, respectively.In (47) we denote with ⟨n 1 l 1 , n 2 l 2 ⟩ a pair of nearest-neighbor Wannier functions, that is, two Wannier functions that overlap on exactly two rim sites.The operator B+ ⟨n1l1,n2l2⟩ creates a pair of particles that are delocalized on these two Wannier functions and whose total spin (the eigenvalue of ( Ŝn1l1 + Ŝn2l2 ) 2 with Ŝnl = ( Ŝx nl , Ŝy nl , Ŝz nl ) T ) is zero, i.e., they form a singlet state.Indeed, this is evident from (47), since the state B+ ⟨n1l1,n2l2⟩ |∅⟩ is symmetric under the exchange of the orbital degree of freedom 1 ↔ 2 and antisymmetric under the exchange of the spin ↑↔↓.
Using the operators introduced above, we can now provide the projected Hamiltonian of the dice lattice.For convenience we break it down into three parts The Hamiltonian Ĥtri.describes the hopping of the onsite pairs (singlets) on the triangular lattice formed by the Wannier function centers, namely the sublattice formed by the hub sites of the dice lattice, and the spin exchange interaction between localized particles on neighboring Wannier functions.It is given by The first term, proportional to the parameter A = 6 + U h ε 4 − /U r , gives the binding energy of an on-site pair.The energy scale used for the projected Hamiltonian is U r c 4 = 1 and A is the only free parameter.The terms in the second and third lines describe the nearestneighbor interaction and the hopping of on-site pairs on the triangular lattice.Taken together, they can be written as an isotropic Heisenberg exchange term Bnl • Bn ′ l ′ for the pseudospin associated to the on-site pair operators Bnl = ( Bx nl , By nl , Bz nl ) T with additional terms.Finally, the last line is the antiferromagnetic Heisenberg exchange interaction of on-site spins, which can be written as Ŝnl • Ŝn ′ l ′ .
The Hamiltonian Ĥtri.takes the same form as the projected Hamiltonian of the Creutz ladder [24] with the only difference being that the Wannier functions in the latter case form a one-dimensional chain instead of a triangular lattice.As in the case of the Creutz ladder, Ĥtri.possesses an extensive number of local integrals of motions given by the operators Ŝ2 nl , which have eigenvalue 3/4 if exactly one particle is present in the Wannier function labeled by nl and zero otherwise.This means that particles can move in the lattice only if they form on-site pairs and are localized otherwise.However, these integrals of motion do not commute with the full projected Hamiltonian due to the additional terms in Ĥkag.+ Ĥtri.−kag., as we will see below.Note that, due to the identity Ŝ2 nl + B2 nl = 3 4 1, the operators B2 nl do not form a set of independently conserved quantities.
The second term Ĥkag. in the projected Hamiltonian describes the hopping of bond singlets and takes the form + cyclic permutations of (1, 2, 3) . ( Here, we have used the abbreviations 1 ≡ n 1 l 1 , 2 ≡ n 2 l 2 and 3 ≡ n 3 l 3 and we denote by ⟨n 1 l 1 , n 2 l 2 , n 3 l 3 ⟩ a triplet of Wannier functions that are two by two nearestneighbors.Therefore their centers form an equilateral triangle whose sides are given by the three vectors v i=1,2,3 in ( 27)- (28).The sum ⟨1,2,3⟩ runs over all such triangles in the triangular lattice.The quantity s(n 1 l 1 |n 2 l 2 , n 3 l 3 ) = ±1 depends on the overlap of the Wannier functions forming a triangle of nearestneighbors and is defined as Note that we have used the fact that the Wannier functions of the dice lattice are real.The quantity s(n 1 l 1 |n 2 l 2 , n 3 l 3 ) is visualized in Fig. 2 as the color of the bonds connecting the sites in a kagome lattice, which is the lattice on which the bond singlets live.Due to the signs of the bond singlet hoppings, the Hamiltonian Ĥkag.
does not have the symmetry of the triangular lattice as in the case of Ĥtri., but instead has the same translational symmetry of the original dice lattice given by the fundamental vectors a i=1,2 (29).The last term in (48) describes the processes by which on-site pairs are converted into bond singlets and viceversa, Due to the sign factor s(1|2, 3), the translational symmetry of this term is the same as Ĥkag. .By allowing the motions of bond singlets, in which two neighboring Wannier functions are occupied each by a single particle, and the conversion of on-site pairs to bond singlets, the projected Hamiltonian of the dice lattice does not possess obvious local integrals of motion as in the case of Ĥtri.discussed above.Our main goal is to investigate whether some sort of conserved quantities associated to localized particles are nevertheless present in the projected Hamiltonian.
An important first step is the solution of the twoparticle problem, which is worked out in Appendix A. For three or more particles, one has to resort to numerical methods.In the next section, we analyse the energy spectrum and the nonequilibrium dynamics of the projected Hamiltonian in the case of three and four particles.

IV. EXACT DIAGONALIZATION RESULTS FOR THE DICE LATTICE
In some previous works, [24] it has been established that unpaired particles are always localized due to the presence of local conserved quantities in specific onedimensional models.These are the same conserved quantities of the Hamiltonian term Ĥtri.(49), as discussed in Section III B. However, they are spoiled by the additional terms Ĥkag.+ Ĥtri.−kag.associated with the bond singlets, which are not present in the one-dimensional case.Therefore, in the case of the dice lattice, it is not clear if local conserved quantities exist.In this section, we address this important question using exact diagonalization.
In Section IV A, we provide numerical results for the energy spectrum of the projected Hamiltonian of the dice lattice in case of an imbalance between spin up and down particles (N ↑ ̸ = N ↓ ).It is shown that the excited states form well-separated and almost degenerate groups of states.This is similar to, but at the same time different from, the perfect degeneracy of the exact eigenstates derived analytically in Section II B, which are found to form the ground state manifold in the case of spin imbalance.It is argued that the lack of degeneracy in the excited states is a finite-size effect.
In Section IV B, we focus on the time evolution from an initial state that includes both on-site pairs and single unpaired particles.It is shown that, in the long time limit, memory of the initial positions of the unpaired particles persists in both the particle density and the spin density.This is a clear signature of nonergodic behavior and many-body localization, which is a consequence of the presence of the local conserved quantities in manybody quantum systems.
All of the numerical results shown here have been obtained using the exact diagonalization package QuSpin.[79,80] Exact diagonalization is performed on a rectangular cluster with periodic boundary conditions composed of N c = N 1 N 2 unit cells, where N 1 and N 2 are the numbers of unit cells in the horizontal and vertical directions, respectively.For instance, a finite cluster with N 1 = 4, N 2 = 2 is shown in Fig. 2.

A. Energy spectrum
In Fig. 3, the lowest eigenvalues of the projected dice lattice Hamiltonian with three particles (N ↑ = 2, N ↓ = 1) are shown for various system sizes.The lowest 2Nc N ↑ −N ↓ = 2N c eigenstates are exactly degenerate and take the form These are precisely the exact eigenstates with localized quasiparticles presented in Sec.II B, see (25).In the balanced case (N ↑ = N ↓ ), it is known [23,25] that the states of the form ( 22) have minimal energy and thus are ground states of the projected Hamiltonian for given number of particles.No analogous statement exists for the imbalanced case, nevertheless one can reasonably expect states of the form ( 25) and (53) to have minimal energy.Our numerical simulations show that this is indeed the case, at least for some values of the parameter A in the projected Hamiltonian (48).In Fig. 3, the most interesting information is contained in the excited state spectrum for which there are no available analytical results.The excited states can be grouped into sets of states with approximately the same energy.Distinct sets are indicated with different colors in Fig. 3.The number of states in these quasidegenerate sets is exactly equal to the degeneracy of the ground state.On  (25).In this specific case, they take the form |Np = 1, On the other hand, the excited states can be grouped in sets of states with approximately the same energy, which have been denoted with different colors.The number of states in these sets of quasidegenerate states is the same as the degeneracy of the ground state.By increasing the system size, the deviation from perfect degeneracy in each set seems to decrease.However, even for the largest system sizes there is no perfect degeneracy in each set of excited states, as shown in the inset in the lower right panel.eigenstates are perfectly degenerate and have an energy E0 = −Ep = −22 for A = 10.These are the exact eigenstates in (25), which for the given number of particles take the form |Np = 1, On the other hand, the excited states can be grouped in sets of states separated by energy gaps, which have been denoted with different colors.Both the number of well-separated sets of excited states and the energy gaps increase with the system size.One such set is visible for N1 = 6, N2 = 3, which is colored in blue.Two such sets are visible for N1 = 8, N2 = 4.The number of states in these sets is exactly four times larger than the degeneracy of the ground state.
increasing the system size, the deviation from perfect degeneracy in each set of states seems to decrease.However, it is noted that even for the largest computed system size, there is no perfect degeneracy in each set of excited states.For example, for system size N 1 = 14, N 2 = 7 the energy width of one of these set of states is roughly 2 × 10 −3 (see inset in Fig. 3).
This peculiar structure in the excited state spectrum is not entirely surprising and can be understood by comparing with what one would obtain in the case of the one-dimensional models in which local integrals of motion are known to exist.[24] In these models, the eigenvalues s nl (s nl + 1) of the operators Ŝ2 nl , that is the total spin on the Wannier function labeled by nl, are good quantum numbers.The lowest lying states have s nl = 1/2 only for a single Wannier function and s nl = 0 otherwise since there is only a single unpaired particle.For any given state of this form, one can obtain a distinct eigenstate with the same energy by applying a translation, due to the translational invariance of the projected Hamiltonian.The new eigenstate is distinct since the localized unpaired particle present in the original state has been moved to a different Wannier function (s nl = 1/2 → s n,l+j = 1/2 for a translation by j unit cells).Thus the lowest lying states have degeneracy which is at least equal to the number of unit cells N c .
The numerical results shown in Fig. 3 suggest that local integrals of motion are present in some approximate sense for the projected Hamiltonian of the dice lattice.The mechanism behind the lifting of the degeneracy is at the present not well understood, except for the fact that it is caused by the motion of the bond singlets and their interaction with localized single particles.As mentioned above, perfect degeneracy is approached in each set of excited states on increasing the system size.Thus, it might by the case that the projected Hamiltonian of the dice lattice possesses exact local integrals of motion only in the limit of infinite system size.We conjecture that these encode the localization of quasiparticles in wave functions that extend over rather large distances.In fact their range, or localization length, should be comparable to the largest size shown in Fig. 3 in order to explain the finite size effects.On the other hand, in the onedimensional models of Ref. 24, the quasiparticle wave functions are compactly localized since they coincide with the Wannier functions.
It is interesting to explore what happens with increasing imbalance.In Fig. 4, for instance, we show the energy spectrum for N ↑ − N ↓ = 2.The ground state manifold is composed of states of the form

) and has degeneracy 2Nc
N ↑ −N ↓ = N c (2N c − 1).In the spectra shown in Fig. 4, the excited states are separated from the ground state manifold by a large energy gap.With increasing system size, sets of excited states well-separated from each other by energy gaps start to appear.Both their number and the energy gaps that separates them increase with the system size.In the figure, one such set is visible for N 1 = 6, N 2 = 3, but two such sets are evident for N 1 = 8, N 2 = 4.In each case, the number of eigenstates in these sets is exactly four times larger than the degeneracy of the ground state manifold.As in the three particle case, the energy width of the well-separated sets of states reduces with increasing system size and it is possible that they approach perfect degeneracy for an infinitely large system.Alternatively, they could break up in several degenerate subsets.
In summary, we observe for N ↑ − N ↓ ̸ = 0 that the excited states form almost degenerate multiplets with a degeneracy that is an integer multiple of the ground state degeneracy.We attribute the deviation from perfect degeneracy to finite-size effects that lead to the breaking of the local integrals of motion present in an infinite lattice.Some rather strong finite-size effects also appear in the long time dynamics, as shown in the next section.

B. Nonequilibrium dynamics
The time evolution of the particle and spin densities obtained by evolving a three-particle initial state of the form with the projected Hamiltonian H int is shown in Figs. 5, 6.This initial state can be equivalently described as either an on-site singlet present at site A with an additional spin-up particle at site B, or as a bond singlet between sites A and B plus a spin-up particle at site A, namely These representations give us insight into the expected behavior of the time-evolved system.During the time evolution, the bond singlet can move away from sites A and B, leaving behind the particle with spin up at site A. Thus, the spin density at site A becomes nonzero (in the initial state the spin density is nonzero only on site B).This differs from when only two particles are present in the system, in which case the spin density remains at zero at all sites where it is initialized to zero.Therefore, no spin transport occurs with two particles, whereas some form of spin transport takes place with three particles, as seen in Fig. 5.
In Fig. 5, the time averaged particle and spin densities are shown for different system sizes.Here, the time average of an observable (computed numerically) is defined as The first panel shows the time averaged particle and spin densities at each site averaged from time t0 to t1 = 10000 according to (57).t0 = 1000 is selected for the first three system sizes while t0 = 3000 is chosen for N1 × N2 = 10 × 5.The error bars indicate the standard deviation of each quantity over the same time interval, as given by (58).The panels in the second column show the initial state of the system |Ψ(0 The diameter of the circle represents the relative particle density at a site and the spin density is represented by a color scale, with red being up spin and blue being down spin.The panels in the third column depict the time averaged particle and spin densities on the lattice.The time averaged particle and spin densities are normalized to one at site B. The spin density is nonzero across all sites of the lattice, but remains largest at site B, where the particle with spin up was initially located.FIG. 6.Time evolution with three particles for system size N1 = 10, N2 = 5.The initial state is shown in the lowest row of Fig. 5.The upper panel shows the particle density from t = 0 until t0 = 3000 for three sites denoted as A, B and C in Fig. 5.The inset shows the rolling average of the particle density in the same sites from t0 = 3000 until t = 7500.The lower panel depicts the spin densities of sites A, B and C from t = 0 until t0 = 3000, with the inset depicting the rolling average of the spin densities from t0 = 3000 until t = 7500.The rolling average in the insets are taken with a time window of ∆t = 3.0, to smooth out fluctuations.The rolling averages are used for the sake of clarity in presentation. where ⟨ Ô(t)⟩ is the expectation value of an observable Ô at time t, t 0 is an initial cutoff time, which is sufficiently large to ensure that the observable is close to its longtime asymptotic value and t 1 is the end time.In the first column of Fig. 5, the standard deviation of the time average of an observable is also shown, and is defined as In Fig. 6, it can be seen that, for a cluster of size N 1 = 10 N 2 = 5, the particle and spin densities approach their respective equilibrium value at around t = 3000, which is chosen as initial cutoff time t 0 .The time scale to reach equilibrium is extremely long and increases with the system size.For instance, we take t 0 = 1000 for the other system sizes shown in Fig. 5. From Fig. 5, it can be observed that the particle and spin densities diffuse throughout the system and becomes nonzero on all lattice sites.However, at site B where the particle with spin up is initially placed, the particle and spin densities are significantly larger than those at the other sites in the lattice, including site A, the initial position of the on-site pair.This behaviour persists as the system size is increased, as seen in Fig. 5.As expected, with an increase in system size, the fluctuations of the particle and spin densities in all lattice sites also decrease, with the exception of site B. A surprising result for system size N 1 = 8, N 2 = 4 is the appearance of three other sites with particle and spin densities comparable to the ones of site B. These three sites are obtained by translating site B by linear combinations of the lattice vectors R 1 /2 and R 2 /2, with R 1 = N 1 a 1 and R 2 = N 2 a 2 the vectors defining the size and shape of the finite cluster with periodic boundary conditions.A similar behavior is observed also for system size N 1 = 4, N 2 = 2, but not in the case of the other system sizes shown in Fig. 5. Thus, the long-time dynamics appears to be very sensitive to finite size effects.The origin of this phenomenon is unclear at present and deserves further investigations.
Analogous results are obtained in the case of four particles N ↑ + N ↓ = 4.The two unpaired particles in the initial state can have either the opposite or same spin.Numerical results for both initial states are shown in Fig. 7. Again, memory of the initial positions of the two unpaired particles persists in the particle and spin densities, which are otherwise approximately uniform over all remaining sites in the long time limit.
The results of Figs.5-7 indicate a peculiar coexistence of ergodic and nonergodic behavior in the same system.On one hand, we have that in the long time limit the spin density distribution is highly nonuniform, in particular it peaks at site B. Thus, memory of the initial position of the particle with spin up is retained during the time evolution.If the time evolution were to be completely ergodic, we would expect the long-time average of the system to coincide with the average value given by the thermal equilibrium state that is the usual Gibbs ensemble.This result is known as the ergodic theorem and is the fundamental postulate of statistical mechanics [81].The equation for the long-time average of an observable stated in (57) is its usual definition in the context of the ergodic theorem.The results of Figs.5-7 do not coincide with the expected thermal equilibrium average, which for a translationally invariant system as studied here, would see the particle density and spin density be uniformly distributed.We can conclude from this that the states do not thermalize and the time evolution displays nonergodic behaviour, which we interpret as a signature of in which the unpaired particles have both spin up.As in Fig. 5, the particle and spin densities retain memory of the initial positions of the unpaired particles (site B and site C) in the long time limit.On the other hand, the particle and spin densities become almost uniform throughout the remaining sites.A difference compared to the three particle case is that the time scale required to reach the long time limit is about ten times smaller and is given approximately by t0.
many-body localization [50,82,83].This characteristic signature of many-body localization has been observed in ultracold gas experiments [51][52][53].On the other hand, any information regarding the initial position (site A) of the on-site pair is completely lost after a sufficiently long time.Therefore, the nonergodic behavior is limited to fermionic quasiparticles, while two-body bound states behave ergodically, as one would expect.
Our results also suggest a violation of the eigenstate thermalization hypothesis (ETH), which asserts that the thermalization of all initial states of a system upon time evolution implies that all many-body eigenstates are thermal as well [84].The initial state (55) used to study the nonequilibrium dynamics is a linear combination of highly excited eigenstates of the system, and the absence of thermalization of the initial state indicates that at least some of the many-body eigenstates are non-thermal as well.The violation of ETH is commonly associated with many-body localization [50,82].
Finally, it is interesting to note that the time scale for the equilibration of the particle and spin densities, which is very long as shown in Fig. 6, seems to be related to the energy scale associated to the lifting of the degener-acy in the spectra shown in Fig. 3 and discussed in Section IV A. The time scale is approximately given by the initial cutoff t 0 used for the time average (57).Indeed, we take t 0 = 1000 for N 1 = 8, N 2 = 4 and t 0 = 3000 for N 1 = 10, N 2 = 5.Correspondingly, the energy width of the quasidegenerate set of states is roughly three times larger for the smaller system size.For instance, in the case of the sets colored in green in Fig. 3, the energy widths are ∆E = 4.8 × 10 −3 and ∆E = 1.9 × 10 −3 , respectively.The orders or magnitude of t 0 and 1/∆E are also in good agreement.Thus, a possible explanation is that the long equilibration times are a manifestation of the local integrals of motion, which are slightly broken by finite size effects.

C. Finite size scaling
In Sec IV B, we calculate the time evolution of the particle density and spin density over long times and observe signatures of many-body localization.However, the phenomena of many-body localization has been debated to be simply a finite size effect [85] and that at infinite sys-tem size, many-body localization would cease to exist.To understand the impact of system size, and as a consequence, finite-size effects on the signatures observed in our model, we perform a scaling analysis.In Fig. 8, we can see the how the asymptotic value of the particle densities and spin densities scale with system size.By observing the trends of these quantities as they approach an infinite system size (N −1 c → 0), we can estimate the impact of the finite size effects.We can see that the particle and spin densities at site A, where the on-site pair was placed, approaches zero as the lattice size increases.At site B, where there is initially a single unpaired spin, the particle and spin densities approach a finite nonzero value as the system size increases.The long-time average of the observables reveal a continued partial breaking of ergodicity.The asymptotic distribution of the spin density is not uniform and peaks on the site in which the unpaired particle was initially located.This is in contrast with the distribution of the spin density at thermal equilibrium, which is uniform as a consequence of translational invariance.The finite-size scaling results indicate that the partial breaking of ergodicity, and memory of the particle and spin density at site B will both remain in an infinite system, supporting the case for true many-body localization in an infinite system.It would be important to perform a finite size scaling analysis by also increasing the particle number, but this is currently beyond reach for exact diagonalization since even for four particles, the system sizes for which we can perform the time evolution is very restricted.

V. CONCLUSION AND DISCUSSION
In this work, we have provided strong evidence through analytic arguments and numerical results obtained with exact diagonalization that propagating bosonic twobody bound states can coexist with the localization of fermionic quasiparticles in lattice models with flat bands.The analytical results, namely exact eigenstates of the projected Hamiltonian of the form (25), are valid for generic lattice models with flat bands under a few assumptions (spin rotational and time-reversal symmetries, uniform pairing condition) and are a simple consequence of the SU(2) symmetry that emerges in the isolated flat band limit.[23] We have also shown how these analytical results can be extended in the case of the dice lattice due to the compact nature of the Wannier functions.
In a previous work, [24] local integrals of motion have been explicitly constructed for a few one-dimensional lattice models with flat bands.As a consequence, fermionic quasiparticles are strictly localized in these models and spin transport is completely suppressed.Similar analytical results are not available for dimension two or higher, therefore we have focused our numerical investigations on the two-dimensional dice lattice.The presence of local integrals of motion in the one-dimensional models of Ref. 24 implies that all the eigenstates are degenerate in the case of spin imbalance, with a degeneracy no less than the number of unit cells N c of the finite cluster.In the case of the two-dimensional dice lattice, we have observed from the energy spectrum computed numerically, that the same degeneracy is approximately present in the excited states.This is in contrast with the ground state multiplet, which is composed of states of the form (25) and is thus perfectly degenerate, see (26).With increasing system size the degeneracy of the excited states seems to be restored.We interpret this finding as a signature of the presence of local integrals of motion in the projected Hamiltonian of the dice lattice, with the slight breaking of the degeneracy in the excited states caused by finite size effects.Further studies are needed to better understand the nature of these finite size effects.
Another clear signature of many-body localization is provided by the time evolution of the particle density and the spin density starting from initial states comprising one on-site pair and one or two unpaired particles.The long time average of the observables reveals a partial breaking of ergodicity.For instance, the asymptotic distribution of the spin density is not uniform and is peaked on the site in which the unpaired particle (or particles) was initially located.This is in stark contrast with the distribution of the spin density at thermal equilibrium, which is expected to be uniform as a consequence of translational invariance.The persisting memory of the initial state in local observables is a characteristic feature of many-body localization and has been observed in experiments with ultracold gases.[52,53] Another interesting observation is that, in the case of the one-dimensional models of Ref. 24, spin transport is strictly forbidden as a consequence of the local conserved quantities, while this is not the case for the dice lattice, since, for initial states with N ↑ + N ↓ ≥ 3, the spin density can become finite at later times even on a Wannier function in which it was initially zero.A possible explanation is that the localized states corresponding to the local conserved quantities have an extension comparable to the size of the system used in the exact diagonalization.
Besides providing interesting information from a theoretical standpoint, our numerical protocol suggests a possible experiment that could be performed using ultracold gases in optical lattices.This might be feasible in the near future since a practical scheme to implement the dice lattice with a finite magnetic flux has been suggested [32] and spin-resolved imaging at the single atom level is nowadays routinely performed.Ultracold gases have the potential to simulate the dice lattice Hamiltonian for much larger system sizes and particle numbers than can be done with exact diagonalization or other currently available numerical methods.
Our work can be extended in several possible directions.For instance, it would be interesting to consider other lattice models to confirm that the numerical results obtained in the case of the dice lattice are generic.Studying other models is probably more difficult numerically than the dice lattice because the Wannier functions are, in general, not compactly localized.This means that the projected Hamiltonian would contain many more terms than in the case of the dice lattice and would be less sparse.Thus, exact diagonalization becomes numerically heavier.In this respect, it is interesting to note that strong finite size effects are still present even for the largest system sizes that we are able to reach at present with exact diagonalization.Thus, improving the efficiency of numerical simulations is of paramount importance.The finite size scaling analysis presented in Sec.IV C provides some evidence that the signature of many-body localization that we observe persists in the limit of infinite size, however it would be important to perform the same analysis with a larger particle number.This is currently out of reach for the numerical methods at our disposal.
From our point of view, the key open question is how to rigorously prove the existence of local conserved quantities, or better yet, how to explicitly construct them.In that sense, it would be interesting to borrow methods from the field of many-body localization, in which local conserved quantities have been investigated extensively, [47,[54][55][56][57] and adapt them to our case.The existence of an extensive number of local conserved quantities is generally considered the defining property of many-body localization.However, even in the most extensively studied models, such as one-dimensional spin chains, there is no consensus on whether a many-body localization transition takes place, precisely because rigorously proving the existence of local integrals of motion in a many-body system is a highly nontrivial task.In-deed, in recent works it has been suggested [86,87] that the critical disorder strength of the many-body localization transition increases with the system size, meaning no transition takes place in the thermodynamic limit.These claims were subsequently criticised [88][89][90] on the basis that it is difficult to extrapolate from the information obtained using small systems.The delicate issue of the existence of many-body localization in the thermodynamic limit is still actively debated and no consensus has been reached yet [85,91].
We believe that lattice models with flat bands provide an interesting platform to tackle the difficult problem of many-body localization and that of its fate in the thermodynamic limit.On one hand, similar to disordered systems, we observe pronounced finite size effects, however, localization is realized in a translationally invariant setting which allows us to derive exact analytical results.A particularly important exact result in the case of the dice lattice is that one term of the projected Hamiltonian possesses an extensive number of local integrals of motion, namely [ Ĥtri. , Ŝ2 nl ] = 0.Moreover, the numerical results suggest that these integrals of motion are not destroyed, but simply deformed by the other terms of the Hamiltonian Ĥkag.+ Ĥtri.−kag., and in fact play an important role in determining the system dynamics.Our hope is that the lessons learned in the future in the relatively controlled setting of translationally invariant lattice models with flat bands can be subsequently transferred to disordered systems.Our work is an initial step in this promising research direction.
We conclude by noting that the idea of localized superconductors, that is, superconductors with localized quasiparticle excitations due to strong enough disorder, has already been proposed long ago [92,93].In this work, we provide evidence that localization of quasiparticles can occur also in two-dimensional lattice models with flat bands even in the absence of disorder.Strictly speaking, the localized quasiparticles are themselves a source disorder [64] and an interesting question is how they affect the supercurrent flow.We believe that the coexistence of antithetic phenomena, such as superfluidity and localization, in the same system is a rather interesting occurrence that deserves further attention.Moreover, superconductors with suppressed quasiparticle transport have a number of interesting technological applications, see for instance, Ref. that lives on the two Wannier functions centered at positions r l1 and r l2 in the triangular lattice.For the purpose of solving the two-body problem, the phase factor χ introduced above is absorbed in the definition of the new bond singlet operators B+ ⟨l1,l2⟩ .Then, we can construct plane wave linear combinations The dispersion of the two-body bound states obtained by the diagonalization of H (2) (k) is shown in Fig. 10.A notable feature is the presence of a flat band at zero energy.The flat band corresponds to the following eigenvector of the two-body Hamiltonian |G f.b. (k)⟩ = 1 A9) which is a linear combination of the bond singlets only, that is, of the states |k, l = 2, 3, 4⟩ in (A7) but not of the on-site singlet states (|k, 1⟩).Therefore, the state (A9) is an eigenstate also of the Hamiltonian that describes the hopping of a bond singlet between nearest-neighbor sites of the kagome lattice, that is the lower 3 × 3 diagonal block in (A8).The dispersion obtained by diagonalizing only this block is shown in Fig. 10 as well.The structure of the states of the kagome lattice flat band has been studied in detail.[95] As in the case of the dice lattice, the flat band subspace is spanned by compact wave func-tions, but with the essential difference that they are not orthogonal since they overlap at most on a single site of the kagome lattice.The presence of the flat band of twobody states in (A9) is not unexpected since it is a consequence of the fact that the flat band two-body Hamiltonian has a rank which is not larger than the number of orbitals per unit cell N orb in the original lattice.
For k = 0 the lower 3 × 3 diagonal block in the twobody Hamiltonian, corresponding to the bond singlets states |k, l = 2, 3, 4⟩, completely decouples from the onsite singlet state |k, 1⟩.This is in agreement with the general exact result presented in Sec.22, according to which the two-body state |N p = 1⟩ = b † |∅⟩ is an exact eigenstate representing a Cooper pair with zero momentum.It is only in the dice lattice that we have the freedom of tuning the energy of this state by varying the parameter A. This is a consequence of the fact that the uniform pairing condition can be relaxed in the case of the dice lattice, see Sec.III A. For k ̸ = 0, the two-body bound states are, in general, linear combinations of onsite and bond singlets, with the exception of the flat band FIG. 1.Schematic of the dice lattice.Six-fold coordinated sites (hub sites, orbital index α = 1, 2) are denoted by hexagons, while three-fold coordinated sites (rim sites) are denoted by triangles.The rim sites can be further divided into two triangular sublattices denoted by up-(α = 3, 5) and down-pointing triangles (α = 4, 6), respectively.The bonds denote nearest-neighbour hoppings which are all real and have the same absolute value t = 1.The color of the bond denotes the sign of the hopping amplitude, grey for +t and red for −t.The rectangular box is the magnetic unit cell, given by the fundamental vectors aj=1,2(29).Two Wannier functions of the two lowest flat bands of the dice lattice Hamiltonian are also shown.The Wannier function w nl , with n = 1, 2, is centered on the hub site α = n in unit cell l and is nonzero only on the same hub site and the adjacent rim sites, therefore it is compactly localized.These "flower states"form the orthonormal basis used in the expansion of the projected Hamiltonian(7).

FIG. 3 .
FIG.3.Lowest energy eigenvalues of the projected dice lattice Hamiltonian Hint (48) obtained from exact diagonalization.Different systems sizes are shown, while the number of particles is fixed to two particles with spin up (N ↑ = 2) and one particle with spin down (N ↓ = 1).The system size is determined by a pair of integers, the number of unit cells in the horizontal N1 and vertical N2 directions, respectively.Only the lowest 16Nc eigenvalues, with Nc = N1N2 the number of unit cells, are shown in each panel.The lowest 2Nc eigenstates are perfectly degenerate and have energy E0 = −Ep = −22 (the only free parameter in the projected Hamiltonian has been fixed to A = 10).These are the exact eigenstates with localized quasiparticles presented in Section II B, see(25).In this specific case, they take the form |Np = 1,I ↑ = {(n, l)}⟩ = b † d † nl↑ |∅⟩.On the other hand, the excited states can be grouped in sets of states with approximately the same energy, which have been denoted with different colors.The number of states in these sets of quasidegenerate states is the same as the degeneracy of the ground state.By increasing the system size, the deviation from perfect degeneracy in each set seems to decrease.However, even for the largest system sizes there is no perfect degeneracy in each set of excited states, as shown in the inset in the lower right panel.

FIG. 4 .
FIG. 4. Same as in Fig. 3 in the case of three particles with spin up N ↑ = 3 and one particle with spin down N ↓ = 1.Three different system sizes are shown: Nc = N1N2 = 8, 18, 32.The lowest 13 2Nc 2 eigenvalues are shown in each panel.The lowest

FIG. 5 .
FIG. 5. Time evolution of the projected dice lattice Hamiltonian for three particles.Different systems sizes are shown in each row (from top to bottom: N1 × N2 = 4 × 2, 6 × 3, 8 × 4 and 10 × 5), while the number of particles is fixed to N ↑ = 2 N ↓ = 1.The first panel shows the time averaged particle and spin densities at each site averaged from time t0 to t1 = 10000 according to(57).t0 = 1000 is selected for the first three system sizes while t0 = 3000 is chosen for N1 × N2 = 10 × 5.The error bars indicate the standard deviation of each quantity over the same time interval, as given by(58).The panels in the second column show the initial state of the system |Ψ(0)⟩ = d †

FIG. 7 .
FIG. 7. Same as in Fig. 5 but for four particles N ↑ + N ↓ = 4 and system size N1 = 6, N2 = 3.The time interval for the time averaged particle and spin densities is taken to be between t0 = 100 and t1 = 5000.The initial state for the top row is |Ψ(0)⟩ = d † and includes two unpaired particles with opposite spin in sites B and C. The initial state for the bottom row is |Ψ(0

FIG. 8 . 1 c
FIG. 8. Finite size scaling of the long time average of the particle and spin densities.The long time average of the particle and spin densities on sites A and B are shown as a function of the inverse number of unit cells N −1 c in the finite cluster.The numerical data shown here are the same as in Fig. 5 and are obtained from the time evolution of the projected Hamiltonian with three particles (N ↑ = 2, N ↓ = 1).

FIG. 10 .
FIG.10.The dispersions En(k) of two-body bound states obtained from the diagonalization of (A8) are shown as blue lines for A = 10.The Brillouin zone corresponding to the triangular lattice(27) and its high-symmetry points are shown in the inset.A notable feature is the presence of a flat band at zero energy, whose corresponding eigenvector is given in (A9).The dashed red lines are dispersions obtained by diagonalizing the Hamiltonian that describes the hopping of the bond singlets between the sites of the kagome lattice, that is the lower 3 × 3 diagonal block in (A8).The green dashed line is the dispersion of the on-site singlets, if these are decoupled from the bond singlets.This dispersion is given by the top diagonal element of the two-body Hamiltonian [H(2) (k)]1,1 = −A − 4 3 i=1 cos(k • vi).