{'fraction_non_alphanumeric': 0.038570618264322176, 'fraction_numerical': 0.012179365979706308, 'mean_word_length': 4.388318193395025, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Driving quantum many-body scars

Periodic driving has been established as a powerful technique for engineering novel phases of matter and intrinsically out-of-equilibrium phenomena such as time crystals. Recent work by Bluvstein et al. [Science 371 , 1355 (2021)] has demonstrated that periodic driving can also lead to a signiﬁcant enhancement of quantum many-body scarring, whereby certain non-integrable systems can display persistent quantum revivals from special initial states. Nevertheless, the mechanisms behind driving-induced scar enhancement remain poorly understood. Here we report a detailed study of the eﬀect of periodic driving on the PXP model describing Rydberg atoms in the presence of a strong Rydberg blockade – the canonical static model of quantum many-body scarring. We show that periodic modulation of the chemical potential gives rise to a rich phase diagram, with at least two distinct types of scarring regimes that we distinguish by examining their Floquet spectra. We formulate a toy model, based on a sequence of square pulses, that accurately captures the details of the scarred dynamics and allows for analytical treatment in the large-amplitude and high-frequency driving regimes. Finally, we point out that driving with a spatially inhomogeneous chemical potential allows to stabilize quantum revivals from arbitrary initial states in the PXP model, via a mechanism similar to prethermalization.

Recent experiments 10 on QMBS utilising Rydberg atom arrays in various geometries, both one-dimensional and two-dimensional, have shown that periodic driving can stabilize the revivals which occur in the static system [49][50][51] . Similar effect was also observed in experiments on a Bose-Hubbard quantum simulator 33 , where QMBS phenomenon was identified in a larger class of initial states. In both cases, the optimal driving frequency was empirically found to be related to the frequency of revivals in the static model. Although the driving protocol is relatively simple, the mechanism behind it remains to be understood. Previous work in Ref. 52 provided an important insight that the experimental driving protocol can be approximated by a simpler two-step (kicked) driving scheme. Such a toy model is easier to treat analytically and it captures some of the observed phenomenology. However, its relation to the experimental driving protocol is not transparent, in particular there is no direct mapping between the parameters of the two protocols (except for the driving frequency, which is equal in the two cases).
Moreover, further questions naturally arise. For example, the only two initial states whose revivals were stabilized by periodic driving (the Néel and polarized states) are precisely the ones that exhibit QMBS in the static version of the same model. This suggests that there are fundamental limitations to the driving, which only "works" for such special states. If the enhancement mechanism is indeed related to QMBS, it would imply that it should be possible to improve the revivals in a variety of other QMBS models that have been studied theoretically. Another question is how to predict the optimal driving parameters without performing a computationally-costly brute force search. To this end, Ref. 10 has identified a window of driving frequencies centered around twice the frequency of revivals in the static model. A subharmonic response to driving was observed in this frequency regime, suggesting a relation to the discrete time crystal, which also warrants further investigation.
In this work we numerically and analytically study the effects of driving on the PXP model describing Rydberg atoms in the regime of strong Rydberg blockade. In Sec. II we introduce the model and investigate in detail the spatially-uniform cosine driving scheme that was used in recent experiments. We identify two regimes with distinct properties (dubbed 'Type-1' and 'Type-2'). In Sec. III we focus on the optimally driven Néel and po-larized state. In Sec. IV we introduce a toy model based on the square pulse driving protocol. We show that this protocol better approximates the cosine-drive dynamics than the previously studied kicked drive 52 . The squarepulse protocol is amenable to analytical treatment in certain limits, as we demonstrate for the high-amplitude and high-frequency drive regimes. Finally, in Section V we discuss the implications of our results on the relation between QMBS and the enhancement by periodic driving, suggesting that QMBS is a necessary ingredient in the Type-1 (typically low-amplitude) regime, while Type-2 driving (typically high-amplitude) can be also applied to non-scarred models. We summarize our findings and identify some future directions in Sec. VI. Appendices A-D provide additional information on system size scaling, a modified (spatially inhomogeneous) cosine driving protocol which can generate revivals from any initial product state, high frequency expansion for the Floquet Hamiltonian, and an alternative way to visualize the trajectory of a state in Hilbert space.

II. PXP MODEL AND COSINE DRIVE
The PXP model 53,54 describes a kinetically constrained chain of spin-1/2 degrees of freedom. The chain serves as a model for the system of Rydberg atoms, each of which can exist in two possible states, |• and |• , corresponding to the ground state and the excited state, respectively. An array of N such atoms is governed by the Hamiltonian where X j = |• j • j |+|• j • j | is the Pauli x-matrix on site j, describing local Rabi precession with frequency Ω. The projectors onto the ground state, P j = |• j • j |, constrain the dynamics by allowing an atom to flip its state only if both of its neighbors are in the ground state. Unless specified otherwise, we use periodic boundary conditions (PBC) and set the PXP term strength to Ω = 1. Before introducing the driving scheme, we briefly review some pertinent properties of the static PXP model. The PXP model is quantum chaotic 49 . However, it was found that preparing the system in the Néel or |Z 2 ≡ |• • • • ... • • state, which has an excitation on every other site, leads to persistent quantum revivals 50,55,56 . Analogously, the translated Néel state, |Z 2 ≡ |• • • • ... • • , also displays revivals. In a conventional thermalizing system, the revivals are not expected as the |Z 2 initial state has predominant support on the energy eigenstates in the middle of the many-body spectrum, thus it corresponds to "infinite" temperature. The presence of revivals in a special initial state in an overall chaotic system was understood to be a manybody analog of the phenomena associated with a single particle inside a stadium billiard, where nonergodicity arises as a "scar" imprinted by a particle's classical periodic orbit 51,57,58 . In QMBS systems, the eigenstates were shown to form tower structures 49 , resulting from the anomalously high overlap of eigenstates with the |Z 2 initial state, and the equal energy spacing between the towers leads to quantum revivals. We mention that the scarring in the PXP model was also shown to persist in higher dimensions 10,59,60 and in the presence of certain perturbations 50,61,62 , including disorder 63 .
The revivals in the PXP model are not perfect, in particular the wave function fidelity, with the initial state chosen to be |ψ 0 = |Z 2 , visibly decays over moderate times 50 . Reference 10 showed that the decay time-scale can be significantly extended by introducing a cosine modulation ∆(t) of the chemical potential, where n i = |• i • i | is the density of excitations on site i, ∆ 0 is the static detuning, ∆ m modulation amplitude and ω the driving frequency. The periodic modulation of the chemical potential in Eq.
(3) was shown to stabilize and enhance the revivals in the quantum fidelity, both in numerical simulations and in experiments on quantum simulators 33 . The enhancement was achieved for two different initial states, the previously mentioned |Z 2 state, as well as the polarized state which has no excitations, |0 ≡ |••• . . . ••• . The optimal driving frequency for the Néel state was found to be close to twice the frequency of revivals in the pure PXP model (∆(t) = 0); for the polarized state, the optimal frequency was instead close to the frequency of revivals in the detuned PXP model (∆(t) = ∆ 0 ). The optimal values of ∆ 0 and ∆ m were found empirically and have different values for these two states.
In order to determine more precisely the optimal driving parameters, we scan the parameter space (∆ 0 , ∆ m , ω) for the highest average fidelity between the times t = 0 and t = τ , where we set τ = 100. We restrict to the case ∆ m ≥ 0 since F t (∆ 0 , ∆ m , ω) = F t (−∆ 0 , −∆ m , ω). In Figs. 1(a) and (b) we plot the average fidelity as a function of ∆ 0 and ∆ m at the fixed frequency ω, for the Néel and polarized initial state. The frequency is fixed to a value close to the revival frequency of the static model, i.e., we choose ω = 2.72 in (a) and ω = 3.71 in (b). These two figures have many similarities, such as parallel diagonal lines consisting of "butterfly-shaped" peaks. However, a closer look reveals an important difference. The Néel state plot contains several additional peaks which are very narrow and typically located between the bright  diagonal lines. For example, one of these peaks is approximately around the point (∆ 0 ,∆ m )=(1.15,2.67) and corresponds to an interesting optimal parameter regime for the Néel state which will be studied below. We will call these narrow peaks "Type-1" peaks. The other bright regions present in both Fig. 1(a) and (b) will be dubbed "Type-2 peaks" (excluding |∆ 0 | 1 and ∆ m |∆ 0 |). These parameter regions are broader, i.e., the parameters do not have to be finely tuned, and they typically occur at higher driving amplitudes than Type-1 peaks. Like Type-1, their position depends on the driving frequency, but the mechanism of revival stabilization is different. We will discuss Type-2 peaks below in relation to the polarized initial state.
It is important to note that the Néel and polarized state are two special cases. The cross sections in Figs. 1(a) and (b) look very different for other initial product states, with low and very narrow peaks of the cost function instead of Type-2 peaks, while the Type-1 peaks are completely absent.
Finally, Fig. 1(c) shows three illustrative cross sections taken at different frequencies, ω ∈ {1, 2, 3}. The initial state is the Néel state and the color scale corresponds to the value of the average fidelity cost function. We see that this cost function has many peaks (bright regions), which move and deform as ω is varied. The position of the fidelity peaks in (∆ 0 , ∆ m ) space depends approximately linearly on the driving frequency. In all three cases, there are two large bright regions in the |∆ 0 | 1 and ∆ m |∆ 0 | regime. When the static detuning ∆ 0 is large, the energy spectrum splits into disconnected bands. This means that the wave function cannot spread through the entire Hilbert space and the fidelity remains relatively high during the evolution, never dropping down to zero. These trivial regimes do not correspond to typical QMBS behavior, hence we exclude those where the minimal fidelity from t = 0 to t = 100 is non-zero. The optimal parameter regimes that will be discussed in detail below are summarized in Table I. Next we determine the optimal window of driving frequencies. In order to do this, we choose one of the narrow peaks of the cost function in Fig. 1 and follow it as we vary the frequency ω. The driving parameters ∆ 0 and ∆ 0 are reoptimized for each frequency. When we leave the optimal frequency window, the peak completely disappears, but other peaks may become more prominent, as summarized in Fig. 2. In order to characterize the quality of revivals, in Fig. 2 we used the following cost function where T is the driving period and M denotes the number of revivals averaged over. The integer l defines the order of the response, e.g., l = 1 for the harmonic and l = 2 for the subharmonic responses we will study below. In both cases we average over lM = 40 driving periods. The value of C l corresponds to the difference between average maximum fidelity and average minimum fidelity over the first M revivals, with the value of 1 corresponding to perfect revivals. In Fig. 2 we use C l instead of the previously introduced average fidelity because it gives better results when comparing different frequencies. Namely, as the driving frequency is increased, the revivals become more dense and the width of an individual revival increases in comparison with the revival period, which then results in increased average fidelity. Therefore, decaying revivals at higher frequencies can lead to higher average fidelity than stabilized and constantly high revivals at lower frequencies. The new cost function does not suffer from this problem and can accurately estimate the quality of revivals independently of the frequency. However, we still use the first cost function for fixed frequency scans, such as those in Fig. 1, since the average fidelity can capture all types of revivals (e.g., harmonic, subharmonic or even incommensurate with the driving frequency), without the need to specify their expected frequency. In Fig. 2(a) we plot the data for the peak which passes through a set of Néel state optimal parameters (∆ 0 = 1.15, ∆ m = 2.67, ω = 2.72). We look at M = 20 revivals in evaluating the cost function C 2 in Eq. (5). The dependence of ∆ 0 and ∆ m on ω is seen to be approximately linear. In Fig. 2(a) we see that there is indeed a window of optimal frequencies, approximately 2.5 < ω < 3.5, where we can obtain a good subharmonic response to periodic driving. This window is centered  at a frequency slightly higher than twice the bare revival frequency of the Néel state in the pure PXP model (ω 0 ≈ 2.66). In Appendix A, we study the system size scaling of this optimal driving frequency, see Fig. 13. For N ≥ 10, the optimal frequency window for the Néel state barely changes with increasing system size. Note that there is a minimal frequency in Fig. 2(a) bellow which it is impossible to obtain revivals. On the other hand, it is possible to enhance the revivals in the high frequency regime, but the lifetime of these revivals is finite and decreases with increasing frequency. We can repeat the same procedure for other Type-1 peaks, see Fig. 2(b) for one example. In this case we start from the optimal parameters which correspond to the "one-loop" trajectory that will be studied in more detail in Sec. V B (∆ 0 = −1.17, ∆ m = 4.81, ω = 2.00). The dependence of optimal ∆ 0 and ∆ m on ω is again linear, but the optimal frequency window is now more narrow and centered around ω = 2, which significantly deviates from ω 0 ≈ 2.66. The results for other Type-1 peaks are found to be similar (data not shown).
Finally, this procedure can also be adapted to the harmonic response, e.g. for the polarized and other initial states with the cost function C 1 . The results for the polarized state and M = 40 in the low-and high-amplitude regime (parameters from Table I) are shown in Fig. 2(c) and (d), respectively. In both regimes, the optimal frequency window is much wider than for the Néel state.

III. NÉEL AND POLARIZED STATE REVIVALS
In this Section we investigate in more detail the optimally driven Néel and polarized initial states.

A. Néel state
Using the procedure described in Sec. II -scanning the parameter space for the highest average fidelitywe found the following optimal parameters for the Néel state: ∆ 0 = 1.15, ∆ m = 2.67, ω = 2.72. This is the lowest-amplitude Type-1 peak in Fig. 1(a). In addition to brute force scanning, we also employed the simulated annealing optimization algorithm, which produced similar optimal values. We note that our driving parameters slightly differ from those considered in Ref. 10, ∆ 0 = 1.00 and ∆ m = 2.00, which were obtained at the fixed frequency ω = 2.66. Fig. 3(a) shows the evolution of quantum fidelity for the optimal parameters quoted above. In contrast to the undriven case, the revivals are clearly enhanced and remain stable over long times. Note that the revival frequency is half the driving frequency. This is an example of period doubling (subharmonic response) which was related to a discrete time crystal in Ref. 10. The growth of the bipartite von Neumann entanglement entropy S A is strongly suppressed by driving, as shown in Fig. 3(b). To evaluate S A , we partition the chain into two equal halves and evaluate the reduced density matrix for one of the subsystems, ρ A (t) = tr B |ψ(t) ψ(t)|, by tracing out the other (B) subsystem. This yields S A = −tr ρ A ln ρ A .
Furthermore, we numerically construct the one period evolution operator of the optimally driven system and diagonalize it in order to obtain the corresponding Floquet modes. The bipartite entanglement entropy versus the quasienergy of all the Floquet modes is plotted in Fig. 3(c). The color represents the expectation value of the number of excitations i n i for each mode. This entropy plot has an interesting "two arc" structure. The origin of these two arcs will be explained in Sec. III C below. The two lowest-entropy modes (marked by the red circles) also have the highest overlap with the Néel state. Two special Floquet modes were also observed in Ref. 52, where the cosine drive was approximated by a delta function pulsed drive. We find that the two modes are approximately (|Z 2 +|Z 2 )/ √ 2 and (|Z 2 −|Z 2 )/ √ 2, which explains the observed period doubling 33 . However, the subharmonic response becomes harmonic if we work in the momentum k = 0 sector and choose a translation invariant initial state (|Z 2 + |Z 2 )/ √ 2. In that case, there will still be two arcs in the entropy plot, but only one high-overlap Floquet mode.
In addition to (∆ 0 , ∆ m ) = (1.15, 2.67) parameters,  there are several other Type-1 peaks in the cost-function scan at fixed frequency ω = 2.72 that we have seen in Fig. 1(a). The entanglement entropy plots of their Floquet modes also display the distinct "two arcs" structure. These other Type-1 peaks also lead to enhanced revivals, although with a finite lifetime at this particular driving frequency. However, some of those peaks are higher at other frequencies and can lead to completely stabilized revivals, see for example Fig. 2(b) and the discussion in Sec. V B. The number, position and height of Type-1 peaks depends on the driving frequency, while the driving amplitude of peaks which result in stable revivals is typically lower than that of Type-2 peaks.

B. Polarized state
For the polarized initial state, particularly in the lowamplitude driving regime, the optimal frequency value is found to be more sensitive to the cost function ( F t or C 1 ) used for the optimization. In particular, the C 1 cost function results in a larger value of the driving amplitude, see Fig. 2(c). Below we will use a set of parameters close to the low-amplitude edge of the optimal frequency window. These parameters were also used to drive the polarized state in a recent experiment on the Bose-Hubbard quantum simulator of the PXP model 33 .
Although the fully polarized state |0 shows no QMBS behavior and quickly thermalizes in the pure PXP model, it was shown that this can be changed by applying a static detuning term and then further enhanced by peri-odic driving 33 . Figure 4(a) shows the evolution of quantum fidelity for the polarized state in the pure PXP model (black), with static detuning (red), and with periodic driving (blue). Similar to the Néel state, the growth of entanglement entropy is strongly suppressed by the static detuning and even more so with periodic driving, see Fig. 4(b). The entanglement entropies and the numbers of excitations for all Floquet modes of the optimally driven system (∆ 0 = 1.68, ∆ m = −0.50, ω = 3.71) can be observed in Fig. 4(c). This driving frequency approximately matches that of revivals in a static system with detuning ∆ 0 = 1.68 (red line in Fig. 4(a)), where QMBS for the polarized state has been previously observed 33 . However, unlike the Néel state, there is now only one arc and a single Floquet mode which has high overlap with the polarized state (≈ 64%, encircled in red).
Apart from this set of driving parameters, which were studied in Ref. 33, there are also other sets of driving parameters that lead to robust revivals from the polarized state. One such example is shown in Fig. 4 55, ω = 2.90). Upon first glance, this situation is completely different from Fig. 4(a)-(c). In the latter case, the revivals essentially stem from static detuning, which creates scarring, and the driving further stabilizes them. By contrast, in the former case, the revivals are generated by periodic driving. The optimal driving frequency was found using simulated annealing and it does not appear to correspond to any special frequency in the static system. Moreover, the Floquet mode entanglement entropy plot also has a different structure, as can be observed in Fig. 4(e). There are now multiple arcs (N/2 + 1 of them), each with a well defined excitation number. This suggests that the Floquet Hamiltonian is fragmented in this parameter regime (see further discussion in Sec. III C below). There is a single Floquet mode which has near-unity overlap (> 90%) with the polarized state, which accounts for the revivals. The mode belongs to its own one-dimensional excitation number sector, since the polarized state is the only state that contains no excitations. The mode is furthermore separated from the neighboring i n i = 1 sector by a gap ∆ , which practically does not depend on the system size -see Fig. 12 in Appendix A.
Based on these insights, in Appendix B we present a modified, spatially inhomogeneous cosine driving scheme, which can enhance the revivals from any initial product state, including generic states that do not display QMBS. The main idea is to change the driving protocol in such a way that the desired initial state becomes the ground state of the new time-dependent term. Driving with optimal parameters then results in similar phenomenology to the polarized state in the high-amplitude regime ( Fig. 4(d)-(f)), including the fragmentation of the Floquet Hamiltonian into sectors with different quantum numbers.
We note that both of the driving regimes presented in Fig. 4 belong to what we referred to as Type-2 peaks. In these cases, the optimal parameters do not require particular fine-tuning, i.e., there is a fairly broad range of parameter values for which the driving stabilizes QMBS revivals, unlike the Type-1 peaks where such regions are more narrow and finding optimal parameters requires more careful optimization. It was not possible to produce a subharmonic response for the polarized state, which is expected since it does not have a translated partner state, such as Z 2 and Z 2 . Since the two Néel states also belong to their own excitation number sector in Fig. 4(f), it is also possible to use the ∆ 0 = 0.64, ∆ m = 7.55, ω = 2.90 parameters to stabilize the revivals from these states, as will be explained in Sec. V.

