Discrete time-crystalline response stabilized by domain-wall confinement

Discrete time crystals represent a paradigmatic nonequilibrium phase of periodically driven matter. Protecting its emergent spatiotemporal order necessitates a mechanism that hinders the spreading of defects, such as localization of domain walls in disordered quantum spin chains. In this work, we establish the effectiveness of a different mechanism arising in clean spin chains: the confinement of domain walls into ``mesonic'' bound states. We consider translationally invariant quantum Ising chains periodically kicked at arbitrary frequency, and discuss two possible routes to domain-wall confinement: longitudinal fields and interactions beyond nearest neighbors. We study the impact of confinement on the order parameter evolution by constructing domain-wall-conserving effective Hamiltonians and analyzing the resulting dynamics of domain walls. On the one hand, we show that for arbitrary driving frequency the symmetry-breaking-induced confining potential gets effectively averaged out by the drive, leading to deconfined dynamics. On the other hand, we rigorously prove that increasing the range $R$ of spin-spin interactions $J_{i,j}$ beyond nearest neighbors enhances the order-parameter lifetime \textit{exponentially} in $R$. Our theory predictions are corroborated by a combination of exact and matrix-product-state simulations for finite and infinite chains, respectively. The long-lived stability of spatiotemporal order identified in this work does not rely on Floquet prethermalization nor on eigenstate order, but rather on the nonperturbative origin of vacuum-decay processes. We point out the experimental relevance of this new mechanism for stabilizing a long-lived time-crystalline response in Rydberg-dressed spin chains.


I. INTRODUCTION
Because of the subtle role played by the temporal dimension, spontaneous breaking of time-translational symmetry has long escaped conclusive theoretical formulations [1][2][3]. A meaningful characterization of an extended many-body system as a time crystal requires robust stationary macroscopic oscillations, without a net exchange of energy with external devices [2]. A leap forward has been taken with the realization that certain nonequilibrium setups [4][5][6][7][8][9] allow to circumvent the obstacles posed by thermal equilibrium [10,11]. Discrete time crystals (DTCs) formed by interacting periodically driven quantum spin systems currently represent the theoretical paradigm of large-scale spatiotemporal ordering [2,12]. Signatures of their stable and robust subharmonic response to the drive have been experimentally observed with several state-of-the-art experimental platforms [13][14][15][16][17][18][19].
A major challenge to realizing a DTC is the fact that the external drive tends to repeatedly inject excitation energy into the system, and the resulting heating generally deteriorates large-scale spatiotemporal ordering.
Protecting order against melting necessitates a mechanism to keep the impact of dynamically generated excitations under control and thus prevent indefinite entropy growth. To date, many-body localization (MBL) [20][21][22][23] represents the single robust mechanism to stabilize a persistent subharmonic DTC response: The strong quenched disorder of a MBL system freezes the motion of local excitations, thereby stabilizing long-range order throughout the many-body spectrum of the system [24,25]. This infinite-time stability can be characterized as eigenstate order. On the other hand, the quest for disorder-free DTCs calls for alternative mechanisms to evade thermalization. The crucial observation that energy absorption from the drive is asymptotically suppressed for large driving frequencies [26,27] allowed Else et al. [7] to formulate the notion of a prethermal DTC as a long-lived (as opposed to permanently stable) dynamical phase exhibiting broken time-translational symmetry. This phase relies on the existence of an effective Hamiltonian governing the transient dynamics in a suitable rotating frame, possessing an emergent Abelian symmetry determined by the driving protocol, and supporting a spontaneous breaking of this symmetry at low enough effective temperature. This condition guarantees that a relevant set of initial states will display quasistationary macroscopic oscillations of the order parameter in the original frame.
In this paper, we explore the efficacy of confinement of excitations to stabilize the DTC spatiotemporal order. Motivated by the main limitations of the theory of prethermal DTC, we perform our analysis at arbitrary driving frequencies and do not require slowly decaying interactions. Our main result is that an extremely longlived DTC response can be stabilized by domain-wall confinement, without relying on a Floquet-prethermal longrange ordered Gibbs ensemble. This occurrence is made possible by the intrinsic slowness of the prethermal dynamics: The meltdown of transient spatiotemporal ordering involves nonlocal processes such as the dynamical generation of unbound domain walls; the stronger domain-wall confinement, the longer such processes take.

II. OVERVIEW OF RESULTS
We provide here an overview of the contents and results of this paper, which is meant as guidance to the reader. Our findings are summarized and illustrated in Figs. 1, 2.
We start (Sec. III) by exactly solving for the orderparameter dynamics of the integrable periodically-kicked transverse-field Ising chain. We establish a possibly long but perturbative decay rate γ ∼ 3 of DTC response, where represents the deviation from a kicking protocol implementing perfect spin flips. Furthermore, we identify the physical mechanism leading to order-parameter meltdown as the spreading of a small density ρ ∼ 2 of dynamically generated reversed spins, as the domain walls delimiting them freely move at velocity v ∼ [ Fig. 1(a)]. In passing, we note that this result clarifies previous contradictory findings on finite-size scaling of the DTC signal in this model [68,69].
Building on this physical intuition, we introduce domain-wall confinement via symmetry-breaking longitudinal fields (Sec. IV), as a mechanism to prevent the spreading of reversed domains. In absence of perfect flips, confinement individually stabilizes both the positively and negatively magnetized states. However, similarly to what happens in MBL and high-frequency driven spin chains [7,8], we show that symmetry-breaking terms are averaged out by the drive, generating deconfined effective dynamics despite domain walls being instantaneously confined at all times. As a result, the orderparameter decay rate is only weakly affected, retaining its perturbative nature γ ∼ 3 [ Fig. 1(b)].
We finally consider extending the range of spin-spin interactions, J i,j = 0 for |i − j| ≤ R, as an alternative route to domain-wall binding, which does not suffer from incompatibility with the (explicit or emergent) Z 2 symmetry (Sec. V). The crucial feature that arises in these systems is the coexistence in the spectrum of both "topologically charged" excitations (kinks and antikinks) and "neutral" confined bound states. During prethermal dynamics, confined excitations only generate vacuum fluctuations, resulting in long-lived coherent oscillations of the order parameter. Dynamical generation of unbound kinks and antikinks triggers the vacuum decay, and hence the decay of DTC spatiotemporal order. We prove that this phenomenon is heavily suppressed by the nonlocal nature of topological excitations. In fact, we rigorously establish an exponential enhancement of the order-parameter lifetime, i.e. γ ∼ 2R+1 , under mild genericity assumptions on the couplings J i,j [ Fig. 1(c)]. In other words, the fastest process leading to DTC order melting occurs at a perturbative order that grows with the interaction range R. Leveraging this result, we finally conjecture that for algebraically decaying interactions J i,j = J/|i − j| α the decay rate is asymptotically suppressed faster than any power of the perturbation, γ ∼ A −1/(α−2) , in the parameter range α > 2 where no prethermal order is possible. (As α approaches 2 from above, however, the lifetime 1/γ is eventually superseded by the heating timescale, see below.) We present numerical simulations that not only confirm our theory predictions, but even point to a much more robust and extended stability than analytically understood.
The consequences of our findings are further discussed in the conclusive Sec. VI: • The stabilization mechanism identified in this work does not rely on Floquet-prethermal finitetemperature order, nor on eigenstate order. Rather, it relies on the long-term metastability of "false vacua", familiar from high energy physics [70,71]: atypical states with finite energy density, that decay through rare macroscopic tunneling phenomena. This idea extends the theory of time crystals beyond previously known mechanisms, circumventing some of their limitations (cf. Fig. 2).
• This work indicates a clear route to observing DTC behavior in a class of quantum simulation platforms of growing importance, namely Rydberg-dressed arrays of neutral atoms [72][73][74]. These systems are characterized by a great degree of control on the interaction range and strength. Compared to conventional prethermal and MBL DTCs, the milder requirements to observe confinement-stabilized DTCs come at the price of reducing the set of initial states  • Finally, as a byproduct on the theory side, our work clarifies a long-standing issue, namely the nature of the apparent anomalous persistence of the order parameter in the quench dynamics of quantum spin chains with algebraically decaying interactions J i,j = J/|i − j| α , α > 2 previously observed in several numerical studies [60,75,76]. Our theory establishes that a strong enhancement of the order parameter lifetime is to be generally expected in the parameter regime 2 < α ∞, and predicts its functional form as a function of the quench magnitude. This has the important consequence, not recognized before, that long-lived nonequilibrium order can be sustained by systems in the "shortrange" regime α > 2, which cannot spontaneously form an ordered state in thermal equilibrium. The meltdown of this metastable order is ultimately triggered at long times by the dynamical generation of deconfined topological excitations (domain walls).

