The fate of quantum many-body scars in the presence of disorder

Experiments performed on strongly interacting Rydberg atoms have revealed surprising persistent oscillations of local observables. These oscillations have been attributed to a special set of non-ergodic states, referred to as quantum many-body scars. Although a significant amount of research has been invested to understand these special states, it has remained unclear how stable scar states are against disorder. We address this question by studying numerically and analytically the magnetization and spatio-temporal correlators of a model of interacting Rydberg atoms in the presence of disorder. While the oscillation amplitudes of these observables decay with time as the disorder strength is increased, their oscillation frequency remains remarkably constant. We show that this stability stems from resonances in the disordered spectrum that are approximately centered at the same scar energies of the clean system. We also find that multiple additional sets of scar resonances become accessible due to the presence of disorder and further enhance the oscillation amplitudes. Our results show the robustness of non-ergodic dynamics in scar systems, and opens the door to understanding the behavior of experimentally realistic systems.


I. INTRODUCTION
Recent progress in the design and coherent control of programmable quantum simulators has made it possible to discover new and striking non-equilibrium phenomena [1][2][3][4]. A particularly remarkable example has been the detection of quantum many-body scars using a 51 qubit cold-atom platform [4]. In this system, strongly interacting Rydberg atoms were initially prepared in a product state and subsequently allowed to evolve under their own unitary dynamics. Unexpectedly, the dynamics of measured observables presented persistent oscillations at finite-energy density. This represented a violation of the Eigenstate Thermalization Hypothesis (ETH) [5], a central tenet of statistical mechanics. The ETH says that states at finite energy density evolve at long times into ergodic states wherein the expectation value of local observables converge to thermal values. The experimental observation of persistent oscillations thus signals that non-ergodic behavior is at play and that physics beyond the ETH was accessed by the Rydberg quantum simulator.
Subsequent analyses revealed that ETH was violated due to the presence of special states located throughout the energy spectrum which were termed many-body quantum scars [6,7]. These states are embedded in an ergodic spectrum and are approximately equally spaced in energy, separated by the characteristic energy scale ω scar . They exhibit less-than volume law entanglement and have high overlap with the initial product state used in the experimental simulation [4]. The observed persistent oscillations at the frequency ω scar thus follow from these properties of quantum many-body scars. These states have garnered intense interest and have been found in a wide range of systems, both static and driven, as well as in higher dimensions [8][9][10][11][12][13][14][15][16][17][18]. Several insights have been obtained by constructing exact scar eigenstates and from analytical results [19][20][21][22][23][24]. Strategies have also been devised to understand and enhance the quality of persistent dynamics by recognizing emergent algebras in scar subspaces [25][26][27]. Furthermore, fundamental questions about the nature of such scar states have been addressed by studying their proximity to integrability, as well as connections with concepts of quantum chaos [28][29][30][31].
An important question that has not yet been addressed is the effect of disorder on the stability of many-body scars. While disorder was discussed in [32,33], it was only in the context of fine-tuned models where disorder was designed not to affect the scar states themselves. The effects of disorder in the PXP model were studied in [34], but the focus was solely on many-body localization physics and not on the non-ergodic dynamics arising from quantum scars. Thus, the broader question of the effects of generic disorder on quantum scars has remained open. One might expect that the randomizing effect of a disordered environment will generically suppress the oscillation amplitude and distort the oscillation frequency of observables. On the other hand, one could also expect that at least some quantum scars will be robust against disorder given that strong kinematic constraints underlie their existence in some models. These constraints introduce strong spatio-temporal correlations, suggesting enhanced robustness with respect to spatial imperfections.
Undertaking this study is important for both practical and fundamental reasons. As a practical matter, near-term quantum simulators, such as those based on superconducting qubits, are naturally affected by qubitto-qubit variations [35][36][37]. If quantum scars are to be probed in such platforms, it is necessary to understand how disorder will affect their detection. At the same time, disorder has played a central role in questions of ergodicity and non-equilibrium dynamics in many-body quantum systems. The fundamental phenomenon of many-body localization (MBL) was predicted for sufficiently strong disorder in a general class of equilibrium interacting electronic systems [38]. Further studies of non-equilibrium phenomena in the presence of disorder have shown that many-body localization and associated phases of matter persist in a wide variety of interacting system even far from equilibrium [39][40][41][42][43][44][45][46][47][48]. In the space of models, systems with scars reside somewhere between the integrable and fully ergodic systems. It is therefore natural to expect their response to disorder to yield a spectrum of unusual features, both in and out of equilibrium. This is indeed what we find in this work.
In this work we study how scars respond to generic disorder. We use the Rydberg Hamiltonian which accurately models the original quantum simulator where quantum scars were discovered [4] (the model system is schematically represented in Fig.1a). We find that eigenstate diagnostics of scar states [6] fail when disorder is turned on. However, clear oscillations can still be identified in observables such as the magnetization of the system as well as spatio-temporal correlators. While disorder indeed leads to the decay of these oscillations with time, their frequency remains close to ω scar for a finite range of disorder strengths, exhibiting significant rigidity to generic perturbations.
We provide a physical picture for these observations by deriving analytically the disorder-averaged dynamics of the magnetization and temporal correlators in the weakdisorder limit. This derivation shows that, instead of remaining discrete eigenstates, scar states become transformed into resonances in the energy spectrum of the disordered system. For sufficiently small disorder strength, these resonances have widths smaller than ω scar , and are centered at the energies of the clean scars states. As a result, the observables we probe continue to oscillate at the scar frequency before decaying irreversibly. Surprisingly, we also uncover multiple additional scar resonances not associated with the previously known scar states; these novel states come into play in the presence of disorder, and, remarkably, enhance the magnetization oscillations at intermediate times. While in the bulk of this work we focus on the regime of strongly interacting Rydberg atoms (the so called "PXP" regime), at the end we fully map out the dynamical regimes of the Rydberg model for arbitrary coupling strength as a function of disorder.
Our results show that quantum many-body scars can be detected in moderately disordered quantum simulators. This in itself can potentially prove a valuable tool, for example, in the calibration of quantum hardware. Furthermore, we find that disorder is able to reveal the multi-tower structure of scar resonances in the Rydberg simulator. This insight serves as a new window into the underlying complexities of systems with scar states that could help further understand the circumstances under which they arise and how they can lead to ETH violation. (1). The excited and ground states of each qubit are separated by an energy Ω. Only when two neighboring qubits are in their excited states will the qubits interact, with energy JR. When JR Ω, a kinematic constraint develops that leads to the formation of quantum many-body scars. (b) Diagram of the dynamical regimes obtained in this work for a chain of Rydberg atoms as a function of the nearestneighbor interaction JR and disorder strengths W : a regime that exhibits scar dynamics due to scar resonances, a regime that is ergodic without scar resonances, a constrained MBL phase, and an MBL without kinematic constraints.