C. Origin of "multiple arcs" structure
In order to explain the structure of entanglement entropy plots in Figs. 3 and 4, we treat the PXP term in Eq. (3) as a perturbation. First we consider the Ω = 0 case when the Hamiltonian is diagonal but still timeperiodic. The Floquet modes are then product states and their quasienergies depend only on the value of static detuning ∆ 0 and their number of excitations i n i , while the quasienergy spectrum has periodicity ω: As the number of excitations increases, the Floquet modes are winding around the spectrum. Different ratios of ∆ 0 /ω will result in different structures. When the PXP term is turned on, these states start to hybridize. However, only those states that are close in excitation number can hybridize, since the PXP term can only change this number by ±1. They also need to be close in quasienergy, but what is "close enough" depends on the magnitude of the PXP term in the Floquet Hamiltonian. Since the cosine drive is difficult to treat analytically, especially for the Néel state optimal parameters which are in the non-perturbative regime, in the following Section IV we introduce a more tractable square pulse approximation of the cosine driving protocol.
In the parameter regime optimal for the Néel state (Type-1), the Ω = 0 Floquet modes form two diagonal lines (even and odd excitation numbers). Only the states belonging to the same diagonal will mix when Ω > 0, which will result in two separate structures at Ω = 1. This is the origin of the two arcs in Fig. 3(c). Note that the Floquet modes in the two arcs no longer consist of the states with only even or only odd number of excitations. No symmetry which distinguishes between the two arcs could be identified.
The situation is different for the polarized state highamplitude optimal parameters (Type-2). At Ω = 0, the Floquet modes which are close in quasienergy are far apart in excitation number, and the magnitude of the PXP term in the Floquet Hamiltonian is small for these parameters, so there is no mixing between the sectors even when Ω = 1. This results in N/2 + 1 separate sectors as shown in Fig. 4(f), each containing only the modes with a certain number of excitations. The low-amplitude regime from Fig. 4(c) is a special case where the smaller arcs merge into one large arc, while the Floquet modes are still arranged according to their numbers of excitations and the i n i = 0 sector is sufficiently separated from the rest.