III. ORDER PARAMETER DECAY IN THE KICKED TRANSVERSE-FIELD ISING CHAIN
Our starting point is the standard Ising Hamiltonian where J > 0 is the ferromagnetic coupling, X j , Y j , Z j are the local spin-1/2 Pauli matrices for the jth spin, and periodic boundary conditions (j +L ≡ j) are assumed. The Floquet dynamics is obtained by periodically intertwining the evolution governed by H J with kicks K = j K j at integer times t n = n = 1, 2, . . . . The resulting singlecycle time-evolution (Floquet) operator reads To observe the simplest realization of time-crystalline spatiotemporal order, the system is initially prepared in the fully polarized state with positive magnetization along theẑ direction, namely, one of the two degenerate ground states of H J . This has a product-state form |+ ≡ | · · · ↑↑↑ · · · , where | ↑ (| ↓ ) denotes the eigenvector of the Pauli matrix Z with eigenvalue +1 (−1). The kick K is taken to rotate each spin by an angle π around a transverse axis: [77] In this case, the time evolved state after n kicks |Φ(n) = U n |+ exhibits a sequence of perfect jumps between |+ and the other ground state of H J , namely |− ≡ |· · · ↓↓↓ · · · . The persistent nonvanishing value of the order parameter m(n) = Φ(n)|Z j |Φ(n) (4) in both space and time, being equal to (−1) n , gives rise to a trivial example of a macroscopic subharmonic response. This behavior, however, relies on a fine-tuning of the kick strength, g = π/2. The existence of a nonequilibrium phase of matter exhibiting time-crystal behavior revolves around the stability of this spatiotemporal order to arbitrary (sufficiently weak) Floquet perturbations in the thermodynamic limit L → ∞.
A. Exact decay rate in the thermodynamic limit and finite-size effects The simplest perturbation to the above Floquet protocol consists in performing imperfect spin flips, i.e., with = 0. Since the perfect kick K π/2 = i L P can be factored out of K and is proportional to the global Z 2 -spin flip operator P = j X j , the expectation value of the local order parameter (4) over n periods can be expressed as where we used the properties P Z j P = −Z j and [V J , P ] = 0, while K = exp(i j X j ). Equation (6) expresses the fact that the absolute value of the magnetization evolves as if it were governed by the Floquet operator K V J with kick strength equal to , where the perfect kick has been completely gauged away by switching to a toggling frame of reference, leading to the multiplicative factor with alternating sign (−1) n .
The persistence of time-crystalline order is related to the preservation of a finite absolute value of the local order parameter |m(n)| for large times n → ∞. This question has been addressed in previous works investigating finite-size chains. In particular, the analysis of Ref. [68] led to a positive answer based on finite-size scaling of the order parameter obtained by exact diagonalization (ED) of short chains L 20. This finding contradicts the generic expectations of the absence of long-range order in excited states of one-dimensional clean short-range interacting systems.
Here, by computing the exact dynamics of the magnetization for an infinite chain using the integrability of the model, as detailed in Appendix A, we establish that the order parameter decays exponentially in time as |m(n)| ∼ e −γn . The rate γ is found to be where the quasiparticle spectrum φ p and the Bogoliubov angle ∆ p resulting from the diagonalization of the quadratic Floquet Hamiltonian are defined by the equations cos(φ p ) = cos(2J) cos(2 ) + sin(2J) sin(2 ) cos(p) and cos(∆ p ) = cos(θ p ) cos(p) + sin(θ p ) cos(2 ) sin(p), respectively. The explicit expression of θ p is given in Appendix A. By Taylor expanding the exact result (7) for small perturbations , we find that the rate γ scales as Figure 3 reports the exact evolution of m(n) for increasing values of , for J = 1 (colored data sets). In all cases, the resulting decay rate excellently reproduces the analytical result (7) (dashed black lines). As it is evident in Fig. 3(a), the decay remains quite slow even for moderate values of . Furthermore, in Fig. 3(b)-(d), we compare the thermodynamic-limit evolution with that of finite chains of length L = 10 ÷ 30. The plot illustrates the dramatic impact of finite-size effects. The apparent persistence of the order parameter in finite systems has been detected in previous works [68,69,78]. For the kicked Ising chain [68,69], these strong finite-size effects can be attributed to a large overlap of the initial state with the magnetized ground state of the integrable Floquet Hamiltonian, in agreement with the numerical findings of Ref. [79] and similarly to the static case [80,81].

B. Physical interpretation of the order-parameter decay as domain-wall spreading
The scaling in Eq. (8) with can be understood in intuitive terms, considering that the system has exact quasiparticles which behave as noninteracting fermions. The dynamics in Eq. (6) is equivalent to a quench from the ground state of the classical Ising Hamiltonian (1), evolving with a kicked Ising chain deep in the ferromagnetic phase | | J. In this case, the free fermions can be interpreted as topologically protected excitations, i.e., domain walls (kinks and antikinks), interpolating between the two degenerate magnetized ground states. To the lowest order in the kick strength , the quench creates a small density ρ = O( 2 ) of spin flips, whose constituent pair of domain walls freely spread along the chain with maximum velocity v = O(| |). The domain of reversed spins extending between a kink-antikink pair grows linearly in time until one of them meets another domain wall initially located far away. The decay rate of the order parameter can be predicted in terms of this semiclassical picture [82]: we find Indeed, the exact result in Eq. (7) obtained from the asymptotic expansion of the determinant of a large block Toeplitz matrix [83,84], precisely takes the form of a product of quasiparticle group velocity |∂ p φ p | at momentum |p| and (for small ) number sin 2 (∆ p ) − ln | cos 2 (∆ p )| of excited quasiparticle pairs with momenta (p, −p) in the initial state, averaged over momenta. This substantiates the intuitive interpretation above. The exact result thus quantitatively confirms the intuitive model of order parameter meltdown by spreading domain-wall pairs, illustrated in Fig. 1(a). More importantly, this picture highlights what makes timecrystalline order doomed to melt in one-dimensional systems: Preventing domain-wall pairs from unbounded separation requires certain microscopic mechanisms, such as 6 Anderson (many-body) localization induced by quenched disorder. In the rest of this article, we will investigate under which circumstances domain-wall confinement can provide a robust stabilization mechanism.

IV. DOMAIN-WALL CONFINEMENT AND DECONFINEMENT IN THE KICKED MIXED-FIELD ISING CHAIN
The previous section unambiguously illustrates how domain-wall spreading underlies the time-crystal melting in clean, locally interacting, spin chains. A celebrated proposal to overcome this occurrence and protect long-range order out of equilibrium hinges upon disorderinduced localization: in such case, domain walls behave like particles moving in a random background, and spatial localization can arise from destructive interference, as first foreseen by Anderson [20], even in the presence of many-body interactions [21][22][23]. This basic mechanism of localization-protected order [24,25] has been proposed to stabilize time-crystalline behavior for arbitrarily long times [5,6]. Here we explore a different mechanism to prevent domain-wall spreading, namely domain-wall confinement.
In this section we briefly review the current understanding of this multifaceted phenomenon (Sec. IV A) and extend it to driven systems, exemplified by the spinchain dynamics (2) subject to weak kicks K 1 about an arbitrary tilted axis. Specifically, we use the tools of many-body perturbation theory to reformulate the order parameter-dynamics in terms of the motion of effective domain walls (Secs. IV B) and demonstrate domain-wall confinement (Sec. IV C). Finally, in Sec. IV D, we come back to our main discrete time-crystal problem with kicks close to perfect spin flips K i L P , and understand the order-parameter melting, illustrated in Fig. 1(b).

A. Confinement in quantum spin chains
For the benefit of readers who may not be familiar with domain-wall confinement in quantum spin chains, in this subsection we provide a brief overview.
Particle confinement is a nonperturbative phenomenon arising in certain gauge theories, which consists in the absence of colored asymptotic states: all stable excitations of the theory above the ground state are colorless bound states of elementary particles [85]. An intuitive picture of this phenomenon is given by the formation of a gauge-field string connecting a quark-antiquark pair, the energy cost of which provides an effective confining potential that grows linearly with the spatial separation between the two particles. As a result, the quark and antiquark bind together into composite neutral particles called mesons. When a large physical separation between them is enforced, the potential energy stored in the string becomes sufficient to produce another pair of particles out of the vacuum, which bind with the old particles to form two mesons, making the observation of isolated quarks impossible.
An analogous confinement phenomenon naturally arises for domain walls in quantum spin chains. Its mechanism was proposed by McCoy and Wu in 1978 [48] and later studied in a variety of theoretical [49-52, 86, 87] and experimental [88][89][90] works. The core ingredient here is a first-order quantum phase transition, i.e., the explicit lifting of a spontaneously broken discrete symmetry. In the ferromagnetic quantum Ising chain (H = −J j Z j Z j+1 − g j X j , with |g| < J), this can be simply realized by introducing a longitudinal field −h j Z j , which generates an energy penalty for the reversed magnetic domain separating a pair of domain walls, analogous to a string tension. The energetic cost for separating the pair thus grows proportionally to the distance, giving rise to a linear confining potential that fully suppresses the spreading at arbitrarily high energies, binding the pair of "topologically charged" particles (kink and antikink) into "topologically neutral" bound states, referred to as mesons by analogy with particle physics. In certain quasi-one-dimensional magnetic insulators, similar effective longitudinal fields are provided at a mean-field level by inter-chain interactions; the resulting tower of mesonic excitations (spinon bound states) has been spectacularly observed with inelastic neutron scattering [88][89][90].
Recently, it has been realized that confinement in quantum spin chains and in (1 + 1)-dimensional lattice gauge theories can be generically mapped onto each other [56] via elimination/introduction of matter degrees of freedom exploiting the local constraints posed by gauge invariance [91,92]. This substantiates intuitive pictures of the nonequilibrium dynamics of spin chains in terms of prototypical phenomena in gauge theories, such as vacuum decay [56,63,93,94], string dynamics [54][55][56][57], and string inversions [91].
Dynamical signatures of domain-wall confinement have been recently attracting a growing interest, starting from Ref. [53]. The suppression of domain-wall spreading stabilizes the order parameter even out of equilibrium. This stabilization also applies to the dynamics starting from the "false vacuum" magnetized against the longitudinal field, which is a very atypical highly excited state without domain walls [55,56]. The basic explanation of this phenomenon is that domain-wall pairs excited on top of the false vacuum are also confined into "anti-mesons". Furthermore, domain-wall pair production out of the false vacuum is strongly suppressed, despite being energetically allowed and entropically favorable, as it requires a locally created virtual pair to tunnel across a high-energy barrier. This effect, akin to the Schwinger mechanism in quantum electrodynamics [95], results in an exponentially long lifetime of the order parameter [56,63].

