Flat-band many-body localization and ergodicity breaking in the Creutz ladder

We study disorder-free many-body localization in the flat-band Creutz ladder, which was recently realized in cold-atoms in an optical lattice. In a non-interacting case, the flat-band structure of the system leads to a Wannier wavefunction localized on four adjacent lattice sites. In the flat-band regime both with and without interactions, the level spacing analysis exhibits Poisson-like distribution indicating the existence of disorder-free localization. Calculations of the inverse participation ratio support this observation. Interestingly, this type of localization is robust to weak disorders, whereas for strong disorders, the system exhibits a crossover into the conventional disorder-induced many-body localizated phase. Physical picture of this crossover is investigated in detail. We also observe non-ergodic dynamics in the flat-band regime without disorder. The memory of an initial density wave pattern is preserved for long times.


Introduction
Localization in non-interacting electron systems has been extensively studied since Anderson discussed the disorder effect on the single-particle electron wavefunction in solids [1]. Presently, what is called Anderson localization (AL), is recognized as a universal phenomenon in various physical systems [2]. In AL quantum system, a single-particle electron wavefunction is exponentially localized with a finite localization length, and an insulating phase forms. Owing to the recent development in the computational power and numerical techniques, study on the effect of the interactions between particles on AL is currently one of the main research topics in condensed matter physics. It is now recognized that AL persists in some cases even if the particles interact. This is called many-body localization (MBL). Mostly by numerical simulations, it has been clarified that the MBL phase exhibits some characteristic properties such as Poisson distribution in the level spacing analysis (LSA) of the energy eigenvalues similar to that of the conventional AL and the logarithmic growth of entanglement entropy. In its glassy dynamics, MBL is closely related with the breaking of eigenstate thermalization hypothesis and ergodicity breaking dynamics [3][4][5][6][7][8]. This means that a closed ergodicity-breaking system does not thermalize for a long time, and if we prepare a non-entangled initial state in such a system, the information of the initial state is conserved for a long time without being lost. Recent experiments on cold-atom gases in optical lattices have reported evidences for the existence of MBL phenomena [9][10][11][12].
Motivated by the above findings, we shall report another type of disorder-free MBL system in this paper. It is a flat-band system with interactions. Certain flat-band structure suppresses particle hoppings effectively and generates a localized Wannier state [30,31] that is similar to the localized states in the conventional AL system. Such a localized Wannier state was theoretically investigated for certain non-interacting flat-band systems with and without weak disorders [32,33]. We are motivated by the existence of such localized wavefunctions and study a flat-band type localization in the Creutz ladder [34]. The Creutz ladder is a simple model and also experimentally feasible in cold atom gases. So far, there are several theoretical proposals for implementation of the model [35][36][37][38], and cold atom experiments realized some related systems [39][40][41], whose the physical properties have been extensively studied [42][43][44].
This paper is organized as follows. In section 2, we introduce the target Creutz ladder model. We focus on the flat-band case, and analytically study properties of the flat-band states. We explicitly reveal the origin of localization and discuss the possibility of MBL with repulsions. Effects of on-site disorders are also discussed, and the global phase structure is given.
In section 3, we present results of the numerical study. We first perform the LSA and also the level spacing ratio (LSR) analysis for the system with weak disorders under the flat-band condition and find that the probability distribution exhibits Poissonian behavior for both the non-interacting and interacting cases, indicating a localization tendency. Interestingly enough as the strength of the disorder is increased, we find that both the LSA and LSR exhibit behavior of Gaussian unitary ensemble (GUE) corresponding to extended (delocalized) states. These results are compared with those of the non-flat case in order to clarify the difference between the flat-band and non-flat-band cases. The above phenomenon is discussed via the analytical study in section 2. Then, we investigate the inverse participation ratio (IPR) to find that its results corroborate the localization tendency of the flat-band Creutz model. In particular, energy-resolved IPR exhibits very interesting behavior, which explicitly clarifies typical properties of the flat-band states as increasing the strength of disorder. We finally investigate distribution of the localization length for typical disorder strengths. Energy-resolved distribution reveals origin of the crossover observed by the LSA and IPR.
In the final subsection of section 3, we study the dynamics in the flat-band Creutz ladder, i.e. we investigate the time-evolution of states in which fermions are periodically put on sites. The result shows ergodicity-breaking dynamics, i.e. the memory of the particle distribution in initial states is preserved for long times. Besides the above important result of the non-ergodicity of the Creutz ladder, we find another interesting phenomenon for the cases of 1/6 and 1/4-particle filling. section 4 is devoted for conclusion. We present the summary and also give future perspective.

