Discrete time-crystalline order enabled by quantum many-body scars: entanglement steering via periodic driving
Nishad Maskara, Alexios A Michailidis, Wen Wei Ho, Dolev Bluvstein, Soonwon Choi, Mikhail D Lukin, Maksym Serbyn
DDiscrete time-crystalline order enabled by quantum many-body scars: entanglementsteering via periodic driving
N. Maskara , A. A. Michailidis , W. W. Ho , , D. Bluvstein , S. Choi , , M. D. Lukin , and M. Serbyn Department of Physics, Harvard University, Cambridge, MA 02138, USA IST Austria, Am Campus 1, 3400 Klosterneuburg, Austria Department of Physics, Stanford University, Stanford, CA 94305, USA Department of Physics, University of California Berkeley, Berkeley, CA 94720, USA and Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA (Dated: March 1, 2021)The control of many-body quantum dynamics in complex systems is a key challenge in the questto reliably produce and manipulate large-scale quantum entangled states. Recently, quench experi-ments in Rydberg atom arrays (Bluvstein et al. Science , 25 Feb 2021) demonstrated that coherentrevivals associated with quantum many-body scars can be stabilized by periodic driving, generatingstable subharmonic responses over a wide parameter regime. We analyze a simple, related modelwhere these phenomena originate from spatiotemporal ordering in an effective Floquet unitary, cor-responding to discrete time-crystalline (DTC) behavior in a prethermal regime. Unlike conventionalDTC, the subharmonic response exists only for N´eel-like initial states, associated with quantumscars. We predict robustness to perturbations and identify emergent timescales that could be ob-served in future experiments. Our results suggest a route to controlling entanglement in interactingquantum systems by combining periodic driving with many-body scars.
Introduction. —Creating and manipulating entangle-ment is a fundamental goal of quantum information sci-ence, with broad implications in computation, metrology,and beyond. At the same time, not all forms of entan-glement are useful. In particular, strongly interactingquantum many-body systems generate large amounts ofentanglement under their intrinsic dynamics, in a processknown as thermalization [1, 2]. However, such dynamicsirreversibly scramble local quantum information, erasingmemory of the initial state. Creating and controlling en-tanglement while at the same time combating thermaliza-tion [3–5] in isolated interacting many-body systems [6–8]is therefore essential for applications of large-scale entan-gled states [9, 10].Experimental studies involving programmable quan-tum simulators based on Rydberg atom arrays [11] havesuggested that interacting quantum systems can exhibit aweak breakdown of thermalization, where certain initialconditions exhibit surprising, persistent many-body re-vivals. This phenomenon can be viewed as resulting fromso-called quantum many body scars (QMBS) [12, 13]–anomalous, non-thermal many-body eigenstates – namedin analogy to non-ergodic wavefunctions in the spectrumof otherwise chaotic single particle Hamiltonians [14]. In-triguingly, in some models with QMBS the system un-dergoes periodic entanglement and disentanglement cy-cles [13, 15–17], providing a potential route to the con-trolled manipulation of entanglement dynamics. In prac-tice, however, QMBS are fragile [12, 15, 18, 19]; sincethey rely on a dynamically disconnected subspace of non-thermalizing eigenstates [5, 15, 20], additional interac-tions generically lead to thermalization [19].Recent experiments [21] demonstrated that periodicdriving can dramatically increase the lifetime of scarredoscillations. This observation is surprising, since theexperiments used driving frequencies resonant with the local energy scales of the system, permitting easy en-ergy absorption and rapid heating towards a featureless,infinite-temperature state. Additionally, the experimentobserved a robust subharmonic response at half of thedriving frequency, suggestive of discrete time-crystalline(DTC) order [22, 23].In this Letter, we propose a theoretical framework forunderstanding these experimental observations by intro-ducing a mechanism whereby driving stabilizes quantumscarred oscillations, prolonging their lifetime and protect-ing them against arbitrary perturbations. Specifically, wefocus on the simplest model describing Rydberg blockade– the so-called PXP model [11, 24, 25] – which is an ideal-ized model for the Rydberg atom array experiment, withthe addition of kicked driving [21]. This model exhibitsrobust subharmonic responses and many-body revivalscoming from an effective many-body spin echo. The de-viation from a perfect echo introduces a small parameter,allowing us to derive an effective prethermal descriptionof the Floquet dynamics, and argue for stability up untilparametrically long times [26].Namely, we construct an effective Hamiltonian in a ro-tating frame, hosting an emergent Z symmetry which isspontaneously broken in its gapped ground state mani-fold. When viewed in the laboratory frame, the systemoscillates between the two spontaneously broken groundstates, resulting in a robust subharmonic response char-acteristic of DTC [26–28]. However, this subharmonic re-sponse is restricted only to N´eel-like initial states whichhave a strong overlap with the ground state of the ef-fective Hamiltonian — a property inherited from QMBS.Our model differs crucially from earlier works on homoge-nous time crystals in 1D [29–32] and mean-field construc-tions [33], in that the trajectory being stabilized is gener-ated by an interacting Hamiltonian, which produces non-trivial entanglement. Therefore, our construction opens a r X i v : . [ qu a n t - ph ] F e b a prospective route towards coherent control of entangle-ment dynamics. Model and phenomenology. —We study a periodicallykicked model H ( t ) = H PXP + θN (cid:80) k ∈ Z δ ( t − kτ ), whichgenerates the following one-period Floquet unitary, U F ( θ, τ ) = e − iθN e − iτH PXP , (1) H PXP = L (cid:88) i =1 P i − σ xi P i +1 , N = L (cid:88) i =1 n i , (2)describing evolution with the PXP Hamiltonian H PXP [11, 24, 25] for time τ , followed by the numberoperator N applied through rotation angle θ . Forsimplicity, the model is defined on a 1D chain of L siteswith periodic boundaries, although much of the analysiscarries over to higher dimensional bipartite lattices.Each site is a two-level system spanned by a ground( ◦ ) and an excited ( • ) state, and periodic boundaryconditions are assumed unless stated otherwise. Op-erators n i = |•(cid:105) (cid:104)•| i and P i = |◦(cid:105) (cid:104)◦| i project a givensite onto the excited and ground states respectively,while σ xi = |◦(cid:105) (cid:104)•| i + |•(cid:105) (cid:104)◦| i generates Rabi oscillations.In the Hamiltonian, σ xi is dressed by projectors onneighboring sites, ensuring that dynamics remain withinthe blockaded subspace where adjacent sites are neversimultaneously excited.For θ = 0 the Floquet dynamics (1) is equivalent toundriven evolution under H PXP . The PXP model isnon-integrable and features rapid growth of bipartite en-tanglement entropy, S ent ( t ) = − tr ρ ln ρ where ρ is thehalf-chain density matrix, from the majority of prod-uct states. In contrast, quenching from the N´eel state | Z (cid:105) = |•◦•◦ . . . (cid:105) leads to coherent oscillations between | Z (cid:105) and its inversion partner | Z (cid:48) (cid:105) , as first seen in [11],with oscillation period τ r ≈ . π that sets an intrin-sic resonant timescale. These oscillations are capturedby the imbalance in excitation number between odd andeven sites, I = (2 /L ) (cid:80) L/ i =1 ( n i − − n i ), see Fig. 1(a).However, dynamics under H PXP still generate entangle-ment, and the coherent many-body oscillations eventu-ally decay, see Fig. 1(a).The addition of strong driving with θ ≈ π almost com-pletely suppresses thermalization at early times, mostclearly seen in the nominal growth of entanglement en-tropy over multiple cycles, see Fig. 1(a). Concomitantly,oscillations of I synchronize to half the drive frequency,a phenomenon known as subharmonic locking. The ori-gin of this response is related to the existence of a spe-cial point at θ = π , because H PXP anticommutes withthe operator C = (cid:81) i σ zi = e − iπN , corresponding toa “particle-hole symmetry” [12]. Indeed, driving with θ = π implements an effective many-body echo, since U F ( π, τ ) = C e − iτH PXP C e − iτH PXP = e iτH PXP e − iτH PXP = ; this implies perfect subharmonic revivals across the entire Hilbert space. However, upon deviating from θ = π , we find that such revivals quickly damp outfor typical initial states without N´eel order, see Fig. 1and [34]. (a) − h I ( t ) i Evolution time t/τ r . . . S e n t ( t ) ( Z ) θ = ( Z ) θ = . π ( Z ) θ = 0 . π (b) Rotation angle θ/ π . . . S e n t . . . f ( ω d / ) | Z i | Z i Figure 1. (a) The density imbalance I and bipartiteentanglement entropy characterize oscillations between thetwo N´eel ordered states ( | Z (cid:105) ), colored orange (undriven)and black (driven), in an infinite size chain simulated viaiTEBD [34]. Adding driving with τ = 0 . τ r / θ = 0 . π to PXP model arrests the growth of entanglement entropy S ent and prolongs the lifetime. In contrast, driven dynam-ics from the | Z (cid:105) = |•◦◦◦•◦◦◦ . . . (cid:105) state thermalize rapidly(blue). (b) Subharmonic weight and average entanglemententropy, computed over 400 cycles ( T = 400 τ ), for an L = 28chain, over a range of θ and τ = 0 . τ r /
2. From the | Z (cid:105) state, these observables form a stable plateau around θ = π .However from the | Z (cid:105) state, which we use as a stand-in for ageneric state, the response quickly disappears away from theeffective echo point at θ = π . In contrast, long-lived oscillations from the N´eel statepersist over a wide range of parameters near θ = π and τ = τ r /
2. To quantify the stability of oscillationsand subharmonic response, we compute the subharmonicweight f ( ω d / ∝ | S ( ω d / | , defined as the normal-ized spectral weight of (cid:104)I ( t ) (cid:105) at half the driving fre-quency ω d = 2 π/τ , rescaled so f ( ω d /
2) = 1 for per-fect subharmonic response at θ = π from the N´eel states,see [34]. The plateaus in the subharmonic weight andtime-averaged entanglement entropy in Fig. 1(b) signal apersistent many-body response at frequency ω d / Many-body echo in su(2) subspace. —The robustness ofthe subharmonic oscillation away from θ = π can be qual-itatively understood in terms of mean-field-like trajecto-ries on an effective Bloch sphere, by invoking the forward-scattering approximation (FSA) introduced in [12, 35].The FSA constructs an L + 1 dimensional subspace thatcaptures dynamics under H PXP from a N´eel initial state,and approximately has the su(2) algebraic structure [12]of a spin- L/ S z op-erator is defined by the difference in excitation numberon odd and even sites, S z = (cid:80) L/ i =1 ( n i − − n i ), so theN´eel state | Z (cid:105) ( | Z (cid:48) (cid:105) ) corresponds to the North (South)pole. The S x operator is approximately proportional to H PXP , and generates a rotation around the x -axis thatexchanges the two N´eel states (blue lines in Fig. 2). Fi-nally, S y is calculated using su(2) commutation relations.Note that we use the weakly deformed PXP model [15]to generate the FSA basis, but consider dynamics under H PXP [34]. ●○●○●○●○○●○●○●○● ○●○●○●○●●○●○●○●○ (a) (b) e i⌧H PXP
2. (a) Thedynamics generated by two periods of U F ( π, τ ) exhibit a per-fect return to the | Z (cid:105) initial state: e − iτH PXP with τ = 0 . τ r under-rotates the | Z (cid:105) state (blue line), then the applicationof e − iπN (orange line) flips the x, y -projections of the spin sothat the second Floquet pulse completes the cycle. (b) Thesame dynamics but for θ = π − .
05 supports a periodic tra-jectory near the N´eel state. Dynamics initialized near theperiodic trajectory precess around it at stroboscopic timesforming cycles depicted for 100 driving periods for three ini-tial states (green ring corresponds to | Z (cid:105) initialization). Computing expectation values of the collective spin op-erators S x,y,z defined above, we visualize the many-bodydynamics from the N´eel initial state under two periodsof Floquet evolution (1) in Fig. 2. As mentioned, H PXP implements a rotation around the x -axis. In contrast,the action of e − iθN pulses is more complex, since the op-erator N does not have a closed form representation inthe su(2) subspace. However, it can be approximated as N ∼ ( S z ) in the vicinity of the N´eel states | Z (cid:105) and | Z (cid:48) (cid:105) ,which accumulate identical phases under e − iθN , see [34]and Fig. 2(b).Figure 2(a) illustrates that at θ = π the second ap-plication of H PXP returns the system to its initial state.Away from the perfect point θ = π , trajectories from | Z (cid:105) are no longer closed, but there exists a nearby closed or-bit with period 2 τ , see Fig. 2(b), explaining subharmonicresponse. In this picture, the existence of periodic trajec-tories is qualitatively similar to mean-field descriptionsof time crystals [33, 36, 37]. However, a key differenceis that the emergent spin- L/ L/ | Z (cid:105) is not closed,the stroboscopic dynamics (with period 2 τ ) exhibit pre-cession around the periodic trajectory, forming islands ofstability similar to Kolmogorov-Arnold-Moser tori in dy-namical systems. The precession leads to characteristicbeatings with an emergent timescale T b , correspondingto the period of motion on the circles around the fixedpoint in Fig. 2(b). The robust subharmonic responseand emergence of a beating timescale T b are a qualita- tive prediction of the spin- L/ L/ L + 1 dimensional subspace.To explain how driving reduces quantum thermalization,we consider the many-body Floquet unitary. Prethermal analysis and effective Hamiltonian. —Weanalyze the many-body dynamics by expanding aroundthe perfect echo point θ = π where the Floquet uni-tary is denoted X τ = U F ( π, τ ). This allows us to write U F ( θ, τ ) = e i(cid:15)N X τ , where (cid:15) = π − θ is a small parameterquantifying the deviation from the perfect point. Since X τ = , this unitary is in the canonical time crystalform [23, 27], which was rigorously analyzed by [26, 38],and we extend their results to the present case. In [34] weshow that the Floquet unitary can be approximated by U F ≈ V e − i(cid:15)H F X τ V † , where V is a unitary frame trans-formation perturbatively close to the identity and H F isan effective Hamiltonian also constructed perturbativelyin (cid:15) . As the DTC phenomenology depends on spectralproperties of the effective Floquet unitary, which are notaffected by V , we base our analysis on the leading ordereffective Hamiltonian and Floquet unitary, H (1) F = −
12 ( N + X τ N X τ ) , U (1) F ( θ, τ ) ≡ e − i(cid:15)H (1) F X τ . (3) H (1) F has an intuitive form, corresponding to the averageHamiltonian in a frame co-rotating with X [34]. Thus,at leading order, (cid:15) sets the timescale of dynamics in therotating frame.A key feature of this result is that the effective Hamil-tonian H F has an emergent Z symmetry [ H F , X τ ] = 0(even beyond the lowest order H (1) F , see [34]) guaran-teed as long as the time-periodicity of the drive is re-spected. H F can rigorously be shown to accurately de-scribe the system at least up to the prethermal timescale T p (cid:38) ( τ /(cid:15) ) e c p /(cid:15) for some constant c p > ε , defined by U (1) F | u (cid:105) = e iε | u (cid:105) . For τ nearan integer multiple of τ r /
2, the Floquet operator has apair of eigenstates characterized by strong overlap with | Z (cid:105) , | Z (cid:48) (cid:105) , and featuring nearly degenerate quasi-energies( τ = 0 , τ r ) or quasi-energies separated by π ( τ = τ r / | Z (cid:105) state. These observations imply the eigenstatescan be well approximated by the long-range correlated“cat” states |±(cid:105) = ( | Z (cid:105) ± | Z (cid:48) (cid:105) ) / √ k ( | + (cid:105) ) and k π ( |−(cid:105) ), and underliespontaneous symmetry breaking (SSB) of the system’stranslation symmetry. However, the emergent symme-tries X τ also play a crucial role, as the π and 0 quasi-energy gaps occur when X τ either exchanges the two N´eelstates ( τ = τ r /
2) or leaves them invariant ( τ = 0 , τ r ).At the level of the effective Hamiltonian H (1) F , these − . − . . . . Q u a s i e n e r g y ( ε / π ) π -paired (a) . . . . . Rotation angle τ/τ r . . . . . E − E g s (b) .
00 0 .
25 0 . N´eel overlap
15 20 25 L − − − E , + − E π , − k , X + k , X − k π , X − k π , X + Figure 3. (a) Eigenspectrum of U (1) F plotted for variousrotation angles τ , L = 16, and (cid:15) = 1, with color intensitythat corresponds to overlap with N´eel state. In the region of τ /τ r ≈ /
2, the states with largest overlap exhibit π -pairing,indicating subharmonic response. However, near τ /τ r ≈ , H (1) F re-veals double degeneracy between ground states from k = 0and π momentum sectors (denoted as k ,π ) with X τ eigenval-ues ± X ± ) in the region that corresponds to π -pairing in (a).The splitting of the ground state manifold vanishes exponen-tially with system size within the π -pairing region, with theinset showing finite size scaling at τ /τ r = 1 / π (0)-paired eigenstates correspond to degenerate groundstates in Fig. 3(b), separated by a finite gap ∆ to excitedstates, and belonging to different (same) symmetry sec-tors of X τ . Hence X τ symmetry breaking in the groundstate is linked to DTC order and the subharmonic oscil-lations of spatial order [23, 27, 34].We argue the observed region with DTC order de-scends from a model with conjectured perfect scars [15,34]. Specifically, if we deform the PXP model as de-scribed in [15], X τ at τ = τ r / |±(cid:105) become true ground states of H (1) F with a constant gap ∆ ≥
1. The PXP model, as well asdriving for τ away from τ r /
2, are weak deformations ofthis drive. However, these deformations do not preservethe emergent symmetry X τ at the level of H (1) F , and coulddestroy the ground state degeneracy. In Ref. [34], we ar-gue that since the emergent symmetry changes slowly aswe deform the drive, the ground states throughout the π -paired region in Fig. 3 can be considered as adiabat-ically connected to |±(cid:105) . Indeed, we confirm the energysplitting in the ground state of H (1) F decreases exponen-tially with system size, see Fig. 3(b) inset, as expectedfor SSB.The above analysis reveals four distinct timescalesemergent in the prethermal regime of Eq. (3). The short-est timescale T s = 2 τ is the subharmonic response. The second timescale, determined by the gap ∆ in the spec-trum of H (1) F , is T b ∝ τ ( (cid:15) ∆) − and comes from over-lap between the N´eel initial state and the lowest lyingexcited states. Semiclassically, T b is the precession pe-riod from Fig. 2(b). Finally, the longest timescale isset by the inverse energy splitting in the ground statemanifold of H (1) F , T g ∝ ( τ /(cid:15) ) e c d L , characteristic of SSB.All phenomenology is ultimately contingent upon thevalidity of the prethermal analysis, which holds until T p (cid:38) ( τ /(cid:15) ) e c p /(cid:15) . If such a bound is saturated and thesystem heats up to an infinite-temperature state beyond T p , the physics associated with T g will become unobserv-able, as for fixed (cid:15) and large enough system sizes T p < T g . Connections to experiments.—
We next demonstratethat the prethermal physics identified above persists be-yond the idealized Floquet model (1). Specifically, wereplace H PXP in Eq. (2) by the Rydberg Hamiltonian H Ry = (Ω / (cid:80) i σ xi − δ (cid:80) i n i + (cid:80) i ( V n i n i +1 + V n i n i +2 ),that includes imperfect Rydberg blockade and next-nearest-neighbor interactions. The PXP Hamiltonian isrecovered from H Ry in the limit V → ∞ , V = 0. Akinto the experiment in [21], we consider a 1D chain with V = V / , V = 10 Ω, and choose δ = V to cancel thestatic background from the next-nearest-neighbor inter-actions [21].Figure 4 illustrates the timescales T s , T b , and T g from stroboscopic dynamics of the revival fidelity F n = |(cid:104) Z | U F ( θ, τ ) n | Z (cid:105)| generated by the kicked Hamilto-nian H ( t ) = H Ry + θN (cid:80) k δ ( t − kτ ). At short times,on the order of tens of driving cycles, we observe a ro-bust subharmonic response at half the driving frequency,and an emergent beating timescale T b . After a few hun-dred driving cycles, the fidelity for even periods F n starts to decrease, while simultaneously for odd peri-ods F n +1 starts to increase. To understand this be-havior, we consider evolution at stroboscopic times andin the rotating frame, where the two nearly degenerate Time t/τ . . . F i d e li t y Time t/τ . . . time T p thermalization T s τ = T g ∝ τ e c g L T b ∝ τ ∆ Figure 4. Dynamics of revival fidelity under the period-ically kicked Rydberg Hamiltonian, and emergent prether-mal timescales. Stroboscopic dynamics of fidelity for θ =1 . π , and τ = 0 . τ r / T s , the beating timescale T b , and Rabi oscillations in thegroundspace, with characteristic timescale T g . Resonant time τ r ∝ / Ω depends on the Rabi frequency. Even (odd) mul-tiples of τ d are colored blue (red). Simulations performed onan L = 14 chain. ground states of H F , |±(cid:105) , form an effective two-level sys-tem with energy splitting ∆ E = E + − E − . The initalstate can be expanded as | Z (cid:105) = ( | + (cid:105) + |−(cid:105) ) √
2, andafter a time T g = π/ (2∆ E ), it evolves into a superpo-sition ( | + (cid:105) + i |−(cid:105) ) / √ | Z (cid:105) − i | Z (cid:48) (cid:105) ) / √ X τ kicks, which exchange the N´eel states ev-ery period. Finally, the prethermal time, when all fideli-ties might be expected to become exponentially small in L and all local observables relax, is not visible for thesystem sizes or times simulated. Discussion. —These considerations demonstrate thatentanglement dynamics associated with quantum many-body scars can be stabilized and steered in the peri-odically kicked PXP model, resulting in an evolutionstrongly reminiscent of prethermal DTC order. Our con-struction relies on the effective many-body π -pulse real-ized through quantum scars, which connect the two N´eelstates via an entangled trajectory, and a driving pulsethat reverses the direction of time. Similar to prethermaltime crystals, the emergent order features a robust, long-lived subharmonic response and spatiotemporal order fora range of parameters. However, an important differenceis that these signatures are present only for eigenstateswhich are perturbatively close to the N´eel initial state,and require sufficiently high fidelity state preparation tobe observed [30]. Nevertheless, we demonstrate that thesignatures of DTC physics survive in an experimentallyrelevant model, thus providing a possible explanation forrecent experimental observations in [21]. Moreover, wetheoretically predict new emergent timescales that couldbe observed in future experiments and suggested the pos-sibility of preparing entangled GHZ state [40] in drivenquench dynamics.The phenomenon described here drastically enhancesthe stability of non-ergodic dynamics thus opening alarge number of exciting directions. Specifically, by ex-tending this construction to the more complicated tra-jectories in the PXP model that connect highly en-tangled states [16] or to quantum scars in other mod- els [5, 17, 41, 42], control over complex entanglement dy-namics could be potentially implemented. From a prac-tical perspective, there remains a number of questionsrelated to experiments in Rydberg arrays [21]. In partic-ular, it is desirable to understand the dynamics in two-dimensional lattices [43], including the situations wheretwo sublattices have different numbers of nearest neigh-bours. In particular, in higher dimensions, there existsan intriguing possibility of realizing a true prethermaltime crystal, with a finite temperature phase transition in H F . It is also desirable to build a theory for higher ordersubharmonic responses observed in experiments [21], andobtain analytical understanding for continuously drivenmodels. Finally, it is important to understand if onecan implement full control over the many-body dynam-ics within the effective spin- L/ Acknowledgments. — We thank Dmitry Abanin, EhudAltman, Iris Cong, Sepehr Ebadi, Alex Keesling, HarryLevine, Ahmed Omran, Hannes Pichler, Rhine Samajdar,Guilia Semeghini, Tout Wang, and Norman Yao for stim-ulating discussions. We acknowledge support from theCenter for Ultracold Atoms, the National Science Foun-dation, the Vannevar Bush Faculty Fellowship, the U.S.Department of Energy, the Army Research Office MURI,and the DARPA ONISQ program (M.L., N.M, W.W.H.,D.B.); the European Research Council (ERC) under theEuropean Union’s Horizon 2020 Research and Innova-tion Programme Grant Agreement No. 850899 (A.M. andM.S.); the Department of Energy Computational Sci-ence Graduate Fellowship under Award Number(s) DE-SC0021110 (N.M.); the Moore Foundation EPiQS ini-tiative grant no. GBMF4306, the National University ofSingapore (NUS) Development Grant AY2019/2020 andthe Stanford Institute for Theoretical Physics (W.W.H.);the NSF Graduate Research Fellowship Program (grantDGE1745303) and The Fannie and John Hertz Founda-tion (D.B.); and the Miller Institute for Basic Researchin Science (S.C.). [1] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol,From quantum chaos and eigenstate thermalization tostatistical mechanics and thermodynamics, Advances inPhysics , 239 (2016).[2] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli,R. Schittko, P. M. Preiss, and M. Greiner, Quantum ther-malization through entanglement in an isolated many-body system, Science , 794 (2016).[3] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Col-loquium: Many-body localization, thermalization, andentanglement, Rev. Mod. Phys. , 021001 (2019).[4] R. Nandkishore and D. A. Huse, Many-body localiza-tion and thermalization in quantum statistical mechan- ics, Annu. Rev. Condens. Matter Phys. , 15 (2015).[5] M. Serbyn, D. A. Abanin, and Z. Papi´c, Quantum Many-Body Scars and Weak Breaking of Ergodicity, arXive-prints , arXiv:2011.09486 (2020), arXiv:2011.09486[quant-ph].[6] I. Bloch, J. Dalibard, and W. Zwerger, Many-bodyphysics with ultracold gases, Rev. Mod. Phys. , 885(2008).[7] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum sim-ulation, Rev. Mod. Phys. , 153 (2014).[8] A. Browaeys and T. Lahaye, Many-body physics with in-dividually controlled rydberg atoms, Nature Physics ,132 (2020). [9] R. Horodecki, P. Horodecki, M. Horodecki, andK. Horodecki, Quantum entanglement, Rev. Mod. Phys. , 865 (2009).[10] L. Pezz`e, A. Smerzi, M. K. Oberthaler, R. Schmied, andP. Treutlein, Quantum metrology with nonclassical statesof atomic ensembles, Rev. Mod. Phys. , 035005 (2018).[11] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Om-ran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres,M. Greiner, V. Vuletic, and M. D. Lukin, Probing many-body dynamics on a 51-atom quantum simulator, Nature , 579 (2017).[12] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn,and Z. Papi´c, Weak ergodicity breaking from quantummany-body scars, Nat. Phys. , 745 (2018).[13] W. W. Ho, S. Choi, H. Pichler, and M. D. Lukin, Pe-riodic orbits, entanglement, and quantum many-bodyscars in constrained models: Matrix product state ap-proach, Phys. Rev. Lett. , 040603 (2019).[14] E. J. Heller, Bound-state eigenfunctions of classicallychaotic hamiltonian systems: Scars of periodic orbits,Phys. Rev. Lett. , 1515 (1984).[15] S. Choi, C. J. Turner, H. Pichler, W. W. Ho, A. A.Michailidis, Z. Papi´c, M. Serbyn, M. D. Lukin, and D. A.Abanin, Emergent su(2) dynamics and perfect quantummany-body scars, Phys. Rev. Lett. , 220603 (2019).[16] A. A. Michailidis, C. J. Turner, Z. Papi´c, D. A. Abanin,and M. Serbyn, Slow quantum thermalization and many-body revivals from mixed phase space, Phys. Rev. X ,011055 (2020).[17] S. Chattopadhyay, H. Pichler, M. D. Lukin, and W. W.Ho, Quantum many-body scars from virtual entangledpairs, Phys. Rev. B , 174308 (2020).[18] V. Khemani, C. R. Laumann, and A. Chandran, Sig-natures of integrability in the dynamics of Rydberg-blockaded chains, Phys. Rev. B , 161101 (2019).[19] C.-J. Lin, A. Chandran, and O. I. Motrunich, Slow ther-malization of exact quantum many-body scar states un-der perturbations, Phys. Rev. Research , 033044 (2020).[20] N. Shiraishi and T. Mori, Systematic construction ofcounterexamples to the eigenstate thermalization hy-pothesis, Phys. Rev. Lett. , 030601 (2017).[21] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Se-meghini, S. Ebadi, T. T. Wang, A. A. Michailidis,N. Maskara, W. W. Ho, S. Choi, M. Serbyn, M. Greiner,V. Vuleti´c, and M. D. Lukin, Controlling quantum many-body dynamics in driven rydberg atom arrays, Science10.1126/science.abg2530 (2021).[22] V. Khemani, A. Lazarides, R. Moessner, and S. L.Sondhi, Phase structure of driven quantum systems,Phys. Rev. Lett. , 250401 (2016).[23] D. V. Else, B. Bauer, and C. Nayak, Floquet time crys-tals, Phys. Rev. Lett. , 090402 (2016).[24] P. Fendley, K. Sengupta, and S. Sachdev, Competingdensity-wave orders in a one-dimensional hard-bosonmodel, Phys. Rev. B , 075106 (2004).[25] I. Lesanovsky and H. Katsura, Interacting Fibonaccianyons in a Rydberg gas, Phys. Rev. A , 041601 (2012).[26] D. V. Else, B. Bauer, and C. Nayak, Prethermal phasesof matter protected by time-translation symmetry, Phys.Rev. X , 011026 (2017).[27] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi,Absolute stability and spatiotemporal long-range orderin floquet systems, Physical Review B , 10.1103/phys-revb.94.085112 (2016). [28] N. Yao, A. Potter, I.-D. Potirniche, and A. Vishwanath,Discrete time crystals: Rigidity, criticality, and real-izations, Physical Review Letters , 10.1103/phys-revlett.118.030401 (2017).[29] B. Huang, Y.-H. Wu, and W. V. Liu, Clean floquet timecrystals: Models and realizations in cold atoms, Physi-cal Review Letters , 10.1103/physrevlett.120.110603(2018).[30] A. Pizzi, D. Malz, G. De Tomasi, J. Knolle, and A. Nun-nenkamp, Time crystallinity and finite-size effects inclean floquet systems, Phys. Rev. B , 214207 (2020).[31] B. Mukherjee, S. Nandy, A. Sen, D. Sen, and K. Sen-gupta, Collapse and revival of quantum many-body scarsvia Floquet engineering, Phys. Rev. B , 245107(2020).[32] H. Yarloo, A. Emami Kopaei, and A. Langari, Homoge-neous floquet time crystal from weak ergodicity breaking,Phys. Rev. B , 224309 (2020).[33] A. Russomanno, F. Iemini, M. Dalmonte, and R. Fazio,Floquet time crystal in the lipkin-meshkov-glick model,Physical Review B , 10.1103/physrevb.95.214307(2017).[34] Supplemental online material.[35] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn,and Z. Papi´c, Quantum scarred eigenstates in a Rydbergatom chain: Entanglement, breakdown of thermalization,and stability to perturbations, Phys. Rev. B , 155134(2018).[36] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya,F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. vonKeyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin,Observation of discrete time-crystalline order in a disor-dered dipolar many-body system, ArXiv e-prints (2016),arXiv:1610.08057 [quant-ph].[37] W. W. Ho, S. Choi, M. D. Lukin, and D. A. Abanin,Critical time crystals in dipolar systems, Phys. Rev. Lett. , 010602 (2017).[38] D. V. Else, W. W. Ho, and P. T. Dumitrescu, Long-lived interacting phases of matter protected by multipletime-translation symmetries in quasiperiodically drivensystems, Phys. Rev. X , 021032 (2020).[39] D. Abanin, W. De Roeck, W. W. Ho, and F. Huve-neers, A rigorous theory of many-body prethermaliza-tion for periodically driven and closed quantum systems,Communications in Mathematical Physics , 809–827(2017).[40] A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T.Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler,S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Mon-tangero, T. Calarco, M. Endres, M. Greiner, V. Vuleti´c,and M. D. Lukin, Generation and manipulation ofSchr¨odinger cat states in Rydberg atom arrays, Science , 570 (2019).[41] K. Bull, I. Martin, and Z. Papi´c, Systematic construc-tion of scarred many-body dynamics in 1d lattice models,Phys. Rev. Lett. , 030601 (2019).[42] K. Mizuta, K. Takasan, and N. Kawakami, Exact Flo-quet quantum many-body scars under Rydberg blockade,Phys. Rev. Research , 033284 (2020).[43] A. A. Michailidis, C. J. Turner, Z. Papi´c, D. A. Abanin,and M. Serbyn, Stabilizing two-dimensional quantumscars by deformation and synchronization, Phys. Rev.Research , 022065 (2020).[44] J. Johansson, P. Nation, and F. Nori, Qutip 2: A python framework for the dynamics of open quan-tum systems, Computer Physics Communications ,1234–1240 (2013).[45] M. B. Hastings, Locality in quantum systems (2010),arXiv:1008.5137 [math-ph].[46] M. B. Hastings and X.-G. Wen, Quasiadiabatic contin-uation of quantum states: The stability of topologi-cal ground-state degeneracy and emergent gauge invari-ance, Physical Review B , 10.1103/physrevb.72.045141(2005). Supplementary material for “Discrete time-crystalline order enabled by quantummany-body scars: entanglement steering via periodic driving”
In this supplementary material we provide additionaldata for the phenomenology of Floquet model. In ad-dition, we describe the procedure used to visualize thedynamics on the Bloch sphere corresponding to collec-tive spin- L/ I. PHENOMENOLOGY OF DRIVEN MODELA. Quantifying the subharmonic weight
In the main text we introduced the subharmonic weight f ( ω ) used to quantify the subharmonic response of thedynamics generated by the Floquet unitary (1) startingfrom the | Z (cid:105) state, defined as f ( ω ) := δω | S ( ω ) | (cid:82) /t /T dω (cid:48) | S ( ω (cid:48) / | (S1)where t is the sampling rate, T is the sampling window, δω is a normalization, and S ( ω ) is the Fourier transformof (cid:104)I ( t ) (cid:105) with the mean subtracted: S ( ω ) := (cid:90) T dt ( (cid:104)I ( t ) (cid:105) − (cid:104)I ( t ) (cid:105) ) e − iωt (cid:104)I ( t ) (cid:105) := 1 T (cid:90) T dt (cid:104)I ( t ) (cid:105) . (S2)In the main text, we chose the normalization δω such thatthe exact subharmonic response in I at θ = π and τ = τ r / f ( ω ). We can decompose any bounded time-varying sig-nal in terms of its Fourier components. We will focus on (cid:104)I ( t ) (cid:105) , the signal given by the excitation imbalance intime. Rigorously, we can look at its so-called (one-sided)power spectrum, defined as P ( ω ) := lim T →∞ T | S ( ω ) | , (S3)which tells us how much ‘power’ (‘energy’ per unit time)is contained in the frequency interval [ ω, ω + dω ]. Notethat this quantity is related to the Fourier transform ofthe autocorrelation function of (cid:104)I ( t ) (cid:105) , by the Wiener-Khinchin theorem. Also, the integral over all frequenciesis the total power of the signal which is assumed finite,can be cast by Parseval’s theorem as (cid:90) ∞ dω P ( ω ) = lim T →∞ T (cid:90) T dt |(cid:104)I ( t ) (cid:105) − (cid:104)I ( t ) (cid:105)| . (S4) Thus, if we had infinite-time knowledge of the signal (cid:104)I ( t ) (cid:105) , we would define the fraction of its spectral weightin the frequency interval [ ω − δω/ , ω + δω/
2] as F ( ω, δω ) := (cid:82) ω + δω/ ω − δω/ dω (cid:48) P ( ω (cid:48) ) (cid:82) ∞ dω (cid:48) P ( ω (cid:48) ) . (S5)A value of F ( ω, δω ) being close to unity implies that theentire power of the signal is contained in the frequencyinterval [ ω − δω/ , ω + δω/ T is finite, since itis required to be smaller or comparable to the prether-mal timescale otherwise thermal fluctuations dominatethe response. In addition, we sample (cid:104)I ( t ) (cid:105) at discretetimes t n = nT /N , thus the integral in the denomina-tor of Eq. (S5) is effectively truncated at some high fre-quency cutoff ω c = 2 N/T set by the sampling rate, andlower cutoff ω = 1 /T set by the finite sampling window.Utilizing relation between P ( ω ) and S ( ω ) we write ourapproximation for F ( ω, δω ) as F ( ω, δω ) ≈ f ( ω ) = | S ( ω ) | δω (cid:82) ω c ω dω (cid:48) | S ( ω (cid:48) ) | . (S6)Finally, introducing discrete Fourier transform explicitly,we write the approximate subharmonic weight as: f ( ω ) = δω (cid:12)(cid:12)(cid:12) t (cid:80) Nn =0 e iωt n ( (cid:104)I ( t ) (cid:105) − (cid:104)I ( t ) (cid:105) ) (cid:12)(cid:12)(cid:12) t (cid:80) Nn =0 ( (cid:104)I ( t ) (cid:105) − (cid:104)I ( t ) (cid:105) ) . (S7) B. Stability to perturbations
In this section we demonstrate the stability againstsmall perturbations of short time dynamics generated bythe Floquet unitary U F ( π − (cid:15), τ ) as defined in Eq. (1)of the main text, in an infinite lattice. In addition, wehighlight the close relation between slowdown of ther-malization, as witnessed by entanglement growth, andthe strength of subharmonic response.To explore the stability of the subharmonic responsewe use the subharmonic weight at half the driving fre-quency, f ( ω d /
2) introduced in Eqs. (S1)-(S2). To quan-tify the rate of thermalization we compute the entangle-ment entropy for a bipartition of the chain S ent ( t ), andcalculate its time-average over c driving cycles,¯ S ent = 1 cτ (cid:90) cτ dt S ent ( t ) . (S8)Averages are computed over a time interval T = cτ , cho-sen to be an integer multiple of the driving period τ .Small values of ¯ S ent correspond to slow thermalization,whereas values of order one signal fast entanglement pro-duction.2 (a) (b) (c) ����������������������� ��� ��� ��� ��� ��������������������� ���������������������� ���� ���� ���� ���� ���������������������� ��������������������� ��� ��� ��� ��� ������������������������ Figure S1. Subharmonic response for the infinite chain remains strong in a broad range of parameters (top panels) and iscorrelated with small entanglement production rate (bottom panels) in the dynamics under Floquet unitary (1). (a) The broadrange of stability near the point (cid:15) = π − θ = 0 originates from the “perfect many-body echo”. The value of τ is fixed to τ r / τ r / θ = 0 . π .(c)Stability with respect to perturbation of the PXP Hamiltonian by next nearest neighbor interactions for θ = 0 . π and τ = τ r / T = cτ with c = 10 for (a) and (c), while c = 20 for (b). To this end, Fig. S1(a) studies the robustness of re-sponse to increasing values of (cid:15) . We observe that for (cid:15)/π < .
3, ¯ S ent remains constant with a small valuewhich primarily comes from entropy generated during themicromotion. For the same (cid:15) -interval the subharmonicresponse is almost maximal, indicating perfect revivalsof the local Rydberg excitations. Next, in Fig. S1(b)we explore the stability when the driving period τ ischanged compared to the intrinsic scar oscillation period τ r /
2. A robust subharmonic response accommodated bya strong suppression of the entropy growth is visible for τ ∈ [0 . , τ r /
2. The figure highlights how the subhar-monic response f ( ω d /
2) and average entropy ¯ S ent movein lockstep.Finally, in Figure S1(c) we demonstrate the broadrange of stability with respect to inclusion of next-nearest-neighbor interactions into PXP Hamiltonian, δH = V (cid:80) i n i n i +2 . Crucially, this perturbation doesnot anti-commute with the symmetry C = (cid:81) i σ zi , ruiningthe perfect many-body echo point even at (cid:15) = 0. Nev-ertheless, similar to the previous cases, we observe thatfor weak enough interactions, thermalization is stronglysuppressed and the subharmonic response is large. Thissupports our argument that the observed subharmonicresponse and slowdown of thermalization are closely re-lated and both are robust against generic weak pertur-bations. These infinite system simulations are in agree-ment and further support the results for the full RydbergHamiltonian on finite systems in the main text, which in-clude imperfect blockade and go out to longer times. C. Simulation Details
Entanglement and density imbalance for infinite chainsare calculated using iTEBD simulations. The initial state is fixed to be | Z (cid:105) product state. The data presented hereand in the main text is obtained using a third order Trot-ter integrator with time step ∆ t = τ / ε iTEBD ≤ − . ED simulations ofboth the PXP model and the full Rydberg Hamiltonian,used the Quantum Toolbox in Python (QuTiP) numeri-cal computing package [44]. II. VISUALIZING DYNAMICS IN SU(2)SUBSPACEA. Construction of the spin- L/ subspace The spin- L/ H ± to the N´eel state | Z (cid:105) (or | Z (cid:48) (cid:105) ). For the Hamiltonian (2), the ladder operators are H ± = L/ (cid:88) i =1 (cid:2) P i − σ ± i P i +1 + P i − σ ∓ i − P i (cid:3) . (S9)The original PXP Hamiltonian is obtained as a sum ofthese operators, H PXP = H + + H − . Repeated applica-tion of the ladder operator H + to the | Z (cid:105) state wouldgenerate the L + 1-dimensional basis of the forward scat-tering approximation. However, the subspace generatedin such way is only an approximate spin- L/ H ± , defined by dressing the pauli oper-ators ˜ σ ± i = σ ± i (cid:32) n max (cid:88) d =2 h d ( σ zi − d + σ zi + d ) (cid:33) (S10)3where h d = h (cid:0) φ ( d − − φ − ( d − (cid:1) − , h ≈ . φ = (1 + √ / n max = 8 in our numerical simulations. The coef-ficients h d decay exponentially with distance and sothis is a quasi-local deformation. As demonstrated inRef. [15], the subspace spanned by the unnormalized vec-tors | Z (cid:105) , ˜ H + | Z (cid:105) , . . . ( ˜ H + ) N − | Z (cid:105) , | Z (cid:48) (cid:105) gives rise to anumerically exact representation of su(2) with spin quan-tum number L/
2. We used this subspace to define spinoperators, but considered dynamics generated by the un-deformed PXP Hamiltonian.The operator S z = (cid:80) L/ i =1 ( n i − − n i ) acts as thespin- L/ S z operator in this space, assuming eigenval-ues − L/ , . . . , L/ L atoms [12, 15].The S x operator is given by ( ˜ H + + ˜ H − ) / (2˜ τ r ), where˜ τ r ≈ . τ r ≈ . π = 4 . τ r with τ r in numerical simulations. This can be done safely sincenone of the observed physics depends sensitively on τ r .Finally, the y -spin operator is obtained from the canon-ical commutation relation, S y = i [ S x , S z ]. Operators S z and S ± = ( S x + iS y ) / L + 1dimensional subspace. Specifically, let ˜ K be a projectoronto the subspace generated by repeated application ofthe deformed ladder operators on | Z (cid:105) . Then, the ladderoperators satisfy, up to numerical precision,˜ K [ S z , S ± ] ˜ K = ± ˜ K S ± ˜ K . (S11)By computing expectation values of these operators,we visualize vectors (cid:104) S α (cid:105) ( t ) as points on the Blochsphere, where N´eel states | Z (cid:105) ( | Z (cid:48) (cid:105) ) correspond to theNorth (South) pole, see Fig. 2. B. Action of N in the subspace In the spin- L/ N behaves as N ∼ ( S z ) in the vicinity of the poles. Toderive this relation, first notice that the operator N be-haves as an Ising interaction within the scar subspace.Lets choose the convention n i = (1 − σ zi ) /
2. Then, re-stricting our attention to the blockaded subspace where n i n i +1 = 0, we can rewrite N as N = L (cid:88) i =1 n i − L (cid:88) i =1 n i n i +1 = L − L (cid:88) i =1 σ zi σ zi +1 . (S12)A similar relation can be easily derived for PXP modelson any lattice. However lattices where sites have differ-ent coordination number, like in the imbalanced latticesstudied in [21, 43], have uncompensated (cid:80) i σ zi terms.Next, we treat the interaction in a mean-field-like ap-proximation. Taking into account the spatial homo-geneity and N´eel ordering, we can associate local den-sities with the global spin expectation value L (cid:104) S z (cid:105) = e i⌧H PXP
4, asthe e − iτH PXP pulse moves the north pole to the equator.The action of the driving pulse e − iθN is to move the statethrough the center of the bloch sphere. As a result, devi-ations from θ = π do not serve to stabilize this trajectory,but instead lead to thermalization, as can be seen by plot-ting the period 2 τ stroboscoipc dynamics. The timescaleof this thermalization is still parametrically slow, occur-ing at a rate set by (cid:15) , but this is qualitatively differentfrom the exponentially slow thermalization observed inthe gapped regime, with τ being close to an integer mul-tiple of τ r / III. EFFECTIVE HAMILTONIAN ANDROBUSTNESS TO PERTURBATIONSA. Derivation of effective description
We consider the time-periodic Hamiltonian H ( t ) = H PXP + θN (cid:88) k ∈ Z δ ( t − kτ ) , (S14)which generates the Floquet unitary U F = e i(cid:15)N X τ , (S15)where X τ = e − iπN e − iτH PXP has the property X τ = 1,and θ = π − (cid:15) . By grouping terms in the Hamiltonian,we can rewrite H ( t ) as H ( t ) = H PXP + πN (cid:88) k ∈ Z δ ( t − kτ ) (cid:124) (cid:123)(cid:122) (cid:125) H ( t ) − (cid:15)N (cid:88) k ∈ Z δ ( t − kτ ) . (S16) The first term H ( t ) generates the propagator U ( t ),which equals X τ after one period, and since X τ = 1 then U ( t ) is time-periodic with double the period U ( t ) = U ( t + 2 τ ). Thus, we can move into a rotating frame,with respect to U ( t ), and analyze the dynamics there.We assume the − (cid:15)N pulse is applied after the πN pulse,so that the Hamiltonian in the rotating frame reads H rot ( t ) = − (cid:15) (cid:104) X τ N X τ (cid:88) k ∈ Z δ ( t − (2 k − τ )+ N (cid:88) k ∈ Z δ ( t − kτ ) (cid:105) , (S17)which is equivalent to the Floquet unitary U F (2 τ ) = e i(cid:15)N e i(cid:15) X † τ N X τ . This Hamiltonian has a twisted time-translation symmetry H rot ( t ) = X τ H rot ( t + τ ) X τ , imme-diately ensuring that time-averaged Hamiltonian is sym-metric under X τ .However, as discussed in the main text, we can makea stronger statement, that the X τ symmetry holds toall orders in a perturbative expansion, using results fromRefs. [38, 39]. The aim is to derive an effective descriptionfor stroboscopic dynamics. To that end we can reparam-eterize the Hamiltonian in the rotating frame to makethe analysis easier, by introducing a new time parameter t (cid:48) defined over one fundamental region t (cid:48) ∈ [0 , (cid:15) ). Theunitary time evolution operator can be written as U ( t (cid:48) ) = U ( t (cid:48) ) U ( t (cid:48) ) , (S18)where U ( t (cid:48) ) is now 2 (cid:15) periodic and is U ( t (cid:48) ) = (cid:40) X τ for 0 ≤ t (cid:48) < (cid:15) (cid:15) ≤ t (cid:48) < (cid:15), (S19)and U ( t (cid:48) ) is generated by the time-dependent Hamilto-nian H ( t (cid:48) ) = (cid:40) −X τ N X τ for 0 ≤ t (cid:48) < (cid:15) − N for (cid:15) ≤ t (cid:48) < (cid:15) (S20)which can be written as H ( t (cid:48) ) = −
12 ( X τ N X τ + N ) −
12 ( X τ N X τ − N ) sgn (cid:20) sin (cid:18) t (cid:48) (cid:15) (cid:19)(cid:21) . (S21)The first term of H ( t (cid:48) ) is the time-independent averageHamiltonian, while the second term is an oscillatory termwhose time-average is 0.Our derivation of the effective description will be basedon the driven Hamiltonian H ( t (cid:48) ), but since the strobo-scopic dynamics are equivalent, the same (stroboscopic)description will apply to the Hamiltonian we first in-tended to analyze, Eq. (S14).The Hamiltonian H ( t (cid:48) ) has some key properties,which enable us to derive the effective Hamiltonian H F H ( t (cid:48) ) isa quasi-local Hamiltonian, it is a sum of bounded lo-cal terms whose amplitudes decay exponentially withsize of support of the terms. The term N = (cid:80) i n i isclearly local, and the norm of n i is one. For the secondterm, X τ N X τ = (cid:80) i X τ n i X τ , notice that X τ is gener-ated by evolution under H ( t ) for a finite time τ . Since H PXP has local bounded operator norm, and N is anon-site operator, X τ only spreads operators a finite dis-tance, with exponentially decaying support beyond thisdistance. This can be rigorously justified using Lieb-Robinson type bounds [45]. Second, H ( t (cid:48) ) still obeys atwisted time-translation symmetry H ( t (cid:48) + (cid:15) ) = X τ H ( t (cid:48) ) X τ . (S22)Therefore, invoking the theorem of Else et al. [38], to-gether with prethermalization bounds from Ref. [39], wecan derive that there is an effective description of U ( t (cid:48) )which reads U ( t (cid:48) ) ≈ U ( t (cid:48) ) V ( t (cid:48) ) e − i(cid:15)H F V (0) † , (S23)where V ( t (cid:48) ) is a 2 (cid:15) -periodic small change of frame (a uni-tary perturbatively close in (cid:15) to the identity, generatedby a quasi-local Hamiltonian), and H F is an effectiveHamiltoninan constructed perturbatively in (cid:15) .To the first couple orders in (cid:15) this reads H F = H + 12 (cid:88) n (cid:54) =0 (cid:15) [ H n , H − n ] nπ + · · · , (S24)where H n are the Fourier modes of H ( t (cid:48) ). The lead-ing order of this expansion, H coincides with the time-averaged Hamiltonian, and gives the expression for H (1) F introduced in the main text, Eq. (3).A key point to note is that H F has the symmetry[ H F , X τ ] = 0. To see this, notice that in Fourier space, H ( t (cid:48) ) = (cid:80) n H n e iπnt (cid:48) /(cid:15) , the twisted time-translationsymmetry reads X τ H n X τ = e iπn H n , and one can readilyverify H F satisfies X τ H F X τ = H F at least up to sec-ond order as given in Eq. (S24) (but we emphasize thesymmetry holds to all orders). Moreover, V ( t (cid:48) ) obeysalso the twisted time translational symmetry, V ( t (cid:48) + (cid:15) ) = X τ V ( t (cid:48) ) X τ . Therefore, our final expression for the effec-tive description of U F reads U F ≈ V (0) X τ e − i(cid:15)H F V (0) † . (S25)This expression implies that stroboscopically, up to asmall static frame change, dynamics are generated byan X τ -symmetric effective Hamiltonian H F interspersedby periodic kicks of X τ .This effective description lasts for at least time t (cid:48) ≥ e c p /(cid:15) , in the sense that evolution under the prethermalHamiltonian well approximates expectation values of lo-cal observables, so it describes the prethermal behavior,see Fig. S3. Since the original time t is related to t (cid:48) by Time t/τ . . . . . . S tr o b o s c o p i c R e v i v a l F i d e li t y Exact H (1) F Figure S3. Stroboscopic dynamics for L = 14, τ = 0 . τ r , (cid:15) = π/
10, plotted using exact floquet dynamics, and firstorder approximation. We see very good agreement up to longtimes. the factor ( τ /(cid:15) ), it means this description lasts in orig-inal time until t ≥ ( τ /(cid:15) ) e c p /(cid:15) . Additionally, all time-scales generated in the time-reparameterized setting aremultiplied by ( τ /(cid:15) ) to convert back to the original timeformulation. This includes the ground-state splitting os-cillation time. B. Robustness to perturbations
Now we argue that the time-crystalline behavior is ro-bust to arbitrary perturbations. Consider a generic time-periodic perturbation V ( t ) = V ( t + τ ) with a finite localoperator norm, and a small parameter δ V (cid:28) /τ . Theresulting perturbed Hamiltonian is given by H ( t ) = H ( t ) + H ( t ) (S26) H ( t ) = H P XP + πN (cid:88) k ∈ Z δ ( t − kτ ) (S27) H ( t ) = − (cid:15)N (cid:88) k ∈ Z δ ( t − kτ ) + δ V V ( t ) (S28)As before, the part H ( t ) generates a time-periodic frametransformation X τ . However now, the local operatornorm of H ( t ), integrated over one period, depends notonly on (cid:15) , but also δ V . We can again write the Floquetunitary as U ( τ ) = U ( τ ) U ( τ ) (S29)where U is given by U ( τ ) = T exp − i (cid:90) τ dt (cid:48) U ( t (cid:48) ) † H ( t (cid:48) ) U ( t (cid:48) ) (cid:124) (cid:123)(cid:122) (cid:125) ˜ H ( t (cid:48) ) (S30)6The Hamiltonian in the rotating frame ˜ H ( t (cid:48) ) is stillquasi-local and has a twisted-time translation symme-try, and therefore the effective Hamiltonian will still bedescribed by Eq. S24 and have an emergent symmetry.The only difference is that the small parameter, the in-tegral of the local norm of H ( t ) in the rotating frame,now goes as (cid:15) (cid:48) = (cid:15) + δ V τ . Hence, additional perturba-tions V ( t ) reduce the prethermal timescale to e c (cid:48) p /(cid:15) (cid:48) . Theunitary U ( τ ) can be approximated as U ( τ ) ≈ e − i(cid:15) (cid:48) H (cid:48) (S31)where to leading order in (cid:15) (cid:48) H (cid:48) = − (cid:15) (cid:15) (cid:48) ( N + X τ N X τ ) + δ V τ (cid:15) (cid:48) ( ¯ V + X τ ¯ V X τ ) , (S32)and ¯ V = 1 τ (cid:90) τ U ( t ) † V ( t ) U ( t ) dt. (S33)From Eq. (S32), we see that whether or not the groundstate of this new effective Hamiltonian is connected to theunperturbed one ( δ V = 0) depends on the ratio of δ V τ /(cid:15) (cid:48) .If there is a gap in the effective Hamiltonian at δ V =0, we expect the symmetry breaking to survive as longas δ V τ (cid:28) (cid:15) . This is why, perhaps counter-intuitively,some deviation (cid:15) from perfect π -pulse of N is required forstability of the subharmonic response against additionalarbitrary small perturbations δ V V ( t ). The necessity ofnon-zero (cid:15) to protect the ground state can be seen inFig. S6b, where the subharmonic response for the fullRydberg Hamiltonian disappears for a range of θ close to π , which corresponds to (cid:15) close to zero. C. Origin of the many-body gap
In this section, we present an argument for why themany-body gap observed in finite size numerics may per-sist in the thermodynamic limit, at τ = τ r /
2. Our anal-ysis is specialized to 1D, as it relies on the conjecture ofa deformed PXP Hamiltonian which exhibits a perfectsu(2) algebra and exchange of N´eel ordered states.We assume there exists a ‘perfect scar’ Hamiltoniangiven by ˜ H PXP = ˜ H + + ˜ H − from Eq. (S10). It can beviewed as a weak quasi-local deformation of original PXPHamiltonian, ˜ H PXP = H PXP + δV (S34)where δV = (cid:80) d ≥ h d P i − σ xi P i +1 σ zi + d are long-rangeterms whose amplitudes decay exponentially fast, see dis-cussion after Eq. (S10).Using the deformed PXP model we define a corre-sponding ˜ X τ operator,˜ X τ = e − iπN e − iτ ˜ H PXP , (S35) which still squares to one ˜ X τ = 1, and rewrite the Flo-quet unitary ˜ U F = e − i(cid:15)N ˜ X τ and effective Hamiltonian˜ H (1) F = − ( N + ˜ X τ N ˜ X τ ). For τ = ˜ τ r / |±(cid:105) = ( | Z (cid:105) ± | Z (cid:48) (cid:105) ) / √ H (1) F .Furthermore, as ˜ X τ exchanges the two N´eel states, the ˜ X τ symmetry is spontaneously broken in this ground space.These ground states are separated from the rest of thespectrum by a constant gap ∆ which is at least one. Con-sider first the − N term. | Z (cid:105) and | Z (cid:48) (cid:105) are the only stateswith L/ L/ − − N/ X τ is unitary, the spectrum of ˜ X τ N ˜ X τ is the same as N .Then, since | Z (cid:105) and | Z (cid:48) (cid:105) also span the two dimensionaleigenspace of ˜ X τ N ˜ X τ with eigenvalue L/
2, the opera-tor norm within the subspace spanned by all remainingstates is at most L/ −
1, and we can conclude that thegap in the second term is also one-half. Therefore, thegap above the ground space in ˜ H (1) F is lower bounded byone, i.e. ∆ ≥ H (1) F is close to ˜ H (1) F . However the two Hamiltonianshave different emergent symmetries, X τ and ˜ X τ respec-tively, so the ‘perturbation’ ˜ H (1) F − H (1) F does not respecteither symmetry. Nevertheless, using the argument insection III E, we can argue the ground space of H (1) F ex-hibits spontaneous symmetry breaking at τ = ˜ τ r /
2. Thisrequires three conditions, which we justify one by one.First, as we saw above, ˜ H (1) F has a gapped ground statewith SSB. Next, we show H (1) F is a deformation of ˜ H (1) F with finite local norm ≈ h = 0 . (cid:28) ∆ that appearssmall enough compared than the gap, that pairing stillappears in H (1) F . Recall H (1) F = − ( N + X τ N X τ ). Incontrast, ˜ H (1) F = − ( N + ˜ X τ N ˜ X τ ). Then,∆ H (1) F = − (cid:16) e − iτH PXP
N e iτH
PXP − e − iτ ˜ H PXP
N e iτ ˜ H PXP (cid:17) (S36)is the difference of the N operator time-evolved underthe two different unitaries. We expect that since thesetwo unitaries are close to one another, ∆ H (1) F should besmall. This is quantified by the Duhamel expansion: | ∆ H (1) F | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) τ dsU ( s ) † [ U † ( τ − s ) N U ( τ − s ) , δV ] U ( s ) (cid:12)(cid:12)(cid:12)(cid:12) , (S37)where U ( t ) = e − itH PXP and U ( t ) = e − it ˜ H PXP . Fromhere we see that, ∆ H (1) F is quasilocal Hamiltonian whoselocal norm is set by δV , and is small as long as h (cid:28) /τ .The final ingredient is to show that two symmetries X τ and ˜ X τ are also close. Specifically, we show the twosymmetries are related by a local unitary transforma-tion, X τ = P ˜ X τ P † , where P is generated by a quasi-local Hamiltonian with small local operator norm. Recall7that ˜ H PXP , H PXP , and δV are all anti-symmetric under C = e − iπN . Thus, we can write the two symmetries as˜ X τ = C e − iτ ˜ H PXP = e iτ ˜ H PXP / C e − iτ ˜ H PXP / (S38) X τ = C e − iτ ( ˜ H PXP − δV ) = e iτ ( ˜ H PXP − δV ) / C e − iτ ( ˜ H PXP − δV ) / Now, we can treat the unitary e − iτ ( ˜ H PXP − δV ) / by mov-ing into the rotating frame w.r.t. ˜ H PXP , specifically e − iτ ( ˜ H PXP − δV ) / = e − iτ ˜ H PXP / P † (S39) P † = T exp (cid:18) i (cid:90) τ dte it ˜ H PXP / δV e − it ˜ H PXP / (cid:19) The frame transformation P is indeed generated by aquasi-local Hamiltonian with small local operator normas long as h (cid:28) /τ , the same dependence as ∆ H (1) F .Then, the argument from Sec. III E apply, and we canconstruct the ground state manifold of H (1) F by an adia-batic deformation from ˜ H (1) F , with a spontaneously bro-ken broken symmetry X τ at τ = ˜ τ r /
2. Thus, thetwo groundstates of H (1) F should be dressed versions of( | Z (cid:105) ± | Z (cid:48) (cid:105) ) / √
2, with an energy splitting exponentiallysmall in L . Note that the difference between ˜ τ r and τ r is small enough that we neglected this technical detail inthe rest of the text, and simply worked at integer multi-ples of τ = τ r / D. Stability with respect to τ deformations In the prior section, we explained why symmetrybreaking occurs in the driven PXP model at τ = τ r / τ near τ r /
2. Consider the driven PXPmodel at two driving periods, τ = τ r / τ (cid:48) . If δτ = τ − τ (cid:48) is small parameter, then the difference be-tween the effective Hamiltonians and emergent Z sym-metries are both perturbatively small in δτ . Then, sincethe ground state at τ is gapped and spontaneously breaks X τ symmetry, we can use the arguments from sec. III Eto argue why the ground state at τ (cid:48) spontaneously breaks X τ (cid:48) symmetry.First, we show the two symmetries are related bya local unitary transformation. Indeed, it is sim-ple to see that X τ (cid:48) = e − iπN e − i ( τ − ( τ − τ (cid:48) )) H PXP = e − iδτH PXP / X τ e iδτH PXP / , where we used the fact that H PXP anti-commutes with C = e − iπN .Then, we need to show the two Hamiltonians H (1) F ( τ )and H (1) F ( τ (cid:48) ) are close. We see the difference H (1) F ( τ ) − H (1) F ( τ (cid:48) ) = 12 ( X τ N X τ − X τ (cid:48) N X τ (cid:48) ) (S40)= 12 e iτH PXP ( N − e − iδτH PXP
N e iδτH
PXP ) e − iτH PXP , is also perturbatively small in δτ using Eq. (S37). Asa result, the Floquet groundstates, corresponding to the ground states of H (1) F , are smoothly connected for nearby τ and remain nearly degenerate, in the regime with largespectral gap ∆. E. Why symmetry breaking survives
It is well known that if a quantum many-body Hamil-tonian H with a symmetry X , has gapped groundstates that spontaneously break X , then spontaneoussymmetry-breaking also persists for a nearby Hamilto-nian H (that is sufficiently close in local norm), as longas the symmetry is preserved, as then the Hamiltonians’ground states can be adiabatically connected [46]. Inthis section we present a slightly modified version ofthe statement: if two quantum many-body Hamil-tonians H , H are (i) perturbatively close, (ii) havesymmetries X and X respectively, which are alsoperturbatively close in a manner defined below, and(iii) H spontaneously breaks the symmetry X in itsgapped ground states, then H exhibits spontaneouslysymmetry-breaking of X in its ground states; moreover,its ground states are adiabatically connected to those of H ’s.We say that a unitary X is perturbatively close toanother unitary X if there exists a perturbativelysmall, quasilocal Hermitian operator S such that X = e iS X e − iS . Now, we can write H = e iS H e − iS + ( H − e iS H e − iS ) . (S41)This shows that H can be written as a deformation of e iS H e − iS . Crucially, by rotating H , we transformed its X symmetry into X , and now all terms are symmetricin X . As e iS is a unitary transformation which preservesspectral properties, the ground states of e iS H e − iS willspontaneously break the X symmetry by virtue of theassumption (iii) that H spontaneously breaks the X symmetry. Now all that remains to argue for is that thedeformation H − e iS H e − iS is perturbatively small sothat H ’s ground states can be constructed adiabaticallyfrom e iS H e − iS .To that end we note that H − e iS H e − iS = ( H − H ) + ( H − e iS H e − iS ) . (S42)By assumption (i) the first parenthesis on the right handside H − H is perturbatively small; by assumption (ii)the second parenthesis is also perturbatively small [thisfollows from using the Duhmael expansion (S38)], andthe statement follows. F. Spatiotemporal ordering in the ground state
Finally, to characterize the ordering in the ground state | ψ (cid:105) of H F via a local order parameter, we calculate a8 . . . . . Rotation angle τ/τ r . . . . . . D e n s i t y o r d e r i n g C n ( q , ω ) ( π, π, π ) , n T = 100 Figure S4. The Fourier transform of the density-density cor-relation function quantifies the order in the ground state of H (1) F , for L = 16 sites. The DTC phase has non-vanishingspatiotemporal ( π, π ) order. Fourier transform of the local density-density correlator C ( q, ω ) = 1 n T L n T (cid:88) n =1 L (cid:88) j =1 e i ( nω + jq ) (cid:104) ψ | n i U nF n i + j U − nF | ψ (cid:105) , (S43)in both spatial and temporal coordinates (restricted tostroboscopic times). The quantity C ( q, ω ) allows us toquantify different kinds of orders: the spatial (N´eel) or-der, where discrete space translation symmetry is spon-taneously broken to two-site translation symmetry, is di-agnosed by a large value of C ( π, C ( π, π ), which includes spontaneous breakingof discrete time translation symmetry. Figure S4 illus-trates that both orders are present in the range of τ where π -pairing of Floquet eigenstates persists. In contrast, invicinity of τ = 0 , τ r only spatial ordering is present, whilethe time translation symmetry remains intact. IV. SIMULATION OF RYDBERGHAMILTONIANA. Fidelity of GHZ state preparation
To quantify GHZ fidelity, we introduce a measurewhich does not depend on the phase of the GHZ state.Assume the quantum state is a pure state, and can bewritten as | ψ (cid:105) = c | Z (cid:105) + c | Z (cid:48) (cid:105) + ... . The GHZ statefidelity we compute in Fig. S5 is F = 12 ( | c | + | c | + 2 | c ∗ c | ) . (S44)When this fidelity F > .
5, the state is verifiably entan-gled [40].To understand the utility of this state preparation formetrology, we also compute the quanum Fisher informa-tion (QFI), with respect to the observable I . Since we fo-cus on pure states, this is simply four times the varianceof the observable, QF I = 4 (cid:104) (∆ I ) (cid:105) . In the ideal case, Time t/ø . . . . . . . G H Z F i d e li t y L 10121416 10 11 12 13 14 15 16
System size F i s h e r I n f o . . . . . G H Z fi d e li t y Time t/ø . . . . . . . G H Z F i d e li t y L 10121416 10 11 12 13 14 15 16
System size F i s h e r I n f o . . . . . G H Z fi d e li t y Figure S5. We look at pulsed driving of the Rydberg Hamil-tonian, with V = 10Ω, θ = 1 . π and τ = τ r /
2. These arethe same parameters of simulation as Fig. 4(a) in the maintext. (top) The GHZ fidelity, as defined in Eq. (S44), reaches amaximum at time T g /
4, which we see is strongly dependent onsystem size. Furthermore, we see (middle) the fidelity of theGHZ state preparation decays with system size. As a result,the quantum Fisher information, quantifying the sensitivityof the prepared state to an applied I field, does not exhibitquadratic Heisenberg scaling, but instead in the regime of in-terest appears best described by a linear scaling, indicatingthe standard quantum limit. However, as we further increasesystem size, this scaling may continue to change. where a GHZ state is prepared, the QFI scales quadrat-ically with the system size. This is known as Heisenberglimit scaling.However, the Rydberg Hamiltonian is a signficant per-turbation away from the ideal pulsed model. Therefore,the groundstate of H F is only perturbatively close to | Z (cid:105) , | Z (cid:48) (cid:105) . As a result, we see the fidelity of GHZ statepreparation drops noticeably with increasing system size.Indeed, we expect the fidelity to decay exponentially withsystem size, as typical for many-body overlap betweentwo quantum states. Furthermore, instead of Heisenberglimit scaling, the QFI exhibits a linear dependence, indi-cating standard quantum limit scaling.These preliminary results suggest that additional de-9 . . Response frequency ωτ / π . . . . . . P u l s e a n g l e − θ / π (a) − − | S ( ω ) | . . Response frequency ωτ / π . . . . . . . P u l s e a n g l e − θ / π (b) − − | S ( ω ) | . . Response frequency ωτ / π . . . . . . . P u l s e a n g l e − θ / π (c) − − | S ( ω ) | Figure S6. Fourier transforms of staggered density (cid:104)I ( t ) (cid:105) , for (a) the pulsed model with PXP Hamiltonian, (b) pulsed modelwith full Rydberg Hamiltonian, and (c) finite-width pulses with full Rydberg Hamiltonian. The subharmonic locking, andemergent beating timescale are very pronounced in the ideal model. However, the signatures near θ = π are robust, and survivefor the Rydberg Hamiltonian, even with finite-width pulses which more closely resembles the experiment [21]. Notice in (b)the subharmonic response disappears in a small region around θ ≈ π . We interpret this as occurring because the gap, setby (cid:15) = π − θ , is insufficiently small to protect against perturbations coming from perturbations to the ideal model. In (c),notice the significant offset of the line that corresponds to the beating timescale (labeled as 1 /T b ) from θ = π as ω →
0. Thislarge offset can be attributed to additional contributions to the effective Hamiltonian, introduced by finite-width pulses, seeEq. (S45). Finally, notice that for the full Ryderbg Hamiltonian, the dynamics are not equivalent under θ → − θ . Data in (a)is for τ = 0 . τ r , (b) corresponds to τ = 0 . τ r , (c) uses τ = 0 . τ r with pulse width τ p = 0 . τ . The system size is L = 24in (a) where we use constrained Hilbert space, whereas L = 14 for panels (b-c) where the full 2 L -dimensional Hilbert space isconsidered. velopments will be required before the time-crystallinebehavior described here can be used to prepare metro-logically useful states in experiments. Nevertheless, ourresults clearly show that for moderate system sizes, thequench dynamics reliably produce states with large over-lap with GHZ states. B. Emergent timescales for finite duration ofdetuning pulses
Finally, we consider an additional deformation of thepulsed driving, by applying the N pulse over a finiteperiod of time. Specifically, we use the following time-dependent Hamiltonian H ( t ) = H Ry + θN for nτ < t ≤ nτ + τ p / H Ry for nτ + τ p / < t ≤ ( n + 1) τ − τ p / H Ry + θN for ( n + 1) τ − τ p / < t ≤ ( n + 1) τ , (S45) where n = 1 , , . . . is a positive integer that correspondsto the current driving period. The system evolves withRydberg Hamiltonian at all times, though for a time τ p ithas an extra contribution from operator N . This drivingprofile is a better approximation for the cosine drivingused in Ref. [21]. Remarkably, we see in Fig. S6 thatdespite additional contributions to the effective Hamilto-nian H F coming from non-commutation of N and H Ry ,the subharmonic timescale T s and beating timescale T b are still clearly visible. Furthermore, T b still exhibits lin-ear dependence on θ − θ , although the base point θ isno longer centered at θ = ππ