B. Effective domain-wall dynamics for weak kicks
Here we show that analogous considerations on domain-wall confinement carry over to driven quantum Ising chains, under the assumption that the periodic kicks are weak. We consider the spin chain evolving with the Floquet operator in Eq. (2), with the simplest Z 2symmetry-breaking periodic kicks, i.e., a tilted rotation axis away from the x − y plane: We defineˆ = max( , h) to conveniently keep track of the global magnitude of the kick. For clarity, we stress that here, we are considering a weak kick K ,h , with 1, in contrast with the almost perfect spin flip K π/2+ , (compare with Eq. (5)) which is relevant for the stability of DTC.
The first problem that we encounter is related to the breaking of integrability for h = 0: domain walls cease to exist as exact quasiparticles throughout the many-body Floquet spectrum, seemingly obscuring the physical interpretation of the magnetization dynamics. However, we can still make sense of an effective picture of domainwall dynamics for weak perturbations, using the rigorous theory of prethermalization of Ref. [96]. Namely, one can construct a sequence of time-periodic, quasi-local unitary transformations {e iˆ m Sm } m=1,...,p , designed to remove from the time-dependent Hamiltonian all the terms that change the total number of domain walls The resulting transformed Floquet operator U p = e −iˆ p Sp · · · e −iˆ S1 K ,h V J e iˆ S1 · · · e iˆ p Sp (12) conserves D 1 up to terms of order p + 1, i.e., U p , D 1 = O(ˆ p+1 ). The perturbative series generated by this construction represents a nontrivial generalization of the well-known Schrieffer-Wolff transformation for static Hamiltonians [97][98][99][100]. The generators {S m } take the form of sums of local operators, whose number and support size grow proportionally to the perturbative order m.
The construction of the transformation e iS ≤p ≡ e iˆ S1 · · · e iˆ p Sp for arbitrarily large p aims at asymptotically producing an exactly domain-wall-conserving Floquet operator. However, the resulting perturbative series is expected to have divergent (asymptotic) character in general, similarly to the static case [100], suggesting a late-time breakdown of the conservation of D 1 and an eventual heating to infinite-temperature. Nevertheless, the transformed picture still contains extremely useful information on the transient dynamics. The analysis of Ref. [96] shows that, when J is strongly incommensurate with the driving frequency 2π, [101] the breakdown of D 1 conservation -and hence heating-must be extremely slow. In fact, the optimal (with respect to suitable operator norms) truncation order depends on the magnitude of the perturbation as p * ∼ C/ˆ 3+δ , where δ > 0 is arbitrary and C = C(δ) is a constant. The very fact that p * scales up with the smallness ofˆ itself yields a nonperturbatively small truncation error. In turn, this translates into a quasi-exponentially long time window [96][102] within which -for the purpose of computing the dynamics of local observables-the nonequilibrium evolution of an initial state |ψ 0 can be approximated as where U p * is the approximate Floquet operator obtained from U p * by truncating terms beyond the orderˆ p * . By construction, U p * exactly conserves the number of domain walls, U p * , D 1 = 0. The physical consequence of this analysis is that the "bare" (i.e., unperturbed) domain-wall occupation number D j,j+1 on the bond (j, j + 1) acquires a perturbative quasi-local "dressing" e −iS ≤p * D j,j+1 e iS ≤p * for small perturbations. The density of such dressed domain walls remains approximately conserved at least for the long timescale T preth in Eq. (13). Note, however, that this heating timescale only represents a rigorous lower bound, and it is not expected to be tight in general. In the integrable limit h = 0, the underlying algebraic structure of the model produces cancellations to all orders which make this perturbative series convergent, leading to exact dressed domain-wall quasiparticles. In this case heating is completely suppressed beyond the above timescale T preth . As soon as h = 0, instead, these emergent domain walls are expected to eventually decay after T preth (see also Refs. [54,103]). In any case, as long as we deal with dynamical phenomena occurring in the long Floquet-prethermal time window 0 ≤ t ≤ T preth , we can switch to the effective picture where states and observables are transformed via the unitary operator e iS ≤p * , and work therein with the effective domain-wall conserving Floquet operator U p * , according to Eq. (14).
As it results from the above discussion, we can analyze the evolution of the order parameter by switching to the transformed picture, where |ψ 0 = e −iS ≤p * |ψ 0 , Z j = e −iS ≤p * Z j e iS ≤p * , and the approximation, due to truncating U p * to U p * , holds up to the long timescale in Eq. (13). In this transformed picture, the number D 1 of domain walls is an exact quantum number, and the Hilbert space fractures into separate blocks labelled by D 1 .
The perturbative construction introduced in Ref. [96] to prove the theorem leading to Eq. (13) is hardly manageable for explicit low-order computations. In practice, we have found it more convenient to resort to a combination of the replica resummation of the Baker-Campbell-Hausdorff (BCH) expansion of Ref. [104] and a standard static Schrieffer-Wolff transformation (see, e.g., Ref. [100]). In both approaches, the existence of a well-defined expansion requires incommensurability of the coupling J with 2π. This condition is necessary to ensure that the domain-wall number uniquely labels the unperturbed sectors of the Floquet operator, and thus remains a good quantum number throughout the perturbative construction. We remark that while the scheme of Ref. [96] generally produces a time-dependent effective Hamiltonian, the combined BCH resummation and Schrieffer-Wolff transformation produce a static effective Hamiltonian, with presumably similar convergence properties.
In Appendix B, we use the latter approach to derive the expression of S and U to lowest order p = 1, reported here: In Eqs. (16), (17), the operators P ↑,↓ j represent local projection operators onto the "up" and "down" states along z of spin j. The off-diagonal processes in U 1 are given by the terms proportional to , which describe nearestneighbor hopping of domain walls to the left or to the right. Processes which create or annihilate pairs of domain walls have been removed from the evolution operator U 1 through the unitary transformation e iS ≤1 ; in this formulation, these processes only show up in the expression of the transformed initial state |ψ 0 = e −iS ≤1 |ψ 0 .
The general structure of the expansion and further details are discussed in Appendix B. Higher-order terms generate further small corrections in the effective Floquet Hamiltonian. In particular, at order p, one has terms in U p that correspond to at most p displacements of one domain wall or more adjacent domain walls to a neighboring bond. For example, at second order one has diagonal terms, plus nearest-neighbor hopping terms analogous to those in U 1 , plus new terms proportional (pair nearestneighbor hopping), and analogous flipped combinations. These additional terms do not modify the conclusions below. Likewise, the higher-order Schrieffer-Wolff generator S ≤p flips at most p neighboring spins.

C. Domain-wall confinement and order-parameter dynamics
Armed with the rigorously established picture of effective domain-wall dynamics for weak kicks, we now discuss domain-wall confinement and its implications for the evolution of the order parameter. Here, as in the previous subsection, we discuss the Floquet dynamics generated by weak kicks U = K ,h V J . In the next subsection we will go back to studying robustness of the DTC signal Equation (15) describes the evolution of the magnetization in the transformed domain-wall picture. The transformed initial state |+ consists of a low density of order 2 of flipped spins, as, to lowest order, this state is obtained by rotating the spins in |+ by an angle / sin(2J) around a transverse axis (note that each spin flip carries two domain walls). Furthermore, since Z j = Z j + O( ), for the purpose of understanding the nature of the evolution of m(n) (i.e., persistent or decaying), we can drop the correction to Z j . Since the initial state is composed of dilute tight pairs of domain walls, we can enlighten the resulting orderparameter dynamics by studying the two-body problem. The intuitive picture of the evolution of m(n) in terms of the motion of domain-wall pairs becomes asymptotically exact in such a low-density limit. The nature of this evolution (persistent or decaying order) depends on the effective Floquet operator U governing the motion of domain walls. Equation (16) describes domain walls of "mass" 2J, hopping to neighboring lattice bonds with amplitude , and experiencing a confining potential V (r) = 2hr which ties them to a neighbor at a distance r.
The two-body problem is obtained by projecting the effective Hamiltonian H 1 , defined by the exponent of U 1 in Eq. (16), onto the sector spanned by two-particle basis states |j 1 , j 2 , where the integers j 1 < j 2 label the positions of two domain walls along the chain. The resulting two-body Hamiltonian is where the hard-core constraint |j 1 = j 2 ≡ 0 is understood. The second term represents nearest-neighbor hops of the domain walls, whereas the first one acts as a linear confining potential V (r) = 2hr as a function of the distance r = j 2 − j 1 > 0. While for h = 0 one has a continuum of unbound traveling domain walls, as soon as h = 0 the spectrum changes nonperturbatively to an infinite discrete tower of bound states. Due to translational invariance of the initial state |+ , domain-wall pairs are only generated with vanishing center-of-mass momentum K = 0. The relative coordinate wavefunction ψ(r) satisfies a Wannier-Stark equation with a hard-wall boundary condition ψ(0) = 0, yield- ing [56,105] the exact mesonic masses where J is the standard Bessel function. [106] For → 0, one finds the energy levels E = 2h , corresponding to a domain of reversed spins, ψ (r) = δ ,r . Figure 4(a) reports a sketch of the low-energy spectrum of the Floquet Hamiltonian in this limit. For finite /h, the eigenfunctions can still be adiabatically labelled by the integer . Panel (b) reports a selection of mesonic eigenfunctions ψ (r) in the center-of-mass frame, for /h = 1.5. Boundary effects are visible for small 2 /h, whereas for larger the wavefunctions become essentially Wannier-Stark localized orbitals, i.e., rigidly shifted copies of each other.
We can formulate a more intuitive analysis of the twobody dynamics, which will turn out to be fruitful later to analyze the time-crystalline behavior. To this aim, we introduce the canonically conjugated operators Q, P defined by which correspond to the position and the momentum in the center-of-mass frame, i.e., the distance between the two domain walls and their relative momentum; one verifies [Q, P ] = i. [107] Using these variables, the centerof-mass frame Hamiltonian becomes where the domain is Q > 0 and a hard-wall boundary condition at Q = 0 is understood. Classical trajectories are bounded in the Q direction and are translationally invariant away from the boundary Q = 0. Indeed, the Heisenberg evolution equations can be integrated exactly in the bulk Q 2 /h (i.e., neglecting the boundary condition), which gives Equation (24) indicates that the relative momentum P revolves freely around the Brillouin zone, whereas the relative coordinate Q performs stable Bloch oscillations, remaining localized near the corresponding initial condition. In physical terms, the mutual confining potential creates an effective Wannier-Stark ladder which pins the distance between the two domain walls.
As anticipated above, the solution of the two-body problem sheds light on the many-body dynamics. Each of the dilute spin flips in the initial state |+ overlaps significantly with the lightest mesonic wavefunctions. [108] The total magnetization thus exhibits a persistent oscillatory behavior, characterized by multiple frequencies associated with the "masses" of the mesonic bound states. As long as the typical separation between distinct initial domain-wall pairs -the inverse −2 of the density in Eq. (18)-exceeds by far the size 1 + /h of the excited bound states, the evolution of the magnetization m(n) can be described as resulting from the coherent superposition of independent small-amplitude meson oscillations, similarly to the undriven case discussed in Ref. [53]. We note that inelastic meson scattering is expected to trigger asymptotic thermalization; however, the inelastic cross section drops rapidly for small [109], making the transient out-of-equilibrium state extremely long lived [110,111].
Analogous considerations can be made for the dynam- ics starting from the "false vacuum" state |− (or, equivalently, taking h → −h). Two consecutive domain walls are subject to an anticonfining potential, decreasing linearly with their separation. Due to lattice effects, however, the domain walls cannot accelerate arbitrarily towards large distances, as their hopping kinetic energy is bounded. Thus, they form a tower of bound states, formally analogous to the true ground-state excitations (21) discussed above. The dynamics of the order parameter thus follows a similar pattern, exhibiting a stable antimagnetization with small-amplitude multiple-frequency oscillations superimposed due to antimesons.
The ultimate decay of the antimagnetization due to resonant domain-wall pair production is expected to occur at very long times, as the high energy cost 4J of creating two domain walls has to be compensated by a gain 2hr associated with a large ground-state bubble of size r extending between them. To realize this, the two locally created virtual domain walls have to tunnel through a distance r * 2J/h, which suggests a nonperturbative lifetime. The more rigorous estimation of Ref. [96] leads to the quasi-exponentially long time in Eq. (13). This phenomenon is closely related to the Schwinger mechanism in quantum electrodynamics, as explained in Ref. [56].
The confinement scenario for the periodically kicked mixed-field Ising chain has been verified numerically, sim-ulating the Floquet dynamics induced by Z 2 -symmetrybreaking periodic kicks K ,h by means of the infinite time-evolving block decimation (iTEBD) algorithm [112]. Figure 5 reports the behavior of the order parameter as a function of the number of n kicks. In fact, either starting from the vacuum |+ [panel (a)], or from the falsevacuum |− [panel (b)], a small value of the longitudinal component h is sufficient to induce a nonperturbative change in the order parameter dynamics: domain walls get confined into (anti)mesons, thus hindering the melting of the order parameter, which remains finite for exponentially long times (continuous colored lines).

