Observation of unconventional many-body scarring in a quantum simulator

,

Coherent manipulation of quantum many-body systems far from equilibrium is key to unlocking outstanding problems in quantum sciences including stronglycoupled 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 manybody systems, leading to the observation of ergodicityviolating 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]. However, due to its high sensitivity to initial states and the fragility of quantum manybody systems, the experimental implementation of quantum many-body scars beyond the Z 2 state in Rydberg atom systems has remained elusive. 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 fundamental understanding of nonergodic dynamics in various research areas ranging from lattice gauge theories to constrained glassy systems.
In this work, 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], see Fig. 1a. We demonstrate that a combination of detuning and periodic driving results in an unconventional manybody scarring regime in the unit-filling state, see Fig. 1b, hitherto believed to undergo fast thermalization [9]. 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. Furthermore, we utilize the interference method to de- velop a protocol for measuring out-of-time correlations without brute-force tomography. We demonstrate this by measuring the single-site fidelity, which shows that scarring brings the system back to the vicinity of the initial state.

OBSERVATION OF Z2 QUANTUM MANY-BODY SCARS
We first benchmark our quantum simulator on Z 2 quantum many-body scars [9]. We start by realizing the PXP model [27,28], which describes a kinetically constrained chain of spin-1/2 degrees of freedom. Each spin can exist in two possible states, |• , |• 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. In our quantum simulator, the PXP model is realized by tuning the on-site Hubbard interaction U to approx-imately match the tilt potential ∆ of the optical latice [29,30], a regime that has been studied extensively in the context of quantum Ising chains [31][32][33]. In the resonant regime U ≈∆ J, where J is the hopping amplitude, three-boson occupancy of any site is strongly suppressed, and doublons can only be created by moving a boson to the left, e.g., . . . 11 . . . → . . . 20 . . ., or destroyed by moving a boson to the right, see Fig. 1a and Methods for details. The unit-filling state |111 . . . maps to the PXP polarized state |0 ≡ |••• . . . , while the state with doublons on every other site, |2020 . . . , corresponds to |Z 2 ≡ |•••• . . . state-the state containing the maximal number of excitations allowed by the constraint in the PXP model.
Remarkably, while the PXP model is quantum chaotic [23], preparing the system in a highly out-ofequilibrium |Z 2 initial state leads to persistent quantum revivals [34][35][36]. The presence of revivals due to 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,37,38]. In many-body scarred systems, eigenstates were shown to form tower structures, illustrated in Fig. 1b. These towers are revealed by the anomalously high overlap of eigenstates with the initial state, such as |Z 2 , and their equal energy spacing is responsible for quantum revivals.
We initialize the experiment by employing a superlattice in the y-direction to prepare an=2 Mott insulator in the left (odd) sites of the double wells with 99.2% fidelity, while removing all atoms on the right (even) sites via site-dependent addressing [39,40]. This gives us the initial state |ψ 0 = |Z 2 = |2020 . . . (see Methods). The superlattice is formed by super-imposing the "short" lattice, with a s = 383.5 nm spacing, and the "long" lattice, with a l = 767 nm spacing. 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. A pancake-shaped trap provides confinement in the zdirection. 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. An external magnetic field gradient can be generated with coils, creating a controllable linear tilting potential together with gravity.
After initial-state preparation, the atoms are isolated in deep lattices with 60E r . E r =h 2 /8ma 2 s is the shortlattice recoil energy, where h is the Planck constant and m is the 87 Rb atomic mass. 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 provided by gravity with U =∆ g ≈16J. After evolution time t, we  . . 2020 . . . , the analog of |Z2 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 ∆g − U ≈ − 2J, the dynamics is ergodic and the system has no memory of the initial state at late times. By constrast, tuning to ∆g≈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. (c)-(d) Periodic modulation of the interaction U (t) = ∆ + U0 + Um cos(ωt) with U0=1.85J, Um=3.71J, ω=3.85J further suppresses the spreading of the wave function in the Hilbert space, leading to an enhancement of scarring. In both the static and driven case, experimental data for M z and P |• are in excellent agreement with TEBD numerical simulations shown by blue and orange solid lines. Gray line shows the modulation U (t).
freeze the dynamics and ramp up the double wells, where we read out the atomic density on the left and right sites of the y-superlattice successively with in-situ absorption imaging [40,41]. This provides access to 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. 2a. Another observable is the density of excitations in the PXP model, which is measured by projecting the number occupancy on each site into even (odd) parity, then reading out the average odd particle density Pn ∈odd (1) [41]. Due to highly suppressed multi-boson 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. 2a. In contrast, in 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. 2b. 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. 2a are visibly damped with a decay rate γ, with γ −1 = 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 manybody wavefunction in the Hilbert space and thereby enhance the scarring effect, as we demonstrate in Fig. 2cd. Crucially, this can be achieved without significantly altering the period of revivals. 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. For the driving parameters in Fig. 2c we find a strong enhancement of the amplitude of the oscillations in staggered magnetization with γ −1 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 Supplementary Material (SM).
We note that the experimental measurement of M z damps slightly faster than the theory prediction, shown by a line in Fig. 2a, 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 scatterring of the lattice lasers. However, with driving we observe a significant increase in coherence time (Fig. 2c).
To avoid the effect of these undesired dephasing or decoherence effects, in the following we limit our investigation up to 60 ms.