Outline
We now outline the structure of this work and summarize the main results. We begin in Section II by introducing the model we use as a case study, namely the Rydberg Hamiltonian which describes interacting Rydberg-atom qubits; the system is shown schematically in Fig.1a. We discuss basic aspects of this model, including how scars arise in the strongly interacting limit, and also the type of disorder we consider in this work.
In Section III, we consider the effects of disorder on the eigenstate properties of the PXP model. We discuss how scar states are no longer detectable in the energy spectrum ( Fig.2) and determine how the ergodicity of the system is affected (Fig.3). We then move on in Section IV to carefully analyze the magnetization of the system. We show that it continues to carry signatures of nonergodic oscillations at the scar frequency for moderate disorder strengths due to the kinematic constraint of the PXP model (Figs.5,6,7).
In order to understand these results, in Section V we derive explicitly the disorder-averaged magnetization dynamics in the weak disorder limit. This derivation shows that the dynamics is controlled by disorder-induced resonances that arise from multiple towers of scar states of the clean system. Although their widths grow with disorder strength and this leads to the decay of the oscillation amplitude, non-ergodic oscillations continue to manifest in the magnetization as long as the width remains smaller than the scar frequency (Fig.9). We further show in Section VI that scar resonances also account for temporal correlations in the system (Figs.11,12).
We bring everything together in Section VII by showing a diagram in the parameter space of interaction and disorder strengths which shows four distinct dynamical regimes (Fig.1a): a regime with scar resonances, a regime that is ergodic and does not exhibit revival dynamics, a constrained many-body localized phase, and another many-body localized phase without kinematic constraints. Finally, in Section VIII we discuss the main conclusions from this work, we comment on how our results can be tested in various quantum simulators, and propose future directions of research.

II. MODEL
In this section, we introduce the model that will serve as a case study to understand how quantum many-body scars respond to disorder. We begin with the standard Rydberg Hamiltonian given by [6,7,11,49] where the Pauli matrices σ x,y,z r act on ground and excited qubit states {| , | } respectively, we defined the projection operator P ± r = 1 2 (I ± σ z r ) , and we assumed periodic boundary conditions. This Hamiltonian models a chain of Rydberg atoms, each driven resonantly with strength Ω (Rabi frequency) and interacting with strength J R when two nearest neighbors are in the excited state [50,51].
In the strongly interacting limit J R Ω, there is a large energy cost J R to go from the state | · · · · · · to the state | · · · · · · , which has two neighboring qubits in their excited state. As a result, when J R Ω, the energy spectrum splits into L bands. The Hilbert space spanned by each of these bands is characterized by a fixed number n exc of pairs of contiguous qubits in their excited state. When the system is initialized within one of these bands, the number n exc is thus preserved by the dynamics.
In particular, the product state |Z 2 = | . . . , which we will use as the initial condition in this work, satisfies n exc = 0. The effective Hamiltonian that generates the dynamics within this subspace, referred to as the PXP model [6], is given by where P = r I − P + r P + r+1 projects into the n exc = 0 sector. The resulting Hamiltonian differs from that of a paramagnet by the presence of projection operators. To understand the impact on the dynamics of these projectors, we can express the Hamiltonian in the equiva- , which holds if we consider the dynamics of initial states that satisfy P|Ψ(0) = |Ψ(0) . This form of the Hamiltonian illustrates that, while in a paramagnetic system the operator σ x r generates rotations of the qubit at the position r without regard to the state of other qubits, in the strongly interacting case this rotation can only happen when the qubits located at the positions (r − 1) and (r + 1) are in their respective ground states.
One striking consequence of this kinematic constraint is the large (although incomplete) and periodic revival dynamics of product states such as the |Z 2 state. Since this is a non-integrable model due to the kinematic constraint, the Eigenstate Thermalization Hypothesis (ETH) dictates that an initial state at finite-energy density should evolve under Eq. (1) into a thermal state at long times. However, the persistent and periodic revival of the Z 2 state constitutes a direct violation of the ETH. Spectrally, this can be understood in terms of a set of L + 1 special eigenstates of the PXP model, referred to as quantum scars, which are embedded in a presumed ergodic spectrum within the n exc = 0 subspace. They are approximately equally spaced by an energy ω scar = 2πν scar = η Ω, η ≈ 0.636.
Scar stats have an exceptionally large overlap with the Z 2 state, and exhibit less-than volume law spatial entanglement even though they reside at finite energy density. The properties of these special eigenstates explain the non-ergodic revival dynamics of the |Z 2 state. Now, the return probability of the Z 2 state obtained by evolving with the PXP model is only partial. Furthermore, its spatial entanglement entropy grows with time, which suggests that the quantum scars of the conventional PXP model exhibit some level of hybridization with ergodic states. To correct for this, a suitable deformation δH can be added to the Hamiltonian that removes this growth in entanglement entropy and further optimizes the return probability of an initial Z 2 to reach unity [25]. One can think of this ideal limit of optimized scars as a reference point for understanding the impact of disorder on its non-ergodic dynamics. We will use the particular deformation where δJ R = −0.026Ω [25]. Longer-ranged terms further improve the quality of scar states, but for our present purposes this deformation is sufficient. The main subject of our work is the effect of disorder on the scar states, the corresponding revivals, and order parameter oscillations. The disorder that we introduce is generic in that it breaks all the spin and translational symmetries of the clean PXP model; namely, we consider the operatorŴ = r,a h a (r)σ a r .
The random fields h a (r) with a = x, y, z are chosen uniformly distributed in the range [−W/2, W/2], where W measures the strength of the disorder. The presence of all Pauli matrices breaks the symmetries that the clean PXP model has, such as chiral and time-reversal symmetries. Throughout most of this work we will mainly study the strongly interacting limit. This effectively projects the disorder operator into the n exc = 0 Hilbert space. As a result, the full Hamiltonian is given by Later, in Section VII, we will lift this strong constraint to study the behavior of the system as a function of interaction strength J R in the full Hilbert space. When disorder is introduced into the system, the fields in Eq. (5) will tend to make the qubits in the system precess around random directions and at different rates. As a result, there will be a competition between the desynchronizing effect introduced by the disorder potential and the clean correlated oscillatory dynamics of the PXP Hamiltonian. Understanding the outcome of this competition is the main goal of this work.

III. SPECTRAL PROPERTIES
In this section, we study the impact of disorder on spectral and eingenstate properties of the Hamiltonian (6) that are commonly attributed to the existence of nonergodic dynamics of scar systems. We present evidence that scars can no longer be clearly identified as eigenstates of the system for moderate disorder strengths. Furthermore, we discuss how disorder initially enhances the ergodicity of the system at low disorder, and eventually induces a many-body localization transition at strong disorder.

