Quantum scars of bosons with correlated hopping

Recent experiments have shown that preparing an array of Rydberg atoms in a certain initial state can lead to unusually slow thermalization and persistent density oscillations [Bernien et al., Nature 551, 579 (2017)]. This type of non-ergodic behavior has been attributed to the existence of"quantum many-body scars", i.e., atypical eigenstates of a system that have high overlaps with a small subset of vectors in the Hilbert space. Periodic dynamics and many-body scars are believed to originate from a"hard"kinetic constraint: due to strong interactions, no two neighbouring Rydberg atoms are both allowed to be excited. Here we propose a realization of quantum many-body scars in a 1D bosonic lattice model with a"soft"constraint: there are no restrictions on the allowed boson states, but the amplitude of a hop depends on the occupancy of the hopping site. 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. We discuss the relation of this model to the standard Bose-Hubbard model and possible experimental realizations using ultracold atoms.


I. INTRODUCTION
Semiclassical 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 5,8,9 in the regime of the Rydberg blockade, where nearest-neighbour 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-of-equilibrium behaviour 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. In the latter case, the Rydberg atoms are prepared in a Néel product state, where every second atom is excited. The collective dynamics from such an initial state was shown to be effectively "semiclassical", 7 i.e., it is effectively restricted to a subset of the Hilbert space and can be interpreted as a nearly-free precession of an emergent spin-1/2 degree of freedom. 16 Another similarity between single-and many-body quantum scars is the presence of non-ergodic eigenstates. As mentioned above, 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 many-body case, non-ergodic eigenstates are broadly defined as those that violate Eigenstate Thermalization Hypothesis (ETH). 17,18 In the Rydberg atom chain describing Ref. 5, such eigenstates appear at evenly spaced energies throughout the spectrum of the system. 6,19,20 Scarred eigenstates violate the ETH in a number of ways: for example, 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, 19 i.e., scarred eigenstates are weakly entangled compared to random vectors in the many-body Hilbert space.
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 non-integrable AKLT model. 21,22 A few of such exact eigenstates are now also available for the Rydberg atom chain model. 23 The collection of models that feature atypical eigen-states is rapidly expanding, including perturbations of the Rydberg atom chain, 19,24,25 theories with confinement, [26][27][28] Fermi-Hubbard model beyond one dimension, 29,30 driven systems, 31 fractional quantum Hall effect in a one-dimensional limit, 32 and models with fractonlike dynamics. [33][34][35][36] In a related development, it was proposed that atypical eigenstates of one Hamiltonian can be "embedded" into the spectrum of another, thermalizing Hamiltonian, 37 causing a violation of a "strong" version of the ETH. 38,39 This approach allows to engineer scarred eigenstates in models of topological phases in arbitrary dimensions. 40 From a dynamical point of view, it has been shown that models with scarred dynamics can be systematically constructed by embedding periodic onsite unitary dynamics into a many-body system. 41 A feature shared by the models currently believed to support quantum many-body scars is the presence of some form of a kinetic constraint. In the Rydberg atom chain, the constraint results from strong van der Waals forces between neighbouring excitations. Formally, this means that all configurations with adjacent Rydberg excitations are removed from the low-energy subspace. The resulting effective Hilbert space, while still exponentially large in the number of atoms, is no longer a tensor product of two-level systems. 42 Such Hilbert spaces occur, for example, in models describing anyon excitations in topological phases of matter [43][44][45] due to a global constraint imposed by an emergent gauge field. Recent numerical investigations of the ETH in such models have argued that they generically thermalize. 46,47 Kinetic constraints emerge naturally in models of lattice gauge theories, [48][49][50] including the Rydberg atom system. 51,52 Recent works on periodically driven optical lattices have started to explore such physics. 53,54 On the other hand, kinetic constraints have been investigated as a possible pathway to many-body localization without disorder. 55 In classical systems, nonthermalizing behavior without disorder is well-known in the context of structural glasses. [56][57][58] The mechanism of this type of behavior is the excluded volume interactions that impose kinetic constraints on the dynamics. 59,60 Similar type of physics has recently been explored in quantum systems where a "quasi many-body localized" behavior was proposed to occur in the absence of disorder. [61][62][63][64][65][66][67][68][69][70][71] In this paper we address the question of the relation between kinetic constraint, slow dynamics and quantum many-body scars. In contrast to previous work, which primarily focused on spins and fermions that are closely related in one dimension, we study kineticallyconstrained models of bosons in one-dimension. This paper is organized as follows. In Sec. II we introduce the models of bosons with correlated hopping and discuss properties of their Hamiltonians when viewed as adjacency matrices of graphs in the Fock space. In Sec. III we investigate thermalization properties of these models by studying their energy level statistics, entanglement entropy of eigenstates, and dynamics under global quench. In Sec. IV 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 such periodic revivals. Finally, our conclusions and a discussion of experimental realizations is presented in Sec. V.

