Slow Thermalization of Exact Quantum Many-Body Scar States Under Perturbations

Quantum many-body scar states are exceptional finite energy eigenstates in a thermalizing system that do not satisfy the eigenstate thermalization hypothesis. We investigate the fate of exact many-body scar states in the so-called PXP model for the Rydberg blockaded atom chain under perturbations. At small system sizes, deformed scar states described by perturbation theory survive. However, we argue for their eventual thermalization in the thermodynamic limit from the finite-size scaling of the off-diagonal matrix elements. Nevertheless, we show numerically and analytically that the nonthermal properties of the scars survive for a parametrically long time in quench experiments. We present a rigorous argument that lower-bounds the thermalization time for any scar state as $t^{*} \sim O(\lambda^{-1/(1+d)})$, where $d$ is the spatial dimension of the system and $\lambda$ is the perturbation strength.

Wigner pioneered the application of random matrix theory to describe quantum chaos [1]. Nowadays, it has even become a definition of quantum chaos: the applicability of the random matrix theory description to the statistical properties of the spectrum and wavefunctions of a quantum mechanical system, in both single-particle and many-body quantum systems. Based on the random matrix theory description, Srednicki and Deutch proposed the Eigenstate Thermalization Hypothesis (ETH), bridging the concept of quantum chaos and the validity of statistical mechanics in closed quantum many-body systems [2,3]. Essentially, ETH implies that the reduced density matrix of a single eigenstate is equal to that of the microcanonical/canonical ensemble, and it is how statistical mechanics emerges in closed quantum systems.
The strong version of ETH proposes that every eigenstate satisfies the above property [4][5][6][7]. However, there can be some exceptional states at a finite energy density that do not satisfy the ETH, while the other states do. Such states are dubbed quantum many-body scar states, in analogy with the single particle scar states [8].
Current analytical understanding about quantum many-body scars relies on the identification of certain exact eigenstates. Moreover, the Hamiltonians with exact scar states are at some special point in some parameter space or have the embedded Hamiltonian structure. An immediate question arises: how robust are the exact quantum many-body scar states under generic perturbations? Are the exact scar states discovered and constructed in several models useful to understand the physics once we add perturbations?
In this paper, we address the above questions by studying the fate of the exact scar states in the PXP model under perturbations. We show that, in the finite-size ED data, there is some apparent robustness in the nonthermal signatures of the exact scar states upon perturbation. The perturbed eigenstates in finite sizes can indeed be understood using standard perturbation theory. However, the finite-size scaling of the matrix elements between the scar states and other eigenstates suggests the eventual thermalization of the scar states at larger system sizes.
Despite the eventual thermalization of the scar states, the thermalization rate is parametrically slow in the strength of the perturbation. In particular, we consider the quench dynamics starting from the exact scar states evolving under the perturbed Hamiltonian. We rigorously lower-bound the time scale for any local observable to thermalize by t * ∼ O(λ − 1 1+d ), where λ is the strength of the perturbation and d is the spatial dimension of the system. Our results suggest that the exact scar states discovered or constructed in different special models can indeed be used to understand the persisting dynamical signatures under perturbations up to some parametrically large time scale and in the thermodynamic limit.
Our results suggest a strong analogy between weakly perturbed scarred systems and weakly perturbed integrable systems. In the latter case, the eigenstates in finite sizes can be understood as perturbed from the special integrable Hamiltonian. Analogously, the deformed scar states in finite sizes are perturbatively related to the ex-act scar states of the corresponding special Hamiltonian. Moreover, in the thermodynamic limit, a weakly perturbed integrable Hamiltonian prethermalizes to a generalized Gibbs ensemble that persists for a parametrically long time scale [36,37]. Similarly, the nonthermal properties of the exact scar states can also survive in quench experiments under perturbations to some parametrically long time even in the thermodynamic limit. We therefore propose that there is a connection between "partially solvable" Hamiltonians and observed quantum many-body scar states.
The paper is organized as follows. In Sec. II, we introduce the PXP model and its properties, and also motivate the specific type of perturbation in our consideration. In Sec. III, we study the signatures of the scars under the perturbation in finite size systems, demonstrating their apparent robustness. In Sec. III A, we apply standard perturbation theory to describe the deformed scar states at small perturbation strength and numerically accessible system sizes. However, in Sec. III B, we argue that the deformed scar states eventually hybridize with the ETH-satisfying states at nearby energies by studying the finite-size scaling of the matrix elements connecting the scar states to other states. Thus, the deformed scar states lose their nonthermal character as L → ∞ and eventually satisfy the ETH. Despite their eventual thermalization, in Sec. IV A, we numerically show that the thermalization is slow under global quenches, and prove a rigorous lower bound on the thermalization time scale in Sec. IV B. In Sec. IV C, we discuss possible scenarios where the bound on the thermalization time can be stronger. We conclude and discuss our work in Sec. V.

II. MODEL
We consider the one-dimensional (1d) PXP model with perturbation. The PXP model is the effective constrained model describing the dynamics of the Rydberg atom chain in the regime of the nearest-neighbor blockade. More specifically, we consider a 1d atom chain with size L and open boundary conditions. There are two degrees of freedom at each site: |0 (atomic ground state), and |1 (atomic excitation). The blockade condition excludes all configurations | . . . 11 . . . with adjacent atomic excitations from the Hilbert space. The dimension of the Hilbert space grows as The perturbed PXP model has the following Hamiltonian:Ĥ whereĤ Here P ≡ |0 0| and X ≡ |0 1| + |1 0|. Before specify-ingV , we summarize some properties ofĤ 0 , in order to motivate the specific perturbation we will consider later. These properties of the PXP model have been discussed in detail in Refs. [10][11][12][13][14]. First,Ĥ 0 has the inversion symmetryÎ : j → L−j +1. It also has the property that if we define the particle-hole transformationĈ where Z ≡ |1 1| − |0 0|, thenĈĤ 0Ĉ = −Ĥ 0 . This implies that, for any eigenstate |E ofĤ 0 with energy E = 0,Ĉ|E is also an eigenstate with energy (−E). The combination ofÎ andĈ guarantees the exponential degeneracy of E = 0 states inĤ 0 . This is due to the difference between the number of states with C = 1 and C = −1 in each inversion symmetry sector I. As a result, in the E = 0 manifold, the eigenstates with a particular inversion symmetry quantum number I = ±1 will also have a definite particle-hole quantum number given by C = I.