IV. SQUARE PULSE DRIVE
The cosine drive in Eq.
(3) is difficult to treat analytically. Therefore, it is desirable to find a simpler model which accurately captures the dynamics, at least at the qualitative level. Previous works 10,52 have considered the following delta function pulse driving model: where T = 2π/ω is the driving period and θ is a parameter of the drive. This driving protocol results in perfect revivals from any initial state when θ = π due to a special symmetry property of the Floquet operator in that case 52 . However, the Néel state is seen to be special when θ is varied away from π: while the revival from other initial states sharply decays, the Néel revival stays robust over a finite range of θ in the vicinity of π. Although this model has been well understood, its relation to the cosine driving scheme is not completely transparent. (b) Red: Cosine drive fidelity, F (t) = | ψ(0)|ψ(t) | 2 . Black (dashed): Delta function pulse fidelity, Néel state, L = 20, optimal driving parameters. (c) Same for square pulse instead of delta function pulse. (d) System size scaling of co-moving fidelity averaged over the first driving period for the square (red) and delta function pulse (blue). Optimal driving parameters: Néel state (squares), polarized state low-amplitude (circles), polarized state high-amplitude (triangles).
We will now show that a simple square pulse model is a better approximation for the cosine drive. The square pulse protocol is defined by as illustrated in Fig. 5(a). We note that similar driving schemes were previously considered in Refs. 47 and 48, however our square pulse model differs from those by a quarter of driving period time shift, which results in very different Floquet dynamics. A comparison of the delta pulse in Eq. (7) and the square pulse in Eq. (8) is shown in Fig. 5(b) and (c), where we plot the quantum fidelity F (t) for the exact evolution under cosine drive, as well as the fidelity for the system driven by the delta or square pulse. Moreover, we directly evaluate and plot the co-moving fidelity, F c (t) = | ψ(t)|ψ (t) | 2 , which represents the overlap between the wave function evolved using the cosine driving protocol, |ψ(t) , and the one evolved using the approximate protocol, |ψ (t) (delta or square pulse). Here, the initial state is the Néel state and the driving parameters are set to their optimal values (∆ 0 = 1.14, ∆ m = 2.08, ω = 2.72 for the square pulse). We note that the optimal driving frequency and static detuning are approximately the same for both the cosine and square pulse drive, but the driving amplitude has to be reoptimized and is approximately 20% lower in the latter case. For the delta function pulse, we use θ = π and the same frequency as the cosine drive.
The co-moving fidelity is seen to be much higher for the square pulse drive, compare Figs. 5(b) and (c). This signals that the square pulse driving protocol is a much better approximation for the cosine drive than the delta function pulse. In particular, while the delta pulse is highly accurate near the revival points, the square-pulse remains accurate throughout the evolution, including away from the revival points. This difference becomes more pronounced in larger system sizes, as can be observed in Fig. 5(d), where we plot the scaling of the co-moving fidelity, averaged over the first driving period with the system size N . In general, the agreement between the cosine drive and square pulse is better when the driving frequency is high and amplitude low. For example, the agreement is particularly good for the polarized state in the low-amplitude regime, where the co-moving fidelity is almost equal to 1 over many driving periods. In addition to accurately capturing the dynamics of the cosine drive, the square pulse drive also reproduces the two special Floquet modes and a "two arc" entanglement entropy spectrum resembling the one in Fig. 3(c). Finally, knowing the optimal driving parameters for the square pulse allows us to estimate the parameters for the cosine drive by simply using the same driving frequency and static detuning while increasing the amplitude by 25%.
The optimal driving parameters identified for the Néel state (∆ 0 = 1.15, ∆ m = 2.67, ω = 2.72) and for the polarized state (∆ 0 = 1.68, ∆ m = −0.50, ω = 3.71) are all of the same order of magnitude as Ω = 1, the characteristic energy scale of the static Hamiltonian. Thus, the stabilization of revivals in these cases due to periodic driving appears to be a non-perturbative effect. Nevertheless, for the square pulse protocol, we can obtain useful analytical insights into the mechanism of the revival stabilization in two extreme limits: the limit of high amplitudes and the limit of high driving frequency. In the remainder of this section, we discuss separately these two cases.

