Long-Range Prethermal Phases of Nonequilibrium Matter

We prove the existence of non-equilibrium phases of matter in the prethermal regime of periodically-driven, long-range interacting systems, with power-law exponent $\alpha>d$, where $d$ is the dimensionality of the system. In this context, we predict the existence of a disorder-free, prethermal discrete time crystal in one dimension -- a phase strictly forbidden in the absence of long-range interactions. Finally, using a combination of analytic and numerical methods, we highlight key experimentally observable differences between such a prethermal time crystal and its many-body localized counterpart.


I. INTRODUCTION
Periodic driving represents one of the most versatile tools for manipulating quantum systems. Classic examples of this abound in magnetic resonance spectroscopy, where it has been used for more than half a century to help narrow spectral line shapes [1][2][3]. More recently, in the context of cold atomic gases, periodic driving has also helped to enable the realization of novel types of many-body interactions [4][5][6][7].
Despite this ubiquity, one place where periodically driven (Floquet) systems have traditionally remained absent is in the study of phases of matter [8][9][10]. Indeed, the usual, statistical mechanical framework for characterizing phases has largely been restricted to the exploration of systems at or near equilibrium. Floquet systems do not fit this category. Rather, they can continuously absorb energy from the driving field, ultimately approaching an infinite-temperature thermal state at late times [11][12][13][14][15][16][17][18][19][20][21][22][23]. As a result, in the thermodynamic limit, the naive conventional wisdom is that all many-body, Floquet systems must behave trivially from the perspective of phases of matter. However, seminal recent works have called this assumption into question.
For example, the presence of strong disorder in one dimension (and possibly higher dimensions) can prevent thermalization by inducing a many-body localized (MBL)  1. (a) Schematic phase diagram for a one-dimensional prethermal time crystal as a function of interaction power law and energy density. The 1D PDTC can only exist for long-range interactions (i.e., J ij ∝ ji − jj −α ) with power law 1 < α < 2 and an energy density that lies in the symmetry-broken phase of the prethermal Hamiltonian D Ã . (b) PDTC Floquet dynamics depicting the magnetization MðtÞ for a system size L ¼ 28. The robust period doubling behavior, which survives for exponentially long times in the frequency of the drive ω, signals prethermal time crystalline order. (c) Table summarizing our analytical results. The star indicates that for this case prethermal phases exist provided that we assume that local observables relax to the Gibbs state of D Ã , which we expect since this is the state that maximizes the entropy subject to the constraint of conservation of energy.
phase [24,25]. When an MBL phase occurs in a Floquet system [13,17,26,27] it can prevent energy absorption from the drive and lead to novel, intrinsically out-of-equilibrium phases of matter [28][29][30][31][32][33][34][35]. However, the dual constraints of strong disorder and low dimensionality significantly limit the scope of both the experiments and models that one can consider, naturally raising the question: can interesting Floquet phase structure survive in periodically driven systems without disorder? An affirmative answer has recently emerged [36] in the context of Floquet prethermalization [37][38][39][40][41][42]. For sufficiently large driving frequencies, a many-body Floquet system can enter a so-called "prethermal regime," where its dynamics are well captured by an effective static Hamiltonian. This static Hamiltonian description necessitates the existence of a conserved energy, which prevents the driven system from heating to an infinite temperature state. Crucially, the lifetime of this prethermal regime has been proven to be exponentially long in the frequency of the drive, providing a parametrically robust mechanism to delay the onset of Floquet heating.
Although such results further cement the power of periodic driving as a technique for Hamiltonian engineering [43][44][45][46], we hasten to emphasize that these results are necessary but not sufficient for proving the existence of intrinsically nonequilibrium, prethermal Floquet phases of matter. Let us unpack this last statement. Our focus in this paper is on phases of matter that cannot exist in equilibrium. This means that the Floquet nature of the system is not simply being used as an engineering trick to stitch two disparate Hamiltonians together, but rather, as a prerequisite ingredient for the existence of a phase with no direct analog in thermal equilibrium. This latter point is most easily summarized as follows: the phase must, at its core, be protected by the discrete time-translation symmetry of the drive [32,33,36].
Thus, in order to prove the existence of prethermal Floquet phases, one must first demonstrate that the prethermal regime can actually preserve the symmetry structure of the driven system. With this in mind, recent progress has precisely demonstrated the existence of emergent symmetries during the prethermal window [36]. The existence of these symmetries can be viewed as a direct manifestation of the discrete time-translation symmetry of the drive. This theoretical framework provides the perfect landscape for realizing prethermal nonequilibrium phases of matter, including prethermal versions of discrete time crystals [28,34], Floquet symmetry protected topological phases [29,32,33,47], and possibly many others [48][49][50][51]. However, this framework leaves open one fundamental challenge, in that it cannot be applied to long-range interacting systems.
More specifically, one cannot ensure that the resulting effective prethermal Hamiltonian possesses any meaningful sense of locality. Without this notion of locality, the evolution of local operators may not be well approximated by the prethermal Hamiltonian. As a result, the usual assumption that the system will evolve to the prethermal Gibbs state and exhibit the phase structure of local and power-law interacting Hamiltonians may not hold. The overarching goal of our work is to tackle this concern, proving the existence of prethermal Floquet phases in manybody systems that exhibit long-range, power-law interactions (i.e., Coulomb, dipolar, van der Waals, etc.) [52][53][54][55][56].
This goal is motivated from two complementary fronts. On the experimental front, many of the platforms most naturally suited for implementing Floquet dynamics exhibit long-range interactions, including dipolar spins in the solid-state, trapped ions, ultracold polar molecules, and Rydberg atom arrays [56][57][58][59][60][61]. Understanding the prethermal properties of this broad class of systems could unlock a myriad of new experimental techniques for Floquet quantum simulation. On the theoretical front, even in equilibrium, it is well known that long-range interactions can lead to symmetry breaking in qualitatively different regimes than that allowed by short-range interactions. This suggests the possibility of finding prethermal Floquet phases that can only be realized in long-range interacting systems.
Our main results are threefold. First, we prove the existence of prethermal Floquet phases of matter ( Fig. 1) in long-range interacting systems, so long as the interactions decay as a power law with exponent α > d, where d is the dimension of the system. Second, we predict the existence of a novel, disorder-free, prethermal discrete time crystal (PDTC) in one dimension. This phase is strictly forbidden in any of the three contexts that we discussed earlier: equilibrium, Floquet MBL, and short-range interacting prethermal Floquet. Indeed, the 1D PDTC can only be realized in a long-range interacting, prethermal Floquet system. Finally, leveraging large-scale Krylov subspace methods, we perform extensive numerics characterizing the emergence of a 1D PDTC in a long-range interacting spin chain. In this context, we highlight one of the key (experimentally observable) differences between the prethermal time crystal and the MBL time crystal, namely, the presence of a phase transition as a function of energy density ( Fig. 1 and Table I).
Our paper is organized as follows. In Sec. II, we lay the framework for understanding Floquet prethermalization TABLE I. Differences between MBL and prethermal discrete time crystalline order in one-dimensional systems. The star next to "Short-range" indicates that the range of the interaction must only be sufficiently short so that MBL is preserved.

MBL TC Prethermal TC
Lifetime τ → ∞ τ ∼ e ω=J local Initial state Any Below T c Requires disorder Yes No Interaction range Short-range* Long-range 1 < α ≤ 2 both with and without an emergent symmetry (although only the former admits nonequilibrium phases of matter). Moreover, we review and contextualize a number of prior results with a particular emphasis on their implications for understanding the dynamics within the prethermal regime. This allows us to formalize the two essential properties for proving the existence of long-range interacting, prethermal phases. Building upon these discussions, in Sec. III, we begin by introducing new machinery to carefully keep track of the spatial structure of the long-range interactions. Leveraging these new tools, we ultimately prove three theorems, which in combination demonstrate the existence of long-lived, nonequilibrium prethermal phases of matter in long-range interacting systems with power laws α > d.
Within this context, we also introduce a novel phase of matter, the 1D prethermal discrete time crystal. In Sec. IV, we perform an exhaustive numerical investigation of a onedimensional Floquet spin chain and demonstrate that it exhibits a PDTC phase, only when the system harbors sufficiently long-range interactions. Using a combination of Krylov subspace methods and quantum Monte Carlo calculations, we identify one of the unique signatures of a PDTC (as compared to an MBL discrete time crystal), namely, that it displays a phase transition as a function of the energy density of the initial state. Finally, we provide a short summary of some of the implications and interpretations in Sec. V.

