Non-ergodic behaviour of clean Bose-Hubbard chains

We study ergodicity breaking in the clean Bose-Hubbard chain for small hopping strength. We see the clear existence of a non-ergodic regime by means of indicators as the half-chain entanglement entropy of the eigenstates, the average level spacing ratio and the dynamics of the Peres Loschmidt echo. We find that this ergodicity breaking falls out the many-body localization paradigm because the average half-chain entanglement entropy of the eigenstates obeys volume law. This ergodicity breaking appears unrelated to the spectrum being organized in quasidegenerate multiplets at small hopping and finite system sizes, so in principle it can survive also for larger system sizes. We find that some imbalance oscillations in time which could mark the existence of a glassy behaviour in space are well described by the dynamics of a single symmetry-breaking doublet and well captured by a perturbative effective XXZ model. We show that the amplitude of these oscillations vanishes in the large-size limit.


I. INTRODUCTION
Thermalization in classical Hamiltonian systems is well understood in terms of chaotic dynamics and the related essentially ergodic exploration of the phase space [1][2][3]. From the quantum point of view the physical mechanism is quite different, involving the eigenstates of the Hamiltonian being fully random strongly entangled states which appear thermal from the point of view of local measurements. This is the paradigm of eigenstate thermalization (ETH), introduced in [4][5][6].
In many cases there is correspondence between classical and quantum thermalization. Usually, when a classically chaotic Hamiltonian is quantized, it gives rise to a Hamiltonian matrix which is a random matrix and its eigenstates have exactly the required properties for local thermalization. This correspondence is nevertheless highly non trivial, because quantum effects can give rise to ergodicity-breaking phenomena with no analog in the classical domain. From a classical point of view, any nonlinear Hamiltonian system with more than two degrees of freedom and no conservation law beyond energy gives rise to chaos and eventually thermalizes. From the quantum point of view the situation is different.
A striking example is the dynamical localization (initially discovered for one degree of freedom [7] and later generalized to the many body case [8][9][10][11]), where a classical ergodic driven system behaves regularly in the quantum regime, leading to a suppression of energy absorption. Another example is many body localization. In contrast with the classical case, a generic disordered and interacting many-body quantum system does not thermalize and spontaneously generates the localized integrals of motion which forbid ergodic behaviour (see an overview of the subject in the reviews [12][13][14]). * On leave It is therefore of the utmost interest to find new cases of complex many body systems which break ergodicity, in order to understand if there are general properties in this lack of correspondence with the classical case. One very interesting advancement in this field has been the proposal of many-body localized systems without disorder. Starting from the early proposals on this subject [15][16][17] there has been a constant interest in understanding if disorder is a necessary ingredient to achieve localization in an interacting system. Most of this activity concerns proposals based on gauge theories which are globally uniform and generate a different realization of an effective disorder in each superselection sector [18][19][20][21][22][23][24]. A uniform Josephson junction chain was studied in [25] where it was shown that at high energies or small Josephson coupling to the system break ergodicity, show a many-body localized phase with a consequent absence of conductivity. The systems considered in these works are not integrable, thus pointing out that non-ergodic behaviour can be found also in non-integrable clean models. Our interest is to further explore these questions and seek for non-ergodic behaviour of a clean quantum many-body system.
Aim of this work is to analyze the the ergodicity properties of the clean Bose-Hubbard model [26]. The equilibrium phase diagram and the dynamical properties of the model were extensively scrutinized in the last three decades also for its importance in the physics of optical lattices [27]. The many-body localization in the presence of disorder was studied in [28][29][30][31][32]. The non-ergodic behaviour of the clean version of this model was already discovered and studied in the two papers [33,34]. Here we confirm and extend the results of these papers. We inquire first of all the relation of the non-ergodicity with the spectral structure of the system. Second, we study if this ergodicity breaking can be interpreted in a manybody localization paradigm. Finally, another goal will be to understand how the glassy behaviour discovered arXiv:2005.09892v1 [cond-mat.quant-gas] 20 May 2020 for this model in [15] is related to non ergodicity.
The Bose-Hubbard model, describing a system of interacting bosons hopping on a d-dimensional lattice, is characterized by two energy scales, the hopping amplitude J and the on-site repulsion U . For J = U the model is known to behave quantum chaotically [35]. Here we extend this analysis and consider several different indicators to examine the (non-)ergodicity as a function of the ratio J/U .
First of all we consider the half-chain entanglement entropies of the eigenstates of the system. We find that their average always obeys a volume law. For large J this volume law is the same obeyed by a fully random state (Page value): In this regime the system obeys eigenstate thermalization [4][5][6]36] and is fully quantum chaotic and thermalizing. For small J ergodicity is broken and the pre-factor of the volume law is significantly smaller than the Page value. Looking at the dependence of this value as a function of J we can find that the ergodic behaviour settles at J * 0.6U . For smaller values of J, this model breaks ergodicity in a way remarkably different from many-body localization, where the eigenstates show instead an area-law behaviour and the averaged entanglement entropy is constant with the system size.
This crossover at J * 0.6 is confirmed by the analysis of the level spacing distribution. Keeping the focus on the spectral properties, we discover that this non-ergodic behaviour appears unrelated to the spectrum being organized in multiplets at small J and small sizes. In this circumstances, in fact, the spectrum can be understood through perturbation theory in J. At J = 0 the spectrum displays large degenerate subspaces. At a finite L, as a small J is switched on, these multiplets acquire a bandwidth of order LJ, remaining well separated between each other by energies of order U . However, we expect the multiplet structure to ultimately disappear for any J, when L U/J. One of our main results is noting that signatures of non-ergodic behavior survive even in cases where the multiplet structure breaks down. This is a crucial result as it disentangles the multiplet structure from non-ergodicity, hinting that non-ergodicity might survive to large system sizes (while the multiplet structure will eventually disappear).
We also study the crossover from ergodicity to non ergodicity by inquiring the dynamics of the Peres Loschmidt echo. We find a wide crossover region, and we see that the value J * 0.8 falls in that region. One possible interpretation is that at system sizes larger than the ones we can access this region shrinks. Nevertheless, for the system sizes we can access with our numerics, we can only clearly state the existence of a regular and an ergodic regime and a transition region between them.
Some words of cautions are necessary at this stage. We do see a non-ergodic regime at small hoppings with characteristics that are clearly different from a many-body localized phase. In this way we confirm and extend the results of [33,34]. However, given the system sizes that we are able to reach, we should keep in mind the possi-bility that this non-ergodic regime is a finite-size effect and we do not know if it can be extrapolated towards the thermodynamic limit. The same limitation does not allow to make a clear-cut observation of the value of the transition point to the ergodic phase.
We conclude our analysis studying the relation of non ergodicity with the glassy behaviour found in [15]. We prepare the system in a number-imbalanced state and we see that, for J small enough, the imbalance shows oscillations with a frequency decreasing with the system size, confirming to longer times the behaviour observed in [15]. We find that the frequency decreases exponentially with the system size. Using perturbative arguments, we show that this phenomenon is described by Rabi oscillations in the symmetry-breaking doublet of an effective XXZ model. The perturbative scheme itself depends on the multiplet structure of the spectrum and both disappear for large system sizes.
The paper is organized as follows. In Section II A we introduce the model and in Section II B the quantities and the observables we use to describe it. In Section III we study the ergodicity breaking by means of the behaviour of the half-chain entanglement entropy whose volume-law behaviour is in sharp contrast with the arealaw behaviour in many-body localized systems. In Section IV we study the ergodicity breaking by means of the spectral properties and we see that it is unrelated to the multiplet structure of the spectrum. In Section V we study the behaviour of the Peres Loschmidt echo. In Section V we study the imbalance dynamics and we interpret the oscillations we find as Rabi oscillations involving a single symmetry-breaking doublet. In Appendixes A and B we discuss in detail the perturbative analysis allowing us to construct an effective XXZ model well describing this phenomenon. In Section VII we draw our conclusions.

