Quantum scars of bosons with correlated hopping

Recent experiments on Rydberg atom arrays have found evidence of anomalously slow thermalization and persistent density oscillations, which have been interpreted as a many-body analog of the phenomenon of quantum scars. Periodic dynamics and atypical scarred eigenstates originate from a “hard” kinetic constraint: the neighboring Rydberg atoms cannot be simultaneously excited. Here we propose a realization of quantum many-body scars in a 1D bosonic lattice model with a “soft” constraint in the form of density-assisted hopping. We discuss the relation of this model to the standard Bose-Hubbard model and possible experimental realizations using ultracold atoms. We find that this model exhibits similar phenomenology to the Rydberg atom chain, including weakly entangled eigenstates at high energy densities and the presence of a large number of exact zero energy states, with distinct algebraic structure. Quantum many-body scarring, a peculiar phenomenon whereby a system thermalizes whilst it keeps returning to its initial state during the time evolution, has recently been observed in experiments on arrays of Rydberg atoms. The authors theoretically investigate the spectral properties of three Hamiltonians using a chain of bosons with density-dependent hopping, providing new insight in the phenomenon of many-body quantum scarring.

S emiclassical studies of chaotic stadium billiards have revealed the existence of remarkable non-chaotic eigenfuctions called "quantum scars" 1 . Scarred eigenfunctions display anomalous enhancement in regions of the billiard that are traversed by one of the periodic orbits in the classical limit when ℏ → 0. It was shown that quantum scars lead to striking experimental signatures in a variety of systems, including microwave cavities 2 , quantum dots 3 , and semiconductor quantum wells 4 .
A recent experiment on a quantum simulator 5 , and subsequent theoretical work 6,7 , have shown that quantum many-body scars can occur in strongly interacting quantum systems. The experiment used a one-dimensional Rydberg atom platform in the regime of the Rydberg blockade 5,8,9 , where nearest-neighbor excitations of the atoms were energetically prohibited. The experiment observed persistent many-body revivals of local observables after a "global quench" 10 from a certain initial state. In contrast, when the experiment was repeated for other initial configurations, drawn from the same type of "infinite" temperature ensemble, the system displayed fast equilibration and no revivals. These observations pointed to a different kind of out-ofequilibrium behavior compared to previous studies of quantum thermalization in various experimental platforms [11][12][13][14][15] .
In both single-particle and many-body quantum scars, the dynamics from certain initial states leads to periodic revivals of the wave function. In the former case, this happens when the particle is prepared in a Gaussian wave packet initialized along a periodic orbit 1 , while in the latter case the revivals can be interpreted as a nearly-free precession of a large emergent su(2) spin degree of freedom 16,17 . Another similarity between single-and many-body quantum scars is the existence of non-ergodic eigenstates. In the single-particle case, such eigenstates are easily identified by their non-uniform probability density that sharply concentrates along classical periodic orbits. In the manybody case, non-ergodic eigenstates are broadly defined as those that violate eigenstate thermalization hypothesis (ETH) 18,19 . Scarred eigenstates violate the ETH in a number of ways: for example, they appear at evenly spaced energies throughout the spectrum 6,20,21 , they have anomalous expectation values of local observables compared to other eigenstates at the same energy density, and their entanglement entropy obeys a sub-volume law scaling 20 .
In recent works, the existence of atypical eigenstates has been taken as a more general definition of quantum many-body scaring. For example, highly excited eigenstates with low entanglement have previously been analytically constructed in the nonintegrable AKLT model 22,23 . A few of such exact eigenstates are now also available for the Rydberg atom chain model 24 . The collection of models that feature atypical eigenstates is rapidly expanding, including perturbations of the Rydberg atom chain 20,25,26 , theories with confinement [27][28][29] , Fermi-Hubbard model beyond one dimension 30,31 , driven systems 32 , quantum spin systems 33,34 , fractional quantum Hall effect in a onedimensional limit 35 , and models with fracton-like dynamics [36][37][38][39] . In a related development, it was proposed that atypical eigenstates of one Hamiltonian can be "embedded" into the spectrum of another, thermalizing Hamiltonian 40 , causing a violation of a "strong" version of the ETH 41,42 . This approach allows to engineer scarred eigenstates in models of topological phases in arbitrary dimensions 43 . From a dynamical point of view, it has been shown that models with scarred dynamics can be systematically constructed by embedding periodic on-site unitary dynamics into a many-body system 44 .
A feature shared by many scarred models is the presence of some form of a kinetic constraint. In the Rydberg atom chain, the constraint results from strong van der Waals forces, which project out the neighboring Rydberg excitations 45 . Such Hilbert spaces occur, for example, in models describing anyon excitations in topological phases of matter [46][47][48][49][50] and in lattice gauge theories [51][52][53] , including the Rydberg atom system 54,55 . Recent works on periodically driven optical lattices have started to explore such physics 56,57 . On the other hand, kinetic constraints have been investigated as a possible pathway to many-body localization without disorder 58 . In classical systems, nonthermalizing behavior without disorder is well known in the context of structural glasses [59][60][61] . The mechanism of this type of behavior is the excluded volume interactions that impose kinetic constraints on the dynamics 62,63 . Similar type of physics has recently been explored in quantum systems where a "quasi manybody localized" behavior was proposed to occur in the absence of disorder [64][65][66][67][68][69][70][71][72][73][74] .
In this paper, we investigate the relation between kinetic constraints, slow dynamics and quantum many-body scars. In contrast to previous work, which focused on models of spins and fermions that are closely related in one dimension due to the Jordan-Wigner mapping, here we study one-dimensional models of bosons with density-assisted hoppings, which realize both "hard" and "soft" kinetic constraints, whilst being non-integrable. Depending on the form of the hopping term, we demonstrate that the models encompass a rich phenomenology, including regimes of fast thermalization, the existence of periodic revivals and many-body scars, as well as the Hilbert space fragmentation that has been found in recent studies of fractonic models [36][37][38][39] . Unlike the experimentally realized Rydberg atom system, we find evidence of many-body scars in a bosonic model without a hard kinetic constraint, i.e., with a fully connected Hilbert space. We identify initial states that give rise to periodic many-body revivals in the quantum dynamics, and we introduce a "cluster approximation" that captures the scarred eigenstates that are responsible for periodic revivals. We discuss possible experimental realizations of these models using ultracold atoms.