II. 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 .

A. Models
We study three different models, defined by the Hamiltonians: All three models contain a free-boson hopping term, b † j b j+1 , which is dressed in various ways by density operators, n j = b † 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 freeboson 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 To avoid clutter, we do not label the vertices in (d), (e) and (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 H1. 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 H2, almost all configurations are well-connected to the rest of the graph. The graphs for H3 show that the Hilbert space is highly reducible: its graph splits into many disconnected components.
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 Sec. II B. 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 labelled 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 Sec. II B below), 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 .

B. Graph structure and bipartite lattice
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 . 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 entanglement entropy. We will provide a more quantitative description and examples that illustrate this in Sec. IV. Finally, the graph of H 3 , due to even stronger constraints, is actually disconnected. This predicts that thermalization and dynamics in the three models will be very different, which we will confirm in Section III. 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, which is discussed in Appendix A.
We note that the relation between H 1 and H 3 is reminiscent of the relation between the quantum East model 72 and the "PXP" model describing the atoms in the Rydberg blockade regime. 6,19,42 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. 72 Hamiltonian H 1 has a bipartite structure, 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 We define n even and n odd as the total numbers of particles at even and odd sites respectively 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 (see Appendix B). 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 , which we define as In this case, the two sets are the configurations with even and with odd distances d a . 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. 73,74 We discuss in detail this zero-energy manifold in Appendix C.

III. DYNAMICS AND ENTANGLEMENT PROPERTIES
We now investigate the phenomenology of the models introduced in Eqs. (1)-(3). We use full 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.

A. Level statistics and entanglement entropy
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) 75 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 ) .  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. 76. For standard random matrix theory ensembles, both P (r) and the mean r are well-known. 77 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 Figs. 2(a), (b) and (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 that we discuss in Appendix C). We exclude the edges of the spectrum because they contain degeneracies which are not symmetry-related. 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. 2(c) 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 L B = L − L A . For a given pure state |ψ , the entanglement entropy is defined as where ρ A = tr B |ψ ψ| is the reduced density matrix of the subsystem A. cent of PXP model. 19,24 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. 78

B. 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 translationinvariant |ψ 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 time-evolved state, we evaluate the quantum fidelity, i.e., the probability for the wave function to return to the initial state. In a thermalizing system, fidelity decays rapidly to zero by the time t * ∼ O(1/J) and remains exponentially suppressed at later times. In scarred models, such as the Rydberg atom chain, 19 fidelity was found to revive to a finite fraction even at long times t 1/J. 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 |...0N 10... . In contrast, a completely uniform configuration |...111... 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 3 lattice sites: Time evolution of the fidelity for the initial state |(210) n for different system sizes L = 3n is shown in Fig. 3(a). The initial state is assumed to be the product state, e.g., |ψ 0 = |210 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., |ψ 0 = 1 √ 3 (|210 + |021 + |102 ). 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. 3(b) we compare the fidelity for the initial state in Eq. (12) 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. 3(c) 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 translationinvariant, it is important how we cut the system, e.g., |...210|210... versus |...2102|10... . 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 con-stant value, which is the case in Fig. 3(c). In Fig. 4 we show the 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 . 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 state 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.

