Long-lived period-doubled edge modes of interacting and disorder-free Floquet spin chains

Floquet spin chains have been a venue for understanding topological states of matter that are qualitatively different from their static counterparts by, for example, hosting π edge modes that show stable period-doubled dynamics. However the stability of these edge modes to interactions has traditionally required the system to be many-body localized in order to suppress heating. In contrast, here we show that even in the absence of disorder, and in the presence of bulk heating, π edge modes are long lived. Their lifetime is extracted from exact diagonalization and is found to be non-perturbative in the interaction strength. A tunneling estimate for the lifetime is obtained by mapping the stroboscopic time-evolution to dynamics of a single particle in Krylov subspace. In this subspace, the π edge mode manifests as the quasi-stable edge mode of an inhomogeneous Su-Schrieffer-Heeger model whose dimerization vanishes in the bulk of the Krylov chain. Floquet physics occurs when a system is driven periodically in time and can be used to uncover phases absent in thermal equilibrium. Here, the authors theoretically investigate heating in a Floquet spin chain, demonstrating the existence and stability of long-lived π edge modes and determining their lifetime from exact diagonalization and Krylov techniques.

P eriodically driven or Floquet systems are promising venues for identifying states of matter that do not exist in thermal equilibrium 1 . Perhaps the most striking among these athermal states are Floquet systems with new topological properties and bulk boundary correspondences that require defining new topological invariants 2,3 . For example, while a static twoband system may be characterized by an integer Z valued or a Z 2 valued topological invariant [4][5][6] , two-band free fermion Floquet systems require at least one more topological invariant so that the system has a Z Z or a Z 2 Z 2 classification [7][8][9] . The additional topological invariant arises due to the periodic nature of the Floquet spectrum. In particular, since the energy in a periodically driven system is conserved only up to integer multiples of the drive frequency Ω, one has quasi-energy bands rather than energy bands. These quasi-energy bands span a Floquet Brillouin zone (FBZ), ϵ 2 ÀΩ=2; Ω=2 Â Ã , with the boundaries of the FBZ being continuous. The additional topological invariant characterizes edge modes that can reside at the zone boundaries ϵ = ±Ω/2.
For two dimensional systems, the additional topological invariant leads to the anomalous Floquet insulator where a bulk has zero Chern number, yet chiral edge modes propagate at the boundary [10][11][12][13] . For one dimensional (1d) systems, the edge modes that reside at the Floquet zone boundary are known as π edge modes [14][15][16] . Since these modes have a periodicity which is twice that of the drive period, they are examples of boundary time crystals [17][18][19] . Of course, the Z Z or a Z 2 Z 2 classification implies that 0 and π edge modes can exist separately, or together, leading to a rich phase diagram [20][21][22][23] .
The π edge modes occurring in most 1d systems [14][15][16]24 also have their origin in what are known as strong π modes (SPM) 25 , whose precise definition we now give. Denoting the π strong mode as Ψ π , these are operators that have the property that they anti-commute with a discrete, say Z 2 symmetry, which we denote by D. Thus, fΨ π ; Dg ¼ 0. In addition, the SPM anti-commutes with the Floquet unitary U that generates stroboscopic timeevolution, {Ψ π , U} ≈ 0. The symbol ≈ is to represent the fact that the anti-commutation is strictly speaking obeyed only in the thermodynamic limit. For a finite size system of length L, {Ψ π , U} ∝ |u| L , where |u| < 1, so that the anti-commutator is suppressed exponentially in the system size. The third feature of the SPM is that it is a local operator with the property Ψ 2 π ¼ Oð1Þ. The existence of a SPM implies an eigenspectrum phase where each eigenstate n j i of a certain parity has a pair Ψ π n j i of the opposite parity, with the quasi-energies of the two states separated by π/T, where T = 2π/Ω is the period of the drive. Note that the above definition of SPMs can be generalized to more complex 2π/k edge modes where k is an integer 26 .
When interactions are included, these operators no longer exactly anti-commute with U in the thermodynamic limit, and therefore acquire a lifetime. Yet a fascinating aspect of these edge modes, which is the central topic of this paper, is that the lifetime far exceeds bulk heating times 25 . These quasi-stable edge modes that almost, rather than exactly anti-commute with U, are referred to as almost strong π modes (ASPM). When disorder is present, the expectation is that many-body localization will make the SPM stable to interactions at all times 20,[27][28][29][30][31] . This paper, in contrast, concerns the stability of π modes for disorder-free chains, and determines how their lifetime depends on interactions. There are of course other contexts, besides ASPMs, where disorder-free Floquet systems can be stable to interactions for long times [32][33][34][35][36][37] .
An analogous definition exists for a strong zero mode (SZM) Ψ 0 38-41 . This is a local operator that obeys all the above properties except that Ψ 0 commutes with the Floquet unitary in the thermodynamic limit Ψ 0 ; U Â Ã % 0. Thus existence of a SZM implies an eigenspectrum phase where the entire spectrum of U is at least doubly degenerate. Adding interactions makes the SZM into a quasi-stable almost strong zero mode (ASZM). Since SZMs and ASZMs have by now been discussed in detail in static systems [42][43][44][45][46][47] , in this paper we will only focus on SPMs and ASPMs. A central goal of this paper is to establish how the lifetime of the π mode depends on interactions. In the process we present Krylov timeevolution as a tool for studying Floquet dynamics. This approach has so far only been employed for static Hamiltonians [46][47][48][49][50] . Here we show how Krylov techniques can be used to extract tunneling times of quasi-stable modes of driven systems.
We will study the stroboscopic dynamics of an open spin-1/2 chain of length L. The dimension of the Hilbert space of the problem is 2 L . The natural numerical method of choice is exact diagonalization (ED), and we are limited to a system size of L = 14. Since we are interested in extracting lifetimes in the thermodynamic limit of L → ∞, ED can be restrictive, and alternate numerical and analytical methods are needed. Thus the results presented here, besides being an explication of the unusual physics of ASPMs, also provide a roadmap for developing methods for understanding slow, system size independent dynamics.

