Emergent strong zero mode through local Floquet engineering

Periodically driven quantum systems host exotic phenomena which often do not have any analog in undriven systems. Floquet prethermalization and dynamical freezing of certain observables, via the emergence of conservation laws, are realized by controlling the drive frequency. These dynamical regimes can be leveraged to construct quantum memories and have potential applications in quantum information processing. Solid state and cold atom experimental architectures have opened avenues for implementing local Floquet engineering which can achieve spatially modulated quantum control of states. Here, we uncover the novel memory effects of local periodic driving in a nonintegrable spin-half staggered Heisenberg chain. For a boundary-driven protocol at the dynamical freezing frequency, we show the formation of an approximate strong zero mode, a prethermal quasi-local operator, due to the emergence of a discrete global $\mathbb{Z}_2$ symmetry. This is captured by constructing an accurate effective Floquet Hamiltonian using a higher-order partially resummed Floquet-Magnus expansion. The lifetime of the boundary spin can be exponentially enhanced by enlarging the set of suitably chosen driven sites. We demonstrate that in the asymptotic limit, achieved by increasing the number of driven sites, a strong zero mode emerges, where the lifetime of the boundary spin grows exponentially with system size. The non-local processes in the Floquet Hamiltonian play a pivotal role in the total freezing of the boundary spin in the thermodynamic limit. The novel dynamics of the boundary spin is accompanied by a rich structure of entanglement in the Floquet eigenstates where specific bipartitions yield an area-law scaling while the entanglement for random bipartitions scales as a volume-law.