A. Destruction of scar eigenstates
A basic question that first arises when disorder is introduced is whether scars can still be identified as eigenstates in the spectrum of the system. One standard way to identify scar states is to single out those eigenstates that have a high overlap the Z 2 state [6,7]. In Fig.2a, we show the distribution of this overlap for W = 0 and W = 0.25ω scar using a fixed disorder realization. While in the clean limit it is easy to identify states that have large overlap with the Z 2 state, it is difficult to do so when disorder is turned on. Disorder has hybridized the cleanlimit scar states with other spectrally near-by states to such an extent that the Z 2 state no longer appears to be concentrated in a small subset of states. As a result, there are no clear signs of scar eigenstates in this overlap diagnostic.
Further evidence that scar eigenstates are absent in the spectrum can be obtained by calculating the spatial entanglement entropy of each energy eigenstate. If we divide the system in two halves A and B, the entanglement entropy of a state |E n is S(E n ) = −Tr ρ (n) where ρ (n) A = Tr B {|E n E n |} and Tr B traces the degrees of freedom in region B. When the entanglement entropy is shown for each eigenstate in the system in the clean limit, as seen in Fig.2b, it is possible to clearly discern isolated scar eigenstates that appear to have less-than volume law entanglement. Upon introducing disorder, however, the states of low entanglement recede quickly into the full set of volume-law ergodic states. Thus, we again see that hybridization of optimized scars with ergodic states has erased signatures of scar states in the entanglement of energy eigenstates.

B. Enhanced ergodicity and localization transition
In addition to obscuring signatures of scar eigenstates as we have seen in the previous section, disorder can also lead to other kinds of non-ergodic behavior such as manybody localization [40,43]. It is thus pertinent that we examine how the ergodicity of the system changes as the disorder strength is varied. To this end, we calculate the mean value of the level spacing ratios [39] where δE n = E n+1 − E n with the energy eigenvalues {E n } sorted in ascending order. Since the disorder operator we have chosen does not respect any symmetries, the disordered Hamiltonian we are studying belongs to the Grand Unitary Ensemble (GUE), implying that [ r n ] ≈ 0.6 when the system is ergodic [52], where · denotes spectral averaging. When the system is localized, such as in an MBL phase, Poisson statistics dictates that [ r n ] ≈ 0.38. We present [ r n ] in Fig.3a as a function of system size. The flow with system size reveals signatures of a localization transition at around W Th-L ≈ 2.13ω scar . When W < W Th-L , the approach to the transition point does not occur monotonically for finite sizes. Instead, the system presents an enhanced approach to the ergodic value [ r n ] = 0.6 at a disorder strength somewhat be-low W Th-L . When W > W Th-L , notice that [ r n ] converges quite slowly to the Poisson value at strong disorder, which is the value that is expected for a manybody localized phase. This is clearly due to the kinematic constraint imposed by the presence of projection operators. Slow convergence to the Poisson value was also obtained in [34] for a similar disordered PXP model, where this type of many-body localization under kinematic constraints was dubbed constrained MBL.
Since level statistics at low disorder can be affected by the symmetries of the clean system, we also compute an additional diagnostic of ergodicity that makes use of the eigenstates of the system to confirm these results. The ETH states that the expectation value of local ob-servablesÔ varies smoothly with energy. As a result, contiguous states satisfy . This operator was used to study signatures of integrability in the clean PXP model [29]. Near infinite temperature, we have S(E) ∼ ln D, where D is the dimension of the Hilbert space. This implies that ∆Ô ∼ D −1/2 if the system is ergodic. We proceed as in [29], and we calculate ∆Ô ∼ D α as a function of D to extract the exponent α. For ergodic systems, the exponent α should thus converge to −1/2, whereas in the non-ergodic regime we expect ∆Ô to decay slower with |α| < 1/2. We show the exponent α extracted by using L = 12, 14, 16 in Fig.3b. The exponent approaches −1/2 below W Th-L , and is significantly suppressed above it, consistent with our findings using energy level statistics.