Results
Models and their Hilbert spaces. A fundamental ingredient of kinetically constrained models is "correlated hopping": a particle can hop depending on the state of its neighbors. In this paper we consider a system of N p bosons on a one-dimensional lattice with L sites. We consider models where the total filling factor, ν = N p / L, is conserved, and we will mainly present results in the dense regime, ν = 1. We have studied models with ν < 1 and ν > 1, but we found them to be either too constrained or not constrained enough, and therefore less interesting. We emphasize that the bosons in our study are not hard-core, i.e., the occupancy of any lattice site can take any value from 0 to N p .
We study three different models, defined by the Hamiltonians: All three models contain a free-boson hopping term, b y j b jþ1 , which is dressed in various ways by density operators, n j ¼ b y j b j . We will show that the position of the density operator n j completely changes the behavior of these models, ranging from fast thermalization to the breakup of the Hamiltonian into disconnected, exactly solvable sectors. For example, note that H 1 and H 2 are related to each other via free-boson hopping, which can be easily proven using bosonic commutation relations.
We will see below that this innocuous free-boson hopping leads to surprisingly different dynamical properties of the two models. The motivation behind introducing three different models in Eqs. (1)-(3) can be summarized as follows. Hamiltonian H 1 describes a model where a particle cannot hop to the left if that site is not already occupied by at least one particle, and cannot hop to the right if it is the only particle left on its initial site. This introduces constraints to the system. Conversely, there are no such constraints in the case of H 2 . Indeed, the hopping coefficients are only modified in intensity by the particlenumber operator. Hamiltonian H 3 introduces additional constraints compared to H 1 . The number of unoccupied sites and their positions remain constant under the action of this Hamiltonian. This leads to different connectivity of the Hilbert space in each of the models, as we explain in the next Section.
We consider periodic boundary conditions (L + 1 ≡ 1) and set ℏ = J = 1. With periodic boundary conditions, all three Hamiltonians H 1 , H 2 and H 3 have translation symmetry, thus their eigenstates can be labeled by momentum quantum number, k, quantized in units of 2π/L. In addition, H 3 has inversion symmetry. We denote by I = 0 and I = 1 the sectors that are even and odd under inversion, respectively.
Without restrictions on the boson occupancy, the Hilbert space of H 1 , H 2 and H 3 grows very rapidly. For L = N p = 12, the Hilbert space size of the k = 0 sector is 112720 (the largest one we will consider for H 1 and H 2 ). As previously mentioned (see also the next Section), the Hilbert space of H 3 splits into many disconnected components, thus it is possible to consider only one connected component at a time and disregard the unoccupied sites whose positions do not change. This is more relevant when looking at properties such as thermalization, than fixing the filling factor. However, the boundary conditions are in that case no longer periodic, and the system does not have translation symmetry. Considering only a system with the size L/2, filling factor ν = 2, open boundary conditions and minimal number of particles per site equal to 1 is completely equivalent to considering the largest component of the full system which has the size L, filling factor ν = 1, periodic boundary conditions and no restrictions on the occupancies. The Hilbert space size of the symmetric invariant sector of the largest connected component of L = N p = 22 is 176484 and this is the largest sector that we will consider for H 3 .
Graph structure of the models. Since we will be interested in the dynamical properties, it is convenient to first build some intuition about the structure of the Hamiltonians of the three models in Eqs. (1)- (3). A Hamiltonian can be viewed as the adjancency matrix of a graph whose vertices are Fock states of bosons, n 1 ; n 2 ; ; n L j i . If the Hamiltonian induces a transition between two Fock states, the corresponding vertices of the graph are connected by a link. The graphs that show how the configuration space is connected have very different structure for the three Hamiltonians H 1 , H 2 , and H 3 , as can be observed in Fig. 1.
The entire graph of H 2 is well connected and it has the same structure as the graph of the standard Bose-Hubbard model: the particle-number operators in H 2 do not introduce any constraints, but only affect the magnitude of the hopping coefficients. In contrast, the H 1 graph shows several clusters of configurations that are weakly connected to the rest of the graph. "Weakly connected" means that there is a small number of connections leading outside the cluster and that their respective hopping coefficients are smaller in magnitude than those of the surrounding connections within the cluster. A state that is initially located inside a cluster is therefore more likely to stay inside during an initial stage of the time evolution, which increases the probability of revivals and slows down the growth of . d-f same as a, b and c but for L = N p = 6. To avoid clutter, we do not label the vertices in d-f. All graphs are weighted, i.e., the line thickness is proportional to the magnitude of the corresponding hopping coefficient. Several different clusters of configurations are visible in the case of H 1 . The clusters start to form already for L = 3 (for example, the configurations 012-021-003 in a) and become more prominent for L = 6 d. In the case of H 2 , almost all configurations are well connected to the rest of the graph. The graphs for H 3 show that the Hilbert space is highly reducible: its graph splits into many disconnected components. entanglement entropy. We will provide a more quantitative description and examples that illustrate this in Section "Quantum scars in H 1 and H 3 models". Finally, the graph of H 3 , due to even stronger constraints, is actually disconnected, which is an example of Hilbert space fragmentation that was previously shown to cause non-ergodic behavior in fracton-like models 37,38 . This predicts that thermalization and dynamics in the three models will be very different, which we will confirm in the following Section. However, we note that the number of connections and the topology of the graph is not the only relevant factor for the dynamics. The magnitude of the hopping coefficients between different configurations is also important (Supplementary Note 1).
We note that the relation between H 1 and H 3 is reminiscent of the relation between the quantum East model 75 and the "PXP" model describing the atoms in the Rydberg blockade regime 6,20,45 . Like H 3 , the PXP model is doubly constrained and inversion symmetric, while H 1 and the quantum East model are asymmetric versions of those two models with only a single constraint. The graph of the quantum East model is similar to that of H 1 , in that it contains bottlenecks which slow down the growth of entanglement entropy 75 .
Dynamics and entanglement properties. We now investigate the phenomenology of the models introduced in Eqs. (1)-(3). We use exact diagonalization to obtain the complete set of energy eigenvalues and eigenvectors, from which we evaluate the level statistics and the distribution of entanglement entropies for the three models. Furthermore, we probe dynamical properties of the models by studying a global quench, simulated via Krylov iteration.
The energy level statistics is a standard test for thermalization of models that cannot be solved exactly. A convenient way to probe the level statistics is to examine the probability distribution P(r) 76 of ratios between consecutive energy gaps s n = E n+1 − E n , r ¼ minðs n ; s nþ1 Þ maxðs n ; s nþ1 Þ : ð5Þ The advantage of studying P(r), instead of P(s n ), is that there is no need to perform the spectrum unfolding procedure-see ref. 77 . For standard random matrix theory ensembles, both P(r) and the mean 〈r〉 are well known 78 . When computing the same quantities in a microscopic physical model, it is crucial to resolve all the symmetries of the model.
The probability distribution P(r) of the ratios of two consecutive energy gaps is shown in Fig. 2a-c for the three Hamiltonians H 1 , H 2 , and H 3 respectively, and two momentum or inversion sectors. In all three cases, the energy levels repel, i.e., the distribution tends to zero as r → 0. For H 2 , the distribution is particularly close to the Wigner-Dyson (non-integrable) line. For H 1 , the distribution is also consistent with Wigner-Dyson when we restrict to the middle 1/3 of the spectrum (and after removing special states with E = 0). We exclude the edges of the spectrum because they contain degeneracies which are not symmetryrelated. However, such states do not appear to have a major effect on the level statistics distribution, which is still closer to the Wigner-Dyson than the Poisson distribution even if they are included. The level statistics of H 3 within the largest connected component of the Hilbert space is shown in Fig. 2c and is also consistent with the Wigner-Dyson distribution without restricting the spectrum. However, we will demonstrate below that the dynamics in some smaller connected components of H 3 can be exactly solved.
As a complementary diagnostic of thermalization, we next compute the entanglement entropy of all eigenstates. We divide the lattice into two sublattices, A and B, of lengths L A and . The black dashed line shows the Poisson distribution, which corresponds to the integrable case, while the red dashed line is the distribution of the Gaussian orthogonal ensemble, which corresponds to the thermalizing case. d-f Entanglement entropies S L/2 of all eigenstates plotted as a function of the eigenstate energy per particle, For a given pure state ψ j i, the entanglement entropy is defined as where ρ A ¼ tr B ψ j i ψ h j is the reduced density matrix of the subsystem A. The scatter plots, showing entanglement entropy of all eigenstates E n j i as a function of their energy E n , are displayed in Fig. 2d-f. Here we take into account the translation symmetry of the system and work in the momentum sector k = 0 for H 1 and H 2 , and consider only the largest connected component and the inversion sector I = 0 for H 3 . The results for other sectors are qualitatively similar.
Entanglement entropy distribution in Fig. 2d, e reveals a striking difference between the Hamiltonians H 1 and H 2 , even though they only differ by a free-boson hopping term, Eq. (4). The model H 1 is constrained, which leads to a large spread of the entropy distribution and many low-entropy eigenstates including in the bulk of the spectrum. From this perspective, H 1 is reminiscent of PXP model 20,25 . By contrast, H 2 has no such constraints and in this case the entanglement entropy is approximately a smooth function of the eigenstate energy. The Hamiltonian H 3 is doubly constrained, and this is reflected in its entanglement distribution, which also shows a large spread and several disconnected bands, reminiscent of an integrable system like the XY model 79 .
Global quenches. The constraints in the models in Eqs. (1), (2), and (3) have significant effects on the dynamics governed by these Hamiltonians. We probe the dynamics by performing a global quench on the system. We assume the system is isolated and prepared in one of the Fock states, ψ 0 , at time t = 0. We restrict to ψ 0 being product states which are not necessarily translation-invariant, as such states are easier to prepare in experiment. However, our results remain qualitatively the same if we consider translation-invariant ψ 0 . After preparing the system in the state ψ 0 , which is not an eigenstate of the Hamiltonian, the system is let to evolve under unitary dynamics, where H is one of the Hamiltonians of interest. From the timeevolved state, we evaluate the quantum fidelity, i.e., the probability for the wave function to return to the initial state. In a general many-body system, fidelity is expected to decay as FðtÞ $ expðÀLðJtÞ 2 Þ. It thus becomes exponentially suppressed in the system size for any fixed time t * , i.e., Fðt Ã Þ $ expðÀcLÞ, where c is a constant. In scarred models, such as the Rydberg atom chain, fidelity at the first revival peak occurring at a time T still decays exponentially, but exponentially slower, i.e., FðTÞ $ expðÀc 0 LÞ, with c 0 ( c. In ref. 20 , for a finite system with L ≲ 32 atoms, the fidelity at the first revival can be as high as~70%, and several additional peaks at times nT are also clearly visible. We first consider the Hamiltonian H 1 . Several configurations exhibit periodic revivals of the fidelity F(t), which can in some cases be higher than 90%. Most of these configurations involve a very dense cluster of bosons such as :::0N10::: j i . In contrast, a completely uniform configuration :::111::: j i thermalizes very quickly. Here we focus on periodically-reviving configurations with density being as uniform as possible. One family of such reviving configurations involves n unit cells made of three lattice sites: Time evolution of the fidelity for the initial state ð210Þ n j i for different system sizes L = 3n is shown in Fig. 3a. The initial state is assumed to be the product state, e.g., ψ 0 ¼ 210 j i for L = 3. The frequency of the revivals in Fig. 3 is approximately the same for all system sizes. We emphasize that similar results are obtained for a translation-symmetric initial state, e.g., Þ . Both cases converge in the large system limit, and the differences are only significant for L = 3 when the revival frequency of the initial state with transition symmetry differs from the frequencies of other system sizes.
In Fig. 3b we compare the fidelity for the initial state in Eq. (9) when it is evolved by all three Hamiltonians in Eqs. (1)-(3). The initial state is fixed to be ð210Þ 5 . We observe that the dynamics with H 3 has very prominent revivals; in fact as we will later show, these revivals are perfect and their period is approximately twice the revival period for H 1 . In contrast, for H 2 the fidelity quickly drops to zero without any subsequent revivals.
Finally, in Fig. 3c we plot the time evolution of entanglement entropy. As expected from the fast decay of the fidelity, the entropy for H 2 rapidly saturates to its maximal value. Moreover, as expected from the perfect revivals in H 3 , the entropy in that case oscillates around a constant value close to zero. For H 1 , we observe a relatively slow growth of entropy, with oscillations superposed on top of that growth, again similar to PXP model 6 . For the initial state that is not translation-invariant, it is important how we cut the system, e.g., :::210j210::: j iversus :::2102j10::: j i . In the first case, the entanglement entropy remains zero for H 3 because no particle can hop from one subsystem to the other, while in the second case the entropy oscillates around a constant value, which is the case in Fig. 3c. In Fig. 4 we show the H 1 evolution of two local observables, density correlations between two adjacent sites 〈n 1 n 2 (t)〉 and density on the first site 〈n 1 (t)〉, starting from the initial state ð210Þ n j i. Unlike fidelity and entanglement entropy, these observables can be easily measured in experiment. Both observables robustly oscillate with approximately the same frequency as the fidelity. The heights of the first few revival peaks are approximately converged for the system sizes ranging from L = 6 to L = 15, which suggests that revivals in such local observables can be observed in the thermodynamic limit. In the following Section, we will show that the oscillations observed in the dynamics from ð210Þ n j istate in Eq. (9) and their frequency can be explained using a tractable model that involves only a small subset of all configurations in the Hilbert space, thus providing a realization of quantum scars in a correlated bosonic system. Our starting point will be the model H 3 , whose graph explicitly separates into disconnected subsets which makes the toy model exact, hence we can analytically calculate the revival frequency. Based on these results, we then introduce an approximation scheme that describes the dynamics from the same initial state under the H 1 Hamiltonian.
Quantum scars in H 1 and H 3 models. The quench dynamics of fidelity and entanglement entropy in Fig. 3 suggest that H 1 and H 3 models are candidate hosts for many-body scarred eigenstates that can be probed by initializing the system in product states ð210Þ n j i. We now analyze the structure of these states using our approach called "cluster approximation" that is introduced in detail in Methods.
The dynamics of H 3 within the sector containing the state ð210Þ n j i can be solved exactly, as shown in Methods. The connected component of the state ð210Þ n j iconsists of all possible combinations of patterns 210 and 120. This means that triplets of sites evolve independently, and dynamically the system behaves as a collection of independent two level systems (spins-1/2). From this observation, it can be shown that revivals will be perfect with a period T 3 = π/2. The same period is obtained for initial product state ð210Þ n j i and its translation-invariant version; if the initial state is both translation-invariant and inversion-symmetric, the period is doubled.
In contrast to the free dynamics in H 3 , the H 1 model exhibits decaying revivals and does not admit an exact description. In order to approximate the quench dynamics and scarred eigenstates in H 1 , we project the Hamiltonian to smaller subspaces of the full Hilbert space. These subspaces contain clusters of states which are poorly connected to the rest of the Hilbert space and thereby cause dynamical bottlenecks. As explained in Methods, the clusters can be progressively expanded to yield an increasingly accurate description of the dynamics from a given initial state.
For our initial state ð210Þ n j i, the minimal cluster is defined as one that contains all the states given by tensor products of 210, 120 and 300 patterns. Similar to the H 3 case, within this approximation, triplets of sites again evolve independently, and the dimension of the reduced Hilbert space is dim H c ¼ 3 L=3 . The time-evolved state within the cluster is given by where the dots denote other configurations. The fidelity is As in the case of H 3 , this result is also valid for the translationinvariant initial state. We see that the period of revivals is T 1 = π/ 4, which is the same as for H 3 with a translation and inversion symmetric initial state. The result of the cluster approximation is compared against the exact result for system size L = 15 in Fig. 5. The frequency of the fidelity revival, shown by the blue line in Fig. 5a, is accurately reproduced in this approximation, however the approximation does not capture the reduction in the magnitude of F(t). Similarly, the dynamics of entanglement entropy, blue line in Fig. 5b, is only captured at very short times. In particular, we observe that the maximum entanglement within the cluster remains bounded even at long times t~10, while the exact entropy continues to increase and reaches values that are several times larger.
To obtain a more accurate approximation, we can expand the minimal cluster with several neighboring configurations on the graph. We define the extended cluster as a set of all states which can be obtained using tensor products of the configurations 210, 120, 300, and 111. The enlarged cluster clearly contains the minimal cluster studied above, but it also includes additional configurations, resulting in a much better prediction for the first revival peak height, while still allowing for analytical treatment.  The dimension of the extended cluster grows as dim Hc ¼ 4 L=3 , and is thus exponentially larger than the minimal cluster approximation. Nevertheless, the extended cluster dimension is still exponentially smaller compared to the full Hilbert space, and within this approximation it is possible to numerically simulate the dynamics of larger systems, L ≲ 30-see Fig. 6a. The revivals are no longer perfect, while their frequency is independent of the system size and closer to the frequency of revivals for the full Hilbert space compared to to the minimal cluster approximation in Fig. 5. The overlap between the eigenstates of the Hamiltonian H 1 reduced to both the minimal and extended cluster and the state ð210Þ 8 is given in Fig. 6b. The eigenstates that correspond to the minimal cluster approximately survive in the extended cluster, where they form a band with the highest overlap.
For the initial product state (210) n , it is possible to analytically obtain the fidelity within the improved approximation for arbitrary system size.  Fig. 6a. It was also found to be a very good approximation for the translation-invariant initial state when L ≥ 9 (data not shown). Figure 7a shows that the logarithm of the fidelity per site, log ðFðTÞÞ=L, at the first peak, saturates at a finite value for large L. In the improved cluster approximation, the first peak height decays as e −0.04L (Supplementary Note 2). For a completely random state, the fidelity would be F $ 1=dim H . In the case ν = 1 and large L, the Hilbert space dimension grows with the system size as This back-of-the-envelope estimate suggests the fidelity of a random state is F~e −1.39L , which decays considerably faster than the first peak height in Fig. 7. The improved cluster approximation correctly reproduces the short-time dynamics, including the first revival peak, and sets a lower bound for the first peak height see Figs. 5 and 7b.
The evolution of the entanglement entropy for the extended cluster approximation is shown in Fig. 5b. Inside the cluster, entropy remains approximately constant with periodic oscillations that have the same frequency as the wave function revivals. Any further growth of the entanglement entropy can be attributed to the spreading of the wave function outside the cluster. To illustrate the "leakage" of the wave function outside the cluster, in Fig. 8 we compute the time evolution of the overlap with a cluster, i.e., the probability to remain inside a cluster at time t, We consider several initial configurations that lie inside or outside the cluster. The configurations initially inside the cluster mostly stay there, and the configuration ð210Þ 4 that has the highest revivals also has the highest overlap. Similarly, configurations initially outside the cluster continue to have negligible overlaps. The overlap starting from the configuration ð210Þ 4 Fig. 6 Cluster approximations. a Fidelity F(t), for the Hamiltonian H 1 and initial states ð210Þ n j i , in the extended cluster approximation for various system sizes. b Eigenstate overlap with the initial state jð210Þ 8 i plotted on a log scale, for both cluster approximations. In the case of degenerate eigenstates the sum of their overlaps is shown. Fig. 7 First peak height. a Logarithm (base 10) of the first revival peak divided by the system size, log ðFðTÞÞ=L, seems to saturate at a finite value in the thermodynamic limit. b Comparison of the logarithm of the first revival peak height for the full dynamics and the improved cluster approximation. The approximation serves as a lower bound. We now summarize the relation between H 3 and H 1 from the point of view of the cluster approximation. For the initial state ð210Þ n j i, the two models yield similar dynamics, compare Eq. (23) and Eq. (12). The only difference is that the revival frequency is doubled in the latter case, which can be easily explained by the symmetry of the initial state and that of the Hamiltonian. Hamiltonian H 3 is inversion-symmetric. If the initial state is also chosen to be inversion-symmetric, the frequency of the revivals doubles. The period is then T inv 3 ¼ π=4, which is equal to the period of revivals T 1 of H 1 in the cluster approximation. This is also proven analytically in Methods, see Eq. (27). For comparison, the revival period for the full Hilbert space is approximately 0.77, which is slightly less than π/4 ≈ 0.79. The Hamiltonian H 1 is not inversion-symmetric, so the frequency does not double for an inversion-symmetric initial state, but the revivals are lower in that case. This shows that it is important for the symmetry of the initial state to match the symmetry of the Hamiltonian.
Finally, the eigenstates of H 1 , projected to the subspace of the minimal cluster approximation, form several degenerate bands whose eigenenergies are equally spaced in integer multiples of 4. Interestingly, some of these eigenstates approximately survive in the full H 1 model, and they are precisely the eigenstates that have the highest overlap with the configurations ð210Þ n j i, see Fig. 9a. In small system sizes, such as L = 6, the surviving eigenstates are also the lowest entropy eigenstates in the middle of the spectrum, which is reminiscent of quantum scars in the PXP model 20 . In larger systems, e.g., L = 12, the surviving eigenstates are slightly lower in entropy than their neighbors, but are far from being the lowest entropy eigenstates, as can be seen in Fig. 9. The lowest entropy eigenstates have high overlaps with other configurations, such as ð3100Þ 3 , as shown in Fig. 9b, c. In the case of ð210Þ n j i, the eigenstates surviving in the full system belong to every other band of eigenstates in the cluster approximation and the number of the surviving eigenstates is n + 1. For even system sizes this counting includes a zero-energy eigenstate. In Methods we discuss in more details the generalization of the cluster approximations to the states of the form ðN10:::0Þ n j i , which were also found to have robust revivals and high overlaps with some low-entropy eigenstates.

Discussion
In this paper, we have introduced three models of bosons with "soft" kinetic constraints, i.e., density-dependent hopping. We have demonstrated that some of these models exhibit similar phenomenology to other realizations of quantum many-body scars, for example the Rydberg atom system 5 . We have studied quantum dynamics of these systems by performing global quenches from tensor-product initial states. We have shown that both the connectivity of the Hilbert space and the relative magnitude of the hopping coefficients have dramatic effects on the dynamics. For certain initial configurations, the constraints can lead to slow thermalization and revivals in the quantum fidelity. The revival frequency can be predicted by considering an exponentially reduced subset of the Hilbert space. For a family of initial configurations of the form ð210Þ n j i, we have derived analytical expressions for the evolution of quantum fidelity within this approximation, which accurately capture the revival frequency obtained from exact numerical data. One notable difference between scarred dynamics in the present bosonic models and the PXP model is that the revivals exist in the absence of a hard kinetic constraint, i.e., in the fully connected Hilbert space. Our cluster approximation also explains the structure of some low-entropy eigenstates in the middle of the many-body spectrum. In addition, we have calculated the evolution of two local observables which are experimentally measurable, density correlations between two neighboring sites and density on a single site, and both of them show robust oscillations over a range of system sizes. We have also shown that the introduced models contain additional special properties, like the exponentially large zeroenergy degeneracy which is related to the bipartite structure of the model. We now comment on the possible experimental realizations of the models we studied. The implementation of a correlated hopping term (n k b y i b j ) in optical lattices has attracted lot of attention due to a possible onset of quantum phases related to high-Tc superconductivity 80 . An early theoretical proposal exploits asymmetric interactions between the two atomic states in the presence of a state-dependent optical lattice 80 . As a result, the obtained effective model corresponds to the inversion-symmetric form of H 1 . In addition, the same term has been found to feature as a higher-order correction of the standard Bose-Hubbard model [81][82][83][84] . Although in this case the term typically represents a modification of the regular hopping term of the order of several percent, its contribution was directly measured 85,86 . More recently, the set of quantum models accessible in cold-atom experiments has been enriched through the technique of Floquet engineering 87 . As a notable example, a suitable driving scheme can renormalize or fully suppress the bare tunneling rate 88 . On top of that, by modulating local interactions an effective model with the density-dependent tunneling term has been engineered 89 . For the models considered in this paper the most promising is a more recent driving scheme exploiting a double modulation of a local potential and on-site interactions 90 . Related . c Entanglement entropy, eigenstates which have the highest overlap with some product states are marked in different colors. sophisticated driving schemes have already enabled a realization of dynamical gauge fields 56,57,91 where both the amplitude and the phase of the effective tunneling are density-dependent. Although these experimental proposals explain how to realize some of the correlated hopping terms present in our models using ultracold atoms in optical lattices, finding a scheme that exactly realizes our models requires further study. We emphasize that other models which would exhibit non-ergodic dynamics and scarred eigenstates as a result of the same mechanism that was explained in this work could be built, for example a linear combination of H 1 and H 2 .
Note added: During the completion of this work, we became aware of ref. 92 which identified non-thermal eigenstates and slow dynamics in the quantum East model. Moreover, a recent study 93 proposed a Floquet scheme for a bosonic model with densityassisted hopping, finding signatures of quantum manybody scars.

Methods
In order to more efficiently describe the dynamics of our models, we introduce a method-"cluster approximation", that is based on Hilbert space truncation inspired by the bipartite graph structure of H 1 . Before providing details about the cluster approximation for H 1 and its generalizations, we present an exact solution for the perfect revivals in H 3 model, which serves as a motivation for the more complicated case of H 1 .
Bipartite lattice and zero modes. The graph of H 1 is bipartite, i.e. all the basis configurations can be divided into two disjoint sets, and the action of the Hamiltonian connects configurations in one set only to the configurations in the other and vice-versa (the Hamiltonian is off-diagonal). One way to sort configurations into these two sets is by parity of the quantity where C = 0 if L is even and C = 1 if L is odd. We define n even and n odd as the total numbers of particles at even and odd sites, respectively, n even ¼ where L 1 = L 2 = L/2 if L is even, and L 1 = (L − 1)/2, L 2 = (L + 1)/2 if L is odd. If only nearest-neighbor hoppings are allowed and if no two odd sites are coupled (if the system has open boundary conditions for any L or periodic boundary conditions for L-even), each hopping either increases n even by one and decreases n odd by one, or vice-versa. This means that each hopping can change Δ a only by ±1. In special cases, like H 1 at filling factor ν = 1, it is also possible to define quantities like Δ a for odd system sizes and periodic boundary conditions. This is a consequence of the constraints imposed by H 1 , i.e., the fact that a particle cannot hop to an empty site to its left (Supplementary Note 3). Note that H 2 in the same geometry is not bipartite.
Another way to sort configurations into two sets is by parity of the distance from the configuration 111:::111 j i , which we define as d a ¼ min n f 111:::111 In this case, the two sets are the configurations with even and with odd distances d a . One hopping can change d a only by ±1 or 0. Changes by other values are not possible by definition if the Hamiltonian is Hermitian (all hoppings are reversible). Both d a and Δ a have the same parity, thus d a must always change after one hopping in even system sizes or in systems with open boundary conditions. As a consequence, d a cannot change by 0 if Δ a can only change by ±1. The graphs of bipartite systems do not contain any loops of odd dimension (triangles, pentagons, heptagons and so on). Moreover, the energy spectra of bipartite systems are symmetric around zero. Their Hamiltonians anticommute with the operator ðÀ1Þ Δ a . The presence of such an operator in a bipartite lattice leads to exact zero energy states in the spectrum 94,95 . It can be shown that the exponentially growing number of zero modes of H 1 is related to the difference between the numbers of elements in the two sets of its bipartite graph (Supplementary Note 4). Additionally, the algebraic structure of zero energy eigenstates can be explained by the structure of the graphsuch eigenstates can be constructed as superpositions of configurations from only one of the sets. Similar properties are found for H 2 for even L, as its graph is also bipartite in that case. The properties of the zero-energy manifold are discussed in more detail in Supplementary Note 4.
Perfect revivals in the H 3 model. We start with a warmup calculation for H 3 acting on L = 3 sites. The connected subspace of 210 contains only two configurations, 120 and 210. The Hamiltonian reduced to this subspace is where the basis vectors are The eigenvalues of H 0 3 are E 1 = −2 and E 2 = 2. The initial state ψ 1 ðt ¼ 0Þ ¼ 210 j i evolves as and the state ψ 2 ðt ¼ 0Þ ¼ 120 j i evolves as ψ 2 ðtÞ ¼ Ài sinð2tÞ 210 j iþ cosð2tÞ 120 j i: Previous results can be straightforwardly generalized to larger systems. Let the length of the system be L = 3n for simplicity. The connected component of the state ð210Þ n j iconsists only of combinations of patterns 210 and 120, which means that triplets of sites evolve independently. From Eq. (20), the initial state ψ n ðt ¼ 0Þ ¼ ð210Þ n j ievolves as ψ L¼3n ðtÞ ¼ cos n ð2tÞ ð210Þ n j i þ ðÀiÞ n sin n ð2tÞ ð120Þ n j iþ ::: where ". . . " denotes contributions of the basis configurations other than ð210Þ n j i or ð120Þ n j i . The fidelity is It follows that the revivals are perfect, with a period T 3 = π/2. This result is also valid for the translation-invariant initial state ð210Þ n j i T , as the connected subspaces of 210, 021 and 102 do not overlap and therefore evolve independently. However, an initial state that is both translation symmetric and inversion symmetric has different dynamics. The inverse of the configuration ð210Þ n j iis the configuration ð012Þ n j i , which is a translation of the state ð120Þ n j ithat belongs to the connected subspace of ð210Þ n j i . The initial state ψ inv n ðt ¼ 0Þ evolves as ψ inv n ðtÞ ¼ cos n 2t þ ðÀiÞ n sin n 2t ð Þ ψ inv n ðt ¼ 0Þ þ ::: ð26Þ and the fidelity is F inv n ðtÞ ¼ jhψ inv n ð0Þjψ inv n ðtÞij 2 ¼ jcos n 2t þ ðÀiÞ n sin n 2tj 2 : The frequency of the revivals is now doubled, so the period is T inv 3 ¼ π=4.
Cluster approximations for the H 1 model. Here we introduce a scheme for approximating the dynamics from initial states (210) n in the H 1 model. As can be observed in Fig. 3, the revival periods are approximately the same for different system sizes. We first focus on the non-trivial case L = 6. Figure 10 shows part of the graph that contains the initial state, 210210 j i . Configurations labeled inside the ellipses denote representatives of an orbit of translation symmetry, i.e., the configurations are translation-invariant such as the one in Eq. (24).
The minimal subcluster of the graph is highlighted in blue color in Fig. 10. This cluster is indeed weakly connected to the rest of the configuration space, as it has only 3 connections that lead outside this cluster (dashed lines) and their hopping coefficients are slightly lower in magnitude than those inside the cluster, meaning that the probability is higher to stay inside the cluster than to leave. The hopping coefficients leading outside are not significantly smaller than the coefficients staying inside, but in combination with the relatively small number of connections this has significant effects on the dynamics. This effect is even more pronounced when the difference in magnitudes is further increased by squaring the particlenumber operators (see Supplementary Note 1).
The minimal cluster from Fig. 10 contains all the states given by tensor products of 210, 120 and 300 configurations. The set of configurations belonging to this cluster could have been chosen differently, but this particular choice has at least two advantages. Firstly, inside this cluster, the evolution of the configuration 210210 j ican be thought of as two subsystems 210 evolving separately. The evolution of all such configurations at different system sizes can be reduced to the evolution of L = 3 subsystems 210, similar to the case of H 3 in the connected subspace of (210) n . Secondly, this definition allows easy generalization to different system sizes L = 3n with initial states (210) n . We would like to emphasize that the cluster was not chosen arbitrarily. The calculations of the probability density distribution starting from the initial configuration 210210 j iand evolving with H 1 have shown that the probability density stays high in this region of the Hilbert space as long as the revivals in fidelity are visible. The configurations important for the dynamics were then identified by analyzing the structure of the graph around the initial configuration.
As an example, consider system size L = 3. The reduced Hilbert space of the cluster H c is spanned by the (non-translation-invariant) configurations The Hamiltonian reduced to this subspace is and its eigenvalues are By generalizing this result to larger systems, it is easy to prove Eqs. (10) and (11). The minimal clusters can be expanded by adding several neighboring configurations. For similar reasons as in the case of minimal clusters, the extended clusters are defined as sets of all states which can be obtained using tensor products of the configurations 210, 120, 300 and 111. In the case of L = 6, the enlarged cluster can be observed in Fig. 10. It contains the minimal cluster studied previously, but it also includes additional configurations shown in green ellipses. Again, the approximation could be improved by including more configurations, but this particular choice is well suited for analytical treatment (Supplementary Note 2) and, as shown above, it gives a good prediction for the first revival peak height.
Generalization to other clusters. Building on the previous observation that some of the low-entropy eigenstates have large weight on ð3100Þ 3 product state, we have investigated periodic revivals from such a larger class of initial states. We find that robust revivals are associated with initial product states of the form ððN À 1Þ1 0::: where N is the length of the unit cell (L = Nn). For example, some of these configurations are ð3100Þ n j i , ð41000Þ n j iand ð510000Þ n j i . Combinations of those patterns such as 310041000 j ialso exhibit similar properties, but we will restrict ourselves to the simpler former cases.
We can construct a generalization of the cluster approximation for configurations of the form in Eq. (31). As in the case of ð210Þ n j i , the dynamics inside one unit cell explains the dynamics of the full system. The generalized clusters can be chosen in such a way that their Hilbert spaces are spanned by N configurations where i takes values 1, 2, …N. If we consider only one unit cell (n = 1), the graph that connects these configurations has a linear structure without any loops, i.e., each configuration i j i is solely connected to the configurations i ± 1 j i, except the two configurations at the edges, 1 j i and N j i, which are only connected to 2 j i and N À 1 j i , respectively. The projection of the Hamiltonian H 1 to this cluster, which we denote by H c 1 , has a very simple structure: it has the form of a tight-binding chain with the only nonzero matrix elements on the upper and lower diagonals: The dynamics within a single unit cell under H c 1 corresponds to density fluctuations between the first and the second site. Following the same procedure as previously, we can now diagonalize H c 1 and compute the fidelity time series for the initial configuration ðN À 1Þ10:::0 j i . This result can be directly generalized to configurations of the form ððN À 1Þ10:::0Þ n j i . The derivation is valid for both translation-invariant and non-translation-invariant initial configurations, as the cluster in Eq. (32) is disconnected from its translated copies. We stress that this disconnection, namely the absence of a hopping term between 1ðN À 1Þ0:::0 j i and 0N0:::0 j i , is a consequence of the constraints imposed by H 1 and it would not hold for H 2 . In this way, we have calculated the time evolution of the fidelity starting from the configurations ð3100Þ n j i(for n = 1, 2, 3, 4), ð41000Þ n j i(n = 1, 2, 3) and ð510000Þ n j i (n = 1, 2), and compared it with the exact numerical results for the full H 1 . The cluster approximation captures both the revival frequency and the height of the first peak. Similar to the ð210Þ n j icase, the approximation can be improved by adding further configurations to the clusters. Moreover, if we want to consider translation-invariant initial states, we can follow the same recipe for ð210Þ n j iby summing translated patterns with the required phase factors given in terms of momenta in multiples of 2π/N. We have checked that revivals appear in these momentum sectors, with roughly the same frequency.
We note that the configurations with larger units cells thermalize more quickly on shorter timescales, but slower at long times. Initially, the states starting from configurations with smaller N have lower entanglement entropies than those with larger N. The Hilbert spaces of large N unit cells are larger, so the entanglement entropy starting from these configurations rapidly grows to the maximal value for that unit cell. However, the only way for the wave function to spread through the entire Hilbert space is that a unit cell reaches a state close to 111. . . 111, so that particles can hop to the other unit cells. This is less likely for large N, and therefore such configurations need long times to fully thermalize. As a result, smaller N entanglement entropies grow faster and after long enough time they overtake those for larger N. For example, in the case of L = 12 and translation-invariant initial states, (210) 4 overtakes (3100) 3 and (510000) 2 around t~2, and (3100) 3 overtakes (510000) 2 around t~80 ( Supplementary Fig. 3).
Finally, we note that non-thermal behavior reminiscent of the one studied here was previously observed in an unconstrained Bose-Hubbard model, for example in the context of "arrested expansion" 96,97 and quenches from superfluid to Mott insulator phase 98,99 . In these cases, the main ingredient is the strong on-site interaction, which causes the energy spectrum to split into several bands. Due to the large energy differences between bands, the dynamics of an initial state from a particular band is at first limited only to the eigenstates that belong to the same band. Additionally, these energy bands are approximately equally spaced, which can lead to revivals in fidelity if several bands are populated. In contrast, our models do not feature on-site interaction, and the mechanism which slows down the spread of the wave function is correlated hopping, which suppresses connections between certain configurations and modifies the hopping amplitudes between others, thus creating bottlenecks that separate different clusters of states.

Data availability
The data that support the plots within this paper and other findings of this study are available at https://doi.org/10.5518/810 100 . Fig. 10 Minimal and extended clusters. Hamiltonian H 1 and system size L = N p = 6. Configurations labeled inside the ellipses are representatives of an orbit of translation symmetry. The minimal cluster is defined by the blue configurations, while green configurations represent the additional components of the extended cluster. Gray arrows connect to configurations outside the extended cluster. The numbers bellow the graph show the distance d a from the configuration 111111 evaluated using Eq. (17).