UNRAVELING THE DETAILS OF SCARRED DYNAMICS VIA QUANTUM INTERFERENCE
Fidelity and entanglement entropy are key observables for characterizing scarring behavior. These observables provide a window to 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 these observables by interfering identical copies in the double wells, analogous to the 50 : 50 beam splitter (BS) interference employed in photonics experiments [42]; see Fig. 3a. 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 Methods). After the interference, a parity projection helps read out the average odd particle density P BŜ n∈odd (1) , which give us access to the second-order Rényi entropy [43]. Due to limitations in our imaging resolution, we were not able to access entropy in extensively large subsystems, however, we measured the entropy of single-site subsystems S (1) = − ln(Tr (1) [ρ(t) 2 ]) = −ln(1 − 2 P BŜ n∈odd (1) ), wherê ρ(t) = |ψ(t) ψ(t)| is the density matrix. Due to the strong suppression of entanglement growth in a scarred system, we find S (1) serves as a good approximation for the half-chain bipartite entropy S L/2 , see Fig. 3b. Both single-site and half-chain entropy are seen to oscillate with the same frequency as P |• , implying that the system returns to the neighborhood of product states |Z 2 and |Z 2 . Our measurement of the Rényi entropy allows to identify a driving protocol that almost fully disconnects the scarred subspace from the thermalizing bulk of the spectrum, trapping the system in a vanishingly small (1) ... Figure 3. Probing many-body scarred dynamics via quantum interference. (a)-(b) After evolution time t, we freeze the dynamics in the y-direction, then by interfering two identical copies in the double wells along the x-direction, we obtain the second order Rényi entropy for a single site, S (1) . The entropy is seen to have robust oscillations with the same frequency as in Fig. 2b, indicating a lack of thermalization at the single-site level. The single-site entropy is a good approximation to the half-chain entropy, S L/2 , evaluated numerically using TEBD (grey line). (c)-(d) We flip all right sites in the x-superlattice to the state |↑ , then by applying an external magnetic field gradient in the y-direction, we create state-dependent linear tilting with ∆ |↑ B = −2∆ |↓ B . Fine tuning of the Hubbard parameters freezes the dynamics of |↑ atoms, which are flipped back to |↓ after evolution time t. We then interfere the left and right copies to read out singlesite fidelity F (1) . The fidelity displays pronounced revivals, indicating periodic returns to the close vicinity of the initial state. Global fidelity, obtained numerically using TEBD, is shown by the grey line. Experimental data is for periodically driven systems with the same parameters as in Fig. 2d. corner of an exponentially large Hilbert space.
Furthermore, we extend the interference protocol to probe unequal-time correlators. Before initiating the evolution, we first transfer atoms on the right sites of the double wells in the x-superlattice to the internal state |↑ = |F =2, m F = − 2 while leaving the atoms on the left sites in the state |↓ = |F =1, m F = − 1 , see Methods for details. By applying the magnetic field gradient, we generate a state-dependent linear tilt with ∆ Fig. 3c. We fine-tune the magnetic field gradient, such that the |↑ atoms are effectively "frozen" in the initial state |ψ 0 , while evolving the chains along the left sites to |ψ(t) . We then interfere these copies of |ψ 0 and |ψ(t) in the x-superlattice, and thus read out the singlesite fidelity with F (1) = Tr (1) [ρ 0ρ (t)] = 1 − 2 P BŜ n∈odd (1) . In Fig. 3d, we observe that F 2 (1) displays persistent revivals at the frequency of M z , revealing the system's periodic return to the vicinity of its initial state. The single-site fidelity F (1) provides an upper bound for the global fidelity F = | ψ 0 | |ψ(t) | 2 for the quenches studied, and the two oscillate with the same frequency, as demonstrated numerically in Fig. 3d (see also SM).

EMERGENCE OF AN UNCONVENTIONAL SCARRING REGIME
Remarkably, we find that a combination of detuning and periodic drive can result in dynamical stabilization of a distinct scarring regime for initial states other than |Z 2 . We highlight this finding in experimental observations of scarring behavior in the polarized state |0 , previously not associated with quantum many-body scars [9].
We start with the unit-filling state |1111 . . . in the lattice along the y-direction [39], which maps to the polarized state in the PXP model (see also Methods). In the absence of detuning or periodic drive, we observe fast relaxation: the density of excitations, single-site entropy, and fidelity all rapidly relax, with no visible oscillations beyond the timescale ∼1/J, see Fig. 4a. 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. 4b. 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. 4c. In particular, both entropy and fidelity now show pronounced oscillations, signaling robust scarinduced coherence at all experimentally-accessible times.
Our experimental observations are explained by theoretical analysis of the PXP model summarized in Figs. 4d, e, f. By computing the overlap of all energy eigenstates |E of the PXP model in Eq. (1) with the polarized state |ψ 0 = |0 , we do not identify any hallmarks of scars, such as ergodicity-violating eigenstates with anomalously enhanced projection on |0 , see Fig. 4d. On the other hand, when we add static detuning, which corresponds to the chemical potential in the PXP model, a band of scarred eigenstates with anomalously large overlap with |0 emerges; see Fig. 4e. Numerical simulations show that the PXP model remains chaotic for the value of static detuning used in experiment (Fig. 4b), and this detuning is not large enough to trivially fragment the entire spectrum into disconnected sectors with the given numbers of excitations (see SM). The band of scarred eigenstates, illustrated by crosses in Fig. 4e, 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. We note that the scarred eigenstates in Fig. 4e are distinct from the known ones associated with the |Z 2 state, but otherwise possess similar scarring properties such as anomalously low entanglement entropy (see SM). Finally, further enhancement of scarring under periodic modulation of the PXP chemical potential is explained by the spectrum of the corresponding Floquet operator, numerically evaluated in Fig. 4f. We observe that a single Floquet mode (denoted by a cross) develops a very large overlap with the |0 state. 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.

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 new route for the investigation of scarring beyond Rydberg atom arrays. By harnessing the effect of detuning and periodic driving, we observed a distinct 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. The versatility of such platforms allows to directly probe the link between many-body scarring and other forms of ergodicitybreaking phenomena, such as Hilbert space fragmentation and disorder-free localization, as the latter can be conveniently studied in our setup by varying the tilt.
Our demonstration of scarring in the |0 state highlights the importance of energy density. While the |Z 2 initial state 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 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 (see SM). Furthermore, one may expect realizations of qualitatively new scarring models by tuning to other resonance conditions and other types of lattices, including ladders and two-dimensional  arrays. While in 1D our Bose-Hubbard lattice with the staggered chemical potential can also be mapped to the quantum link model providing a different PXP realization [44], such a mapping does not directly extend to higher dimensions and potentially results in yet another class of scarred models.
Finally, our methods for probing unequal-time correlators allow for state-of-the-art monitoring of nonequilibrium dynamics and its applications in quantum technology. Notably, this protocol can be used to probe the global fidelity, as well as the out-of-time correlations between arbitrary quantum states, e.g., | ψ(t 1 )| |ψ(t 2 ) | 2 , with the help of single-atom resolution quantum gas microscopes [43,45]. This would empower detailed experimental studies of exotic quantum phenomena such as dynamical quantum phase transitions [46]. 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 [47]. Moreover, the dynamical manipulation of a many-body system employed in this work can be directly used to prepare states with extensive multipartite entanglement [48], thus lending itself to novel protocols for phase estimation and quantum metrology.

Mapping between PXP and Bose-Hubbard models
The Bose-Hubbard model in a tilted optical lattice is described by the Hamiltonian where J is the hopping amplitude,b,b † are the standard Bose annihilation and creation operators, the interaction L denotes the number of sites in the chain and we assume open boundary conditions. Moreover, we fix the filling factor to ν=1, i.e., the total number of bosons is also equal to L.
In order to realize the PXP model in the Bose-Hubbard quantum simulator, we tune the parameters to the resonant regime U ≈∆ J [29,30]. 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 unitfilling, identifying the bonds that carry PXP excitations and replacing the corresponding sites in the Mott state with 11→20, see SM for further details. Applying this rule across the chain allows to map any basis state of the PXP model to a corresponding Fock state in the Bose-Hubbard model, e.g., the |Z 2 state maps to the Fock state |. . . 2020 . . . . Fig. M1 illustrates the profound change in the connectivity of the Fock space near the resonance U ≈∆ J, with en emergent dynamical subspace isomorphic to the PXP model in the sector containing the |Z 2 state.