A. The model
The one-dimensional Bose-Hubbard model [26] is characterized by the Hamiltonian the first term describes the on-site repulsion while the second the hopping between neighboring sites. In Eq.(1) L is the system size,â j are bosonic operators andn j ≡ a † jâ j are number operators, and j labels the site. In the rest of the paper we will fix the interaction strength U = 1; we will explore the different dynamical behaviours modifying the hopping strength J. The total number of bosons operatorN ≡ jn j is conserved. Defining the filling factor ν ≡ N /L, it is possible to see that the Hilbert space dimension is dim H = L(ν+1)−1 L−1 (if Lν is an integer). This implies that, in order to see if there is thermalization, we have to compare the behaviour of the observables with a finite T = ∞ value. Here and, for small system size, we do not need any truncation. In our analysis we will restrict to the case ν = 1 (one boson per site) and label the state using the Fock basis, i.e. the basis of the simultaneous eigenstates of all then j operators (we will denote it as {|n }).
Throughout the work we will consider periodic boundary conditions (unless otherwise specified), so the system has the translation and the inversion symmetries. Therefore, choosing initial states invariant under these symmetry operations, the dynamics restricts to the Hilbert subspace H S (L) fully symmetric under these symmetry operations. This fact restricts the dimension of the interesting Hilbert subspace, so we can perform full exact diagonalizations for system sizes up to L = 11. In particular, the fully-symmetric subspace has a value of the dimension smaller than the full Hilbert space dimension of a factor ∼ 2L.
The existence of an ergodicity breaking phase in this model has been yet argued in [33,34].

B. Quantities and observables
Our analysis will include both observing the dynamics of the system after an initial preparation in a a given non-equilibrium (high-energy state) and by a statistical analysis of the properties of the spectrum. Correspondingly we will consider different quantities.

Eigenvalues and eigenstates
We define the Hamiltonian eigenvalues as E α and the corresponding eigenvectors |ϕ α . In order to see if there is eigenstate thermalization we focus on the properties of eigenstates.
We consider an important basis independent quantity, the half-chain entanglement entropy, which is defined in the following way. We divide the system in two partitions A and B. When L is even they are both long L/2. When L is odd, one of them is long L/2−1 and the other L/2+1. So the full Hilbert space (not the symmetrized one) has a tensor product structure H(L) = H A ⊗ H B . Once we have done that we define for each eigenstate where Tr B is the partial trace over H B . Notice that now the states are taken living in the full Hilbert space. In order to pass from |ϕ α S to |ϕ α one must apply a linear transformation (an isometry). We will see that the entanglement entropy is a very precious quantity to see if there is ETH and thermalization. We study also the properties of the energy eigenvalues, density of states Eq. (15) and average level spacing ratio Eq. (14) (see Sec. IV). Throughout all the text we define the average over the eigenstates as