Certain one-dimensional many-body Hamiltonians are endowed with a quasi-local operator called a strong zero mode (SZM), which exists on the boundary of the system and is conserved for long times [39,40].This operator can be utilized to encode stable quantum information away from the ground state and does not require low temperatures.The existence of SZM can be proven in integrable models while in the presence of weak integrability-breaking perturbations, these quasi-local operators are expected to exhibit prethermal behavior [41][42][43].Models with Hilbert space fragmentation can host a statistically localized SZM, without the presence of integrability [44].The connection between integrability and SZMs is of significant interest for the breakdown of thermalization.In this context, symmetries play a crucial role in fragmenting the Hilbert space and forming quantum scars, which stabilize SZM.This class of phenomena falls under the rubric of partially integrable models, where local dynamical properties exhibit non-thermal behavior.Under a global periodic drive that breaks integrability while preserving the global symmetry, the SZM destabilizes and has a finite lifetime [45,46].These boundary observables take a parametrically long time to thermalize while the bulk of the system thermalizes quickly.These approximate zero modes are continuously connected to exact zero modes in the integrable limit of the models.Indeed, in most cases, the presence of a global drive appears to be crucial for realizing prethermal or athermal behavior due to the emergence of conservation laws.
Generically, an interacting Floquet system heats up to an infinite temperature state, however, prethermal or slow dynamics can be induced due to the separation of scales between the frequency of the drive and the system's internal energy scales [47][48][49][50][51].The distinct stroboscopic dynamical regimes of a periodically driven system are governed by the Floquet Hamiltonian H F , defined through the propagator U (τ ) as where H(t) = H(t + τ ).H F plays a central role in understanding the slow dynamics or heating to infinite temperature [52][53][54].A perturbative description of the prethermal timescales can be developed in terms of an expansion of H F in τ , where τ is the time-period of the Floquet drive.In general, the perturbative corrections in a periodically driven nonintegrable system become non-local and the series diverges in the thermodynamic limit, signaling rapid thermalization.Therefore, in finite size systems the choice of drive protocol which FIG. 1.A schematic figure showing the model interaction and drive protocol adapted in this work.Red (blue) color is used to denote ferro-(antiferro-) magnetic exchange in consecutive links.A stationary global field is applied in the z-direction whereas a time-periodic field in the x-direction is applied at a few specific sites.
increases the radius of convergence of such perturbation series is of importance for Floquet engineering.In the presence of a high-frequency drive, nonintegrable systems can exhibit prethermal behavior where emergent conservation laws are applicable for exponentially long times.In another class of phenomena global drives can result in the formation of Floquet quantum scars [55][56][57][58][59][60] and dynamical many-body freezing [61][62][63], where local observables remain athermal under unitary evolution for a specific choice of drive protocols -an exclusively drive-induced non-thermal phenomenon at intermediate frequencies.The frozen states are characterized by local conservation laws and lead to scarring of the Floquet spectrum.
A question naturally emerges, what happens when a system is driven locally?Unlike the global drive which can preserve translational symmetry, a local drive explicitly breaks it.Recently, the real-space profile of thermalization for local or spatially modulated drives has been shown to exhibit rich behavior where the drive can disentangle spins leading to large variations in thermalization times and even steady-state properties [64,65].The disentangling effect of the drive can generate cold spots which persist for significantly longer than thermalization times in the rest of the system.
In this work, we consider the fate of thermalization for a boundary-driven nonintegrable spin chain.In the absence of driving, the model hosts a set of scarred eigenstates that act as a large spin given by the SU(2) symmetry of the model, showing persistent oscillations for certain unentangled initial states.In the presence of a local drive, this manifold is destroyed by a novel form of dynamical disentangling of the spin due to the emergence of a discrete symmetry at the dynamical freezing frequency.We unravel a new mechanism for local ergodicity breaking due to the formation of an approximate SZM at a certain discrete set of frequencies.Furthermore, using Floquet engineering of a cluster of sites (see the illustration in Fig. 1) we are able to control the exponential time scales of boundary relaxation characterized by the emergence of an asymptotically exact SZM.The athermal dynamics is accompanied by a rich entanglement structure of Floquet eigenstates which exhibit both thermal and athermal properties.
The rest of the paper is structured as follows.In Sec.II we introduce our model, and discuss its symmetries and nonintegrable nature.This is followed by the definition of the local Floquet protocol.In Sec.III we discuss the boundary-driven protocol where the drive acts on a single site.Using higher-order Floquet-Magnus (F-M) expansion, we demonstrate the slow dynamics of the edge spin and the resulting properties of the approximate SZM.In Sec.IV, we give the full phenomenology of two and multi-site driving, discovering the optimal drive protocol which leads to the freezing of the boundary spin.We discuss how this novel dynamical property is accompanied by the emergence of a SZM in Sec.IV C. In order to interpret our results on the dynamics, we analyze the entanglement structures of the Floquet eigenstates in Sec.V, and show their thermal and athermal properties for specific bipartitions.Sec.VI summarises the salient results and their implications for athermal behavior in Floquet systems.
We consider the one-dimensional "staggered" (alternating ferro-antiferromagnetic exchange interaction) nearest neighbor spin-1/2 Heisenberg model with a globally applied homogeneous magnetic field, as the undriven part of our system.Its Hamiltonian is, where i refers to a site index and S i ≡ (S x i , S y i , S z i ) is the vector comprised of the usual spin-1/2 operators, and h is the magnetic field strength.N is the number of spins and we always consider open boundary conditions (OBC).As previously mentioned elsewhere by some of us [64], the model in Eq. ( 2) has a long history partly due to its relevance to the Haldane spin-1 chain [66][67][68][69], which can be realized in the limit of the ferromagnetic bonds being much stronger than the antiferromagnetic ones.
For h ̸ = 0, the magnetic field term reduces the symmetry of the Hamiltonian from SU(2) to U (1), thus the total magnetization along the z axis (S z tot = i S z i ) is conserved.A subset of U (1) is the Z 2 symmetry, i.e.H 0 commutes with the D z operator, where with α = x, y, z.In terms of spatial symmetries, H 0 is reflection symmetric (w.r.t. the central bond) only for even N .For odd N , H 0 is found to anticommute with the operators for α = x, y, where we have defined, R is the reflection operator w.r.t. the central site and SWAP(i, j) is the SWAP gate, defined as SWAP ≡ |↑↑⟩⟨↑↑| + |↑↓⟩⟨↓↑| + |↓↑⟩⟨↑↓| + |↓↓⟩⟨↓↓|.Consequently, H 0 has a spectral reflection symmetry (i.e. if there is an eigenstate |ψ⟩ of energy E, then there is also an eigenstate A |ψ⟩ with energy −E) for odd N .
We find that there are exponentially many (with system size) exact zero energy states at h = 0 for odd N , intriguingly such a feature has also been reported in other models that harbor a superspin structure, which translates to the existence of QMBS [24,[70][71][72][73][74].For any eigenstate |ψ⟩ with magnetization m z , there is always an eigenstate R |ψ⟩ (D α |ψ⟩, for α = x, y) with magnetization m z (−m z ).For even N , both the spectral reflection (E → −E) symmetry and exact zero energy states are lost but the latter property still remain valid.
Despite its superficial resemblance to the integrable uniform Heisenberg chain, the staggered model is nonintegrable.Focusing on even N , here N = 16, we work in the zeromagnetization sector (S z tot = 0) and switch off the magnetic field at the last site to remove the presence of conventional symmetries (e.g. S 2 conservation, inversion symmetry) [75], which must otherwise be accounted for in the determination of level statistics.The middle 60% of the many-body spectrum is used for the purpose of computing level statistics, this is done to focus on properties at infinite temperature.Due to the limited statistics of eigenvalue spacings available for every bin, it is instead useful to monitor the integrated level spacing distribution, where P (s ′ ) is the probability density of the unfolded energy level spacings (spacings divided by their mean).Note that I(s) is simply the cumulative distribution function of the unfolded spacings s.
In Fig. 2 we plot I(s) for the case of h = 1, and find that it closely resembles the results of a Gaussian orthogonal en-semble (GOE) of the random matrix theory [76], a signature of chaotic, nonintegrable systems [77].In comparison, when the bonds are made uniform (instead of staggered), we find that I(s) follows the Poisson distribution, as expected for integrable models.

B. Strong zero mode
Our work here is inspired by the aim of engineering a strong zero mode, a quasi-local conserved operator, localized at the edge of a one-dimensional quantum system which gives rise to long coherence times of the edge spin [39-41, 44, 78].Given previous work on the search and characterization of such modes, it has been established that the SZM operator Ôszm in a spin-1/2 system must satisfy the following conditions: 1.There exists a discrete symmetry operator D in the system, i.e. [D, H] = 0, implying the eigenstates are split into two sectors, labeled by the eigenvalues of D.

4.
Ôszm asymptotically commutes with the Hamiltonian H of the system, i.e. ||[H, Ôszm ]|| ∼ exp(−aN ) with a > 0. A consequence of this is that Ôszm is a conserved quantity in the thermodynamic limit and that the gap falls off exponentially in the system size.
We, therefore, expect that a system endowed with a SZM is characterized by a spectrum-wide (quasi)-degeneracy and relaxation of the edge that can be slowed down by increasing the system size.It is also important to satisfy all these conditions simultaneously to have SZM in a system, for example, one can construct an operator which satisfies the last condition but does not correspond to a SZM.The undriven staggered Heisenberg chain (H 0 ) is not expected to possess a normalizable SZM and, consequently, any edge spin should quickly relax to zero starting from any randomly chosen ordered state.
In a bid to effectively decouple the edge spin from the rest of the system, we adopt the strategy of driving it locally with a protocol we describe next.

C. Drive protocol
Throughout this work, we adopt a Floquet drive protocol that involves applying a time-periodic magnetic field in the xdirection on a subsystem that includes the boundary spin.The total time-dependent Hamiltonian is given by where sgn(x) is the sign function, and thus its profile is that of a "square pulse", and s d is the set of driven sites (see Fig. 1).γ and ω are the drive amplitude and frequency, respectively and τ ≡ 2π/ω is the time duration of one drive cycle.The time-averaged field in the x-direction is zero at stroboscopic times (integer number of drive cycles).We note that some aspects were previously considered for the special case of a single driven bulk site [79][80][81].
The time evolution operator over one drive cycle can be written as Since we focus on stroboscopic dynamics, t = nτ where n is an integer, analyzing the properties of U (τ ) or equivalently H F is adequate for long time dynamics.A key contribution of our work is to explicitly derive H F , to some approximation, in a bid to understand the exact stroboscopic dynamics seen in our numerical calculations.The explicit form of H F was computed to the lowest order in the F-M expansion for the case of a single driven site in the bulk [64], a result we revisit in the next section for a boundary driven case.

D. Observables
Both from an experimental and theoretical point of view it is advantageous to monitor the time dependence of local operators, here we focus on ⟨ψ(nτ )|S x i |ψ(nτ )⟩, where is the time-dependent state starting from an initial state |ψ 0 ⟩.The initial state is typically chosen to be "simple", for example, a product state in the x basis, with the hope that it can be easily realized in experiments [16,[82][83][84].
Further theoretical insights can be obtained from the timedependent von Neumann entanglement entropy, where A refers to a single site or collection of sites.The reduced density matrix of this collection of sites, ρ A , is computed by tracing out all degrees of freedom not part of A (labeled by Ā) in the full density matrix In this work, we choose A to include either a single site, or a collection of driven sites.Entanglement entropy is a useful metric for distinguishing between low-entangled states (area law or lower) or those with high entanglement (typically, volume law).Thus, it reveals the existence or absence of (Floquet) ETH and is also useful for the detection of QMBS [24,29].

III. BOUNDARY-DRIVEN APPROXIMATELY-CONSERVED STRONG ZERO MODE A. Overview of salient features of dynamics
Given the setup in the previous section, we study the stroboscopic relaxation of the spins when only the left boundary spin is driven i.e. s d = {1}.The main purpose is to explore the existence of parameters that freeze the dynamics of the edge spin via the emergence of a conserved quantity (for example, S x 1 ) which can lead to realizing a strong SZM.Instead, we find an approximate SZM for the boundary-driven protocol, i.e., no operator satisfies the criteria outlined previously.However, the lessons learned from understanding this case facilitate a series of drive protocols that lead us to a SZM, demonstrated in the following sections.
We summarize our key physical findings in this section, which relies on constructing H F using a resummed F-M expansion to different levels of accuracy.First, we revisit the result for H F to zeroth order (l = 0) in the F-M expansion, by transforming to the rotating frame.We show that this approximation, while qualitatively correct for many situations, fails to capture the dynamics evaluated from exact numerics near certain special drive frequencies which we dub as "freezing condition" [62].This failure motivates the technically challenging computation of H F to second order which contains multiple terms, we classify as "local" and "non-local".We find that the local terms, treated accurately to high order in the inverse drive strength, renormalize the freezing frequency in excellent agreement with the exact results.Finally, we demonstrate that describing the slowdown of the boundary spin at and near the freezing frequency requires the incorporation of additional non-local terms.We initialize our system to a simple product state with all spins pointing along the −x axis, In the absence of an applied magnetic field (h = 0), this state, being the maximally polarized ferromagnet along the −x direction, is an exact eigenstate of H 0 with energy E = 0 for odd N and E = 1/4 for even N .When decomposed in terms of S z tot projected eigenstates, we get When a nonzero z-field (h ̸ = 0) is turned on, the projected S z tot states (on the R.H.S.) remain eigenstates of H 0 , but their degeneracy is split.Thus |ψ 0 ⟩ no longer remains an eigenstate of H 0 .These projected states were shown to have low entanglement [29], in violation of the expected volume law at infinite temperature.Thus, this state serves as perhaps the simplest example of an embedded superspin (at high energy) in an otherwise chaotic spectrum, this is the basic idea in the core of perfect QMBS [29,85].
Though the system initialized with |−X⟩ exhibits perfect oscillations (in the Loschmidt echo and various observables) under the action of H 0 , the situation is starkly different when the system is subject to a local periodic drive [Eq.(7)].In previous work [64], we found that there are at least two parametric regimes depending on the kick strength and frequency.Typically, the driven spin thermalizes collectively with the rest, however, when a freezing condition involving the drive strength and frequency is satisfied, the driven spin is almost (but not strongly) decoupled from the rest of the system.We will elaborate on this in the next subsection.
In what follows we have worked primarily with the |−X⟩ initial state, however, we emphasize that all our results remain insensitive to this choice and are valid for any product state in S x -basis.

B. Failure of the zeroth order Floquet Hamiltonian
We perturbatively construct and analyze the Floquet Hamiltonian, H F , using a resummed F-M expansion (as charted out in detail in Appendix A).In this scheme, where the various orders are denoted by H (l) F and, in principle, can be calculated for each l.
A computation of the l = 0 term gives, H where λ ≡ γτ /2 = πγ/ω.Additionally, the calculation for l = 1 shows that H F vanishes for all drive parameters.When λ = 2πk, where k is an integer, the only non-zero term involving the driven boundary site in This suggests that the edge spin can be completely frozen if one tunes the drive frequency to ω k c = γ/(2k) since the following commutation relation F + H is exactly satisfied.We focus on the k = 1 case, though all the discussions remain valid for any k.Fig. 3 (a) presents the stroboscopic dynamics of S vN i , for all i with s d = {1}.We find that the relaxation rate of the boundary spin is indeed very small at ω ∼ ω k c (compared to other sites which thermalize very quickly in spite of being undriven), however, it is always nonzero.This can be seen in the inset of Fig. 3(a), where we plot S vN 1 averaged over the first 2000 cycles with ω near ω 1 c .The frequency which gives the slowest relaxation (ω = ω k m ) shifts to a somewhat higher value compared to ω k c (i.e.ω k m > ω k c ).This result can not be explained by the lowest-order (l = 0) Floquet Hamiltonian.
We also monitor the size dependence of the relaxation rate of S vN 1 in Fig. 3(b) using a combination of exact diagonalization for N ≤ 14, and matrix product state based time-evolving block decimation (TEBD) technique [86,87] for N > 14, setting the bond dimension to its maximal value of 2 N/2 .A time step of δt = τ /50 was used for the TEBD calculations [88], and we checked there was almost no difference in our results for δt ∼ τ /30 to δt ∼ τ /100.While increasing the system size, we find that the relaxation of the edge spin slows down, but most likely not exponentially.Moreover, the relaxation dynamics tends to be saturated beyond a threshold system size (≈ 16 as can be seen in Fig. 3(b)), an observation that is consistent with the phenomenology of an approximate SZM [40].

C. Higher order corrections and freezing condition
To explain our numerical observations, we need to go to higher terms in the F-M expansion.However, before we dive into this calculation, we note that the F-M expansion to obtain the Floquet Hamiltonian for interacting systems is plagued by various convergence issues [49].Though these issues can be circumvented to some extent by systematic resummations in certain parameters (e.g.drive frequency) [52,89,90], the calculation of the resummed expansion becomes increasingly difficult in higher order in other parameters (e.g.drive amplitude) due to the presence of nested commutators of manybody terms.In spite of this, the calculation of higher-order terms in F-M has been fruitful, as it can reveal physical effects not accessible in the lower order expansion.

As anticipated, the calculation of H (2)
F is somewhat complicated by the existence of multiple terms, for completeness the various contributions have been described in Appendix A.
Here we categorize the resulting expression into "local" and "non-local" parts, H The "local" part, consisting of only one-and two-body (n.n.) operators (including the boundary site), can be written as where p ≡ (h, γ, τ ) represents the set of drive parameters (see Appendix A for the full expression of all the a i coefficients).The ellipsis designates terms, which commute with S x 1 .It is important to note that our classification of "local" is based on the fact that the terms have the same functional forms as those appearing in H (0) F , hence only renormalizing their strength.Interestingly, we find that a i s contain terms with similar magnitude as the zeroth (l = 0) term [i.e.O(1/γ)], some of which can be nonzero even at the special freezing frequencies ω k c .Collecting only the O(1/γ) terms from H F, loc , we rewrite H F as: F + H (2) With this expression it becomes clear that In fact, there is not a single p for which all S x 1 non-commuting terms in Eq. ( 18) simultaneously go to zero.
In other words, it is impossible to dynamically freeze the edge spin via the emergence of a corresponding local conserved quantity.This result must be contrasted with the case of H F for the case of global driving with a square pulse [62], where the contributions are expected to be at least O(1/γ 2 ) or smaller (see Appendix C).
The local terms, to order 1/γ, do not significantly affect the freezing condition obtained from H (0) F , they alone are quantitatively inadequate for explaining the observed shift in the freezing frequency.Thus we retain higher order local terms in F,loc and define F + H (2) To demonstrate that the Hamiltonian in Eq. 19 captures the shift of the point of minimum relaxation to a frequency slightly higher than ω k c , we calculate a matrix element of H F in the S x -basis, H F (2, 1) where |1⟩ = |ψ 0 ⟩ and |2⟩ = σ z 1 |1⟩.Fig. 4 shows a comparison of different levels of approximation for H F (2, 1) with the exact numerical calculation.One can see that |H F (2, 1)| obtained from Eq. ( 18), though always nonzero (supporting the absence of freezing of the edge spin), does not match the exact results, particularly it still shows the minimum at a frequency which is extremely close to ω k c , F yields substantially good agreement with the exact numerics, in particular the location of the minima of matrixelements (like H F (2, 1)) which determines the freezing frequency.
Going beyond individual matrix elements, we compare the exact stroboscopic dynamics with that generated by the H F,loc [O(1/γ 3 )] in Fig. 5.This approximation is excellent for capturing the dynamics of all spins away from the freezing frequency.At and near the freezing frequency, it captures the dynamics of all, but the boundary spin, accurately.Thus the local H F approximation can capture not only the individual matrix elements but also the stroboscopic dynamics quite accurately at least for short times.

D. Role of non-local terms to slowdown the edge-spin
Interestingly, the failure of the local F only for the dynamics of the boundary spin at and near ω k m reveals the importance of the non-local terms, that we have ignored up to this point.We have numerically checked that near the freezing frequency many matrix elements of H F , generated by long-range and multi-spin interaction terms involving the boundary spin are quite significant relative to the local matrix elements [like H F (2, 1) which shows a minimum at this parameter regime].We find that such non-local terms, appearing in H (2) F , can at most be four-body operators having support on sites 1-4.We have calculated all such terms (dubbed as H   The non-local terms, though very weak, will proliferate in further higher orders; a careful inclusion of them may increase the timescale of the agreement with exact numerics.In fact we find that the non-local and multipsin terms already present in H (2) F,nonloc are sufficient to capture the slow dynamics of the edge spin quite accurately up to n ∼ 200.The strength of the nonlocal terms calculated in 2nd order (l = 2) will be renormalized in higher order in the F-M expansion which in turn will extend the agreement with exact dynamics to longer times.This is quite nontrivial for two reasons.First, any nonthermal (or prethermal) phenomenon in a Floquet system is naively expected to be linked with the existence of a local H F whereas a non-local H F usually promotes heating.But here we observe the exact opposite, the thermalizing dynamics of the bulk sites is well captured by a local H F (even at ω k m as shown in Fig. 5) but not the freezing dynamics of the boundary site which necessitates the presence of certain longrange terms for its onset.Here we also note that H (2) F,nonloc is completely off-diagonal in the x-basis and the proliferation of these terms is expected to dephase the initial state |ψ 0 ⟩.Second, from a practical point of view, it establishes the freezing phenomenon reported in this work as a genuine drive-induced many-body effect, which is challenging to engineer with finetuned long-range terms in a time-independent Hamiltonian.
As noted previously, the rate of relaxation of the boundary spin at the special frequencies (ω k m ) decreases further with increasing system size, as can be seen in Fig. 3 (b).This is at odds with the usual expectation that increasing the system size generically increases the phase space available for thermalization and should thus speed up the process.This apparent anomalous behavior already suggests that boundary driving alone can be enough to trigger the onset of an approximate SZM [40,43] in the system.However, the relaxation time of the boundary spin is significantly shorter than the time associated with a SZM, in which case it is exponentially large in system size.
The inability of boundary-site driving to fully freeze the boundary spin and the presence of non-local terms motivates the next natural question: what happens when multiple sites (including the boundary one) are driven together?Is there any requirement on the number of driven sites as well as the specific locations where the drive should be applied to maximize the lifetime of the edge spin?We address these questions in the next section.

IV. EMERGENCE OF THE STRONG ZERO MODE VIA MULTI-SITE DRIVING A. Optimal drive protocol
We now consider a situation when multiple spins in the bulk are driven simultaneously, in addition to the boundary.Intuitively, the increase in the number of driven sites would lead to more rapid thermalization of all the spins.As discussed in Sec.III D, each driven site generates a larger collection of non-local terms which slow down the edge spin near the freezing frequency.In this section we first explore the landscape of few-site protocols and investigate their effects on the edge spin.This subsequently leads us to the development of a multi-site driving protocol, resulting in the further slowdown of the edge spin.
When the boundary and a generic bulk site are driven together, our results indicate that the slowdown of the boundary site is not improved.We demonstrate this behavior with an example for s d = {1, 3} in Fig. 6(a,c), where we observe that the relaxation of the boundary site is faster than the case of s d = {1}.The undriven sites are found to thermalize quickly, with the exception of the intermediate spins, which exhibit slower relaxation to thermal equilibrium.However, remarkably we find that a particular two-site driving protocol results in strong freezing of the edge spin: when in addition to the boundary site the fifth site is also driven (s d = {1, 5}).Indeed, the first site shows almost no relaxation of S x 1 for the drive cycles plotted in Fig. 6(b,d).The intermediate sites assume a different character now, showing a strong temporal oscillation before eventually decaying out to zero.This particular driving protocol also causes a widened frequency window which admits slow (but nonzero) relaxation of the edge spin.Moreover, we find that the relaxation of the edge spin, analogous to the boundary-driven case, slows down with increasing the system size [see the inset in Fig. 7].
This freezing behavior is necessarily a high-order manybody effect, as the drive itself only produces interactions with a range of three sites away from the driven one (in the 2nd term in F-M).Therefore, driving site 5 does not change the structure of the effective Hamiltonian from Eq. ( 18), only adding terms supported on sites 2-8.This implies that the origin of the suppression of the boundary relaxation is a higher- order many-body process, where the driven site 5 impacts sites 2-4, which in turn impact the boundary.This complex collective behavior seems to be a unique property of only this particular two-site protocol (s d = {1, 5}), and will become the key ingredient in our subsequent proposal for the Floquet engineering of the SZM.
We also note that the additional driven site is itself not frozen.This can be understood by noticing that the Hamiltonian now includes newly generated three-body interaction terms (e.g. S x 4 S z 5 S x 6 , see Appendix B), which renormalizes the field strength at site 5.Note that such a term cannot occur for the boundary-driven case.To illustrate this one can, for example, compare the aforementioned matrix element H F (2, 1) with the matrix element representing single spin-flip events at site 5, ⟨ψ 0 | H F σ z 5 |ψ 0 ⟩.Indeed, we find the latter to be larger in magnitude, contributing to a much faster dephasing of site 5 as compared to site 1.

B. Freezing of edge spin on increasing the number of driven sites
Motivated by our observation of the slowdown of the edge due to the two-site protocol, we consider the driving of additional sites.Increasing the number of driven sites further, we find that generically there is a rapid relaxation of all sites, including the boundary.However, when the driving sequence is specifically chosen as s d = {1, 5, 9}, the freezing of the boundary site is significantly enhanced compared to the twosite protocol with s d = {1, 5}.We can conclude that the further slowdown has a similar origin to the two-site protocol: a high-order many-body process, this time involving site 9. Based on the spatial periodicity of the few-site protocols, we propose the following Floquet engineering protocol to further slowdown the edge spin, where the driving sequence is 1 vs n for s d = {1, 5} for different system sizes.All parameters are the same as in Fig. 6. composed of the cluster s d = {1, 5, 9, ..., 4k + 1}.Subsequently, we will use the term "every fourth site" to describe this cluster, and designate the total number of driven sites as n 4th d .The existence of this protocol does not preclude other protocols which can stabilize a SZM, and we consider this as a particular case in a general class of locally driven protocols.
In Fig. 7 we show how the growth of entanglement entropy of the boundary site decreases while increasing the size of the cluster s d : we find a hierarchy of freezing behaviors with the growing number of driven spins.This conclusion is supported by the exact diagonalization results allowing us to reach large number of drive cycles, see Fig. 7. Furthermore, the inset in Fig. 7 shows that although increasing the system size alone improves the freezing behavior, the slowdown can be most effectively achieved by also increasing the number of driven sites.
Intuitively, the extremely slow dynamics in the proposed protocol is directly related to the gaps between quasi-energy levels in the Floquet spectrum.This can be seen by considering thermalization of the expectation value of a local observable for a driven state |ψ(nτ )⟩ [91][92][93]: where ϵ r and |Φ r ⟩ are the r-th quasienergy and eigenstate in the first Floquet Brillouin zone (BZ).The first sum in the equation controls the relaxation speed, while the second sum gives the long-time expectation value in the thermalized state (diagonal ensemble average).For our protocol, the zerothorder effective Floquet Hamiltonian H F along with the 2 nd order corrections, are invariant under S x → −S x transformation, which constraint the diagonal ensemble average of S x 1 to be zero.Since our exact calculations of the expectation value of S x 1 in Floquet eigenstates also yield zero up to machine precision for a wide range of parameters, the invariance is expected to exist at all orders of H F .
For generic non-integrable systems the off-diagonal contribution decays to zero quite rapidly.But here we find that the energy gaps between the states (say r and q) which dominate the matrix elements are extremely small.For those states |ϵ r − ϵ q | ranges from ∼ 10 −5 for s d = {1, 5} dropping to ∼ 10 −8 for s d = {1, 5, 9, 13}.This is further supported by the level statistics of the Floquet Hamiltonian, which shows a strong Poisson distribution and the absence of level repulsion, putatively indicating the emergence of a symmetry or fragmentation of the Hilbert space (see Appendix D).We provide a coherent analysis of this extensive quasi-degeneracy throughout the entire Floquet spectrum, which shows that the freezing ultimately has its origin in the emergence of the SZM.

C. Emergence of the strong zero mode
In this section, we elucidate the step-by-step emergence of the strong zero mode which is the crux of the exponentially slow dynamics of the edge mode operator.We expect that a system endowed with a SZM, as defined in Sec.IIc, is characterized by a spectrum-wide (quasi)-degeneracy and relaxation of the edge that can be slowed down by increasing the system size.
First, let us consider the existence of a discrete symmetry operator.The undriven system has a Z 2 symmetry of ), yet the drive in the x-direction generally destroys it.However, we can exploit the freedom in the driving parameters and note that near the special frequencies ω k m the terms in Eq. ( 18) that do not commute with D z vanish.We test this for the full Floquet Hamiltonian numerically in Fig. 8(a), where we plot |⟨D z ⟩| = 2 N r=1 |⟨Φ r |D z |Φ r ⟩|/2 N ; this quantity approaches 1 when each Floquet eigenstate is simultaneously an eigenstate of the proposed symmetry operator.We find that although D z is not an exact symmetry, it is nearly conserved at the special freezing frequencies.The formation of SZM requires a special pairing structure of the Floquet eigenstates and quasi-degenerate eigenvalues, along with the global discrete symmetry D z .
We now label the Floquet eigenstates by the eigenvalues of the symmetry operator D z as We show that the eigenstates occur in pairs in the Floquet spectrum where the operator which toggles between the states |Φ + r ⟩ and |Φ − r ⟩ is localized at the edge.This operator anticommutes with the global symmetry operator D z .We propose that an approximate operator, which satisfies these properties, is O szm ≈ S x 1 .This is shown by the numerical calculation of the matrix elements between |Φ + r ⟩ and |Φ − r ⟩.In Fig. 8(b) the average matrix element over all pairs of eigenstates is plotted as a function of ω.This quantity measures how S x 1 connects the two symmetry sectors, and should approach 1/2 when the pairing is exact, when the matrix element is non-zero for a single pair of Floquet eigenstates.Indeed, we find that increasing the number of driven sites leads to a more accurate m which implies the emergence of SZM induced quasidegeneracy between the Floquet eigenvalues from different D z sectors.N = 10, γ = 15 for all the plots.We have performed a moving average over the raw data (in an ω interval of 0.008) for panel (b).The vertical dashed line is a guide to the eye to the position of strongest freezing.
pairing of the states.We propose an ansatz for the Floquet eigenstates with the following structure: where states |ξ ± r ⟩ and | ξ± r ⟩ obey conditions ⟨ξ This is a general entangled state between the edge spin and the rest of the system with the edge spin constrained to satisfy ⟨S x 1 ⟩ = 0.The structure of entanglement between them is encoded in the overlap properties of |ξ ± r ⟩ and | ξ± r ⟩, which will be highlighted in Sec V.These states can be utilized to encode the SZM operator with the leading order term dominated by S x 1 .
The corrections to the leading order term are controlled by the extent of deviation from the complete overlap between |ξ + r ⟩ and |ξ − r ⟩ states.Despite the rapid thermalization in the bulk, on increasing the cluster of driven sites we find a proliferation of the number of Floquet eigenpairs satisfying the condition R r ≈ 1, ∀r (see Appendix E).Therefore, the form of the strong zero mode is indeed well-approximated by Eq. 23.The existence of this normalizable operator explains the degeneracy present in the spectrum and relates it directly to the pairing of Floquet eigenstates.
Finally, we turn our attention to the commutation properties of Ôszm with the Floquet Hamiltonian.The commutator assumes the following form in terms of eigenstates and eigenvalues of H F , (24) which implies that all spectral gaps must simultaneously vanish exponentially with N for the SZM to be a quasi-conserved operator.In Fig. 8(c) we present results for the average pairing gap where indeed the gap decreases significantly with increasing size of the driven cluster (n 4th d ) in the proximity of the freezing frequency.Interestingly, when we fix n 4th d and start increasing the system size, the gap decreases weakly with N [see Fig. 9(a)].We associate this slow decay of ∆(N, n 4th d ) with N to the convergence properties of the corrections in Eq. ( 23), which appear to converge for smaller N , with ||[H F , Ôszm ]|| decreasing in magnitude.Although, subsequently for larger N the corrections dominate Ôszm with the commutator reaching a finite threshold signalling the presence of an approximate SZM (as witnessed previously for the boundary driven case in Sec.III).On the other hand, if the system size is fixed, the gap exhibits an exponential fall with n 4th d as shown in Fig. 9(b) [see Appendix F for more details].This behavior of the pairing gap explains our results for the freezing dynamics of the boundary site in Fig. 7, where we found a similar dependence of the relaxation times on the system size and the driven cluster size.Therefore, we conclude that increasing both N and n 4th d leads to an asymptotic exponential decrease of the average pairing gap and hence the commutator of Ôszm with the Floquet Hamiltonian when n 4th d scales linearly with N .
In summary, we have demonstrated that the protocol involving driving every fourth site in the chain leads to the emer- gence of a SZM.We have proposed an approximate form of an operator that satisfies all the features of SZM which is directly related to the spectrum-wide quasi-degeneracy and slow relaxation rate of the boundary spin.We note that any bulk site can not be frozen in this manner which suggests that a strictly convergent SZM, localized in the bulk, can not be constructed.

V. ATHERMAL ENTANGLEMENT STRUCTURE OF FLOQUET EIGENSTATES
In this section, we discuss the non-equilibrium properties of the system through the lens of the entanglement entropy of Floquet eigenstates.The athermal and thermal behavior represented in the dynamics and spectral properties of the boundary and bulk degrees of freedom respectively, have an imprint on the entanglement features of the eigenstates as well.In fully thermal eigenstates at infinite temperature, the reduced density matrix of any subsystem exhibits maximal entanglement.The disentangling nature of the periodic drive introduces several novel features in entanglement which we categorize between thermal and athermal depending on the partitioning of subsystems.
The von Neumann entanglement entropy of the boundary site with the rest of the system in the Floquet eigenstates (Eq.22) can be written as In the protocol with s d = {1}, |δ ± r | is found to be large for many eigenstates, thus giving rise to many weakly entangled states, nonetheless in a large majority of states the edge spin remain maximally entangled with the bulk as shown in Fig. 10(a).In the other limit, for s d = {1, 5, 9, 13} the boundary spin is maximally entangled in all the Floquet eigenstates [see Fig. 10(b)] which also implies |δ ± r | is vanishingly small.This is counterintuitive since, the enhanced freezing of the boundary spin is expected to lead to disentangling of the edge spin in the Floquet eigenstates, which instead appears fully ergodic for single-site observables.Therefore, the dynamics of the boundary spin, however slow it may be for finite systems, is only a prethermal phenomenon with the athermal behavior emerging in the thermodynamic limit.The boundary spin takes an exponentially long time in N to dephase and finally reaches a featureless state for any finite system size.
However, the scaling of the Floquet eigenpair gap in Fig. 9 suggests that it should vanish in the thermodynamic limit.Consequently, the decoherence time of the boundary spin diverges as well as the states |Φ + r ⟩ and |Φ − r ⟩ become degenerate which allows us to choose 2 as the Floquet eigenstates, where the boundary spin is disentangled from the rest of the system.Thus, the apparent paradox between the stroboscopic dynamics of the edge spin and the ergodic properties of the Floquet eigenstates is resolved in the thermodynamic limit.
Interestingly, we do find athermal entanglement signatures in the Floquet eigenstate when every fourth site is driven and the system is bipartitioned into driven and undriven sites [see Fig. 10(e)].For many eigenstates, the entanglement remains close to ln 2, as can be seen in Fig. 10(d).Note that typical states in general should satisfy volume-law scaling for such bipartitions, an example being the protocol s d = {1} as shown in Fig. 10(c).Therefore, our most efficient protocol (s d = {1, 3, 5, 9}) generates exponentially many Floquet eigenstates with subextensive (athermal) entanglement, for a bipartition separating the driven and undriven sites.However, any random bipartition such as half-chain entanglement appears to produce a volume law scaling which is reminiscent of the behavior of rainbow states [94,95].Our entanglement results showcase a rich and complex structure based on the geometry of the driven sites, where we observe signatures of both prethermal and athermal behaviors.

VI. CONCLUSION AND DISCUSSION
The phenomena of SZM where a quasi-local conservation law localized on the boundary, is proven to exist in onedimensional integrable models with a discrete symmetry.In non-integrable models, this operator usually develops a finite lifetime and ceases to be a conserved quantity in the thermodynamic limit.In this paper, we have shown the emergence of SZM in a nonintegrable spin chain through local Floquet engineering.We begin with a time-independent non-integrable Hamiltonian with a continuous U (1) symmetry which is broken by the periodic drive.We consider two classes of protocol, one where only the boundary spin is driven while the multisite drive also includes a subset of spins in the bulk.Through a second-order F-M expansion of the Floquet Hamiltonian, we show that the edge spin can undergo a significant slowdown for the boundary drive at certain freezing frequencies, which arise due to the non-local corrections in the Floquet Hamiltonian.In contrast, such higher-order processes were previously reported to destroy the dynamic localization in interacting systems [96].At the freezing frequency, we argue that the system develops an approximate global discrete Z 2 symmetry which is responsible for the boundary slowdown because of the formation of an approximate SZM.We find the boundary relaxation slows down with increasing system size but the relaxation time scale does not appear to grow exponentially.
We propose a multi-site drive protocol that further slows down the boundary relaxation.A special sequence of driven sites leads to the slowdown which scales exponentially with the number of driven sites.The emergence of a SZM is shown through signatures in spectral properties and eigenstate overlaps.The significantly enhanced slowdown of the boundary spin coexists with the rapid thermalization of the bulk degrees of freedom which realizes a novel regime that exhibits thermal and athermal behavior.The intermediate athermal behavior also leaves an imprint in the entanglement structure of the Floquet eigenstates.Although the dynamics of the boundary spin relaxes on a time scale exponentially large in system size, the entanglement of the boundary spin is maximal in the Floquet eigenstates.Along with maximal single-site entanglement, random bipartitions also exhibit volume law scaling.Perhaps not surprisingly the bipartition involving the driven sites is significantly athermal.We leave the complete description of the entanglement structure as an interesting future problem.The eigenstates also show thermal and athermal properties reflecting the non-trivial consequences on the thermalization of the driven system.Our description provides a novel realization of an emergent SZM in a non-integrable system realized through local Floquet engineering.
The stability of local athermal dynamics on a background of a thermalizing environment is of relevance for protecting quantum information.The emergence of quasi-local conservation laws with exponentially long relaxation times could potentially play a role near the many-body localization transition [97,98].The interplay between the local Floquet drive and emergent symmetries can provide a path towards realizing novel states with a rich entanglement structure far from the ground state.A general understanding of the entanglement structure of states which are generically volume-law entangled, yet contain elements of athermal character are potentially important for stabilizing non-equilibrium quantum orders away from the ground state.Our prescription for realizing local athermal dynamics utilizes the non-local multi-spin interactions generated by local Floquet engineering.This can be a valuable tool for engineering models with novel athermal character including QMBS [99].The use of optical tweezers in Rydberg atom arrays opens the opportunity for implementing local drives where the physics of SZM can be investigated [16,100].Majorana zero modes in the ground state have been studied experimentally on solid state platforms [3,101,102].Floquet driving provides another means of preparing such zero modes which can be robust to noise.
The stroboscopic Floquet Hamiltonian in the F-M expansion is given by We calculate this expression analytically up to the 2nd term (l = 2) in the subsections below.
, (A19) , (A20) where we have used the values of the integrals, Values of these integrals substituted in Eq. (A10) directly give the H F,loc .The interesting fact is that even H F,loc contains terms which are O(1/γ).This modifies the freezing conditions which were believed to be true at least in O(1/γ).

d. Full H F (local)
As mentioned before, the resummed H F is a series expansion in 1/γ and for many purposes, we can neglect the higher order (which are at least 1/γ 2 or smaller) terms.Therefore, summing up all the one and two site terms up to 2nd order in F-M expansion which are O(1/γ), we get the local Floquet effective Hamiltonian [Eq.(18) in the main text] F + H For concreteness, we also write the full local H F (considering the full contribution from the l = 2 term in F-M expansion) F + H (2) Hamiltonian is now given by W (t) gives a global rotation of all the terms about the x-axis.Consequently, the Heisenberg part remains intact because of the SU (2) symmetry.The zeroth term (l = 0) in the F-M is given by where λ = γτ /2.We get H F = 0 as before.To calculate the 2nd term (l = 2) in F-M, we find where Thus, using the integrals in Eq. (A25), we get where Note that, unlike the local driving, the l = 2 terms here are at least O(1/γ 2 ) or smaller.Finally, the full H F can be obtained by adding all these contributions : H F .Thus, we see, only the field terms are renormalized but no additional interacting term is generated.Therefore, in this case, the global driving only induces a fully coherent oscillation from the initial x-polarized (ψ 0 ) state at any drive parameter The driven state being a globally rotated version of ψ 0 always remains an eigenstate of the SU (2) symmetric Heisenberg interaction part and hence does not suffer any dephasing.In this case, the field terms can get completely canceled at fine-tuned drive frequencies resulting in true freezing of the wavefunction.The exact numerical data including the shift of the freezing frequencies (from ω k c to ω k m ) can be accurately captured by H F [O(1/γ 3 )] as shown in Fig. 11.

Appendix D: Level statistics of the Floquet Hamiltonian
In this appendix, we study the level statistics of the Floquet Hamiltonian (H F ).At generic drive frequencies (ω ≲ γ) H F hardly has any symmetry and hence we increasingly order the Floquet eigenvalues (ϵ i ) in the entire first Floquet BZ.At the freezing frequencies, D z is an emergent conserved quantity and we confine ourselves in a specific D z sector.We then calculate the gap between consecutive eigenvalues: g i = ϵ i+1 − ϵ i (we are using a different symbol for the gap to distinguish it from the spectral gap defined in the main text) and their ratios defined as The distribution of r i gives important information about the integrability of the system.For example, chaotic (nonintegrable) systems exhibit a GOE distribution, whereas integrable systems admit Poisson level statistics.In our case, we find the distribution to be more like Poissonian at generic drive frequencies (see Fig. 12).Here we note that in addition to integrability (an unlikely possibility in our case), Poissonian statistics can appear in several other situations, for example, if there is some hidden unresolved symmetry or if the Hilbert space is fragmented.We find the following operator to be a good candidate of an approximate discrete symmetry at any frequency F ] = 0 for any λ.Note that M converts to D z when λ = 2πk.We also numerically find that for any driving protocol with every fourth site driven, ⟨Φ r |M|Φ r ⟩ ≈ ±1 on the full range of λ.The existence of such operators, which may remain weakly conserved, is most  Here we corroborate the emergence of SZM (presented in the main text) with more numerical results.If a Floquet system described by the Floquet Hamiltonian H F possesses a SZM with S x 1 as the leading part, then the whole Floquet spectrum can be decomposed into quasi-degenerate eigenpairs |Φ ± r ⟩ (with D z |Φ ± r ⟩ = ±|Φ ± r ⟩ and H F |Φ ± r ⟩ = ϵ ± r |Φ ± r ⟩) where |Φ ± r ⟩ can be obtained from |Φ ∓ r ⟩ under the action of S x 1 up to a global phase.We increasingly order both ϵ + r and ϵ − r s within the first Floquet-BZ, and also define ⟨ξ − r |ξ + r ⟩ = ⟨ ξ− r | ξ+ r ⟩ = R r e iθ r , with R r < 1 being the absolute value and θ r being the phase.Such labeling of eigenpairs may lead to some subtleties in some cases which we clarify now.We find that some Floquet eigenpairs (say, ϵ ± p ), obtained in this way, may yield a super low value of |⟨Φ − p |S x 1 |Φ + p ⟩| instead of ∼ 0.5.Such ghost eigenpairs are found to appear always in pairs (say, p, p+1) and a cross pairing may restore the value of the matrix element i.e. |⟨Φ − p(p+1) |S x 1 |Φ + p+1,(p) ⟩| ∼ 0.5.Such cross-pairing may also reduce one of the pairing gaps but then the other one will definitely increase.This kind of ghost pairs get rarer with the emergence of SZM.But, when one of the Floquet eigenvalues of a normal pair crosses the Floquet BZ boundary, that results in a cascade of mismatch between the increasingly ordered ϵ ± r which generates ghost pairs throughout the spectrum.This leads to a sharp drop in |⟨Φ − |S x 1 |Φ + ⟩| at some ω which are mostly away from the special frequencies (ω k m ) where even the conservation of D z , the first criteria of having a SZM, doesn't satisfy.We discard such rare pathological points from our consideration while plotting Fig. 8(b) in the main text.Thus the correspondence between the increasingly ordered ϵ ± r may get ill-defined at drive frequencies away from the slowest relaxation points (ω k m ).But as we increase the number of optimally chosen driven sites (n 4th d ) and tune the drive frequency towards ω m , the correspondence gradually becomes more and more accurate signaling the emergence of SZM.
We demonstrate this in Fig. 13 where we show the behavior of |⟨Φ − r |S x 1 |Φ + r ⟩| in parallel with the behaviors of R r with increasing n 4th d .One can see that there are many pairs with |⟨Φ − r |S x 1 |Φ + r ⟩| < 0.5 including the ghost pairs with |⟨Φ − r |S x 1 |Φ + r ⟩| ∼ 0 for s d = {1}.The number of such pairs decreases with increasing n 4th d and totally disappears when every fourth site is driven in a chain of fixed length N [as shown in Fig. 13(d,h)].This follows from the consistent behavior of R (as shown in Fig. 13 (e)-(h)) which was conjectured in the main text as a requirement for the emergence of SZM.We note here that the lower value of |⟨Φ − |S x 1 |Φ + ⟩| for s d = {1} compared to the value for other higher n 4th d is not an artifact of just the presence of ghost pairs.There are many pairs with |⟨Φ − |S x 1 |Φ + ⟩| a bit less than 0.5 for direct pairing albeit much higher than the corresponding value for cross pairing.Such pairs also reduce in number and disappear with increasing n 4th d , ramping up the value of |⟨Φ − |S x 1 |Φ + ⟩| towards 0.5.The parallel approach of all the R r towards one enable us to arrive at the condition: S x 1 |Φ ± r ⟩ = 1 2 e ±iθ r |Φ ∓ r ⟩ via the ansatz used in Eq. 22 in the main text.
In Fig. 14 we show the behavior of S vN 1 for all the Floquet eigenstates in parallel with the behaviors of δ ± r = ⟨ξ ± r | ξ± r ⟩.There are lots of low-entangled Floquet eigenstates for s d = {1} but all of them get strongly entangled for s d = {1, 5}.This is accompanied by a huge reduction in the average value of |δ ± |, we denote by |δ ± |.Further increase of driven sites (n 4th d ) is found to reduce |δ ± | only a little bit.Consequently, the entanglement spectrum also becomes only slightly narrower, though by now they are already saturated at their maximal value to a good extent.Thus in the asymptotic limit, when every fourth site is driven, the boundary spin forms Bell like pairing with the rest of the system in all the Floquet eigen-

FIG. 2 .
FIG.2.Integrated level spacing distribution I(s) for H 0 (in blue) and H 0 without staggering (in yellow), both with h = 1, N = 16.Black lines show a comparison with Poisson statistics, the signature of integrability (dashed), and GOE, valid for chaotic systems (dotted).

1 FIG. 3 .
FIG. 3. Stroboscopic dynamics of (a) single site entanglement entropy S vNi for all the sites (i) for N = 14.The red (dashed) line is used to denote the maximum single-site entropy (ln 2).S vN 1 (averaged S vN 1 over first 2000 cycles) is plotted with ω in the inset showing the absence of dynamic freezing of the edge spin even for the most appropriate choice of the drive frequency.A vertical dashed line is also drawn at ω 1 c = 7.5 to clearly show the minima in S vN 1 (ω 1 m ) is shifted from ω 1 c .(b) entanglement entropy of the boundary site (S vN 1 ) for different system sizes (N ).h = 1, γ = 15, ω = ω 1 m = 7.53 for both the panels.The N = 16 and N = 18 calculations were done with matrix-product-state-based TEBD calculations for a time step of δt = τ /50 and a maximum bound dimension of 2 N/2 .

( 2 )
F,nonloc ) with their respective strengths in terms of the drive parameters (see Appendix A).Such terms when included in H F are found to improve the agreement with exact numerics for the boundary spin (see Fig.5(b)).