State preparation and detection
Our experiment starts out with a two-dimensional Bose-Einstein condensate of 87 Rb atoms prapared 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 super-imposing two standing waves with laser frequency λ s =767 nm and λ l =1534 nm, which can be described by where V x(y) s(l) is the depth of short (long) lattice in x(y)direction, k=2π/λ s is the short lattice wave number, and θ x(y) the relative phase between the short and long lattices in 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 [39].
Atoms in even sites are removed by performing siteselective addressing. This is done by first setting θ y =0 to form double wells, 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 in even sites to |↑ and remove them with the imaging laser [40]. 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 [39]. This results in the unit-filling state |1111 . . . which corresponds to the polarized state |0 in the PXP model.

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. (3) by setting θ x =0. In the non-interacting limit, indistinguishable bosonic particles coming into the the interference at t=0 interfere according to the bosonic bunching. Therefore, equal number of atoms coming into the two ports at t=0 results in P BŜ n∈odd =0 at t BS , while different number of atoms interfering results in P BŜ n∈odd =0.5. Each copy of atoms coming into the interference are prepared individually, hence no global phase between them, resulting in the equivalence between the two output ports [43].
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 Starting from the initial |Z2 state, chains with |↑ atoms remain largely at Pn ∈odd (1) ≈0 indicating they remain "frozen" in the initial state, while |↓ atoms evolve with the scarred dynamics. (b) Starting from the initial polarized state |0 , chains with |↑ atoms remain largely at Pn ∈odd (1) ≈1, while the |↓ atoms exhibit scarred dynamics in the presence of static detuning.
kHz, achieving an interaction of U ≈360 Hz. Two examples are shown here in Fig. M2, where we interfere product states |1, 1 (Fig. M2a) or |2, 0 ( Fig. M2b) 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 10 doublewell units. We find good agreement at later times, while the earlier times are affected by the finite time in the lowering and rising of lattice potentials, which takes 50 µs. The finite interaction strength and inter-double-well tunneling results in about 1% error in the beam splitter operation in the simulation, but this is beyond the precision of our absorption imaging.

Internal-state-dependent evolution
We use site-selective addressing in the x-superlattice to transfer all right sites to the internal hyperfine state |↑ with a microwave pulse, while leaving all left sites in |↓ [40]. By applying an external magnetic-field gradient in the y-direction ∂ y B, the magnetic energy shift between nearest-neighbor short lattice sites of different hyperfine states can be described by the anomalous Zeeman effect where µ B is the Bohr magneton, the Landé factor g F is 1/2 for |↑ , and -1/2 for |↓ . Therefore, we create a state-dependent linear potential with ∆ |↑ B =−2∆ |↓ B . We fine-tune the gradient strength for different initial conditions such that the total linear potential ∆ U would avoid resonant processes while suppressing the direct tunneling. Hence, atoms in the |↑ state are "frozen" in the initial states during the entire evolution. To minimize the "leak" from the initial states, the Hubbard parameters are optimized experimentally by measuring the odd particle density Pn ∈odd (1) . For the initial |Z 2 state, Pn ∈odd (1) should remain close to 0, while for the polarized state |0 , Pn ∈odd (1) should remain close to 1. We find the optimum magnetic field gradient for |Z 2 state to be ∂ y B≈11.5 G/cm, corresponding to ∆ |↑ ≈200 Hz, meanwhile, the |↓ atoms feel an increased linear potential with ∆ |↓ = ∆ |↓ g + |∆ |↓ B |≈1120 Hz, where we adjust the interaction U and tunneling J such that ∆ |↓ ≈U ≈16J, which give rise to the scarred dynamics, see Fig. M3a. For the polarized state, however, ∂ y B≈14.3 G/cm and ∆ |↑ ≈50 Hz results in better freezing of the |↑ atoms, therefore we set U ≈∆ |↓ ≈1200 Hz, see Fig. M3b. In both cases we find around 3% "leaking" from the initial states within the investigation time we measure single-site fidelity. This leaking contributes to the small deviation in single-site fidelity we find between experiments and TEBD simulations in Fig. 3 and Fig. 4.
After evolution time t, chains with |↓ atoms have evolved to |ψ(t) , while chains with |↑ atoms have largely remained in |ψ 0 . We ramp up the short lattice in the y-direction to freeze the dynamics of |↓ atoms along the chains, and flip all |↑ atoms back to |↓ . Then interfere copies of |ψ 0 and |ψ(t) in the double wells formed by the x-superlattice. After the parity projection, we read out single-site fidelity with F (1) = Tr (1) where J is the hopping amplitude, U is the interaction strength, ∆ is the tilt and L the number of sites. 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.
∆ ≈ U resonance: first order terms In the U, ∆ J limit, the energy spectrum of the Hamiltonian in Eq. (S1) 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 can therefore is described by an effective Hamiltonian S2], as will be shown below (see also Ref. [S3] for the original derivation of the mapping and a recent review [S4]).
The connected component of the Hilbert space contains only certain types of two-site configurations (20,11,12,02,01), while all other two-site configurations are forbidden (22, 21, 10, 00). If we consider the configuration 20 to be an excitation, all allowed configurations can be mapped to those of the PXP model as follows: .
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 2 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 a neighboring site on the left corresponds to creating an excitation, moving to the right to annihilating, while the delta functions act as constraints.
Finally, the effective Hamiltonian is therefore equivalent to the PXP Hamiltonian when

Higher order terms in the mapping
The effective Hamiltonian of Eq. (S2) results from the first-order Schrieffer-Wolff transformation [S7] wherê H 0 =Ĥ U +Ĥ ∆ . In this section we look at the relevant terms that arise in the effective Hamiltonian at second order. To simplify the notation we write these terms as sums of range-3 operators, where |111 120| j denotes the operator changing sites j − 1, j and j + 1 from 120 to 111 while leaving all other sites unaffected.
First, we can identify the matrix elements that take the system out of the PXP sector. This happens by the appearance of sites with 3 bosons via the operator There are also off-diagonal matrix elements connecting states within the PXP sector, given bŷ There are also additional off-diagonal matrix elements connecting states outside of the PXP sector, but as they do not directly influence the dynamics out of it we do not describe them here.
Finally, the diagonal operator in this sector is given bŷ where the two-site operator |11 12| j acts on sites j and j + 1. As bulk terms get added, the overall diagonal factors are extensive in the system size in the Fock basis.
The state with the lowest on-site potential is |111 . . . 11 with a value of − (L−1)J 2 U . The maximum is ≈ 4J 2 L 3U for the state |120120 . . . 120 , which corresponds to the Z 3 state in the PXP terminology.
In order to see how these second-order terms change the effective model we can rewrite Eqs. (S7)-(S8) for the PXP model with N = L − 1 sites. To do this we introduce the single-site projector on the excited statê Q j = |• •| = 1 −P j . We then obtain and respectively. We notice that the off-diagonal correction has the form of a constrained XY term.