Dynamics
We will study dynamics with different initializations. First of all, in order to explore space localization, we will follow a protocol introduced by [15,37] and focus on the evolution of the imbalance. In order to do that, we will initialize the system in the state and study the evolution of the imbalance operator In this analysis we can restrict to a portion of the Hilbert space (the one fully symmetric under translations of two sites) which is larger than the fully symmetric one, but still numerically affordable. For L ≤ 10 we will use full exact diagonalization. For L = 12, on the opposite, we will resort to Krylov technique (implemented in Expokit [38]) and we will truncate the Hilbert space so that the number of bosons per site will be smaller than some threshold m th (we will always consider m th ≤ 8). We will take a Krylov subspace of dimension M K ≤ 32. With this technique we can address only values of J 0.25: in this case, for the sizes we consider, the time-evolved state does not deviate too much from the initial one. In order to understand if the Krylov technique is applicable with our numerical resources, for each of the considered cases we consider different values of m th and M K and verify that the result converges if we increase these values. Inspired by [39], we will also consider the following protocol. We will initialize the system in the state and we will perform two different evolutions, one with HamiltonianĤ[J] and another with HamiltonianĤ where the hopping strength has been increased byĤ δ ≡ H(J + δJ). Also here we use exact diagonalization and Krylov technique, as we better specify later. We will study the evolution of the Peres Loschmidt echo It is important to state the definition of the infinite-time average of a time-dependent quantity O(t),

III. BEHAVIOUR OF THE ENTANGLEMENT ENTROPY
We start our analysis by considering the half-chain entanglement entropy of the eigenstates Eq. (2). We plot some examples of scatter plots of S (α) L/2 versus E α for different values of L and J in Fig. 2. As a reference we plot also the line corresponding to the "Page value" S (Page) L/2 . With that phrase we mean the value obtained with the following construction, performed for a 1/2-spin chain model by Don N. Page [40]. We take a state |ψ S such that with α n,ψ and β n,ψ drawn from from a normal distribution with zero mean and N ψ the appropriate normalization constant [41]. We project it in the full Hilbert space and we evaluate the entanglement entropy. We repeat this procedure over N rand realizations and we take the average. We see that for J above some threshold, the graph becomes a more-or-less continuous curve which touches the Page value. This marks the setting-on of the ETH. Therefore we can qualitatively see that the system breaks ergodicity at small J and obeys ETH at large J. Now we are going to discuss this feature more quantitatively.
In order to do that, we move to study the average over the eigenstates of the entanglement entropy. We show some examples of scaling with the system size L in the inset of Fig. 1. We notice that there is always a linear increase, no matter the value of J, but the difference in slope between small and large J is striking. In the former case (J = 0.01, J = 0.0225) there is ergodicity breaking and the slope is smaller. In the ergodic cases (J 2.92, J 9.8) we show there is still linear increase but with a larger slope, very near to the one of the fully-random Page value.
Indeed, if ETH holds, then the slope of the average entropy should coincide with the slope of the Page entropy. To show this, we recall that, assuming ETH and considering a subsystem of l consecutive sites with l L, the reduced density matrix of an eigenstate is equivalent to a thermal density matrix. In particular, it follows that the entanglement entropy of the reduced density matrix is equal to the thermal entropy computed over the appropriate ensemble. This fact has been shown to be valid up to l = L/2 (with better approximation for increasing L) for a spin 1/2 Heisenberg chain in Ref. [42].
We assume the reduced density matrix being thermal for l = L/2 also in our case (we will see that it leads to conclusions confirmed by our numerical analysis). With this assumption, we can compute the S L/2 for eigenstates at energy E as the corresponding thermal entropy S(E) in the microcanonical ensemble. In the definition of S(E) an average over an energy shell of width ∆E is implicit. We can also average S L/2 over eigenstates in the energy shell [E − ∆E/2, E + ∆E/2] in order to get more regular S L/2 (E) curves. In the microcanonical ensemble S(E) = log (ρ S (E)), where ρ S is the density of states with energy E with E α running over all the energy eigenvalues. (Also ρ S (E) is regularized by averaging over the energy shells). At this point it is convenient to express entropies S and energies E in terms of the corresponding energy densities ε = E/L and entropy densities s = S/L. If we take an energy-shell width ∆E much smaller than the width of the energy spectrum, we can write the entanglement entropy averaged over the eigenstates as Since ρ S (ε) = exp(Ls(ε)), when L 1, we are justified in computing the integrals using a saddle-point approximation. Thus replacing where ε ∞ is the energy density of the maximum entropy states (corresponding to T = ∞), s(ε ∞ ) coincides with the Page value, and W denotes a non-universal manybody bandwidth. The saddle point calculation then gives with the O(1) correction being non-universal, as they depend on W . We can confirm the analysis above by doing a linear fit S L/2 ∼ A + B S L. We show the dependence of the slope B S on J in the main panel of Fig. 1. We can notice that the linear coefficient increases until reaches a plateau around J * ∼ 0.6. This plateau coincides with the Page value inside the errorbars. Looking at Fig. 3 we see that around J * the average level spacing ratio r attains the Wigner-Dyson value. So, we can say that for J * 0.6 the system has reached ergodicity.
So, there is a clear crossover from non-ergodicity to ergodicity, but in both regimes the half-chain entanglement entropy obeys a volume law with the system size. This is a relevant result and clarifies that the system in the ergodicity-breaking regime behaves in a way different from many-body localization where the half-chain entanglement entropy obeys area law. Of course these statements apply only to the system sizes we can access numerically and no statement can be done for larger system sizes. Nevertheless, we are quite confident that at J * there is a strong modification in the dynamical properties of the system, as we are going to further confirm in the next section by means of the spectral properties of the system.