IV. SYSTEM DYNAMICS
The results from the previous section raise the question of whether non-ergodic dynamics can still be detected in the presence of disorder. To explore this question, we turn to the study of the magnetization of the system. We define the on-site magnetization as where |Ψ(t) is the time-evolved state obtained from solving the Schrodinger equation with the initial condition |Ψ(0) = |Z 2 . We further define the Fourier In what follows, we will first calculate the magnetization in the clean limit in order to understand how it depends on the presence of scar eigenstates. We will then investigate its behavior when disorder is introduced. spectrally ordered spectrally ordered ordered by towers ordered by towers . Although the behavior of these matrix elements is not known for arbitrary wave vector k, progress can be made for the particular case when k = π and the initial state is |Ψ(0) = |Z 2 [25,49]. Given that the Z 2 state is approximately spanned by optimized scar states, we only need to consider matrix elements φ l∈Λs , where Λ s is the set of indices of optimized scar states. In particular, approximate ladder operators σ ± s can be defined within the scar subspace where which connect scar states that differ by an energy ±ω scar (see Appendix A for details). The matrix elements φ l∈Λs that contribute to the magnetization are those for which ω l l = ±ω scar . Plugging this into Eq. (11) with k = π leads to the magnetization components which exhibit oscillations at the scar frequency with fixed amplitude. These magnetization oscillations capture the non-ergodic dynamics induced by the presence of scars in the spectrum, and could thus be useful observables to track any remaining signatures of scar physics when the system becomes disordered. Now, when disorder is introduced, other states of the clean system outside of the scar subspace will inevitably become populated. As a result, in the following sections we will need to also understand the behavior of matrix elements φ l ∈Λs . Not much has been said in the literature, however, about their behavior. These states have been largely assumed to be ergodic, which might suggest that they do not contribute to oscillations at the scar frequency. Consider, however, Fig.4a, where we show a density plot of all the matrix elements | φ at the same disorder strengths. For weak disorder, oscillations occur at (k, ω) = (π, ωscar). As the disorder strength is increased, the distribution broadens, and eventually concentrates at (k, ω) = (π, 0).
performs an appropriate re-ordering of the energy basis by grouping them into sets that are connected via the same ladder operators σ ± s . This re-ordered basis reveals that the energy eigenstates can be organized into towers labeled by an index J, with each tower containing a number D J of states labeled by an index m. There are a few states near the middle of the spectrum that do not seem to form towers with any other state, so they form sets with a single state each; such sets do not contribute appreciably to the magnetization.
Using this reordering, we again show | φ Fig.4b (details of how we performed this re-ordering are presented in Appendix B). The matrix is approximately block-diagonal with respect to J, and each block is tri-diagonal, so we approximately write for a = y, z. In this expression, we defined Γ a,κ Jm = φ Jm . Similar to the case of optimized scars (which in our notation corresponds to the J = 1 tower), this tri-diagonal structure selects frequencies in Eq. (11) that are close to ω scar for most towers. As we will see in the following section, upon introducing disorder, these additional towers will become populated and their nonergodic oscillations near the scar frequency will be revealed.
We note that it is remarkable that the PXP model is comprised of multiple towers of scar states, of which the set of optimized scar states is just one example. One of the reasons that such towers have remained hidden is, in part, because the Z 2 state is spanned largely by only the J = 1 tower. The multi-tower structure of the PXP model can lead to non-ergodic SU (2) dynamics with a variety of product states as we discuss elsewhere [53]. A pronounced broadening occurs in the paramagnetic case, whereas in the PXP model the dynamics is consistently centered and localized at the scar frequency.
distribution [|M z (k, ω)|] for a few values of the disorder strength, as shown in Fig.5. Here, and throughout this work, all disorder-averaged quantities shown are averaged over 500 disorder realizations and are denoted by square brackets [·]. As can be seen in Fig.5, for the weak disorder case W = 0.1ω scar , there is a clear periodic pattern in both space and time. Temporally, each qubit oscillates as a function of time at the scar frequency, initially between the maximum amplitudes [M z (r, t)] = ±1 and gradually decaying with time. Spatially, the even and odd sites appear anti-ferromagnetically ordered throughout the dynamics, as expected from the kinematic constraint of the PXP model. This leads to sharp peaks in [|M z (k, ω)|] at (k, ω) = (π, ±ω scar ) as shown in the right column of Fig.5, consistent with Eq. (13). As the disorder strength is increased, the magnetization amplitude decays faster in time, until there is no longer any clear oscillation at the scar frequency for sufficiently strong disorder. Correspondingly, the Fourier peaks broaden, especially in the frequency domain. It is clear, however, from both [M z (r, t)] and [|M z (k, ω)|], that the frequency of oscillation continues to be dominated by ω = ω scar for a finite range of disorder strengths even though the qubits are experiencing random and biased fields. This indicates that the dynamics of the system continues to exhibit signatures of many-body scars. For strong enough disorder, the distribution eventually focuses at (k, ω) = (π, 0), indicating that the dynamics of the system has been rendered temporally trivial while maintaining antiferromagnetic ordering.
In order to better understand how stable the oscillations remain at the scar frequency, it is useful to compare [|M z (π, ω)|] for both the deformed PXP and the paramagnetic case J R = 0. The paramagnetic limit is equiva-lent to removing the projection operators from the PXP model. Because of this, the contrast between the two systems can shed light on the possible stabilizing effects that the kinematic constraint has on the oscillations of the scar system. In Fig.6a,b we show [|M z (π, ω)|] for both limits as a function of disorder strength. A sharp difference is revealed in their response to the presence of disorder. In the paramagnetic case, the dominant frequency of oscillation at the lowest disorder strength W = 0.1 ω scar we show in this figure is centered at ω = Ω, although it is already visibly broadened at this disorder strength. As W is increased, the broadening grows continuously, to the point that there is no clear dominant frequency of oscillation. As expected, the ease with which the dynamics develops a broad range of frequencies is the manifestation of each qubit oscillating at its own local disorder-induced Zeeman field.
By contrast, in the strongly interacting case there is a clear persistent dominant frequency of oscillation centered around ω = ω scar . This occurs in spite of the qubits experiencing the same random fields as those used in the paramagnetic case. The strong kinematic constraint that correlates the rotation of a given qubit with respect to its neighbors continues to robustly enforce oscillations at the scar frequency. As a check of the stability of this observation, we find that the oscillations at the scar frequency persist when the system size is increased, as we exemplify in Fig.7a, where [|M z (π, ω)|]/L does not change appreciably for L = 12, 14, 16 with W = 0.1ω scar . We thus find clear evidence that the interacting system continues to sustain oscillations at the scar frequency for a significant range of disorder strengths, even when conventional eigenstate diagnostics fail.
Oscillations at the scar frequency terminate at a characteristic disorder strength W c . This can be seen in Fig.6a, where there appears to be a marked change in dynamical regimes wherein a strong signal in [| M z (π, ω)|] eventually yields to a clear peak at zero frequency when the disorder strength is increased. To better visualize this, in Fig.7b we show [| M z (π, ω scar )|] and [| M z (π, 0)|] as a function of disorder strength. At weak but non-zero disorder, the dynamics is dominated by [| M z (π, ω scar )|]. As the disorder strength is increased, the peak at ω scar gradually decreases, while at the same time the zero frequency component increases, eventually rendering the time evolution trivial at strong disorder. The change between both regimes can be taken to occur when both frequency components have the same magnitude, namely round W c ≈ 0.63ω scar . As can be seen in this figure, the crossing between both curves does not appear to change as the system size is increased.
To capture this transition more clearly, we make use of the normalized frequency distribution [47,54]  as a function of disorder strength, for three system sizes. There is a clear and stable transition between a regime with scar dynamics and a dynamically trivial regime. In this figure, we chose q = 6 for convenience as it accentuates the position of the peak, but other values of q also work.
to study the spatial localization of wave functions.
As we show in Fig.7b, Λ (q=6) ω is suppressed at both weak and strong disorder, as expected since the system is strongly dominated by ω = ω scar and ω = 0, respectively. At intermediate disorder strengths, Λ (q) ω develops a peak at W c that does not appear to shift significantly with system size. This peak is present for other values of q = 6 as well, so our particular choice for q is somewhat arbitrary. It is used in this work because it yields a sharper peak than those obtained with lower values of q (higher values can make it even sharper). The peak signals a spread of the frequency components at W c , allowing us to identify two distinct dynamical regimes, one with clear oscillations at the scar frequency and another that is manifestly featureless.
These results thus show that oscillations at the scar frequency persist over an appreciable time scale even as dis-order is making the system more ergodic. It is clear that the increased ergodicity induced by disorder is consistent with the temporal decay of the magnetization amplitude. What is nontrivial is that these oscillations remain close to ω scar , indicating temporal rigidity.
Given that the characteristic disorder strengths W c and W Th-L appear to be distinct (as indicated in Fig.3), there are three distinct dynamical regimes that arise as a function of disorder strength: a regime with oscillations at the scar frequency, a fully thermal phase, and a constrained MBL phase. In the following section, we will present a quantitative interpretation for these observations that will show that oscillations at the scar frequency occur due to the presence of scar resonances in the spectrum.

V. SCAR RESONANCES
The behavior of the magnetization poses a number of questions regarding the presence of scars in the disordered system. Scars no longer exist as eigenstates in the spectrum, and yet the magnetization shows multiple oscillations in an increasingly ergodic system before decaying at long times. Central to the underlying dynamics is the fact that disorder couples optimized scars with the rest of the spectrum. This brings into focus, on the one hand, the manner in which scars hybridize with nearby ergodic states and, on the other, the role that states outside of the scar subspace play in the decay of the magnetization signal.
As we will see in this section, the robustness in the oscillations occurs in fact via two mechanisms. First, we will present evidence that scars continue to exist in the system in the form of resonances with a nonzero width. As time goes by, the finite life time of these resonances leads to their decay into other states in the spectrum. Additionally, although the states outside the scar subspace are more highly entangled and ergodic, they also lead to oscillations close to the same scar frequency, contributing appreciably to the magnetization of the system.
To begin our description of the problem, suppose we initialize the system in the Z 2 state which can be accurately expressed as |Ψ(0) = l∈Λs ψ l |φ (0) l , where Λ s is the set of indices of optimized scar states. We express the time-evolved state in the basis of clean eigenstates in the form , where the sum runs over the full set of energy eigenstates of the clean system since disorder will inevitably induce transitions to states outside of the scar subspace. To satisfy the initial conditions, we must have that A l∈Λs (0) = ψ l , and A l ∈Λs (0) = 0. The disorder-averaged magnetization is then The magnetization is thus determined by the dynamics of [A * l (t)A l (t)] as well as the structure of the matrix ele-   . Any deviations from the clean limit will arise from the factor [A * l (t)A l (t)] in Eq. (14). This factor can change the oscillation amplitude as well as the oscillation frequency of the magnetization [ M a (π, t)].