NUMERICAL METHODS
In the main text and this Supplementary Material, we use two types of numerical methods for modeling the experiment. For small Bose-Hubbard chains with L 12 sites, we use exact diagonalization techniques to obtain the full energy spectrum of the Hamiltonian and directly access the system's eigenstate properties. Unless specified otherwise, we restrict the occupancy of any site to be at maximum 3 bosons, as our results are found to be insensitive to allowing more than 3 bosons on any site.
To access dynamics in much larger systems, L 50 sites, we use TEBD variational method [S8] implemented in TenPy package [S9]. We employ the second order Trotter decomposition with time step 2.5 × 10 −5 /J and maximum bond dimension χ max = 3000. Such a small time step was necessary because some of the quantities we are interested in, e.g., the fidelity density, are sensitive to otherwise negligible fluctuations in the revival peak heights that appear for longer time steps.
For numerical simulations via either one of these methods, it is convenient to work in natural units =1. We adopt this convention for presenting all numerical results in this Supplementary Material.

Numerical demonstrations of the mapping between PXP and tilted Bose-Hubbard model
In this Section, we numerically corroborate the mapping between PXP and Bose-Hubbard models introduced in Sec. . Specifically, we use exact diagonalization to demonstrate the consistency between dynamics and eigenstate properties in the PXP model and the Bose-Hubbard model tuned to the resonance U =∆.
In analogy with the PXP model, the system initialized in the state 2020 . . . 201 is expected to oscillate between this state and the state 12020 . . . 20. This is not only the case for the effective model (S4) which is exactly equivalent to PXP, but also for the full tilted Bose-Hubbard model (S1) at the U = ∆ resonance, as can be observed in Fig. S1.  In Fig. S2 we show the evolution of the bipartite von Neumann entanglement entropy, S vN (t) = −Tr A (ρ A lnρ A ), whereρ A is the reduced density matrix for subsystem A of length L A . The system is initially prepared in the state 2020 . . . 201 or the completely homogeneous state 111 . . . 111. As in the PXP model, the entanglement entropy for the 2020 . . . 201 state exhibits slow and approximately linear growth in time. In contrast, the entanglement entropy for the state 111 . . . 111 rapidly saturates, implying that the system quickly thermalizes. The evolution of density imbalance between the even and odd sites M z = ( n odd − n even )/( n odd + n even ), which corresponds to staggered magnetization in the PXP model, is shown in Fig. S3. This is one of the quantities that was experimentally measured in the main text. Here we again compare the evolution with the full tilted Bose-Hubbard Hamiltonian (S1) and the effective Hamiltonian (S4), the latter being equivalent to the PXP model, and we find excellent agreement between the two.
The overlap of eigenstates in the full model in Eq. (S1) for a chain of L = 9 sites with the initial 202020201 state is given in Fig. S4. For parameter values U =∆=12, and J = 1, the energy spectrum is split into multiple bands with approximately constant expectation value of the sum of interaction and tilt terms Ĥ U +Ĥ ∆ , as indicated by different colors. The inset shows the top part of the highest-overlap band (around E = 202020201|Ĥ|202020201 = 432). This band is described by the effective Hamiltonian (S2), which preserves the expectation value Ĥ U +Ĥ ∆ and is equivalent to the PXP Hamiltonian. A band of L = 9 scarred eigenstates is visible in the inset, as expected from the analogy with the PXP model. These scarred eigenstates are responsible for the revival dynamics in Fig. S1. As the two Néel states have the maximal number of dou-blons at filling factor ν = 1, this type of dynamics also leads to oscillations in doublon number, which can be experimentally measured. As a side note, the system is also described by PXPlike effective models at other integer filling factors. The reviving initial states are of the form |(n + 1)(n − 1)(n + 1)(n − 1) . . . (n + 1)(n − 1)n for ν = n, e.g. |3131 . . . 312 for ν = 2 and |4242 . . . 423 for ν = 3, as shown in Fig. S5. Revival frequency increases with n as n(n − 1), but the revivals decay faster for larger n.

GLOBAL VS. SINGLE-SITE FIDELITY
In this section we discuss the relation between two measures of similarity of pure states used in this work. One measure is the quantum fidelity i.e., the global overlap between two pure states. This measure is very convenient for numerical simulations and theoretical analysis, but hard to measure in experiment. For this reason, we also consider a different measure consisting of an average of local measurements: is the density matrix obtained by performing the partial trace on all sites except sites j to j − 1 + r.
Both quantities, F and F (r) , are real and obey and F (r) (|φ , |ψ ) = F (r) (|ψ , |φ ), (S15) F(|φ , |ψ ) = F(|ψ , |φ ). (S16) It is also important to note that While for arbitrary states F (r) is neither an upper bound nor a lower bound of F, it does not mean that this is never the case. We are now limiting our study to the case where the state φ is a product state. The consequence of that is that the reduced density matrixρ φ j,j+r−1 will correspond to a pure state for any r. We can then chose a basis for each site such that |φ is a product of local basis states, and so a Fock basis state for the whole Hilbert space. Let us then denote the orthonormal states of this basis by |α .
This allows us to rewrite the reduced density matrix asρ φ j,j+r−1 = |φ j,j+r−1 φ j,j+r−1 |, where |φ j,j+r−1 corresponds to the state φ for sites j to j + r − 1 (remember that we can only do this because |φ is a product state).
This formulation implies the following simplification Tr |φ j,j+r−1 φ j,j+r−1 |ρ ψ j,j+r−1 where α|φ j,j+r−1 denotes the product of |α j,j+1−r and |φ j,j+1−r , which is always either 0 or 1. This means that α ψ 2 contributes to F (r) (|φ , |ψ ) with a weight of 1 N +1−r each time r consecutive sites are in the same state in |α and |φ . This simple rule allows us to not only derive this inequality between F and F (r) , but also to compare the effect of r on F (r) . Indeed, if with r a basis state |α has a weight of n L+1−r , then if n > 0 with r−1 it has a minimum weight of n+1 L+2−r . As n ≤ L+1−r, it implies that n L+1−r ≤ n+1 L+2−r . If n = 0 then for r − 1 the same state cannot contributes less, and so for any n it contributes more or the same amount. Hence we can conclude that Finally, it is important to note that the inequality n L+1−r ≤ n+1 L+2−r is saturated if and only if n = L + 1 − r, meaning that |α and |φ are the same. This is important as F (r) (|φ , |ψ ) is a weighted sum of all the α ψ 2 with weights equal or smaller to 1. In order for F (r) (|φ , |ψ ) to be equal to F (r−1) (|φ , |ψ ), all weights corresponding to a non-zero | α|ψ | 2 must stay the same. But the only weights that are not increasing are either the ones that stay equal to zero or the one of |α = |φ which stays equal to one. This implies that all inequalities of Eq. S19 are simultaneously saturated if and only if |ψ = |φ (in which case they are all equal to one) or F (1) (|φ , |ψ ) = 0 (in that case they are all equal to 0). It is also possible for some of them to be saturated. This can only happen if any consecutive subset of length m of |φ and |ψ are orthogonal, meaning the all F (r) are equal to 0 for all r ≥ m.
In the experimental setup we only have access to the single-site fidelity F (1) , which already mimics the behavior of the real fidelity F (see Fig. S6). While it bounds F from above, that bound is fairly loose. If we instead look at its square F 2 (1) , we can see that it approximates F much better as it takes a lower value when F is close to zero. While F 2 (1) is not guaranteed to be an upper bound of F, our theoretical simulations indicate that it still effectively acts as one for the conditions we study. Our simulations also show that the single-site second Renyi entropy S (1) shows a very similar behavior to the bipartite half-chain second Renyi entropy S. While it is limited in the range of values it can take, S (1) is clearly able to distinguish between the two regimes we are seeing in our setup: rapid entropy growth until a plateau is reached, and very slow entropy growth with oscillations on top.