II. PRETHERMALIZATION
In an interacting, many-body quantum system, one generally expects dynamics to push the local state of the system toward equilibrium via a process known as thermalization [62][63][64][65]. However, in certain cases, the timescale τ Ã at which thermalization occurs can be significantly larger than the timescale associated with the intrinsic local interactions of the Hamiltonian 1=J local [66]. In such cases, before full thermalization actually occurs (i.e., for times t < τ Ã ), the system can first approach a different equilibrium state determined by an effective Hamiltonian-this process is called prethermalization; the time interval associated with it is known as the prethermal regime, while the effective Hamiltonian is referred to as the prethermal Hamiltonian.
Systems exhibiting prethermalization generally have two distinct energy scales. In static systems, this typically requires the underlying Hamiltonian to exhibit two very different couplings which lead to both "fast" and "slow" degrees of freedom (d.o.f.). Prethermalization can then be understood as the equilibration of the fast d.o.f. with respect to a slowly varying background arising from the dynamics of the slow d.o.f. In this case, τ Ã is expected to depend algebraically on the ratio of the energy scales [67].
Exponentially long Floquet heating time.-Unlike static systems, Floquet systems always exhibit two distinct energy scales: the local energy scale J local and the frequency of the drive ω. To this end, a Floquet system can almost naturally be expected to exhibit a long-lived intermediate prethermal regime when these two energy scales are sufficiently different; our focus is, of course, on the case in which ω ≫ J local . In that case (typically referred to as Floquet prethermalization), τ Ã scales exponentially with the ratio of these two energy scales, ω=J local , rather than algebraically [37][38][39][40][41].
The physical intuition for this exponential scaling is simple. Given a local energy scale J local , the many-body system requires ω=J local rearrangements in order to absorb a single quantum of energy from the drive. When interactions are local, the system cannot efficiently make a large number of correlated local rearrangements. Thus, the associated rate of energy absorption (i.e., Floquet heating) is exponentially small in ω=J local , leading to a heating timescale, τ Ã ∼ e ω=J local . This physical picture also helps to explain why long-range interacting Floquet systems with power laws α < d cannot exhibit a prethermal regime. In such systems, the energy scale associated with a single local rearrangement diverges as a function of the system size (i.e., the system exhibits a superextensive many-body spectrum), implying that a single local rearrangement can, in principle, absorb an energy quantum from the drive regardless of the magnitude of the driving frequency.
Approximation of local Floquet dynamics.-While we focused above on the existence of an exponentially long Floquet prethermal regime, as we alluded to earlier (while emphasizing the importance of locality), this is not the only constraint that one needs to worry about. Rather, just as important is whether one can prove that there actually exists a local prethermal Hamiltonian D Ã that approximately generates the dynamics of the Floquet system during the prethermal regime. A bit more precisely, the unitary time evolution operator U f that generates the exact Floquet dynamics during a single driving period T should be approximated by: And, more importantly, one hopes that this approximation correctly captures the dynamics of local observables until the Floquet heating timescale. A priori, this need not be the case, and, in fact, the exact Floquet dynamics might not have any effective Hamiltonian description. Indeed, the difference between proving the existence of a conserved energy (i.e., measured with respect to the prethermal Hamiltonian) versus proving that the prethermal Hamiltonian correctly generates the local dynamics is stark. For example, although the Floquet heating time τ Ã has been proven to be exponentially long in generic systems with extensive energy scales (including long-range interacting systems [36,[38][39][40][41] and even classical systems [71]), proving that the associated prethermal Hamiltonian describes the dynamics of local observables has only been achieved for a significantly smaller class of systems [36,38,40,72]. In fact, in certain systems it has been shown that the prethermal Hamiltonian does not generate the actual Floquet dynamics [71].
Generalizing to the case of an emergent symmetry.-Up to now, we have focused on how an effective static description of the Floquet system (governed by the prethermal effective Hamiltonian) can emerge during the prethermal regime, both in the context of a conserved energy as well as in the context of generating local dynamics. While powerful in and of itself, this description limits Floquet systems to mimicry of equilibriumlike physics within the prethermal regime. This is because, at the moment, our effective static description has forgotten about the structure of the original time periodic drive. Luckily, this need not be the case.
Before formalizing this last statement, let us illustrate it with a simple example. Consider an S ¼ 1=2 spin undergoing a π=2 rotation every period T. In the absence of any perturbing field, the spin will return to its original orientation every four periods. Crucially, it turns out that even in the presence of small interactions (with respect to the driving frequency ω ¼ 2π=T), this picture remains true for an extremely long timescale. One can gain some intuition for this by noting that all of the interactions which fail to commute with the π=2 rotation get "echoed out" (i.e., they average to zero in the toggling frame that rotates by π=2 each Floquet period), which means that at leading order in the inverse frequency, they do not contribute to the dynamics. We emphasize, however, that the general results we eventually consider will hold not just at leading order but also at higher orders.
Armed with this simple example, let us now formalize how extra symmetry structure can emerge in the prethermal regime of Floquet systems. In particular, if U f contains a large rotation X that returns to itself after N periods, X N ¼ 1 (in our example with the π=2 rotation, N ¼ 4) and generic interactions (whose strength is much smaller than the driving frequency), then U f can be exponentially well approximated by a much simpler evolution [36]: where D Ã is the effective prethermal Hamiltonian that commutes with the rotation X and U is a time-independent unitary change of frame, which is close to the identity. Note that we will often choose to work directly in the rotated frame given by U, so that the evolution is (approximately) given byŨ app f rather than U app f . The above discussion encodes a few important consequences. First, since D Ã commutes with X, it remains an exactly conserved quantity under this approximate evolution. Taking into account the exponentially small error terms (which track the differences between this approximate evolution and the exact Floquet evolution) leads to D Ã being exponentially well conserved. Second, while X was not a symmetry of the original evolution, it has become a Z N symmetry of the approximate time evolutionŨ app f ; this emergent symmetry is protected by the underlying discrete time-translation symmetry of the Floquet evolution operator. As we discuss later, one can leverage this emergent symmetry to realize novel Floquet phases within the prethermal regime, including phases like the time crystal, which break the discrete time translation symmetry of the underlying drive. Third, let us emphasize that the presence of X withinŨ app f ensures that for every period, the system undergoes a nontrivial rotation that remains finite even in the high-frequency limit, ω → ∞; this corresponds to the remnant "Floquet structure" that remains within the prethermal regime. However, when one considers the evolution every N periods, one finds that the dynamics are simply generated by the static prethermal Hamiltonian D Ã : Finally, we emphasize that the emergent Z N symmetry is relevant only within the prethermal regime, where the total energy is also exponentially well conserved.
A. Prethermal emergent symmetry as a framework for nonequilibrium phases of matter In this section, we further elucidate the role of the emergent symmetry and how it provides a natural framework for realizing nonequilibrium phases of matter. Since the time evolution every N periods is captured by the prethermal Hamiltonian D Ã [Eq. (3)], there exists a timescale τ pre after which the system has "prethermalized" into a Gibbs state of D Ã and, thus, is locally described by ρ ∝ e −βD Ã , with a temperature β −1 determined by the system's initial energy density.
Let us now examine the evolution of this equilibrium state under a single period ofŨ app f . In general, ρ will evolve trivially because the equilibrium state respects the emergent symmetry X: However, if D Ã exhibits a spontaneously symmetry-broken (SSB) phase with respect to X, ρ can instead approach the equilibrium state within a particular symmetry-breaking sector; let us refer to such a spontaneously symmetrybroken state as ρ SSB . In this case, although ρ SSB evolves trivially under D Ã , the action of X is to rotate ρ SSB into a distinct symmetry-breaking sector ρ 0 SSB : During each period, the state rotates between the different symmetry-breaking sectors, only coming back to its original FRANCISCO MACHADO et al. PHYS. REV. X 10, 011043 (2020) 011043-4 sector after N periods (X N ¼ 1). The subharmonic nature of this behavior becomes transparent by measuring the order parameter, which is a local observable whose expectation value is different in each of the symmetry sectors.
In the language of time crystals, the fact that the underlying Floquet evolution has a period of T, while observables exhibit an enlarged periodicity NT, precisely corresponds to the discrete breaking of time-translation symmetry [28,31,[34][35][36]57]. For the remainder of this section, we continue to use the example of time crystalline order to highlight some of the unique features of prethermal nonequilibrium phases (Table I).
First, in order to meaningfully label the prethermal time crystal as a phase of matter, one needs to show that it remains stable under small perturbations. This is guaranteed so long as the discrete time-translation symmetry of the drive is not broken; in particular, this symmetry protects the emergent Z N symmetry, and we know that a phase that spontaneously breaks a Z N symmetry should be stable with respect to perturbations that do not explicitly break the symmetry.
Second, because our construction requires the system to prethermalize to an SSB state of D Ã , the observation of a prethermal time crystal depends on the choice of initial state (Table I). In particular, the initial energy density must correspond to a temperature below the critical temperature of the SSB phase transition. We emphasize that because the underlying transition of D Ã is sharp in the thermodynamic limit, there is an equally sharp transition between the prethermal time crystal and the trivial prethermal regime as a function of energy density (as long as τ Ã ≫ τ pre [73]).
Third, as the system begins absorbing energy from the drive at τ Ã , the temperature of the system will eventually cross the critical temperature of the SSB transition, leading to the loss of time crystalline order-the prethermal time crystal phase will always have a finite (but large) lifetime. To this end, depending on the energy density of the initial state, the lifetime of the time crystalline behavior can exhibit two distinct behaviors. If the energy density is below the critical SSB temperature, the system prethermalizes to the SSB phase and the timescale τ TC at which the time crystalline order parameter decays is similar to the heating timescale: τ TC ∼ τ Ã ∼ e ω=J local . If, on the other hand, the energy density is above the critical SSB temperature, the system will simply prethermalize to the symmetry preserving (trivial) phase and any transient time crystalline order can only occur before prethermalization, τ TC ≲ τ pre ∼ OðJ −1 local Þ. Differences between the many-body localized and prethermal discrete time crystal.-We end this section by juxtaposing the above discussions about the prethermal discrete time crystal with its many-body localized counterpart. Our focus is on highlighting the key differences between the two phases, as summarized in Table I. These differences can be divided into two categories: (1) the stability of the time crystal and (2) the restrictions on systems that can host a time crystal. Concerning the former, in contrast to the exponentially long lifetime of the PDTC, the ergodicity-breaking properties of Floquet manybody localization enable the MBL time crystal to persist to infinite times. Moreover, while the stability of the MBL time crystal can be independent of the initial state, the PDTC can only occur for a finite range of initial energy densities.
Let us now turn to the restrictions on systems that can realize an MBL versus a prethermal time crystal. In the MBL case, such systems are required to have strong disorder [74] and are unstable to the presence of an external bath [75], long-range interactions [76,77], and higher dimensions [77]. By contrast, the prethermal time crystal suffers from none of these restrictions and requires only two ingredients: a Floquet frequency that is larger than the local bandwidth and the existence of a static Hamiltonian D Ã with a spontaneously symmetry-broken phase. Crucially, in one dimension, this latter ingredient requires us to consider longrange interacting systems with power law 1 < α < 2 [78]; for such power laws, it is known that even a 1D system can exhibit a finite-temperature SSB phase, skirting the conventional Landau-Peierls argument that discrete symmetry breaking is forbidden for short-range interacting systems in 1D.

B. Prethermalization in long-range interacting systems
Before proving the existence of long-range interacting, prethermal phases of matter, we briefly contextualize a number of prior results with a particular emphasis on their implications for understanding the dynamics within the prethermal regime.
In particular, we now formalize the two different properties (for which we previously gave intuition) that U app f should satisfy in order to be of the broadest interest and most useful. We simplify the following discussion by focusing on the case without an emergent symmetry, Eq. (1), but our analysis carries over to the case with an emergent symmetry [Eq. (2)] by rotating into the frame U: to be a good approximation to U f , a naive first requirement is that the difference between the two unitaries be small. This can be encoded in a bound of the form where Λ is the volume of the system. Such a result would ensure that the error associated with the approximation in Eq. (1) is exponentially small in the frequency of the drive. However, owing to its volume dependence, this bound, at first, suggests that U app f is not meaningful in the thermodynamic limit, Λ → ∞. In particular, if one simply computes the overlap between wave functions evolved under the approximated and the true evolution, it would go to zero: But, of course, one is typically not interested in capturing the dynamics of the full quantum wave function (which cannot be measured), but rather in the dynamics of local observables. Unfortunately, by itself, Eq. (6) is insufficient to analyze the error in the evolution of generic local observables. Nevertheless, it can still be used to prove important results on the dynamics of extensive quasiconserved quantities. Of particular interest is the dynamics of the energy density D Ã =Λ. Since it remains constant under U app f , bounding the error growth of this observable provides an immediate upper bound on the heating rate under the true evolution.
To this end, by combining knowledge of the structure of the approximate unitary [Eq. (1)] with the error in the unitaries [Eq. (6)], one can immediately conclude that D Ã =Λ remains exponentially well conserved under the evolution: As promised, this formalizes the statement that the energy of the system is conserved up to an exponentially long timescale τ Ã and, thus, that the infinite temperature state cannot be reached before τ Ã . Note that for other extensive quantities conserved by D Ã , similar bounds can also be derived. (ii) Approximation of local dynamics. At this point, we have not yet formalized the statement that U app f is the correct "effective" generator of the true Floquet dynamics, only that the energy density remains conserved [79]. By filling in this gap, we would be able to rigorously connect the prethermal regime with the equilibrium properties of D Ã . This can be achieved by bounding the error in the dynamics of a generic local observable O as for some finite δ. Crucially, this result is independent of the volume of the system, meaning that it remains applicable even in the thermodynamic limit. This formalizes the intuition that, even if the global wave function is not perfectly captured by U app f [Eq. (7)], the local properties remain correct. Supplementing this result with an understanding of the equilibrium properties of D Ã as well as the structure of the unitary evolution (i.e., the emergent symmetry) will ultimately enable us to prove the existence of long-range, prethermal phases of matter.
Having formalized these two properties, we are now in a position to contextualize prior results on prethermalization in long-range interacting systems, without an emergent symmetry. In the case of an exponentially long thermalization time [property (i) above], the approximate unitary U app f has been proven to satisfy Eq. (8) for power laws α > d [38,39]. For approximating local dynamics [property (ii) above], the approximate unitary U app f has been proven to satisfy Eq. (9) for power laws α > 2d [38,39]. The discrepancy between these two regimes arises from the fact that Lieb-Robinson bounds with power-law light cones have been proven only for α > 2d [80][81][82][83]. When attempting to extrapolate to the case with an emergent symmetry in the prethermal regime, the above prior techniques do not appear readily generalizable [38,39].
Indeed, even for short-range interactions [36], generalizing to the case of an emergent symmetry requires the use of a different construction [40]. Curiously, although not explicitly discussed, many of the arguments found in this construction [40] generalize directly to the long-ranged case with little modification. In particular, the construction depends on the number of lattices sites each interaction term couples, which remains small even for long-range interactions (e.g., the long-range Ising interaction found in trapped ion experiments only couples pairs of sites [57]). As a result, one can directly use this construction for any power law α > d to create the approximate Floquet unitary U app f and to prove that it satisfies property (i), i.e., that it exhibits an exponentially long thermalization timescale. Extending to the case of an emergent symmetry then naturally follows by using the arguments found in Ref. [36].
Key challenge.-Unfortunately, since the construction found in Ref. [40] retains no spatial information about D Ã , one is unable to prove that U app f satisfies property (ii), i.e., that the dynamics of local observables are accurately captured.
Crucially, the lack of spatial information about D Ã prevents the application of Lieb-Robinson bounds, implying that any bound on the error of local observables diverges with the system size. To better understand the essential role of the Lieb-Robinson bounds, let us recall that the Floquet unitary is given by the exact expression [40]: where T denotes time ordering and V Ã ðtÞ is a timedependent interaction such that the sum of terms acting on any one site is exponentially small in frequency. One then builds the approximate unitary evolution U app f by disregarding the role of the exponentially small V Ã ðtÞ.
To understand how much error is accrued in this approximation, it is crucial to understand how a local operator O "spreads" under the evolution generated by D Ã . The bigger the volume of O, the larger the number of terms in V Ã ðtÞ it can overlap with and whose contribution we are missing when we disregard the role of V Ã ðtÞ. As such, the rate of error growth is simply bounded by the sum of the local terms of V Ã ðtÞ within the support Λ OðtÞ of the operator OðtÞ, while the total error δOðtÞ is the integral: δOðtÞ ∼ e −ω=J local R t 0 dt 0 Λ Oðt 0 Þ . The role of the interaction range is now apparent. If the original Floquet evolution is short-ranged, both the resulting D Ã and V Ã ðtÞ are also short-ranged and the evolution exhibits a finite Lieb-Robinson velocity v LR . The volume of the operator OðtÞ is then bounded by ∝ ðv LR tÞ d , and the error δOðtÞ ∼ t dþ1 e −ω=J local remains small for an exponentially long time in the frequency.
In contrast, when the original Floquet evolution is longrange, the volume of the operator O can grow much faster than Oðt d Þ. For example, for interactions decaying with power laws α ≤ 2d, only an exponential light cone has been proven, Λ OðtÞ ∼ e dηt [84]. In this case, the error δO ∼ e −ω=J local þdηt remains small for only a short time proportional to the frequency of the drive. For α > 2d, a power-law light cone has been proven [80][81][82][83], suggesting that if D Ã can be shown to exhibit an α > 2d spatial decay, one can immediately apply current Lieb-Robinson bounds. Of course, we hasten to remind the reader that in order to apply these long-range Lieb-Robinson bounds, one must first extend prior results (in the context of an emergent symmetry [36,40]) to determine the spatial decay of D Ã which, a priori, may be quite different from the decay of HðtÞ.
Prethermal phases in finite-size systems.-Up to now, our discussion has focused on the thermodynamic limit, where Lieb-Robinson bounds are required to prove that local dynamics are captured by U app f . However, in finite system sizes, Eq. (6) can actually be enough to guarantee that the prethermal Hamiltonian properly captures the dynamics. In particular, by setting the frequency of the drive large enough, i.e., ω ≫ log Λ, the approximate Floquet unitary is close to the full unitary evolution and the global wave function of the system is well approximated, regardless of the locality of the interactions. In this case, any observable (local or not) is well captured by the prethermal Hamiltonian until a timescale τ O ∼ Λ −1 e ω=J local (which remains smaller than the thermalization timescale τ Ã by a factor of Λ). Nevertheless, as long as τ pre is smaller than τ O , the system is guaranteed to approach the Gibbs state of D Ã , and this intermediate window (τ pre < t < τ O ) can host prethermal phases of matter.

C. Summary of key analytical results
Our main analytical results are twofold. First, we present a new construction for D Ã that explicitly retains information about the spatial locality of the interactions. Our construction naturally addresses the case where D Ã hosts an emergent Z N symmetry, extending prior results [36] to the case of long-range interactions. Second, using this novel construction, we are able to apply appropriate longrange Lieb-Robinson bounds to ensure that the prethermal Hamiltonian captures the local dynamics within the prethermal regime [property (ii)] and, thus, to prove the existence of long-range prethermal phases of matter.
For α > 2d, the existence of power-law light-cone Lieb-Robinson bounds allows us to prove that the local dynamics are accurately captured by U app f up to the Floquet heating timescale, τ Ã ∼ e ω=J local [third row of table in Fig. 1(c)]. This ensures that within the prethermal regime, the system will approach the equilibrium state of the prethermal Hamiltonian D Ã ; combined with the existence of an emergent symmetry (protected by the time-translation symmetry of the drive), this proves the existence of prethermal phases of matter [fourth row of table in Fig. 1 For d < α < 2d, we are not be able to directly invoke such power-law light-cone Lieb-Robinson bounds. In this case, the equilibration dynamics within the prethermal regime are less clear. Nevertheless, one expects that the approximate conservation of energy density means that local observables still relax to the Gibbs state of D Ã , since this is the state that maximizes the entropy subject to the constraint of conservation of energy. Under this assumption, we show that the robustness of prethermal phases of matter extends to power laws d < α < 2d as well [fourth row of table in Fig. 1(c), where the star indicates this additional assumption]. Moreover, in finite-size systems, one can prove rigorous statements without making this assumption, as discussed in the previous section.
In summary, our work demonstrates that prethermal phases of matter exist for all extensive power-law interacting systems (α > d).

III. RIGOROUS STATEMENT AND PROOF OF PRETHERMALIZATION RESULTS IN LONG-RANGE INTERACTING SYSTEMS
In this section, we describe our novel analytic construction, which extends prior results on prethermal phases [36,40] to the long-range interacting case. At its heart, this construction exactly transforms the initial time-dependent Hamiltonian into a new Hamiltonian composed of a static term D Ã (with an emergent Z N symmetry) in addition to small error terms. Crucially, this transformation captures two complementary properties. First, it ensures that the error terms are exponentially small in the frequency of the drive. Second, it guarantees that D Ã and the small error terms inherit the same locality properties as the original Hamiltonian; if the original Hamiltonian is long-ranged, the transformed Hamiltonian will also be long-ranged.
As discussed in Sec. II B, the first property allows us to prove an exponentially long thermalization timescale, in agreement with previous bounds [36,[38][39][40]. Meanwhile, the second property enables us to prove a much stronger statement, namely, that local observables remain well approximated by the long-range prethermal Hamiltonian throughout the prethermal regime (for power laws α > 2d)-a statement which has not been addressed in any prior literature for longrange interacting, prethermal systems with an emergent symmetry.
To guide the reader through this rather technical section, we present a short road map below. We begin by providing a careful treatment of previous results on prethermalization (Sec. III A). This introduces the necessary context to discuss the novel ideas required for our construction (Sec. III B). Next, we precisely state the key result of our construction in the form of Theorem 1 (Sec. III C). Finally, we discuss three immediate consequences of our construction (Sec. III D): (1) that local observables are well captured by the approximate Floquet unitary for α > 2d (Theorem 2), (2) how prethermal phases of matter arise even for α > d (Theorem 3), and (3) how our ideas can be directly generalizable to static systems with a near integer spectrum.

A. Previous results
Analyzing the Magnus expansion.-In Refs. [38,39], the main theoretical tool used to analyze the prethermal regime is the formal Magnus expansion of the single period time evolution operator U f . This procedure defines the Floquet Hamiltonian H F as a formal series expansion in the period of the drive T: with K m being operators and m the order of the Magnus expansion. Although such a series will, in general, not converge (otherwise there is a quasi-local Hamiltonian H F which is conserved under the dynamics of the system), understanding its truncation remains very useful. First, by truncating the Floquet Hamiltonian at the correct order, m¼0 T m K m , one obtains an exponentially good approximation to the full F . This implies that, over a single period of the evolution, the energy density hH ðn 0 Þ F i=Λ remains exponentially well conserved in the frequency of the drive; this corresponds to property (i) of Sec. II B. Because this analysis relies only on the few-bodyness of the interaction and the existence of a finite local energy scale, it holds for both short-and long-range interacting systems with α > d.
Second, for power laws α > 2d, one can use Lieb-Robinson bounds with power-law light cones [80][81][82][83] to prove that H ðn 0 Þ F is also the approximate generator of the dynamics of local observables for exponentially long times; this corresponds to property (ii) of Sec. II B. Combining these two conclusions, one proves the existence of a longlived prethermal regime whose dynamics are well captured by the prethermal Hamiltonian for short-and long-range interacting systems with power law α > 2d [first and second rows of the table in Fig. 1(c)]. Again, we emphasize that this construction does not prove the existence of an emergent symmetry in the prethermal regime; to obtain this result requires (to the best of our knowledge) a different approach.
Rotating into an appropriate frame.-To this end, a different approach [40] was pursued which enabled the proof of an emergent symmetry in the prethermal regime [36]. The main idea is to find a sequence of frame rotations where each rotation reduces the magnitude of the driven part of the evolution. Stopping the iteration at the correct step minimizes the driven component and proves the existence of a long-lived prethermal regime.
In more detail, one begins by separating the Hamiltonian HðtÞ ¼ H 0 ðtÞ into two components: a static D 0 and a driven V 0 ðtÞ term. Performing a rotation into a new frame, one obtains a new Hamiltonian H 1 ðtÞ that exactly describes the evolution, but where the norm of the driven term V 1 ðtÞ is reduced (while the static component D 1 is slightly modified); repeating such a process for n steps reduces the magnitude of the drive V n ðtÞ exponentially in n. However, much like the Magnus expansion result, this process cannot continue indefinitely or the system would be described by a static quasi-local Hamiltonian and thus fail to thermalize to the infinite-temperature state. The optimal iteration step is given by n Ã ∼ Oðω=ln 3 ωÞ, leading to the final Hamiltonian H n Ã ðtÞ: Since the local terms of the driven part V n Ã ðtÞ are exponentially small, the full evolution is approximately generated by the static component, Analogous to the Magnus expansion approach, one can prove that D n Ã =Λ remains exponentially well conserved under a single period: for some volume and frequency independent constant C; the thermalization timescale is then exponentially long in the frequency of the drive. Using this approach, one can also prove that the prethermal Hamiltonian can approximate the dynamics of local operators provided that the original evolution is governed by a Hamiltonian with short-range interactions. The source of this additional restriction is that, unlike the Magnus expansion approach, this construction cannot keep track of the range of interactions due to the way it accounts for the size of the Hamiltonian terms. More specifically, the proof ensures that any one operator does not grow to act on too many sites, without bounding the distance between the sites it acts on. In short-range interacting systems, this distinction is unimportant because the two measures of size are proportional; it is then guaranteed that D n Ã remains short-ranged and that the appropriate Lieb-Robinson bounds can be used to show it approximately generates the dynamics of local operators. However, this distinction becomes crucial in long-range interacting systems where these two measures can be very distinct leading to the breakdown of the proof, as explained in more detail in Sec. III B.
Generalizing to a prethermal emergent symmetry.-Understanding the limitations of this construction [40] is crucial because it provides the only path (to our knowledge) to prove the emergence of symmetries in the prethermal regime [36]. The main insight for this generalization is that the previous construction can be slightly modified to preserve the structure of the original Floquet unitary. Consider a Floquet unitary of the form where where E 0 corresponds to the static terms of the evolution that do not commute with the symmetry X. In this case, E 0 and V 0 ðtÞ are both the error terms we wish to minimize (in this language, the original construction corresponds to the specific case when N ¼ 1, X ¼ 1, and E 0 ¼ 0 [40]). To adapt their construction, one first rotates the system such that E 0 becomes time periodic, while keeping D 0 unchanged; the system is now fully characterized by D 0 and a new drive V 0 0 ðtÞ. One can now directly employ the previous construction to reduce the magnitude of the newly defined driven part [36]. The resulting new Hamiltonian contains terms E 1 and V 1 ðtÞ whose magnitude is reduced and a static D 1 whose magnitude slightly increases. Applying this procedure n Ã times reduces the size of E n Ã and V n Ã ðtÞ optimally, such that the unitary evolution is well approximated by the action of X and an evolution under the final static term D n Ã ¼ D Ã [Eq. (2)]. Let us emphasize that this picture is exact in a slightly rotated frame, U ≈ 1 þ Oðω −1 Þ, arising from the small rotation necessary to transform each E n into a driven term.
Because this analysis follows the results of Ref. [40], the results have the same scope with regards to the range of the interactions. In particular, the heating rate of the system is exponentially slow in frequency for both short-and longrange interactions with power law α > d; however, local observables are only provably well captured by the prethermal Hamiltonian in short-range interacting systems. Proving this result in full generality is the goal of the next few sections and will open up an entirely new landscape for investigating nonequilibrium phases of matter and their quantum simulation in long-range interacting quantum optical platforms.

B. Main ideas of proof for long-range generalization
In this section, we outline the novel ideas required to extend prior results [36,40] to long-range interacting systems; our main result is summarized in Theorems 1 and 2. For more details, see the Appendix B for the complete proof.
The main hurdle in generalizing the previous results to long-range interacting systems is to understand how the spatial structure of the interactions changes as one performs the necessary frame rotations.
We highlight, with a simplified example illustrated in Fig. 2, the importance of the range of interactions to the spread of operators. Although this example uses time evolution, the intuition carries over to the case of a frame rotation generated by some short-or long-range operator. Consider an operator O ¼ σ x i and a short-range interacting At early times, the spread of the operator is given by Crucially, the growth of the operator can happen only where it fails to commute with the Hamiltonian. Because the Hamiltonian is short-ranged, the range (spatial extent R) of the time-evolved operator is proportional to the size of the support of the operator (number of sites k it acts nontrivially on). This distinction may not seem meaningful for short-range interacting systems, but in long-range systems it becomes crucial. For example, if we consider long-range interactions such as then the spread of the operator is given instead by In this case, the time-evolved operator immediately becomes a sum of terms that connect two very distant points. While each term is two-bodied-i.e., the size of the support remains small with k ¼ 2-it can connect two points that are arbitrarily far away-i.e., the range R is arbitrarily large. We now connect this intuition to a careful analysis of the prethermal Hamiltonian. Starting from two-body interactions [such as Eq. (18)], the usual construction performs a rotation (informed by the driven part of the Hamiltonian) that generates a new Hamiltonian with higher-body and further extended terms [36,40]. To properly characterize the resulting final prethermal Hamiltonian, it is crucial to account for both the support size k and the spatial extent R of the terms, as these two properties play different roles in our result.
In particular, we need to ensure that terms that have either a large support size or a large range have a small magnitude. More precisely, if their magnitude decays exponentially with support size k, one can prove that there is a prethermal Hamiltonian exhibiting an exponentially long heating timescale. If their magnitude also decays with R with sufficiently large power law, one can employ the necessary Lieb-Robinson bounds to prove that the prethermal Hamiltonian is the approximate generator of the dynamics. In our work, we prove that this condition holds even when there is an emergent symmetry.
This latter point has eluded previous results [36,40] because their construction was unable to keep track of the spatial structure of interactions; in particular, a distinction is not made between an operator that acts on many sites (large k) and a few-body interaction that acts on sites far apart (large R).
To overcome this issue, our strategy is to imbue the construction with extra structure that enables us to keep track of the range and the size of the operator separately. To this end, we introduce the notion of an "R-ranged set" and use it to build "R-ranged operators." By representing the Hamiltonian in terms of R-ranged operators, we will ultimately be able to keep track of both the range R as well as the size k of the rotated Hamiltonian throughout the construction.
Let us begin by defining an R-ranged set. Schematically, an R-ranged set is a union of "clusters," each separated by distance at most R. As a result, any two of its sites are connected via a sequence of "jumps" of size at most R through the set, as shown in Fig. 2(c). Formalizing this picture, we define an R-ranged set as a set Z R of sites of our system, such that for x; x 0 ∈ Z R , there exists a sequence of elements ( At first sight, this definition appears more involved than simply characterizing a set based on its diameter (i.e., largest distance between two of its elements). This is on purpose. Indeed, our definition of an R-ranged set has the following crucial property: if two R-ranged sets have a nontrivial intersection, then their union is itself an R-ranged set. The same is not true for two sets with diameter at most R.
To see the importance of this property, let us first define an R-ranged operator as an operator whose support is an R-ranged set. The previous property of R-ranged sets immediately manifests in the following: if one takes two R-ranged operators A R 1 , B R 2 , then e A R 1 B R 2 e −A R 1 will be a maxðR 1 ; R 2 Þ-ranged operator. If we consider an operator written as a sum of R-ranged terms, then we can easily keep track of the range of each term as we perform a frame rotation (here, corresponding to e A R 1 ). When applied to the construction of the prethermal Hamiltonian, we can easily keep track of the R rangeness of each term of the original Hamiltonian throughout the different rotations.
The idea now is that we will consider potentials made up of a hierarchy of different-ranged interactions, decaying in an appropriate way with range. Specifically, we introduce a parameter σ > 0 (the value of which we will choose later), and define a sequences of ranges R l ¼ e σl . Then we will define a range-indexed potential to be a formal sum:  (19)]. In the short-range case (a), the operator remains close to its original location. For the operator to spread to a far away location, it requires many actions of H sr , which leads to a correspondingly large increase in its support; the range and support are closely related notions of size. In the long-range case (b), this need not be the case. The operator can very quickly spread across the system without a significant increase to its support; the range and the support of the operator capture very different notions of size. (c) An R-ranged set is a set where any two elements can be connected via a sequence of "jumps" (within the set) of size no greater than R. We illustrate this concept with the gray, green, and orange sets, each representing a different R-ranged set. Crucially, this definition is closed: when two R-ranged sets have a nonempty intersection, their union is also an R-ranged set (e.g., the gray and green sets). If they do not intersect, the union of two Rranged sets need not form an R-ranged set (e.g., the green and orange sets).
where Φ Z;l is supported on the R l -ranged set Z. Here we have introduced Z R l , the collection of all possible R l -ranged sets. Now we introduce a norm whose finiteness ensures our desired condition, namely, that the strength of the interactions decays exponentially in the size of their support k and as a power law in the range R. Specifically, we define a norm that depends on two parameters κ; γ > 0 according to where γ characterizes the power law of the long-range decay. This is a generalization of the norm used in Refs. [36,40], which did not keep track of the decay with range.
As an example, we note that for a two-body long-ranged Hamiltonian such as Eq. (18), our new norm Eq. (21) is finite in the thermodynamic limit provided that γ < α − d.
To see this, note that we can set where lðrÞ is the smallest l such that R l ≥ r. Then we have that On a d-dimensional lattice, we have X i;j∶ r<distði;jÞ≤r 0 for some constant C, and hence we find provided that γ < α − d.
However, we emphasize that our results also hold for Hamiltonians that are not just two-body. The only condition is that they decay fast enough with distance such that the norm in Eq. (21) is finite.
C. Statement of the prethermalization theorem for long-range interacting systems We have now set up all of the requisite tools. Our key contribution is developing the techniques required to analyze the range of the Hamiltonians produced by the aforementioned iterative construction, which leads to the following two main results (for details, see the Appendices).
First, we show that, by revisiting systems with shortrange interactions, we can obtain stronger bounds by simply replacing the particular sequence of numbers "κ n " chosen in Ref. [40] with a more optimized version. Second, by leveraging the properties of R-ranged operators and our particular choice of the sequence R l , we encode the information of the two-parameter norm Eq. (21), which captures the long-range nature of the interactions, back into the original one-parameter norm Eq. (22). This enables us to make use of the exact same analysis as in the short-range case, while keeping track of the long-range nature of the interactions via this encoding. Our final result is as follows.
If λ ≤ C 1 (the high-frequency regime), then there is a unitary transformation U which transforms the evolution to where LONG-RANGE PRETHERMAL PHASES OF NONEQUILIBRIUM … PHYS. REV. X 10, 011043 (2020) 011043-11 and Moreover, U is locality preserving and close to the identity in the following precise sense: for any range-indexed potential Φ.
D. Consequences of Theorem 1

Approximate form of the Floquet unitary
The end goal of Theorem 1 is to prove that the discussion in Sec. II A for realizing prethermal phases of matter (e.g., the prethermal time crystal) carries over to systems with power-law decaying interactions.
To this end, we build the approximate Floquet unitary evolution, It then immediately follows that property (i) from Sec. II B is satisfied: the energy density hD Ã i=Λ remains approximately conserved until the heating time τ Ã ∼ 2 n Ã . At this point, this just recovers an already obtainable result (even for long-range interactions) directly from the arguments of Ref. [36], albeit with an improved bound on the heating time since n Ã now lacks any logarithmic corrections in λ. Crucially, however, our choice of norm also guarantees that the interactions in D Ã [as well as E Ã and V Ã ðtÞ] remain power-law decaying in space. This allows us to consider how well U app f approximates the dynamics of local observables [property (ii) in Sec. II B], which requires the use of Lieb-Robinson bounds.

Approximation of local observables
As previously discussed in Sec. II B, proving that local dynamics are well captured by the prethermal Hamiltonian requires the existence of Lieb-Robinson bounds with power-law light cones. However, such bounds, in turn, require the prethermal Hamiltonian to exhibit the correct locality properties; its terms must decay, at most, as a power law of their range.
In our construction, this is guaranteed by the finiteness of our two-parameter norm [captured in Eqs. (34)- (36)], where the power-law decay of each term is characterized by the parameter γ Ã . Crucially, Theorem 1 guarantees that γ Ã can be chosen arbitrarily close to the parameter γ that characterizes the power-law decay of the original Hamiltonian of the system. This ensures that the prethermal Hamiltonian exhibits the same locality properties as the original Hamiltonian. Let us emphasize, however, that in the case where the original Hamiltonian contains two-body interactions, γ does not correspond to the exponent α that appears directly in the magnitude of each individual term [as in Eq. (18)]; rather, as we found in Eq. (27), γ must be smaller than α − d.
This language also enables us to immediately use Lieb-Robinson bounds available in the existing literature for multibody long-range interacting Hamiltonians [83]. In particular, as we show in Appendix C, any long-range interacting Hamiltonian H with bounded norm kHk κ;γ and γ > d satisfies the assumptions of Ref. [83], and therefore obeys a power-law light-cone Lieb-Robinson bound. We emphasize the requirement of a Lieb-Robinson bound for interactions with arbitrary k-bodyness, since our construction does not guarantee that the k-bodyness of the original Hamiltonian is preserved by the prethermal Hamiltonian.
Combining our knowledge of the locality of the prethermal Hamiltonian with the necessary Lieb-Robinson bounds, we prove the second main result of our work: all local observables are accurately captured by the approximate unitary U app f throughout the entire prethermal regime. This statement is formalized into the following theorem (see Appendix C for the proof).
Theorem 2. Approximation of local observables: Consider the scenario described in Theorem 1. Definẽ where U is the rotation constructed in Theorem 1, and define the corresponding approximate unitaryŨ app f ¼ Xe −iD Ã T by discarding the E Ã and V Ã terms in Eq. (33). Suppose that γ Ã > d, where d is the spatial dimension. Then for any η satisfying ðdþ1Þ=ðγ Ã þ1Þ< η<1, and for any local observable O supported on a set S, we have where τ ¼ ðC 6 λÞm, where C 6 is a constant that depends only on κ Ã and γ Ã , and C is a constant that depends only on the geometry of the system (but not its volume), the spatial dimension d, the size of the set S, and on η.
Before concluding this section, we hasten to emphasize that if novel multibody Lieb-Robinson bounds can be extended to power laws γ > 0, the construction presented in this work will immediately carry over. Such improvements would be in agreement with previous numerical and experimental results [85][86][87][88], as well as a recent proof for the particular case of two-body long-range interacting systems in one dimension [89].

Prethermal phases for power laws d < α < 2d
Unfortunately, we cannot prove a result as strong as Theorem 2 for 0 < γ Ã < d (corresponding to initial twobody Hamiltonians with d < α < 2d). Nevertheless, we can at least show that the dynamics of local observables are well approximated byŨ app f at short times (see Appendix D for the proof).
Theorem 3. Approximation of local observables (for short times): Consider the scenario described in where U is the rotation constructed in Theorem 1, and define the corresponding approximate unitaryŨ app f ¼ Xe −iD Ã T by discarding the E Ã and V Ã terms in Eq. (33). Then for any local observable O supported on a set S, we have, for any positive integer m satisfying mλ ≤ C 7 , where C 7 is a constant that depends only on κ Ã , and C 8 is a constant that depends only on κ Ã and the size of the set S. The assumptions of Theorem 3 differ from Theorem 2 in that Theorem 3 does not require γ Ã > d, but has an upper bound on the number of periods m which can be considered. For small enough λ (that is, high enough frequency), m max ¼ bC 7 =λc > 1, so one can at least accurately describe the dynamics of local observables during a single driving period.
The consequence of this result is as follows. Suppose that at some time t ¼ nT, the local observables are approximately described by the Gibbs ensemble of D Ã , or some spontaneous symmetry-broken sector thereof, which we call ρ. As mentioned in Sec. II C, we reemphasize that is a somewhat nontrivial assumption in the absence of a proof that the approximate unitary accurately describes the dynamics of local observables during the whole approach to thermal equilibrium; however, it follows if we assume that the system maximizes its entropy subject to the constraint of conserving energy density (which remains true for exponentially long times). Then, after one more driving period, the local state is approximately described by the rotated Gibbs ensemble,Ũ app f ρðŨ app f Þ † ¼ XρX † (using the fact that ½ρ; D Ã ¼ 0). This is all we need to repeat the arguments of Sec. II A about nonequilibrium prethermal phases of matter.

Extension to static systems
The long thermalization timescale of driven systems can also be generalized to static systems whose dynamics are dominated by an operator P with integer spectrum [36,40]: where ½D; P ¼ 0, while ½V; P ≠ 0 and u is the largest energy scale. In this setup, there is a change of frame where P becomes quasiconversed. To intuitively understand how this conservation emerges, it is simplest to consider a infinitesimal evolution under Δt ¼ δt=u: where the integer spectrum of P ensures that X ¼ e −iδtP with N ¼ 1=δt. However, we can make δt to be as small as possible, increasing the size of the emergent symmetry. In the δt → 0 limit, where Eq. (43) becomes exact, N → ∞ and we can think of the emergent symmetry as a continuous Uð1Þ symmetry, generated by the "number" operator P.
Analogously to the driven case, a time-independent change of frame U ensures that this emergent symmetry is approximately conserved until an exponentially long time in 1=u. This was proven in Ref. [40], closely following their techniques for driven systems. In a similar fashion, our construction immediately adapts to the proof of the longlived prethermal regime in static systems, allowing its extension to long-range interactions. As an application, we note that the existence of a prethermal continuous time crystal in an undriven system [36] can now be generalized to systems with long-range interactions.

IV. LONG-RANGE PRETHERMAL DISCRETE TIME CRYSTAL IN ONE DIMENSION
We now turn to the example of a nonequilibrium prethermal phase, where long-range interactions are essential to its stability-the disorder-free one-dimensional prethermal discrete time crystal. In particular, we study a one-dimensional periodically driven spin-1=2 chain with long-range interactions decaying with a power law d < α < 2d. Using massively parallel matrix-free Krylov methods [90][91][92][93], we compute the late time Floquet dynamics for system sizes up to L ¼ 28. This enables us to highlight many of the features of prethermal phases discussed in Sec. II A. First, by directly comparing shortand long-range interactions, we evince the crucial role of power-law interactions for stabilizing a 1D PDTC (Fig. 3). Second, by varying the energy density of the initial state, we access the aforementioned transition between the PDTC and the trivial phase (Fig. 4). These two phases can be easily distinguished by the different scaling behavior of the time crystal's lifetime τ TC : in the PDTC phase it follows the heating timescale. τ TC ∼ τ Ã ∼ e ω=J local , while in the trivial phase it is bounded by the prethermalization timescale, τ TC ≲ τ pre ∼ Oð1=J local Þ. We corroborate that our observed finite-size crossover matches the location of the phase transition independently computed via quantum Monte Carlo simulation of the corresponding equilibrium finite-temperature phase transition. These results provide insight into the experimental signatures of the PDTC, as well as direct measures of the relevant energy and timescales.

A. Model and probes
To generate Floquet dynamics that host a PDTC, the evolution must satisfy two properties: first, it must lead to a prethermal Hamiltonian D Ã with a robust emergent Z N symmetry, and second, D Ã must exhibit a spontaneous symmetry-breaking phase. We engineer a drive, motivated by current generation trapped ion experiments, that exhibits both.
To ensure that the emergent symmetry exists in the prethermal regime, we design a Floquet evolution that matches the form of Eq. (29)  and H x , for times T and T x , respectively. By choosing and σ ν i the Pauli operator acting on site i, the second part of the evolution flips all spins around thex direction (in the language of NMR, this portion of the evolution corresponds to a global π pulse): The resulting Floquet evolution then reads: matching Theorem 1, with N ¼ 2 and drive frequency ω ¼ 2π=T [95]. We emphasize that ½X; H ≠ 0; X is not a symmetry of the evolution. Next, to ensure that the associated prethermal Hamiltonian D Ã exhibits a spontaneous symmetry-breaking phase with respect to X, it must include long-range interactions with a power law d < α < 2d. However, D Ã results from the construction in Theorem 1 and thus corresponds to a complicated, frequency-dependent object. Fortunately, as part of Theorem 1 we saw that D Ã remains close (at high frequencies) to D, the original static symmetry respecting component of H, as defined by Eq. (29). Since H is time independent [Eq. (46)], D has a very simple form: it precisely contains the terms of H that are even under X. Thus, by including a long-range Ising interaction (which commutes with X) directly in H, one can guarantee that both D and D Ã exhibit a finite-temperature paramagnetic to ferromagnetic symmetry-breaking phase transition [78].
Combining the long-range Ising interaction with additional generic terms (that help to break integrability) leads to the following long-range Hamiltonian H: When we compare to the "short-range version" of this Floquet evolution, we will simply truncate the Ising interaction in H to nearest and next-nearest neighbor; we denote this corresponding short-range Hamiltonian as H s . For the remainder of this work we consider units where J ¼ 1 and use the parameters d < α ¼ 1.13 < 2d and fJ x ; h x ; h y ; h z g ¼ f0.75; 0.21; 0.17; 0.13g in a spin chain of size L with periodic boundary conditions [96]; we have verified that the observed phenomena are not sensitive to this particular choice of parameters. We note that, due to our choice of an antiferromagnetic coupling J > 0, the ferromagnetic phase occurs at the top of the spectrum of D Ã .
Finally, let us emphasize the role of the field term h x σ x i and nearest neighbor interactions J x σ x i σ x iþ1 to the thermalization properties of D Ã . While favoring the disordered phase, they also ensure that, to zeroth order in ω −1 , D Ã is not trivially diagonal and that, at large frequencies, the dynamics under D Ã are generic and thermalizing; as a result, both J x and h x control the timescale at which the system approaches the prethermal state τ pre .
Having described our model, we now introduce the diagnostics used to characterize its Floquet evolution. First, we consider the energy density of the system. Naively, one wishes to compute the energy density with respect to the full prethermal Hamiltonian D Ã ; however, its numerical construction and evaluation is very costly. Therefore, we will instead measure the energy density with respect to D, which remains close to D Ã at high frequencies. Second, we consider the half-chain entanglement entropy, S L=2 ¼ −Tr½ρ L=2 log ρ L=2 , where ρ L=2 ¼ Tr 1<i≤L=2 jψihψj, as a probe of the prethermalization and thermalization dynamics of the system.
To probe time crystalline behavior, we wish to consider an observable that can exhibit a subharmonic response to our driving protocol. From our discussion in Sec. II A, a suitable probe should be related to the order parameter of the paramagnetic to ferromagnetic transition in our model's prethermal Hamiltonian; for example, hσ z i ðtÞi for some site i. However, to reduce fluctuations owing to the small support of hσ z i ðtÞi, we find it convenient to average over the different sites of the system; let us then define It might have seemed more natural to consider instead the average magnetization σ z ðtÞ ¼ L −1 P L i¼1 hσ z i ðtÞi, but MðtÞ, which is related to a two-time correlation function, provides a clearer window into the early time decay of the period doubling behavior. Since we consider initial product states of σ z , Mðt ¼ 0Þ is guaranteed to be 1, its maximal value.
After the system prethermalizes to D Ã (for t > τ pre ), MðtÞ approaches a plateau whose sign will change every other period in the PDTC phase. Crucially, at this point and for translationally invariant systems (like our model), MðtÞ becomes proportional to the average magnetization σ z ðtÞ, which itself matches σ z i (for any i). As a result, MðtÞ is equally sensitive to the late time decay of the time crystalline behavior (provided that the initial magnetization is nonzero).
While MðtÞ is nonzero in the PDTC phase, it can also remains nonzero in the absence of a PDTC, e.g., in the ferromagnetic phase of a "static" Hamiltonian. The true order parameter for the PDTC phase must then measure the subharmonic (i.e., period doubling) response of MðtÞ. To this end, we introduce the PDTC order parameter: In the PDTC phase, MðtÞ will remain finite and sign changing every period and thus ΔMðtÞ will be nonzero. By contrast, in the symmetry-unbroken phase, all observables [including MðtÞ] quickly become T periodic and ΔMðtÞ approaches zero.

B. Exponentially long-lived PDTC
Before addressing the long-range PDTC, we begin by exploring the Floquet evolution of its short-range counterpart, H s , where previous results have proven the existence of an exponentially long-lived prethermal regime [36,[38][39][40][41]. As shown in Fig. 3(a), this is indeed borne out by the numerics: the energy density remains approximately constant until a late time τ Ã D Ã when hDi=L approaches its infinite temperature value of zero. By increasing the frequency of the drive, one observes an exponential increase in τ Ã D Ã , in agreement with analytic expectations [36,[38][39][40][41] and previous numerical studies [87]. These observations are mirrored in the evolution of the entanglement entropy S L=2 [ Fig. 3(d)]. There, the approach to the infinite temperature value, S T¼∞ L=2 ¼ ½ðL logð2Þ − 1=2 [94], In both the short-range model (g) and the "hot" long-range initial state (i), any period doubling behavior of the magnetization quickly decays as the system approaches, independently of frequency, the prethermal state atτ pre . By contrast, in the "cold" long-range initial state (h), the magnetization exhibits a robust period doubling behavior for as long as the energy density remains conserved; the decay of both quantities occurs at τ Ã ¼ Oðe ω=J local Þ and the prethermal time crystal is robust. This distinction is even clearer when considering the ω → ∞ limit of our Floquet evolution, where the magnetization shows no signs of decay.
occurs at τ Ã S L=2 , which is also exponentially controlled by the frequency of the drive. The agreement between τ Ã D Ã and τ Ã S L=2 corroborates the existence of a single thermalization timescale τ Ã that controls the approach to the infinitetemperature state. For the remainder of this work we quantify τ Ã using τ Ã S L=2 . Furthermore, S L=2 also informs us about the equilibration with respect to the prethermal Hamiltonian D Ã ; as the system evolves and approaches the prethermal state, the entanglement entropy approaches a plateau that remains constant until the drive begins heating the system at τ Ã . The timescale when S L=2 approaches this plateau value is frequency independent. In fact, the system's prethermalization is well captured by the ω → ∞ Floquet evolution (black dotted line in Fig. 3(d)). In this limit, we have , so the evolution for even periods is exactly generated by the static Hamiltonian D; for odd periods the wave function must be rotated by X (which does not affect S L=2 or hDi=L). This agreement with the ω → ∞ limit highlights that the dynamics within the prethermal regime are indeed well approximated by the prethermal Hamiltonian, D Ã ≈ D.
Finally, we turn to MðtÞ, our diagnostic for time crystalline order. From the discussion in Sec. II A, the lack of a spontaneous symmetry-breaking phase in shortrange interacting one-dimensional systems is expected to preclude the existence of the PDTC phase. In particular, any transient period doubling behavior should quickly decay as the system approaches the prethermal state at τ pre . This is precisely what is observed in the dynamics of MðtÞ, as shown in Fig. 3(g); while at very early times, even and odd periods exhibit almost opposite MðtÞ, by the timescale τ pre , MðtÞ has decayed to zero and the system no longer exhibits any time crystalline behavior. Thus, the transient signatures of a time crystal "melt" as the system equilibrates to the prethermal Hamiltonian D Ã , clearly demonstrating the system's lack of a true PDTC phase.
We now contrast this behavior to the long-range case using the same initial state, as evinced in Figs. 3(b), 3(e), and 3(h). With respect to the thermalization dynamicscaptured by hDi=L and S L=2 as shown in Figs. 3(b) and 3(e), respectively-the short-range and long-range models exhibit qualitative agreement; an increase in the frequency of the drive leads to an exponential increase in the thermalization timescale τ Ã . We note, however, an important quantitative difference. In particular, the value of J local extracted from the scaling τ Ã ∼ e ω=J local is larger in the long-range interacting system. This increase is due to the greater number of interactions terms in the Hamiltonian and is in agreement with previous numerical results [87]. In addition, τ pre remains frequency independent and the prethermal dynamics are in excellent agreement with the ω → ∞ time evolution [ Fig. 3(e)].
The difference between the short-and long-range interacting systems becomes apparent when considering the PDTC order. In particular, in the long-range model, the subharmonic response of MðtÞ survives well beyond τ pre and lasts until the heating timescale τ Ã . This behavior is robust. By increasing the frequency of the drive, the lifetime of the time crystal increases, mirroring the exponential growth of the thermalization timescale; the decay of time crystalline behavior is no longer determined by dynamics within the prethermal window, but rather by heating toward infinite temperature.

C. Role of the initial state
Another distinct feature of the PDTC is its sensitivity to the energy density of the initial state. Unlike the MBL time crystal [28,30,31,34,35,57], which can exhibit period doubling for all physically meaningful initial states, the stability of the prethermal time crystal relies on the prethermal state's approach to the symmetry-broken phase of D Ã . As a result, its stability is intimately related to the phase diagram of D Ã . Because hD Ã i=L remains approximately conserved until τ Ã , the energy density of the initial state is equal to the energy density of the prethermal state. With this in mind, one can then translate the initial energy density into the temperature β −1 of the prethermal state via the relation hD Ã ðt ¼ 0Þi ¼ Tr½D Ã e −βD Ã =Tr½e −βD Ã . By choosing initial states with different energy densities, one can effectively vary the temperature of the prethermal state across the phase transition; the resulting MðtÞ dynamics display qualitatively distinct behaviors.
This difference is manifest when we compare the dynamics of a "cold" state (near the top of the many-body spectrum [97]), Figs. 3(b), 3(e), and 3(h), with the dynamics of a "hot" state (near the center of the many-body spectrum), Figs. 3(c), 3(f), and 3(i). Despite exhibiting the same thermalization behavior to infinite temperature, the period doubling behavior of the hot state decays significantly faster; indeed, the decay of MðtÞ [and thus ΔMðtÞ] is frequency independent and occurs as the system approaches the prethermal state at t ≲ τ pre , well before the heating timescale τ Ã . This behavior is directly analogous to that of the shortrange model.
To directly connect the stability of the prethermal time crystal to the equilibrium phase diagram of D Ã , we study the decay timescale τ TC of the PDTC order parameter ΔMðtÞ across the spectrum of D Ã (Fig. 4). (For details on the extraction of these timescales, see Appendix F.) Crucially, τ TC exhibits important differences between the short-and long-range cases [Figs. 4(a) and 4(b) respectively]. In the short-range case, the frequency of the drive has no discernible effect on the lifetime of ΔMðtÞ (except for the highest energy state, which we discuss below).
In the long-range case, the behavior is significantly richer and modifying the driving frequency has a different effect depending on the energy density [ Fig. 4(b)]. The most distinct behaviors occur deep in the paramagnetic phase (near the center of the spectrum) and deep in the ferromagnetic phase (near the top of the spectrum). In the former, we observe the same frequency independent behavior of τ TC that characterized the short-range model-the decay timescale of ΔMðtÞ is simply determined by the prethermalization dynamics. In the latter, the behavior is dramatically distinct: τ TC increases exponentially with the drive frequency, following the thermalization timescale τ Ã ; in fact, the two timescales approach one another with increasing frequency-this is the key signature of the PDTC phase, namely, that the decay of the time crystalline order is limited only by the late time Floquet heating dynamics.
Having understood the behavior deep within each phase, we now turn to the transition between the two. At first glance, it appears that the onset of the exponential frequency scaling (and thus the transition to the PDTC phase) occurs at a lower energy density than what is expected [dark shaded region of Fig. 4(b)]. This expectation is based on an independent quantum Monte Carlo calculation for the transition in D (see Appendix H). As we explore below, this apparent inconsistency instead corresponds to a small finite frequency effect arising from the slow thermalization dynamics of D Ã near the phase transition, as schematically depicted in Fig. 5.
As a system approaches a phase transition, critical slowing-down causes its thermalization timescale to diverge. As a result, even in the paramagnetic phase, the decay of ΔMðtÞ can occur at very late times; we refer to this decay timescale as τ mag . In the paramagnetic phase, τ mag is finite, while in the ferromagnetic phase, it is infinite.
At low frequencies, if the system is near the phase transition on the paramagnetic side, τ mag can be finite but much larger than τ Ã . In this case, the decay of ΔMðtÞ is set by heating rather than the prethermal dynamics of D Ã even though the system is in the trivial phase. The situation is resolved upon increasing the frequency of the drive, at which point τ Ã and τ TC will both increase exponentially until they reach the magnetization decay time τ mag ; then, τ TC again becomes bounded by τ mag , losing its frequency dependence, while τ Ã keeps increasing exponentially with frequency. Thus, at large enough frequencies, it is always the case that, in the paramagnetic phase, the decay of ΔMðtÞ arises from the dynamics of D Ã .
This behavior is evinced in Fig. 4(b) in two distinct ways. First, by directly simulating the decay of ΔMðtÞ in the ω → ∞ limit (where heating cannot occur), we observe a significant increase of the decay time near the transition. In particular, in the paramagnetic phase, we observe a decay timescale which diverges around the transition at hDi=L ≈ 2.0-this is direct evidence for the presence of slow prethermalization dynamics near the transition. Second, near the transition to the ferromagnetic phase, the disagreement between τ mag (as measured by the decay of the magnetization in the ω → ∞ evolution) and τ TC occurs deeper in the trivial phase for smaller frequencies.
Interestingly, the above discussion also explains the long thermalization time found in the edgemost state of the short-range model, Fig. 4(a). In this case, the initial state is close to the zero temperature ferromagnetically ordered FIG. 5. Schematic explanation of the behavior near the transition of the long-range model (Fig. 4). There are two competing timescales: the heating time τ Ã and the magnetization decay time τ mag of the prethermal Hamiltonian D Ã [captured by the red squares in Figs. 4(a) and 4(b)]. As the system approaches the phase transition into the ferromagnetic phase (shaded region) from the paramagnetic side, τ mag diverges (red dashed line). The relaxation time τ TC is given by the smaller of these two timescales. In (most of) the paramagnetic phase, τ mag is smaller and approximately frequency independent; while in the ferromagnetic phase, τ Ã is smaller; τ TC shares its strong frequency dependence. , an analogous behavior occurs near the center of the spectrum. However, as one moves to higher energies across the paramagnetic to ferromagnetic phase transition (red shaded region), τ TC becomes exponentially dependent on the frequency of the drive and τ TC approaches τ Ã . In this regime, τ TC is set by the exponentially slow heating rather than the prethermal dynamics for all frequencies-the prethermal time crystal is stable.
state, leading to a finite, but very large prethermalization timescale. This very long prethermal equilibration time might also underlie the recent observations of long-lived period doubling behavior in the prethermal regime of shortrange interacting systems [98][99][100], where no finite-temperature phase transition or stable PDTC should occur.

V. CONCLUSION
Using a combination of analytical and numerical results, we demonstrate the existence of prethermal nonequilibrium phases of matter in long-range interacting systems with power laws α > d. This prethermal approach contrasts with recent MBL-based studies of Floquet phases, since it does not require disorder, nor is it limited by the dimensionality of the system. We emphasize the generality of our analytic construction, whose limitations arise only from the lack of an appropriate Lieb-Robinson bound for d < α < 2d. However, even in this regime, on quite general grounds, we expect the system to approach the Gibbs state with respect to the prethermal Hamiltonian and, thus, for prethermal phases of matter to be well defined. Finally, we predict the existence of a novel, disorder-free, prethermal discrete time crystal in one dimension. This phase is strictly forbidden in equilibrium, Floquet MBL, and shortrange interacting prethermal Floquet systems. Note added.-Recently, we became aware of a related complementary work on locality and heating in periodically driven, power-law interacting systems [101].

APPENDIX A: SHORT-RANGED PROOF
In this Appendix, we prove an improved version of the prethermalization theorem for short-ranged Hamiltonians. This improved version will eventually be the key to extending to the case of long-range power-law interactions.
Consider a finite set of sites Λ that characterize our system. Each site is assigned a finite Hilbert space, so the total Hilbert space becomes the tensor product of these local Hilbert spaces. One can then define any operator as a sum of terms acting on different parts of the system: where Q Z is an operator that acts on Z ⊆ Λ. The collection of Q Z is often referred to as a potential [40]. Despite this decomposition not being unique, our result constructs new potentials from an initial input potential, so this ambiguity does not affect our proof. We begin by introducing a one-parameter norm [40]: The finiteness of this norm in the limit of infinite volume indicates that the interactions are decaying exponentially with the size of their support. We can extend this definition to time-periodic potentials QðtÞ by considering the time average of the instantaneous norms: The statement of our theorem is as follows. Theorem 4: Suppose we have a time-periodic Hamiltonian HðtÞ ¼ Hðt þ TÞ which induces a Floquet evolution over a period T: such that D and E are time independent and Fix some κ 0 > 0, and define Now fix any 0 < C < 1. Then there exist constants C 1 ; …; C 5 > 0, depending only on C and κ 0 , with the following properties.
If λ ≤ C 1 (high-frequency regime), then there is a unitary transformation U which transforms the evolution to FRANCISCO MACHADO et al.
PHYS. REV. X 10, 011043 (2020) 011043-18 where and Moreover, U is locality preserving and close to the identity, in the following precise sense, for any potential Φ. Note that this is very similar to Theorem 1 of Ref. [36]. It differs, however, in two important ways. First, scaling of n Ã lacks the logarithm corrections with λ (which is proportional to the inverse frequency) found in Ref. [36]; as a result, the bound on the size of the residual "error" terms (V Ã and E Ã ) scales more stringently with frequency. Second, the norm k · k κ Ã with respect to which the final bounds are obtained has a parameter κ Ã which does not depend on λ. Roughly, the κ Ã for which a finite bound can be obtained can be thought of as setting an upper bound on the locality of the Hamiltonians, so the second condition means that D Ã , V Ã , and E Ã do not become more nonlocal as the frequency increases (whereas the theorems of Refs. [36,40] did not exclude this possibility).

Iteration
Following Ref. [36], the idea is to construct the necessary rotations iteratively. At step n of the iteration, there is a slightly rotated frame where the Floquet evolution operator U f is in the form with X N ¼ 1: We are interested in performing a unitary transformation, such that H n becomes closer to a time-independent term which commutes with the symmetry X. We begin by writing H n ðtÞ as the sum of two different contributions, D n and B n ðtÞ. D n corresponds to the time-independent part of H n ðtÞ which commutes with X-the "good" part-and it is given by where h·i T corresponds to the time averaging across a period: and h·i X corresponds to the symmetrization with respect to X, defined as Together, time averaging and symmetrization guarantee that D n is both time independent and commutes with X. B n ðtÞ is then the remaining "bad part" of H n ðtÞ and is composed of a time-independent term E n , which does not commute with X, and a time-dependent term V n ðtÞ: where V n ðtÞ is chosen such that At each step of the iteration we reduce the norm of B n ðtÞ by performing a transformation informed by H n . The construction for the iteration is exactly the one described in Ref. [36], and we do not repeat it here. We only differ from Ref. [36] in how we analyze the bounds satisfied by the iteration, as we describe in the next section.

Analysis of bounds
Now we prove bounds on the result of the iteration. Our first result is Lemma 1, a slightly modified form of Theorem 4 (Theorem 4 itself will eventually arise as a collorary), in which the constants are more explicitly stated.
Lemma 1: There are order 1 constants u and v (not depending on any other parameters) with the following properties.

011043-19
such that D and E are time independent and ½D; X ¼ 0: Fix some κ 0 > 0, and define Now fix any 0 < C < 1. Then suppose that Then there is a unitary transformation U which transforms the evolution to where and Moreover, U satisfies for any potential Φ.
Proof.-To prove Lemma 1, following Refs. [36,40], we introduce a decreasing sequence of numbers κ n > 0. The key difference between our analysis and that of Refs. [36,40] is in how we choose this sequence κ n . In particular, we choose this sequence in a way that is frequency dependent, meaning that it depends on the parameters λ and μ that appeared in the statement of the lemma. The higher the frequency (i.e., the smaller λ and μ), the slower κ n will decrease, which allows us to run the iteration to a larger order n Ã .
First of all, let us define dðnÞ ¼ kD n k κ n ; vðnÞ ¼ kV n k κ n ; eðnÞ ¼ kE n k κ n ; We recall the following bounds from Appendix A.4 of Ref. [36] (note that these bounds are independent of the choice of κ n ): where ε n ¼ 2TmðnÞv 0 ðnÞ½dðnÞ þ 2v 0 ðnÞ; ðA38Þ Note that there is an extra factor of 2 in Eq. (A38), which corrects an error [102] in Ref. [36]. These bounds hold provided that These results can be recast in a more intuitive manner as follows. Our eventual goal is to argue by induction. Suppose our induction hypothesis is that, given some h that is independent of the iteration order, vðnÞ; eðnÞ ≤ 1 Then we will make sure to choose κ nþ1 in terms of κ n such that the following conditions are satisfied: The point is that Eq. (A45), combined with Eq. (A43), ensures that Eq. (A41) is satisfied, and then Eq. (A44) combined with the induction hypothesis ensures that vðn þ 1Þ; eðn þ 1Þ; 2δdðnÞ ≤ 1 which, in turn, ensures that Eq. (A43), one of our induction hypotheses, is satisfied for n → n þ 1 (we consider the other one later). One way to ensure Eqs. (A44) and (A45) is to define for some ϵ > 0 that we will choose later. Then, where u < 1=2 and v are new constants introduced such that ffiffiffiffiffiffiffiffiffiffi ffi Computing explicitly for v, one obtains Equation (A44) is then satisfied provided that Meanwhile, for Eq. (A45) to be satisfied, we note that Therefore, Eq. (A45) is satisfied provided that In summary, the conditions on ϵ are that 6ðN þ 3Þ max 12 We choose to continue the iteration only while Accordingly, we will set ϵ ¼ bκ 2 0 ; then Eq. (A59) requires only that With this choice, we see that κ n ¼ κ 0 ffiffiffiffiffiffiffiffiffiffiffiffiffi 1 − bn p . Finally, we can complete the argument. The main missing piece is to show that the induction hypothesis Eq. (A42) is satisfied. Indeed, from Eq. (A46) we have that and, thus, Therefore, if we set h ¼ λ þ f½4ðN þ 3Þ þ 1=2gμ, then given the assumptions of Lemma 1, we can continue the induction up to the maximum iteration order n Ã . Finally, we need to prove Eq. (A35). From the form of the iteration (see Ref. [36]), we have where kA n k κ n ≤ NeðnÞT. Let us define Φ n ¼ e iA n Φ n−1 e −iA n , Φ 0 ¼ Φ. Then, from Lemma 4.1 of Ref. [40] and Eqs. (A43) and (A44), and the fact that h ≥ λ, we obtain and, hence, Then, we also have LONG-RANGE PRETHERMAL PHASES OF NONEQUILIBRIUM … PHYS. REV. X 10, 011043 (2020) 011043-21 ≤ e μ=2λ μ 4λ 1 2 n kΦk κ 0 ; ðA74Þ from which we conclude by summation and the triangle inequality that This completes the proof of Lemma 1. ▪ Now let us state how to prove Theorem 4. Lemma 1 (with μ ∼ λ) already takes us most of the way there, but it does not give the Oðλ 2 Þ scaling of kD − D Ã k κ Ã nor the OðλÞ scaling of kUΦU † − Φk κ Ã . The idea to fix this gap is that one should first do a single iteration of the procedure of Ref. [36], with κ 0 − κ 1 held fixed independently of λ (rather than the prescription above, for which κ 1 − κ 0 → 0 as λ → 0). In that case, we see from Eq. (A38) that ϵ 0 ¼ Oðλ 2 Þ. Now we apply Lemma 1 to the D 1 , V 1 , E 1 that result from the first iteration. We see that we can set the μ appearing in the statement of Lemma 1 to be Oðλ 2 Þ. Theorem 4 immediately follows.

APPENDIX B: PROOF OF THEOREM 1
In this Appendix, we prove our main theorem, Theorem 1 from Sec. III C. One of the principal ingredients is a new version of the prethermalization theorem for shortrange interactions, which we describe in Appendix A.
Here we extend this proof to range-indexed potentials, as introduced in the main text; recall that these are formal sums, where we have introduced a sequence R l ¼ e σl , and Z R l is the set of all R l -ranged subsets of sites (recall the definition of R-ranged set from Sec. III B). We define the formal commutator of two range-indexed potentials according to The idea is that we take the commutator of ½Φ Z 1 ;l 1 ; Θ Z 2 ;l 2 to be supported on Z 1 ∪ Z 2 , and then we observe that, in fact, if Z 1 and Z 2 are nondisjoint R l 1 -and R l 2 -ranged sets, respectively, then indeed Z 1 ∪ Z 2 is a maxfR l 1 ;R l 2 g¼R maxfl 1 ;l 2 granged set. This is true because an R 0 -ranged set is also an R-ranged set for R > R 0 , and the union of two nondisjoint R-ranged sets is also an R-ranged set. Then, we define the exponential action of one potential on another according to Recall from the main text that we introduced a twoparameter norm for range-indexed potentials, according to We will find it convenient to fix some κ 0 , γ and define a one-parameter norm for range-indexed potentials: We emphasize that this is not the same norm as Eq. (A2) for a potential Φ which does not keep any information regarding the range. Now we can prove the following key lemma. Lemma 2: Let Φ,Θ be range-indexed potentials, and let 0 < κ 0 < κ. Then, Proof.-This is analogous to Lemma 4.1 in Ref. [40]. Indeed, the proof carries through in exactly the same way, line by line, just replacing sums over Z with sums over (Z; l). The key fact for that proof was that for a collection of sets S 0 ; …; S m which is connected (i.e., it cannot be separated into nondisjoint subcollections), the size of their union P ¼∪ m k¼0 S k can be bounded by the sum of the sizes of each S i as For us, the analogous fact is as follows. Let S 0 ; …; S m be a connected collection of sets, and let l 0 ; …; l m ≥ 0. Then we have that ▪ In fact, Lemma 2 is already sufficient to allow us to extend Theorem 4 to range-indexed potentials. The reason is that the only two things we needed to prove Theorem 4 were the bounds Eq. (A36) and Lemma 4.1 of Ref. [40]. However, the only only nontrivial property of potentials that was used in deriving Eq. (A36) in Refs. [36,40] was Lemma 4.1 of Ref. [40] itself.
By generalizing Lemma 4.1 of Ref. [40] to Lemma 2 (which applied to range-indexed potentials), all of the argumentation in Theorem 1 from Sec. III C immediately carries over.

APPENDIX C: LIEB-ROBINSON BOUNDS FOR LONG-RANGE INTERACTIONS AND THE APPROXIMATION OF LOCAL OBSERVABLES
In this Appendix, we give the proof of Theorem 2 from Sec. III D 2.
We restrict our attention to sets of sites Λ that can be embedded in a Cartesian space R d , such that for any x ∈ Λ there exists r x ∈ R d such that distðx; yÞ ¼ jr x − r y j. We also assume that there is a smallest distance min x;y distðx; yÞ ¼ a, which we normalize to be 1.
The important result that we will use is that there is a Lieb-Robinson bound for time evolution by range-indexed potentials with bounded norm k · k κ;γ , so long as γ > d.
Lemma 3: Lieb-Robinson bounds for generic graded potentialsLet ΦðtÞ be a (time-dependent) graded potential with kΦk κ;γ < ∞ for some κ > 0 and γ > d. Let A be an operator supported on the set X ⊆ Λ, and let B be an operator supported on Y ⊆ Λ. Define the time evolution τ t ðAÞ as the time evolution of A according to ðd=dtÞτ t ðAÞ ¼ i½τ t ðAÞ; ΦðtÞ. Then for any η with ½ðdþ1Þ=ðγ þ1Þ< η<1, there is a Lieb-Robinson bound: and K 1 and K 2 are constants that depend only on the geometry of the system and on η, and we have defined Proof.-This is a corollary of Theorem 1 in Ref. [83]. To show that the theorem applies, we need only ensure that the assumptions of Sec. I of Ref. [83] are satisfied. First, observe that there is always a rescaling of time (which might be nonlinear) such that kΦðtÞk κ;γ becomes independent of t and equal to kΦk κ;γ . Now define Φ Z ¼ P ∞ l¼0 Φ Z;l (where we take Φ Z;l ¼ 0 if Z is not an R l -ranged set). Then we have, for any x ∈ Λ, s ∈ ½0; t: X where we used the fact that any R l -ranged set Z ∈Z R l satisfies diamðZÞ≤R l jZj, and the fact that max x¼½0;∞Þ e −κx ðκxÞ γ ¼ e −γ γ γ . Moreover, for any x ∈ Λ: X y∈Λ X Z∋x;y Hence, we see that the assumptions of Theorem 1 of Ref. [83] are satisfied with J ¼ e −γ ðγ=κÞ γ kΦk κ;γ ; ðC17Þ Therefore, the Lieb-Robinson bound follows from Ref. [83]. ▪ where in Eq. (C39) we used Bernoulli's inequality. Finally, we have X z f 3 ðdistðz; xÞ; tÞ ≤ K 4 ðτ þ τ β ÞjXj n Ã þ2 ; where which is finite in the thermodynamic limit provided ηγ > d. ▪ In this Appendix, we will deal only with potentials (not range-indexed potentials). Starting from a range-indexed potential we can construct a potential just by defining Φ Z ¼ P ∞ l¼0 Φ Z;l . We define the Heisenberg evolution of a (timeindependent) potential Θ by a (time-dependent) potential ΦðtÞ through the Dyson series for Heisenberg evolution, i.e., where ad Φ Θ ¼ ½Φ; Θ. This satisfies Our key result is as follows. Lemma 5: Consider numbers 0 < κ 0 < κ, and suppose that 3tkΦk κ 0 ≤ κ − κ 0 . Then, Here we defined Proof.-This is basically a time-dependent version of Lemma 4.1 from Ref. [40]. The proof proceeds in a nearly identical way. Indeed, we have where we defined kΦ Z k ¼ ð1=tÞ R t 0 kΦ Z ðtÞk. The rest of the proof proceeds identically to Lemma 4.1 of Ref. [40]. ▪ A corollary of this (or, in fact, of Lemma 4.1 of Ref. [40]) is as follows.
Lemma 6: For any potential W, we have Proof.-Just use the fact that ▪ Now we can prove a result about approximation of local observables.
Lemma 7: Define λ ¼ maxfkΦk κ ; kΦ 0 k κ g. Suppose that 12λt ≤ ðκ − κ 0 Þ. Then, where we defined ΔðtÞ ¼ ΦðtÞ − Φ 0 ðtÞ, kΦk κ ¼ ð1=tÞ R t 0 kΦðsÞk κ ds (and similarly for Φ 0 , Δ), and and therefore, ≤ CMkΔðsÞk κ 1 kE Φ ðsÞΘk κ 1 ðD16Þ where we have invoked Lemmas 5 and 6. This then gives Finally, we obtain where we invoked Lemma 5 once more. ▪ An immediate corollary is as follows. Lemma 8: Define λ ¼ maxfkΦk κ ; kΦ 0 k κ g. Suppose that 24λt ≤ κ. Let O be an observable supported on a set S. Then where we defined ΔðtÞ ¼ ΦðtÞ − Φ 0 ðtÞ, kΦk κ ¼ ð1=tÞ R t 0 kΦðsÞk κ ds (and similarly for Φ 0 , Δ), and Proof.-We define κ 0 ¼ κ=2 and treat O as a potential with a single term O S ¼O. Then kOk κ ¼e κjSj kOk. Moreover, we observe that δ ≔ E Φ O − E Φ 0 O, considered a potential, only takes nonzero values on sets Z that contain S. Therefore, given some s ∈ S, we have  Fig. 3 of the main text, we studied the late time Floquet dynamics of different initial states. The main feature that underlies much of our results is the existence of a long-lived prethermal plateau, where the system approaches an equilibrium state with respect to the prethermal Hamiltonian D Ã . In the main text, we studied the system's equilibration via the dynamics of energy density, entanglement entropy, and global magnetization (where the latter two exhibit long-lived plateaus consistent with the evolution under D, the zeroth term of D Ã ). In this Appendix, we supplement this analysis with the dynamics of local observables where we observe the approach of the dynamics to that of the prethermal Hamiltonian. Curiously, by studying the dynamics of the σ x operator, we observe evidence of the small, but finite, rotation of frame U that appeared in the statement of our theorem.
Our results are summarized in Fig. 6, where we consider the dynamics of σ z 4 , σ z 10 , σ x 4 , σ x 10 for the initial states considered in the main text, Fig. 3. We focus on the dynamics of even (full lines) and odd periods (thin dashed lines) independently in order to highlight any time crystalline behavior the local observables might possess (indeed this behavior is clear in the dynamics of σ z ). We also consider the time evolution in the ω → ∞ limit, where U f ¼ Xe −iDT (thin dashed line). This evolution enables us to see how well the full Floquet dynamics is captured by D Ã within the prethermal regime.
In particular, we wish to emphasize three different features in the dynamics of local observables. First, for the initial states that fail to approach the symmetry-broken prethermal phase, first and third column of Fig. 6, we observe that the dynamics of local observables under the Floquet evolution closely follows the dynamics of local observables under D until a late time approach to their infinite temperature value. By increasing the frequency of the drive, we observe this agreement extending to longer and longer times, emphasizing that D Ã is indeed the generator of the local dynamics of the system in the prethermal regime and that deviations occur due to the heating at a timescale τ Ã ∼ e ω=J local .
Second, this picture is not so clear when considering the initial state which approaches a symmetry-broken state in the prethermal regime, second column of Fig. 6. While the dynamics of σ z in this case are also very well described by D, the same is not true when considering σ x . We can attribute this to the effect of the small change of frame U; in the original lab frame, the system is really evolving under UD Ã U † rather than D Ã . Hence, measuring σ x in the lab frame is equivalent to measuring Uσ x U † in the rotated frame (where the evolution is governed by D Ã ). The latter has some overlap with σ z , which has large expectation value in the spontaneous symmetry-broken phase of D Ã (but zero expectation value in the symmetry-unbroken phase). Hence, since U is Oð1=ωÞ close to the identity, one finds that there is an Oð1=ωÞ contribution to the expectation of σ x in the lab frame, which disappears as ω → ∞, as can be observed in the numerics. [Note that other observables in principle could display the same effect, both inside and outside of the prethermal time crystal phase, but one can check by explicitly computing the perturbative expansion for U that the Oð1=ωÞ corrections happen to be much smaller in those cases.] These Oð1=ωÞ corrections also differ between odd and even periods (i.e., they exhibit time crystalline behavior), which is consistent with the picture that they arise from the overlap of Uσ x U † with σ z .
Finally, by comparing the dynamics of σ z 4 and σ z 10 , we can directly observe the local prethermalization of the system. In our choice of states, these two observables take opposite initial values, yet the translation invariance of our system implies that they must prethermalize to the same value. In particular, in the symmetry-broken phase, the thermal value of σ z is large, and so the sign of one of the local observables must change. Since the chain is mostly pointing up, σ z 10 , which started with a negative value, must prethermalize to a finite positive value, matching the magnetization of the remaining spins (including σ z 4 ). This is indeed what we observe, supporting the claim that the system approaches the prethermal state and that we are indeed observing the prethermal time crystalline phase.

APPENDIX F: EXTRACTION OF THE THERMALIZATION TIMESCALES
In order to better understand the thermalization dynamics of our Floquet evolution, we quantify the timescale at which different quantities approach their late time thermal values. In particular, we focus on the following quantities: the energy density of the system hDðtÞi=L, entanglement entropy S L=2 ðtÞ, time crystalline order parameter ΔMðtÞ, and the average magnetization in thex direction M x ðtÞ, where the latter is defined as We define the associated decay times as τ Ã D Ã , τ Ã S L=2 , τ Ã TC , and τ Ã S x , respectively. Although the complete dynamics of each quantity OðtÞ is nontrivial, at late times the system is in a local thermal state with respect to D Ã and their dynamics become much simpler. In particular, we observe that they exhibit an exponential approach to their infinite temperature value O T¼∞ : FIG. 6. Analysis of the evolution of different single spin operators-σ z 4 , σ z 10 , σ x 4 and σ x 10 -for the different conditions considered in Fig. 3: the short-range model (a),(d),(g),(j), a "cold" initial state in the long-range model (b),(e),(h),(k), and a "hot" initial state (e),(f),(i),(l). On the different single spin observables, we observe the approach to a position-independent constant within the prethermal regime, consistent with the plateau observed in the ω → ∞ limit Floquet evolution, further suggesting that the system has approached a thermal state of the prethermal Hamiltonian D Ã . By increasing the frequency of the driven system, we observe this agreement extending to later time, highlighting that the disagreement occurs due to the late time heating which becomes meaningful τ Ã ∼ e ω=J local . We also note that this simple picture is more complex in the case of σ x . In this case, one needs to account for the small frame rotation U which can induce a finite overlap between Uσ x U † and an observable that fails to commute with X.
Although this prescription is not exact and small deviations are observed, it provides a simple and robust way of extracting the thermalization timescale associated with each quantity.
This functional form motivates the following fitting procedure.
(i) We consider the evolution dynamics at every other period, so as to avoid any systematic effects of the period doubling behavior on the fits. The only observable where this effect is significant is thex magnetization M x ðtÞ (as discussed in Appendix E). Nevertheless, we observe that the extracted timescales are consistent regardless of the parity of the period considered. (ii) We restrict the data for the fit to the regime where jO T¼∞ − OðtÞj > ϵ for some small ϵ (ϵ ¼ 0.05 for energy density, ϵ ¼ 0.1 for time crystalline order parameter and entanglement entropy, and ϵ ¼ 0.015 forx magnetization). We found this cutoff necessary to ensure that the fitted curves captured the correct approach and were not dominated by the very small late time fluctuations close to the thermal value. (iii) We fit the linear relation y ¼ x=a þ b to log jO T¼∞ − OðtÞj as a function of t. The decay timescale is immediately given by the extracted value of a. (iv) Finally, we estimate the error of the procedure by partitioning the data in five regions and performing the same fitting procedure. The error is given by the weighted standard deviation of these results with respect to the global fit. Before moving on, let us note a small detail regarding the entropy timescale. Near infinite temperature β −1 , the entanglement entropy scales as β 2 as opposed to β like the other observables. As a result, to ensure that τ Ã S L=2 is capturing the same heating timescale τ Ã , the extracted value must be multiplied by a factor of 2 (for more details, see the Appendix of Ref. [103]).
Finally, the time evolution of the entanglement entropy also provides one more timescale: the time at which the system has approached the prethermal state τ pre .  7. Example of the fitting procedure for extracting the decay times for a particular initial state evolved with the long-range Floquet evolution. We apply the same procedure to all initial states in both the short-and the long-range model. We observe that a simple exponential decay captures the approach of different observables to their thermal values: (a) energy density hDðtÞi=L, (b) time crystalline order parameter ΔMðtÞ, (c)x magnetization M x ðtÞ (here plotted with a moving average over five points for clarity), and (d) half-chain entanglement entropy S L=2 ðtÞ. (e) The entanglement entropy provides an extra timescale τ pre which captures the approach to the prethermal state. The x axis in the shaded region is linear with time to emphasize the early time entanglement entropy behavior.
(f) Comparison of the different decay times. The decay time of the energy density τ Ã D Ã, entropy τ Ã S L=2 , andx magnetization τ Ã S x provide different estimates of the true thermalization timescale of the system τ Ã . Because this particular initial state is a "cold" state of the longrange model, it hosts a prethermal time crystal; the decay of the time crystalline order parameter also occurs at τ Ã . The agreement of all these timescales further corroborates the existence of a prethermal time crystal and the existence of a single thermalization timescale. Finally, we observe that τ pre occurs at a much earlier, frequency-independent timescale.
Unfortunately, the entropy dynamics are much more complex, so the above detailed fitting procedure does not apply. As a result, we follow a different procedure. Using the evolution of the initial state under the static Hamiltonian D, we obtain an approximation to the prethermal entanglement entropy value S pre L=2 , averaging the entanglement entropy value at late times. The time at which the driven system reaches 0.9S pre L=2 provides an estimate for τ pre . The error of this procedure is estimated by measuring the times at which the evolution reaches ð0.9 AE 0.05ÞS pre L=2 . We summarize both fitting procedures in Figs. 7 and 8, where we consider an initial state evolved under the longand the short-range model, respectively. The resulting decay times are plotted in the bottom right-hand panel, where we see agreement between all measures of the heating timescale τ Ã , as well as the existence of a much earlier, frequency-independent, decay time associated with the approach to the prethermal regime τ pre .

APPENDIX G: FURTHER EVIDENCE OF CRITICAL SLOWING-DOWN
As we approach the phase transition of D Ã from the paramagnetic side, we begin to observe the extension of the lifetime of the time crystalline order parameter, despite the system being in the trivial phase. This does not correspond to the breakdown of the prethermal phase, but rather extra physics in the equilibration dynamics under the prethermal Hamiltonian D Ã . In particular, this corresponds to the known phenomena of critical slowing-down. When one is close to the phase transition, small fluctuations in energy alter significantly the system's tendency to order or not; the system is unable to efficiently "choose" which side of the transition it actually is and equilibration takes a long time. This results in significant fluctuations in the dynamics and an enhancement in the timescale at which the system approaches the prehtermal state τ pre .
We can corroborate this hypothesis by investigating the dynamics of different initial product states evolving under the static Hamiltonian D. We focus on the entanglement entropy as its behavior has the simplest expectation; starting from zero, we expect the entanglement entropy to monotonically increase and approach a well-defined plateau corresponding to the equilibrium state. This is exactly what we observe for initial states far away from the phase transition, blue curves in Fig. 9. For initial states near the phase transition (on either side), red curves in Fig. 9, we observe a slower rate of entropy growth, plagued by much larger fluctuations. Moreover, these states also exhibit a very late approach to a well-defined plateau; some curves have yet to approach such a plateau although we are considering very late time dynamics, t ≳ 1000=J.  Fig. 7, but considering an initial state time evolved with the short-range Floquet evolution. As in Fig. 7, we observe that a simple exponential decay captures the broad features of the approach of the different quantities to their thermal values. Moreover, we also observe a good agreement between the τ Ã D Ã , τ Ã S L=2 , and τ Ã S x as measures of the thermalization time τ Ã . However, unlike the long-range case, the time crystalline order parameter (b) decays at a much faster, frequency-independent, timescale. This time is on the same order of τ pre , further corroborating that, in this case, the decay of the time crystalline order arises from the dynamics of the prethermal Hamiltonian.

APPENDIX H: QUANTUM MONTE CARLO CALCULATION
One of the requirements for a prethermal time crystal is a spontaneous symmetry-broken phase of the prethermal Hamiltonian; as long as the system thermalizes to a spontaneous symmetry-broken phase of D Ã , it will exhibit long-lived time crystalline behavior. As such, whether the system is in the prethermal time crystal phase is dependent on the temperature β −1 of the system as it prethermalizes to D Ã . In particular, as the system crosses the critical temperature T c , the system transitions from the prethermal time crystal phase to the prethermal trivial phase.
In order to estimate T c and by extension the critical energy density of the initial state ϵ c , we perform a quantum Monte Carlo simulation to understand the transition temperature of D Ã . Unfortunately, the full D Ã depends on the frequency of the drive. Fortunately, since we are working in the large frequency regime, we expect the transition to be dominated by the zeroth order term of D Ã , given by D: For ease of the numerical methods, for this analysis we invert our Hamiltonian by taking J to be negative, inverting the spectrum of the system. In this case, the bottom of the spectrum corresponds to the ferromagnetic ordered regime we observe at the top of the spectrum in the numerical calculations of Sec. IVof the main text. We note that h x and J x are kept positive to ensure that the Hamiltonian is signproblem-free. Since we expect the nature of the transition to be classical, we believe the difference of sign in these couplings does not significantly change the position or properties of the transition. In fact, when comparing our quantum Monte Carlo results to the classical model with J x ¼ h x ¼ 0, the location of the transition does not change; we believe flipping the sign of these couplings will not alter the stability and location of the phase.
To accommodate the periodic boundary condition of our problem, we modify the simple power-law behavior to the closest periodic function that describes a long-range decay, as it avoids any discontinuity in the derivative of the interaction.
For the case of this numerical investigation, we are interested in the finite-size crossover regime between the ferromagnetic and paramagnetic phases. This is of particular importance to correctly estimate the critical temperature, as long-range interacting systems often exhibit significant finite-size effects.
To diagnose the crossover, we make use of the heat capacity of the system which should present a divergence in the thermodynamic limit. In the finite system case, no true divergence occurs, but the presence of a peak in C V corresponds to a finite-size transition or crossover. The position of such a peak can then be used for estimating the critical temperature of the finite system T L¼22 C .
Using the information about the energy density of the system, as illustrated in Fig. 10, we numerically We observe that, for states away from the phase transition (blue lines), the evolution is characterized by a fast approach to a well-defined constant plateau. However, for initial states near the phase transition (red lines), the approach takes a very long time, displaying a slowly growing entropy for very long times and displaying large fluctuations. The initial states, marked in red, correspond to the state lying in the transition region in Fig. 4 of the main text.
differentiate the data with respect to temperature to obtain the heat capacity of the system. The position of the transition is then obtained by fitting the top of the peak in heat capacity to a Gaussian distribution. We estimate the uncertainty region associated with T L¼22 c as the region where the Gaussian distribution remains above 90% of its peak value (blue shaded region), leading to the estimation Finally, we can use the energy density curve to translate between critical temperature T L¼22 C and the critical energy density ϵ L¼22 C (red shaded region):