A. Decay of optimized scar amplitudes
Insight into the behavior of the amplitudes A l∈Λs (t) can be gained by calculating their dynamics in the weakdisorder limit W ω scar . When the disorder strength is weak, scar eigenstates primarily hybridize with nearby states within a spectral neighbourhood of size ∆ ω scar , as we schematically represent in Fig.8a. Since the number of states in this spectral neighbourhood is exponentially large in the system size, the coupling of the scar state to its spectral surrounding leads to the irreversible decay of the amplitude A l∈Λs (t). The dynamics corresponds to the decay of a discrete state that is embedded in a quasi-continuum. Standard perturbative calculations of this type of decay [55] lead us to derive the approximate expression for the amplitude (see Appendix C for details): where we defined The expression Eq. (15) represents the dynamics induced by a general disorder operatorŴ , although below we will use the explicit form defined in Eq. (5). The factor e −iδE l (t)t has fixed magnitude and encodes perturbative shifts in the scar energy due to the disorder potential. By contrast, the factor e − 1 2 λ l (t)t leads to the decay of the scar amplitude. In fact, at sufficiently long times, the factor (1 − cos ω ll t)/(ω 2 ll t) is highly peaked at low frequencies and converges to ∼ πδ(ω ll ). The sum over states in Eq. (16b) can then be turned into an integral that effectively yields the decay rate of the scar state as determined by Fermi's Golden rule. We will not take this limit here since we are dealing with finite sized systems with discrete spectra and, as we shall see, the expressions in terms of discrete sums account very well for the observed dynamics. The next step is to calculate the average [A * l ∈Λs (t)A l∈Λs (t)] over realizations of the random fields h a (r) in Eq. (5). To perform this average explicitly, we make use of a Gaussian probability distribution with the same first and second moments, [h a (r)] = 0 and [h 2 a (r)] = W 2 /12 respectively, of the box disorder distribution we used numerically. We are thus led to the explicit disorder-averaged expression , h n = h a(n) (r(n)) with n = 1, . . . , 3L, a(n) = (n − 1)mod3, r(n) = n−1 3 , we made the identification σ 1,2,3 r = σ x,y,z r , and G l (t 1 , t 2 ) = l =l,l ∈Λs |ω ll |<∆ ρ * ll · ρ T ll f mn (t 1 , t 2 ), (19) where f mn (t, t ) = iω −2 mn (ω mn (t − t ) − i(exp (iω mn t ) − exp (iω mn t))). We defined this function in terms of two temporal arguments for convenience, as it will arise again in this form in the following section. For the same reason, although the sum in Eq. (19) has the redundant constraints l = l and l ∈ Λ s , it will be useful later when we discuss the amplitudes A l ∈Λs (t) The expression Eq. (17) shows us that the magnetization decays on account of both the square root and exponential factors, encoding two sources of decay. The exponential factor represents dephasing induced by random energy shifts from Eq. (16a) which are linear in the disorder strength. The factor det{α l l (t)} arises in part from the exponential factor in Eq.15, which encodes information about the finite life time of the scar state as it decays into its spectral neighbourhood.
On the other hand, the magnetization oscillation frequency can also change if A * l ∈Λs (t)A l∈Λs (t) develops a time-dependent imaginary component. Inspection of Eq. (17) reveals that this can happen via the secondorder energy shift in Eq. (16a). Since this shift is suppressed with respect to the first order term in the energy shift, we expect deviations from ω scar to be small. This explains why, for the disorder strengths in Fig.6a, which are small with respect to ω scar , there appears to be some level of rigidity in the frequency of oscillation of the magnetization.
We now use Eq. (17) to compute the contribution of the scar states to the magnetization using the restricted sum where we defined Q Jm (t) = [A * Jm+1 (t)A Jm (t)], and Ω Jm = ω Jm+1,Jm . In Fig.9a we present [ M a (π, t)]| Λs for the numerically-obtained values and our explicit disorder-averaged result. In calculating the analytic result, the matrix elements ρ ll were calculated numerically using the eigenstates of the clean system. Both match very closely, which confirms that our basic picture of the decay of scar amplitudes is accurate. When compared, however, with the full magnetization signal, this contribution is not sufficient to account for the full magnetization of the system, even though it clearly is sufficient in the clean limit. The inevitable conclusion is that states outside of the scar subspace are contributing to the dynamics of the system. What is perhaps most remarkable is that this additional contribution must oscillate at the same frequency ω scar , even though it involves states that are known to be more ergodic and more highly entangled than the optimized scar states. This requires us to investigate how such states can coherently contribute to the observed magnetization, as we do so in the following section. The key to identifying the missing contribution to the magnetization can be gleaned from Eq. (14). It is clear that as the scar states become de-populated, amplitudes for states outside the scar subspace will inevitably begin to grow. In order to obtain the dynamics arising from the other towers, the decay of the amplitudes A l ∈Λs (t) can be calculated similar to the case of optimized scars, as we schematically represent in Fig.8b. Following the same main steps we took for the set of optimized scars, the solution for A l ∈Λs (t) can be written in the integral form (see Appendix C for details) where we defined .
Here, n 0 (l) refers to the optimized scar state that is spectrally closest to the l-th state, as illustrated in Fig.8b. These functions are straightforward generalizations of δE l (t) and λ l (t) in Eq. (16). To get an intuitive sense of how the solution Eq. (22) behaves, we can momentarily take δE l (t, τ ) → δE l and λ l (t, τ ) → λ l to be constant since they are expected to be slowly varying. We similarly do this for δE n0(l) (t) → δE n0(l) and λ n0(l) (t) → λ n0(l) . We can then perform the time integral in Eq. (22) explicitly to obtain This expression illustrates the overall behavior of the amplitudes A l ∈Λs (t). Its absolute value vanishes at t = 0 as well as at long times, reaching a maximum at intermediate times. The time scales in this evolution are the decay rate λ n0(l) of the optimized scar states that play the role of sources for the J = 1 towers; and the decay rate λ l of the J = 1 towers themselves which transition further to their spectral neighbourhood. The balance between these two time scales determines the maximum value reached by |A l ∈Λs (t)|. Note from Eq (22) that this maximum value is also controlled by the random matrix element φ n0(l) . The final step we need to take is to perform the disorder average A * l ∈Λs (t)A l ∈Λs (t) . Since it is somewhat more complicated and not very illuminating, we leave its explicit expression in the Appendix C. Using this expression, we can then calculate the contribution to the magnetization from the J = 1 towers as To make it easier to visualize the shape of these distributions, we binned the energy axis and summed the probabilities inside each bin. As can be seen, for sufficiently weak disorder, the scar states are broadened in the disordered setting but remain distinguishable. Eventually, at a disorder strength Wc ≈ 0.63ωscar (indicated here by the green frame), the width of these states increases so much that they can no longer be spectrally resolved.
In Fig.9b, we compare our explicit expression (red dots) with the numerically-obtained contribution from the lower towers (yellow solid line), finding good agreement between them. The close match confirms our assumptions about the existence of lower towers as well as the manner in which they become populated. In Fig.9c we add the analytic result when summed over all towers and again find good agreement with the numerical values for the full magnetization. Although it is more subtle in this case to infer that the J = 1 towers do not change the oscillation frequency, in Fig.9c we see that the magnetization from lower towers continues to oscillate at a frequency close to the scar frequency.