D. Deconfinement by driving and DTC lifetime
In the last section we established that domain-wall confinement induced by a Z 2 -symmetry-breaking kick component stabilizes both the magnetization when quenching from the "true vacuum" and the antimagnetization when quenching from the "false vacuum" (cf. Fig. 5).
We are now ready to come back to our main problem of time-crystalline order, and discuss how generic (non-Z 2symmetric) kick imperfections impact the order parameter lifetime determined in Sec. III.
We consider kicks K in Eq. (2) of the form of imperfect spin flips: To make progress, generalizing the approach of Sec. III, we switch to the toggling frame, i.e., rewrite the two-cycle Floquet operator reabsorbing the perfect kick: where, as in Eq. (6), we have exploited the fact that (−i) L K π/2 = P flips the Z axis, leaving the Z 2symmetric interactions V J invariant. Here, however, the toggling frame makes the symmetry-breaking longitudinal component h of the kick periodically flip in sign. The dynamics can thus be seen as generated by strong interactions and weak kicks only (without perfect flips), alternating the sign of the longitudinal component of the kick at each period. The theory of Ref. [96] discussed above for U , straightforwardly applies to U 2 as well: it guarantees the existence of a close-to-identity time-periodic unitary transformation e iS ≤p such that, in the transformed frame, the two-cycle Floquet dynamics described by Eq. (26) approximately conserves the number of domain walls over a long prethermal timescale analogous to Eq. (13). In the previous section we have discussed the transformation e iS ≤p for U = K ,h V J [see Eqs. (12) and (14)]. Unfortunately, there is generally no simple relation between e iS ≤p and e iS ≤p , because the generators S ≤p and the effective Floquet operator U p depend on both and h. In particular, h and −h produce generally different operators S ≤p . This fact prevents us from straightforwardly combining the two transformations for K ,h V J and K ,−h V J into a single one for U 2 .
However, the lowest-order result in Eq. (17) shows that S ≤1 is actually independent of h. This simplification allows us to directly combine the two transformations into a single one for U 2 right away: substituting into Eq. (26), we find that the resulting lowest-order transformation e iS ≤1 coincides with e iS ≤1 in Eq. (17): where U 1 (h) is expressed in Eq. (16). In contrast with U 1 (±h), U 2 1 is not expressed as the exponential of a time-independent local Hamiltonian, but it results from two time steps where the longitudinal field switches between h and −h. On the other hand, the occurrence that the lowest-order transformationS ≤1 = S ≤1 is timeindependent here relies on the special form (25) chosen for the kick perturbation, which makes the present derivation especially simple.
Thus our problem amounts to study the driven dynamics of a dilute gas of domain-wall pairs. To the level of approximation considered above, domain walls are subject to the unitary dynamics expressed by Eq. (16) where, crucially, the confining string tension h regularly flips in sign at integer times, as dictated by Eq. (27). Generalizing the analysis of Sec. IV C, the nature of the evolution of the order parameter -and hence the fate of timecrystal behavior in the presence of confinement-will be essentially determined by the solution of the two-body problem.
Solving the two-body dynamics amounts to composing the evolution map in Eq. (24) with +h and −h. Even though domain walls are completely bound into mesons or antimesons within each individual period, we demonstrate that the periodic switching between the two leads to deconfinement, and thus meltdown of the system magnetization. The exact two-cycle Floquet map restricted to the two-body space in the "bulk" (i.e., for Q 2 /h) is equivalent to the composition of two maps given by Eq. (24) for t = 1, with +h and −h, respectively. The result of this composition is This two-cycle map is equivalent to one generated by an effective static Hamiltonian, which reads as can be readily verified.
Remarkably, H cm is a pure-hopping Hamiltonian, without interaction potentials. Its eigenstates are no longer bound states localized around a finite value of Q, but deconfined plane waves with a definite momentum P . The periodic switching between +h and −h averages out the Z 2 -breaking confining potential, and effectively restores the symmetry, similarly to what happens in highfrequency-driven discrete time crystals [7].
A semiclassical description of the effective domain-wall dynamics is portrayed in Fig. 1(b), where it is highlighted how the periodic switching of the sign of the confining potential h leads to an effective ballistic spreading of the reversed bubble delimited by two domain walls, with a renormalized maximum effective velocity v eff = v eff ( , h). From Eq. (30), we find the approximate dispersion rela- As we neglected all terms in H 2-body beyond the first order in and h [cf. Eq. (19)], we cannot expect the correction to be quantitatively accurate. However, numerical simulations shown in Fig. 6  by a factor which is compatible with a relative decrease of order h 2 for both the density ρ of spin flips in the initial state and the spreading velocity, according to the formula γ ∼ ρv eff . The bottom line of this analysis is that the symmetrybreaking perturbation, leading to confinement of excitations in the static case, does enhance the lifetime of timecrystalline behavior, although only by parametrically decreasing the prefactor a = a(J, h) of the perturbative decay law γ ∼ a 3 . The above arguments can be easily modified to account for generic symmetry-breaking perturbations of the perfect kick. Such perturbations, thus, do not qualitatively modify the picture of the order parameter meltdown found in the integrable case in Eq. (8). Although we derived this result in the lowest perturbative order, we expect that the Z 2 symmetry is restored to all orders, similarly to what happens in MBL and high-frequency driven prethermal DTCs [7,8]. Our simulations of Fig. 6 show no signature of slowdown of the computed exponential decay, clearly confirming the expectation.
We finally remark that, similarly to the integrable case, the order-parameter dynamics is strongly affected by the finite size of the chain: our ED data for L = 10 ÷ 30, shown in Fig. 6 (shaded lines), display a deceptive persistence of the order parameter, as previously observed in different contexts [113], whereas it eventually decays to zero for systems in the thermodynamic limit (thick straight lines). Note that both operators K and V can be exactly applied to the many-body wavefunction, which allows us to efficiently simulate dynamics of 30 spins with a reasonable amount of resources.

E. Enhancement of DTC lifetime at dynamical freezing points
The theoretical analysis above assumed the Ising coupling J = O(1) to be the dominant scale in the Floquet operator, and the transverse and longitudinal fields , h 1 to be weak perturbations to the perfect spin flip. In this Subsection, building on the recent Refs. [114][115][116], we point out that the DTC lifetime can be strongly enhanced if the kick perturbation is tuned to particular large [O(1)] angles, referred to as dynamical freezing points. The argument goes as follows.
We first of all rewrite our kicking protocol in Eq. (25) in the equivalent form considered in Ref. [116]. The twocycle Floquet operator U 2 = P KV P KV can be turned into the form PKV PKV with a modified kick and finally perform a unitary change of frame, to obtain our claim above, upon posing˜ ≡ β and h ≡ α − γ. Note that the initial state and the magnetization observable are unaffected by the last unitary transformation e iαZ , so the two Floquet operators are fully equivalent. Thus, we will consider the equivalent form of the Floquet operator with the modified kick in Eq. (33).
The key observation of Ref. [115] -exploited in Ref. [116] in the context of DTC response in NMR experiments -is that certain suitably chosen longitudinal field driving amplitudes proportional to the driving frequency may cause an averaging effect on the noncommuting transverse field perturbation, such that in the rotating frame co-moving with the strong drive the perturbation is heavily suppressed. This idea is implemented in our scheme similarly to what done in Ref. [116]: We rewrite our two-cycle Floquet operator as As in Refs. [115,116], here we consider the regime where J and˜ are small, so that we can apply a standard highfrequency expansion to the transformed Floquet operator U 2 . In our case, this expansion is equivalent to the plain BCH expansion. To lowest order, we obtain The lowest-order dynamics are solvable by mapping to free fermions, similarly to our treatment in Sec. III: rotating around the Z axis by a suitable angle we realize that the exponent in Eq. (38) describes a standard quantum Ising chain with coupling 2J and a transverse field generally proportional to˜ . Thus, excitations are domain-wall-like for˜ small enough. As we proved in Sec. III, the magnetization of this driven spin chain will decay exponentially with a rate γ ∼˜ 3 . Nevertheless, we may cancel the perturbation to lowest order by settingh = π(k + 1/2) with k integer: at these dynamical freezing points the two-cycle Floquet operator reads where the last equality holds to lowest order. As it is evident, in this case the decay rate of the magnetization becomes of higher order. Furthermore, unlike the general case, for these specific choices ofh the exact Floquet dynamics is solvable to all orders in J,˜ , as is evident from the first equality above: all four unitaries are mapped to quadratic fermions by the same Jordan-Wigner transformation. Computing the exact dispersion relation of the fermionic quasiparticles, we have verified that the excitations' bandwidth grows proportionally to˜ 2 for small /J -in other words, there is no accidental cancellation of O(˜ 2 ) corrections.
The argument above shows that our Floquet model with longitudinal field set to one of the dynamical freezing pointsh = π(k + 1/2) becomes an exactly solvable Floquet model, Eq. (39) above, similar to the one we considered in Sec. III, but with the strength of the perturbation˜ being replaced by˜ 2 due to a smart choice of the driving protocol. Applying our analysis of Sec. III, we conclude that the magnetization decays exponentially in this model, with a strongly suppressed rate γ ∼ (˜ 2 ) 3 =˜ 6 . We note that this lifetime is (much) longer than that found in the XXZ-type spin chain considered in Ref. [116]. Lastly, we remark that ifh is slightly detuned from a dynamical freezing point, the Floquet dynamics may still be analyzed using the theory developed in the previous Subsections, showing that the expected decay rate remains qualitatively unaltered (γ ∼˜ 6 ).