IV. PERIODIC REVIVALS
In this Section, we introduce a toy model that explains the dynamics and predicts the frequency of revivals in H 1 from the initial state |(210) n in Eq. (12). 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.

A. H3 model: Perfect revivals
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 3 are E 1 = −2 and E 2 = 2. The initial state |ψ 1 (t = 0) = |210 evolves as and the state |ψ 2 (t = 0) = |120 evolves as 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 consists 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 Eq. (15), the initial state |ψ n (t = 0) = |(210) n evolves as |ψ L=3n (t) = cos n (2t)|(210) n + (−i) n sin n (2t)|(120) n + ...
where "..." denotes contributions of the basis configurations other than |(210) n or |(120) n . 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 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 is the configuration |(012) n , which is a translation of the state |(120) n that belongs to the connected subspace of |(210) n . The initial state evolves as |ψ inv n (t) = (cos n 2t + (−i) n sin n 2t) |ψ inv n (t = 0) + ... (21) and the fidelity is The frequency of the revivals is now doubled, so the period is T inv 3 = π/4. B. H1 model: Cluster approximation In Sec. II, we have demonstrated that the Hilbert space of H 1 consists of several weakly-connected clusters of configurations. We will now show that this cluster structure is responsible for the H 1 dynamics observed in Fig. 3. Some initial configurations are located in parts of the Hilbert space that are weakly connected to the rest and this leads to the wave-function "getting stuck" in some of those parts instead of spreading through the full Hilbert space. Here we identify such clusters that contain the configurations (210) n , which exhibit interesting dynamics with periodic revivals, as found in Fig. 3.

Cluster approximation
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. Fig. 5 shows part of the graph that contains the initial state, |210210 . Configurations labelled inside the ellipses denote representatives of an orbit of translation symmetry, i.e., the configurations are translation-invariant such as the one in Eq. (19). The minimal subcluster of the graph is highlighted in blue color in Fig. 5. 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 lower in magnitude than those inside the cluster, meaning that the probability is higher to stay inside the cluster than to leave. Inside this cluster, the evolution of the configuration |210210 can 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 . The cluster from Fig. 5 contains all the states given by tensor products of 210, 120 and 300 configurations.
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 E 1 = −4, E 2 = 4, E 3 = 0. The initial configuration |210 evolves according to From this, we can read off the period of revivals: T 1 = π/4. Similar to the H 3 case, we can now generalize to larger systems. Within the cluster approximation, triplets of sites again evolve independently, and the dimension of the reduced Hilbert space is The time-evolved state within the cluster approximation is given by |ψ c n (t) = cos n (4t)|(210) n + ...
As in the case of H 3 , this result is also valid for the translation-invariant initial state. 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. 6. The frequency of the fidelity revival, shown by the blue line in Fig. 6(a), 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. 6(b), 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 takes values that are several times larger.