Creutz ladder model and flat-band localization
In this work, we study an interacting Creutz ladder model with the Hamiltonian [34] [ V n n n n n n n n n n i h . c .
where ( †) a j and ( †) b j are the fermion annihilation (creation) operators on the upper and lower chains, respectively, and subscript j denotes a unit cell. ( ) n a b j , is the number operator of the particle on the upper (lower) chain. t 1 and t 0 are the intra-chain and inter-chain hopping amplitudes, respectively. V is the intra-chain and inter-chain repulsions, as depicted in figure 1(a), which is one of the simplest interactions suitable for the present study as we explain shortly. There are two possible ways to implement this type of interactions in real experiments: (I) Method to use electric or magnetic dipole-dipole interactions between atoms [45,46], (II) To use natural overlap of Wannier functions between neighboring sites connected by horizontal and diagonal links induces to this type of interaction. The case (I) may induce vertical interactions, but we ignore them in this work. We verified that the vertical interactions do not change the subsequent numerical results substantially. Obviously, the repulsive Vinteraction prefers the density-wave configurations in the ladder direction. ( ) m a b j , is a random disorder chemical potential, which has a uniform distribution, such as , and breaks the chiral symmetry 4 . This choice of the disorder plays a significant role in the localization phenomenon in the present model as we explain shortly.
The energy spectrum of the non-interacting case of H in equation (1), with V=μ=0 is given as where k is the wave number and the bandwidth is | ( )| t t 2 1 0 . As shown in figures 1(b) and (c), the band is flat for t 0 =t 1 with E(k)=m2t 0 , whereas it is dispersive for ¹ t t 0 1

5
. The noninteracting case of H in equation (1) with μ=0 belongs to the BDI class in the topological classification theory 4 The disorders act independently in each site and they effect destroying the localized Wannier state in equation (2). If one employs chiral symmetric disorders instead, the system is fairly robust against the disorders. 5 Such a hopping amplitude ratio can be easily realized in real experiments [41]. By controlling a modulation amplitude of a driving optical lattice, t 0 =t 1 condition can be achieved. [47][48][49][50]. Hence, the model has chiral, time-reversal and particle-hole symmetries. In particular, the chiral symmetry makes the energy spectrum symmetric around zero energy. In addition, at the flat-band point, t 0 =t 1 , a localized Wannier state exists in the system, whose wavefunction for the lower spectrum is given by [30,31,37] where | ñ 0 is the empty state. The state |Y ñ w j spans over two adjacent unit cells, i.e. it is a four-site localized state, and there are two |Y ñ w j ʼs per site. It is quite useful to study analytically the flat-band case of the present system for the forthcoming numerical investigation. In that case, the hopping part of the Hamiltonian reduces to the following one, H flat , Then, we introduce the following operators . This transformation is a kind of detangling for a lattice system [32]. Under this transformation, the Creutz ladder is detangled into a simple lattice system where each lattice site is completely decoupled each other. In terms of w Aj and w Bj , H flat is expressed as and straightforward manipulations show Equations in equation (6) reveal very important properties of the Creutz ladder mode with the flat-band coupling, i.e. in terms of {w A , w B }-'particles', w A(B) -particle hops only left (right)-hand site and changes to w B(A) -particle. Therefore, the {w A , w B }-particles strictly localize on two adjacent rungs of the ladder. It is obvious that the Wannier state in equation (2) is nothing but a static state composed of a pair of nearest-neighbor {w A , w B } such as Similarly, the upper-spectrum eigenstates can be constructed easily as ( Therefore, the flatband Hamiltonian, H flat , can be expressed in terms of the following operators, †  W j , that create energy eigenstates   One may wonder how the original fermion a j (b j ) behaves. Obviously, they do not create an eigenstate of the Hamiltonian. However, a j and b j are a simple superposition of w Aj and w Bj , i.e.
( ) The above state in equation (9) is obviously localized. The analytical study in the above gives the following important observations on the Creutz ladder model in equation (1): (i) in the clean and non-interacting flat-band case, the Creutz ladder system is strictly non-ergodic and all eigenstates are localized; (ii) the localization 'length' is four lattice sites. The Wannier state in equation (2) resides on four sites. In the state expressed by equation (9), a particle resides on a single site and six sites with equal probability. Such a localized particle can be regarded as a concrete example of a flat-band compactons. More general discussion and construction for the flat-band compactons have been given in [51,52]; (iii) under a disorder such as μ a,j =μ b, j , the w-particle picture is robust, i.e. no on-site mixing of the w A and w B -particles takes place, and therefore the above localization properties are intact. On the other hand, a disorder such as m m ¹ a j b j , , , which we employ in the present work, tends to break the w-particle picture as it induces an on-site mixing; (iv) similarly, the interaction term in equation (1) is expressed by the w-particle in the diagonal form and therefore, the w-particle picture is robust even in the presence of the interaction.
Before going into the practical calculations, we shall give some comments. (1) In the following section, we consider the 1/8-filling case. In such a low commensurate filling, particles described by equation (2) do not overlap substantially [53]. Then, it is expected that the w-particle picture is preserved even for rather strong V-interactions under weak disorder, and the system exhibits localization. This is nothing but a new kind of MBL. The conventional disorder-induced MBL needs sufficiently strong disorders [6]. On the other hand, our considering MBL is induced by the flat-band, i.e. distractive interference of hoppings. (2) In the ordinary AL systems, localization length depends on the disorder strength. On the other hand in the above MBL regime, the Wannier state in the flat-band has finite components in definite lattice sites. We note that this properties give certain suggestions on the set up of an initial state for observing MBL dynamics in simulations that we shall give in later section. (3) Increasing the disorder strength μ, the w-particle picture is getting unstable, and the genuine flat-band localization is expected to be destroyed. We expect that a crossover takes place from the flatband localized states to a new kind of states at a critical disorder strength, μ c .

Numerical studies
In this section, we shall study the Creutz model by the numerical methods. As a hallmark of localization and (non)ergodicity, we investigate the level spacing, the IPR and the temporal evolution of inhomogeneous states.
Obtained results all support the picture of the flat-band localization given in section 2. Furthermore, the numerical results show interesting behavior of the model, in particular at relatively high fillings, which come from the interplay between the locality of the flat-band regime and the repulsion. In what follows, we employ t 1 as a unit of energy.

Level spacing analysis
We first perform the LSA by full-diagonalization of the Hamiltonian H in equation (1), under the periodic boundary condition. In the LSA, we employ the usual unfolding analysis [54]. In the unfolding method [21], we first prepare a set of energy-eigenvalue spectrum {E i } (i=1, 2, L, N D ; N D is the Hilbert space dimension) in ascending order, and then calculate the average level spacing of the original spectrum By using ΔE, we define a new level spacing set {s i } as s i =(E i+1 −E i )/ΔE. From the set {s i }, we obtain the statistical distribution P(s), which is to be compared with the level statistics of the random matrix theory. When we use multiple realizations (samples) of the disorder, we average P(s) with respect to them to obtain the final result of P(s).
On performing the LSA for the disorder-free case (μ=0), it is important to note that the system has the translational symmetry. This symmetry generally leads to numerous degeneracies in the energy eigenvalues. Because of the degeneracies, it is not simple to obtain the probability distributions of the level spacing without ambiguities [14,15]. To avoid this difficulty, we consider the cases with small but finite disorders. In the presence of disorders, even those that are extremely weak, the degeneracies of the energy eigenvalues are solved. In practical calculations, we consider the upper and lower chains with length L=16 and number of particle N=4 6 . From the LSA, one can examine the localization properties of the system. In general, for an ensemble of localized states, the probability distribution exhibits Poisson statistics, such as ( ) ( ) µ -P s s exp P , where s denotes the unfolded level spacing. Contrastingly, for an ensemble of delocalized (extended) states, the probability distribution is to be GUE, with characteristics such as ( ) [55-59, 6, 7, 60, 61]. 7 Figure 2(a) shows the obtained probability distribution for various disorder strengths for the noninteracting flat-band (V = 0, t 1 =t 0 ). We find that for a weak disorder (μ=1), the probability distribution is extremely similar to Poisson statistics. This result indicates the existence of localized states even in a weak disorder. With increasing disorder strength, we observe an interesting phenomenon, i.e. first, the statistics changes from Poisson to GUE-like, and then it returns to the Poisson statistics. Calculations for μ=6 and μ=30 shown in figure 2(a) clearly exhibit this behavior: Poisson→GUE→Poisson. The above behavior of the Creutz ladder model is similar to that in other flat-band models in [60][61][62]. The previous studies focus on a single-particle spectrum, however the Creutz model here includes interaction. The novelty of the results in figure 2 is that even for interacting many-body cases, the level statistical changes first from Poisson to GUE-like, and return to the Poisson. We understand our findings as follows. The Poisson statistics for the μ=30 ensemble originates from the conventional AL that is induced by disorder. Contrastingly, the Poisson-like statistics for the μ=1 ensemble arises from the flat-band properties of the model. Crossover takes place from the flat-band localization to the disorder-induced AL as the disorder increases 8 . This conclusion is in good agreement with the observation in section 2 and will be corroborated by the subsequent IPR calculation. Figure 2(b) shows the LSA of the interacting cases with a weak disorder, μ=1. We find that even for finite interactions V=1 and 6, the Poisson-like statistics persists. This result is indicative of the disorder-free MBL induced by the flat-band structure. This is again in good agreement with the observation in section 2 We also study the non-flat-band case (t=6t 0 ), which we regard as a reference system with respect to the AL in finite-size systems. Figure 2(c) shows the LSA of a non-interacting non-flat-band for various μʼs. The μ=1 and μ=6 results are close to GUE, whereas for a larger disorder, μ=60, the conventional disorder-induced AL occurs. This delocalization-like behavior is robust to the interaction, as shown in figure 2(d). The obtained result, in particular for the non-interacting case, seems to contradict the common belief that all the states are localized in 1D random-potential systems. Probably, this is a finite-size effect, i.e. for a weak disorder, μ=1, localization lengths of certain part of states are larger than the system size. By comparing the results in figures 2(a) and (b) with those in figures 2(c) and (d), we find that the localization in the flat-band case is obviously stronger than that in the non-flat-band case, indicating that their mechanisms are different as we discussed in section 2. We will confirm this observation by calculating other quantities. The level spacing ratio in separate energy sectors is numerically studied in section 3.2 to complement the above LSA. In addition, we investigate finite-size effects for the LSA in figure 2(b). It is displayed in appendix A.

Averaged level spacing ratio (LSR)
The LSR is often used for study of localization, which is a kind of numerical analysis of the LSA [7,64]. In this section, we study the energy-resolved LSR to see the localization tendency of various energy sectors. To this end, we introduce a normalized energy scale ò i , which is defined by where E 1 and E N D are the ground state and maximum excitation energies as before. By definition, 0ò i 1. LSRs of the energy eigenvalues {E i } (in ascending order) are defined as . To obtain average value ( ) á ñ  r as a function of ò, we average r k over 1000 energy eigenstates in the vicinity of ò and 20 disorder realizations. The value of ( ) á ñ  r gives us an estimate of the (non-)localization tendency of the states around the energy density ò. For the Poisson random matrix ensemble (localized state), á ñr 0.386. On the other hand, for an ergodic state (extended state), á ñr 0.600 (GUE). As we show, ( ) á ñ  r in the present system varies from 0.4 to 0.55. This result indicates that coexistence of extended and localized states is realized.
For the flat-band case (t 0 =t 1 ) in figure 3, we display V-dependence of ( ) á ñ  r with the strength of the disorder μ=1 and μ=6. Let us see V=0, μ=1 data first. All ( ) á ñ  r s are close to the value of the Poisson distribution (∼0.386), but in the intermediate energy region (ò∼0.6), the upward deviation from the Poisson distribution exists. This tendency increases for the weak interaction V=1, whereas in the larger interaction cases V=3 and 6, the tendency is weakened. Therefore, even though there is a small ò-dependence in ( ) á ñ  r , the whole states tend to localize in the weak disorder and flat-band case. This result supports the result in figure 2 in section 3.1. In passing, V-dependence in ( ) á ñ  r in figure 3(a) may imply a V-induced weak spectral transition [65].
Let us turn to the μ=6 case in figure 3 has larger values in all cases compared with the μ=1 case. Maximum value of ( ) á ñ  r is 0.55, which is close to the GUE value. Therefore, we expect that extended states exist in the region of μ=6, and they are located in the center of the energy spectrum. This observation is in good agreement with the studies of the IPR and the dynamical behavior of the Creutz ladder given in the subsequent sections.

IPR and crossover
We calculate the IPR, which is often used for the study of localization. By diagonalizing the Hamiltonian in equation (1), we obtain all the eigenvectors, | | . For these eigenstates, the IPR is defined as In particular for the AL with N particles, the localization length, R ℓ [in units of the lattice spacing] is given by (IPR) ℓ ;1/(R ℓ ) N . 9 We average (IPR) ℓ over all the states for fixed μ and V. The averaged IPR is denoted by á ñ IPR . Figure 4(a) shows the μ-dependence of á ñ IPR in the non-interacting case (V = 0). For a sufficiently weak disorder (μ  1), the obtained á ñ IPR in both the flat-band (t 0 =t 1 ) and non-flat-band (t 0 =6t 1 ) is small compared with that in the strong-disorder regime (μ10), where the value of á ñ IPR is large owing to the existence of the conventional disorder-induced AL. In the weak-disorder regime, there exists a clear difference in the á ñ IPR of the flat-band and non-flat-band cases 10 , i.e. the value of the á ñ IPR of the flat-band is obviously much larger than that of the non-flat-band, as shown in the inset of figure 4(a). This means that the flat-band system tends to localize more strongly than the non-flat-band system 11 .The origin of this difference is clearly explained in section 2. It is intriguing to see that  á ñ IPR 0.02 gives an estimation of the localization length, R ℓ ;4.0, which is close to the estimation of the localization length given in section 2.
It is interesting to observe that in the vicinity μ∼6, á ñ IPR decreases in the flat-band system, as shown in the inset of figure 4(a). This behavior is in good agreement with the results of the LSA presented in figure 2(a) and the LSR in figure 3. In fact for μ=6, the LSA of the flat-band shows a GUE-like behavior. Again this behavior of á ñ IPR is an evidence of the crossover, and we estimate μ c ∼6.
As our main concern is the MBL state in the flat-band, we study the interacting cases with finite Vʼs. Calculations of the IPR for the case, V=1, are shown in figure 2(b). We find that the value of á ñ IPR of the flatband increases in the weak-disorder regime compared with the V=0 case, and it again decreases considerably near μ∼6 as in the V=0 case. We investigated cases for other values of V and found similar behavior of á ñ IPR . We therefore conclude that MBL exists in the flat-band Creutz ladder model in the weak-disorder regime, reflecting the flat-band structure. Moreover, a crossover from flat-band MBL to disorder-induced MBL takes place as the In the μ10 regime, the conventional disorder-induced AL/MBL phase appears. The system size is L=12, and the particle number is N=3. 9 More precisely for the AL without inter-particle interactions, ( ) , where n labels quantum states occupied by Nparticles, and R n is the localization length of n-th state. See [21]. 10 The small values of á ñ IPR are partly owing to the degeneracy originating from the translational symmetry of the model. In the smalldisorder regime, the breakdown of the translational symmetry is weak. Accordingly, there exist numerous quasi-degenerate states. 11 Simple estimation of the average of R ℓ , á ñ R , is obtained by using the obtained results of á ñ IPR . It gives  á ñ R 8 for the non-flat-band and  á ñ R 4 for the flat-band. This result for the flat-band is reminiscent of the Wannier state in equation (2).
disorder increases. This is one of the main conclusions of this work. In section 3.5, we shall give a physical picture of the above crossover that is obtained by calculating energy-resolved localization lengths.

Detailed study of IPR: energy-resolved analysis
In figure 4, we showed the mean value of the IPR obtained by averaging all eigenstates. We observed that the IPR exhibits a very interesting behavior as a function of the disorder strength μ, i.e. it substantially decreases in the region μ=1.0-10. In section 3.3, we emphasized that this behavior of the IPR is consistent with the LSA and LSR. In this subsection, we investigate the energy dependence of the IPR, (IPR) ℓ , as we studied the energyresolved LSR ( ) á ñ  r in section 3.2. We also study effects of the interactions. Figure 5 shows the disorder (μ) and interaction (V ) dependence of the IPR for states with various energies. Results of the flat-band cases (t 1 =t 0 ) are in figures 5(a)-(d). There, for all Vs except for V=6, the IPR decreases in the region 1μ  10 in all energy eigenstates. In particular, in the central region of ò, this behavior is remarkable. This indicates that all states tend to extend in the region 1μ  10 in the flat-band system. We think that this peculiar behavior (see the results of the non-flat-band case below) stems from the fact that in 'weak disorder' below μ;1, all the states sustain properties of the flat-band localization although energy splitting takes place as a result of the on-site disorder. In other words for 'strong disorder' (μ>10), genuine localization due to disorder takes place as the disorder is strong enough to dominate the flat-band effects. Therefore, a crossover takes place in the intermediate regime 1μ  10, as we explained in the previous sections.
By close look at V=3 case in figure 5(c), we find that the data for ò=0, 0.1, 0.8, 0.9 and 1.0 (i.e. areas of the tail of the energy spectrum) exhibit only a slight decrease in the IPR in 1μ  10. This tendency is stronger for the V=6 case in figure 5(d). There, the data for ò=0, 0.1, 0.8, 0.9 and 1.0 shows almost no decrease in the value of the IPR in 1μ  10. Accordingly, a 'quasi-mobility edge' seems to exist in for V3.
In summary, the IPR of the flat-band regime shows that for small V, as increasing the disorder μ from the flat-band localization, there exists a crossover regime (in 1μ  10) from the flat-band localization to the disorder-induced genuine MBL. In this crossover regime, all states tend to extend, and for larger μ, all states are strongly localized. On the other hand for large V, such a crossover is blown away, and the direct transition from the flat-band localization to the disorder-induced MBL takes place. What states are realized in the crossover regime is an interesting problem. Coexistence of localized and extended states may occur there as ( ) á ñ  r implies. It is also important to study if the above properties of the Creutz ladder are common to other flat-band systems. These are future works.
We also studied the non-flat-band case (t 1 =6t 0 ). Obtained results of the energy-resolved IPR are shown in figures 5(e)-(h). For all Vs, the IPR for μ  1 is much smaller than the IPR of the flat-band case in figures 5(a)-(d). This result is consistent to the result in figure 4. The behaviors of the IPR for each ò in 1μ  10 are almost the same with the different Vs. Furthermore contrary to the flat-band case, states only located in the central region of the energy spectrum tend to extend and low and especially high-energy states tend to localize there. This behavior comes from the fact that in the non-flat case, quantum states have different features with each other depending on their energy. Close look at the data reveals that in the weak disorder regime 10 −2 μ  1.0, the band-edge states (low and high energy states) start to localize. This is a common picture of the weak localization. Data seem to indicate that there exists a transition from the weak disorder to strong disorder as μ is increased. However, location of this transition may depend on each state. This behavior is similar to that in the conventional MBL transition in a typical random Heisenberg spin chain [64].

Distribution of localization length and emergence of 'critical edge'
In the previous two subsections, we calculated the IPRs to study AL and MBL by varying the strength of the disorder. The results showed the sharp contrast between the flat and non-flat cases. In this subsection, we study the distribution of the localization length as a function of energy by using the relation between the IPR and localization length, ( This investigation is important to examine the finite-size effect, and to verify that the results of IPR obtained in the previous subsections for L=12 are reliable. To this end, we study the systems with L=12 (N = 3, 24 sites) and L=16 (N = 4, 32 sites) focusing on some interesting disorder strengths, μʼs. Besides the finite-size effect, this investigation reveals very important properties of the present system, as we see later in the present subsection.
We first show the distributions of the localization length averaged over the entire energy eigenstates, which correspond to the IPR in figure 4. We consider the flat-band cases with V=1, μ=0.8 and V=1, μ=6, and the non-flat-band case with V=1, μ=6. The results for the system size L=12 are displayed in figure 6(a). For the flat-band case, the localization length for μ=6 is larger than that for μ=0.8, which agrees with the calculations of the IPR in figure 4. More important observation is that the majority of the localization lengths in the distribution in both cases are fairly small compared with the system size, i.e. {R ℓ }<9. This result seems to indicate that the system size L=12 is large enough to calculate the localization length for the flat-band case. On the other hand for the non-flat band case, typical localization length R ℓ ∼9, and therefore, the localization length may not be estimated correctly.
To examine the above observation for the flat and non-flat cases, we studied the L=16 system. Obtained results are shown in figure 6(b). For the flat-band case, the maximum of the localization length { } ℓ  R Max 9, which is the same with that in the L=12 case. For the non-flat band case, the the maximum of the localization length is slightly larger than that in the L=12 case, but { } ℓ  R Max 10.5. These results seem to indicate that the estimations of the localization length are reliable for both the flat and non-flat cases with the above parameters.
We also studied the energy-resolved localization length for the L=16 system with the above parameters and obtained very important observations. The calculations are shown in figures 7(a) and (b). For the non-flat-band system in figure 7(a), the distribution is dominated by a sharp peak and there is a moderate peak very close to the sharp peak. For the flat-band case with a small chemical potential μ=0.1, the distribution has a single moderate peak centered at R ℓ =3.
On the other hand for the flat-band case with μ=6 in figure 7(b), the distribution has a different shape depending on eigenenergy. States in the band center have a fairely large localization length, whereas at the band edges, the states are localized. From the calculations of the non-flat band case with μ=6 and flat-band case with μ=0.1 in figure 7(a), we observe that the states far from the band center are localized due to the flat-band localization, which is one of the properties of the genuine Creutz ladder system. On the other hand for the states in the band-center regime, AL caused by the disorder potential is the main mechanism of localization as in the non-flat system. In other words, there exists a critical strength of the disorder at which the flat-band structure shown in figure 1(b) is destroyed, and the upper and lower bands merge. In this sense, there exists a 'critical edges' separating the flat-band localized and AL regimes. (We estimate them as ò=0.3 and ò=0.85, respectively.) Schematic picture is shown in figure 8, which displays an intuitive understanding of the crossover observed by the LSA and IPR. Anyway, more detailed study on this kind of crossover is a future problem.

Ergodicity-breaking dynamics
The above results of the LSA and IPR indicate that disorder-free single-particle localization and MBL occur in the flat-band Creutz ladder. This motivates us to simulate the dynamics of the Creutz ladder. In the conventional disorder-induced AL and MBL, information of an initial density wave pattern is stored for long times [3][4][5][6]8]. This behavior is a hallmark of ergodicity breaking and indicates the breaking of the eigenstate thermalization hypothesis [3-6, 8, 66]. Here, we focus on the disorder-free cases and investigate whether the flat-band Creutz ladder exhibits ergodicity breaking dynamics. To this end, we employ the time-dependent exact diagonalization method with the periodic boundary condition [67,68].
As discussed in section 2, a particle wave function in the flat-band regime tends to have a non-vanishing amplitude only on definite adjacent finite sites. Therefore, we expect that the localization of the flat-band system exhibits different behavior depending on the particle fillings. This expectation obviously comes from the observation that the Pauli exclusion principle and the repulsion work substantially at large fillings but less effectively at low fillings. In fact at large fillings, the repulsions between particles come to effective, and they suppress movements of the particles. As a result, localization is enhanced.
To verify the above expectation, we investigate three cases of particle filling, 1/8, 1/6 and 1/4-fillings in our numerics. To see their dynamics, we prepare a specific initial state for each filling such as where q is taken as follows for each filling, q=1/8, 1/6 and 1/4, respectively. This initial state is a totally nonentangled Fock state, and therefore it is quite suitable for detect the localization dynamics [8]. To characterize the localization dynamics, we measure the long-time average of the return probability [14,69] ( )  where P(t) is a return probability at t. Here, as shown in [5,14,69], if the initial state is given by | | ℓ ℓ ℓ where ò ℓ is the eigenenergy of | ℓ y ñ. Accordingly, P is related to the level spacing of the eigenenergy. That is, the above expression of P indicates that states with small level spacings contribute more to P. Since the Poisson distribution (realized in localized regimes) has a small level spacing regime, the value of P tends to be large. Therefore, localization enhances the value of P. In a conventional localization state, entanglement of eigenstates is fairly suppressed, and each eigenstate | ℓ y ñ tends to be close to the Fock state | ñ F m . Then, a finite (IPR) ℓ implies a finite P although P is indirectly related to IPR defined in section 3.3.
If the value of P is finite, the memory of the initial state is preserved. This implies that an ergodicity breaking takes place and the system exhibits localization. In our practical numerics, we set the unit of time by  t 1 , set the long-time limit as t=10 3 [ ]  t 1 in equation (11), and use the time slice, dt=10 −3 [ t 1 ]. We put μ=0 for all the calculations.
To begin with, let us verify a single particle localization dynamics. The initial state is set to | | † y ñ = ñ = a 0 j ini 8 with the system size L=15. Figure 9(a) shows the dynamical behavior of The numerical result exhibits a clear localization since oscillates, and also its oscillating period agrees with the analytical result of equation (9). Under the flat-band condition t 0 =t 1 , the single particle certainly localizes. The detailed density dynamics for rung j=8, 8 , 8 is also plotted in figure 9(b) for the flat-band case and (c) for non-flat band case. For the flat band case, the initial single particle is localized with an oscillation between j=7 and j=9 rungs, corresponding to the analytical result of equation (9). On the other hand, see figure 9(c) for the non-flat band case ¹ t t 1 0 with V=0, the oscillation of decays immediately. Let us turn to the multiple-particle system. We calculated P(t) for various filling cases and interaction strengths V. First, we consider the 1/8-filling case. We expect that the inter-particle distance of the initial state is sufficiently large there, and the particles do not substantially interact with each other. However, this does not necessarily mean that the repulsion does not influence the dynamics of each particle at all. The numerical result for the L=12 three-particle system is displayed in figure 10(a). For various Vs in the flat-band case, P(t) takes a finite large value for long times, i.e. P∼0.37. This indicates the strong localization of particles and the ergodicity . The green line is the analytical oscillation solution of obtained by equation (9). Density distribution dynamics for the rung, n a,j +n b, j (b) for the flat-band condition, (c) for the non-flat band condition, t 0 =2t 1 .
breaking. The independence of the value of V in the dynamics originates from the large inter-particle distance. On the other hand, the results for the non-flat band case with t 0 =2t 1 show that the value of P(t) suddenly decays, and therefore the dynamics of the non-flat band system is ergodic. These numerical results are consistent with the results of the LSA shown in figure 2(b).
Second, let us turn to the 1/6-filling case. For the initial state of equation (10) in the non-interaction case V=0, we expect that each particle starts to oscillate around the adjacent rungs as described by the singleparticle solution of equation (9) in section 2. There, although each single particle wave function spreads a little, the overlap and interference between particles are not substantial, and therefore we expect the multiple-particle system for V=0 exhibits strong localization similar to the 1/8-filling case above. However, once V is switched on, oscillating single particles start to interact with each other, and the single-particle localization picture may be affected by the existence of the interaction V. Figure 10(b) is the result of P(t) for the L=9 three-particle system. For V=0, as we have expected, P(t) exhibits strong localization P∼0.37. Remarkably for a finite V, P(t) remains a definitely finite value, P∼0.15. Even for a finite V with the 1/6 filling, the system exhibits the ergodicity-breaking dynamics, but the value of P is a little smaller than that of the V=0 case and the 1/8 filling, i.e. the interacting system is moderately localized. For this moderate-localization regime, it is difficult to judge whether the LSA and LSR obey Poisson or GUE ensembles. The calculations of the LSA and LSR shown in appendix B prove this expectation. Properties of the moderate localization are interesting and warrant deep study as a future work.
Third, we focus on the 1/4-filling case. The inter-particle distance is small and the overlap of the single particle oscillating wave functions is so large that we expect the interaction V drastically changes the localization properties of the system. Figure 10(c) is the result of P(t) for the L=8 four-particle system. Interestingly enough, depending on the value of V, the dynamical behavior of P(t) drastically changes. For the non-interacting V=0, the moderate localization appears since P∼0.13. As increasing V from V=0, for weak but finite V cases, the localization is highly suppressed, i.e. the system tends to be extensive since P<0.1. However for large V6, the P increases to P>0.2. That is, the interaction V suppresses the localization tendency first, but it starts to enhance localization as V exceeds a critical value. We expect that in the localized regime for large V, the particles repel each other strongly, and then particles are squeezed and localized moderately. The localization length of the moderate localization for large V may be a little larger than that of the strong localization with P∼0.37. Calculations of the averaged LSR, ( ) á ñ  r , in appendix B suggest that the band-edge eigenstates in the moderate-localized regime have stronger tendency of delocalization compared with the strong-localized states.
We conclude that for small size systems, a disorder-free MBL exists in the flat-band Creutz ladder both with and without interactions, and it exhibits ergodicity-breaking dynamics. In figure 11, we summarize the results of the numerical calculations and show the qualitative dynamical properties of the system as a function of the interaction V for various particle fillings.

Conclusion
In this work, we have clarified disorder-free single-particle localization and MBL phenomena induced by the flat-band structure of the Creutz ladder model. We found that the flat-band localization originates from a localized Wannier state (FB compacton), and the localization length is quite short compared to that of ordinary AL. The localization length of the flat-band Creutz ladder system is also insensitive to the strength of the disorder although the localization length of AL is strongly influenced by the disorder strength. As a result, effects of interactions in flat-band localization depends on particle filling (inter-particle distance) significantly.
In section 3, we extensively studied the flat-band localization properties by using some conventional numerical methods. We extracted the localization properties from the statistical properties of the static spectrum and eigenstates in the Creutz ladder Hamiltonian. In the flat-band regime, the LSA exhibits Poisson distribution in the weak-disorder regime with or without interactions. This indicates that the flat-band model exhibits a (many-body) localization induced by the flat-band nature not by disorder as in AL. After that, we calculated the LSR from the spectrum and also the IPR from eigenstates of the model in order to capture the localization tendency in the real space. We found that they support the LSA result.
We also studied the flat-band localization from the view point of the dynamical aspect. We found that flatband localization tends to prevent the system from thermalization. The single-particle localization picture was analytically given in section 2. If we put on a single particle on a single site on the flat-band system, the single particle localizes with oscillating. To estimate this dynamical localization with or without interactions, we performed exact dynamical simulations for small size systems. To judge the thermalization and ergodicity breaking, we employed the return probability, which quantifies how much information of initial state wave function remains. By calculating the long-time average of the return probability, we characterized an ergodicitybreaking dynamics similar to the conventional disorder-induced AL and MBL dynamics, and also found rich localization properties as varying particle filling and the repulsive interaction. In summary, even in the interacting cases, the system exhibits localization and ergodicity-breaking dynamics. Our numerics is only for small system sizes but exact, and therefore our results can be a benchmark for future simulations with large system sizes, e.g. by using Krylov subspace method. We also expect that the findings in the present work are useful for future real experiments on cold atoms such as [39,41].
We also expect similar phenomena in other flat-band models, such as the saw-tooth, diamond and Lieb lattice models, which are to be realized in experiments [70][71][72][73][74], and also some studies in the Creutz ladder in the clean limit [75,76] pointed out the presence of conserved quantities. Such conjecture may support our results.