V. STABILIZATION OF DTC RESPONSE BY INTERACTIONS BEYOND NEAREST NEIGHBORS
A weaker version of domain-wall confinement also arises in the ordered phase of spin chains with an in-teraction range extended over multiple sites, even in the absence of explicit symmetry-breaking fields. The basic mechanism was identified in Ref. [60]: the separation of two domain walls involves an increase of the configurational energy, due to the increase in the number of frustrated bonds between pairs of spins beyond the nearest neighbors. As discussed in Refs. [60,61] and experimentally verified in Ref. [62], this gives rise to an effective attractive potential v(r) between two domain walls. The resulting physics is thus reminiscent of that generated by a longitudinal field. Here, however, the interaction tail tunes the shape and depth of the effective potential well, which, for the kind of interactions relevant to this work, grows sublinearly at large distances.
Below we show that this general domain-wall binding mechanism also arises in the Floquet context. Crucially, as it does not rely on explicitly breaking the symmetry, it does not suffer from time-averaging effects demonstrated in Sec. IV, opening the door to a true enhancement of the time-crystal order lifetime, as pictorially illustrated in Fig. 1(c). As a core result of this section, we establish that an increase of the interaction range R (i.e., J i,j = 0 only if |i − j| ≤ R) leads to an exponential enhancement of the order-parameter lifetime, γ ∼ a 2R+1 . We further discuss the practical conditions for observing this confinement-stabilized DTC behavior, and the implications for experiments with Rydberg-dressed interacting atomic chains, recently realized in the lab [72,74]. Based on this result, we will finally argue that long-range algebraically decaying interactions with a generic exponent α (not necessarily smaller than 2) stabilize the order parameter over timescales longer than any inverse power of the kick perturbation .
This section is organized as follows. In Sec. V A, we generalize the derivation of an effective domain-wall conserving Floquet Hamiltonian of Sec. IV B to a chain with arbitrary Ising couplings J i,j beyond nearest neighbors. Hence, in Sec. V B we study the two-body problem, which is richer than the corresponding one studied in Sec. IV C due to the coexistence of bound states and deconfined continuum, and determine the conditions for domain-wall binding and their implications for the evolution of the order parameter. Building on the intuition from the twobody problem, in Sec. V C we switch to a more scrupulous analysis and prove that, in the asymptotic regime of weak kick perturbation , the order-parameter decay rate gets heavily suppressed as γ ∼ a 2R+1 . This and related theory predictions are numerically verified in Sec. V D. Finally, in Sec. V E we discuss the limit of long-range interactions and argue that the decay becomes nonperturbatively small.

A. Effective domain-wall dynamics for weak kicks
We generalize Eq. (2) by considering a kicked Ising chain with arbitrary longer-range couplings, defined by 14 the Floquet operator JrZj Zj+r (40) with periodic boundary conditions. We have denoted by J ≡ (J 1 , J 2 , . . . , J R ) the array of coupling strengths at increasing distances, with π/2 > J 1 ≥ J 2 ≥ · · · ≥ J R > 0, and we implicitly assumed R < L/2. We consider the range R fixed and independent of the system size L (the long-range limit R ∝ L will be discussed in Sec. V E). For simplicity, we take the kick K as in Eq. (5), i.e., As shown in the previous Sec. IV, an explicit symmetrybreaking component h of the kick is not expected to qualitatively enhance the time-crystal lifetime. Thus, to lighten our arguments, we will discard it. We will later numerically verify that, indeed, those terms do not impact the lifetime of the DTC response. Note that, despite h = 0, the couplings beyond nearest neighbors break integrability. As in Eq. (6), transforming to the toggling frame gives so we focus on weak kicks K close to the identity. Along the lines of the first part of this paper, we want to build intuition on the evolution of the order parameter in terms of the dynamics of domain walls. The construction of the effective domain-wall conserving Floquet Hamiltonian of Sec. IV B can be straightforwardly generalized to the present case of the Floquet operator U = K V J with interactions beyond nearest neighbors. To this aim, we first work in the regime of "weak confinement" J 2 , . . . , J R J 1 . As the dominant scale J 1 in the problem couples to the number of domain walls D 1 [defined in Eq. (11)], we can set up a perturbation theory similar to that of Sec. IV B, where we had h J. The derivation closely parallels that of Sec. IV B, based on the general theory of Ref. [96]. In this case the perturbative parameter isˆ = max( , J 2 , . . . , J R ). The rigorous bounds of Ref. [96] ensure that the density of perturbatively dressed domain walls remains accurately conserved over the long prethermal timescale in Eq. (13), where the numerical constant c is here adjusted to account for the longer range R of the perturbation operator.
Similarly to the case discussed in Sec. IV B, for loworder explicit computations it is more practical to follow a different scheme from Ref. [96] and aim for a static effective Floquet Hamiltonian H F by combining the two exponentials of the product K V J using the replica resummation of the BCH expansion [104], order by order in the kick imperfection , and hence perform a conventional static Schrieffer-Wolff transformation on H F , order by order inˆ . The presence of arbitrary Ising couplings here requires a nontrivial generalization of the calculation in Ref. [104]. The structure of the resulting expansion is worked out in Appendix B.
The lowest-order result reported here is a simple generalization of Eqs. (16), (17): Due to the arbitrary range of the perturbation, the appearance of higher order terms in the effective Hamiltonian and in the Schrieffer-Wolff generator is more cumbersome than that discussed in Sec. IV B. Appendix B presents the general hierarchical structure of the expansion and explicitly reports the second-order result, for illustration. Terms of order p in H F contain at most p spin-flip operators S ± j ≡ 1 2 (X j ± iY j ) separated by a distance at most R, i.e., they contain strings of the form S µ1 j1 · · · S µq jq with µ 1 , . . . , µ q = ±, q ≤ p, and j 1 < · · · < j q with j n+1 − j n ≤ R for all n. Such a product of spin-flip operators (which we occasionally refer to as a quasilocal "cluster") is multiplied by complicated products of diagonal Z j operators -with coefficients depending on the couplings {J r } -located not farther away than R sites from the cluster. Finally, projectors P ↑ j , P ↓ j applied to spins adjacent to the position of the spin-flip operators ensure that the spin flips only move domain walls without creating or annihilating them. Likewise, the Schrieffer-Wolff generator at order p contains clusters of at most p spin-flip operators, with the same locality properties as discussed above.

B. Domain-wall binding
Similarly to Sec. IV C, we now analyze the dynamics of the order parameter in terms of the motion of conserved domain walls in the Schrieffer-Wolff-transformed picture. The initial state in the transformed picture |+ = e −iS ≤p * |+ can be viewed as a dilute gas of quasilocal clusters of p = 1, 2, . . . , p * flipped spins, the respective density being suppressed as 2p . To lowest order p = 1, these are isolated spin flips, i.e., adjacent pairs of domain walls. Thus, to this order of approximation, the evolution of the order parameter can be understood starting from the two-body problem.
Let us begin with qualitative considerations. Interactions beyond the nearest neighbors lead to the presence of a discrete set of bound states, coexisting with a continuum of unbound domain walls for larger energy. Furthermore, interactions between distant spins also favor the formation of more structured "molecular" bound states out of larger clusters of domain walls. The rich nonequilibrium dynamics of the system, including the anomalously slow decay of the order parameter, results from the coexistence of topological and nontopological excitations in the spectrum (i.e., unbound domain walls and bound pairs) which we now turn to quantitatively analyze.
The effective domain-wall conserving Floquet Hamiltonian (43) projected onto the two-particle sector gives the following first-quantized two-body problem, analogous to Eq. (19): where The two-body potential v J (r) grows as a function of the distance up to r = R, then flattens out. Thus, the potential well hosts a finite number of bound states, which grows to R−1 upon decreasing → 0. In field-theoretical language, this discrete set of energy levels forms the mass spectrum of nontopological particles. Above this, a continuum of scattering states appears, built out of two unbound domain walls; in field-theoretical language, the spectrum contains stable topologically charged particles, i.e., kinks and antikinks. The topological nature of these excitations stems from the fact that they can be locally created or destroyed in globally neutral pairs only, not individually. The bound (nontopological) and unbound (topological) excitations can be distinguished by being labelled by a real-space or momentum-space quantum number. To understand this, let us transform to the center-of-mass frame Ψ(j 1 , j 2 ) = e iK(j1+j2) ψ(j 2 − j 1 ) and set K = 0 due to translational invariance of the nonequilibrium initial state (cf. Sec. IV C). The reduced wavefunction ψ satisfies the Schrödinger equation in the domain r > 0, subject to the boundary condition ψ(0) ≡ 0. This equation defines the center-ofmass frame Hamiltonian H cm . For → 0, the eigenfunctions in the center-of-mass frame ψ (r) = δ ,r correspond to contiguous reversed domains of spins, with eigenvalues E = v J ( ). In this limit, the discrete label = 1, . . . , R−1 thus has the physical meaning of distance between the two domain walls. Due to the discreteness of the spectrum, this labelling can be adiabatically continued to finite , where eigenfunctions feature quantum fluctuations of the physical distance. On the other hand, the degenerate levels E R = E R+1 = · · · = v J (∞) split into a continuous band E k = v J (∞)−4 cos k, with eigenfunctions labelled by the relative momentum k of the two domain walls. The binding potential v J only affects these eigenfunctions via the relative scattering phase e iδ k between incoming and outgoing waves. We can estimate the stability range of the most excited [( = R − 1)-th] bound state by the condition that the continuum band E k does not overlap the discrete energy level E R−1 . This yields the range J R ; for larger , the bound state hybridizes with the unbound continuum. Energy levels and eigenfunctions are further renormalized by higher-order processes of orderˆ 2 ,ˆ 3 ,. . . to be added to H 2-body , not included in Eq. (45). These additional terms consist of diagonal terms and longer hops of domain walls by at most 2, 3, . . . lattice sites, respectively. For small enough , the resulting quantitative corrections do not alter the qualitative structure of the spectrum discussed here.  Fig. 4, both discrete "mesonic" bound states and unbound "domainwall-like" states appear. The latter are highlighted by blue lines in the pictorial sketch, and their energy levels are marked by thick blue bars, indicating that they form continuous bands of width ∝ . Panel (b) reports a selection of eigenfunctions of the two-body problem (47) in the center-of-mass frame, for algebraically decaying coupling J r = 1/r α truncated to the range r ≤ R = 10, α = 3, and = 0.1139. In this case, due to small couplings beyond the nearest neighbors and comparatively large hopping , the number of bound states (red curves) is only 3. As → 0, it grows to R −1 = 9. Higher excited eigenstates are unbound plane waves (blue curves).
The two-body problem already gives us important hints on the nonequilibrium evolution of the order parameter. As a matter of fact, confined domain walls only produce a weak oscillatory behavior of m(n) with frequencies related to the excitation energies (masses) of the bound states. Within the prethermal time window 0 ≤ t ≤ T preth in Eq. (13), the order parameter decay is only ascribed to the dynamical production of unbound domain walls. As discussed above, domain-walllike excitations exist as higher energy excitations. While in thermal equilibrium such domain-wall excitations are finitely populated, imperfect kicks K will only excite the domain-wall continuum very weakly, precisely with amplitude R , as generating unbound domain walls requires to flip a critical contiguous domain of at least R spins. This key insight suggests that the destruction of longrange order in the Floquet-prethermal Gibbs ensemble e −βH F /Z is a very slow process, leaving large room for nonequilibrium time-crystalline behavior.  (47) in the center-of-mass frame, for Jr = 1/r α for r ≤ R = 10, α = 3, and = 0.1139.