Results and discussion
Model. We will study an open chain of length L whose stroboscopic time-evolution is generated by the following Floquet unitary where Above σ x,y,z are the Pauli matrices, and in what follows we set J x = 1. Equation (1) is an example of a ternary drive where during one period, a magnetic field of strength Tg is applied, followed by the application of a nearest neighbor exchange interaction of strength T along the x direction, and this is followed by the application of a nearest neighbor exchange interaction of strength TJ z along the z-direction. The Floquet unitary has a discrete symmetry as it commutes with D ¼ σ z 1 σ z L . When J z = 0, the model can be mapped to free fermions through a Jordan-Wigner transformation. In this free limit any operator can be expanded in a basis of 2L Majorana fermions or Pauli string operators. Thus when J z = 0, the dimension of the space needed to diagonalize the problem is only 2L rather than the dimension of the full Hilbert space 2 L . The existence of this free limit has been valuable for identifying strong modes. In fact, at the special point of Tg = π, the Floquet unitary for J z = 0 reduces to U ¼ ðÀiÞ L De Ài T 2 J x H xx 20 . At this special point, the SPM is simply the Pauli spin operator on the first site Ψ π ¼ σ x 1 . It is straightforward to check that σ x 1 has all the properties of a SPM as summarized in the Introduction. It is a local operator with Ψ 2 π ¼ 1, it anti-commutes with D, and it anti-commutes with U. Even away from this special point, although the SPM is a more complex operator, it continues to have the property that it has a non-zero overlap with σ x 1 15,25 . Including interactions, i.e, J z ≠ 0, as anticipated in the Introduction, the SPM changes to an ASPM. Moreover, this operator continues to have an overlap with σ x 1 for weak interactions. Therefore, an efficient way to determine whether the system hosts these special edge modes is through the study of the following autocorrelation function We only study the dynamics at stroboscopic times, so n is an integer. Moreover, A ∞ is an infinite temperature average as the trace is over the entire Hilbert space. Thus this quantity is a highly out of equilibrium measure of the system dynamics. When an ASPM or an ASZM is non-existent or if we replace σ x 1 by an operator deep in the bulk of the chain, A ∞ will decay to zero within a few drive cycles. Throughout this paper, all numerical results will be for g = 0.3 and T = 8.25, but for different values of J z . For this value of g, T, when J z = 0 we have a SPM. When J z is non-zero, the SPM changes to an ASPM.
Exact diagonalization (ED). We first present results for the stroboscopic time-evolution of A ∞ from ED. Figure 1 plots A ∞ for J z = 0.01, for three different system sizes L = 6, 8, 10. The x-axis which denotes the stroboscopic times is on a logarithmic scale. One finds that A ∞ flips sign between neighboring stroboscopic times, thus we have an ASPM. Moreover, for small system sizes, as L increases, the lifetime of the ASPM increases (compare L = 6 with L = 8 in the figure), but eventually the lifetime reaches a system size independent value. For the chosen parameters, the system size independent lifetime is reached by L = 8 as the plots for L = 8 and L = 10 nearly coincide. We refer to this system size independent lifetime as the lifetime in the thermodynamic limit. Moreover, for J z = 0.01, this lifetime is around 10 4 drive cycles. Note that the high frequency limit requires T ≪ 1, and since we have T = 8.25, we are far from the high frequency limit. Thus, the bulk is in fact heating, as expected for a periodically driven interacting, and disorder-free system [51][52][53][54][55] . The evidence from ED for bulk heating was already presented in ref. 25 , where it was shown that the autocorrelation function for bulk operators decays to zero within a few drive cycles, and that the entropy density rapidly approaches the maximum possible value (accounting for finite size effects). Further below we discuss the signature of bulk heating in Krylov subspace. We also extend the results of ref. 25 by extracting the interaction dependence of the lifetime of the πmode using several different approaches: ED, Krylov dynamics, and domain wall counting.
We now discuss how the lifetime in the thermodynamic limit depends on J z . Figure 2 shows the autocorrelation function accounting for its period-doubled behavior A p 1 ðnTÞ ¼ ðÀ1Þ n A 1 ðnTÞ. Since the sign-flipping of A ∞ is absorbed by the (−1) n factor, A p 1 has a smoother behavior in time. We plot A p 1 for different J z and for system sizes L = 8, 10, 12, 14. The smallest possible J z we can study is J z = 0.001 because for J z values smaller than this, the system size needed for the lifetime to become L independent, is larger than L = 14. While the stroboscopic times in Fig. 1 were linearly separated, the stroboscopic times in Fig. 2 are logarithmically separated as the lifetimes increase dramatically with decreasing J z , and linearly separated points are numerically too costly to compute. Figure 2 clearly shows that as J z increases (panels (a-g)), the lifetime of the ASPM decreases. In fact the thermodynamic lifetimes are already reached for L = 8 when J z = 0.01 as the plots  The autocorrelation function accounting for the sign-change between neighboring strobosocopic times A p 1 ðnTÞ ¼ ðÀ1Þ n AðnTÞ is plotted for stroboscopic and logarithmically separated times nT, where T is the period and n is an integer. The results are for a magnetic field of strength g = 0.3, period T = 8.25, and for different interaction strengths J z . The xaxis is on a logarithmic scale. From panels a-g J z increases. Panel g corresponds to J z = 0.01 and the corresponding time evolution for A ∞ (nT) for linearly separated times is shown in Fig. 1. Also shown is a fit (purple line) to an exponential A p 1 ðnTÞ ¼ e ÀΓnT A p 1 ð10TÞ, where Γ is a decay rate.
for all the four system sizes lie on top of each other (panel (g)).
Recall that Fig. 1 is a more detailed version of A ∞ for precisely this value of J z , where the stroboscopic times are linearly separated, and one smaller system size, L = 6, is shown in order to highlight the system size dependence. We employ the ansatz that A p 1 ðnTÞ ¼ e ÀΓnT A p 1 ðn 0 TÞ, where we choose n 0 = 10 as it takes about 10 drive cycles for the initial transients to decay. The decay rate Γ is extracted from determining the time at which A p 1 ðΓ À1 Þ % e À1 A p 1 ð10TÞ. This ansatz is plotted in Fig. 2 and captures the time-evolution very well. The decay rates Γ are plotted in Fig. 3 where now it is the yaxis that is plotted on a logarithmic scale. The almost linear slope for 1/J z ≫ 1 suggests that (restoring J x ) Above c is a O(1) number that depends on g, J x . Thus for small enough J z , ED indicates that the decay rates are non-perturbative in the strength of the interactions J z . As J z is further increased, we do not expect the edge mode to be an ASPM, and its decay rate will be determined entirely by perturbative processes Γ / OðJ α z Þ, where α is a power of O(1).
Krylov subspace dynamics. In order to understand the origin of the non-perturbatively long lifetimes for J z ≪ 1, and its relation to bulk heating, we map the stroboscopic time-evolution of σ x 1 to dynamics of a free particle in a Krylov subspace employing a recursive Lanczos scheme [48][49][50]56,57 . This mapping to single particle physics will allow us to develop a tunneling picture for the lifetime of the ASPM.
Let us define the Floquet Hamiltonian as H F ¼ i lnðUÞ=T. The stroboscopic time-evolution after m periods of a Hermitian operator O can be written in terms of H F as follows where we define LO ¼ ½H F ; O. To employ the Lanczos algorithm, we recast the operator dynamics into vector dynamics by defining O j Þ ¼ O. Since we are concerned with infinite temperature quantities, we have an unambiguous choice for an inner product on the level of the operators, ðAjBÞ ¼ Tr A y B Â Ã =2 L . The Lanczos algorithm iteratively finds the operator basis that tri-diagonalizes L.
We begin with the seed "state", . The recursive definition for the basis . It is straightforward algebra to check that the above procedure will yield a L of the form The basis spanned by O n Á lies within the Krylov subspace of L and O 1 Á . We refer to this tri-diagonal matrix as the Krylov Hamiltonian H K , H K ¼ ∑ n b n ðc y n c nþ1 þ c y nþ1 c n Þ, and the 1d lattice it represents, as the Krylov chain.
For free systems, the operation L O n Á can be efficiently solved in the Majorana basis. If the starting operator is a single Majorana, then the dimension of the Krylov subspace of that operator will scale as 2L, as free system dynamics can only mix the individual Majoranas among themselves. Outside of free problems, the size of the full set of |O n ) will be large. For example, a system of size L will have~2 2L possible basis operators. Since we are interested in the thermodynamic limit for the lifetime of the edge operator, in what follows, we will treat the Krylov chain to be a semi-infinite chain. The Krylov chain of interest to us is the one where the seed operator Thus the dynamics of A ∞ has been transformed into that of a semiinfinite, single-particle problem where A ∞ (nT) is now the probability that a particle initially localized at the end of the Krylov chain, stays localized at the end at time nT.
As a point of orientation, let us discuss the details of the Krylov subspace in the free limit. In the Majorana basis, the stroboscopic time evolution of an operator a ! ¼ a 1 ; a 2 ; a 2L Â Ã T is (see where K is a 2L × 2L orthogonal matrix and the a i are Majoranas with a 1 ¼ σ x 1 . For studying SM dynamics, our seed operator is a ! ¼ a 1 ; 0; 0; 0 Â Ã T . The components of K can be determined analytically. On comparing Eqs. (5) and (7), we identify the operator iTL with ln K. Since L is an operator, whose precise form depends on the basis, we have argued that the Krylov Hamiltonian H K is related to i ln K by a simple basis rotation. The form of i ln K becomes particularly simple close to the exactly solvable point gT = π and in the high frequency limit T ≪ 1. Denoting s 1 ¼ sinðgTÞ, in the first order in s 1 and T we find (see Supplementary Note 1) The analytic result in Eq. (8) shows that i ln K is like a Su-Schrieffer-Heeger (SSH) model 58,59 with a topologically non- Fig. 3 Decay rates of the π mode extracted from exact diagonalization (ED). The decay rate from the ED data in Fig. 2 plotted on a logarithmic scale vs 1/J z where J z is the interaction strength. The constant slope for 1000 ≥ 1/J z ≥ 500 indicates that for 1/J z ≫ 1, the decay rate is Γ % e Àc=J z , where J z is in units of J x and c is a O(1) number.
trivial dimerization for |s 1 | < T. Moreover the overall shift of π ensures that the edge mode of the SSH model is pinned at π rather than at zero energy. The SSH model is a band insulator with a band-gap controlled by the strength of the dimerization ||s 1 | − T|. In contrast, when the dimerization is zero, the model is a trivial metal. We will see below that switching on interactions leads to inhomogeneities such that insulating regions of non-zero dimerization coexist with metallic regions of zero dimerization.
In the limit where Eq. (8) is valid, we can derive the Krylov Hamiltonian analytically (see Supplementary Note 1). We find that b odd = |s 1 |, b even = T, with a constant term ± π along the diagonals. Thus when |s 1 | < T, the Krylov Hamiltonian is a topologically non-trivial SSH model that hosts a zero mode. The constant term along the diagonal shifts its energy to π.
The b n s for g = 0.3 and T = 8.25, and for the free case J z = 0 are shown in Fig. 4a. For this case, |s 1 | and T are no longer small. Thus there are differences in the Krylov parameters between this case, and the exact solution around gT ≈ π just discussed. One is that the Krylov Hamiltonian has zeros on the diagonals away from the exactly solvable limit. The second is that the hopping on the very first site is large. However, as suggested by the analytic form in the exactly solvable limit, the Krylov chain is a SSH model with a uniform dimerization after n ≿ 4. Figure 4b shows the corresponding spectrum of the Krylov chain, where the three horizontal red lines are at ϵ = 0, ±π. Modes at ±π that are also separated from the bulk modes by a gap, are clearly visible. Thus we see that even though the diagonal term of the Krylov chain is zero, it is the initial large hopping of b 1 ≈ π that ensures that the edge modes of the SSH model are pinned at ±π.
In fact the effective model for the Krylov chain for J z = 0 can be written as H K = H SSH + H E , where H SSH represents a SSH model and captures the behavior from sites n ≿ 4, while H E is an edge Hamiltonian that captures the physics on the first few sites. The SPM ψ π is a zero mode of H SSH , while H E ψ π = πψ π . Thus H K ψ π = πψ π . To obtain a π edge mode, the parameters of H E are finely tuned, while H SSH only requires its dimerization to be topologically non-trivial to ensure a zero mode. The role of H E is to raise the energy of the zero mode to π.
Note that this mapping from the Floquet unitary U to the Krylov chain has lost information about the periodic nature of the spectrum of U, and this manifests as finely tuned b n at the edge of the Krylov chain when a π mode exists. Nevertheless this mapping to an effectively free model helps to arrive at a tunneling estimate for the lifetime of the π mode when interactions are nonzero. We discuss this below.
We now switch on interactions. We expect the dynamics to explore larger regions of the Hilbert space, resulting in more complicated b n . These are shown in Fig. 5b-e for system size L = 12 and with different J z . Figure 5g-j show the corresponding spectra. For easy comparison between the free and interacting cases, Fig. 5a, f correspond to J z = 0. The blue circles (all panels Fig. 5) correspond to carrying out the Lanczos procedure in the spin basis, which is the natural choice when interactions are present. In contrast, the free case involved performing Lanczos in the Majorana basis Fig. 5a, f. The periodicity of U is lost in the Lanczos approach, and the resulting b n are sensitive to the choice of the branch of lnðUÞ. This leads to b n s (blue circles Fig. 5b-e), which do not bear much of a resemblance to the b n s of the free case Fig. 5a, making them harder to interpret. In particular, the free b n s have a perfectly dimerized form for n > 3, and therefore a periodicity of 2. In contrast, the b n s shown by the blue circles (Fig. 5b-e) have a longer periodicity, close to 3. The spectra for the spin basis are shown in Fig. 5g-j (blue circles). These spectra are not restricted to the FBZ. In addition, the periodicity of 3 manifests as 3 gaps for the spectra shown by the blue circles ( Fig. 5g-j), with the gaps located at ±π, 0. These gaps are most clearly visible for the smallest J z = 0.001 (Fig. 5g). In contrast, the spectra of the free b n s (Fig. 5f) have only two gaps. The additional gap in the spectra shown by the blue circles arises because the system has lost information that the quasi-energy spectra are continuous with π being the same as −π.
One may map the dynamics to an alternate Krylov subspace using an Arnoldi iteration scheme 60 that works directly with the Floquet unitary, rather than its logarithm, and therefore bypasses some of the ambiguities of the Lanczos iteration. Alternatively, below we devise a scheme that can extract the relevant physics from Lanczos by a suitable gauge choice.
Since the spectra are periodic, a physically more suitable gauge choice for the Krylov Hamiltonian is the one where the spectra are folded back to lie within the FBZ (orange crosses in Fig. 5g-j). This folding requires transforming the Krylov Hamiltonian H K ! U KεFBZ U y K , whereε FBZ is a diagonal matrix where all the energies lie in the FBZ, and U K is the unitary matrix that diagonalizes H K before the folding. The folding procedure is nonlocal, and therefore gives a new Hamiltonian, which is no longer tri-diagonal. Therefore, a second Lanczos iteration is carried out to recover the tri-diagonal form, resulting in a new set of b n s that are shown by orange crosses in Fig. 5b-e. After this transformation, the new b n s bear a closer resemblance to the b n s of the free case, thus making them easier to interpret.
In comparison to the free case, one notices that a dimerization persists even with interactions, but is non-uniform, and gradually decreases into the bulk of the chain. This is visible in both gauges, i.e., blue circles and orange stars in Fig. 5a-e. The larger J z is, the more rapidly the dimerization decays into the chain. The contrast is most visible between J z = 0.001 (Fig. 5b) and J z = 0.05 (Fig. 5e). The region of the chain where there is no dimerization, represents a metallic state. Thus we have a spatially inhomogeneous system in the presence of J z where a disordered insulating region (represented by spatially fluctuating but non-zero dimerization) is separated from a metallic bulk. An operator with weight in the metallic bulk will spread rapidly, and its autocorrelation function will decay to zero within a few drive cycles. The existence of the metal is the signature of bulk heating because the metal has no localized states. The emergence of the metal is especially clear after the folding procedure where the gaps at zero quasi-energy begin to fill up after the folding, compare folded spectra represented by orange stars with the unfolded spectra represented by blue circles in Fig. 5g-j. Recall that for the free case the dimerization exists throughout the bulk. Thus the structure of the b n s for J z ≠ 0 gives us evidence that a quasi-localized edge mode can exist despite bulk heating.
The above picture also clarifies how the ASPM acquires a lifetime. Essentially the edge mode is localized initially at the left end of the chain, and is separated by a finite region of dimerization from the metallic bulk. Therefore, it has a non-zero probability of tunneling into the metallic region. Below we estimate the lifetime of the ASPM by determining this tunneling probability.
In order to make the discussion more quantitative, in Fig. 6a we plot the dimerization, i.e., absolute value of the nearest-neighbor b n s, M(n) = b n+1 − b n , corresponding to the data represented by the orange stars of Fig. 5. We plot this quantity after performing a moving average over 4 sites, and denote it as 〈|M(n)|〉 4 . (See Supplementary Note 2 and Supplementary Fig. 1 for the data without the averaging, and with only 2-site averaging for comparison.) We note that 〈|M(n)|〉 4 does not change with J z for the first couple of sites (provided J z < 0.05), while away from the edge, 〈|M(n)|〉 4 decreases with n when J z ≠ 0. In contrast, 〈|M(n)|〉 4 stays constant for the free case. The fact that the first few sites of the Krylov chain do not change with J z implies that H E is not sensitive to J z . We therefore adopt a simple model of a Krylov chain for the sites n ≿ 4, with two slowly varying parameters, a nearest-neighbor average hopping (b n + b n−1 )/2, and the dimerization M(n) 46,47 (see  Supplementary Note 3).
We now emphasize some important differences between the b n s for the ASZM in static systems 46,47 and the same for the ASPM for Floquet systems. One is the sensitivity to gauge choice for the latter due to the freedom in shifting the quasi-energies by integer multiples of 2π (see detailed discussion above). The second difference is that when the Floquet spectrum is bounded by the FBZ, the average b n s do not increase unboundedly with n, unlike in static systems. Thus we derive a continuum model under the assumption that the nearest-neighbor average hopping is spatially uniform, and that the dimerization is slowly varying in space. These assumptions map the Krylov chain onto a Dirac model with a spatially inhomogeneous mass (see Supplementary Note 3) i∂ tΨ ¼ mðXÞσ y þ σ x i∂ X Â Ã Ψ, where m(n) = M(2n). For J z = 0, the mass is spatially uniform and topologically non-trivial, m > 0, with M(2n) = − M(2n + 1) = m. For J z ≠ 0, this mass vanishes into the bulk. While the precise model for how it vanishes is complicated, we adopt a simple ansatz where m(X) = M 0 θ(X − X 0 ). Then a WKB treatment shows that the lifetime of the edge mode is (see Supplementary Note 3) Γ % 4M 0 e À2M 0 X 0 . The fact that the edge mode is at π energy rather than at zero energy enters in the boundary condition via H E , where a strong local hopping pins the edge mode to π.
We now discuss the J z dependence of M 0 . Figure 6b shows that hjMðn ¼ 24Þji 4 $ J À1 z . Since the decay-rate Γ depends on the mass M 0 exponentially, and M 0 ∝ 1/J z , we conclude that Γ $ e Àc=J z .
Bound for lifetime from domain wall counting argument. We now present an alternate argument for the non-perturbatively Fig. 5 Interaction dependence of the parameters of the Krylov chain. The Krylov chain hopping parameters b n (panels a-e) and the corresponding spectra (panels f-j), where n labels the sites of the Krylov chain. The results are for a magnetic field of strength g = 0.3, period T = 8.25, and for the interaction J z increasing from panels a-e and panels f-j, with the J z = 0 case shown in panels a, f. The system size for J z = 0 is L = 20 (same as Fig. 4), while that for the interacting cases is L = 12. Blue circles are for the unfolded spectra, while orange crosses for non-zero J z are for a gauge choice where the spectra have been folded back into the Floquet Brillouin zone (FBZ). The folded spectra (orange crosses (g-j)), and the corresponding b n (orange crosses (b-e)) bear a closer resemblance to the free problem (a, f). The horizontal red lines in panels f-j indicate energies at 0, ±π. long lifetime in Eq. (4). We show below that despite the low frequency driving, the energy required to flip the spin on the very first site is highly off-resonant, and requires rearranging many domain walls in the bulk. This phenomena was also noted in previous studies on static systems 42,45,46,61 . Thus the boundary is in a prethermal state [62][63][64][65][66] , despite the thermalized bulk. We argue this physics by deriving a Floquet Hamiltonian H F in the limit of J z ≪ 1 and |gT − π| ≪ 1. Recall gT = π and J z = 0 is the exactly solvable limit where the SPM is Ψ π ¼ σ x 1 . Let us defineĴ z ¼ J z T=2;ĝ ¼ Tg=2 À π=2 andĴ x ¼ TJ x =2. We will work in the limitĴ x ) 1 andĴ z ;ĝ ( 1. We cannot perform a high-frequency expansion 67 to construct H F asĴ x is not small. Nevertheless, H F to first order inĴ z ;ĝ but to arbitrary orders inĴ x may be derived from an infinite resummation of the Baker-Campbell-Hausdorff formula 68,69 , leading to the following nonlocal perturbed Ising model 25 Above where h E denotes the contributions from the edge spins σ α 1;L and h B denotes the bulk spins. Let us consider the energetics involved in flipping σ x 1 for g ¼Ĵ z ¼ 0. Since the boundary spin has only one neighbor, this flip costs energyĴ x . However, the energy cost for creating a domain wall in the bulk is 2Ĵ x due to the two neighboring sites. Thus there is a mismatch ofĴ x ) 1 between a bulk and an edge excitation, making the flipping of an edge spin impossible. However this simple argument does not account for processes that can make the domain wall hop from site to site, resulting in a lowering of its energy. Thus, we have to revisit the energy argument by accounting for the kinetic energy of the domain walls.
In order to develop our argument, we first note that H xx counts the number of domain-walls N ¼ ∑ i σ x i σ x iþ1 . Therefore, we write H F as a part that commutes with the number of domain walls, and a part that changes the number of domain walls. In particular, . Both D and V can be written in the form similar to Eq. (9), i.e., as sums over local strings of Pauli matrices. The precise forms of D and V are not needed for our qualitative argument. We only assume that the operator norms of each of the local terms of D and V are Oðĝ;Ĵ z Þ. If these assumptions are valid we can apply the arguments of refs. 42,45,46,61 to find a lower bound on the lifetime of the ASPM. Here, we briefly summarize the physical picture.
First, consider the HamiltonianĴ x N þ D. The spectrum of thê J x -term are states that are separated by multiples ofĴ x because N counts domain walls. The D-term causes the domain walls to move without changing their number. Diagonalization ofĴ x N þ D results in domain wall "bands" with a typical bandwidth ϵ $ jjDjj $ Oðĝ;Ĵ z Þ which is much smaller than the separationĴ x between the bands.
The V-term of the total Hamiltonian does change the number of domain walls, but only by an even number due to the parity symmetry of the total Hamiltonian. A single application of V therefore changes the energy by about 2Ĵ x and is off resonant with the costĴ x of flipping the boundary spin. It is impossible to absorb the energyĴ x within few orders of perturbation theory in V. However, the creation and annihilation of a pair of domain walls would lead to the change of the energy by the order of the bandwidth ϵ (Ĵ x . Therefore, we estimate that one needs of the order ofĴ x =ϵ powers of V to offset the energyĴ x required to flip a boundary spin. The probability corresponding to the required order of perturbation theory goes as jjVjj ½ Ĵ x =ϵ $ jjVjj ½ J x =OðJ z ;gÞ where ||V|| denotes the typical size of the matrix element that creates a domain wall. The above expression is a lower bound for the lifetime. For example, in the two integrable limits (which is a property of the exact rather than the approximate H F ) J z → 0, g ≠ 0 and g → 0, J z ≠ 0, the lifetime should diverge. Empirically combining this observation with the rough estimate above we expect 1/ϵ = O(1/J z , 1/g). When J z ≪ g this empirical formula replaces ϵ → J z in the above estimate making it consistent with Eq. (4).

Conclusions
ASPMs are fascinating objects which have lifetimes that far exceed bulk heating times. Besides presenting evidence for this, we developed a method for extracting their lifetimes by mapping their dynamics to single-particle quantum mechanics in Krylov subspace. While we studied the lifetime for J z ≪ g, determining the lifetime when g, J z are comparable is left for future studies. Our Krylov method for determining lifetimes is generalizable to any spatial dimension, to closed and open systems, and to static and driven systems. In addition, the resistance to heating of the π mode is promising for its experimental realization 70 .

Data availability
All relevant data are available from the corresponding author upon reasonable request.