<figure><img src="image_1.png" /></figure>

HAL
open science

# Quantum scars of bosons with correlated hopping

Ana Hudomal, Ivana Vasić, Nicolas Regnault, Zlatko Papic

► To cite this version:

Ana Hudomal, Ivana Vasić, Nicolas Regnault, Zlatko Papic. Quantum scars of bosons with correlated hopping. Communications Physics, 2020, 3 (1), pp.99. 10.1038/s42005-020-0364-9. hal-03020106

HAL Id: hal-03020106

https://hal.science/hal-03020106v1

Submitted on 28 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

COMMUNICATIONS
PHYSICS

ARTICLE

https://doi.org/10.1038/s42005-020-0364-9

OPEN

<figure><img src="image_4.png" /></figure>

Check for updates

# Quantum scars of bosons with correlated hopping

Ana Hudomal<fnref n="1" />, Ivana Vasić<fnref n="1" />, Nicolas Regnault²,³ & Zlatko Papić⁴

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.

<footnote n="1">¹Institute of Physics Belgrade, University of Belgrade, 11080 Belgrade, Serbia. ²Joseph Henry Laboratories and Department of Physics, Princeton University, Princeton, NJ 08544, USA. ³Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot, Sorbonne Paris Cité, 75005 Paris, France. ⁴School of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, UK.</footnote>

<footnote>✉email: ana.hudomal@ipb.ac.rs</footnote>
Semiclassical studies of chaotic stadium billiards have revealed the existence of remarkable non-chaotic eigenfunctions called “quantum scars”¹. Scarred eigenfunctions display anomalous enhancement in regions of the billiard that are traversed by one of the periodic orbits in the classical limit when $h \to 0$. It was shown that quantum scars lead to striking experimental signatures in a variety of systems, including microwave cavities², quantum dots³, and semiconductor quantum wells⁴.

A recent experiment on a quantum simulator⁵, and subsequent theoretical work⁶,⁷, 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⁵,⁸,⁹, where nearest-neighbor excitations of the atoms were energetically prohibited. The experiment observed persistent many-body revivals of local observables after a “global quench”¹⁰ 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 behavior compared to previous studies of quantum thermalization in various experimental platforms¹¹⁻¹⁵.

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¹, 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¹⁶,¹⁷. 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 many-body case, non-ergodic eigenstates are broadly defined as those that violate eigenstate thermalization hypothesis (ETH)¹⁸,¹⁹. Scarred eigenstates violate the ETH in a number of ways: for example, they appear at evenly spaced energies throughout the spectrum⁶,²⁰,²¹, 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²⁰.

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²²,²³. A few of such exact eigenstates are now also available for the Rydberg atom chain model²⁴. The collection of models that feature atypical eigenstates is rapidly expanding, including perturbations of the Rydberg atom chain²⁰,²⁵,²⁶, theories with confinement²⁷⁻²⁹, Fermi-Hubbard model beyond one dimension³⁰,³¹, driven systems³², quantum spin systems³³,³⁴, fractional quantum Hall effect in a one-dimensional limit³⁵, and models with fracton-like dynamics³⁶⁻³⁹. In a related development, it was proposed that atypical eigenstates of one Hamiltonian can be “embedded” into the spectrum of another, thermalizing Hamiltonian⁴⁰, causing a violation of a “strong” version of the ETH⁴¹,⁴². This approach allows to engineer scarred eigenstates in models of topological phases in arbitrary dimensions⁴³. 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⁴⁴.

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⁴⁵. Such Hilbert spaces occur, for example, in models describing anyon excitations in topological phases of matter⁴⁶⁻⁵⁰ and in lattice gauge theories⁵¹⁻⁵³, including the Rydberg atom system⁵⁴,⁵⁵. Recent works on periodically driven optical lattices have started to explore such physics⁵⁶,⁵⁷. On the other hand, kinetic constraints have been investigated as a possible pathway to many-body localization without disorder⁵⁸. In classical systems, non-thermalizing behavior without disorder is well known in the context of structural glasses⁵⁹⁻⁶¹. The mechanism of this type of behavior is the excluded volume interactions that impose kinetic constraints on the dynamics⁶²,⁶³. 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⁶⁴⁻⁷⁴.

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³⁶⁻³⁹. 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, $v = N_p/L$, is conserved, and we will mainly present results in the dense regime, $v = 1$. We have studied models with $v < 1$ and $v > 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:

$$H_1 = -J \sum_{j=1}^{L} (b_j^\dagger b_{j+1} n_j + n_{j-1} b_j^\dagger b_{j-1}), \quad (1)$$

$$H_2 = -J \sum_{j=1}^{L} (n_j b_j^\dagger b_{j+1} + b_j^\dagger b_{j-1} n_{j-1}), \quad (2)$$

$$H_3 = -J \sum_{j=1}^{L} (n_{j+1} b_j^\dagger b_{j+1} n_j + n_{j-1} b_j^\dagger b_{j-1} n_j). \quad (3)$$

All three models contain a free-boson hopping term, $b_j^\dagger b_{j+1}$, which is dressed in various ways by density operators, $n_j = b_j^\dagger 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,

