Observation of many-body scarring in a Bose--Hubbard quantum simulator

The ongoing quest for understanding nonequilibrium dynamics of complex quantum systems underpins the foundation of statistical physics as well as the development of quantum technology. Quantum many-body scarring has recently opened a window into novel mechanisms for delaying the onset of thermalization by preparing the system in special initial states, such as the $\mathbb{Z}_2$ state in a Rydberg atom system. Here we realize many-body scarring in a Bose-Hubbard quantum simulator from previously unknown initial conditions such as the unit-filling state. We develop a quantum-interference protocol for measuring the entanglement entropy and demonstrate that scarring traps the many-body system in a low-entropy subspace. Our work makes the resource of scarring accessible to a broad class of ultracold-atom experiments, and it allows one to explore the relation of scarring to constrained dynamics in lattice gauge theories, Hilbert space fragmentation, and disorder-free localization.


I. INTRODUCTION
Coherent manipulation of quantum many-body systems far from equilibrium is key to unlocking outstanding problems in quantum sciences including strongly coupled quantum field theories, exotic phases of matter, and development of enhanced metrology and computation schemes. These efforts, however, are frequently plagued by the presence of interactions in such systems, which lead to fast thermalization and information scrambling: the behavior known as quantum ergodicity [1][2][3]. A twist came with recent advances in synthetic quantum matter, which enabled detailed experimental study of thermalization dynamics in isolated quantum many-body systems, leading to the observation of ergodicity-violating phenomena in integrable [4] and many-body localized systems [5,6].
More recently, quantum many-body scarring has emerged as another remarkable ergodicity-breaking phenomenon, where preparing the system in special initial states effectively traps it in a "cold" subspace that does not mix with the thermalizing bulk of the spectrum [7,8]. Such behavior hinders the scrambling of information encoded in the initial state and suppresses the spreading of quantum entanglement, allowing a many-body system to display persistent quantum revivals. Many-body scarring was first observed in the Rydberg atom experimental platform [9,10], and subsequent observations of weak ergodicity breaking phenomena have attracted much attention [11][12][13]. On the other hand, theoretical works have unearthed universal scarring mechanisms [14][15][16][17], pointing to the ubiquity of scarring phenomena in periodically driven systems [18][19][20] and in the presence of disorder [21,22]. Given that many-body scarring in Rydberg atom systems has previously been reported in a single initial state, the Z 2 -ordered state, many questions remain about the overall fragility of this phenomenon and its sensitivity to the initial condition. It is thus vital to extend the realm of scarring to a greater variety of experimental platforms and more accessible initial conditions, which would empower a fundamental understanding of nonergodic dynamics in various research areas ranging from lattice gauge theories to constrained glassy systems.
In this paper, we observe many-body scarring in a largescale Bose-Hubbard quantum simulator, where we employ a tilted optical lattice to emulate the PXP model, a canonical model of many-body scarring [23][24][25][26]. We demonstrate that many-body scarring can result from a larger set of initial states, including the unit-filling state at finite detuning, hitherto believed to undergo fast thermalization [9]. Furthermore, we demonstrate that periodic driving can be used to enhance scarring behavior. Taking advantage of spin-dependent optical superlattices, we measure the system's entanglement entropy by interfering identical copies in the double wells. We show the average entropy of single-site subsystems to be a good approximation of half-chain bipartite entropy, revealing a key property of scarring: the "trapping" of the quantum system in a low-entropy subspace, which prevents its relaxation into the exponentially large Hilbert space.
The remainder of this paper is organized as follows. In Sec. II we introduce our experimental setup and show how it can realize the PXP model. In Sec. III we benchmark our quantum simulation by observing many-body scarring from the previously known Z 2 initial state. We also demonstrate the enhancement of scarring under periodic driving. In Sec. IV we present our measurements of entanglement entropy, providing deeper insight into the slow thermalization dynamics associated with scarred initial states. Finally, in Sec. V we extend the scarring phenomenon to a regime at moderate detuning for the unit-filling initial state. Our conclusions are presented in Sec. VI, while Appendixes A-D contain the derivation of the PXP mapping, further details on state preparation and measurement techniques, and a numerical study of other scarred initial conditions.