A. Signatures of the exact scar states
In Ref. [14], four exact scar states were discovered inĤ 0 for even system sizes L, labeled as |Γ αβ , where α, β ∈ {1, 2}. Here we assume that the states are normalized. For the readers' convenience, we summarize the wavefunctions and some essential properties of these states in Appendix A. In particular, |Γ 11 and |Γ 22 have energies E = 0; the states |Γ 12 and |Γ 21 have energies E = √ 2 and E = − √ 2 respectively. Despite being in the middle of the spectrum at energy density corresponding to infinite temperature, the states |Γ αβ have constant bipartite entanglement entropy scaling with the subsystem length ("area law"), instead of the volume law scaling required by the ETH. Moreover, the expectation values of some local observables in these states do not agree with the thermal ensemble values at infinite temperature, therefore violating ETH. One intriguing feature is their valence bond solid (VBS) order, which is also used to identify their translation symmetry breaking. The order parameter is defined aŝ whereD j,j+1 ≡ |01 10| + H.c.
detects the dimer (bond) strength. In the thermodynamic limit, in the bulk, Γ αβ |D j,j+1 |Γ αβ = 0 if j is odd and −2/9 if j is even. On the other hand, an infinite temperature thermal ensemble would give the thermal value 1 D L Tr[D j,j+1 ] = 0 for all j, where D L is the dimension of the Hilbert space of system size L.