Improved cluster approximation
To obtain a more accurate approximation, we now expand the cluster with several neighboring configurations: |111300 , |111210 , |111120 and |111111 . This enlarged cluster contains the minimal cluster studied previously, but it also includes additional configurations shown in green ellipses in Fig. 5. Similar clusters can be defined in larger systems if L = N p is divisible by 3, by using tensor products of any of the configurations 210, 120, 300 and 111.
The dimension of the extended cluster grows with the system size as and is thus exponentially larger than the previous cluster approximation. Nevertheless, its dimension is still exponentially smaller compared to the full Hilbert space, and within this approximation it is possible to simulate the dynamics of larger systems, L 30 -see Fig. 7. The revivals are no longer perfect and their frequency slightly changes with the system size, but it converges for large L. The frequency is closer to the frequency of revivals for the full Hilbert space than in the original cluster approximation, as previously shown in Fig. 6. For the initial product state (210) n , it is possible to analytically obtain the fidelity within the improved approximation for arbitrary system size. Similar to the previous methods, it can be shown where α = 9 + √ 57 ≈ 4.06815, β = 9 − √ 57 ≈ 1.20423, b ≈ 0.694113 and d ≈ 0.134933. These equations were obtained by solving the eigenproblem for L = 3 and generalizing the results to larger systems. The derivation can be found in Appendix D.
Eq. (30) is in excellent agreement with the numerical results from Fig. 7(a). We also found it to be a very good approximation for the translation-invariant initial state when L ≥ 9 (data not shown). 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. 6 and 8. Fig. 8(a) shows that the negative logarithm of the first peak height per site − log(F (T ))/L saturates at a finite value for large L. In Appendix D we show that the first peak height in the improved cluster approximation decreases with the system size as e −0.04L . 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 results in the fidelity of a random state F ∼ e −1.39L , which decays significantly faster than the first peak height e −0.04L . The evolution of the entanglement entropy for the extended cluster approximation is shown in Fig. 6(b). 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. 9 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 approximately predicts the revival peak heights for the full dynamics in Fig. 6. We conclude this Section by discussing the relation between H 3 and H 1 reduced to the cluster approximation. For the dynamics from the initial state |(210) n , the two models yield similar results, compare Eqs. (18) and (30). 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 proved analytically in Sec. IV, see Eq. (22). 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 that the symmetry of the initial state matches the symmetry of the Hamiltonian.