A. High-amplitude regime
When the driving amplitude ∆ m is high, the square pulse allows us to describe the dynamics perturbatively. At first order in the small parameter λ ± ≡ Ω/(∆ m ± ∆ 0 ) 1, it straightforward to obtain the Floquet Hamiltonian using the Pauli algebra. We start from the Hamiltonian whereX j ≡ P j−1 X j P j+1 denotes the Pauli x-matrix dressed by projectors, and ∆(t) is given by the square pulse in Eq. (8). We ignore the energy shift −∆(t)/2 and assume ∆ m ± ∆ 0 Ω > 0. The evolution operator from t = 0 to t = T /4 and from t = 3T /4 to T is denoted U + and from t = T /4 to t = 3T /4 by U − . At first order λ ± 1, we can factorize these evolution operators since the commutator between λX j +Z j and λX j +Z j on different sites j and j is second order in λ. We can then expand the exponentials using the properties of Pauli operators (Z 2 = I as well asX 2 = I) in terms of sine and cosine functions, and evaluate the total evolution operator for one period U = U + U − U + . At first order in λ, the final expression can be easily converted back into exponential form, from which the first-order Floquet Hamiltonian can be directly read off. While the algebra is tedious, the calculation is straightforward and results in the following first-order Floquet Hamiltonian The first-order Floquet Hamiltonian contains a PXP term whose magnitude is controlled by the prefactor g, which depends on the values of driving parameters ∆ 0 , ∆ m and the drive period T = 2π/ω. The g prefactor explains several features of the fidelity phase diagram, as can be seen in Fig. 6. When g is small, the Floquet Hamiltonian is almost diagonal and the Floquet modes are close to product states, which have well defined numbers of excitations. As long as g is small enough (bright regions in Fig. 6(a)), there is no mixing between different excitation-number sectors. This leads to fragmentation and the polarized state ends up in its own one-dimensional sector, while the two Néel states belong to another two-dimensional sector. These states are therefore close to a Floquet mode and keep reviving in the driven system, which leads to high average fidelity (bright regions in Fig. 6(b)). This explains the origin of Type-2 peaks. Conversely, when g is large (for example around the diagonals ∆ 0 = ±∆ m ) the sectors are mixed and the Floquet modes have large variance of the excitation number, meaning that they are superpositions of many product states with different excitation numbers. We do not expect revivals in these regions and Fig. 6(b) confirms that this is mostly the case. However, the optimal driving parameters for the Néel state (Type-1 peaks) are actually in the dark region along the diagonal.
Some features of the phase diagram in Fig. 6(b) cannot be explained only by Eq. (12). Irrespective of the value of g, there are dark lines (no revivals) around ∆ 0 = nω/2, where n is an integer (this is more visible for the cosine drive in Fig. 1). The explanation is related to the quasienergies of the Floquet modes when the PXP term is turned off (Ω = 0), see the discussion in subsection III C. In this case, the Hamiltonian is diagonal and the Floquet modes are product states. Their quasienergies depend only on the value of static detuning ∆ 0 and the Floquet spectrum has periodicity ω, see Eq. (6). When ∆ 0 = nω, all Ω = 0 modes have the same quasienergy = 0. For Ω = 1, any finite value of g will lead to mixing between all sectors. Consequently, no single mode will be close to a product state and there will be no revivals. This underlines the importance of non-zero static detun-ing. The situation is similar when ∆ 0 = (2n + 1)ω/2, but in this case the Ω = 0 Floquet modes split into two groups with same quasienergies (even and odd numbers of excitations). Finite g then causes mixing within these two groups, leading to the appearance of two arcs in the entanglement entropy plots (similar to Fig. 3(b)).
In conclusion, the square pulse drive provides analytical understanding of a large part of phase diagram corresponding to high driving amplitudes. However, we note that the agreement between the square pulse and cosine drive is not as good in the high amplitude regime for the polarized state with parameters ∆ 0 = 0.64, ∆ m = 7.55, ω = 2.90, recall Fig. 5(d). Nevertheless, the square pulse still accurately captures the revivals and the co-moving fidelity is higher than for the delta function pulse protocol.