C. Exponential suppression of the order-parameter decay rate
Generalizing the reasoning of Sec. III B, the solution of the two-body problem suggests that the decay rate is severely suppressed: The density of critical domains in the initial state is of order ρ ∼ ( R ) 2 , and the spreading velocity of their constituents domain walls is v ∼ | |, leading to the estimate For R = 1, domain walls are the only stable excitations in the spectrum, and we recover the exact result γ ∼ 3 of Sec. III A. For R > 1, the appearance of bound states is expected to significantly slow down the order parameter decay. The argument based on the two-body problem is, however, too naive, as it completely neglects all multibody processes and interactions between confined domain walls. In fact, not only adjacent domain walls attract each other via the two-body potential v J (r), but also mesonic bound states themselves are subject to an effective attraction when their distance is less than R. [117] Such attractive forces can be thought of as residual interactions between bound states of elementary particles, physically analogous to nuclear forces that keep together protons and neutrons (i.e., bound states of quarks), or to molecular forces that keep together atoms (i.e., bound states of electrons and nuclei). The meson-meson attractive potential can be defined as (49) where E J ( 1 , r, 2 ) is the configurational energy of two reversed domains of length 1 and 2 , separated by r spins in between. One finds with n s ≡ min(s − r, 1 , 2 , r + 1 + 2 − s). Considering fixed 1 and 2 , this two-body potential w J, 1 , 2 (r) grows monotonically as a function of the distance r, from w(r = 1) < 0 to w(r ≥ R) ≡ 0. For instance, two isolated spin flips ( 1 = 2 = 1) experience an attractive potential proportional to the bare coupling, w J,1,1 (r) = −4J r . It is easy to imagine a process where a cluster of few flipped spins or small domains, initially distant from each other, agglomerate more tightly, making the gained reciprocal meson-meson interaction energy available to progressively enlarge a domain and finally release a domainwall from the attraction of the rest of the cluster at distance > R, such that it can freely travel away and melt the order parameter. In Fig. 9 we sketch such an example, where a finite (R-independent) cluster of initial spin flips possesses sufficient energy to release a traveling domain wall.
It is a-priori conceivable that such processes could trigger a comparatively fast decay of the order parameter at low perturbative order. Remarkably, however, it is possible to rigorously exclude such a scenario, and prove that the fastest decay process occurs at order R.
To establish this result, we need to modify the perturbative Schrieffer-Wolff transformation of Sec. IV B to explicitly account for the fact that all couplings J 1 , . . . , J R are large compared to when the asymptotic regime → 0 is considered. The unperturbed Floquet operator JrZj Zj+r defines highly degenerate sectors of the many-body Hilbert space, identified by the energy levels where E GS = −L r J r is the unperturbed ground-state energy of the fully polarized state |+ , and the nonnegative integer n r ∈ N has the meaning of total number of frustrated bonds at distance r. Under the assumption of strong incommensurability of the couplings {J r } in Eq. (B14) -necessary to derive a well defined static Floquet Hamiltonian, as shown in Appendix B 1-each degenerate sector is in one-to-one correspondence with the set {n r }. For → 0, transitions between such sectors are energetically suppressed and can be adiabatically eliminated. In other words, we can dress the effective Hamiltonian within each sector, order by order in , to account for all resonant processes occurring via virtual transitions between different sectors. The strong incommensurability condition on the couplings guarantees that each such transition is accompanied by a finite energy denominator. The construction can thus be formally carried out to all orders in [96]. Let us illustrate how this procedure works within our time-independent approach. Starting from the Floquet operator U = K V J , we combine the two exponentials into a Floquet Hamiltonian H F p by generalizing the replica calculation of Ref. [104] as detailed in Appendix B 1. Hence, we seek a modified Schrieffer-Wolff unitary transformation e i m Sm , iteratively for m = 1, . . . , p, which eliminates from H F p all terms of order m that violate the conservation of the any of the operators [H F R,p , D r ] = 0 for all r = 1, . . . , R.
where e iS ≤p ≡ e i S1 · · · e i p Sp . We remark that this scheme should be distinguished from the perturbation theory of Sec. V A, aimed at the conservation of the quantity D 1 only for the effective Floquet Hamiltonian H F 1,p ; on the contrary, the effective Floquet Hamiltonian H F R,p obtained here conserves all quantities D r , r = 1, . . . , R. To emphasize this distinction, here we use a calligraphic notation for the generator S of the perturbative scheme in powers of the kick , to avoid confusion with the previous generator S of the perturbative scheme in powers ofˆ .
As usual in the many-body context (cf. Sec. IV B), the proliferation of possible processes as well as the decrease of energy denominators at high perturbative orders is expected to lead to a divergent (asymptotic) character of the series, corresponding to a mutual hybridization of the many-body energy bands arising from the splitting of the highly-degenerate unperturbed levels. However, the perturbative series provides valuable information on the slowness of the dynamical delocalization in Hilbert space. The timescale T preth over which the approximate conserved quantities {D r ≡ e −iS ≤p * D r e iS ≤p * } R r=1 appreciably depart from their initial value can be estimated by finding the optimal truncation order p * of the series. [118] Within the time-dependent construction of Ref. [96], p * is found to depend on the magnitude of the perturbation as C −1/(2R+1+δ) where δ > 0 is arbitrary and C = C(δ) is a constant. This result leads to a stretched exponential lower bound on the prethermal timescale, which generalizes Eq. (13): Within the long Floquet-prethermal window 0 ≤ t ≤ T preth , the dynamics is guaranteed to take place within the Hilbert space sectors defined by the set of eigenvalues n r of the dressed operators D r , i.e., the (dressed) total numbers of frustrated bonds at distance r. Since the time-independent approach followed here (replica resummation + standard static Schrieffer-Wolff) basically produces the same set of conserved quantities as the timedependent scheme of Ref. [96], it is natural to assume that this alternative scheme yields similar nonperturbatively long heating timescales (see also the relative discussion in Refs. [96,104]).
Since H F R,p * conserves all operators {D r }, an initial configuration can only evolve within the corresponding resonant subspace with fixed number of frustrated bonds n r at distance r. The effective Hamiltonian H F R,p * is thus much more constrained than H F 1,p * in Eq. (43). The transformed-picture initial state e −iS ≤p * |+ is a lowdensity superposition of isolated quasilocal clusters of flipped spins. At pth order in , these clusters may comprise at most p flipped spins.
We are now in a position to clearly formulate our question on the order parameter evolution: Can an initial cluster of p < R flipped spins evolving via H F R,p * resonantly excite unbound domain walls? Under the single assumption that the array (J 1 , . . . , J R , 2π) is strongly incommensurate [as specified by Eq. (B14) in Appendix B 1], we prove that such a configuration with p < R is never resonant in energy with order-melting configurations, i.e., configurations possessing a contiguous domain of R reversed spins.
Referring to the illustration in Fig. 10, the claim is that the top and bottom configurations are necessarily separated by an energy mismatch. By the incommensurability condition, the existence of this energy mismatch is equivalent to the occurrence that n r = n L r + n R r for some 1 ≤ r ≤ R (note that the process in Fig. 9 is a particular case of Fig. 10).
The proof follows from the Lemma. A domain-wall-like configuration has n r ≥ r.
To show this, we consider a domain-wall-like configuration, with all spins pointing up (down) for j ≤ j L and down (up) for j > j R (where j R ≥ j L ). We focus on the r sublattices {i + jr|j ∈ Z}, labelled by i = 1, . . . , R. Each sublattice exhibits a domain-wall-like configuration, and hence, it has at least one frustrated nearest-neighbor bond. Each such bond maps to a bond at distance r in the original lattice. Thus, the original configuration has at least r frustrated bonds at distance r, QED. Using this lemma, it is easy to prove the main claim: considering again Fig. 10, the bottom configuration is energetically equivalent to two isolated domain-wall-like configurations, defined by the content of the regions denoted L and R, as the reversed bubble in between can be made infinitely large without frustrating any further bonds. By the lemma, n L r ≥ r and n R r ≥ r. In particular, n L R + n R R ≥ 2R. On the other hand, the top configuration in Fig. 10 has p flipped spins, each of which can frustrate at most two bonds at distance r, for all r. In particular, n R ≤ 2p. Hence, clearly, the resonance condition n R = n L R + n R R is impossible to satisfy when p < R, which concludes the proof, QED.
The result proved here implies that processes involving at least R spin flips are responsible for the order parameter decay. In particular, we already identified above the fastest such process, arising from terms . These terms give rise to isolated domains of R contiguous reversed spins, with density ∝ 2R , appearing in the transformed-picture initial state |+ = e iS ≤p * |+ at order R. The two domain walls separated by R sites are free to hop away from each other with amplitude , thus spreading the reversed domain and melting the system magnetization at a rate γ ∼ 2R+1 , cf. Eq. (48). Figure 2(a) reports a sketch of the timeline derived here; we reiterate that the heating timescale is only a lower bound, and might be non-tight in general.
We finally observe that tuning the longitudinal kick component to a dynamical freezing point as outlined in Sec. IV E, the decay rate of the order parameter may be further strongly suppressed, as the magnitude of the perturbation gets effectively replaced by 2 . Our theory shows that this results in a decay rate γ ∼ 4R+2 . As interactions beyond the nearest neighbors are quite generic in experiments, our results indicate that, upon tuning the driving parameters to a dynamical freezing point, DTC response might appear perpetually stabilized for all practical purposes.