B. Particle-hole odd perturbation
We are interested in the fate of the above scar states and their signatures under perturbation. Moreover, we also want to examine if the E = 0 degenerate manifold has any relevance to the robustness of the scar states. Thus, we examine the states |Γ 21 and |Γ I ≡ (|Γ 11 is the normalization factor (see Appendix A). Both states have inversion quantum number I = −(−1) L/2 , and |Γ I additionally has the particle-hole quantum number C = −(−1) L/2 . To compare their behavior against some thermal or chaotic state, we will also examine the state |Γ th , which is picked as the eigenstate with eigen-index three more than the index of |Γ 21 in the I = −(−1) L/2 symmetry sector.
The simplest inversion-symmetric perturbation that has the propertyĈVĈ = −V iŝ This perturbation was first studied in Ref. [12]. It was identified that at λ ≈ −0.02,Ĥ is close to some unknown integrable point. Moreover, in the periodic boundary condition version ofĤ, at λ ≈ −0.053, a set of nearly perfect scar states was numerically found to lead to nearly perfect revivals [28] [38].

III. SIGNATURES OF THE SCAR STATES UNDER PERTURBATION
To see how robust the scar states are under the perturbation Eq. (6), we first examine the loss of the fidelity as we increase λ in an open chain of length L = 20. In Fig. 1 we examine the overlaps of exact eigenstates of the perturbed Hamiltonian with the unperturbed states, | E n (λ)|Γ | 2 , where |E n (λ) runs over eigenstates ofĤ = H 0 + λV , while |Γ is |Γ 21 , |Γ I , or |Γ th in panels (a), (b), or (c) respectively. In Fig. 2, we examine the bipartite entanglement entropy of |E n (λ) and the VBS order parameter E n (λ)|M |E n (λ) (associated with the unperturbed scar state |Γ 21 ) measured in the exact eigenstates of the perturbed Hamiltonian as a function of the perturbation strength.
In Fig. 1(a), | E n (λ)|Γ 21 | 2 clearly shows some apparent robustness under the perturbation [the numerical values can be seen in Fig. 4(a)]. In the region λ > 0, there is a single perturbed eigenstate which can be traced back to |Γ 21 , despite multiple avoided level crossings when the perturbation strength is increased. On the other hand, on the λ < 0 side, |Γ 21 has significant overlap with several states near the approximate integrable point λ ≈ −0.02 [12]. The spread in overlap makes it difficult to identify a single perturbed state associated with the unperturbed state |Γ 21 . As a comparison, in Fig. 1(c), we show the fidelity loss of |Γ th . We can see that the overlaps of |Γ th spread over multiple states at a smaller value of λ as compared to those of |Γ 21 . This suggests that |Γ th hybridizes more strongly with the other states in the spectrum upon perturbation.
In Fig. 1(b), |Γ I appears to be even more robust under the perturbation than |Γ 21 . Here, the overlaps between the |E n (λ) = 0 states and |Γ I are summed up and displayed in the figure. That is, we show the weight of |Γ I projected into the perturbed E = 0 manifold. The numerical values can be seen more easily in Fig. 4(b).
In Figs. 2(a) and (b), we examine the bipartite entanglement entropy and the VBS order of the exact perturbed eigenstates near the unperturbed scar state |Γ 21 under the perturbation. Again, on the λ > 0 side, the perturbed eigenstates corresponding to the highest overlap with |Γ 21 also shows low entanglement entropy and significant VBS order. Note that at λ ≈ −0.02, the entire spectrum shows features of low entanglement, which is again a manifestation of the proximity to some integrable point [12].

A. Perturbation theory
Here we show that the eigenstates which have high overlaps with the unperturbed states |Γ 21 or |Γ I are perturbed version of |Γ 21 or |Γ I . The fact that such a perturbation theory is controlled for our system sizes is a manifestation of the effective weakness of the perturbation as far as the scar states are concerned, which we will further quantify below. However, we expect this perturbation theory to fail for large enough L for any λ = 0 (see Sec. III B).
The standard non-degenerate perturbation theory to N -th order gives the perturbative corrections to the unperturbed energy E For |Γ 21 , we numerically calculate the perturbed energies E [N ] pert. , N = 1, 2 and 3 and compare them to the ED energy E exact as shown in Fig. 3. The ED energy E exact is the energy of the eigenstate that has the maximum overlap with |Γ 21 . From the figure, we see that the perturbed result agrees very well with the ED result for small the unperturbed scar states, | En(λ)|Γ21 | 2 and | En(λ)|ΓI | 2 respectively, and (c) the unperturbed (presumably) thermal state, | En(λ)|Γ th | 2 . The horizontal axis is the perturbation strength λ while the vertical axis shows eigenstate energies. It appears that the scar states |Γ21 and |ΓI hybridize with other states relatively weaker compared to the thermal state |Γ th , and exhibit some robustness in the finite-size ED spectrum. Note that for this system size when going from λ = 0 to λ = 0.01, the scar state |Γ21 crosses roughly 40 states while still maintaining its fidelity. At λ ≈ −0.02, the system is near some approximate integrability point. perturbation strength λ. Note that as mentioned previously, the weight of |Γ 21 is spread over several states on the λ < 0 side, making it hard to single out a particular eigenstate for the energy comparison, resulting in some non-smoothness in E exact .
For |Γ I , we need to use degenerate perturbation theory since |Γ I resides in the E = 0 degenerate manifold. However, if we choose the eigenbasis with definite inversion quantum number, then this basis is already appropriate for non-degenerate perturbation theory. Indeed, an important corollary of the inversion symmetry ofV and of the propertyĈVĈ = −V is that E n = 0|V |E m = 0 = 0 for any eigenstates n and m in the zero-energy manifold. The reason is thatV changes the particle-hole quantum number but preserves the inversion symmetry. As discussed in Sec. II, the E = 0 states have definite pairs of particle-hole and inversion quantum number. This also implies that, to first order, the E = 0 states are not hybridizing with each other. Furthermore, using this property, we can prove that the perturbed wavefunction |Γ (n) I will have definite particlehole quantum number C and hence E (n) Γ I = 0 to any order. We present the proof in Appendix. B. However, we emphasize that the proof does not imply the convergence of the perturbation theory in the thermodynamic limit.
Having established good agreement between the perturbed energy and the ED energy, we further show that the perturbed wavefunction is a good description for the ED wavefunction. The (unnormalized) perturbed wavefunction to first order is given as We use this wavefunction to provide quantitative understanding of the fidelity loss obtained in ED. Specifically, we calculate and compare this with the ED result F exact (λ; Γ) ≡ | E n (λ)|Γ | 2 (where the ED state is selected by maximizing the overlap). In Figs. 4(a) and (b), we compare the perturbation theory with the ED results, for |Γ 21 and |Γ I respectively. In Fig. 4(a), the ED result has some sudden jumps in the overlap on the λ > 0 side due to accidental events of (avoided) "level crossings." When there is another eigenstate which happens to be very close in energy, the weight of |Γ 21 will be spread over these states-known as accidental hybridization. However, once the level crosses, the overlap recovers. In fact, this phenomenon stems from the fact that Γ 21 |V |Γ 21 behaves differently from the nearby states n (0) |V |n (0) . This is also only possible when there are scar states present, since ETH would predict n (0) |V |n (0) to have the same value (up to some finite size correction) for |n (0) 's with the same energy density. On the other hand, for λ < 0, the weight of |Γ 21 is almost always spread over several states. It is therefore less justified to single out a particular special state for the overlap comparison.
Turning to Fig. 4(b), since there is no level-crossing through the E = 0 manifold, the ED result of the overlap to |Γ I is smooth in λ. We therefore see that, up to the accidental hybridization, the perturbation theory provides a good understanding for the perturbed ED scar states in this system size for both |Γ 21 and |Γ I .
We note that technically, the Rayleigh-Schrodinger perturbation theory is not applicable when n (0) |V |Γ It is expected to already be not valid when there are (avoided) level crossings, as observed in the case of |Γ 21 . The good agreement between the eigenenergy E exact and the perturbation theory E pert. given the multiple level crossings is therefore indeed remarkable. It reflects the smallness of the off-diagonal matrix elements while the crossing levels are "pushed through" (as a function of λ) by the differing diagonal matrix elements. To better explain the relative accuracy of the perturbation theory for this system size, we present a slightly modified version in Appendix C.

B. Finite-size scaling and eventual thermalization
Despite the success of the perturbation theory in describing L = 20 ED results, here we study the finitesize scaling of the off-diagonal matrix elements and argue that the exact scar states will hybridize with the FIG. 3. The energy difference between the perturbation theory Epert. (1st to 3rd order) and the ED result Eexact, where Eexact is the energy of the eigenstate with the maximum overlap with |Γ21 . The perturbation theory gives accurate predictions of the energy of the perturbed eigenstate in the range of small perturbation strength λ. Inset: The energies obtained from ED compared to the perturbation theory predictions up to 3rd order. The dominant shift in the energy is already captured by the 1st order, i.e., the "diagonal" part of V , which also suggests relative weakness of the off-diagonal matrix elements. See text for the discussion of the accuracy of the perturbation theory and also Appendix C for an improved treatment incorporating the diagonal part of V in the unperturbed part.
other ETH-satisfying states and thermalize eventually. First, we examine the distribution of the amplitudes of the off-diagonal matrix elements | n (0) |V |Γ | for each of the three states |Γ = |Γ 21 , |Γ I , and |Γ th . This is plotted on the left side in Fig In the left panels, we see that for | n (0) |V |Γ 21 | and n (0) |V |Γ th , there is a state in each case with zero overlap (showing up at numerical error threshold ∼ 10 −14 ) withV |Γ 21 orV |Γ th respectively. These states are just C|Γ 21 andĈ|Γ th respectively. One can easily see that Γ|ĈV |Γ = − Γ|VĈ|Γ = − Γ|ĈV |Γ * = 0, where in the last equality we used time-reversal-like symmetry, or the fact that Γ|VĈ|Γ is a real number. On the other hand, for |Γ I , we see two states with zero matrix elements n (0) |V |Γ I . These states are in fact |Γ 21 and |Γ 12 . However, these matrix elements are zero due to the special structure of the states, instead of some symmetry reasoning. Also note that this is not reflected in the plot of | n (0) |V |Γ 21 |, because the |E Comparison of the ED and the 1st-order perturbation theory results for the overlap between the descendant states and the unperturbed scar states, for (a) |Γ21 and (b) |ΓI (in ED, the descendants are the eigenstates with the largest such overlaps) at L = 20. The sudden drops in (a) on the λ > 0 side in the overlap in ED are caused by accidental hybridization that occurs very close to the avoided level crossings. On the λ < 0, the weight of |Γ21 is spread over several states, resulting in low overlap. The perturbation theory provides a good understanding for the eigenstate at small λ.
are obtained through numerical diagonalization, and |Γ I is a superposition of these states.
For |Γ th , we see that the distribution of | n (0) |V |Γ th | is fairly uniform and does not have any features within the energy window |E On the other hand, we see that for |Γ 21 and |Γ I , there is some "horn" structure at |E 6. This is due to the strong connection between the exact scars and other quasiparticlelike excitation states on top of the exact scar states. (See Ref. [14,16] for such a "quasiparticle picture" of the other nonexact scar states in the PXP model in the spirit of single mode approximation and its multi-mode generalization.) In all three cases, there is a rapid dropoff of the matrix elements for |E Γ th | 5, which reflects the locality of the Hamiltonian-see Appendix D for details.
We can see some suppression in the amplitudes of the matrix elements | n (0) |V |Γ 21 | and | n (0) |V |Γ I | near zero |E We think that for these system sizes, this suppression in the off-diagonal matrix element amplitudes further assists the validity of the perturbation theory description in Fig. 4 discussed earlier. However, by examining the distributions in more narrow energy windows, we find that this only affects the connection to the thermal states quantitatively and not qualitatively.
We further use the statistics of | n (0) |V |Γ | to argue for the eventual thermalization of the original unperturbed eigenstates. For an ETH system under a perturbation V , within some energy window |E The distribution of the matrix elements | n (0) |V |Γ |, with |Γ chosen as |Γ21 , |ΓI , and |Γ th , shown for system sizes L = 12 and L = 22. In all cases, there is a rapid falloff outside some energy window because of the locality ofĤ0 and V . In the case of the scar states, we see strong "horns" at |E Γ | ≈ ±2.6 (attributed to matrix elements to other scars in the spectrum) and also some suppression for small |E Γ |, while in the thermal state case the distribution is more uniform. Nevertheless, the typical values of the matrix elements decrease with L comparably for the scar and thermal states. In the right panels, we zoom in to show more detailed features of the distribution. expects the off-diagonal matrix elements n (0) |V |m (0) to scale as D −1/2 L , where D L is the dimension of the Hilbert space for the chain of length L (D L ∼ φ L for the Rydberg-blockaded chain). The reason is that if |n (0) is thermal or chaotic, than |n (0) behaves essentially like a random vector in the Hilbert space when considering offdiagonal matrix elements. We can viewV |m (0) as some fixed "direction" in the Hilbert space, and the overlap between a random vector and a vector in some direction will be of order D −1/2 L . Note that we needed only the state |n (0) to be "thermal" (i.e., essentially "random") while the other state |m (0) can be either thermal or non- Γ | is among the smallest 400. For all statesthermal or scar-the behavior of the off-diagonal matrix elements satisfy the scaling predicted by the random matrix theory.
thermal as long asV |m (0) is not somehow special, e.g., |m (0) does not happen to be an exact eigenstate ofV . On the other hand, the density of states at energy density corresponding to infinite temperature grows as ∼ D L . Therefore the ratio between the typical off-diagonal matrix elements and the level-spacing will grow as D 1/2 L , which signals the strong hybridization or thermalization.
A well known counterexample which circumvents thermalization is many-body localization [39][40][41][42][43][44][45]. In such systems, the off-diagonal matrix element for nearby states scales as e −L/ξ , where ξ is some localization length. In a crude estimate, if ξ < L/ ln D L , the ratio between typical off-diagonal matrix element and level-spacing will in fact decrease with L [46]. The eigenstates of manybody-localized systems are thus perturbatively accessible at any L.
Returning to our setting of following the fate of the scar states under the perturbation, we have already observed that they generically have nonzero matrix elements to all eigenstates with the same quantum numbers, and that they appear to lose fidelity upon increasing the perturbation strength or the system size. We also expect that the scar states are exceptions in the unperturbed PXP model while most of the eigenstates are ETH-satisfying. The above scaling of the matrix element between a scar state and an ETH state with L should therefore hold. Nevertheless, we would like to check this explicitly and examine more detailed quantitative features of such con-nections. Accordingly, we next study the finite-size scaling of the off-diagonal matrix elements.
In Figs. 6(a) and (b), we show the statistics of the off-diagonal matrix elements. In Fig. 6(a), we averaged | n (0) |V |Γ | in the energy window |E n = 0 states in the |Γ I case. Not surprisingly, | n (0) |V |Γ th | shows the φ −L/2 scaling predicted by the random matrix theory description. Importantly, we also see the same scaling for the scar states | n (0) |V |Γ 21 | and | n (0) |V |Γ I |, but with smaller amplitude than for the thermal state. To avoid the arbitrariness of specifying the energy window for averaging, in Fig. 5(b), we take 400 matrix elements corresponding to 400 states with the smallest |E (We do not show L = 12 and L = 14 data because 400 states is larger than the Hilbert space dimension for L = 12 and larger than half of the Hilbert space dimension for L = 14.) For both averaging protocols, we see that the off-diagonal matrix element scaling of the thermal |Γ th and scar states |Γ 21 and |Γ I all satisfy the random matrix theory prediction, while the overall amplitudes for the scar states are smaller. This also explains the apparent weaker hybridization observed in Fig. 1: the apparent robustness of the scars compared to the thermal state is due to the smaller amplitude of the off-diagonal matrix elements.
The relevant quantitative difference is even stronger than suggested by the averages in Fig. 6 because the number of states being averaged covers a wider window than the "dip" close to zero |E Fig. 5, while the states inside this dip play a larger role in the perturbed wavefunction at these sizes. As mentioned earlier, looking more closely at the matrix elements inside the dip (e.g., averaging over smaller number of states, all the way down to just few closest states), we do not see any faster decay with L than expected in the random matrix theory.
To conclude, from the above scaling argument comparing the off-diagonal matrix elements and the levelspacing, we expect the scar states will eventually thermalize in the thermodynamic limit. Finally, we also note that comparison between the |Γ 21 and |Γ I cases suggests that the degenerate E = 0 manifold (to which the latter state belongs) is in fact irrelevant to the hybridization of the scar states.

IV. SLOW THERMALIZATION OF LOCAL OBSERVABLES
Despite the eventual thermalization of the scar states under perturbations, in the following, we show numerically and analytically, that the signatures of the scars can survive up to some long time set by the perturbation strength even in the thermodynamic limit. This slow thermalization is analogous to that in systems that are weakly perturbed from integrability [36,37] or that have weak breaking of some conservation laws [47]. Specif-ically, we consider the following general global quench setting: we consider the initial state as an exact scar state |Γ of some HamiltonianĤ 0 , and let it evolve under the perturbed HamiltonianĤ =Ĥ 0 + λV . Since |Γ is a scar state, there is some local observablem which is nonthermal. We then examine how quickly such an observable behaves at late times, and if it reaches its thermal value.
A. Slow decay of VBS order in numerics First, we numerically show the slow thermalization using time-evolved block decimation (TEBD) method [48] [49]. We choose the initial state as |Γ 21 or |Γ I and evolve it underĤ =Ĥ 0 + λV , whereĤ 0 andV are given in Eqs. (2) and (6)  At short times, the VBS order deviates from the initial value as t 2 . This can be easily understood from the timereversal-like symmetry. Considering the Taylor series the coefficient of t 1 is given by Now, Γ(t)|m|Γ(t) is a real number by the hermiticity of the observablem, while Γ|[Ĥ,m]|Γ is a real number because the operatorsĤ andm have real-valued matrix elements in the Rydberg atom basis, and the wavefunction |Γ has real-valued amplitudes in this basis. Hence, we conclude that d dt Γ(t)|m|Γ(t) | t=0 = 0. The coefficient of the t 2 growth is this is generically not zero, and its sign determines if Γ(t)|m|Γ(t) curves up or down initially. The initial t 2 behavior stops at around t ≈ 1, and beyond this time the VBS order appears to relax almost linearly. (The oscillations that are clearly visible on top of the overall decay are roughly at frequency ω ≈ 3 and can be traced to the "horn" features in the distribution of the matrix elements discussed in Fig. 5.) At small λ = 0.02, the VBS order has an almost negligible relaxation rate. As λ increases, the relaxation of the VBS order becomes visible within the time window shown in the figure. Note that in Figs. 1 and 2, we study the perturbation strength |λ| ≤ 0.05. Here we also study larger perturbation strengths to have noticeable VBS order relaxation.
Given the almost linear relaxation behavior of the VBS order, we empirically extract the relaxation rate by fitting the data via the least-squares method to in the time interval t ∈ [1,30]. Fig. 7 (e) plots the relaxation rate W λ . The relaxation rate appears to vanish faster than linearly at small λ, which corresponds to a relaxation time that diverges faster than λ −1 . However, since the slope obtained from the linear fit at the smallest |λ| ≤ 0.05 are numerically very small and may not be very reliable, we cannot extract the precise functional form.
In Fig. 7 (f), we show the snapshots of the dimer strength pattern Γ 21 (t)|D j,j+1 |Γ 21 (t) . It clearly shows that the dimer strength pattern is uniform on even/odd sites, som is indeed representative of the VBS order.
To understand the effects of finite size, in Fig. 8, we show the VBS order decay Γ 21 (t)|m|Γ 21 (t) for different system sizes. We can see that the traces of the small system sizes L = 18 and L = 30 are converging to the trace of L = 48 up to the time of about 10 or slightly larger. We think that this time is what is often referred to as "recurrence" time in ED studies of quantum many-body dynamics, which is roughly the time it takes for information to propagate across the entire system [50,51]. In this picture, the L = 48 data up to such time t 10 is already representative of the thermodynamic limit. We therefore think that the almost linear relaxation behavior observed up to this time is already representative of the thermodynamic limit. The relaxation is clearly visible for λ ≥ 0.1, and we expect eventual thermalization for such perturbation strengths; the relaxation is not visible for the smaller λ = 0.02 and 0.05. Although we do not expect any qualitative changes as we vary λ, the characteristic relaxation time scale can grow as λ approaches zero. We therefore argue that our exact scar states eventually thermalize for any nonzero λ.
It is interesting to note that all sizes in Fig. 8 show qualitatively similar relaxation behavior that continues well beyond the above estimated recurrence times for these sizes. (For the largest perturbation strength |λ| = 0.2, the systematic relaxation for the smallest size L = 18 stops around t 50, beyond which time the trace wanders non-systematically.) At present, we do not have a good understanding of this observation. Nevertheless, it suggests that the presented time range in Figs. 7 and 8 may be representative of the thermodynamic limit behavior beyond the recurrence time t 10.
Finally, we note that the dynamical signature of the VBS order is essentially the same starting from either |Γ 21 or |Γ I for large enough system sizes. This is therefore another piece of evidence suggesting that the E = 0 manifold does not have significant effects on the observable dynamical signatures of the scar states.

B. Rigorous bound on the thermalization time
We have shown numerically that small perturbation strength indeed gives slow thermalization of the VBS order. However, as in most numerical calculations, the results may be strongly affected by finite-size effects, and one may worry that the slow thermalization may not be a true phenomenon in the thermodynamic limit. Remarkably, in the following, we give a rigorous lower bound on the thermalization time for the scar states, even in the thermodynamic limit, that diverges when the perturbation strength goes to zero. The following theorem in fact is valid in any dimension, though we will specialize it to one dimension first.
Recall that we consider the global quench setting aŝ We also further assume that bothĤ 0 andV are local Hamiltonians, i.e., lattice sums of local terms, e.g.,V = jv j wherev j acts only on degrees of freedom near site j.
Consider a local observablem, which we assume to be localized near the origin j = 0 (at the middle of the chain), and its expectation value, Γ(t)|m|Γ(t) . Consider the time derivative where in the last line we used that |Γ is an eigenstate of H 0 . For t = 0, the right hand side is given by the expectation value of [λV ,m] in the initial state |Γ . Since this is a local operator, the expectation value is λ times an O(1) number, and the initial rate of change of the observable is thus of order λ. (In the specific quench setting considered in the previous subsection, this term is actually zero because of the time reversal invariance ofĤ, but we will not use this in the developments below.) The small parameter λ is present as a factor in Eq. (18) at any time t, but we have to consider the possibility that it multiplies a function of t that grows with time.
A conservative bound can be obtained by noting that the time-evolved operator e iĤtm e −iĤt is significantly spread only over region |j| ≤ v LR t, where v LR is the Lieb-Robinson velocity for the Hamiltonian H [52][53][54]. Therefore, this time-evolved operator can have a significant commutator only with local terms in V that are inside this region, i.e.,v j with |j| ≤ v LR t. Now, for any j, the norm of the operator [λv j , e iĤtm e −iĤt ] is bounded by 2 λv j m , which is λ times an O(1) number. Hence, the total contribution to Eq. (18) from |j| ≤ v LR t is bounded by λ(c 0 + c 1 t), where c 0 and c 1 are fixed O(1) numbers (and for concreteness we take the spatial dimension to be d = 1). On the other hand, using the Lieb-Robinson bounds [v j ,m(t)] ≤ c exp[−a(j − v LR t)], the total contribution from |j| > v LR t can be bounded by λc 0 , where c 0 is also a fixed O(1) number. Thus, we have Γ|[λV , e iĤtm e −iĤt ]|Γ ≤ λ(c 0 + c 1 t) .
A minor point: The Lieb-Robinson bounds are used here with respect to the full HamiltonianĤ =Ĥ 0 + λV , which depends on the parameter λ. However, typically used estimates-effectively, upper bounds-of the Lieb-Robinson velocity use operator norms of local terms in the Hamiltonian and will depend smoothly on λ, which will only introduce a smooth λ dependence in the parameters c 0 and c 1 . Hence, the leading λ dependence of the above bound is essentially unchanged. One can provide a slightly different argument using the Duhamel formula [55] where this minor point does not arise at all.
The analysis so far is completely general and applies to anyĤ 0 and any initial eigenstate |Γ . Let us now see how we can use it to lower-bound the persistence time of nonthermal properties when |Γ is a scar eigenstate ofĤ 0 . In this case, Γ|m|Γ differs by a finite amount from the expectation value ofm in the nearby "thermal" eigenstates, which we will denote as m th.;Ĥ0, E0= Γ|Ĥ0|Γ . We assume that |Γ(t) will eventually "thermalize" with respect toĤ =Ĥ 0 + λV . (This is the "worst case" scenario, as suggested by the finite-size scaling analysis in Sec. III B. Otherwise, the nonthermal persistence time is infinite.) The eventual expectation value is then By perturbation theory in equilibrium quantum statistical mechanics, barring the unlikely situation whereĤ 0 happens to have a first order transition at the temperature corresponding to average energy E 0 , we expect that m th.;Ĥ, E= Γ|Ĥ|Γ is close to m th.;Ĥ0, E0= Γ|Ĥ0|Γ in smallness in λ. Hence, we expect that lim t→∞ Γ(t)|m|Γ(t) is a finite amount away from Γ|m|Γ . Given the bound in Eq. (18) on the time derivative, we have Hence, at least until time t (λ) ∼ λ −1/2 the observable will still be a finite amount away from the thermal value. Some remarks are in order here. First and most importantly, while the bound Eq. (19) on the derivative of the observable is valid for any initial eigenstate, if |Γ were a "thermal" state (i.e., satisfying ETH with respect toĤ 0 ), the above argument would not give us a divergent relaxation time for small λ since in this case Γ|m|Γ = m th.;Ĥ0, E0= Γ|Ĥ0|Γ and is close to lim t→∞ Γ(t)|m|Γ(t) = m th.;Ĥ, E= Γ|Ĥ|Γ . Thus, to obtain the divergent t (λ), it was crucial to use the nonthermal property of the scar states, namely that local observables have expectation values different from the thermal ones at the same energy density. Second, the above argument generalized to d dimensions would replace the bound Eq. (19) on the derivative of the observable by c 0 t d−1 + c 1 t d , which would yield thermalization time at least as long as t (λ) ∼ λ −1/(d+1) .

C. Possibility of stronger bounds on the thermalization time
We now ask if we can obtain stronger bounds on the thermalization time than the above t (λ) ∼ λ −1/2 in d = 1. The bound we obtained is indeed very general and used almost no information aboutĤ 0 andV . The only assumptions we made were the locality ofĤ 0 andV and that |Γ violates ETH.
On the other hand, we suspect that in the problem of the PXP scar states, the thermalization time diverges with a stronger power law, perhaps even as λ −1 in our specific model. Such a suspicion can already be supported from Fig. 7, where we see that | Γ(t)|m|Γ(t) − Γ|m|Γ | grows almost linearly in t instead of the quadratic growth that appeared in the upper bound in Eq. (21). Equivalently, examination of the derivative d Γ(t)|m|Γ(t) /dt of the measured observable in Fig. 7 shows that the derivative remains bounded at least up to the time that the numerical calculation is reliable. Namely, there is no evidence of the linear in time growth that appears in the rigorous upper bound in Eq. (19), which suggests significant over-estimation in the bound.
More generally, the reason for this suspicion is as follows. Consider the step in the argument in Sec. IV B where for |j| ≤ v LR t, we used the bound [v j , e iĤtm e −iĤt ] ≤ 2 v j m . While this is probably the best bound for the operator norm of the commutator, we are actually interested in the expectation value of the commutator evaluated in the initial state |Γ . This expectation value is essentially a retarded Green's function between the local operatorsm andv j , where the time dynamics is determined byĤ and the initial ensemble is given by the pure state |Γ . The retarded Green's function in turn can be related to a dynamical time-ordered correlation function between the operators at space-time points (0, t) and (j, 0) (we take t > 0 throughout), where we have used hermiticity ofm andv j . It is likely to be very crude to use the operator norm to bound this expectation value. Instead, we suspect that the above retarded Green's function decays both in space and in time, and the sum over |j| ≤ v LR t can be bounded by a slower tdependence than t 1 in Eq. (19), perhaps even by a tindependent number in many cases. This suspicion is based on the physical intuition that correlation functions decay at large spatial or temporal separation (note that Γ|e iĤtm e −iĤt |Γ and Γ|v j |Γ are both real numbers, so the above expression equals the imaginary part of the connected correlation function). While for |j| > v LR t this expectation is formalized by the Lieb-Robinson bounds, here we need the regime t > |j|/v LR , i.e., inside the "light cone." If |Γ were a thermal state andĤ a thermalizing Hamiltonian, it would be natural to expect exponential decay of such correlations in time also inside the "light cone." This would completely suppress the t 1 piece in Eq. (19).
As a less favorable example, suppose we know that the above correlation function is bounded by ∼ t −p for all |j| ≤ v LR t. Then in Eq. (19) we would be able to replace the right-hand-side bound by λ(c 0 + c 1 t 1−p ) and in Eq. (21) the right-hand-side bound by λ(c 0 t+c 1 t 2−p /(2− p)). If p ≥ 1, then we could argue that the observable would still be a finite amount away from the thermal value at least until time t (λ) ∼ λ −1 , while if p < 1 we would only be able to argue that t (λ) ∼ λ −1/(2−p) . These conclusions would hold also if the retarded correlation function is bounded inside the light cone by |G R (x, t)| ≤ A/(v 2 LR t 2 − x 2 ) p/2 , which is the form encountered at zero-temperature quantum critical points described by conformal field theories (strictly speaking, we also need p < 2 for the spatial integral over |x| < v LR t to be convergent, where for concreteness we specialized to d = 1).
From the above discussion, we see that if we allow arbi-traryĤ and |Γ , we probably cannot improve the bound in Eq. (19) just on general grounds: the details of the system at hand are important. However, for many systems it is likely that this bound significantly overestimates the actual rate of change of the observable, and the thermalization time will actually be significantly longer than suggested by the rigorous argument. Our direct numerical study suggests that this is the case for the exact PXP scar states under the PXPZ perturbations.

V. CONCLUSIONS
In this paper, we examine the fate of exact quantum many-body scar states under perturbations. In particular, we consider the PXP model with two exact scar states at energies E = ± √ 2, and two exact scars that are in the E = 0 degenerate manifold. We consider the perturbationV in Eq. (6) that preserves the particle-hole property, the inversion symmetry, and hence the E = 0 manifold (but not the exact scars), to investigate the effects of the perturbation on the exact scar states. We also compare the scar states to some thermal states under the perturbation.
In finite sizes, we find robust signatures of deformed scars in the perturbed Hamiltonian. In particular, the smallness of the bipartite entanglement entropy and the presence of the VBS order seem to survive to some degree under the perturbation. The hybridization of the scar states to other states also seems to be weaker compared to the hybridization of the thermal states, which is seen in the loss of fidelity relative to the unperturbed states. Furthermore, we use Rayleigh-Schrodinger perturbation theory to construct perturbed scar states and find good agreement with the ED states.
Nevertheless, by examining the finite-size scaling of the off-diagonal matrix elements connecting the scar and thermal states, we conclude that the aforementioned robustness will be lost in the thermodynamic limit. In particular, we find that the off-diagonal matrix elements between the scar states and the nearby thermal states scale as φ −L/2 , in agreement with the ETH picture of the thermal states. Although these matrix elements have relatively smaller amplitudes compared to the off-diagonal matrix elements between thermal states, the difference appears to be quantitative and not qualitative. Such a scaling of the matrix elements is not enough to beat the exponential decrease of the level-spacing ∼ φ −L , and accordingly we expect the eventual thermalization of the scar states. In addition, while not presented in this work explicitly, we also performed similar finite-size scaling analysis on the embedded Hamiltonian constructed in Ref. [29], and reached the same conclusions. Accordingly, we expect that the eventual thermalization of the exact scar states under perturbations holds more generally.
On the other hand, despite the eventual thermalization, we show that the nonthermal properties of the exact scar states can survive for some parametrically long time even in the thermodynamic limit. Specifically, we show numerically that the VBS order in the exact scar states in the PXP model can survive for a long time after a quench to the perturbed Hamiltonian. We also present a general theory that rigorously lower-bounds the thermalization time as t * ∼ O(λ −1/(1+d) ), where d is the dimension of the system. In practice, depending on the details of the system, the thermalization time scale can even be longer. By comparing the actual time evolution of the derivative of the observable in the numerical simulation with the bounds used in the rigorous argument, we propose that in the specific perturbed PXP model the thermalization time diverges at least as λ −1 .
Our work provides an important foundation for understanding how generic the exact scar states found and constructed in various special models are, and possible relevance of such exact scar states for understanding apparent "scarness" of nearby models where exact scars are not known. We are considering the "worst case" scenario, where an exact scar state in the spectrum is surrounded by chaotic states. In this case, the random matrix theory description of the off-diagonal matrix elements scaling perhaps is unavoidable and leads to eventual thermalization. However, the thermalization time can be large as long as the perturbation is small. This means that in experiments, some nonthermal signatures in dynamics can indeed be understood from some specialĤ 0 and its exact scar states. Furthermore, Ref. [23] proposed an "embedded Hamiltonian" formalism as a way to engineer Hamiltonians which have nonthermal states inside the spectrum. Our general lower bound on the thermalization time suggests that the nonthermal signatures of these states can survive for some long time even when perturbations break the exact embedded structure. This also establishes the possibility of engineering embedded Hamiltonians to protect quantum information.
If the above scenario of thermalization of scar states under generic perturbations is true, we are led to the following speculation about the PXP model itself. Original ED studies found a number of approximate scar states in the model. We speculate that the PXP model may be proximate to some model that has a larger number of exact scars. Under the perturbation that takes this unknown model to the PXP model, some of these become the approximate numerical scars, while others remain exact for special reasons. From our systematic study of the Schmidt numbers of all eigenstates in the PXP model, we conjecture that only the presented exact scar states are truly exact, in the sense that they can be expressed analytically for any system size, while the other scars are no longer exact in the PXP model. In this case, we would predict that these other scar states do not survive in the thermodynamic limit, but their apparent scarness can still persist up to large sizes/long times. It would be interesting to look for such a tractable "mother model" that could explain all scars that are very prominent in the PXP model. Reference [14] discovered four exact scar states in the PXP modelĤ 0 , Eq. (2), with open boundary conditions. We denote the states as |Γ αβ , where α, β ∈ {1, 2}. The wavefunctions for |Γ αβ are most economically expressed as matrix product states.
Defining the boundary vectors v 1 = (1, 1) T and v 2 = (1, −1) T and 2 × 3 and 3 × 2 matrices we can express the scar states as where N αβ = 2 3 L/2 +(−1) L/2+α+β is the normalization factor. The states |Γ 11 and |Γ 22 are not orthogonal but have overlap The normalization factor of |Γ I = 1 √ N (|Γ 11 − |Γ 22 ) is therefore as stated in the main text. These scar states have the following symmetry properties. For the inversion, we have Hence I|Γ I = −(−1) L/2 |Γ I , i.e., |Γ I has the same inversion quantum number as |Γ 21 . For the particle-hole transformation, we have In this appendix, we prove that the perturbation energy for |Γ I ≡ |Γ is an eigenstate of the particle-hole transformation with eigenvalue C = −(−1) L/2 to any order. (We caution, however, that this does not imply that the perturbation theory converges.) First, we briefly review the formal development of the Rayleigh-Schrodinger perturbation expansion. We closely follow the notations and settings in Ref. [56]. In the perturbation theory, we are trying to solve the eigenvalue equation Γ I is the energy shift. Note that we can immediately see Γ We further assume the (formal) series expansion of |Γ I (λ) = (B4) Moreover, from Eq. (B3) and plugging in the series expansion of |Γ I (λ) , equating the same order of λ, we have, at N -th order Therefore, C|Γ order. This indeed also implies ∆ (n) Γ I = 0 for any order n.

Appendix C: Diagonally improved perturbation theory
To further elaborate on the use of the perturbation theory in finite sizes in Sec. III A when there are avoided level-crossings, particularly for the scar state |Γ 21 , we modify the perturbation theory in the following way. We include in the "unperturbed" Hamiltonian the diagonal part ofV , and treat the off-diagonal part ofV as perturbation. More specifically, in the eigenbasis ofĤ 0 , namely |n (0) , the matrix element of the diagonal part of V is [V diag. ] nm = n (0) |V |n (0) δ nm and the off-diagonal part isV off-diag. =V −V diag. . The unperturbed Hamiltonian is nowĤ 0 =Ĥ 0 + λV diag. and the perturbation is λV =V off-diag. . Clearly, the unperturbed eigenstates are still |n (0) .
Based on this regrouping, we can use the standard perturbation theory Eqs. (7) and (9) with the replacements In Figs. 9(a) and (b), we show the comparison between the diagonally improved perturbation theory and the ED results for the state |Γ 21 , similar to Figs. 3 and 4. Note that on the λ < 0 side, the perturbation theory also shows a drop of the overlap around λ ≈ −0.02, though the recovery at λ ≈ −0.04 is spurious. On the λ > 0 side, the diagonally improved perturbation theory also shows the accidental hybridization phenomenon. It happens when E We can now appreciate the improvement provided by this perturbation theory but also its eventual failure. For an isolated level crossing driven by increasing λ and the difference between λ n (0) |V |n (0) and λ Γ 21 |V |Γ 21 , this perturbation theory is accurate both well before the crossing and well after the crossing, which is why we see improvement in the 2nd-order estimate of the energy in Fig. 9(a) compared to Fig. 3, particularly on the λ > 0 side. However, once the separation between successive energy level crossings (set by the density of states) becomes smaller than the duration of the avoided level crossings (set by the off-diagonal matrix elements), the improved perturbation theory also fails. Upon increasing the system size, the density of states and hence the frequency of level crossings increases very fast and cannot be compensated by the decrease of the typical off-diagonal matrix elements, as discussed in Sec. III B. Hence the failure of the perturbation theory and the eventual thermalization of the scar states appear inevitable in the thermodynamic limit. In the perturbed PXP model here, for the perturbation strength λ 0.05, our largest (b) Comparison of the fidelity obtained in the diagonally improved perturbation theory and the ED result. The perturbation theory now also accounts for the phenomenon of accidental hybridization.
ED system sizes are just reaching the regime where this starts to happen. 10. System size scaling of the variance ofV in eigenstates |Γ21 , |ΓI , and |Γ th . Since the thermal state |Γ th is inherently random, its system size dependence is "noisy." First, we note that where all states are assumed normalized (here and below, we use the same notation as in Sec. III A). This allows us to think about the squared amplitudes of the matrix elements between the |Γ and other states |n (0) as some "weights," where the total weight is given by the variance of the perturbationV in the state |Γ . SinceV is a local Hamiltonian,V = jv j , the variance ofV grows linearly with the system size L: Here we have assumed translational invariance and that connected correlation functions of local observables are short-range, which is true for our exact scar states with the finite bond dimension. We expect this to be true also for short-range-correlated thermal states (e.g., away from any finite-temperature critical point or critical phase), in particular for the thermal states used in the main text. Figure 10 shows measured variance ofV in the scar and thermal states used in the main text as a function of L and confirms the expected linear in size scaling. The plot for the thermal state has noisy L dependence, which is expected since the thermal states are inherently "random," so even our "deterministic" procedure of picking the thermal state whose index is given by the index of the exact scar state plus three, when going from one size to the next has effective randomness in it. On the other hand, the plots for the scar states use the same MPS but for different systems sizes, and there is no randomness in this.
We next observe that the mixing weights | n (0) |V |Γ | 2 are significant only for states with |E The reason is again locality ofV : Eachv j when acting on |Γ can "add" or "remove" only O(1) energy as measured byĤ 0 , sov j |Γ expanded over the eigenstates ofĤ 0 is concentrated on |n (0) with |E . This observation can be formalized as follows. Consider n =Γ HereŴ is an operator which is a sum of local termsŵ j (the factor of i = √ −1 makesŴ hermitian to simplify expressions). Hence, similarly to Eq. (D2), the variance ofŴ in the state Γ grows linearly with the system size, We then conclude that for large L The left hand side can be interpreted as an average of over the distribution of n's with weights proportional to | n (0) |V |Γ | 2 , and we see that this average is an L-independent number, as claimed. Figure 5 in the main text indeed shows that the matrix elements | n (0) |V |Γ | become very small once |E (0) n −E (0) Γ | exceeds some characteristic energy scale. This is true for both the scar and thermal states |Γ studied.
1. Simple-minded application to the time-independent 2nd-order perturbation theory and the large L limit As an application, consider the formal 2nd-order perturbation theory expression for the ratio of the probability of being in any state other than Γ to the probability of remaining in the state Γ [cf. Eq. (9) in the main text]: We can bound this as follows: Here we used Eq. (D1) and introduced some fixed energy scale ∆, which will be chosen below. We can now use Eq. (D3) to bound the second term in the square brackets as follows: var(Ŵ ; Γ) ≥ n, |E We finally obtain where in the very last equation we picked ∆ to make this lower bound as large as possible. Thus, we see that the probability for the perturbed eigenstate to remain in the initial state |Γ decreases to zero with L, at least in this formal treatment. Of course, it is known that the fidelity of a many-body state under a generic perturbation goes to zero in the thermodynamic limit, and the above is not intended as any serious proof but only as a simple illustration of thinking about the distribution of the matrix elements. Note that the bound does not use any information about this distribution except the variance; in particular, it applies also, e.g., for a gapped ground state that is separated from the rest of the states by a gap-the fidelity under perturbation still goes to zero. In the case of states at finite energy density that are surrounded by many thermal states, for very large L the above lower bound is actually a gross underestimate of how poorly the formal static second-order perturbation theory performs: With the level spacing decreasing as ∼ D −1 L (where D L is the dimension of the total Hilbert space which grows exponentially with L), and the matrix elements decreasing as ∼ D 2. Application to the time-dependent 2nd-order perturbation theory for fidelity after a quench As another application, consider the following quantity: which arises when we apply the time-dependent 2ndorder perturbation theory to the quench setting described in Sec. IV. Specifically, the system starts in the state |Γ at time t = 0 and evolves under the perturbed Hamil-tonianĤ =Ĥ 0 + λV . The above P (t) is the perturbative result for the probability of being in any state other than Γ at time t; hence, the probability of remaining in the state Γ, or fidelity at time t after the quench, is 1 − P (t). Clearly, P (t) also lower-bounds the quantity R in Eq. (D7) that arises in the time-independent 2nd-order perturbation theory. However, the dynamical quench setting is more interesting in that the stated results are potentially more accurate and useful on reasonable time or length scales because of the effective control over the denominators provided by the sin 2 (E (0) n −E For very small t such that for all significantlyparticipating n the quantity (E   (1), so we actually expect the above formula to be an accurate description of the formal perturbation theory result up to t ∼ O(1). In particular, this formula is valid for t ∼ 1/ √ λ 2 αL, beyond which the perturbation theory would give the probability P (t) exceeding unity. Beyond this time, we clearly cannot apply such a formulation of the perturbation theory, but until a somewhat smaller time of the same order such an application actually appears sensible.
More precisely, let us pick a number y ∈ (0, π) and set f y = sin 2 (y)/y 2 . We have the following bounds: where in the third line we used Eqs. (D1) and (D8), while in the last line we used Eqs. (D2) and (D5). For simplicity, let us keep O(1) number y and hence f y fixed, i.e., we do not try to further optimize this degree of freedom. We see that, e.g., for t < y α/β and t > 1/ 3λ 2 f y Lα/4, which can be satisfied simultaneously for large enough L, we have P (t) > 1. Hence, the 2nd-order perturbation theory already fails beyond a time that scales as L −1/2 , as claimed earlier.
To summarize, the above analysis suggests that the time scale for the fidelity loss goes to zero for large L but only as a power law L −1/2 . We remark that the above arguments are true for both scar and thermal states, and the difference in such short-time fidelity loss is only quantitative. Nevertheless, the above analysis helps, e.g., when we want to understand ED results for the fidelity loss in the quench setting of Sec. IV. We did not present such fidelity results as they do not allow defining slow thermalization in the thermodynamic limit. Instead, as presented in the main text, measuring local observables shows that the nonthermal signatures of the scar in the initial state persists to nonzero time even when L → ∞, and this thermodynamic-limit time diverges when the perturbation strength goes to zero.