B. High-frequency regime
Although the optimal set of parameters for the Néel state is deeply in the non-perturbative regime with Ω = 1, ∆ 0 = 1.15, ∆ m = 2.67 and ω = 2.72 of the same order of magnitude, in Fig. 2 we have seen that it is also possible to somewhat stabilize the revivals by driving with relatively high frequencies. If we follow the same peak by reoptimizing the static detuning and amplitude while increasing the frequency, as in Fig. 2(a), the resulting driving parameters still result in a subharmonic response and its Floquet modes are split into two arcs. This suggests that we can obtain useful insights into the revival stabilization mechanism by studying the dynamics using the high-frequency expansion.
In Appendix C we derive the following third order Floquet Hamiltonian for the square pulse driving protocol in the limit of high driving frequency: In addition to the renormalized PXP term and time averaged detuning ∆ 0 , this Hamiltonian also contains a diagonal PZP term and constrained nearest-neighbor hopping terms. Note that the Floquet Hamiltonian only describes stroboscopic evolution at times t = nT , where n is an integer. We have numerically confirmed that the first three orders are enough to correctly reproduce the dynamics of the first few revivals at ω = 10. Further increase in frequency makes the agreement even better. In Fig. 7 we plot the subharmonic cost function, Eq. (5), for the Néel state at ω = 10, and compare the evolution using H correctly reproduces most features of this phase diagram and approximately locates the narrow Type-1 optimal parameter regime (around ∆ 0 ≈ 5, ∆ m ≈ 13). At the optimal parameters, the dominant term in H (3) F,ω is the static detuning, while PZP and hopping are relatively small. The energy spectrum thus resembles that of a strongly detuned PXP model, with several disconnected energy bands defined by their numbers of excitations and the gaps between them approximately equal to ∆ 0 . If we want to obtain a subharmonic response, this gap should be half the driving frequency, so we expect the optimal detuning to be around ∆ 0 = ω/2, which is indeed the case here. The best revivals will be obtained when the energy bands are very narrow, which will make the eigenstates which have the highest overlap with the Néel state almost exactly equidistant. The widths of individual bands are controlled by the other three terms. Assuming that P XP , P ZP and P σ ± σ ∓ P are approximately of the same order of magnitude and minimizing the sum of their coefficients, T 2 (3∆ 0 ∆ m − ∆ 2 m + 12∆ m )/96, at Ω = 1 and fixed T and ∆ 0 , results in the optimal driving amplitude which gives ≈ 13.5 -close to the numerically obtained optimal value. We can reproduce Fig. 2(a) for the square pulse drive (not shown) and compare the dependence of optimal ∆ 0 and ∆ m on the driving frequency with our analytical predictions from the previous paragraph. The optimal static detuning indeed closely follows the relation ∆ 0 (ω) = ω/2, even at relatively low frequencies (ω 5). As expected since the high frequency expansion is no longer valid in that regime, there are some deviations from this rule at lower frequencies, for example ∆ 0 ≈ 0.42ω at ω = 2.72. Regarding the relation for the optimal driving amplitude ∆ m (ω) = 3ω/4 + 6 in Eq. (14), the numerical results start to approach this line at higher frequencies, ω > 10. Since we have already established that the optimal values of static detuning are approximately equal for the cosine drive and square pulse, while the driving amplitude is around 20% lower for the square pulse, we can use these analytical results to estimate the Type-1 optimal driving parameters for the cosine driving protocol.