IV. SPECTRAL PROPERTIES
The quantum chaoticity properties of a system appear also through the properties of the spectrum of its Hamiltonian. An important role is played by the level spacings λ α = E α+1 −E α . Indeed, the distribution of the (normalized) level spacings takes a universal form in case of quantum chaos. If the dynamics is fully chaotic and thermalizing, the Hamiltonian behaves as a Gaussian-Orthogonal-Ensemble random matrix [43][44][45][46] and the level spacings obey the Wigner-Dyson (WD) distribution. On the opposite, a classically integrable system generically shows a Poisson distribution [47] of the normalized λ α . We remark that the Poisson level spacing distribution is just a sufficient condition for integrability, not a necessary one [47,48]. A very convenient tool to distinguish these extreme cases and all the intermediate ones is the average level spacing ratio This quantity attains a value r WD 0.5295 for a fullychaotic Wigner-Dyson level-spacing distribution and r P 0.386 for a Poisson distribution. In many-body localization phases there is a superextensive number of localized integrals of motion, a situation closely resembling classical integrability, and the level spacing ratio attains the Poisson value. In this context, the transition to ergodicity is marked by a crossover to the Wigner-Dyson value which becomes sharper and sharper as the system size increases [49].
We evaluate the average level spacing ratio for our model and we plot it versus J for different system sizes in Fig. 3. We see that it attains the Wigner-Dyson value for J 0.8 and that this feature is stable when increasing the system size. The crossover does not become sharper and sharper as the system size is increased as it occurs in contrast for many-body localization systems (see for instance [49]). The value of r intermediate between Poisson and Wigner-Dyson marks that the system is not ergodic (the Hamiltonian is significantly different from a random matrix) and does not thermalize. From this fact we also know that this ergodicity breaking is most probably not associated to any space localization. Indeed, from the existing numerical evidence [12] we know that whenever there is a strong form of localization (Anderson or manybody localization) the average level spacing ratio gets the Poisson form.
We remark that the level spacing ratio moves to the Wigner-Dyson value at a value of J consistent with the value of J * we have found for the crossover in the behaviour of the entanglement entropy (see Section III). We also to remark that the nonergodic behaviour of the level spacing ratio is independent of the fact that the spectrum is organized in multiplets when J and L are small. Indeed, the spectrum is massively degenerate for J = 0; when J U one can apply second order perturbation theory in J/U (see A for details) and see that the degenerate levels move to quasidegenerate multiplets. We can see an example of this fact in Fig. 4 (an analogous figure can be found in [50]). Here we show the density of states for different values of J and L. We can see that for small J the density of states shows a series of spikes, each one corresponding to a quasi-degenerate multiplet.
Nevertheless we see a tendency of this structure to vanish and get diluted in a smooth continuum as the system size L is increased (this fact is quite apparent for J 0.17 and J 0.25). Moreover, already at J 0.38 there is no trace of this structure when the system size is L = 9. On the opposite, the average level spacing ratio attains the Wigner-Dyson value around J ∼ 0.6 and the entanglement-entropy slope attains the Page value at J * ∼ 0.6. So, the organization of the spectrum in multiplets and the ergodicity breaking appear to be independent phenomena.
We can get a further confirmation of this fact by looking at the energy gap between two nearby multiplets. We choose in particular the gap between the multiplet around the energy density 0.5U and the one immediately above it. The 0.5U is the energy density of this multiplet when J = 0 and it is degenerate. As we show in the inset of Fig. 5 this degeneracy is resolved for J = 0, and for J large enough nearby bands merge with each other. If we fix L and we change J we see that the number of states in each multiplets does not change. So, in order to evaluate the gap between the two nearby multiplets we are interested, we need to order the eigenstates in increasing energy order and compute the difference between the two eigenvalues on the two sides of the vertical black line in the inset of Fig. 5. More formally, this operation defines the gap as the difference of the minimum energy of the upper multiplet and the maximum energy of the lower multiplet. We show the dependence of this gap on J in the main panel of Fig. 5. We see that for L = 10 the gap decreases below 10 −3 already at J 0.17 and for larger values of J one cannot speak about multiplets separated by a gap anymore.
Nevertheless, the results of the average level spacing ratio r should be taken with caution. What we have said is certainly true for the system sizes we can numerically address, but there might be unexpected developments for larger system sizes. For instance, the slight increase with L that one can see between J ∼ 0.1 and J ∼ 0.4 in Fig. 3 could lead to a convergence to the Wigner-Dyson value in the thermodynamic limit. A similar very slight increase with L can be found in the average level spacing ratio resolved in energy, as one can see in [50]. What we can say for sure is that at finite system sizes there is a clear distinction in the behaviour of r at large J where it attains the Wigner-Dyson value and at small J where it does not; the two regimes seem quite robust when one varies the system size and respectively correspond to an ergodic and a non-ergodic behaviour in the half-chain entanglement entropy; they seem to have no relation with the multiplet structure of the spectrum. In the next section we will see that the transition between these two regimes can also be seen in the dynamics by means of the Peres Loschmidt echo.