ALTERNATIVE MAPPING BETWEEN PXP AND TILTED 1D BOSE-HUBBARD WITH STAGGERED DETUNING
There is another mapping between the tilted 1D Bose-Hubbard and the PXP model which is based on the additional staggered potential term added to the the model in Eq. (S1): The parameter δ determines the energy offset between even and odd lattice sites. The odd ones are now "plus" sites where the δ term is positive, while the even one are "minus" sites where it is negative. This means that for an odd chain of length L there are N = L−1 2 minus sites and N + 1 = L+1 2 plus sites. The model in Eq. (S20) has been experimentally studied (see [S10] and references therein). Its mapping to the U(1) quantum link model has been already established in the literature, as is the mapping between the U(1) quantum link model and the PXP model [S11]. However, the equivalence of the tilted Bose-Hubbard model with staggered detuning and the PXP model was to our knowledge never explicitly stated, so we will briefly explain it here. This mapping is valid in the regime U ≈ 2δ J and at filling factor ν = 1/2.
When U ≈ 2δ J, the second-order process 101 ↔ 020 becomes resonant. Nonzero tilt ∆ makes other relevant second-order processes such as 100 ↔ 001 offresonant. In this regime and for odd system size L with filling factor ν = L+1 2L , the effective Hamiltonian at second order of the model in Eq. (S20) is fragmented. One of these fragments can be mapped to the PXP model up to some diagonal boundary terms. To find the corresponding state in PXP, one only needs to look at the "minus" sites. Doublons on these sites are mapped to PXP excitations, as the resonant processes cannot create two doublons on two adjacent minus sites. Due to the nature of the resonant process they can never be singly-occupied, and empty minus sites are mapped to non-excited atoms. This means that the corresponding PXP model has length N = L−1 2 , which is just the number of minus sites.
The   S7 shows the wave function fidelity over time for the model in Eq. (S20) with U =2δ > ∆ J, and for the PXP model it can be mapped to. To directly compare the two models, we have rescaled the time axis by the constant κ, which takes the value √ 2J (2) = √ 2×4J 2 U /(U 2 −4∆ 2 ) in the Bose-Hubbard model and Ω in the PXP model. The staggered Bose-Hubbard model oscillates between two product states, 02000200 . . . 02001 and 100200020 . . . 0020, which are the analogs of the Néel states in the PXP model. Overall, the dynamics is seen to be very similar in the two models, with the slight difference between the two being likely due to the boundary terms defined in Eq. (S21) below.
To derive the mapping to the PXP model rigorously, we can separate the Hamiltonian in Eq. (S20) asĤ = H 0 −JV and perform the Schrieffer-Wolff transformation [S7]. HereĤ 0 encompasses all the diagonal terms while V simply corresponds to hopping, which assumes that δ, ∆, U J. Furthermore, we will only focus on the regime U = 2δ, in which case there are no first order terms. If U is close but not equal to 2δ, then the effective Hamiltonian at first order will contain diagonal terms proportional to |U − 2δ|. Finally, we only focus on the connected component of the second order Hamiltonian that can be mapped to the PXP model, meaning that resonant processes like 02010 ↔ 11001 are ignored as these configuration cannot appear in the Hilbert space component of interest.
In the relevant part of the Hilbert space, the only offdiagonal resonant process at second order is 101 ↔ 020, which appears with a weight of There are also two allowed second-order di-agonal processes in the bulk of the chain: 010 ↔ 010 and 020 ↔ 020. They have a respective weight of J (2) and 2J (2) . However, as creating a new doublon means emptying two singly-occupied sites, the diagonal matrix elements do not change under the off-diagonal process.
The only exception to this is hopping at the boundaries of the chain. For the leftmost site, only hopping to the right and then back is possible, leading to a contribution of 2J 2 U −2∆ instead of J (2) . For the rightmost site only hopping to the left is possible and the contribution is 2J 2 U +2∆ . This means that not all diagonal matrix elements are the same but they vary between J (2) N = L−1 2 J (2) and (N + 1)J (2) = L+1 2 J (2) . However the differences between the diagonal elements are O(1) and do not scale with L, so they become negligible for large system sizes. All together, the second order Hamiltonian can be mapped to the following model.
where Z denotes the usual Pauli z matrix [P j andX j were defined previously below Eq. (S5)]. It is worth mentioning that for U = 2δ, the third order effective Hamiltonian is identically zero due to the absence of diagonal elements in the perturbation V (which is simply the hopping) and the next correction only happens at fourth order.

EFFECT OF PERIODIC DRIVING ON Z2 SCARS
In this section we numerically explore the effect of driving on the stabilization of many-body scars and revival dynamics, both in the PXP and tilted Bose-Hubbard models.