V. RELATION TO QUANTUM MANY-BODY SCARS
Finally, we come to the important set of questions concerning the relation between QMBS and the enhancement of fidelity revivals by periodic driving. We focus on two questions in particular: (i) is it still possible to stabilize the revivals from some initial state if the PXP part of the Hamiltonian in Eq. (3) is replaced by some other, fully ergodic Hamiltonian? (ii) how do we prove that the optimally-driven system is "scarred"? The second question, in particular, is non-trivial as many probes of QMBS in static systems do not directly translate to the driven case. For example, QMBS are typically diagnosed as a band of atypical eigenstates at equidistant energies. Since the driven model is time-dependent, the Floquet modes assume the role of the eigenstates. However, in all the cases studied here, there is no band of special Floquet modes in the quasienergy spectrum of the Floquet Hamiltonian. Instead, there is only one or in some cases two special modes which have unusually low entanglement entropy and very high overlap with the initial state. Also, it could be said that these modes are typically located at the edges of the quasienergy spectrum rather than the middle of it, see Fig. 3(c). However, one must keep in mind it is hard to rigorously define the "edge" when the quasienergy spectrum is periodic. In light of this difference between driven and static systems, it is pertinent to ask if the dynamics in the driven case is still "scarred" in the sense that it represents merely an enhancement of the effect already present in the static case. In this section we address these questions.
A. The effect of perturbation QMBS can be destroyed by applying suitable perturbations to the PXP model 50 . The destruction of QMBS was found to coincide with the overall level statistics of the model converging faster to the Wigner-Dyson distribution, signalling faster thermalization. Here we study the PXPXP perturbation and apply the usual spatially homogeneous cosine driving protocol, Eq. (3), to the perturbed PXP Hamiltonian: The PXPXP perturbation is off-diagonal, thus it changes the connectivity of the Hilbert space. In Fig. 8 we scan the (∆ 0 , ∆ m ) parameter space for three different perturbation strengths, while the driving frequency is fixed at ω = 2. Similar to Fig. 1, the color corresponds to the average fidelity F (t) t between t = 0 to t = 100, such that the bright regions correspond to parameter regimes where we expect robust revivals. The perturbation destroys some of the optimal parameter regimes, mostly those close to the diagonals ∆ 0 = ±∆ m . For example, the main Type-1 optimal regime for the Néel state belongs to this region, as shown in Fig. 1(a) for ω = 2.72. The higher-amplitude optimal parameter regions (Type-2) are affected by the perturbation and their position, shape and size all depend on the perturbation strength. However, these bright regions are still present even when the perturbation is fairly strong (W = 0.5). Those regions are related to Hilbert space fragmentation, see Fig. 4. This allows us to conclude that the Floquet Hamiltonian fragmentation mechanism is not specific to the PXP model, which is in line with the discussion from Section III C. However, a possibility remains that the original, weakly non-ergodic PXP model is in some way important for the revival enhancement in the low-amplitude driving regime for the Néel state (Fig. 3).
Similar behavior was also observed with different types of perturbations, such as next-nearest neighbor interaction W 1 i n i n i+2 , or constrained nearest neighbor hop- It is, however, interesting to note that, while destroying scars for the Z 2 state, the next-nearest neighbor interaction can actually create scars for the polarized state, which is reminiscent of the effects of static detuning 33 . The stabilization of revivals from the polarized state using periodic driving works even better (quantified by higher C 1 ) in a certain range of nonzero values of W 1 that coincides with the existence of quantum many-body scars for this state in the static system.

B. Trajectory in the Hilbert space
Another way to relate the reviving dynamics in the driven system with QMBS in the static case is to examine the Hilbert space trajectory of the Néel state, with and without driving. We can visualize this trajectory by plotting the expectation values of the (normalized) numbers of excitations on the even and odd sublattice, n A and n B respectively, see Fig. 9. The two Néel states only have excitations on one of the sublattices, so they correspond to the bottom right (1,0) and top left (0,1) corner, while the polarized state with no excitations is in the bottom left corner (0,0). All other states on the main diagonal or above it are forbidden due to the PXP blockade constraint.
Up to this point, for the Néel state we have mostly used the driving frequency ω = 2.72 which is approximately twice the frequency of revivals in the static PXP model (ω ≈ 2.66). Similarly, the optimal driving frequency for the polarized state in the low-amplitude regime is close to the frequency of revivals in the detuned PXP model. However, we were also able to stabilize the revivals at driving frequencies which are not particularly close to the original revival frequency. Upon first glance, this fact does not seem to support the claim that the driven system hosts QMBS, however we will show bellow that driving with frequencies close to ω = 2.72 leads to Hilbert space trajectories with some special properties.
It is possible to obtain high and long-lived revivals from the Néel state at frequencies significantly below ω = 2.72. Two examples of such trajectories are shown by the red lines in Fig. 9(a) and (b), for driving frequencies ω = 1 and ω = 2 respectively, while the trajectory in the static model is plotted in blue. For comparison, the trajectory for ω = 2.72 is shown in Fig. 9(c). In this special case, the trajectory precisely repeats the approximate first revival period of the undriven case. It could therefore be said that driving with optimal parameters stabilizes the scarred dynamics of the pure PXP model. This suggests that the true optimal driving frequency is indeed related to quantum many-body scars in the static model. The most significant difference between this trajectory and the previous two is that the low-frequency trajectories contain loops. The number of loops grows as the driving frequency decreases. Although the trajectory can still be made to oscillate between the two Néel states at lower frequencies, the intermediate dynamics differs from that in the static PXP model. However, the revivals cannot be stabilized for any arbitrarily chosen driving frequency. When the frequency is too high, for example ω = 5 in Fig. 9(d), the revivals decay and the trajectory explores only one corner of the Hilbert space, staying in the vicinity of the initial state. This could be better understood using the square pulse approximation of the driving protocol from Section IV. Close to integer multiples of the driving period, the drive is in the high-detuning regime (∆ 0 +∆ m ) where the spectrum is split into bands and the Néel state is approximately one of the eigenstates. Frequent returning to this regime therefore prevents the wavefunction from moving away from the Néel state. The intermediate low-detuning regime (∆ 0 − ∆ m ) needs to be long enough in order for the trajectory to reach the other side of the Hilbert space, and its detuning low enough so that the system is closer to the pure PXP model. In other words, the frequency needs to be low enough for the driven system to oscillate between the two Néel states.
Two remarks are in order. First, driving with highamplitude optimal parameters for the polarized state can also produce high revivals from the Néel state. This is not surprising since in this case there are two Floquet modes with high overlap with the Néel states (see Fig. 4(f)), but the wavefunction is in that case "stuck" in one corner of the Hilbert space. The trajectory is similar to Fig. 9(d), albeit more narrow. Second, representing the Hilbert space trajectory in the n A -n B plane is not suitable for the polarized state, since both that state and the full (uniformly driven) Hamiltonian are symmetric with respect to the two sublattices, thus n A = n B at all times, i.e. the trajectory starting from the polarized state always stays on the diagonal. In Appendix D we discuss an alternative way to plot the trajectory for the polarized state.