$$H_2 = H_1 - J \sum_j (b_j^\dagger b_{j+1} + b_j^\dagger b_{j-1}), \quad (4)$$

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 particle-number 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 \equiv 1$) and set $\hbar = 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\pi/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 $v = 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 $v = 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 adjacency matrix of a graph whose vertices are Fock states of bosons, $|n_1, n_2, \dots, n_L\rangle$. 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

<figure><img src="image_5.png" /><figcaption>**Fig. 1** Connectivity of the Hilbert space. Adjacency graph for a $H_1$, b $H_2$, c $H_3$, all for $L = N_p = 3$. 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.</figcaption></figure>
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 = \frac{\min(s_n, s_{n+1})}{\max(s_n, s_{n+1})} \qquad (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 $\langle r \rangle$ 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 \to 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 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. 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

<figure><img src="image_6.png" /><figcaption>**Fig. 2** Level statistics and entanglement. **a–c** Probability distribution of the ratios of two consecutive energy gaps. **a** $H_1$ (middle third of the spectrum without $E = 0$ states, $L = N_p = 12$), **b** $H_2$ (full spectrum, $L = N_p = 12$) and **c** $H_3$ (largest connected component of $L = N_p = 22$). 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, $E/N_p$. **d** $H_1$ ($L = N_p = 12$, $L_A = 6$, $k = 0$), **e** $H_2$ (same) and **f** $H_3$ in the largest connected component of $L = N_p = 20$, $L_A = 10$, $I = 0$. The inset shows all connected components for $L = N_p = 12$, $L_A = 6$, $k = 0$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 3 Dynamics of quantum fidelity and entanglement entropy for initial configurations in Eq. (9). a Time evolution of fidelity $F(t)$ in Eq. (8) for system sizes $L = 3n$. The evolution is governed by the Hamiltonian $H_1$, different colors represent different system sizes $L$. b Fidelity evolution $F(t)$ for the Hamiltonians $H_1$, $H_2$ and $H_3$ and system size $L = 15$. c Entanglement entropy evolution $S_{L,A}(t)$ for the same cases as in b.</figcaption></figure>

$L_B = L - L_A$. For a given pure state $|\psi\rangle$, the entanglement entropy is defined as

$$S_A = -\mathrm{tr}_A(\rho_A \ln \rho_A), \qquad (6)$$

where $\rho_A = \mathrm{tr}_B|\psi\rangle\langle\psi|$ is the reduced density matrix of the subsystem A. The scatter plots, showing entanglement entropy of all eigenstates $|E_n\rangle$ 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, $|\psi_0\rangle$, at time $t=0$. We restrict to $|\psi_0\rangle$ 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 $|\psi_0\rangle$. After preparing the system in the state $|\psi_0\rangle$, which is not an eigenstate of the Hamiltonian, the system is let to evolve under unitary dynamics,

$$|\psi(t)\rangle = \exp\left(-\frac{i}{\hbar}Ht\right)|\psi_0\rangle. \qquad (7)$$

where $H$ is one of the Hamiltonians of interest. From the time-evolved state, we evaluate the quantum fidelity,

$$F(t) = |\langle\psi_0|\psi(t)\rangle|^2, \qquad (8)$$

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) \sim \exp(-L(Jt)^2)$. It thus becomes exponentially suppressed in the system size for any fixed time $t^*$, i.e., $F(t^*) \sim \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) \sim \exp(-c'L)$, with $c' \ll c$. In ref. $^{20}$, for a finite system with $L \lesssim 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...\rangle$. In contrast, a completely uniform configuration $|...111...\rangle$ 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:

$$|210210...\ 210\rangle \equiv |(210)^n\rangle. \qquad (9)$$

Time evolution of the fidelity for the initial state $|(210)^n\rangle$ for different system sizes $L=3n$ is shown in Fig. 3a. The initial state is assumed to be the product state, e.g., $|\psi_0\rangle = |210\rangle$ 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., $|\psi_0\rangle = \frac{1}{\sqrt{3}}(|210\rangle + |021\rangle + |102\rangle)$. 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\rangle$. 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., $|...210|210...\rangle$ versus $|...2102|10...\rangle$. 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.
<figure><img src="image_8.png" /><figcaption>Fig. 4 Evolution of local observables for the Hamiltonian $H_1$. a</figcaption></figure>

Correlations between adjacent sites $\langle n_1n_2(t) \rangle$ for different system sizes and the initial state $|(210)^n\rangle$. b Density on one site $\langle n_1(t) \rangle$.

In Fig. 4 we show the $H_1$ evolution of two local observables, density correlations between two adjacent sites $\langle n_1n_2(t) \rangle$ and density on the first site $\langle n_1(t) \rangle$, starting from the initial state $|(210)^n\rangle$. 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\rangle$ state 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\rangle$. 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\rangle$ can be solved exactly, as shown in Methods. The connected component of the state $|(210)^n\rangle$ 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 this observation, it can be shown that revivals will be perfect with a period $T_3 = \pi/2$. The same period is obtained for initial product state $|(210)^n\rangle$ 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.

<figure><img src="image_9.png" /><figcaption>Fig. 5 Comparison of the full dynamics against the minimal cluster (1) and extended cluster (2) approximation schemes. We consider the system size $L=15$ with the initial state $|(210)^5\rangle$. a Time evolution of the fidelity. The frequency of revivals is approximately the same in both cases, but the results for the extended cluster show better agreement with the results for the full Hilbert space. b Time evolution of the entanglement entropy.</figcaption></figure>

For our initial state $|(210)^n\rangle$, 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 \mathcal{H}^c = 3^{L/3}$. The time-evolved state within the cluster is given by

$$|\psi_n^c(t)\rangle = \cos^n(4t)|(210)^n\rangle + \dots, \quad (10)$$

where the dots denote other configurations. The fidelity is

$$F_n^c(t) = |\langle\psi_n^c(0)|\psi_n^c(t)\rangle|^2 = |\cos(4t)|^{2n}. \quad (11)$$

As in the case of $H_3$, this result is also valid for the translation-invariant initial state. We see that the period of revivals is $T_1 = \pi/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 \sim 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.
<figure><img src="image_10.png" /><figcaption>**Fig. 6** Cluster approximations. **a** Fidelity $F(t)$, for the Hamiltonian $H_1$ and initial states $|(210)^n\rangle$, in the extended cluster approximation for various system sizes. **b** Eigenstate overlap with the initial state $|(210)^8\rangle$ plotted on a log scale, for both cluster approximations. In the case of degenerate eigenstates the sum of their overlaps is shown.</figcaption></figure>

The dimension of the extended cluster grows as $dim \mathcal{H}^c = 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 \lesssim 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 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\rangle$ 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. Similar to the previous methods, it can be shown (see Supplementary Note 2)

$$F_{L=3n}^{\bar{c}}(t) = 4^n |b^2 \cos(\alpha t) + d^2 \cos(\beta t)|^{2n}, \quad (12)$$

where $\alpha = \sqrt{9 + \sqrt{57}} \approx 4.06815$, $\beta = \sqrt{9 - \sqrt{57}} \approx 1.20423$, $b \approx 0.694113$ and $d \approx 0.134933$. Eq. (12) is in excellent agreement with the numerical results in Fig. 6a. It was also found to be a very good approximation for the translation-invariant initial state when $L \ge 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 \sim 1/\dim \mathcal{H}$. In the case $v = 1$ and large $L$, the Hilbert space dimension grows with the system

<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8** Evolution of the probability to remain inside the minimal cluster. $O_{\text{Cluster}}$ as defined in Eq. (14). Initial configurations are indicated in the legend. Solid lines: configurations initially inside the cluster. Dashed lines: configurations initially outside the cluster (all except $|(111)^4\rangle$ are randomly chosen). Similar results are obtained for the extended cluster (not shown). System size $L = 12$.</figcaption></figure>

size as

$$\dim_{\mathcal{H}} = \begin{pmatrix} 2L-1 \\ L \end{pmatrix} \approx \begin{pmatrix} 2L \\ L \end{pmatrix} \approx \frac{4^L}{\sqrt{\pi L}}. \quad (13)$$

This back-of-the-envelope estimate suggests the fidelity of a random state is $F \sim 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$,

$$O_{\text{Cluster}} = \sum_{a \in \text{Cluster}} |\langle a | \psi(t) \rangle|^2. \quad (14)$$

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\rangle$ 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\rangle$
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Non-ergodic eigenstates. **a** Overlap of the configuration $|(210)^4\rangle$ with all the eigenstates of $H_1$, $H_1^\text{c}$ and $H_1^\text{e}$ versus the eigenstate energy for sector $k=0$ and system size $L=N_p=12$. **b** Same for $|(3100)^3\rangle$. **c** Entanglement entropy, eigenstates which have the highest overlap with some product states are marked in different colors.</figcaption></figure>

approximately predicts the revival peak heights for the full dynamics.

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\rangle$, 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_3^{\text{inv}} = \pi/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 $\pi/4 \approx 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\rangle$, 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\rangle$, as shown in Fig. 9b, c. In the case of $|(210)^n\rangle$, 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\rangle$, 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\rangle$, 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 zero-energy 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_i^\dagger 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-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
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 density-assisted hopping, finding signatures of quantum many-body 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

$$ \Delta_a = \frac{|n_{\text{even}} - n_{\text{odd}} + C|}{2}, \qquad (15) $$

where $C=0$ if $L$ is even and $C=1$ if $L$ is odd. We define $n_{\text{even}}$ and $n_{\text{odd}}$ as the total numbers of particles at even and odd sites, respectively,

$$ n_{\text{even}} = \sum_{l=1}^{L_1} n_{2l}, \quad n_{\text{odd}} = \sum_{l=1}^{L_2} n_{2l-1}, \qquad (16) $$

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_{\text{even}}$ by one and decreases $n_{\text{odd}}$ by one, or vice-versa. This means that each hopping can change $\Delta_a$ only by $\pm 1$.

In special cases, like $H_1$ at filling factor $v=1$, it is also possible to define quantities like $\Delta_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\rangle$, which we define as

$$ d_a = \min_n \{ \langle 111...111 | H_1^n | a \rangle \neq 0 \}. \qquad (17) $$

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 $\pm 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 $\Delta_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 $\Delta_a$ can only change by $\pm 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)^{\delta_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 graph – such 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

$$ H'_3 = - \begin{pmatrix} 0 & 2 \\ 2 & 0 \end{pmatrix}, \qquad (18) $$

where the basis vectors are

$$ \begin{pmatrix} 1 \\ 0 \end{pmatrix} = |210\rangle, \quad \begin{pmatrix} 0 \\ 1 \end{pmatrix} = |120\rangle. \qquad (19) $$

The eigenvalues of $H'_3$ are $E_1 = -2$ and $E_2 = 2$. The initial state $|\psi_1(t=0)\rangle = |210\rangle$ evolves as

$$ |\psi_1(t)\rangle = \cos(2t)|210\rangle - i \sin(2t)|120\rangle, \qquad (20) $$

and the state $|\psi_2(t=0)\rangle = |120\rangle$ evolves as

$$ |\psi_2(t)\rangle = -i \sin(2t)|210\rangle + \cos(2t)|120\rangle. \qquad (21) $$

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\rangle$ consists only of combinations of patterns 210 and 120, which means that triplets of sites evolve independently. From Eq. (20), the initial state $|\psi_n(t=0)\rangle = |(210)^n\rangle$ evolves as

$$ |\psi_{L=3n}(t)\rangle = \cos^n(2t)|(210)^n\rangle + (-i)^n \sin^n(2t)|(120)^n\rangle + \dots \qquad (22) $$

where “...” denotes contributions of the basis configurations other than $|(210)^n\rangle$ or $|(120)^n\rangle$. The fidelity is

$$ F_{L=3n}(t) = |\langle\psi_n(0)|\psi_n(t)\rangle|^2 = |\cos 2t|^{2n}. \qquad (23) $$

It follows that the revivals are perfect, with a period $T_3 = \pi/2$. This result is also valid for the translation-invariant initial state $|(210)^n\rangle_T$,

$$ |(210)^n\rangle_T = \frac{1}{\sqrt{3}} \left( |(210)^n\rangle + |(021)^n\rangle + |(102)^n\rangle \right), \qquad (24) $$

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\rangle$ is the configuration $|(012)^n\rangle$, which is a translation of the state $|(120)^n\rangle$ that belongs to the connected subspace of $|(210)^n\rangle$. The initial state

$$ |\psi_n^{\text{inv}}(t=0)\rangle = \frac{1}{\sqrt{2}}|(210)^n\rangle_T + \frac{1}{\sqrt{2}}|(120)^n\rangle_T \qquad (25) $$

evolves as

$$ |\psi_n^{\text{inv}}(t)\rangle = (\cos^n 2t + (-i)^n \sin^n 2t) |\psi_n^{\text{inv}}(t=0)\rangle + \dots \qquad (26) $$

and the fidelity is

$$ \begin{aligned} F_n^{\text{inv}}(t) &= |\langle\psi_n^{\text{inv}}(0)|\psi_n^{\text{inv}}(t)\rangle|^2 \\ &= |\cos^n 2t + (-i)^n \sin^n 2t|^2. \end{aligned} \qquad (27) $$

The frequency of the revivals is now doubled, so the period is $T_3^{\text{inv}} = \pi/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\rangle$. 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 particle-number 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\rangle$ 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$. 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\rangle$ and evolving with $H_1$
<figure><img src="image_14.png" /><figcaption>**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 below the graph show the distance $d_o$ from the configuration 111111 evaluated using Eq. (17).</figcaption></figure>

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

$$ \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} = |300\rangle, \quad \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix} = |210\rangle, \quad \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} = |120\rangle. \qquad (28) $$