D. Numerical simulations for R > 1
The above theoretical analysis is supported by extensive numerical simulations that we performed for R = 2 (namely, for the Ising chain with next-to-nearestneighbor interactions). In our numerics, the couplings have been fixed as J 1 = 1/ζ(α) and J 2 = (1/2) α /ζ(α), where ζ(x) = ∞ r=1 1/r x denotes the Riemann zeta function. We fix either α = 2.25 or 3 (thus larger than 2), to be consistent with the next section concerning long-range systems. This choice is largely arbitrary at this level; we anticipate that it allows a direct comparison with the data presented in the next section. Note that the incommensurability of the couplings {J r } would be guaranteed for irrational α's, but we take rational values to further test of the robustness of our analytical predictions.
We found that the Floquet evolution of the absolute value of the order parameter is always compatible with the theoretically expected exponential decay for all kick imperfections = 0. To validate our theory prediction on the suppression of the decay rate γ compared to the deconfined case of nearest-neighbor interactions, we performed thermodynamic-limit iTEBD simulations for a sequence of small values of . Results are reported in panels (a) and (b) of Fig. 11. Note that the curvature of the exponential decay is hardly visible on the accessible timescale for such small values of ; however, γ can be accurately extracted. The resulting scaling of γ vs is reported in panel (c), which clearly shows how the numerical data points follow the theory prediction (dashed black lines). The nearest-neighbor analytical result, given by Eq. (7) and valid for R = 1, is also reported for comparison (dotted black line). We finally report in panel tudinal kick component h = 0.4. As expected from our analysis of Sec. IV, the decay rate is essentially unaltered.
In Fig. 12 we compare the thermodynamic data with the finite-size ED results. (Note that for small the entanglement entropy growth is slow enough to push iTEBD simulations to unusually large numbers of driving periods [119].) Also in this case, the absolute value of the order parameter gets stuck to a nonzero value for L < ∞. However, both iTEBD data and finite-size results show the same frequency of oscillations.
The frequencies of such oscillations can be extracted from the Fourier power spectrum of the time series for |m(n)| (Fig. 13). We observe that data in the thermodynamic limit and at finite size are consistent, and manifest the same main peaks symmetrical around ω = π [panel (a)]. (Note that the finite-size spectrum shows some spurious frequencies due to time-recurrence effects.) Further confirmation of the validity of our perturbative analysis is provided by the scaling of the position of the main peak, which is approaching the "classical" value of the single spin-flip excitation 4J 1 + 4J 2 for → 0 [panel (b)].

E. Long-range limit
The results of the previous section demonstrate that the order parameter decay rate is suppressed as γ ∼ a R 2R+1 , as → 0. In this last section, we discuss the long-range limit R → ∞. In particular, we focus on the experimentally relevant case of algebraically decaying interactions J r = J/r α . This model is relevant to the dynamics of effective qubits in quantum simulators based on trapped ions (tunable 0 < α < 3) [120,121] and Rydberg atoms (α = 6) [122,123]. In these setups, nontrivial quantum dynamics are generated by additional magnetic fields acting on the spins, whose spatiotemporal variations can be efficiently controlled in the experiment. We thus hereafter interpret J = {J 1 , . . . , J R } as a finite-range truncation of the parent sequence of cou- [Note that for a generic (irrational) value of α, these truncated sequences are expected to satisfy the strong incommensurability condition in Eq. (B14).] Within this perspective, it is interesting to shift our viewpoint to the functional dependency of the rate γ on R.
Extracting the scaling of the prefactor a R would involve keeping track of the magnitude of the subset of processes of order R which trigger the order parameter decay in our combined replica + Schrieffer-Wolff transformations. This is in principle straightforward but practically unfeasible, due to the rapid growth of the complexity of high-order perturbative computations. However, as a crude conservative estimate, we can bound a R from above by the total magnitude of all terms of order R . This type of bounds are worked out in the related analysis of Ref. [96], as well as in many previous works on rigorous prethermalization theory [7,26,27,100,124], to obtain estimates of the thermalization timescales, like Eq. (55). The ubiquitous scenario resulting from these works is that the total magnitude (measured by a relevant operator norm) of all terms perturbatively generated at order p first decreases exponentially with p, before plateauing at p = p * and finally diverging rapidly. In the case of interest here, Ref. [96] finds p * ∼ −1/(2R+1+δ) . Since any finite range R is largely superseded by p * for small enough perturbation , the exponential suppression 2R+1 dominates over the prefactor a R , and the decrease of γ upon increasing R is effectively exponential.
However, taking the limit R → ∞ is subtle, as it does not commute with the asymptotic perturbative limit → 0: Setting heuristically R = ∞, the heating bound in Eq. (55) trivializes. Taken literally, this occurrence suggests that a fast violation of the effective conservation laws of {D r } has to be expected for the long-range interacting system, which could in principle lead to a fast order parameter meltdown. In particular, this would preclude any meaningful extrapolation of the results of the previous section to long-range interactions.
Numerical simulations, however, suggest the opposite behavior: we find that increasing the truncation radius R of long-range interactions J r = J/r α to the maximum R = L/2, for arbitrary α, leads to a dramatic increase of the order parameter lifetime, as clearly visible from the data shown in Fig. 14  In these simulations, we have chosen Kac rescaled interactions, i.e., J = 1/ζ(α) (cf. Fig. 11), so that results for different α can be fairly compared. The reported data point at a robust stabilization of the DTC signal, well beyond our analytical theory of Sec. V C: in fact, one can see that the kick strength taken there, = 0.114 and 0.171, correspond to quite big rotations of the spins at each kick, by angles ≈ 13 • and ≈ 20 • , respectively. These perturbations are actually much larger than the considered couplings beyond nearest neighbors. In spite of this, the order parameter decay is extremely suppressed upon increasing R, and hardly visible in the long-range limit; moreover, this occurrence is not a finite-size effect.
To resolve this apparent contradiction, we observe that the bound (55) on prethermalization is unnecessarily pretentious for our purposes: it expresses the expected timescale of quasiconservation of a large number of operators {D r } R r=1 , uniformly in the many-body spectrum. As we take R = L/2, energy levels become infinitely dense in the thermodynamic limit away from the band edges, thanks to the strong incommensurability condition which prevents them from being degenerate. In fact, the preservation of an extensive number of commuting operators {D r } L/2 r=1 would make the system effectively many-body localized, contrary to conventional delocalization scenarios for translationally invariant models [37,42]. The slow dynamics of highly excited states resulting from this long-range limit is thus nonstandard, and the actual heating timescales (or thermalization timescales for time-independent systems) are presently unclear, even for static (undriven) systems [60,61,75,76].
On the other hand, the evolution of the order parameter relevant to this work takes place in a particular corner of the many-body Hilbert space, corresponding to the low-energy sector of an approximate Floquet Hamiltonian. While the perturbative series might be severely divergent at low orders in the long-range limit when measured by uniform operator norms, the same bounds are far too loose when the construction is restricted to a low- energy sector with dilute excitations, relevant for the purpose of this work. [125] In other words, while low-order perturbative transitions from highly excited configurations with extensively many frustrated bonds are very likely to hit resonances, this becomes extremely unlikely when the initial state is the polarized state |+ considered in this work, since quasilocal clusters of flipped spins generated by the weak kicks are far away from each other and hence can hardly cooperate to produce resonant transitions. This argument suggests that the lower bound on the timescale T preth in Eq. (55) is far too conservative for the dynamics originating from the state |+ , for arbitrary R. A tighter bound would be needed to correctly account for the decrease of the density of states at low energy. In the same regime, a long-lasting suppression of heating has to be expected in the limit of long-range interactions R → ∞, as well. We reiterate that, even for finite-range interactions, our numerical simulations indicate a stability well beyond the analytical theory presented abovecf. Figs. 1, 14. Along the lines of the discussion above, our intuition is that while plenty of resonances are bound to occur when the range is long, in the dilute quench dy-namical setup of our work many of these resonances are inconsequential, as the only "dangerous" processes for the order parameter meltdown are those that cause an inner rearrangement of a cluster of domain walls involving the relocation of the leftmost (or rightmost) of them to distance R from the next. It is likely that such specific resonances form a much smaller set that is further dynamically obstructed compared to the naive expectation from brute-force norm bounds.
Setting up a generalized Schrieffer-Wolff perturbative scheme aimed at estimating the timescales involved in the intricate slow dynamics of the long-range interacting chain appears as a formidable problem, which we leave to future investigations. Here, building on the insight of Sec. V C and above, we formulate the conjecture that the timescale of the fastest process leading to the order parameter decay can be estimated in terms of the initial density of bubbles of reversed spins, whose walls are free to spread away from each other. Remarkably, such plausible scenario leads to a decay rate γ beyond perturbation theory, i.e., smaller than any power of .
As suggested in the discussion above, even though the elimination of domain-wall nonconserving processes is formally valid throughout the entire many-body spectrum only for α 1, in the low-energy sector one can naively perform several perturbative steps to preserve D 1 (and a few other operators D r , as well) for arbitrary α. Thus, we reconsider the two-body problem of Eq. (45), and set J r = J/r α and R = ∞. We obtain where v α (r) = 4 The potential grows from v α (r = 1) = 4Jζ(α) to v α (r = ∞) = 4Jζ(α − 1) (for α > 2) or ∞ (for α ≤ 2). The asymptotic behavior at large r is for α > 2.
(59) Let us now discuss the result in Eq. (59). For α ≤ 2, the binding potential is confining at large distances. Hence, the hierarchy of domain-wall bound states exhausts the excitation spectrum [60]. The infinite energetic cost of isolated domain walls for α ≤ 2 underlies the persistence of long-range order at finite temperature [126,127], circumventing the standard Landau argument against the existence of thermal phase transitions in one dimension. Indeed, at a prethermal level, the effective Floquet Hamiltonian resulting from the resummation of the BCH expansion supports long-range order at finite temperature for small but finite . On the other hand, the two-body potential is bounded at large distances for α > 2. In this case, freely travelling domainwall states appear, similarly to finite-range systems. Due to their nonlocal nature, isolated domain-wall-like excitations cannot be locally created or destroyed; thus, they give rise to topologically protected quasiparticles, which form a continuum in the excitation spectrum above a discrete sequence of nontopological bound states. Unlike the case R < ∞, however, the potential only flattens out asymptotically for r → ∞. Consequently, the number N α of such bound states critically depends on the hopping amplitude . The finite statistical weight of domain-wall-like quasiparticles in thermal equilibrium is what prevents long-range order at finite temperature for α > 2. Figure 15 reports an illustration of the two-body spectrum, obtained within a semiclassical approximation Here we take = 1, so that h ≡ 2π. The continuous function vα(Q) has been taken as in Eq. (59). For this choice of parameters, the center-of-mass potential well hosts four bound states, marked by red trajectories, bounded in Q. All higher excited eigenstates (blue trajectories) are unbound plane waves, and form a continuum labelled by the asymptotic momentum P for Q → ∞.
(which becomes quantitatively accurate in the continuum limit, i.e., for highly excited states). Here we set α = 3 and = 0.1139. Quantized trajectories undergo a transition between spatially localized (red) and delocalized (blue), representing a discrete sequence of nontopological confined bound states below a continuum of topological unbound domain walls. For this choice of parameters, N α ( ) = 4; however, this number grows as → 0. A clear signature of these bound states is further given by the presence of pronounced peaks in the power spectrum of the absolute value of the order parameter time series. In Figure 16 we show the Fourier transform of finite-size data for α = 2.25. In the long-range limit (R = ∞) the two main peaks correspond to the lowerenergy bound states, namely the single and the double spin-flip excitation. In the classical limit ( → 0) their energies are given respectively by ω 1 = 4 L/2 r=1 r −α /ζ(α) and ω 2 = ω 1 + 4 L/2 r=2 r −α /ζ(α). On the contrary, the main signatures in the order parameter evolution of the presence of unbound domain walls in the spectrum is the (very slow) overall decay of the signal in Fig. 14.
Generalizing our argument of Sec. V B, the number of bound states N α can be estimated starting from the observation that isolated domain walls can freely hop to neighboring sites with amplitude , which gives the dispersion law E k = 1 2 v α (∞) − 2 cos k. The unperturbed bound-state wavefunctions for → 0, ψ (r) = δ ,r , ∈ N, are precluded from hybridizing with the domain-wall continuum when their unperturbed energy E = v α ( ) is below the "ionization threshold" v α (∞) − 4 . The equation thus identifies the highest stable bound state ≡ N α . Using the asymptotic expansion in Eq. (59), we obtain where c α = 4(α − 2)(α − 1). This result expresses how the number of bound states diverges as → 0 for all α > 2. Accordingly, the physical size of a critical reversed bubble triggering the order parameter meltdown grows unbounded in the asymptotic regime of weak perturbation. In other words, our conjecture that the formation of a critical bubble is the fastest process leading to the decay of the order parameter, suggests that this decay has a nonperturbative origin in the long-range interacting spin chain, where we have defined A ≡ 2(c α /J) −1/(α−2) . We note that as α approaches 2 from above, the lifetime γ −1 from Eq. (62) diverges. The lack of an appropriate heating bound in this regime (as discussed above) prevents us from estimating the location of a presumable crossover region α ≈ α * ≥ 2 between a vacuumdecay driven (α α * ) and a heating-driven (α α * ) order-parameter decay. In any case, we reiterate that the order-parameter lifetime is expected to be nonperturbatively long for all α's.
An explicit comparison between the various timescales can be drawn in the high-frequency driving limit of our Floquet model, i.e., taking (J/r α )Zj Zj+r (63) with τ small. In this case, rigorous bounds on the heating timescale are available for all values of α. We refer to Fig. 2(b) for a comprehensive illustration. For 1 < α ≤ ∞, an exponential lower bound T preth ≥ exp(C/τ ) applies uniformly in α [64] (assuming Kac rescaling of J as above). Within this long time window, heating is prevented by the quasiconservation of an effective Hamiltonian emerging in the toggling frame, of the form Z j Z j+r r α + X j +τ (. . . )+τ 2 (. . . )+. . .
(64) For α > 2, this Floquet Hamiltonian governs the evolution of local observables [65]. The arguments of this Section apply to the dynamics generated by H F , leading to the long lifetime 1/γ of the confinement-stabilized DTC order parameter in Eq. (62). Knowledge of the heating timescale allows us to keep track of a presumable crossover value α * = α * ( , τ ) > 2 between this and the vacuum-decay timescale. Upon decreasing → 0, the value of α * is pushed to ∞; conversely, as τ → 0, it is pushed towards 2. (We reiterate that the heating timescale is only a lower bound, and might be non-tight in general.) For 1 < α ≤ 2, the Floquet Hamiltonian supports long-range order at finite temperature. Thus, the mechanism of Floquet prethermalization suffices to stabilize DTC order. Its lifetime is only limited by heating processes [7,65]. Finally, for 0 ≤ α ≤ 1, dynamics starting from a fully polarized initial state (such as our |+ ) is governed by an emergent semiclassical description over a timescale that diverges with system size [66,67], and the underlying asymptotic "decoupling" between the few collective degrees of freedom coupled to the drive and the extensive set of microscopic degrees of freedom, prevents the system from absorbing energy and heating up [128,129]. Because of this nontrivial mechanism, the emergent DTC behavior has mean-field character in this regime [9,130].
To summarize, while the Floquet-prethermalized state features a nonvanishing density of travelling domain walls precluding large-scale spatiotemporal order, the dynamical production of these excitations by kick imperfections goes through extremely slow vacuum-decay processes, delaying the Floquet prethermalization itself, ultimately paving the way for a long window of genuinely nonequilibrium time-crystalline behavior compatible with numerical simulations. Our argument that the fastest ordermelting process is the fractionalization of a critical-size reversed bubble into a deconfined kink-antikink pair leads to the nonperturbatively long lifetime in Eq. (62).