C. From scarred eigenstates to Floquet modes
Finally, we would like to explore the connection between scarred eigenstates in the static PXP model and the two π-paired (separated by ∆ = ω/2 = π/T in quasienergy) Floquet modes which have high overlap with the Néel state. To this end, we interpolate between the static (∆(t) = 0) and optimally driven system (∆(t) = ∆ 0 + ∆ m cos(ωt)) by first increasing the static detuning from 0 to ∆ 0 while keeping ∆ m = 0, and then increasing the driving amplitude to its optimal value. The driving frequency is kept constant throughout this process (ω = 2.72). The Floquet modes can be formally computed even in the ∆ m → 0 limit. In this case the Floquet modes should be exactly equal to the eigenstates of the static system. However, the Floquet quasienergy spectrum is periodic with periodicity ω, and we therefore need to "unfold" it in order to match the Floquet modes with their corresponding eigenstates. This situation is shown in Fig. 10(a) where we plot the overlap of the PXP eigenstates and Floquet modes with the Néel state. For clarity, we use a small system size N = 10. The shaded region marks the original [−ω/2, ω/2] interval of the Floquet spectrum which is translated by ω several times so that the full energy spectrum is covered.
Next, the Floquet modes are computed for various values of static detuning ∆ 0 and their quasienergies and overlap with the Néel state are plotted in Fig. 10(d).
The highest overlap states of the bare PXP model are continuously deformed into the highest-overlap states of the detuned model at ∆ 0 = 1.15, which are also shown in Fig. 10(b). Similarly, these two states are then transformed into the two special Floquet modes by increasing the driving amplitude ∆ m from 0 to ∆ m = 2.67, see Fig. 10(e). To summarize, it is possible to establish a connection between the scarred eigenstates of the static PXP model (for example, = 0 and ≈ ω/2 at ∆ 0 /∆ 0 = 0 in Fig. 10(d)) and the π-paired modes which are responsible for near-perfect revivals in the optimally driven model (two highest overlap modes at ∆ m /∆ m = 1 in Fig. 10(e)). Driving leads to scarred Floquet modes which are more separated from the bulk than in the original bare PXP model, see Fig. 10(c). The π-pairing was already present in the static model and was never broken during the interpolation procedure.
The interpolation between the detuned and optimally driven system was repeated for the polarized state and the results are shown in Figs. 10(f) and (g). Unlike for the Néel state, it is not possible to continuously deform the eigenstate of the static detuned model into the the high overlap Floquet mode of the driven system in the high-amplitude regime, see Fig. 10(f). However, as shown in Fig. 10(g), the low-amplitude driving regime is different. In that case, one of the scarred eigenstates indeed transforms into the highest overlap Floquet mode.

VI. CONCLUSIONS
In this work we studied the periodically driven PXP model in order to gain a better understanding of how driving can enhance the revivals that are present in the static model due to QMBS. We introduced a square pulse approximation for the cosine driving protocol which was used in earlier experiments, and demonstrated that this approximation much more accurately reproduces the dynamics of the cosine drive than the previously studied delta function pulse drive, while still being amenable to analytical treatment. This allowed us to analytically predict the optimal driving parameters in the highamplitude and high-frequency regimes. We were also able to modify the cosine driving protocol in such a way that any initial product state can be made to revive until late times.
We identified two distinct classes of optimal parameter regimes, which were dubbed Type-1 and Type-2. The second class is related to the splitting of the Floquet Hamiltonian into disconnected sectors, where a single Floquet mode belongs to a one-dimensional sector and has high overlap with one of the product states. The sectors consist only of Floquet modes with a certain number of excitations, thus the drive, being the largest energy scales, produces an emergent U(1) symmetry, similar to prethermalization scenario 64 . Type-2 mechanism of revival creation and stabilization is well understood, relatively trivial and likely not connected with QMBS. The optimal driving parameter regimes were found for completely arbitrary initial product states, and were shown to be insensitive to perturbations which destroy quantum many-body scars in the original PXP model.
The other class of optimal driving parameters (Type-1) is more interesting, but also more challenging to study analytically due to it being in the non-perturbative regime. This revival enhancement mechanism was so far only observed for the Néel state and is destroyed by perturbations which also destroy QMBS. Visualising the trajectory of the evolved state in the Hilbert space and comparing the static and driven cases reveals that optimal driving stabilizes this trajectory -the system is made to precisely repeat the dynamics of the first revival period of the undriven model. When the system is initialized in one of the two Néel states, the response to driving is subharmonic, meaning that the revival frequency is half the driving frequency 52 . However, the subharmonic response becomes harmonic when the initial state is a translationinvariant superposition of the two Néel states. The origin of the subharmonic response is related to the presence of two π-paired Floquet modes in the quasienergy spectrum of the optimally driven system. These two atypical modes have very high overlap with the Néel states and can be thought of as equivalents of scarred eigenstates in a periodically driven system.
An interesting direction for future work is the extension of driving protocols to other models which host inexact QMBS states and revival dynamics. This would involve finding a suitable driving term whose ground state is the scarred initial state. While Type-2 parameters (high-amplitude regime) are expected to trivially work even in the absence of scars, a more interesting question is whether Type-1 phenomenology exists beyond the PXP model. In this Appendix we examine the system size scaling of several important quantities. We fix the driving frequency and scan the (∆ 0 , ∆ m ) parameter space in the vicinity of two Type-1 peaks, (∆ 0 , ∆ m , ω) = (−1.17, 4.81, 2.00) and (∆ 0 , ∆ m , ω) = (1.15, 2.67, 2.72), for different system sizes N . The results are presented in Fig. 11(a)-(d), showing that the optimal driving parameters are relatively stable for N ≥ 10. Additionally, we plot the maximal average fidelity F t as a function of system size in Fig. 11(e)-(f). Based on these results, we predict that the same parameters will still be optimal and produce revivals of finite height in relatively large systems. Such system sizes are out of reach of numerical simulations based on exact diagonalization, but could be realized experimentally. For example, a recent experiment has realized the PXP model with approximately 50 sites using a Bose-Hubbard quantum simulator 33 . In Fig. 4(f) we showed the quasienergy spectrum of the Floquet modes which corresponded to the optimally driven polarized state in the high-amplitude regime. This spectrum is split into separate bands according to the number of excitations, with a single Floquet mode approximately equal to the polarized state belonging to a one-dimensional sector (0 excitations) and separated from the neighboring sector (1 excitation) by a gap ∆ . In Fig. 12 we show that this gap is almost independent of system size. The inset zooms in on the data points in order to show that weak exponential decay is nevertheless present. If this trend continues, we estimate a finite gap even at sizes N ≈ 300, implying that it would still be possible to create fidelity revivals in such large systems by applying the same driving parameters found in smaller systems.  We also investigated the dependence of the optimal frequency on system size. We use the same procedure as in Fig. 2 and plot the cost function C 2 as a function of driving frequency in Fig. 13. The optimal frequency window is quite stable for N ≥ 10, although it becomes slightly narrower as the system size increases.