The Hamiltonian reduced to this subspace is

$$ H_1^c = - \begin{pmatrix} 0 & 2\sqrt{3} & 0 \\ 2\sqrt{3} & 0 & 2 \\ 0 & 2 & 0 \end{pmatrix}, \qquad (29) $$

and its eigenvalues are $E_1 = -4$, $E_2 = 4$, $E_3 = 0$. The initial configuration $|210\rangle$ evolves according to

$$ |\psi_1^c(t)\rangle = -\frac{i}{2}\sin(4t)(\sqrt{3}|300\rangle + |120\rangle) + \cos(4t)|210\rangle. \qquad (30) $$

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)^n\rangle$ 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

$$ \left| \left( (N-1) 1 \underbrace{0...0}_{N-2} \right)^n \right\rangle, \qquad (31) $$

where $N$ is the length of the unit cell ($L = Nn$). For example, some of these configurations are $|(3100)^n\rangle$, $|(41000)^n\rangle$ and $|(510000)^n\rangle$. Combinations of those patterns such as $|310041000\rangle$ 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. (31). As in the case of $|(210)^n\rangle$, 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

$$ |i\rangle = \left| \left( (N+1-i)(i-1) \underbrace{0...0}_{N-2} \right)^n \right\rangle, \qquad (32) $$

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\rangle$ is solely connected to the configurations $|i \pm 1\rangle$, except the two configurations at the edges, $|1\rangle$ and $|N\rangle$, which are only connected to $|2\rangle$ and $|N-1\rangle$, respectively.