V. PERES LOSCHMIDT ECHO
We focus here on the Peres Loschmidt echo P(t) defined in Eq. (7). We consider δJ = 10 −3 and plot some instances of time traces in Fig. 6. We see that, for J 1, P(t) oscillates around 1 and changes very slightly with the system size (inset of Fig. 6). For larger J it decays and then oscillates around an average value (main panel of Fig. 6). The last behaviour is a marker of quantum chaos in small quantum systems [39].
In a many-body context like the present one, in order to understand if there is quantum chaos we consider also the scaling with the system size L. We consider the time-average of the Loschmidt echo P and we plot it versus J for different system sizes in the main panel of Fig. 7. We qualitatively see that for J 0.1 P is almost size independent (see an example in the inset of Fig. 7). We have also checked that for J 0.8 there is an exponential decay with the size (see examples in the inset of Fig. 7). Between J0.1 and J = 0.8 there is a fluctuating behaviour in L without clear features, which possibly develops in a clear trend for larger system sizes. We remark that the curve for L = 7 crosses the one for L = 8 for J 0.6, while the curves for larger sizes cross only in J = 0. So, we have identified a clear regular regime for J 0.1 and an ergodic regime for J 0.8. In this regime, indeed, the Peres Loschmidt echo exponentially decays to zero with increasing system size, as appropriate for a fully chaotic system. In this regime, we can fit log(P) versus L with the curves log(P) = A P + B P L and we plot the slope B P versus J in Fig. 8. We consider two values of δJ (10 −3 and 10 −4 ) and we see that the curves do not differ significantly from each other, marking the robustness of the Peres Loschmidt echo as a detector of chaos. We have also checked that the curves in Fig. 6 do not change significantly if we change δJ from 10 −3 to 10 −4 : The behaviour of the Loschmidt echo is robust in δJ in all the regimes we find.
To summarize, due to the limited system sizes we can access to, we cannot pinpoint a precise crossover point, but a crossover region between regularity and chaos for J ∈ [0.1, 0.8]. We can just notice that the crossover point J * 0.6 that we found in the sections above is contained in the crossover region we have found here, so the results found with entanglement entropy, spectral analysis and Peres Loschmidt echo are in agreement.

VI. DYNAMICS OF THE IMBALANCE
We start considering some instances of evolution of the imbalance I(t) = ψ(t)|Î |ψ(t) withÎ defined in Eq. (5). We consider also the evolution of the overlap in Fig. 9 we plot some time traces of the imbalance (upper panel) and of the overlap (lower panel). For J = 0.01 we qualitatively see that there are large-amplitude longperiod regular oscillations with a period strongly increasing with the system size (we will see that the increase is exponential). On the opposite, for J = 10 there are very small irregular oscillations resembling a random noise (as appropriate for local thermalization, as we will have seen above the system obeys ETH in this regime). Indeed, for different values of J, we can see two very different regimes. In order to better understand them and their connection with the crossover to thermalization we have discussed above, let us move to a more quantitative analysis of the amplitude and the period of the imbalance oscillations. We quantify the amplitude of the oscillations by means of the infinite-time fluctuations of the imbalance Energy eigenvalue densities Eα/L in increasing order versus progressive number α for L = 8. We focus on the band around energy density 0.5U . We see that the number of states in the band does not depend on J. We can evaluate the gap between the two nearby bands we are interested in as the difference of the two energy eigenvalues across the vertical black line. states and eigenvalues of the Hamiltonian and assume that there are no degeneracies in the spec- trum, we easily obtain the formula where we have defined I α γ ≡ ϕ α |Î |ϕ γ . Notice that Λ is the IPR of the initial state |ψ 02 in the basis of the Hamiltonian eigenstates. We show ∆I 2 versus J in Fig. 10 upper panel and Λ versus J in Fig. 10 lower panel. We consider different system sizes and we see that for J around 0.17 the curves drop down to a plateau. In this plateau we see that ∆I 2 and Λ decrease by an order of magnitude when the size is increased by two sites. Moreover, the curves drop to the plateau at a value of J which steadily decreases as L increases. We remark that the curves decrease with increasing L also for small J (although in a much slighter way). We can quantitatively estimate also the period of these oscillations by looking at the Fourier transform of the imbalance and the overlap. We define the Fourier transform as (we consider for instance the case of the imbalance) where t f is the total evolution time. Performing the Fourier transform for I(t) and Λ(t), we consider the main peak at non-vanishing frequency (Λ(ω) has a huge peak at vanishing frequency due to the non-vanishing time average of Λ(t)). We find that these two quantities coincide and we plot them in the upper panel of Fig. 11. This peak frequency is the frequency of the dominant imbalance oscillations and we term it ω peak . It increases quadratically with J up to around J ∼ 0.17 and decreases exponentially with the system size. More precisely, in this interval of J, ω peak is well described by the formula We explicitly show this in the lower panel of Fig. 11 where we plot log ω peak /J 2 versus L. We see that the curves are actually straight lines and the curves for different values of J overlap when J 0.17. We also perform a minimum square fit of the curve for J = 0.01 with the formula log ω peak /J 2 = log(B) − αL. We find a good agreement for the fit, as we show in Fig. 11, and we get as parameters log(B) = 1.153 ± 0.148, α = 0.496 ± 0.016.
In the plot we show also a curve for a reduced model indicated as "XXZ effective model". We can use secondorder perturbation theory in J/U to interpret the dynamics at small J (J U ) as the dynamics of an effective XXZ model. This model breaks the Z 2 symmetry in its ground state for infinite size. For finite size, there is a symmetric ground state and an antisymmetric firstexcited state separated by an exponentially small splitting. Preparing the system in a symmetry-breaking state, it will show Rabi oscillation with a frequency given by this splitting and described by a formula like Eq. (22). It is important to remark that in general the overlap of the initial symmetry-breaking state with the two smallest- energy states giving rise to the Rabi oscillations is smaller than 1 and exponentially decreases with the system size. So in the thermodynamic limit the Rabi oscillations tend to have vanishing amplitude. This fact is true also for the state |ψ 02 which in the XXZ language is represented as a symmetry-breaking state. So, we do not need to extrapolate our curves to infinite size in order to know that the amplitude of the imbalance oscillations vanish in the thermodynamic limit (as also Fig. 10 suggests) and so there is no many-body localization. We are going to describe our analysis in detail in the following and in Appendixes A, B.