C. Generalization to other clusters
The eigenstates of H 1 , projected to the subspace of the minimal cluster approximation (as defined in Sec. IV B), 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 . 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. 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. 10. The lowest entropy eigenstates have high overlaps with other configurations, such as |(3100) 3 , as shown in Fig. 10(b). In the case of (210) n , 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.
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 where N is the length of the unit cell (L = N n).
For example, some of these configurations are |(3100) n , |(41000) n and |(510000) n . Combinations of those patterns such as |310041000 also 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. (33). As in the case of |(210) n , 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 is solely connected to the configurations |i ± 1 , except the two configurations at the edges, |1 and |N , which are only connected to |2 and |N − 1 , 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 . This result can be directly generalized to configurations of the form |((N − 1)10...0) n . The derivation is valid for both translation-invariant and non-translation-invariant initial configurations, as the cluster in Eq. (34) is disconnected from its translated copies. We stress that this disconnection, namely the absence of a hopping term between |1(N − 1)0...0 and |0N 0...0 , 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 (for n = 1, 2, 3, 4), |(41000) n (n = 1, 2, 3) and |(510000) n (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 case, the approximation can be improved by adding further configurations to the clusters. Moreover, if we want to consider translationinvariant initial states, we can follow the same recipe for |(210) n by 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.

V. CONCLUSIONS
In this paper, we have introduced three models of bosons with "soft" kinetic constraints, i.e., densitydependent 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 , 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. 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. Finally, we have also shown that the introduced models contain additional special properties, like the exponentially large zero-energy degeneracy which is related to the bipartite structure of the model.
Finally, we comment on the possible experimental realizations of the models we studied. The implementation of a correlated hopping term (n k b † i b j ) in optical lattices has attracted lot of attention due to a possible onset of quantum phases related to high-Tc superconductivity. 79 An early theoretical proposal exploits asymmetric interactions between the two atomic states in the presence of a state-dependent optical lattice. 79 As a result, the obtained effective model corresponds to the inversionsymmetric 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. [80][81][82][83] 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. 84,85 More recently, the set of quantum models accessible in cold-atom experiments has been enriched through the technique of Floquet engineering. 86 As a notable example, a suitable driving scheme can renormalize or fully suppress the bare tunneling rate. 87 On top of that, by modulating local interactions an effective model with the density-dependent tunneling term has been engineered. 88 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. 89 Related sophisticated driving schemes have already enabled a realization of dynamical gauge fields 53,54,90 where both the amplitude and the phase of the effective tunneling are density-dependent.
Note added: During the completion of this work, we became aware of Ref. 91 which identified non-thermal eigenstates and slow dynamics in the quantum East model. These two Hamiltonians have the same constraints as the Hamiltonian H 1 and therefore the same graphs as in Fig. 1. However, the particle number operators n j are squared in H 1a and replaced with delta functions (1 − δ nj ,0 ) in H 1b . This makes the clusters from Section IV B even less connected to the rest of the configurations in the case of H 1a and more connected in the case of H 1b .

VI. ACKNOWLEDGMENTS
As anticipated, the revivals become more prominent for H 1a , with fidelity peaks reaching more than 95%, while the peaks almost disappear for H 1b , as illustrated in Fig. 11(a). In addition, the entanglement entropy quickly saturates in the case H 1b , while the growth is significantly suppressed in the case of H 1a , as can be observed in Fig. 11(b). The distribution of entanglement entropy across all eigenstates is also affected by the change of coefficients (not shown). H 1a has a spectrum with many low-entropy eigenstates, while the spectrum of H 1b is almost thermal and resembles that of H 2 . The probability distribution of consecutive gaps in the energy spectrum of H 1a is close to the Poisson distribution, which implies that H 1a is almost integrable. On the other hand, the distribution for H 1a is Wigner-Dyson, like in H 1 .

Appendix B: Bipartite structure of H1
We have shown in Sec. II B that the H 1 model is bipartite for open boundary conditions irrespective of the system size L parity or for periodic boundary conditions when L is even. We prove in this appendix that this property of H 1 holds true when L is odd and ν = 1.
Due to the constraints imposed by H 1 , a particle cannot hop to an empty site to its left. At the filling factor ν = 1, all configurations except 111...111 contain at least one empty site. These configurations can be connected to 111...111 by hoppings only to the right, which is also the shortest possible path d a , defined in Eq. (7). The empty sites can be filled only with particles that come from the site on their left, as hopping from the other side is forbidden. This implies that at least for one pair of adjacent sites (an empty site and the filled one on its right) there will be no particles hopping between them on the path to 111...111. We can then redefine the numbering of sites to start from the filled site in this pair. This is equivalent to setting the right (filled) site as the first and the left (empty) site as the last site in the chain and imposing open boundary conditions. In this way, no two odd sites will be coupled and the argument that the absolute difference between the numbers of particles on even and odd sites can only change by ±1 will still be valid.
Unlike H 1 , H 2 in the same geometry is not bipartite. The reason for this is that there are no constraints in the case of H 2 , so the shortest path to 111...111 can include hoppings both to the right and to the left, which means that it is not always possible to choose the numbering in such a way that no two odd sites are coupled. Because of the open boundary conditions, the Hamiltonian H 3 in its largest connected component is also bipartite for all system sizes.

Appendix C: Zero modes
An interesting feature of H 1 model is the large number of zero energy states in the middle of its spectrum. The number of these zero modes, found by brute force diagonalization, is listed in Table I for different system sizes and momentum sectors. Similar property is found for H 2 -see Table II, with the notable difference that there are no zero modes when the number of sites L is odd.    The origin of the zero modes is the underlying bipartite structure of the Hamiltonian. 73 H 1 does not connect configurations within the same class. For example, a graph that shows how the configurations for L = N p = 4 are connected is displayed in Fig. 12.
Here we will refer to the two classes as the "green" (even) and the "red" (odd) configurations. Each basis configuration can be uniquely assigned a green or a red label according to the parity of its distance d a , Eq. (7), from the configuration |111...111 . If this number is even, the configuration is green, and if it is odd, the configuration is red. This separation into two classes is a consequence of the constraints present in the Hamiltonian H 1 . Hamiltonians without such constraints, for example H 2 or the standard Bose-Hubbard model, do not exhibit this bipartite structure for odd system sizes, see Section II B and Appendix B more details. In these cases, it is not possible to uniquely determine whether a particular configuration is green or red. The lack of bipartite structure is the reason for the absence of zero energy eigenstates of H 2 in odd dimensions, which was observed in Table II. However, the configuration space of H 2 is still bipartite in even dimensions, allowing for the existence of some zero modes in those cases.
Low-entropy zero-energy states can be constructed as superpositions of either only green or only red configurations. For example, in the case of L = N p = 4, the simplest and therefore the lowest-entropy zero mode can be constructed using only two green product states (encircled by a dashed line in Fig. 12) where | . . . T was defined in Eq. (19). There is another zero mode in this case, and it can be formed by adding more green configurations to the superposition. The number of zero-energy eigenstates is related to the difference between the numbers of green and red configurations, 73 as we will now explain.
As the Hamiltonian H 1 only connects green configurations to red configurations and red to green, we can rewrite it in the following way: where |R i are the red product states and |G j are green. Its square, H 2 1 , connects green configurations to green and red to red, and it is therefore block diagonal. The blocks are CC † and C † C, where C is a matrix with the elements c ij . The dimensions of C are r × g, where r is the number of red configurations and g of green. C and C † can be factorized using singular value decomposition. From this structure we can see that the energy spectrum is symmetric around zero and that the minimal number of zero-energy states is |g −r|. The zero-energy eigenvectors can also be obtained as the ground states of H 2 1 . Similar analysis and counting of the zero modes in PXP model was performed in Ref. 6 and 92. Table III shows the difference between the numbers of red and green states g − r for different system sizes and the number of zero-energy states N 0 in those systems. The number of zero-energy states, found by exact diagonalization, in all cases satisfies the anticipated inequality N 0 ≥ |g − r|. In fact, the bound is almost always saturated, N 0 = |g −r|, except when L = N p = 4n+2, n ∈ Z. Interestingly, the minimal number of zero modes |g − r| for L = N p = 4n is equal to the Hilbert space dimension for L = N p = 2n (both total and for k = 0 sector only). This leads to the conclusion that the number of zero-energy states grows exponentially with the system size. It can also be noticed that the total difference g − r for L = N p = 2n + 1 is always twice the difference for L = N p = 2n.  Table III. The difference between the number of green and red configurations g − r and the number of zero-energy states N0 (determined by exact diagonalization) for different system sizes. Overall, the derived bound for the number of zero modes is found to be very tight in finite systems where it can be independently confirmed by explicit diagonalization ("NA" denotes cases where this was not possible).
Appendix D: Derivation of fidelity in the extended cluster approximation for |(210) n For system size L = 3, the Hilbert space of the extended cluster is spanned by only four configurations: The Hamiltonian reduced to this subspace is Its eigenvalues are and its eigenvectors where α = 9 + √ 57 ≈ 4.06815, β = 9 − √ 57 ≈ 1.20423, a ≈ 0.591050, b ≈ 0.694113, c ≈ 0.388150 and d ≈ 0.134933. There are no simple analytical expressions for the coefficients a, b, c and d.
The period of revivals is approximately T ≈ π/α ≈ 0.772241, and the first peak height exponentially decreases as Eqs. (D6) and (30) are exact for non-translation-invariant initial states, but just an approximation for the translation symmetric case. This is due to the fact that different translations of 300, 210 and 120 no longer evolve independently in that case, as they are connected to each other through the configuration 111. However, this approximation becomes better with increasing the system size, as the configuration (111) n becomes further away from the initial state (210) n and the probability that this configuration will be reached decreases.