The projection of the Hamiltonian $H_1$ to this cluster, which we denote by $H_1^c$, 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:

$$ H_1^c_{i,i+1} = H_1^c_{i+1,i} = (N-i)\sqrt{i(N+1-i)}. \qquad (33) $$

The dynamics within a single unit cell under $H_1^c$ corresponds to density fluctuations between the first and the second site. Following the same procedure as previously, we can now diagonalize $H_1^c$ and compute the fidelity time series for the initial configuration $|(N-1)10...0\rangle$. This result can be directly generalized to configurations of the form $|(N-1)10...0\rangle^n$. 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\rangle$ and $|0N0...0\rangle$, 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\rangle$ (for $n = 1, 2, 3, 4$), $|(41000)^n\rangle$ ($n = 1, 2, 3$) and $|(510000)^n\rangle$ ($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\rangle$ case, 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\rangle$ by summing translated patterns with the required phase factors given in terms of momenta in multiples of $2\pi/N$. We have checked that revivals appear in these momentum sectors, with roughly the same frequency.

We note that the configurations with larger unit 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 \sim 2$, and $(3100)^3$ overtakes $(510000)^2$ around $t \sim 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/810100.
## Code availability

Code is available upon reasonable request.

Received: 6 November 2019; Accepted: 5 May 2020;

Published online: 01 June 2020

## References

1. Heller, E. J. Bound-state eigenfunctions of classically chaotic Hamiltonian systems: Scars of periodic orbits. *Phys. Rev. Lett.* **53**, 1515-1518 (1984).

2. Sridhar, S. Experimental observation of scarred eigenfunctions of chaotic microwave cavities. *Phys. Rev. Lett.* **67**, 785-788 (1991).

3. Marcus, C. M., Rimberg, A. J., Westervelt, R. M., Hopkins, P. F. & Gossard, A. C. Conductance fluctuations and chaotic scattering in ballistic microstructures. *Phys. Rev. Lett.* **69**, 506-509 (1992).

4. Wilkinson, P. B. et al. Observation of 'scarred' wavefunctions in a quantum well with chaotic electron dynamics. *Nature* **380**, 608–610 (1996).

5. Bernien, H. et al. Probing many-body dynamics on a 51-atom quantum simulator. *Nature* **551**, 579–584 (2017).

6. Turner, C. J., Michailidis, A. A., Abanin, D. A., Serbyn, M. & Papić, Z. Weak ergodicity breaking from quantum many-body scars. *Nat. Phys.* **14**, 745–749 (2018).

7. Ho, W. W., Choi, S., Pichler, H. & Lukin, M. Periodic orbits, entanglement, and quantum many-body scars in constrained models: Matrix product state approach. *Phys. Rev. Lett.* **122**, 040603 (2019).

8. Schauß, P. et al. Observation of spatially ordered structures in a two- dimensional Rydberg gas. *Nature* **491**, 87–91 (2012).

9. Labuhn, H. et al. Tunable two-dimensional arrays of single Rydberg atoms for realizing quantum Ising models. *Nature* **534**, 667–670 (2016).

10. Calabrese, P. & Cardy, J. Time dependence of correlation functions following a quantum quench. *Phys. Rev. Lett.* **96**, 136801 (2006).

11. Kinoshita, T., Wenger, T. & Weiss, D. S. A quantum Newton's cradle. *Nature* **440**, 900–903 (2006).

12. Schreiber, M. et al. Observation of many-body localization of interacting fermions in a quasirandom optical lattice. *Science* **349**, 842–845 (2015).

13. Smith, J. et al. Many-body localization in a quantum simulator with programmable random disorder. *Nat. Phys.* **12**, 907–911 (2016).

14. Kucsko, G. et al. Critical thermalization of a disordered dipolar spin system in diamond. *Phys. Rev. Lett.* **121**, 023601 (2018).

15. Kaufman, A. M. et al. Quantum thermalization through entanglement in an isolated many-body system. *Science* **353**, 794–800 (2016).

16. Choi, S. et al. Emergent SU(2) dynamics and perfect quantum many-body scars. *Phys. Rev. Lett.* **122**, 220603 (2019).

17. Bull K., Desaules J.-Y. & Papić Z. Quantum scars as embeddings of weakly broken Lie algebra representations. *Phys. Rev. B* **101**, 165139 (2020).

18. Deutsch, J. M. Quantum statistical mechanics in a closed system. *Phys. Rev. A* **43**, 2046–2049 (1991).

19. Srednicki, M. Chaos and quantum thermalization. *Phys. Rev. E* **50**, 888 (1994).

20. Turner, C. J., Michailidis, A. A., Abanin, D. A., Serbyn, M. & Papić, Z. Quantum scarred eigenstates in a Rydberg atom chain: entanglement, breakdown of thermalization, and stability to perturbations. *Phys. Rev. B* **98**, 155134 (2018).

21. Iadecola, T., Schecter, M. & Xu, S. Quantum many-body scars from magnon condensation. *Phys. Rev. B* **100**, 184312 (2019).

22. Moudgalya, S., Rachel, S., Bernevig, B. A. & Regnault, N. Exact excited states of nonintegrable models. *Phys. Rev. B* **98**, 235155 (2018).

23. Moudgalya, M., Regnault, N. & Bernevig, B. A. Entanglement of exact excited states of Affleck-Kennedy-Lieb-Tasaki models: exact results, many-body scars, and violation of the strong eigenstate thermalization hypothesis. *Phys. Rev. B* **98**, 235156 (2018).

24. Lin, C.-J. & Motrunich, O. I. Exact quantum many-body scar states in the Rydberg-blockaded atom chain. *Phys. Rev. Lett.* **122**, 173401 (2019).

25. Khemani, V., Laumann, C. R. & Chandran, A. Signatures of integrability in the dynamics of Rydberg-blockaded chains. *Phys. Rev. B* **99**, 161101 (2019).

26. Michailidis, A. A., Turner, C. J., Papić, Z., Abanin, D. A. & Serbyn, M. Slow quantum thermalization and many-body revivals from mixed phase space. *Phys. Rev. X* **10**, 011055 (2020).

27. Kormos, M., Collura, M., Takács, G. & Calabrese, P. Real-time confinement following a quantum quench to a non-integrable model. *Nat. Phys.* **13**, 246–249 (2016).

28. James, A. J. A., Konik, R. M. & Robinson, N. J. Nonthermal states arising from confinement in one and two dimensions. *Phys. Rev. Lett.* **122**, 130603 (2019).

29. Robinson, N. J., James, A. J. A. & Konik, R. M. Signatures of rare states and thermalization in a theory with confinement. *Phys. Rev. B* **99**, 195108 (2019).

30. Vafek, O., Regnault, N. & Bernevig, B. A. Entanglement of exact excited eigenstates of the Hubbard model in arbitrary dimension. *SciPost Phys.* **3**, 043 (2017).

31. Iadecola, T. & Žnidarič, M. Exact localized and ballistic eigenstates in disordered chaotic spin ladders and the Fermi-Hubbard model. *Phys. Rev. Lett.* **123**, 036403 (2019).

32. Haldar, A., Sen, D., Moessner, R. & Das, A. Scars in strongly driven Floquet matter: resonance vs emergent conservation laws. Preprint at http://arxiv.org/abs/1909.04064 (2019).

33. Schecter, M. & Iadecola, T. Weak ergodicity breaking and quantum many-body scars in spin-1 XY magnets. *Phys. Rev. Lett.* **123**, 147201 (2019).

34. Iadecola, T. & Schecter, M. Quantum many-body scar states with emergent kinetic constraints and finite-entanglement revivals. *Phys. Rev. B* **101**, 024306 (2020).

35. Moudgalya, S., Bernevig, B. A. & Regnault, N. Quantum many-body scars in a Landau level on a thin torus. Preprint at http://arxiv.org/abs/1906.05292 (2019).

36. Pai, S. & Pretko, M. Dynamical scar states in driven fracton systems. *Phys. Rev. Lett.* **123**, 136401 (2019).

37. Khemani, V. & Nandkishore, R. Local constraints can globally shatter Hilbert space: a new route to quantum information protection. Preprint at http://arxiv.org/abs/1904.04815 (2019).

38. Sala, P., Rakovszky, T., Verresen, R., Knap, M. & Pollmann, F. Ergodicity breaking arising from Hilbert space fragmentation in dipole-conserving Hamiltonians. *Phys. Rev. X* **10**, 011047 (2020).

39. Khemani, V., Hermele, M. & Nandkishore, R. M. Localization from shattering: higher dimensions and physical realizations. Preprint at http://arxiv.org/abs/1910.01137 (2019).

40. Shiraishi, N. & Mori, T. Systematic construction of counterexamples to the eigenstate thermalization hypothesis. *Phys. Rev. Lett.* **119**, 030601 (2017).

41. D'Alessio, L., Kafri, Y., Polkovnikov, A. & Rigol, M. From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics. *Adv. Phys.* **65**, 239–362 (2016).

42. Gogolin, C. & Eisert, J. Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems. *Rep. Prog. Phys.* **79**, 056001 (2016).

43. Ok, S. et al. Topological many-body scar states in dimensions one, two, and three. *Phys. Rev. Res.* **1**, 033144 (2019).

44. Bull, K., Martin, I. & Papić, Z. Systematic construction of scarred many-body dynamics in 1D lattice models. *Phys. Rev. Lett.* **123**, 030601 (2019).

45. Lesanovsky, I. & Katsura, H. Interacting Fibonacci anyons in a Rydberg gas. *Phys. Rev. A* **86**, 041601 (2018).

46. Feiguin, A. et al. Interacting anyons in topological quantum liquids: The golden chain. *Phys. Rev. Lett.* **98**, 160409 (2007).

47. Trebst, S. et al. Collective states of interacting Fibonacci anyons. *Phys. Rev. Lett.* **101**, 050401 (2008).

48. Chandran, A., Schulz, M. D. & Burnell, F. J. The eigenstate thermalization hypothesis in constrained Hilbert spaces: a case study in non-Abelian anyon chains. *Phys. Rev. B* **94**, 235122 (2016).

49. Lan, Z. & Powell, S. Eigenstate thermalization hypothesis in quantum dimer models. *Phys. Rev. B* **96**, 115140 (2017).

50. Chandran, A., Burnell, F. J. & Sondhi, S. L. Absence of Fibonacci anyons in Rydberg chains. *Phys. Rev. B* **101**, 075104 (2020).

51. Lan, Z., van Horssen, M., Powell, S. & Garrahan, J. P. Quantum slow relaxation and metastability due to dynamical constraints. *Phys. Rev. Lett.* **121**, 040603 (2018).

52. Smith, A., Knolle, J., Moessner, R. & Kovrizhin, D. L. Absence of ergodicity without quenched disorder: From quantum disentangled liquids to many-body localization. *Phys. Rev. Lett.* **119**, 176601 (2017).

53. Brenes, M., Dalmonte, M., Heyl, M. & Scardicchio, A. Many-body localization dynamics from gauge invariance. *Phys. Rev. Lett.* **120**, 030601 (2018).

54. Surace, F. M. et al. Lattice gauge theories and string dynamics in Rydberg atom quantum simulators. Preprint at http://arxiv.org/abs/1902.09551 (2019).

55. Magnifico, G. Real time dynamics and confinement in the $Z_n$ Schwinger-Weyl lattice model for 1+1 QED. Preprint at http://arxiv.org/abs/1909.04821 (2019).

56. Görg, F. et al. Realization of density-dependent Peierls phases to engineer quantized gauge fields coupled to ultracold matter. *Nat. Phys.* **15**, 1161–1167 (2019).

57. Schweizer, C. et al. Floquet approach to $Z_2$ lattice gauge theories with ultracold atoms in optical lattices. *Nat. Phys.* **15**, 1168–1173 (2019).

58. Abanin, D. A., Altman, E., Bloch, I. & Serbyn, M. Colloquium: many-body localization, thermalization, and entanglement. *Rev. Mod. Phys.* **91**, 021001 (2019).

59. Binder, K. & Young, A. P. Spin glasses: experimental facts, theoretical concepts, and open questions. *Rev. Mod. Phys.* **58**, 801–976 (1986).

60. Berthier, L. & Biroli, G. Theoretical perspective on the glass transition and amorphous materials. *Rev. Mod. Phys.* **83**, 587–645 (2011).
61. Biroli, G. & Garrahan, J. P. Perspective: the glass transition. *J. Chem. Phys.* **138**, 12A301 (2013).

62. Fredrickson, G. H. & Andersen, H. C. Kinetic Ising model of the glass transition. *Phys. Rev. Lett.* **53**, 1244–1247 (1984).

63. Palmer, R. G., Stein, D. L., Abrahams, E. & Anderson, P. W. Models of hierarchically constrained dynamics for glassy relaxation. *Phys. Rev. Lett.* **53**, 958–961 (1984).

64. Carleo, G., Becca, F., Schiró, M. & Fabrizio, M. Localization and glassy dynamics of many-body quantum systems. *Sci. Rep.* **2**, 243 (2012).

65. De Roeck, W. & Huveneers, F. Asymptotic quantum many-body localization from thermal disorder. *Commun. Math. Phys* **332**, 1017–1082 (2014).

66. Schiulaz, M. & Müller, M. Ideal quantum glass transitions: many-body localization without quenched disorder *AIP Conf. Ser.* **1610**, 11 (2014).

67. Yao, N. Y., Laumann, C. R., Cirac, J. I., Lukin, M. D. & Moore, J. E. Quasi-many-body localization in translation-invariant systems. *Phys. Rev. Lett.* **117**, 240601 (2016).

68. Papić, Z., Stoudenmire, E. M. & Abanin, D. A. Many-body localization in disorder-free systems: the importance of finite-size constraints. *Ann. Phys.* **362**, 714–725 (2015).

69. van Horssen, M., Levi, E. & Garrahan, J. P. Dynamics of many-body localization in a translation-invariant quantum glass model. *Phys. Rev. B* **92**, 100305 (2015).

70. Veness, T., Essler, F. H. L. & Fisher, M. P. A. Quantum disentangled liquid in the half-filled Hubbard model. *Phys. Rev. B* **96**, 195153 (2017).

71. Smith, A., Knolle, J., Kovrizhin, D. L. & Moessner, R. Disorder-free localization. *Phys. Rev. Lett.* **118**, 266601 (2017).

72. Kim, I. H. & Haah, J. Localization from superselection rules in translationally invariant systems. *Phys. Rev. Lett.* **116**, 027202 (2016).

73. Yarloo, H., Langari, A. & Vaezi, A. Anyonic self-induced disorder in a stabilizer code: Quasi many-body localization in a translational invariant model. *Phys. Rev. B* **97**, 054304 (2018).

74. Michailidis, A. A. et al. Slow dynamics in translation-invariant quantum lattice models. *Phys. Rev. B* **97**, 104307 (2018).

75. van Horssen, M., Levi, E. & Garrahan, J. P. Dynamics of many-body localization in a translation-invariant quantum glass model. *Phys. Rev. B* **92**, 100305 (2015).

76. Oganesyan, V. & Huse, D. A. Localization of interacting fermions at high temperature. *Phys. Rev. B* **75**, 155111 (2007).

77. D'Alessio, L. & Rigol, M. Long-time behavior of isolated periodically driven interacting lattice systems. *Phys. Rev. X* **4**, 041048 (2014).

78. Mehta, M. L. *Random Matrices* Vol. 142 (Elsevier, 2004)

79. Alba, V., Fagotti, M. & Calabrese, P. Entanglement entropy of excited states. *J. Stat. Mech.* **2009**, P10020 (2009).

80. Eckholt, M. & García-Ripoll, J. J. Correlated hopping of bosonic atoms induced by optical lattices. *New J. Phys.* **11**, 093028 (2009).

81. Mazzarella, G., Giampaolo, S. M. & Illuminati, F. Extended Bose Hubbard model of interacting bosonic atoms in optical lattices: From superfluidity to density waves. *Phys. Rev. A* **73**, 013625 (2006).

82. Bissbort, U., Deuretbacher, F. & Hofstetter, W. Effective multibody-induced tunneling and interactions in the Bose-Hubbard model of the lowest dressed band of an optical lattice. *Phys. Rev. A* **86**, 023617 (2012).

83. Lühmann, D.-S., Jürgensen, O. & Sengstock, K. Multi-orbital and density-induced tunneling of bosons in optical lattices. *New J. Phys.* **14**, 033021 (2012).

84. Dutta, O. et al. Non-standard Hubbard models in optical lattices: a review. *Rep. Prog. Phys.* **78**, 066001 (2015).

85. Jürgensen, O., Meinert, F., Mark, M. J., Nägerl, H.-C. & Lühmann, D.-S. Observation of density-induced tunneling. *Phys. Rev. Lett.* **113**, 193003 (2014).

86. Baier, S. et al. Extended Bose-Hubbard models with ultracold magnetic atoms. *Science* **352**, 201–205 (2016).

87. Eckardt, A. Colloquium: Atomic quantum gases in periodically driven optical lattices. *Rev. Mod. Phys.* **89**, 011004 (2017).

88. Eckardt, A. et al. Exploring dynamic localization with a Bose-Einstein condensate. *Phys. Rev. A* **79**, 013611 (2009).

89. Meinert, F., Mark, M. J., Lauber, K., Daley, A. J. & Nägerl, H.-C. Floquet engineering of correlated tunneling in the Bose-Hubbard model with ultracold atoms. *Phys. Rev. Lett.* **116**, 205301 (2016).

90. Zhao, H., Knolle, J. & Mintert, F. Engineered nearest-neighbor interactions with doubly modulated optical lattices. *Phys. Rev. A* **100**, 053610 (2019).

91. Barbiero, L. et al. Coupling ultracold matter to dynamical gauge fields in optical lattices: From flux attachment to $Z_2$ lattice gauge theories. *Sci. Adv.* **5**, eaav7444 (2019).

92. Pancotti, N., Giudice, G., Cirac, J. I., Garrahan, J. P. & Bañuls, M. C. Quantum East model: localization, non-thermal eigenstates and slow dynamics. Preprint at https://arxiv.org/abs/1910.06616 (2019).

93. Zhao, K., Vovrosh, J., Mintert, F. & Knolle, J. Quantum many-body scars in optical lattices. *Phys. Rev. Lett.* **124**, 160604 (2020).

94. Sutherland, B. et al. Localization of electronic wave functions due to local topology. *Phys. Rev. B* **34**, 5208–5211 (1986).

95. Inui, M., Trugman, S. A. & Abrahams, E. Unusual properties of midband states in systems with off-diagonal disorder. *Phys. Rev. B* **49**, 3190–3196 (1994).

96. Heidrich-Meisner, F. et al. Quantum distillation: dynamical generation of low-entropy states of strongly correlated fermions in an optical lattice. *Phys. Rev. A* **80**, 041603 (2009).

97. Ronzheimer, J. P. et al. Expansion dynamics of interacting bosons in homogeneous lattices in one and two dimensions. *Phys. Rev. Lett.* **110**, 205301 (2013).

98. Greiner, M., Mandel, O., Hänsch, T. W. & Bloch, I. Collapse and revival of the matter wave field of a Bose-Einstein condensate. *Nature* **419**, 51–54 (2002).

99. Kollath, C., Läuchli, A. M. & Altman, E. Quench dynamics and nonequilibrium phase diagram of the Bose-Hubbard model. *Phys. Rev. Lett.* **98**, 180601 (2007).

100. Hudomal, A., Vasić, I., Regnault, N. & Papić, Z. Supporting data for Quantum scars of bosons with correlated hopping. https://doi.org/10.5518/810 (2020).

## Acknowledgements

The authors thank Thomas Iadecola for fruitful discussions. A.H. and I.V. acknowledge funding provided by the Institute of Physics Belgrade, through the grant by the Ministry of Education, Science, and Technological Development of the Republic of Serbia. N.R. was supported by the Department of Energy Grant No. DE-SC0016239, the National Science Foundation EAGER Grant No. DMR 1643312, Simons Investigator Grant No. 404513, ONR Grant No. N00014-14-1-0330, the Packard Foundation, the Schmidt Fund for Innovative Research, and a Guggenheim Fellowship from the John Simon Guggenheim Memorial Foundation. Z.P. acknowledges support by EPSRC grant EP/R020612/1 and the National Science Foundation under Grant No. NSF PHY-1748958. Part of the numerical simulations were performed on the PARADOX-IV supercomputing facility at the Scientific Computing Laboratory, National Center of Excellence for the Study of Complex Systems, Institute of Physics Belgrade. The authors would also like to acknowledge the contribution of the COST Action CA16221.

## Author contributions

A.H., I.V., N.R., and Z.P. contributed to developing the ideas, analyzing the results and writing the manuscript. A.H. performed the calculations and designed the figures.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s42005-020-0364-9.

Correspondence and requests for materials should be addressed to A.H.

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_15.png" /></figure>

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020