A. Interpretation of the imbalance oscillations
As we show in detail in Appendix A, for J 1, whenever the system size of the system is so small that the gaps between the different quasidegenerate multiplets are still open (see Sec. IV), we can perform a second-order pertur- We see that the overlap is non-vanishing mainly for two quasidegenerate states. As we show in table I, the splitting of these two states is exponentially small in the system size and coincides with the corresponding Rabi-oscillation frequency shown in Fig. 11. bation theory in J/U to understand the level structure in the different multiplets. In particular, if we focus on the multiplet to which our initial state |ψ 02 belongs, we find that the dynamics is induced by the effective Hamiltonian z j (23) witht = (1/2)J 2 /U ,Ṽ = 4J 2 /U andμ = −J 2 /U . This is a XXZ Hamiltonian and |ψ 02 corresponds in this representation to the Néel state |↑↓↑↓ · · · . The remarkable thing about the XXZ Hamiltonian is that it has just one single pair of eigenstates (call them |ϕ ± ) which break the Z 2 symmetry in the thermodynamic limit. As we show in detail in Appendix B, these two states are separated by a splitting ∆ exponentially small in the system size and proportional to J 2 . This is the key leading to the Rabi oscillations. The initial state has square overlap mainly with the two states in this doublet (see some example in Fig. 12), and we have checked that the splitting ∆ between these two states coincides with the frequency ω peak shown in Fig. 11 (see some examples of this fact in table I). Therefore, one can see oscillations with a frequency exponentially small in the system size just because the splitting between |ϕ + and |ϕ − is exponentially small in the system size. Nevertheless, as we show in detail in Appendix B, the overlap of |ψ 02 with the symmetry-breaking doublet tends to 0 in the thermodynamic limit and so the Rabi oscillations disappear in that limit.
This can be also seen in a simpler way, considering that the gaps between the different multiplets tend to vanish L J = 0.01 J = 0.11389 ∆ ω peak ∆ ω peak 6 1.5 · 10 −5 1.5 · 10 −5 2.08 · 10 −3 2.08 · 10 −3 8 6.24 · 10 −6 6.14 · 10 −6 7.68 · 10 −4 7.68 · 10 −4 10 2.53 · 10 −6 2.39 · 10 −6 3.12 · 10 −4 3.12 · 10 −4 TABLE I. Comparison of the splitting ∆ of the two maximumoverlap states of Fig. 12 with the corresponding value of ω peak (lower panel of Fig. 11) for two values of J. The agreement is excellent, up to errors due to the finite time over which the Fourier transform is performed. As the lower panel of Fig. 11 shows, both the quantities decay exponentially in L, as confirmed by the theoretical analysis of Appendixes A and B.
in the thermodynamic limit and the perturbation theory leading to Eq. (23) is no more valid for L larger than some threshold. In Fig. 5 we show the closing with the increasing system size for a gap directly relevant for our analysis. These effects are just the ones we neglect in deriving Eq. (23), and it is because of this neglect that our effective model predicts the expression Eq. (22) correctly, but with incorrect values of B and α. We see this in the lower panel of Fig. 11 where we compare ω peak /J 2 of the Bose-Hubbard model and ω peak /J 2 coming from the simulation of the effective model and we see that they are different (see Appendix B for details).
So, to summarize, the Rabi oscillations are a phenomenon related to the behaviour of just two states in the spectrum and disappear in the thermodynamic limit. They have therefore no relation with the crossover to non ergodicity we have studied in Secs. III, IV and V. The latter phenomenon involves all the spectrum and we have argued it to be independent to the multiplet structure. On the opposite, the perturbative model describing Rabi oscillations strongly relies on this multiplet structure. Its disappearance at large L leads to the destruction of the Rabi oscillations.