Driven PXP model
Periodic driving has been shown to enhance and stabilize the revivals in the PXP model [S12, S13]. The optimal driving frequency was found to be close to twice that of revivals in the pure PXP model without driving. We consider the following spatially-uniform cosine driving scheme which was also experimentally implemented in Ref. [S12] Here, µ 0 is the static detuning, µ m modulation amplitude and ω driving frequency. For simplicity, when work-ing with the PXP model, we impose periodic boundary conditions (PBC). We determine the optimal values of µ 0 and µ m by scanning the parameter space for the highest time-averaged fidelity. To make sure we found the globally-optimal values of the driving parameters, we have also checked the results against the simulated annealing algorithm implemented in GNU Scientific Library. The evolution of quantum fidelity starting from the Néel state |Z 2 = |• • • • • • . . . can be seen in Fig. S8(a), both without driving (red) and driven with optimal driving parameters µ 0 /Ω = 1.15, µ m /Ω = 2.67 and ω/Ω = 2.72 (blue). Driving leads to high revivals whose amplitude remains close to 1 over very long times. Additionally, the driving also strongly suppresses the growth of entanglement entropy, as can be observed in Fig. S8(b).
In order to understand the mechanism of revival enhancement, in Fig. S9   The states on the main diagonal (except the Néel states) and in the triangle above it are forbidden due to the PXP constraints.
In the undriven case (purple), the trajectory at first oscillates between the Néel and anti-Néel states in the corners, while passing through a region with a lower number of excitations (bottom left). However, as the time progresses, the wavefunction starts to thermalize and the trajectory drifts towards towards the average numbers of excitations ( n A , n B ) = (0.276, 0.276). When the driving is turned on (green), the trajectory continues to approximately repeat the first revival period in the undriven case and does not seem to thermalize even at very late times. In this way the revivals are stabilized and enhanced. Another effect of driving is that the overlap with the anti-Néel state is now lower, but its peaks do not decay with time.
Finally, in Fig. S10 we studied the Floquet modes of the driven PXP model. The Floquet modes are a generalization of eigenstates for periodic time-dependent HamiltoniansĤ(t + 2π ω ) =Ĥ(t). Unlike the eigenstates, the Floquet modes evolve in time, but they are timeperiodic with the same periodicity as the driven Hamltonian, Φ n (t + 2π ω ) = Φ n (t). We have computed all the Floquet modes Φ n (t = 0) for the driven case by numerically constructing the evolution operator over one period T = 2π ω ,Û (T ), and diagonalizing it. In Fig. S10 we plot the bipartite von Neumann entanglement entropies of all the Floquet modes for the optimal driving parameters. There are two symmetric "arcs" in the entropy plot, which suggests that the Floquet Hamiltonian might be fractured into two components. The expected numbers of excitations for each mode are represented by different colours. The two lowest entropy modes have the highest overlap with the Néel state |Z 2 . One of them is approximately Φ 1 (0) = (|Z 2 + |Z 2 )/ √ 2 and the other is close to Φ 2 (0) = (|Z 2 − |Z 2 )/ √ 2, while the quasienergy separation between them is ∆ = 1 − 2 ≈ ω/2. This provides a simple explanation for the revival dynamics starting from the Néel state, as will be outlined below.
Let us assume that the two idealized states Φ 1 (0) and Φ 2 (0) are indeed Floqeut modes. The initial state |Z 2 will then be a superposition of only these two modes and will evolve as After one driving period, the two Floqeut modes will return to their initial states, but the relative phase will be e i ω 2 2π ω = e iπ = −1. The wavefunction after one period will therefore be in the anti-Néel state (with an unimportant phase prefactor), ψ(T ) = e −i 1T |Z 2 . It will take two driving periods for the relative phase to again become +1 and the wavefunction to return to the initial |Z 2 state. This is the origin of the period doubling (subharmonic response to periodic driving) which was observed in previous works [S12].
We note that the period doubling will disappear if we resolve the translation symmetry and work only in the k = 0 momentum subspace. The initial state (|Z 2 + |Z 2 )/ √ 2 will in that case have high overlap with only a single Floquet mode and will trivially oscillate with the same frequency as the periodic drive. As the tilted Bose-Hubbard model can be mapped to the PXP model in the U ≈ ∆ J limit, see Section , we also expect to be able to enhance many-body scarring via periodic modulation of the term corresponding to the number of excitations. In the Bose-Hubbard case, such a term is conveniently provided by the on-site interaction strength U . However, we cannot use periodic boundary conditions due to the linear tilt which would be discontinuous at the boundary. We therefore consider the Bose-Hubbard model with open boundary conditions and periodically modulate the interaction strength U (t): with the driving given by The driving parameters, U 0 , U m and ω, are the static detuning and the modulation amplitude of the interaction strength and the driving frequency, respectively. The modulation of interaction strength indeed leads to enhanced revivals in the Bose-Hubbard model, see Fig. S11. In particular, the slope of entanglement growth is significantly reduced, with scarred oscillations becoming more pronounced. However, in local observables, such as the density of doublons, the effects of driving are less striking than in the pure PXP model. The reason for more modest enhancement of revivals in the Bose-Hubbard model is the competition between stabilization of revivals within the PXP subspace and the processes which destroy the mapping to PXP model, such as the terms creating 3 or more bosons on a site. Additionally, the optimal driving parameters are not the same as those for the PXP model (up to the trivial rescaling by Ω = √ 2J to match the normalization of off-diagonal matrix elements). Increasing the tilt parameter ∆ brings the tilted Bose-Hubbard model closer to the PXP model, but it is still necessary to perform a separate optimization of driving parameters.