VI. CONCLUSIONS
In this paper, we have established a framework to understand and compute the order-parameter evolution in periodically driven quantum spin chains, hinging upon the effective dynamics of emergent domain-wall excitations. Within this framework, we have analyzed the impact of domain-wall confinement on the order-parameter decay, and established that a slight increase of the interaction range can result in a dramatic extension of the lifetime of DTC response, despite the absence of a longrange ordered Floquet-prethermal states. The results of this paper delimit and characterize the theory of time crystals for disorder-free, finite-range interacting quantum spin chains driven at arbitrary frequencies, extending the current state of the art.
Observing confinement-stabilized DTC response does not require disorder, high-frequency drives, or fat-tailed interactions, which may not be easily accessible in many experimental setups. A naturally suited platform for implementation of the physics discussed here is given by Rydberg-dressed spin chains [72][73][74]. The tunable-range Ising interactions realized in these systems, described in our notations as have an almost flat core for r r c , and cross over to a quick decay in the intermediate range r ≈ r c . The value of the effective range r c can be efficiently tuned in the experiment, making this setup ideal to observe the confinement-stabilized DTC response identified in this work. Indeed, necessary ingredients such as preparation of fully polarized states, application of global pulses, long coherence times and monitoring of the collective magnetization, have already been demonstrated in these experiments.
An important remark, however, is that the set of initial states which give rise to such a response is more limited than for the prethermal DTC. The flexibility in perturbing the initial state by arbitrary local operations is set by the same parameter scale that controls the duration of the signal. This is, in a sense, reminiscent of DTC behavior associated with quantum many-body scars [30,131]. More generally, this work provides further evidence that increasing the range of interactions may generate nonthermal behavior in certain regimes [60,61,66,67,[132][133][134][135], and may help realizing genuinely nonequilibrium phases [65,128].
The results reported here clarify some confusion on the role of the system size in clean short-range interacting "time crystals", consistent with the numerical analysis of Ref. [79]. Furthermore, a crucial byproduct result of this work is the clarification of the nature of the apparent anomalous persistence of the order parameter observed in several numerical investigations of long-range quantum Ising chains with α 2 after a global quench of the transverse field from a ferromagnetic ground state [60,75,76]. Due to pronounced finite-size effects and severe slowdown of the dynamics in this regime, purely numerical calculations face significant challenges. The theory developed in this paper for the more general Floquet setting (Sec. V) applies equally well to these quench dynamics (of course, incommensurability of the couplings with 2π need not be assumed in this case). This work provided solid analytical evidence that a strong enhancement of the order parameter lifetime is to be generally expected in the parameter regime 2 < α ∞, and predicted its functional form as a function of the quench magnitude. This has the important consequence, seemingly not recognized before, that a long-lived nonequilibrium order can be sustained by systems that cannot exhibit long-range order in equilibrium, via the suppression of the dynamical creation of deconfined topological excitations (domain walls), ultimately responsible for melting the metastable nonequilibrium order. This occurrence is related to familiar macroscopic tunnelling phenomena in high-energy physics, such as the Schwinger mechanism [70] and Coleman's false vacuum decay [71], where atypical excited states may only decay (and hence thermalize) through slow nonperturbative processes.
On a technical level, a few points remain open and are left to future work, including a more rigorous derivation of Eq. (62), and a better understanding of the role of resonances in the unperturbed spectrum.
with D, j Z j Z j+1 = 0 and V purely off-diagonal between sectors with different number of domain walls. Explicitly we have Now we fix S 1 in such a way to exactly cancel V in the transformed Floquet operator, i.e., This condition for S 1 is solved by which corresponds to Eq. (17) in the main text. With such a choice, one recovers the Floquet operator U 1 = e −iH ≤1 in Eq. (16).