C. Distribution and stability of scar resonances
Our description of the disordered system suggests thinking about the PXP spectrum as being comprised of multiple towers of scars states that exhibit varying levels of ergodicity and entanglement. Upon disordering the system, these approximate scar states acquire a finite life time, making them effectively resonances in the disordered spectrum. To illustrate the formation of scar resonances in the system, we can calculate the distribution of the clean scar states in the basis of disordered eigenstates. In Fig.10, we show the disorder-averaged and coarse-grained probability [| φ l |φ (0) l∈Λs | 2 ] as a function of the energy E l of the disordered eigenstate |φ l . To make it easier to visualize the overall shape of these distributions, we binned | φ l |φ (0) l∈Λs | 2 in discrete ranges of energy and summed them over each bin. For weak disorder relative to W c , clear spectrally-localized distributions are obtained that are centered at each of the scar energies of the clean limit. As the disorder strength is increased, their widths increase, which is what leads to the decay of the magnetization. Since their width is less than ω scar , the distribution corresponding to each scar state is clearly distinguished from each other. As long as the width of these resonances is small compared to the scar frequency, oscillations at the scar frequency will continue to be resolved in the dynamics. Importantly, the distributions remain centered at the scar energies, which is consistent with the rigidity in the oscillation frequency we found in Fig.6a, which we also concluded from our analytic result Eq. (17). When the disorder strength becomes sufficiently strong, the distributions overlap significantly (as illustrated by the green box in Fig.10), and thus we expect oscillations at the scar frequency to disappear.
This picture of resonances allows us to infer an approximate value for W c . Given that the contribution from J = 1 towers is smaller and appears at finite times, we focus on the main scar tower to determine this value. The width of these distributions in the frequency domain determines the decay of the scar amplitudes A l∈Λs (t). As long as this broadening is sufficiently small, A l∈Λs (t) will evolve slowly with respect to the scar period 2π/ω scar . When the disorder strength is increased, their widths continue to grow, up until when W approaches W c , at which point the distributions overlap significantly and are no longer distinguishable. We expect that W c is the disorder strength for which, after a time T = 2π/ω scar , the probability |A l∈Λs (T )| 2 = e −λ l (T )T is reduced to 1/e of its value at t = 0. The disorder averaged result, obtained from Eq. (17), leads to (24) This produces the value W c ≈ 0.34ω scar for scar states near the middle of the spectrum for L = 14, which is of the order of the value W c ≈ 0.63ω scar we found numerically in Fig.7b.
There is an important and subtle question regarding the scaling to the thermodynamic limit. Even though we found numerically that the magnetization of the system did not change appreciably with system size (see Fig.7a), it was argued in [56] that perturbations that are added to a system with optimized scar states must inevitably thermalize for a sufficiently large system. Furthermore, using Lieb-Robinson bounds [56], it was argued that it would still be possible to observe non-ergodic physics for a time scale t * ∼ O − 1 2 , where is the strength of the perturbation. While in [56] the PXP Hamiltonian and the perturbation were clean systems, their argument has some bearing on our results. The reason that scars might be expected to thermalize in the disordered system is that matrix elements of the disorder operator between a scar state and another state that is ergodic must scale as D −1/2 , whereas the density of states scales as D (D being the dimension of the Hilbert space). This suggests that the factor λ l (t) in Eq. (16b) can potentially remain finite, so that non-ergodic dynamics occurs for a time scale t * ∼ O W −2 , consistent with [56]. However, higher order terms in the perturbative calculation of the amplitude might get overwhelmed by the exponential growth of the density of states. Although we do not find numerical evidence of this in this work, we cannot rule out this possibility. We do note, however, that the arguments in [56] assume that all states outside of the scar subspace are exactly ergodic. Instead, as we have found here, they realize additional scar towers which, although more highly entangled, are not completely volume-law states. This adds another subtle layer to the problem, which warrants further investigation.