QUANTUM MANY-BODY SCARS IN THE POLARIZED STATE
In the main text we reported the observation of manybody scarring associated with the state that contains no doublons, |111 . . . , or equivalently the fully-polarized state | • • • . . . in the PXP model. In this section we provide extensive theoretical evidence for many-body scarring in the polarized state. While the polarized state does not exhibit many-body scarring in the pure PXP model, consistent with previous work [S14], it does display weak signatures of non-ergodicity in local observables for sufficiently small systems. In this section, we show that static detuning and its periodic modulation can be used to stabilize the scarring from this initial state. As we will demonstrate below, the many-body scarring in the polarized state is distinct from the previously studied "dynamical freezing" regime associated with | • • • . . . state in the PXP model driven by a square pulse protocol [S15].  are no significant revivals in wave function fidelity, some oscillations in the number of excitations are still visible.
A closer look at the eigenstates of the PXP Hamiltonian and their overlap with the polarized state reveals the underlying reason for this behaviour, see Fig. S13. In Fig. S13(a) we plot the overlap of all PXP eigenstates with the Néel state, showing the well known band [S5] of scarred eigenstates marked by the red crosses and corresponding tower structures. In contrast, there is no such band of high-overlap eigenstates for a randomly chosen state, see Fig. S13(c).The polarized state is between these two cases, as can be observed in Fig. S13(b). Although there is no well defined band of scarred eigenstates as for the Néel state, there is still a number of unusually high-overlap eigenstates which are marked by the black dots. Finally, in Fig. S13(d) we show the entanglement entropies of all eigenstates. The lowest-entropy eigenstates are the Néel state scars (red crosses), but the eigenstates with the highest overlap with the polarized state (black dots) also have lower then average entanglement entropies. Thus, we conclude that the polarized state is poised to develop many-body scarring by a suitable perturbation of the PXP model. We next show that this can be achieved by applying static detuning.  In addition to having the highest overlap with the polarized state, the special states are also approximately equidistant in energy and have lower entanglement entropy than most other eigenstates. These are all paradigmatic properties of quantum many-body scars. However, one striking difference compared to the Néel state scars is that the highest-overlap states are not concentrated in the middle of the spectrum. Instead, most of them are located at one edge of the energy spectrum, but the band of atypical states still continues well into the higher energy densities, see Fig. S14(c). The fact that special eigenstates are biased towards one end of the spectrum is expected since the detuning potential breaks the particlehole symmetry of the PXP Hamiltonian [S6]. This emergence of scarred eigenstates significantly affects the revival dynamics, as illustrated in Figs. S12 and S15. For the Néel state in Fig. S15(a), the detuning monotonically destroys the revival, until we reach the regime of very large detuning µ 0 3 which places the Néel state in its own fragment of the Hilbert space. By contrast, for the polarized state in Fig. S15(b) we see the revivals start to emerge at moderate detuning µ 0 /Ω ≈ 1. The frequency of the revival is found to match the energy separation between the scarred states in Fig. S14(c). The oscillations in the number of excitations are also enhanced and their frequency has changed to the frequency of fidelity revivals. This is the regime that corresponds to the many-body scarring observed in experiment. We note that the revivals from the polarized state also persist in the trivial large-detuning limit (µ 0 3) where the polarized state is effectively in its own fragment of the Hilbert space, similar to the Néel state.
The addition of detuning not only affects the shorttime dynamics, but also infinite-time expectation values. After a quench, the value of any observable will reach the value predicted by the diagonal ensemble However we also expect the observable to thermalize towards the value predicted by the canonical ensemble O th = Tr ρ thÔ , whereρ th = 1 Z e −βĤ with Z = Tr e −βĤ and β the inverse temperature. Note that we also restrictĤ to the symmetry sector invariant under translation and spatial inversion as it is the only one compatible with the |0 state. A large difference between the predictions of these two ensembles for a given initial state is an indicator of the violation of the Eigenstate Thermalization Hypothesis [S16, S17]. For the PXP model we will use the operatorn = 1 N jn j , which counts the average number of excitations in the system [S18], and denote the difference between the ensemble predictions by δn. The Néel state is most athermal at zero detuning, while the peak for the polarised state occurs around µ 0 /Ω = 1.7, see Fig. S16. For larger values of the detuning these two states become respectively the topmost and ground states, meaning that the temperature is ±∞ and both ensembles agree exactly.

Periodic driving in the PXP model
Finally, in order to stabilize revival and many-body scarring in the polarized state at late times, we need to modulate the detuning amplitude, in addition to the static detuning. Using the same driving protocol as for the Néel state in Eq. (S23), we can enhance and stabilize the revivals from the polarized state at late timessee Fig. S12 (blue lines). The optimal driving frequency is found to be close to the frequency of revivals in the undriven case with static detuning.
In Fig. S17 we plot the entanglement energies of all the Floquet modes Φ n (t = 0) for the optimal driving parameters. As in Fig. S10, the colour scale represents the expected number of excitations for each mode. The translational symmetry is now resolved and we show only the modes inside the k = 0 momentum subspace. There is a single mode that has high overlap with the polarized state, which explains the revival dynamics in Figs. S12 and S15(b). Note that there is no period doubling in this case.
Detuning and periodic driving in the tilted Bose-Hubbard model Finally, we confirm that our conclusions about manybody scarring associated with the polarized state also hold in the full tilted Bose-Hubbard model in the regime U ≈ ∆, where we expect the effective description to be close to the PXP model. We will show that the driving leads to a strong suppression of entanglement growth and makes off-resonant the processes that cause leakage out of the PXP subspace.
In Fig. S18 we compare the dynamics at the resonance U = ∆ (black lines, corresponding to the pure PXP model), at U = ∆ + U 0 (red lines, corresponding to the PXP model with static detuning), and for U (t) = ∆ + U 0 + U m cos(ωt) (blue lines, corresponding to the periodically driven PXP model). Due to the very fast growth of the Hilbert space size, we restrict the maximal number of bosons per site to 3. The results are consistent with those for the PXP model shown in Fig. S12. Note that the frequency of fidelity revivals in Fig. S18(a) is the frequency of PXP revivals multiplied by a factor of √ 2 which comes from the off-diagonal matrix elements in the Bose-Hubbard model. The expected number of doublons, which is related to the number of PXP excitations is shown in Fig. S18(b).
The growth of entanglement entropy is suppressed by the addition of static detuning and even more by periodic driving, see Fig. S18(c). There are two factors that contribute to this behaviour. One is the dynamics inside the PXP subspace. The other is related to the leakage out of this subspace, which is represented by the number of sites with 3 particles in Fig. S18(d). The static detuning by itself significantly decreases this quantity, while the periodic driving does not seem to result in a substantial further improvement for the polarized state.

OTHER QUANTUM MANY-BODY SCARRED STATES
In addition to the Néel state and the polarized state, we also find other initial states which revive in the PXP model with static detuning, Eq. (S28). These initial states are the ground states ofĤ(µ i ) and they exhibit revivals when the detuning is quenched to a different value, H(µ i ) →Ĥ(µ f ). This setup generalizes the quench protocols studied in the main text. For example, setting µ i →−∞, the ground state is simply the Néel state and then quenching to µ f =0 (pure PXP model) leads to the appearance of Z 2 quantum many-body scars. Similarly, if we set µ i →∞ then the ground state is the |0 state and quenching to µ f =1.68Ω also leads to scarring, as shown in the main text.
However, we observe similar scarring phenomenology in a larger set of initial conditions by varying the parameters µ i and µ f . Unlike the Néel state and the polarized state, the ground states ofĤ(µ i ) are far from product states for |µ i | < 2. However, they also have low entanglement and can be prepared experimentally. In Fig. S19 we illustrate this with an example for µ i =−0.76Ω and µ f =1.6Ω. For this set of parameters, we recover a similar scarring phenomenology as shown in the main text for the polarized state with µ 0 =1.68. However, we stress that the initial state considered here, i.e., the ground state ofĤ(µ i ), is now far from both the Néel and polarized states (the overlap with these states is on the order 10 −5 ). We emphasize that these values of µ are not fine tuned, and we find large regions of µ i and µ f leading to scarring.
In Fig. S19 (a) we recover a dynamics close to what could be observed for the polarized state 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Ĥ(µ f ) is shown in Fig. S19(b). Comparing with the top states having overlap with the polarized state at the same value µ f (red crosses) we see a similar pattern. More than that, the same atypical eigenstates have a high overlap for both states, but the phase is different in the two cases. This is similar to what is observed with the Néel state and its translated version for µ f =0. Both have the same overlap magnitude with each eigenstate, but the phases are different.