VII. CONCLUSION
In conclusion we have studied the ergodicity breaking in the clean Bose-Hubbard model.
We have seen that for small hopping strength J the model breaks ergodicity. We have seen this fact looking at the entanglement entropies of the eigenstates. In the non-ergodic regime the average over the eigenstates of their half-chain entanglement entropy linearly increases with the system size with an angular coefficient smaller than the fully-ergodic one. So the system is non-ergodic. The volume-law behaviour of the entanglement entropy is in strong contrast with the area law behaviour seen in many-body localization. This ergodicity breaking is confirmed by the spectral properties: The average level spacing ratio is not Wigner-Dyson in the same range of J found for the non-ergodicity of the entanglement entropy and this property seems to be stable when size is increased. Most importantly, it appears to be independent of the spectrum being organized in quasidegenerate multiplets at small J and L. The multiplets are doomed to disappear for the system size beyond some threshold, but the ergodicity breaking possibly survives. We have also studied the time behaviour of the Loschmidt echo and we have found that it is in agreement with these findings.
Then we have moved to study the dynamics of the imbalance. We have found that it oscillates with a period exponentially large in the system size when J U , suggesting freezing (and then some sort of space localization) in the infinite-size limit. We have interpreted this phenomenon through a perturbative theory in J. We have constructed an effective XXZ model describing this dynamics. We have found that this effective model shows a doublet breaking the Z 2 and the states in this doublet have a splitting exponentially small in the system size. The imbalance oscillations with a period exponentially large in the system size were actually the Rabi oscillations in this doublet. Therefore these oscillations are a phenomenon involving just a two-dimensional subspace of the Hilbert space, opposite to the ergodicity breaking which involves all the spectrum. The overlap of the initial state with this doublet vanishes in the thermodynamic limit and so also the imbalance oscillations vanish in this limit. We have argued that releasing the approximations leading to the XXZ effective model, the decay of the imbalance-oscillations amplitude with the system size gets stronger. Beyond some value of L they must disappear because the gaps between different multiplets close up. Therefore, in the thermodynamic limit we expect to see no imbalance oscillations and the apparent glassy behaviour we see occurs only at finite system sizes.
Future research will focus on the study of the relevance of our results for the many body localization occurring in the Josephson-Junction chain [25]. Naively one could expect some relation because the Bose-Hubbard model at high energies maps to the Josephson junction chain [51]. Nevertheless, this mapping is valid at equilibrium while here we are inquiring non-equilibrium properties. Another direction of research will be to inquire if it is possible to construct for the clean Bose-Hubbard chain an extensive number of local integrals of motion, as usually occurs for quantum integrable models [52][53][54]. We apply degenerate perturbation theory in J to Eq. (1) when J U 1. Before we go on with this analysis we have to specify some details about the structure of the spectrum. When J = 0, the Hamiltonian Eq. (1) behaves as the operatorV , therefore it shows massively degenerate eigenspaces, as elucidated in [10]. If we switch on a small value of J U , we find a strong mixing inside the eigenspaces, while different eigenspaces interact starting from second perturbative order in J/U . As a result the eigenspaces are transformed into multiplets. The different multiplets are separated by gaps of order U while the levels inside each multiplet have a separation order (J/U ) 2 . We can see the existence of these multiplets for small J by looking at the density of states We show some examples of ρ(E) for different J and different L in Fig. 4. We see that for J ∼ 0.25 it shows a series of δ peaks centered around integer values of the energy E (upper left panel of Fig. 4) to a continuum. Each point corresponds to a multiplet. Increasing the system size L, at some point mixing occurs also across the degenerate subspaces. We see this fact in Fig. 5 where we show the closing with L of the gap between the band at energy around L/2 (this is the energy expectation of |ψ 02 ) and the next one. Nevertheless, we see in this figure that for J small enough the gap is open at the system sizes we can access to, and we can exploit this feature of the spectrum in order to perform a perturbative-analysis interpretation.
We are choosing as initial state of our dynamics the imbalanced state |ψ 02 [Eq. (4), so we are interested in studying the energy splitting inside the unperturbed eigenspace withV -eigenvalue V = L/2.
We can show that the eigenvalues have no corrections at first order in J, as follows. At first order in J the effective Hamiltonian is given by where δ Va,V b means that states |a and |b are unperturbed eigenstates in the sameV -eigenspace. Given a state |b = |n 1 n 2 · · · , we can see that H (1) ef f |b = 0 iff. ∃ j s.t. n j = n j+1 ± 1. It immediately follows that the initial state |ψ 02 belongs to a subspace where H (1) ef f ≡ 0, so that there still is no dynamics at first order in perturbation theory.
Non-trivial dynamics starts to appear at second order in perturbation theory. We apply degenerate perturbation theory and consider coupling only within the samê V -eigenspace. The effective Hamiltonian at second order is then a H (2) So we need to hop two bosons and get back to a state with the same V = L/2. In order to describe the situation, we consider |b = |n 1 n 2 · · · and focus on a nearest neighbour pair (i, j). We can reach the virtual state |c with modified occupation of the pair (i, j) given by n (c) i = n i + 1 and n (c) j = n j − 1 (the discussion of the case n (c) i = n i − 1 and n (c) j = n j + 1 is essentially the same). We now distinguish three cases depending on which sites we act to go from |c to |a : (a) (i, j): In this way n i → n i + 2 and n j → n j − 2. We obey the condition V a = V b iff. n i = n j − 2. The net result is that we moved two bosons from site j to site i.
(b) (j, i): In this way we go back to the initial state. This will produce a nearest-neighbour coupling term diagonal in the number basis.
(c) none of the above: Let's call (l, m) the new pair of nearest neighbour sites on which we acted upon. The term thus produced in H (2) ef f is highly nonlocal. H (2) ef f can move a boson from i to j and from l to m if the overall eigenvalue of V is conserved.
In order to diagonalize the effective Hamiltonian H In an approximate way, we could say that the subspace H (1) ef f = 0 is given by the linear combination of product states |a satisfying the condition H (1) ef f |a = 0. This is not exact, as e.g. |20012 has an overlap with states in the sector H bosonic states |020202 · · · and |202020 · · · . In finitesize systems the true eigenstates will be the even and odd superposition of these states which have the form |ϕ ± = (|↑↓↑↓ · · · ± |↓↑↓↑ · · · ) / √ 2. These states are called cat states, being superpositions of macroscopically ordered classical states, and the elements of each superposition are related with each other by the Z 2 symmetry.
The states |↑↓↑↓ · · · and |↓↑↓↑ · · · are separated by a gap of order J 2 from the rest of states and are degenerate at order 0 int/Ṽ . Applying degenerate perturbation theory int/Ṽ in the subspace generated by these two states, one sees that they are connected at order L. Therefore one sees that the two eigenstates are the cat states and that there is a splitting between them ∆(L) ∼ J 2 (t/Ṽ ) L , which is therefore exponentially small in the system size.
This result is at the roots of the imbalance oscillations. Preparing the system in the imbalanced state |ψ 02 is equivalent to prepare the model in the state |↑↓↑↓ · · · . The evolution amounts to Rabi oscillations of the system between |↑↓↑↓ · · · and |↓↑↓↑ · · · with period ∆(L) ∼ J 2 exp(−| log(t/Ṽ )|L), giving rise to the imbalance oscillation frequency formula Eq. (22) These conclusions are valid for the limitṼ /t 1 but in our caseṼ /2t = 3.5 we have a strictly similar picture. We still have a quasi-degenerate symmetry-breaking doublet at an extremum of the spectrum. These two states are separated from the rest of the spectrum by a gap of order J 2 and they are even and odd superposition of macroscopically ordered symmetry-breaking states, separated by a splitting ∆(L) ∼ J 2 exp(−α XXZ L), for some α XXZ > 0. These two states generate the so-called symmetry-breaking manifold.
The main difference from the limitṼ /t 1 is the following. Before, in the limitṼ /t 1, the state |ψ 02 was an element of the symmetry-breaking manifold, so the dynamics was fully described by this two-level system. Now, |ψ 02 will have an overlap smaller than one with the symmetry-breaking manifold and this overlap tends to 0 in the limit L → ∞ [58]. In this limit the overlap with the rest of the spectrum dominates. Because the only states able to generate a Rabi-oscillation dynamics between two different symmetry sectors are the ones in the symmetry-breaking manifold, we see that the amplitude of the Rabi oscillations goes to 0 in the thermodynamic limit.
We can clearly see this effect directly in the XXZ effective model. We prepare the system in the state |↑↓↑↓ · · · (which corresponds to |ψ 02 in this representation) and we study the equivalent of the imbalance for that model, that's to say the z-staggered magnetization m z S (t) = 1 L L j=1 (−1) j ψ(t)|σ z j |ψ(t) .
First of all we can see qualitatively the existence of the Rabi oscillations in the time traces of Fig. 13 upper panel. These plots are performed for a generic J U because the choice of J is immaterial, being J 2 just a global coefficient of Eq. (A6) giving rise to an overall rescaling of the frequencies. So it is easy to see that our model predicts Rabi oscillations with a frequency proportional to J 2 , as we observed for the Bose-Hubbard Hamiltonian in Fig. 11. Moreover, performing the Fourier transform and finding its main peak as we did above for the imbalance, we find that the Rabi oscillation frequency ω peak exponentially depends on the system size, as we show in the lower panel of Fig. 11. Therefore, our effective model correctly predicts for the frequency of the imbalance oscillations a dependence given by Eq. (22), ω peak = B XXZ J 2 exp(−α XXZ L), and the results of the Bose-Hubbard Hamiltonian obey this rule as soon as J U and we are in the regime where the XXZ effective model is valid.
Looking at the curves in the lower panel of Fig. 11, we see that this prediction is correct qualitatively but not quantitatively. From one side, we see that the quantitative agreement is not extremely bad, because from the numerical fits we get log(B XXZ ) = 0.735 ± 0.020 and α XXZ = 1.106 ± 0.002 which are the same order of magnitude as the log(B) and α found for the imbalance oscillations. From the other side, in getting the effective model Eq. (A6) we have neglected many terms At least for the system sizes we can numerically reach, neglecting these terms has the effect of accelerating the exponential decay of the oscillations frequency, thereby decreasing the splitting in the symmetry-breaking manifold. Moreover, we see that neglecting these terms we have transformed a non-integrable model, like the Bose-Hubbard one, in an integrable one like the XXZ one. We can see this effect especially in the way the amplitude of the Rabi oscillations decays in the two cases, much faster in the Bose-Hubbard model (see lower panel of Fig. 13, where we evaluate the oscillation amplitude as time fluctuations and we define ∆(m z S ) 2 ≡ (m z S ) 2 − (m z S ) (N ∈ N) we can apply an analysis strictly similar to the one performed above and show that also in this case there are Rabi oscillations as in the |ψ 02 case. Applying the perturbation theory at order 2N between nearby sites (similarly to what is done at order 2 in [15,55,56]) we find the effective XXZ model where A 2N and D 2N are positive real numbers. Also here we get a symmetry-breaking two-state manifold separated from the rest of the spectrum by a gap of order J 2 .
[58] (), it is not difficult to show that the overlap of the state |ψ02 and any state of the symmetry-breaking manifold vanishes for L → ∞. Indeed, both of them can be seen as ground states of some non-critical Hamiltonian and then both of them can be written as matrix-product states. The overlap of two matrix-product states in a system with L sites can be written as the trace of some transfer matrix to the power L. Because the transfer matrix has eigenvalues smaller than 1 for the overlap of two different states, the overlap vanishes for L → ∞.