VI. SPATIO-TEMPORAL CORRELATIONS
Our quantitative perturbative description of the magnetization dynamics can also be extended to provide physical understanding of spatio-temporal correlations starting from different initial states, such as clean scar eigenstates or the Z 2 product state. It was shown in Ref. [49] that the equal-time connected correlator φ (0) n | (Pσ z r (0)P) Pσ z r+R (0)P |φ (0) n c evaluated with respect to an optimized scar eigenstate converges to a nonzero value for large R, suggesting that the states are long-range spatially ordered. Furthermore, it can be argued that scars are temporally ordered as well [49]. We here consider the spatio-temporal correlator Temporal correlations can be explicitly calculated for the spatial Fourier transform C a (t) = L R,r0=1 dte iπR C a (r 0 , r 0 + R, t). Following the same steps as in the calculation of the magnetization, we find that (see Appendix D for details): where T a (t) ≈ In this expression, we defined ψ a m = s,d ψ J=1,m−sd Γ a,sd 1,m−sd with s = ±1, d = 1, 3. For convenience, in these sums we set ψ J,m = Γ a,sd J,m = 0 whenever m is out of range in the sum. In contrast to the magnetization, for the calculation of temporal correlators we need to include matrix elements beyond the tri-diagonal that are small but nevertheless contribute to the final temporal correlator. The expression Eq. (26) is completely determined, since we have already calculated [A * l (t)A l (t)] in previous sections. Similar to the magnetization, the dynamics of temporal correlators is thus also controlled by the existence of multi-tower scar resonances.
The expression Eq. (26) depends on both the choice of axis index a of the Pauli operators, and the initial state |Ψ . As a first example, we calculate the correla-tor using the clean scar eigenstate at E = 0 for magnetization component along a = z axis. In Fig.11 we show the contributions from the J = 1 and J = 1 towers for both the numerical (dots) and analytic (solid lines) values, showing again very good agreement between both. The corresponding spatio-temporal correlator [|C z (r 0 , r 0 + R, t)|] in the (R, t) plane (with r 0 an odd site) underlying these results is shown in Fig.11a. The zero-energy scar state clearly leads to a periodic pattern in the (R, t) plane indicating that spatially separated qubits are periodically correlated in time. While correlations are reduced by disorder as the temporal separation is increased, we have not observed decay as a function of R. This is counter intuitive, as disorder usually leads to a decay of spatial correlations with a characteristic length scale determined by the disorder strength.
In this case, this decay appears to be absent, suggesting that long-range order is maintained. This is quite remarkable, and deserves more detailed study as a function of system size and disorder strength. It is also reminiscent of spatio-temporal correlations in time crystals, the non-equilibrium phases of matter that spontaneously break the time-translational symmetry of the Hamiltonian [48,57,58]. Time-translation symmetry is broken when lim R→∞ C a (r, R, t) oscillates with a frequency that differs from the temporal period of the Hamiltonian and is robust to small perturbations, signifying rigidity [46]. In the scar case, the Hamiltonian has full (in contrast to discrete) time translation symmetry, and we indeed find evidence of spatial long-range order with a rigid frequency of oscillation that remains close to the clean scar frequency in spite of the presence of disorder.
While evaluating correlators starting with the scar eigenstates is illuminating, these states may not be easily accessible experimentally. For this reason, we also consider evaluating correlations starting from the |Ψ = |Z 2 initial state. In this case, connected correlators with respect to the σ z r operators vanish, so instead we analyze the a = y component. In Fig.12 we again show the contributions from the J = 1 and J = 1 towers for both the numerical (dots) and analytic (solid lines) values, showing again very good agreement between both. Notably, although oscillations continue to occur at the scar frequency, the underlying spatio-temporal pattern is very different from the pattern that we obtained starting from a scar eigenstate. Now the correlation function begins concentrated at R = 0, and spreads out as a function of time until it covers the full extent of the system. These correlations can in principle be probed in experiments. However, note that these results involve projected Pauli operators. To make even closer connection with future experiments, in the following section we evaluate them for full Rydberg Hamiltonian without any projection operators.