II. MAPPING THE PXP MODEL ONTO THE BOSE-HUBBARD MODEL
The PXP model [27,28] describes a kinetically constrained chain of spin-1/2 degrees of freedom. Each spin can exist in two possible states, |• and |• , corresponding to the ground state and excited state, respectively. An array of N such spins is governed by the Hamiltonian whereX = |• •| + |• •| is the Pauli x matrix, describing local spin precession with frequency . The projectors onto the ground state,P = |• •|, constrain the dynamics by allowing a spin to flip only if both of its neighbors are in the ground state.
A remarkable property of the PXP model is that it is quantum chaotic and yet it exhibits persistent quantum revivals from a highly out-of-equilibrium |Z 2 ≡ |•••• · · · initial state [23,[29][30][31]. The presence of revivals from a special initial state in an overall chaotic system was understood to be a many-body analog of the phenomena associated with a single particle inside a stadium billiard, where nonergodicity arises as a "scar" imprinted by a particle's classical periodic orbit [16,32,33]. In many-body scarred systems, eigenstates were shown to form tower structures [23]. These towers are revealed by the anomalously high overlap of eigenstates with special initial states, and their equal energy spacing is responsible for quantum revivals. While previous experiments on Rydberg atoms [9,10] have primarily focused on the |Z 2 initial state, we will demonstrate that the PXP model can effectively emerge in the Bose-Hubbard model, allowing us to identify scarred revivals from a larger set of initial conditions, including the polarized state |0 ≡ |••• · · · .
Our experiment begins with a 87 Rb Bose-Einstein condensate, which is compressed in the z direction and loaded into a single layer of a pancake-shaped trap. We then perform the superfluid-to-Mott-insulator phase transition with optical lattices in the x-y plane. In both x and y directions, we have a superlattice that is formed by superimposing the "short" lattice, with a s = 383.5 nm spacing, and the "long" lattice, with a l = 767 nm spacing [34,35], each of which can be individually controlled. We realize independent one-dimensional (1D) Bose-Hubbard systems in the y direction by ramping up the short-lattice depth in the x direction over 40E r , with E r = h 2 /8ma 2 s being the short-lattice recoil energy, where h is the Planck constant and m is the 87 Rb atomic mass. The short lattice in the y direction makes an approximately 4 • angle with gravity, which results in a static linear tilt per site of g = 816 Hz; see Fig. 1(a). An external magnetic field gradient B may be further added to create a tunable linear tilting potential = g + B . The effective Hamiltonian describing our simulator iŝ where J is the hopping amplitude,b andb † are the standard Bose annihilation and creation operators, the interaction energy isÛ = (U/2) L i=1n i (n i − 1), and the tilt potential iŝ = L i=1 in i . L denotes the number of sites in the chain with open boundary conditions, and we restrict ourselves to a total filling equal to 1, i.e., with the same number of bosons as lattice sites.
In order to realize the PXP model in the Bose-Hubbard quantum simulator, we tune the parameters to the resonant regime U ≈ J [36,37], which has been studied extensively in the context of quantum Ising chains [38][39][40]. In this regime, three-boson occupancy of any site is strongly suppressed, and doublons can only be created by moving a particle to the left, e.g., · · · 11 · · · → · · · 20 · · · , or destroyed by moving a particle to the right. The states of the PXP model are understood to live on the bonds of the Bose-Hubbard model. An excitation in the PXP model • j, j+1 , living on the bond ( j, j + 1), corresponds to the creation of a doublon 2 j 0 j+1 on site j in the Bose-Hubbard chain. We identify the unit-filling state |111 · · · with the PXP polarized state, |0 . Any other configuration of the PXP model can be mapped to a Fock state in the Bose-Hubbard model by starting from the unit filling, identifying the bonds that carry PXP excitations and replacing the corresponding sites in the Mott state with 11 → 20. Applying this rule across the chain allows us to map any basis state of the PXP model to a corresponding Fock state in the Bose-Hubbard model; for example, the |Z 2 state maps to the Fock state |· · · 2020 · · · . Figure 1  where the linear tilting potential is applied. Spin-dependent superlattices consisting of two standing waves in each direction can be individually controlled for state preparation and measurement. At the resonance U ≈ J, the dominant hopping process is 11 ↔ 20. The PXP excitations, •, live on the bonds between the lattice sites. The doublon configuration 20 in the Bose-Hubbard model maps to an excitation in the PXP model, while all other configurations are mapped to an empty site, •. For example, the given state |· · · •••••••• · · · maps to the Fock state |· · · 120201120 · · · . (b) Emergence of the PXP subspace in the Bose-Hubbard model at the resonance U ≈ J. Dots represent Fock states of the tilted Bose-Hubbard model with five bosons on five sites (restricting ourselves to at most three bosons on any site). Lines denote the allowed hopping processes. The color scale shows the sum of interaction and tilt energies Û +ˆ for each Fock state, and this value is conserved by resonant processes. The PXP dynamical subspace and its Fock states are explicitly labeled.

III. OBSERVATION OF Z 2 QUANTUM MANY-BODY SCARS
To prepare the initial states, we first employ an entropy redistribution cooling method [34] with the superlattice in the y direction to prepare ann = 2 Mott insulator in the left (odd) sites, while removing all atoms on the right (even) sites via site-dependent addressing [35]. This gives us the initial state |ψ 0 = |Z 2 = |2020 · · · (see Appendix B). In the region of interest, we have prepared 50 copies of the initial state |ψ 0 isolated by the short lattice along the x direction. Each copy extends over 50 short lattice sites along the y direction.
We quench the system out of equilibrium by abruptly dropping the y-lattice depth to 11.6E r , which corresponds to switching J from 0 to 51(1) Hz. This is done while simultaneously adjusting the lattice depth in the x and z directions accordingly, such that the interaction strength matches the linear tilt with U = = g ≈ 16J. After evolution time t, we freeze the dynamics by ramping up the y-lattice depth rapidly and read out the atomic density on the left ( n Left ) and right ( n Right ) sites of the double wells formed by the y superlattice successively [35,41]. This provides access to the density imbalance, M z = ( n Left − n Right )/( n Left + n Right ), an observable corresponding to the staggered magnetization in the PXP model; see Fig. 2(a). Another observable is the density of excitations in the PXP model, which is measured by projecting out the even atomic number occupancy on each site, then reading out the average odd particle density Pn ∈odd (1) [41]. Due to highly suppressed multiboson occupancy, we have P |• = n doublon (1) ≈ (1 − Pn ∈odd (1) )/2.
Away from the resonance, the dynamics is ergodic, and the staggered magnetization present in the initial |Z 2 state quickly decays with time; see Fig. 2(b). In contrast, by tuning to the vicinity of the resonance, = U , we observe distinct signatures of scarring: The system approximately undergoes persistent oscillations between the |Z 2 ≡ |•••• · · · configuration and its partner shifted by one site, |Z 2 ≡ |•••• · · · , as can be seen in the staggered magnetization profile and the density of excitations in Fig. 2(b). The density of excitations does not distinguish between |Z 2 and |Z 2 states; hence there is a trivial factor of 2 difference between the oscillation frequencies of P |• and M z .
The scarred oscillations in Fig. 2(b) are visibly damped with a decay time τ = 49.6 ± 0.8 ms. Nevertheless, as shown in Ref. [10], by periodically driving the system it is possible to "refocus" the spreading of the many-body wave function in the Hilbert space and thereby enhance the scarring effect, as we demonstrate numerically in Fig. 2(c) and experimentally in Fig. 2(d). Our driving protocol is based on modulating the laser intensity of the z lattice, which translates into periodic modulation of the interaction energy, U (t ) = + U 0 + U m cos(ωt ), while is kept fixed. This results in a modulation of the density of doublons in the chain, acting as the analog of the chemical potential in the PXP model.
Numerical simulations of the PXP model with the driven chemical potential, shown in Fig. 2(c), demonstrate the dynamical stabilization of the Hilbert space trajectory. We visualize the trajectory by plotting the average sublattice occupations, n Left and n Right , normalized to the interval [0,1]. The |Z 2 and |Z 2 states are located at the coordinates (1,0) and (0,1), which are the lower right and upper left corners of this diagram, respectively. The polarized state |0 is at the origin (0,0).
In the undriven case [left panel of Fig. 2(c)], the trajectory at first oscillates between |Z 2 and |Z 2 states, while passing through a region with a lower number of excitations. However, as the time passes, the trajectory drifts, exploring progressively larger parts of the Hilbert space. By contrast, when the driving is turned on [right panel of Fig. 2 Observation of Z 2 quantum many-body scars in a Bose-Hubbard quantum simulator. (a) Starting from the state |ψ 0 = |· · · 2020 · · · -the analog of the |Z 2 state in the PXP model-we utilize gravity to provide linear tilt = g . We characterize quench dynamics by measuring density imbalance and the number of doublons, corresponding to staggered magnetization M z and density of excitations P |• in the PXP model. In the detuned regime − U ≈ −2J, the dynamics is ergodic, and the system has no memory of the initial state at late times. (b) Tuning to U ≈ , we observe persistent oscillations in both M z and P |• . This memory of the initial state is a signature of weak ergodicity breaking due to quantum many-body scars. trajectory approximately repeats the first revival period of the undriven case, even at late times. Thus the driving stabilized the scarred revivals without significantly altering their period.
Experimental measurements on the driven Bose-Hubbard model in Fig. 2(d) find a strong enhancement of the amplitude of the oscillations in staggered magnetization with the decay time τ increasing to 208 ± 10 ms, while the period remains nearly the same as in the static case. Optimal driving parameters were determined numerically using a combination of simulated annealing and brute-force search; see Supplemental Material [42].
We note that the experimental measurement of M z damps slightly faster than the theory prediction, shown by a curve in Fig. 2(b), at late times (t > 60 ms). We attribute this to an inherent residual inhomogeneity across the lattice, which results in dephasing between different parts of the system, as well as possible decoherence induced by scattering of the lattice lasers. To avoid the effect of these undesired dephasing or decoherence effects, in the following we limit our investigation to times up to 60 ms.

IV. UNRAVELING THE DETAILS OF SCARRED DYNAMICS VIA QUANTUM INTERFERENCE
Entanglement entropy is key for characterizing scarring behavior. Entropy provides a window into the evolution of the system's wave function and the spreading of quantum entanglement. For a system trapped in a scarred subspace, thermalization is inhibited, and the system exhibits suppressed entropy growth and periodic fidelity revivals. Measuring these observables usually requires brute-force state tomography, but for our 50-site Bose-Hubbard system with a Hilbert space dimension exceeding 10 28 , this approach is generally impossible.
However, the superlattice in the x direction allows us to probe entanglement entropy by interfering identical copies in the double wells, analogous to the 50 : 50 beam splitter (BS) interference employed in photonics experiments [43]; see Fig. 3(a). This is done by freezing the dynamics along the chains in the y direction after evolution time t; then we interfere copies of |ψ (t ) in the double wells formed by the x superlattice (see Appendix C). After the interference, a parity projection helps read out the average odd particle FIG. 3. Probing many-body scarred dynamics via quantum interference. (a) After evolution time t, we freeze the dynamics in the y direction, and then by interfering two identical copies in the double wells along the x direction, we obtain the second-order Rényi entropy. (b) The entropy for a single site, S (1) , is seen to have robust oscillations with the same frequency as in Fig. 2(b), indicating a lack of thermalization. The slow growth of entropy in the absence of driving (upper panel) is strongly suppressed when we drive the system using the same parameters as in Fig. 2(d) (lower panel). In both cases, the single-site entropy is a good approximation to the half-chain entropy, S L/2 , evaluated numerically using TEBD (gray curve).

density P BŜ
n∈odd (1) , which gives us access to the second-order Rényi entropy [44]. We measure the entropy of single-site subsystems S (1) = −ln(Tr (1) [ρ(t ) 2 ]) = −ln(1 − 2 P BŜ n∈odd (1) ), whereρ(t ) = |ψ (t ) ψ (t )| is the density matrix. Entanglement entropy S (1) , shown in Fig. 3(b), grows much more slowly than expected in a thermalizing system. The growth is accompanied by oscillations with the same frequency as P |• in Fig. 2(b), implying that the system returns to the neighborhood of product states |Z 2 and |Z 2 . Furthermore, the entropy growth becomes almost fully suppressed by periodic driving, indicating that the scarred subspace disconnects from the thermalizing bulk of the spectrum. Numerical timeevolving block decimation (TEBD) simulations confirm that this lack of thermalization at the single-site level provides a good approximation for the behavior of larger subsystems, as demonstrated by the half-chain bipartite entropy S L/2 plotted in Fig. 3(b). This shows that scarring traps the system in a vanishingly small corner of an exponentially large Hilbert space.

V. EMERGENCE OF DETUNED SCARRING IN THE POLARIZED STATE
Up to this point, we have provided extensive benchmarks of our quantum simulator against the previously known case of Z 2 quantum many-body scars [9]. In this section we demonstrate that our quantum simulator also hosts distinct scarring regimes for initial states other than |Z 2 , which are enabled by detuning and further stabilized by periodic drive. We highlight this finding by observation of scarring behavior in the polarized state |0 , previously not associated with scars.
We first prepare the unit-filling state |1111 · · · by transferring |2, 0 to |1, 1 states in the superlattice [34], which maps to the polarized state in the PXP model (see also Appendix B). In the absence of detuning or periodic drive, we observe fast relaxation: Both the density of excitations and single-site entropy rapidly relax, with no visible oscillations beyond the timescale ∼1/J; see Fig. 4(a). Interestingly, when we bias the system by a static detuning, U 0 = −2.38J, we observe the emergence of oscillations in all three observables, accompanied by a slight decay; see Fig. 4(b). Finally, if we also periodically modulate the interaction with amplitude U m = 1.54J and frequency ω = 4.9J × 2π , we find a dramatic enhancement of scarring [ Fig. 4(c)]. In particular, entropy now shows pronounced oscillations, signaling robust scar-induced coherence at all experimentally accessible times.
The intuitive picture behind our observations is summarized as follows. In the absence of detuning or periodic drive, the system initialized in the polarized state undergoes chaotic dynamics and rapidly explores the entire Hilbert space. By biasing the system via static detuning, thermalization can be suppressed over moderate timescales. Finally, by periodically driving the system it is possible to "refocus" the spreading of the many-body wave function in the Hilbert space and thereby enhance the scarring effect, similar to the findings of Ref. [10] for the |Z 2 state. In the remainder of this section, we present our theoretical analysis of the experiment that supports this interpretation of the dynamics. Figure 5 shows the results of exact diagonalizations of the PXP model in Eq. (1) in the presence of static detuning, H (μ) =Ĥ PXP + μ 0 in i , wheren i takes a value equal to 1 if site i contains an excitation and 0 otherwise. The static chemical potential μ 0 is proportional to the Bose-Hubbard detuning parameter U 0 in Fig. 4. Figure 5(a) plots the overlap of all energy eigenstates |E of the pure PXP model (μ 0 = 0) with the polarized state |ψ 0 = |0 . As expected, we do not see any hallmarks of scars, such as ergodicity-violating eigenstates with anomalously enhanced projection on |0 . Moreover, the lowest-entropy eigenstates, denoted by squares in Fig. 5(b), are the known Z 2 scarred eigenstates [31] which are hidden in the bulk of spectrum when the overlap is taken with the |0 state.
On the other hand, when we add the static chemical potential μ 0 = 1.68 , corresponding to the detuning value in Fig. 4, a band of scarred eigenstates with anomalously large overlap with |0 emerges; see Fig. 5(c). The band of scarred eigenstates, highlighted by star symbols in Fig. 5(c), spans the entire energy spectrum, but their support on |0 is biased towards the ground state due to the breaking of particle-hole symmetry by detuning. The detuned scarred states also have anomalously low entanglement entropy, as seen in Fig. 5(d).
A few comments are in order. We note that exact diagonalization confirms that the PXP model remains chaotic for the value of detuning used in Fig. 5(c), and this detuning is not large enough to trivially fragment the entire spectrum into disconnected sectors with the given numbers of excitations [42]. Moreover, we confirmed that the scarred eigenstates in Fig. 5(c) are distinct from the ones associated with the |Z 2 state in Fig. 5(a). Thus it remains to be understood if these eigenstates can be described within the su(2) spectrum-generating algebra framework developed for the |Z 2 state in Ref. [45].
Nevertheless, similar to the |Z 2 case, the scarring from the |0 state can be further enhanced by periodic modulation of the PXP chemical potential, μ(t ) = μ 0 + μ m cos(ωt ). By evaluating the corresponding Floquet operator, we find that a single Floquet mode develops a very large overlap with the |0 state [42]. The existence of a single Floquet mode, whose mixing with other modes is strongly suppressed, gives rise to robust oscillations in the dynamics well beyond the experimentally accessible timescales.
To probe the ergodicity of the dynamics from the polarized state, we compare the difference between the predictions of the diagonal and canonical ensembles for an observable such as the average number of excitations; see Fig. 6(a). These two ensembles are expected to give the same result if the strong eigenstate thermalization hypothesis (ETH) holds [46] and all eigenstates at a similar energy density yield the same expectation value for local observables. Figure 6(a) shows that the discrepancy between the two ensembles is the strongest around μ 0 ≈ 1.68 , where we observe strong scarring. For μ 0 / close to 0, the polarized state For the optimal value μ 0 = 1.68 , identified in (a) and also used in experiment, the trajectory avoids the high-leakage region and approximates well the quantum dynamics, while it is not limited to a small corner of the many-body Hilbert space. thermalizes quickly towards the thermal value expected for a state whose expectation value of the energy is near the middle of the many-body spectrum. For very large μ 0 / , we enter a trivial regime where the polarized state is close to the ground state and only a few eigenstates at low energies are relevant for the dynamics. Hence, in this regime, quenching from the polarized state is similar to quenching from a thermal state at a very low temperature, and the agreement between the two ensembles is again very good. However, in this regime, only a very small part of the many-body Hilbert space is explored by the dynamics. This is not the case in the scarred regime that we investigate experimentally, and this can be demonstrated by studying the relevant classical limit, as shown next.
In the single-particle case, scarred quantum dynamics originates from an unstable periodic orbit in the classical limith → 0 [47]. In a many-body system, one approach to establishing a quantum-classical correspondence is to project the Schrödinger dynamics into a variational manifold, e.g., spanned by matrix product states [48], a method known as the "time-dependent variational principle" (TDVP). It was shown that the scarred dynamics of the |Z 2 state in the PXP model is well captured by the TDVP approach, allowing one to identify a classical orbit [16]. In Figs. 6(b) and 6(c) we utilize the TDVP approach to gain a semiclassical understanding of the detuned scarred dynamics from the |0 state. We parametrize the TDVP manifold using translation-invariant, spin-coherent states compatible with the Rydberg blockade constraint [49]. The states are defined by the Bloch sphere angles θ and φ, where sin(θ ) is proportional to the density of excitations, while φ describes the phase. In the thermodynamic limit, we can obtain classical equations of motion for θ and φ (see Ref. [50] for a detailed derivation). Figure 6(b) demonstrates that this classical dynamical system provides an excellent approximation of the quantum trajectory for sufficiently large values of μ 0 , including μ 0 = 1.68 .
To quantify the accuracy of the TDVP approach in capturing the quantum dynamics, we use "quantum leakage": the instantaneous norm of a component of the state vector that lies outside the TDVP manifold, γ 2 ≡ (1/N )|| |ψ − iĤ |ψ || 2 [16]. For the initial state |0 , the leakage has a simple analytic expression γ 2 = 2 sin 6 θ/(1 + sin 2 θ ) [50]. The leakage is higher as θ is increased, corresponding to a larger density of excitations. In this regime, i.e., for small values of μ 0 / , the PXP constraint has a strong effect, and the spin-coherent state ansatz does not faithfully capture the dynamics. On the other hand, for large values of μ 0 / , the leakage is low, but θ is confined to values near zero; thus the trajectory does not explore much of the Hilbert space. This corresponds to the trivial case where the dynamics is confined to very low densities of excitations, rendering the constraint unimportant. Finally, in the intermediate regime of μ 0 / where we observe the scarring, the TDVP dynamics is able to "avoid" the highleakage area, as seen in Fig. 6(c), while at the same time θ is not pinned to zero and the dynamics is not confined to one corner of the Hilbert space.

VI. DISCUSSION AND OUTLOOK
We performed a quantum simulation of the paradigmatic PXP model of many-body scarring using a tilted Bose-Hubbard optical lattice. We demonstrated the existence of persistent quantum revivals from the |Z 2 initial state and their dynamical stabilization, opening up a route for the investigation of scarring beyond Rydberg atom arrays. By harnessing the effect of detuning, we observed a scarring regime associated with the polarized initial state. As the latter state is spatially homogeneous, its preparation does not require a superlattice, which makes further investigations of scarring phenomena accessible to a large class of ultracold-atom experiments.
Moreover, we have demonstrated that periodic driving can lead to a striking decoupling of the scarred subspace from the rest of the thermalizing bulk of the spectrum, as revealed by the arrested growth of entanglement entropy. The mechanism of this enhancement is a subject of ongoing investigations. On the one hand, Ref. [51] used a kicked toy model to argue that the scarring enhancement originates from a discrete time crystalline order. On the other hand, Ref. [52] studied the cosine drive employed in experiment, finding two distinct regimes with long-lived scarred revivals. In the regime corresponding to the parameter values in Fig. 2 above, the driving parameters need to be fine-tuned to match the intrinsic revival frequency of the undriven scarred system. Moreover, the stabilization was no longer possible when the system was perturbed by terms which destroy scarring in the undriven case. This suggests that driving indeed acts as an enhancement mechanism, preventing dynamics from "leaking" into the thermalizing bulk.
Our demonstration of scarring in the |0 state highlights the importance of energy density. While the |Z 2 has predominant support on the eigenstates in the middle of the spectrum, i.e., it constitutes an "infinite temperature" ensemble, the support of the |0 state is biased towards one end of the spectrum as a result of particle-hole symmetry breaking via the detuning potential. This suggests that, depending on the effective temperature, one can realize scarring from a much larger class of initial states with a suitable choice of detuning and periodic driving protocols. We illustrate this in Appendix D by simulating the quench of the chemical potential in the PXP model (see also Ref. [50]).
The versatility of optical lattice platforms allows one to directly probe the link between many-body scarring and other forms of ergodicity-breaking phenomena, such as Hilbert space fragmentation and disorder-free localization, as the latter can be conveniently studied in our setup by varying the tilt. In this context, we note that Ref. [12] has recently used the tilt potential to demonstrate Hilbert space fragmentation in the Fermi-Hubbard optical lattice. By contrast, in this paper we explored ergodicity breaking due to many-body scars occurring within a single fragment of the Hilbert space. While many-body scarring can be induced in the Fermi-Hubbard model by tuning to a similar resonance condition [53], the underlying mechanism is an approximate dimerization of the chain, which is conceptually different from the PXP-type scarring considered here.
In future work, it would be interesting to explore realizations of new scarring models by tuning to other resonance conditions and other types of lattices, including ladders and two-dimensional arrays. Indeed, it is known that the U(1) quantum link model (QLM) [54,55] can be exactly mapped to the PXP model [56]. As such, recent large-scale experiments realizing the U(1) QLM [57,58] can in principle also probe our results. A proposal has recently been introduced to extend these setups to (2 + 1)D [59], where a mapping between the U(1) QLM and PXP model does not hold, which would allow one to probe how the scarring regimes discovered in this paper would behave in higher spatial dimensions. Finally, the observation of long-lived quantum coherence due to scarring and its controllable enhancement via periodic modulation lays the foundation for applications such as quantum memories and quantum sensing [60,61]. The Hamiltonian describing our 1D Bose-Hubbard model is given in Eq. (2) of the main text, with J denoting the hopping amplitude,Ĥ U denoting the corresponding interaction term, andĤ denoting the tilt potential. We denote by L the number of lattice sites and assume open boundary conditions (OBCs). Unless specified otherwise, we fix the filling factor to ν = 1, i.e., the number of bosons is equal to the number of sites in the chain.
In the U, J limit, the energy spectrum of the Hamiltonian in Eq. (2) splits into bands with approximately constant expectation value of the diagonal terms, Ĥ U +Ĥ ≈ const, and the Hilbert space becomes fragmented. At the U ≈ J resonance, the only process which conserves Ĥ U +Ĥ is 11 ↔ 20, i.e., doublons can only be created by moving a particle to the left and destroyed by moving a particle to the right. In the connected component of the Fock state |111 · · · 111 , the system in the resonant regime is described by an effective Hamiltonian which results from the first-order Schrieffer-Wolff transformation applied to Eq. (2) [62]. In the Supplemental Material [42] we discuss the effect of higher-order terms in the Schrieffer-Wolff transformation.
Note that excitations live on the bonds between sites and this mapping also includes links to the two surrounding sites. For example, the configuration · · · 2020 · · · maps to • • In this equation, the index i labels the sites, while j labels the bonds between sites. The Kronecker delta functions have been expressed in terms of projectors,P j = |• j • j |, and the bosonic hopping terms correspond to the spin raising and lowering operators,σ ± j , on the bond j. We can use delta functions because there are no configurations with more than two particles per site in this connected component and the only possible values ofn i (2 −n i ) are 0 and 1. Moving a particle to the neighboring site on the left corresponds to creating an excitation, and moving to the right corresponds to annihilating, while delta functions act as constraints.
Finally, the effective Hamiltonian is equivalent to the PXP Hamiltonian when we set = − √ 2J and N = L − 1, withX j ≡ |• j • j | + |• j • j | denoting the usual Pauli x matrix. In the case of OBCs, the two boundary terms becomeX 1P2 and P N−1XN . Note that the effective bosonic model in Eq. (A3) for a system size L is equivalent to the PXP model for size N = L − 1 since the number of bonds is the number of sites minus 1.
In the PXP model, the initial states which lead to pronounced quantum revivals are the two states with the maximal number of excitations: the Néel states, |••••· · ·•• and |••••· · ·•• [23,31]. The equivalent states in the tilted Bose-Hubbard model are |2020 · · · 201 and |12020 · · · 20 , for odd system sizes, and |2020 · · · 20 and |120 · · · 201 for even sizes. In our experimental setup, it is not possible to exactly prepare the |2020 · · · 201 state due to the inability to independently control single sites. Instead, our experiment realizes the |2020 · · · 20 state, which corresponds to the Néel state |••••· · ·••• in the PXP model with an odd number of sites and OBCs. terms Ĥ U +Ĥ , as indicated by different colors. The inset shows the top part of the highest-overlap band, around the energy E = 202020201|Ĥ |202020201 = 432. This band is described by the effective Hamiltonian (A1), which preserves the expectation value Ĥ U +Ĥ and is equivalent to the PXP Hamiltonian. A band of scarred eigenstates is magnified in the inset, and indeed resembles similar plots for the PXP model [31]. As the two Néel states have the maximal number of doublons at filling factor ν = 1, this type of dynamics also leads to oscillations in doublon number, which was experimentally measured in Fig. 2.

APPENDIX B: STATE PREPARATION AND DETECTION
Our experiment starts out with a two-dimensional Bose-Einstein condensate of 87 Rb atoms prepared in the hyperfine state |↓ = 5S 1/2 |F = 1, m F = −1 . By applying a microwave pulse, atoms can be adiabatically transferred to the state |↑ = 5S 1/2 |F = 2, m F = −2 , which is resonant with the imaging laser and thus can be detected. The atoms are initially confined to a single layer of a pancake-shaped trap with 3 µm period. In both x and y directions, we have an optical superlattice that can be controlled separately. Each superlattice potential is generated by superimposing two standing waves with laser frequency λ s = 767 nm and λ l = 1534 nm, which can be described by s cos 2 (kx) − V x l cos 2 (kx/2 + θ x ), V (y) = V y s cos 2 (ky) − V y l cos 2 (ky/2 + θ y ), where V x(y) s(l ) is the depth of the short (long) lattice in the x (y) direction, k = 2π/λ s is the short-lattice wave number, and θ x(y) is the relative phase between the short and long lattices in the x (y) direction.
We first perform a cooling technique by loading the atoms into a staggered superlattice in the y direction at θ y = π/4, meanwhile ramping up only the short lattice in the x direction. We tune the y-superlattice potential to create a Mott insulator withn = 2 filling in odd sites, while even sites form an = 1.5 superfluid, serving as a reservoir for carrying away the thermal entropy [34].
Atoms at even sites are removed by performing siteselective addressing. This is done by first setting θ y = 0 to form double wells and then tuning the polarization of the short-lattice laser along the y direction to create an energy splitting between even and odd sites for the |↓ -to-|↑ transition. We transfer the atoms at even sites to |↑ and remove them with the imaging laser [35]. In this way we have prepared the initial |Z 2 state |2020 · · · . The same site-selective addressing procedure is also utilized to read out atomic density on even and odd sites separately in experiment. Inside each isolated double-well unit, we can perform state engineering that transfers the state |2, 0 to |1, 1 [34]. This results in the unit-filling state |1111 · · · which corresponds to the polarized state |0 in the PXP model.

APPENDIX C: QUANTUM INTERFERENCE IN THE DOUBLE WELLS
The beam splitter (BS) interference is realized in the balanced double wells formed by the superlattices in the x direction, expressed in Eq. (B1) by setting θ x = 0. In the noninteracting limit, indistinguishable bosonic particles coming into the interference at t = 0 interfere according to the bosonic bunching. Therefore having equal numbers of atoms coming into the two ports at t = 0 results in P BŜ n∈odd = 0 at t BS , while having different numbers of atoms interfering results in P BŜ n∈odd = 0.5. Each copy of atoms coming into the interference is prepared individually, and hence there is no global phase between them, resulting in the equivalence between the two output ports [44].
To implement the quantum-interference protocol, we quench the x-lattice potentials to V x s = 6E r and V x l = 5E r , resulting in the intra-double-well tunneling at J ≈ 740 Hz and inter-double-well tunneling J ≈ 35 Hz. Simultaneously, we lower the lattice depth in the x direction to 25E r and trapping frequency in the z direction to 1.4 kHz, achieving an interaction of U ≈ 360 Hz. Two examples are shown here in Fig. 8, where we interfere product states |1, 1 [ Fig. 8(a)] or |2, 0 [ Fig. 8(b)] in the double wells and read out the average odd particle density. At t BS = 0.14 ms we identify the beam splitter operation, where |1, 1 gives P BŜ n∈odd (1) = 0.01(3), while |2, 0 gives P BŜ n∈odd (1) = 0.48 (3). We simulate the interference dynamics with a 20-site chain consisting of ten double-well units. We find good agreement at later times, FIG. 9. Emergence of many-body scarring by quenching the chemical potential in the PXP model from μ i = −0.76 to μ f = 1.6 . (a) The dynamics of quantum fidelity (blue solid curve) is similar to that of the polarized state for the same value of μ f (red dashed curve). The overlap between the time-evolved state and |0 (black dash-dotted curve) shows that a significant state transfer occurs between them. (b) The overlap of the prequench ground state with the eigenstates ofĤ (μ f ) displays characteristic scar tower structures. Red crosses denote the highest overlaps with the |0 state in each scarred tower. The overlap of the prequench ground state with |0 and |Z 2 states is given in the inset. All data are for N = 32 spins in the zero-momentum, inversion-symmetric sector of the Hilbert space.
while the earlier times are affected by the finite time in the lowering and raising of lattice potentials, which takes 50 µs. The finite interaction strength and inter-double-well tunneling together result in about 1% error in the beam splitter operation in the simulation, but this is beyond the precision of our absorption imaging.

APPENDIX D: OTHER SCARRED STATES
In addition to the |Z 2 and |0 states, we find other reviving states in the PXP model with static detuning,Ĥ (μ), introduced in Sec. V. These initial states are the ground states ofĤ (μ i ), and they exhibit revivals when the detuning is quenched to a different value,Ĥ (μ i ) →Ĥ (μ f ). This setup generalizes the quench protocols studied in the main text. For example, setting μ i →−∞, the prequench ground state is simply the |Z 2 state, and then quenching to μ f = 0 (pure PXP model) gives rise to scarred many-body revivals. Conversely, if we set μ i →∞, the ground state is |0 , and quenching to μ f = 1.68 also leads to scarring, as this value corresponds to the Bose-Hubbard detuning value in Fig. 4.
We numerically identify similar scarring phenomenology in a larger set of initial conditions by varying the parameters μ i and μ f . In Fig. 9 we present an illustrative example for μ i = −0.76 and μ f = 1.6 . Unlike the |Z 2 and |0 states, the ground state ofĤ (μ i ), for general values of |μ i | < 2, is not a product state. Nevertheless, such ground states have low entanglement entropy and can be prepared experimentally, while at the same time they are nearly orthogonal to the |Z 2 and |0 states (the overlap with the latter is on the order 10 −5 ). We emphasize that this does not require fine-tuning: We find large regions of μ i and μ f leading to scarring.
The dynamics in Fig. 9 is similar to that of the polarized state evolved withĤ (μ = 1.68 ). During the evolution, the state periodically transfers to the polarized state and then returns to itself. The frequency of revivals is approximately the same as that for the polarized state evolved with the same static detuning μ f , but the revivals are more prominent. The overlap of theĤ (μ i ) ground state with all the eigenstates of H (μ f ) is shown in Fig. 9(b). These overlaps exhibit a similar pattern to the overlap of eigenstates with the polarized states (red crosses). Furthermore, the atypical eigenstates appear to be the same in the two cases, up to a difference in phase. This is similar to what we find for the |Z 2 and |Z 2 states at μ f = 0: Both states have the same magnitude of overlap with each eigenstate, while the phases are different.