Appendix B: Spatially inhomogeneous drive
Based on previous results, we now introduce a generalized driving protocol which can stabilize revivals from any product state. A key observation is that the two initial states whose revivals can be stabilized by driving (Néel and polarized) are the ground and the highestexcited state of the spatially homogeneous driving term in Eq. (3), since these two states have the highest and lowest possible number of excitations allowed by the blockade. In order to achieve a similar scenario for other states, the driving term should be chosen as spatiallydependent, where n 0i = ψ(0)|n i |ψ(0) is determined by the initial state. In this protocol, each site is driven by ±∆(t), where the sign depends on whether the site initially contains an excitation or not, see Fig. 14 for an illustration. In this way, the initial state becomes the (nondegenerate) ground state of the driving term. Figure 14. Modified driving protocol for an arbitrary product state. The protocol is now site-dependent, i.e., each site is driven by a term (−1) n i ∆(t) whose sign depends on whether or not the site hosts an excitation.
Using the inhomogeneous driving protocol, we were able to generate very high and long lived revivals for any initial state, even randomly chosen ones with no symmetries. The driving parameters ∆ 0 , ∆ m and ω need to be separately optimized for each initial state. To this end, we use the simulated annealing optimization algorithm, which limits the system sizes we can reach. However, when the initial state is sufficiently symmetric, for example Z 2 , Z 3 , Z 4 and combinations of such states, it is sufficient to find the optimized parameters in a smaller system and the same driving parameters can be applied to the same state in larger systems, as was the case for Z 2 with uniform driving.
Some optimal parameter values obtained by simulated annealing are listed in Table II. It is interesting to note that these values are similar in most cases, with |∆ 0 | < 1, ∆ m 1 and ω ≈ 2. From these values, we see that the phenomenology is similar to the polarized state in the high-amplitude driving regime (cf. Fig. 4(d)-(f)): the quasienergy spectrum is split into sectors with different values of i (−1) n0i n i and there is a single Floquet mode which has very high overlap with the initial state. Scanning the cost function over the parameter space at fixed frequency results in similar plots as in Fig. 1, but with Type-2 peaks only (no subharmonic response). Ultimately, the mechanism of scar enhancement in this case is very special case of prethermalization 64 : ∆ m assumes the role of the large energy scale which gives rise to an approximate U(1) symmetry. The initial state belongs to a one-dimensional sector and its mixing with other states is suppressed over long times, in accordance with general results in Ref. 64 . To illustrate the spatially inhomogeneous drive in Eq. (B1), we consider a completely random state |•••••••••••••••• . This state has no special properties and quickly thermalizes in the pure PXP model, but driving with optimal parameters results in stable revivals, see Fig. 15(a). It might be hard to see due to a large number of different sectors and relatively small system size, but the quasienergy spectrum is also fragmented in this case. As can be observed in Fig. 15(b), there is again a single high-overlap mode which belongs to a one-dimensional sector, since it is the only product state which has the minimal possible value of i (−1) n0i n i . The driving amplitude is the largest energy scale (∆ m = 6.36). (C4) Using the expression for the commutator, after some algebra, we obtain which finally leads to the expression Eq. (14) in the main text. We see from Eq. (C6) that, in addition to the PXP and detuning terms, two new terms are generated in the third order expansion: a diagonal PZP term and an offdiagonal hopping term. The "time-crystal-like" dynamics from the Néel state is generated by the interplay of such terms. The resulting Hamiltonian is sufficiently local and captures the subharmonic fidelity revival quite well. At higher frequencies, the agreement between BCH expansion and full dynamics under H(t) is systematically improved. At lower frequencies, one needs to go higher order in the expansion to ensure a good agreement. We note that although new non-local terms (with support over larger clusters of sites) will proliferate at higher orders, the terms in Eq. (C6) will also occur repetitively, thus renormalizing their strength. For example, if we expand the high amplitude expression of H (1) F from Eq. (12) we get Note that the first two terms are exactly same as in Eq. (14). Thus to get a better Floquet effective Hamiltonian, one can just re-sum the terms in Eq. (14). Numerically, one can tune the coefficients of the terms in Eq. (14) manually to obtain a better matching between the fidelity dynamics calculated using H F and H(t). Ultimately, however, the BCH expansion is expected to diverge at lower frequencies, which is at odds with the numerics that shows the existence of long-lived subharmonic response at frequencies ω ∼ 1. We leave the interesting question how to adapt the BCH expansion to describe this non-perturbative regime for future work.
Appendix D: A different way of visualising the Hilbert space trajectory In the main text, we analyzed the trajectory of a driven QMBS system in the Hilbert space by evaluating the av-erage number of excitations of the two sublattices for the |Z 2 initial state. This method is less useful for the |0 state since in this case the Hilbert space trajectory in just a diagonal line n A = n B . A more generic method for visualising the trajectory is to compute the reduced density matrix for a small subsystem, e.g., ρ 2 for a 2-site subsystem, and plot its eigenvalues. Due to the PXP constraints, there are only three possible configurations for two sites (••, •• and ••), thus ρ 2 has three eigenvalues λ i . However, only two of those are independent since tr(ρ 2 ) = 1, which allows us to plot the trajectory in 2D.
The trajectories of the polarized state in two different optimal driving regimes can be observed in Figs. 16(a) and (b). Here we use the largest and the smallest eigenvalue of ρ 2 and compare the cases without driving (bare PXP model), with static detuning only, and with periodic driving. The polarized state is located at the (1, 0) point (bottom right corner), since its reduced density matrix corresponds to a pure state (λ 1 = 1, λ 2 = λ 3 = 0). Similar to the Néel state, the driven trajectory approximately repeats the first part of the static trajectory. There are no revivals from the polarized state in the bare PXP model, so the trajectory at later times ends up in the region where all three eigenvalues are of approximately equal magnitudes, implying thermalization. In the lowamplitude regime, the static detuning alone is enough to produce revivals, as evidenced by the trajectory which periodically returns to the vicinity of the initial state, while periodic driving further stabilizes the detuned trajectory, see Fig. 16(a). In the high-amplitude regime, the revivals are created purely by driving. In that case the driven trajectory is even more narrow and returns closer to the initial state, see Fig. 16(b).
For comparison, we use the same method to plot the trajectories for the Néel state in Fig. 16(c)-(d). In these two figures we use two different sets of driving parameters, one of which corresponds to the ideal trajectory with no loops in Fig. 9(c) and the other to a "one-loop" trajectory from Fig. 9(b). However, it is now much harder to see the difference between these two cases. In both cases the largest eigenvalue stays close to 1.