VII. DIAGRAM OF DYNAMICAL REGIMES FOR THE RYDBERG HAMILTONIAN
We now return to the full Rydberg Hamiltonian to examine the dynamics arising from scar resonances in the full parameter space (W, J R ). This amounts to removing all projection operators in Eq. (6). Since the Hilbert space of the Rydberg Hamiltonian grows as 2 L , we are constrained in this section to smaller system sizes L = 8, 10, 12. We focus on the connected temporal correlator This correlator differs from Eq. (26) by the absence of projection operators. This makes it a directly experimentally accessible quantity, although it has the downside that we can no longer evaluate it analytically. In Fig.13a, we show a contour plot of its Fourier transform [| C Ryd y (ω scar )|] evaluated at the scar frequency in the (W, J R ) parameter space. Oscillations at the scar frequency occur for interaction strengths greater than J R ≈ Ω and for disorder strengths up to around W ∼ W c . To pinpoint the boundary for this region, we again resort to the frequency participation parameter adapted to this correlation function Λ We find a line of parameters along which Λ (q) ω reaches a maximum, as presented in show Fig.13a with a dashed blue line. This boundary encompasses a region that is consistent with high values of [| C Ryd y (ω scar )|]. There is thus an ample regime of parameters for which it is possible to observe non-ergodic dynamics, even for moderate interaction strengths and finite disorder. Interestingly, signatures of scar resonances are present even for low interactions when J R ∼ Ω.
We further examine the ergodicity of the system using the mean level spacing ratio [ r n ]. Since the energy spectrum shifts significantly throughout the (W, J R ) parameter space, we average over the spectral range [ is the energy spread of the Z 2 state with respect to the Rydberg Hamiltonian, since this is the dynamically relevant part of the spectrum. In Fig.13b, we show a contour plot of [ r n ] in the (W, J R ) plane. A clear region is discernible for which [ r n ] is close to the ergodic value. A rough boundary for this region can be estimated to occur at the contour [ r n ] ∼ 0.45, which is the approximate value obtained in Fig.3, and is shown by a magenta dashed in line in Fig.13b. The value of [ r n ] generally increases to the left of this line as the system size is increased through L = 8, 10, 12, and decreases to the right of this line. We note that due to the limited system sizes, there are some parameters points for which it is difficult to clearly discern a direction of flow for [ r n ] as a function of L. Notwithstanding this, the boundary [ r n ] ∼ 0.45 is consistent at strong interactions with the transition at the disorder strength W T h−L we found for the PXP model. It is also consistent with the weaklyinteracting limit, where [ r n ] indeed is very close to the Poisson value, as expected for an MBL phase which corresponds to a weakly interacting disordered paramagnet that is expected to be localized. Note that, similar to what we found in Fig.3, there is a region to the right of the ergodic boundary for which, even though [ r n ] flows to lower values as a function of system size, it decreases very slowly. This is consistent with the constrained MBL phase we found in the PXP limit.
Overall, it is clear that the ergodic region in Fig.13b is larger than the region of scar resonances in Fig.13a. We thus obtain the overall diagram of dynamical regimes presented in Fig.1a. There are four regimes as a function of disorder and interaction strengths: a regime with scars resonances, another with a fully ergodic spectrum, a con- strained MBL phase, and finally an MBL phase without kinematic constraints. Of note is an unexpected behavior within the ergodic region around J R ≈ 4ω scar , where [ r n ] appears to indicate weak ergodicity. Notably, there is a slight improvement in the strength of oscillations at the scar frequency in Fig.13a in this region, perhaps suggesting that the weaker ergodicity is less detrimental for the existence of anomalous dynamics for these interaction strengths.

VIII. DISCUSSION AND CONCLUSIONS
In this work, we studied the impact of disorder on systems with quantum many-body scars. We found that scar eigenstates become resonances that have a finite life time and remain centered at the scar energies of the clean system. As a result, observables exhibit non-ergodic oscillations at the same scar frequency of the clean limit for a time scale t * ∼ O W −2 . We confirmed this picture of scar resonances by calculating explicitly the disorderaveraged magnetization and temporal correlators of the system, which match closely with the values obtained numerically. We further examined the range of interaction and disorder strengths for which scar resonances are present in the Rydberg model, and mapped out a diagram of other dynamical regimes that are accessed when scar resonances are lost at strong disorder and weak interactions.
A natural system where these results can be probed is of course the Rydberg-atom simulator used to originally discover scar states. In this system, scars were observed for J R ∼ 2π × 24 MHz and Ω ∼ 2π × 2 MHz [4]. As a result ω scar ∼ 2π × 1.3 MHz and J R ≈ 19ω scar , which places this experimental system well within the regime that has scar eigenstates along the W = 0 line in Fig.1a. The type of Zeeman disorder we considered can be implemented using drives with spatially-varying Rabi frequencies. This platform should thus make it possible to observe the changes in dynamical regimes that arise both as a function of disorder as well as interaction strengths.
Superconducting qubit systems constitute another possible quantum simulator in which to study scar resonances. In charge qubit systems, for example, it is possible to realize the model [59] which constitutes the transverse Ising model in the presence of a longitudinal field. The τ a r are Pauli matrices that act on the superconducting qubit charge states {|0 , |1 } . Up to a constant term and a unitary rotation, this model is equivalent to the Rydberg Hamiltonian Eq. (1) when Ω x = 4J xx = J R , which is what leads to the kinematic constraint in the system. In order to reach the regime of scar resonances, according to the diagram in Fig.1a, we must further have Ω z < J xx . In these systems, it is possible to reach J xx ∼ 2π × 100 MHz, as well as have a flux-tunable qubit frequency Ω z [59], so reaching the scar resonance regime is feasible. Interestingly, Ising-type models of the form Eq. (29) can be realized in trapped ions as well [60]. However, an inevitable new ingredient in such a system is the presence of power-law Ising interactions [60][61][62]. Power-law interactions can lead to nontrivial dynamics that need to be examined more closely in the context of quantum manybody scars since they could lead to qualitative changes in the diagram Fig.1a.
Finally, one could alternatively probe the effects of disorder in scar systems with digital quantum simulators, although there are some challenges for currently available hardware [63]. However, there has been significant progress in the fidelity of single and two-qubit gates in both ion traps [64,65] and SC qubits [66]. Furthermore, there has been progress in algorithms for quantum simulation beyond the coherence time of quantum devices [67]. As a result, it might be possible to access the physics of scar resonances in some devices in the near future. Interestingly, in addition to having the potential to model the Rydberg Hamiltonian, digital simulators can directly simulate the strongly interacting limit described by the PXP model using Toffoli gates that effectively impose the kinematic constraint on sets of three contiguous qubits.
The results presented in this work invite further investigation into the non-ergodic nature of scar systems along several directions. For example, the rigidity we found in the spatio-temporal correlations of the PXP model poses intriguing questions regarding a possible framework for understanding quantum scar systems as realizations of continuous time translation symmetry breaking in static (undriven) systems. It would also be interesting to explore further the transitions at strong disorder that eventually lead to the fully thermal and constrained MBL phases. Experimentally, our findings could be used to devise strategies to calibrate quantum simulators that suffer from spatial imperfections such as in SC qubits.
Finally, note that while our work has focused on the one dimensional Rydberg model as a case study, scar resonances can arise in other models and also in higher dimensions. This is clear from the derivation of Eq.17, which did not make assumptions about the specific model at hand. The results in this work thus apply quite generally to other systems and consequently lay the groundwork to study other aspects of non-ergodic dynamics arising from scar states in the presence of disorder. , where we defined ω l l = E l . Thus, the problem reduces to determining the functional form of φ (0) l | σ a π |φ (0) l . As was found in [25], angular momentum operators can be defined within the scar subspace as S y = 1 2η r (−1) r P s (σ y r + δσ y r ) P s , where P s = l∈Λs |φ (0) l φ (0) l | is the projector into the subspace of optimized scars, and δσ a r = 2δJ R Ω σ z r+2 + σ z r−2 σ a r . These operators approximately satisfy angular momentum commutation relations within the subspace of scar states, so that the L + 1 states span the tower of maximum total angular momentum S = L/2. Now, note that S x and S y are slight deformations of the operators σ x 0 and σ y π , respectively, since δJ R /Ω is small. As a result, we approximate We emphasize that, while we have disregarded the deformation term proportional to δJ R /Ω in these angular momentum operators, the eigenstates {|φ (0) l } are obtained with the deformed Hamiltonian. This reduces the deviation from the ideal su(2) algebra found in [25].
Next, the optimization of scar states implemented in [25] is designed so that the L + 1 states {|K n } obtained by the Forward Scattering Approximation (FSA) proposed in [6] become exact eigenstates of S z so that K n |S z |K n ≈ m n δ nn , where m n = −L/2, . . . , L/2. As a result, the projection operator into the scar subspace can be written in the unitarily equivalent form P s = n |K n K n |, since the {|K n } must span this subspace. These FSA states are constructed by repeated application of the operator H + = r (σ x r + δσ x r ) + i r (−1) r (σ y r + δσ y r ) on the Z 2 state [25]. Since [σ z π , H + ] = 2H + and |Z 2 is an eigenstate of σ z π , then K n |σ z π |K n = 2m n δ nn . Thus, we can write approximately Having established the approximate su(2) algebra satisfied by σ x 0 and σ y,z π , we now express the latter in the form P s σ z π P s = P s σ + s + σ − s P s (A4) P s σ y π P s = ηP s σ + s − σ − s P s .
Here, the σ ± s act as ladder operators on the scar eigenstates: φ (0) l ∈Λs |σ ± s |φ (0) l∈Λs ∝ δ l ,l±1 where the notation l ±1 refers to the index of the scar state that is spectrally separated by ±ω scar from the l-th scar state. Thus, the matrix elements φ (0) l ∈Λs | σ a π |φ (0) l∈Λs precisely pick out the phases that oscillate at the frequency ω l±l = ω scar . We then obtain the magnetization components M y (π, t) ≈ ηL sin (ω scar t) , M z (π, t) ≈ −L cos (ω scar t) , To arrive at this expression, we have used that where Λ s excludes the highest-energy scar state. These are the expressions presented in Eq.13 of the main text.

Appendix B: Reordering of basis
In this appendix, we discuss how to re-order the energy eigenstates into towers. Suppose we take an energy eigenstate |φ (0) m1 from the bottom of the spectrum and we apply the ladder operator once. We use the operators defined in the main text As we discussed in Appendix A1, these are approximate ladder operators within the space of optimized scar states. Here, however, we will examine their behavior in the full PXP Hilbert space. In general, applying σ + s on any given state can result in a linear combination over all eigenstates of the clean system. However, numerically one finds that this linear combination is highly concentrated at a particular state |φ m2 . By repeating this procedure, we eventually arrive at a state |φ (0) mmax for which the further application of the operator does not produce a peaked distribution. Thus, σ + s acts like a ladder operator for the set of states {|φ (0) mi }, even in cases where the states do not belong to the scar subspace. By iteratively applying this procedure to all eigenstates, we can systematically group all energy eigenstates into sets {|φ (0) Jm }, where J labels the set and m = 1, . . . , D J labels the state within that set. In this notation, we can include the scar space in the set J = 1, since this set is generated by starting with the ground state. We found that a small number of states do not have strong matrix elements with any other state in the basis, so they form sets by themselves according to our algorithm. Clearly, such states do not contribute to the magnetization. This is the set of steps we took to re-order the basis of energy eigenstates. For example, the case L = 8 presented in the main text leads to 13 towers of sizes: D 1 = 9, D 2 = D 3 = 7, D 4 = D 5 = 6, D 7−13 = 1. ρ n0(l )n0(l ) τ 2 − ρ n0(l),n(l) τ 1 . Note that we can write