FIG. 5 .
FIG. 5. Comparison of stroboscopic dynamics from exact numerics and analytical H F obtained from the resummed F-M expansion at (a) ω = 5 (b) ω = 7.53 (c) ω = 10.γ = 15, h = 1, N = 8.Away from the freezing frequency the stroboscopic dynamics is well captured for all spins within a "local" approximation for H F , in comparison to exact numerics.At or near the freezing frequency, non-local terms are required to accurately capture the behavior of the driven boundary spin.

14 FIG. 7 .
FIG. 7. Long time stroboscopic dynamics of S vN1 for different multisite driving showing the optimal choice of drive locations (every fourth site) to freeze the boundary spin.Inset shows S vN 1 vs n for s d = {1, 5} for different system sizes.All parameters are the same as in Fig.6.

FIG. 9 .
FIG. 9. Scaling of averaged pairing gaps (|∆|) (a) with system size (N ) for a fixed number of driven sites and (b) with increasing the size of the driven cluster n 4th d for a fixed system size (N = 14).All other parameters are the same as in Fig. 6.The dashed green line in subfigure (b) is a best-fit curve to an exponential ansatz, |∆| ∼ exp(−2.9n4th d ).

FIG. 10 .
FIG. 10.Entanglement entropy of the boundary site with the rest of the system (S vN 1 ) in the upper panels, and sites {1,5,9,13} with the rest of the system (S vN 1,5,9,13 ) in the lower panels for N = 14.The driven sites are mentioned at the top of each column.All other parameters are the same as in Fig. 6.Page values of entanglement entropies are shown by the red dashed line.The bipartition used to calculate the entanglement entropy in panels (c) and (d) is shown below in subfigure (e).

FIG. 12 .
FIG. 12. Level statistics of the Floquet eigenvalues, (a) at a generic drive frequency ω = 5.6 (b) at the freezing frequency ω = ω 1 m = 7.53.The blue line corresponds to the Wigner surmise (an approximation to the GOE statistics of chaotic systems), while the red line is the Poisson distribution (integrable systems).N = 14, γ = 15.
Appendix E: More results on the emergence of SZM and entanglement of the Floquet eigenstates