EFFECT OF DETUNING ON THE SPECTRAL STATISTICS OF THE PXP MODEL
In this section we show that the addition of finite detuning to the PXP model does not make this model integrable. We study the energy level spacings s n = E n+1 − E n , which we normalize to have s n = 1. For an integrable model, {s n } should follow the Poisson distribution, while for a chaotic model we expect to see the Wigner-Dyson distribution. A convenient way to probe level statistics is by computing the so-called r parameter [S19], defined as the average of level spacing ratios: r n = min(s n , s n−1 ) max(s n , s n−1 ) .
For the Poisson statistics, we expect r ≈ 0.39, while r ≈ 0.53 for Wigner-Dyson. In Fig. S20 we show that r tends towards 0.53 as N increases, for all values of µ.
In general, as µ becomes larger, the convergence is slower because the detuning approximately conserves the number of excitations. Beyond that, one can also notice two dips in r at µ = 0 and µ ≈ 1.6, hinting that near these values PXP is close to another integrable model. For pure PXP this had been noted and previously investigated with various other perturbations [S20]. The full distribution of the s n is shown in Fig. S21 for µ = 0, 1, and 1.6829 for N =32 spins. In all cases, we see that the distribution resembles Wigner-Dyson, even though in the last case it is skewed towards zero.
In conclusion, for any finite value of µ, the PXP model is non-integrable and its level statistics follow the Figure S19. Emergence of many-body scarring by quenching the PXP model from µi = −0.76Ω to µ f = 1.6Ω for N = 32 spins in the k=0, p=1 symmetry sector. The dynamics is very similar to that of the polarized state at the same µ f , with the overlap between the two suggesting that state transfer happens between them. The overlap between these two states and the eigenstates ofĤ(µ f ) also shows similar tower structures.  Figure S20. r for the PXP model with various system sizes N and detuning µ. For all values of µ shown, the spectral statistics flows towards Wigner-Dyson value, as the r parameter increases with system size. However the convergence is slower near µ = 0, µ = 1.68, and in general as µ becomes larger.
Wigner-Dyson distribution in a large enough system size. Interestingly, the level statistics suggests a proximity to an integrable model at the points where we find good revivals due to scars: at µ = 0 for the Néel state and near µ = 1.68 for the polarized state. These results are in accordance with the discrepancies observed between the diagonal and canonical ensembles in Fig. S16.

SYSTEM-SIZE SCALING OF THE REVIVAL FIDELITY
An important question concerns the stability of revivals in the thermodynamic limit. In particular, due to the cost of non-linear optimization, the driving parameters were obtained in relatively small systems, therefore it needs to be checked whether the same parameters work as well in large systems. For the calculations in this section, we perform time evolution in large systems using the TEBD algorithm implemented in TenPy package [S9], as explained in Sec. .
Figs. S22 and S23 show the system size scaling of the first three revival peaks for different initial states, both with and without driving. The results were obtained from numerical simulations of the tilted Bose-Hubbard model, Eq. (S26), with open boundary conditions, ∆/J = 16, and maximally 3 particles per site. This particle number limit is an reasonable assumption since the periodically driven interaction strength U (t) = ∆+U 0 +U m cos(ωt) is large compared to the hopping amplitude J. In the case of global fidelity F (t), we plot the so-called fidelity density − ln(F (nT ))/L, where T is the revival period and n ∈ {1, 2, 3}. The single-site fidelity F (1) (t) is a local quantity, so it does not need to be rescaled by the system size L. We therefore simply plot the peak heights F (1) (nT ). For the Néel state 2020 . . . 20, the fidelity density is expected to converge to a constant value in the limit of large L. This is consistent with our results in Fig. S22(a), where we plot the fidelity density after one, two and three driving periods. The driving parameters are the same for all system sizes, ω = 3.85, U 0 = 1.85 and U m = 3.71. Due to the minus sign in the definition, lower fidelity density corresponds to higher revival peaks and vice versa. As can be observed in Fig. S22(a), periodic driving leads to increased revivals over a broad range of system sizes and there is no indication that this will change for L > 50. The revivals are decaying with time, but the decay is significantly slower when the driving is turned on. We can thus conclude that periodic driving with these parameters both enhances and stabilises the revivals, even in relatively large systems.
The scaling of the single-site fidelity can be observed in Fig. S22(b). This experimentally measurable quantity represents a tight upper bound for the global fidelity when the system is initialized in a product state, see Sec-tion . The results are similar to those for the global fidelity. In all cases, the revival heights are rapidly converging towards a constant value. Again, the revivals in driven systems are significantly higher than those without driving and the difference between them increases with time. The effects of periodic driving are even more striking with the polarized state 111 . . . 111 as the initial state, as shown in Fig. S23. There are no notable revivals in global fidelity when the driving is turned off. The dashed lines in Fig. S23(a) correspond to irregular minor local maxima which are present in smaller systems. Even these local maxima disappear with increasing system size, which explains why some data points are missing. In contrast, driving with parameters ω = 4.90, U 0 = −2.38 and U m = 1.54 produces very high revivals which do not decay, either with time or with system size. The singlesite fidelity tells a similar story, see Fig. S23(b), however in this case there are revivals even in the absence of driving, consistent with dynamics of local observables in Fig. S12(b).
Finally, we note that the Néel and polarized states are the only two initial product states for which we were able to find optimal driving parameters that lead to robust revivals at late times. This is true both for the tilted Bose-Hubbard model, Eq. (S26) in the ∆ ≈ U regime, and for the PXP model with a spatially uniform driving protocol. For other initial states, such as Z 4 state with an excitation on every fourth site or, equivalently, 20112011 . . . 2011 in the tilted Bose-Hubbard model, it is possible to stabilize a small number of revivals at short times. In contrast to the Néel and polarized states, these revivals are found to decay quickly with time as well as with system size.