Higher-order and fractional discrete time crystals in clean long-range interacting systems
HHigher-order and fractional discrete time crystals in clean long-range interacting systems
Andrea Pizzi, Johannes Knolle,
2, 3, 4 and Andreas Nunnenkamp Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United Kingdom Department of Physics, Technische Universit¨at M¨unchen, James-Franck-Straße 1, D-85748 Garching, Germany Munich Center for Quantum Science and Technology (MCQST), 80799 Munich, Germany Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom
Discrete time crystals are periodically driven systemscharacterized by a response with periodicity nT , with T the period of the drive and n > . Typically, n is an integerand bounded from above by the dimension of the local (orsingle particle) Hilbert space, the most prominent exam-ple being spin- / systems with n restricted to . Here weshow that a clean spin- / system in the presence of long-range interactions and transverse field can sustain a hugevariety of different ‘higher-order’ discrete time crystalswith integer and, surprisingly, even fractional n > . Wecharacterize these (arguably prethermal) non-equilibriumphases of matter thoroughly using a combination of exactdiagonalization, semiclassical methods, and spin-wave ap-proximations, which enable us to establish their stabilityin the presence of competing long- and short-range inter-actions. Remarkably, these phases emerge in a model withcontinous driving and time-independent interactions, con-venient for experimental implementations with ultracoldatoms or trapped ions.Introduction — Due to its foundational and technologicalrelevance, the study of condensed matter systems out of equi-librium has attracted growing interest in recent years, account-ing among others for the discovery of dynamical phase transi-tions [1, 2], quantum scars [3] and, particularly, discrete timecrystals (DTCs) [4–12]. A DTC is a non-equilibrium phaseof matter breaking the discrete time translational symmetry ofa periodic (i.e., Floquet) drive. In the thermodynamic limit,the defining feature of a n -DTC is a subharmonic response at /n -th of the drive frequency ( n > ), which is robust to per-turbations and which persists up to infinite time [12]. Follow-ing the first seminal proposals [4–8], DTCs have been widelyinvestigated both theoretically and experimentally [12–23].In this context, most work has focused on spin- / sys-tems, which have largely been shown to exhibit a -DTCwhere at every Floquet period each spin (approximately) os-cillates between the states |↑(cid:105) and |↓(cid:105) leading to period dou-bling ( n = 2 ). This fact naturally emerges from the dimension of the local Hilbert space of the spins [6], and can be gener-alized to n -DTCs in models of n -dimensional clocks [24, 25].Another well-studied setting is that of bosons in a gravita-tional field bouncing on an oscillating mirror [9], where thesingle-particle Hilbert space dimension is infinite (as the parti-cle’s position is continuous) and where n -DTCs with arbitraryinteger [26] and fractional [27] n have been shown.In these systems, heating to a featureless ‘infinite tempera-ture’ state is typically avoided by introducing disorder, whichleads to a (Floquet) many-body-localized (MBL) phase [5, 8]. FIG. 1.
Higher-order and fractional discrete time crystals. (a) Aspin- / chain with long-range interactions and initial z -polarizationis driven with a monochromatic transverse magnetic field of strength h , inducing a spin precession around x . (b) The time crystallinity isprobed by the Fourier transform | ˜ m ( ν ) | of the magnetization along z . The spectrum fragments in a multitude of plateaus with constantfrequency /n for a magnetic field strength h in a finite range ≈ /n ,each of which signals a higher-order n -DTC robust to perturbationsof the drive ( n is indicated in blue font for some of the resolvedDTCs). Especially remarkable are fractional n -DTCs, with n = p/q and p and q some coprime integers. This spectrum refers to the LMGlimit ( α = 0 , λ = 0 ), at fixed interaction J = 0 . , restricting to thefirst frequency Brillouin zone − . ≤ ν ≤ . , for and drive periods in the top and bottom panels, respectively. Alternatively, in clean (i.e., non-disordered) systems, heatingcan be escaped with all-to-all interactions [12, 25, 28], or sig-nificantly slowed down with long-range interactions [29–32].Very recently, Ref. [32] has provided the theoretical frame-work to study Floquet, clean, long-range interacting systems,in which novel prethermal phases of matter are expected.While their framework allows for the possibility of n -DTCswith n larger than the size of the local (or single-particle) a r X i v : . [ c ond - m a t . o t h e r] J a n Hilbert space, their concrete examples are limited to n = 2 .From our analysis below, we see that part of the difficulty innumerically observing what we call ‘higher-order’ DTCs maylie in their emergence at system sizes which are typically be-yond the reach of exact diagonalization.Here, we overcome this limitation by considering a sys-tem amenable to a set of complimentary methods, which en-able us to discover an unusually rich dynamical phase dia-gram hosting a zoo of novel, exotic, (and arguably prether-mal) non-equilibrium phases of matter. More specifically, weshow that a clean spin- / chain in the presence of long-range interactions (Fig. 1a) can sustain robust higher-order n -DTCs with integer and, remarkably, even fractional n > (e.g., n = 3 , , / and beyond). These novel dynamicalphases give rise to a peculiar fragmentation of the magnetiza-tion spectrum, which is intriguingly reminiscent of the plateaustructure of the Fractional Quantum Hall Effect.In the following we present a rather general model of long-range interacting spins, thoroughly study its semiclassical(that is, mean-field) limit, and finally show that the physicsobserved extends far beyond the fine-tuned limit. We notethat our work is distinct from traditional MBL DTCs, whichdo not have such a semiclassical limit. On a conceptuallevel our analysis is closer to that of equilibrium statisticalphysics, where for example the ferromagnetic phase in theIsing model is best understood in a mean-field descriptionwhich is exact in the limit of all-to-all interactions. In ourout-of-equilibrium and clean setting, the existence of a con-ceptually simple mean-field limit is particularly valuable, andhighlights profound differences between the clean DTCs con-sidered here and the pioneering works on MBL DTCs. Results — We consider a one-dimensional chain of N spinsin the thermodynamic limit ( N → ∞ ), driven according to thefollowing time-periodic Hamiltonian H ( t ) = J N N,α N (cid:88) i,j =1 i (cid:54) = j σ zi σ zj ( r i,j ) α + λ N (cid:88) j =1 σ zj σ zj +1 − πh [1 + sin(2 πt )] N (cid:88) j =1 σ xj , (1)where σ ( x,y,z ) j denote the standard Pauli operators for the j -thspin, periodic boundary conditions are assumed, and both (cid:126) and the drive frequency have been set to . J measures thestrength of a power-law interaction with characteristic expo-nent α , λ is the strength of a nearest-neighbor interaction, and πh is the average over one drive period of the monochromatictransverse magnetic field. The Kac normalization N N,α = (cid:80) Nj =2 1( r ,j ) α guarantees extensivity, and conveniently allowsto stretch the model to the Lipkin-Meshkov-Glick (LMG)limit of all-to-all interactions ( α = λ = 0 ), in which the un-derlying complex physics is reduced to its essence and mosteasily interpreted.The dynamics from an initially z -polarized state | ψ (0) (cid:105) = |↑ , ↑ , . . . , ↑(cid:105) is integrable in the non-interacting limit J = λ = 0 , for which the magnetization m ( t ) = (cid:104) σ zj (cid:105) ( t ) at strobo-scopic times t = 0 , , , . . . reads m ( t ) = cos(2 πht ) , that is h is the system’s characteristic frequency. The essential ques-tion to diagnose a n -DTC is whether, upon switching on theinteractions, there exists a finite range of h for which the sys-tem’s characteristic frequency ν remains instead locked to aconstant value /n < , signalling the stability of the DTC toperturbations of the drive. In the following we answer thisquestion affirmatively not only for the well-known n = 2 case, but, if the interactions are sufficiently long-range, alsofor integer and even fractional n > , corresponding to thehigher-order DTCs.For the sake of clarity, we first focus on the LMG limitof all-to-all interactions ( α = λ = 0 ), which allows for aconceptually simple semiclassical interpretation of the variousdynamical phases. The dynamics of the system is in this casedescribed by a semiclassical Gross-Pitaevskii equation (GPE)for the complex fields ψ ↑ and ψ ↓ (details in SupplementarySection I) dψ ↑ dt = πh [1 + sin 2 πt ] ψ ↓ − J | ψ ↑ | ψ ↑ ,dψ ↓ dt = πh [1 + sin 2 πt ] ψ ↑ − J | ψ ↓ | ψ ↓ , (2)where we can identify | ψ ↑ | − | ψ ↓ | → m = (cid:104) σ zj (cid:105) and ψ ∗↑ ψ ↓ = | ψ ↑ || ψ ↓ | e iθ → (cid:104) σ xj (cid:105) + i (cid:104) σ yj (cid:105) . While in the limit α → the dynamics described by the GPE (2) is indeed -dimensional and lacks any sense of locality, it carries thesignature of many-body interactions in its nonlinearity ratherthan in an exponentially large number of degrees of free-dom (similarly, e.g., to the paradigmatic mean-field equation m = tanh[ Jm/k B T ] of the Ising model in equilibrium).Note, the presence of such a limit highlights qualitative differ-ences between clean long-range DTCs and MBL DTCs, andis at the heart of their much richer phenomenology.The dynamics of the magnetization m is obtained integrat-ing the GPE (2) from an initially z -polarized state ( ψ ↑ (0) =1 , ψ ↓ (0) = 0 ), and the corresponding Fourier transform | ˜ m ( ν ) | versus the magnetic field strength h is plotted inFig. 1b. As it is well-known [12], the -DTC results in thesystem characteristic frequency ν being locked to / for h ≈ / . Surprisingly, the same locking occurs at frequencies /n with integer and fractional n > (e.g. n = 3 , , / ),giving rise to a fragmentation of the spectral line of ˜ m ( ν ) ina sequence of plateaus of constant frequency for a finite rangeof h ≈ /n . Each of these plateaus signals a higher-order(possibly fractional) DTC, the width of the plateau being asignature of the DTC’s robustness to drive perturbations. The‘halos’ surrounding the plateaus in Fig. 1 correspond to in-commensurate (non-subharmonic) frequencies adding a time-glassy aspect to the DTCs. The magnitude of these secondarypeaks is in the order of a few percent compared to the domi-nant subharmonic peak, resulting in weak aperiodic fine fea-tures on top of the subharmonic response.The plateau at ν = 0 for h ≈ signals the tendency of thespins to remain aligned along z in a dynamical ferromagnetic FIG. 2.
Phase space structure of the dynamical phases.
Poincar´e maps of the semiclassical dynamics (2) for various magnetic field strengths h and a fixed interaction J = 0 . . Red markers highlight the trajectory starting in the z -polarized state ( m = 1 , θ = 0 , green asterisk). (a)dynamical ferromagnet (F): the magnetization m remains ≈ at all times; (b) stroboscopic ferromagnet (sF): the magnetization m changessign during the micromotion and yet it remains positive at stroboscopic times; (c) -DTC: the system alternatively visits two islands of thephase space – one with m ≈ (numbered as ) at even times, and the other with m ≈ − (numbered as ) at odd times; (d,e) higher-order n -DTCs with integer n = 4 , , respectively: the system visits cyclically n islands of the phase space (accordingly numbered in red), with onetour of the islands corresponding to one complete revolution of the spins around the Bloch sphere. (f) Higher-order n -DTC with fractional n = q/p = 8 / : it takes p revolutions of the spins for the system to tour q islands of the phase space, resulting in a sharp magnetizationoscillation frequency ν = p/q . The insets on the right zoom on the island visited at times t = 8 k + 5 , k = 0 , , , . . . for the -DTC (top)and the / -DTC (bottom). phase (F). This corresponds to macroscopic quantum self-trapping of weakly driven bosons in a double well [33], whichcan in fact be exactly mapped to the LMG limit (details inSupplementary Section I). For h ≈ , , , . . . , the spins com-plete approximately , , , . . . revolutions around the Blochsphere at each drive period, respectively, and yet maintain apreferential alignement along z at stroboscopic times, in whatmay be called a stroboscopic-ferromagnetic phase (sF).Our results are confirmed by exact diagonalization studies.Thanks to the all-to-all coupling of the LMG limit, the dy-namics is in fact confined to the symmetric sector, whose sizegrows only linearly with the number of spins N . This allowsa scaling analysis extended up to large system sizes, show-ing a progressive emergence of the spectral line plateaus foran increasing number of spins N . For the standard -DTC,the plateau is clearly visible already for N (cid:39) , whereas,crucially, for the -DTC it appears only for N (cid:39) (see de-tails in Supplementary Section II). This observation stronglysuggests that signatures of the higher-order n -DTCs arisefor larger system sizes as compared to the standard -DTC,making them generally elusive to exact diagonalization tech-niques. This fact might explain the difficulties in observinghigher-order DTCs in the past and motivates the choice ofmodel (1) in the first place. The stroboscopic dynamics generated by the GPE (2) canbe conveniently described with Poincar´e maps, popular toolsin dynamical systems theory that here provide an immedi-ate interpretation of much of the underlying physics, whichto some extent characterizes the dynamical phases also whendeviating from the LMG limit. In Fig. 2, the trajectory start-ing in the z -polarized state (green asterisk) is highlighted withred markers. For a weak drive h ≈ , the spins tend to remainaligned along z in a dynamical ferromagnetic phase (a), givingrise to a Poincar´e map which closely resembles the phase por-trait of undriven bosons in a double well [33]. For h ≈ , themicromotion consists of approximately an entire revolutionof the spins around the Bloch sphere per period (not shown),with a preferential z alignment restored at stroboscopic timesdespite the detuning in the magnetic field strength (b). For h ≈ /n and n = 2 , , in (c), (d) and (e), respectively, the n -DTC results in the presence of n ‘islands’ in the phase spacewhich the system visits sequentially jumping from one to thenext at each drive period. In n drive periods, the system visitsall the n islands once, and the magnetization m completes oneoscillation.Furthermore, for h ≈ / the system behaves as a n -DTCwith fractional n = q/p = 8 / (f). In this case, the systemcyclically visits q islands of the phase space in q drive pe- m ( t ) | m ( (cid:1) ) | ∼ m ( t ) | m ( (cid:1) ) | ∼ Frequency (cid:1)
20 40Time t600 10080(a)(b) J = . J = . M agne t i c f i e l d s t r eng t h h - D T C
10 0.5Interaction J 10(d) Thermala bt = 4k t = 4k+3t = 4k+1 t = 4k+2 Trivial0 0.25 0.5|m(1/4)| ∼ (c) - D T C 〈 d 〉 t a bTrivialTrivial Thermal FIG. 3.
Many-body nature of the higher-order discrete time crys-tals.
The robustness of the higher-order time crystals is inducedby the interactions, justifing their classification as non-equilibriumphases of matter. For concreteness, we show this for the -DTC inthe LMG limit. (a,b) Magnetization m ( t ) at stroboscopic times (left)and respective Fourier transform | ˜ m ( ν ) | (right) for a slightly detunedmagnetic field strength h = 1 / . . For a weak interaction J =0 . in (a), the system trivially oscillates at frequencies ν ≈ ± . ,whereas a larger interaction J = 0 . in (b) re-establishes a robustsubharmonic response at frequency ν = 1 / . (c) Amplitude of thesubharmonic peak | ˜ m (1 / | in the ( J, h ) plane. The -DTC phaseopens up from the integrable point J = 0 , h = 0 . , that is the in-teraction makes the -DTC robust. (d) A thermal region of the phasespace is characterized by a finite value ∼ of the decorrelator time-average (cid:104) d (cid:105) t , corresponding to semiclassical chaos. Both (c) and(d) are computed over drive periods. riods. Differently from a q -DTC, however, during this timethe magnetization m completes p oscillations, resulting in acharacteristic frequency p/q . Finally, larger interactions J lead to chaotic Poincar´e maps (not shown, but analogous toRef. [12]), signaling thermalization [34].Since the island-to-island hopping that underpins the sub-harmonic response holds for any point of any island, the is-lands themselves can be interpreted as stability regions of theDTCs with respect to perturbations of the initial state. ThePoincar´e maps also provide an interpretation of the subdomi-nant frequencies visible as halos around the plateaus in Fig. 1,that are in fact associated with the revolution period of theintra-island orbits (e.g., those shown in the insets of Fig. 2).This also explains why these frequencies are sensitive to per-turbations of the drive, which deform the shape of the islandsand thus the orbits revolution periods, but are only weaklysensitive to perturbations of the initial state. It is well-established for the standard -DTC that the robustsubharmonic response hinges on the interaction being suffi-ciently strong. The fact that interactions are necessary for therobustness of DTCs is critical, as it underpins the many-bodynature of the DTCs and it justifies their classification as non-equilibrium phases of matter [10]. It becomes thus of primaryimportance to assess the role of the interactions also for thehigher-order DTCs. To this end, as a concrete example, inFig. 3 we investigate the effects of the interaction strength J on the -DTC. If the interaction is weak, a slightly mistakenmagnetic field strength h = 1 / (cid:15) , with (cid:15) (cid:28) , originatesin envelopes (that is, beatings) with period ∼ /(cid:15) , resultingin the Fourier transform ˜ m being peaked at ν ≈ h and intrivial dynamics (a). Crucially, stronger interactions can com-pensate the mistake in the flipping field (b): the envelopes in m ( t ) disappear, the peak in ˜ m is set back to the subharmonicfrequency ν = 1 /n , and the discrete time symmetry is bro-ken. The time-glassy character of the DTC is observed in asmall aperiodic modulation of the magnetization on top of thesubharmonic response.The subharmonic peak magnitude | ˜ m (1 / | can be used totrace out the -DTC phase in the ( J, h ) plane (c). The -DTCphase opens up from the integrable point J = 0 , h = 1 / for increasing interactions, in analogy with the opening of thestandard -DTC from J = 0 , h = 1 / [16]. This open-ing, which in dynamical system theory would be referred toas Arnold’s tongue, confirms that larger interactions J allowthe higher-order DTCs to bear larger detunings in the field h . However, at even larger J (cid:39) . semiclassical chaossets in and the time crystalline order is broken irrespectivelyof h . To see this, we introduce a decorrelator (cid:104) d ( t ) (cid:105) t (seeMethods), measuring the average distance between two ini-tially very close copies of the system evolving under Eq. (2). (cid:104) d (cid:105) t ∼ corresponds to sensitivity to the initial conditions,that is, to classical chaos, which in turn signals quantum ther-malization [34].As shown, the DTCs rely on the interactions being suffi-ciently (but not too) strong. Crucially, in contrast to the stan-dard -DTC, higher-order DTCs also necessitate the interac-tions to be sufficiently long-range. We now probe the robust-ness of the higher-order DTCs along yet a different directionin the drive space, exploring the effects of non-all-to-all in-teraction on higher-order DTCs, particularly assessing theirstability upon breaking the mean-field solvability of the dy-namics with power-law ( α > ) and nearest-neighbor ( λ > )interactions. In this case, the system is no longer described asa collective spin, and spin-wave excitations are rather gener-ated. To account for them, we adopt a spin-wave approxima-tion (see Methods), in which the central dynamical variable isthe density of spin-wave excitations (cid:15) ( t ) (cid:15) = 2 N N (cid:88) k (cid:54) =0 (cid:104) b † k b k (cid:105) , (3)where b † k and b k are bosonic creation and annihilation opera-tors for the spin-waves excitations with momentum k . n . n . i n t e r a c t i on (cid:1) S p i n w a v e den s i t y (cid:2) ( t ) Time t 1 1.2 1.410 -2 -1 〈 (cid:2) 〉 t (a)10 -2 -1 (cid:3) (cid:3) T h e r m a l -3 N(b) (cid:3) c N1.410 Power exponent (cid:3) (d) 〈 (cid:2) 〉 t -2 -1 a ThermalPre -thermal
LMG (c) |m(+1/4)|0 0.25 0.5 ∼ a Power exponent (cid:3) FIG. 4.
Stability and prethermalization with power-law andnearest-neighbor interactions . The higher-order and fractionalDTCs survive, most likely in a prethermal fashion, when deviatingfrom the LMG limit. For concreteness, we focus on the -DTCat h = 0 . and J = 0 . , and consider the effects of power-law( α > ) and nearest-neighbor ( λ > ) interactions. (a) If the interac-tions are sufficiently long-range (that is α is small enough, here for afixed λ = 0 . ), the density of spin-wave excitations (cid:15) remains smallthroughout several time decades. Conversely, shorter-range interac-tions lead to the proliferation of spin-wave excitations that makes thesystem quickly thermalize destroying any time crystalline order [15].(b) A sharp transition between these two regimes is highlighted bythe time-average (cid:104) (cid:15) (cid:105) t over periods versus α (at a fixed λ = 0 . ).The critical α c at which (cid:104) (cid:15) (cid:105) t crosses the threshold . (inset), growsand possibly saturates with the system size N , suggesting the stabil-ity of the -DTC in the thermodynamic limit N → ∞ . (c,d) Thestability of the -DTC for a whole region of the parameter space sur-rounding the LMG point α = λ = 0 is highlighted plotting the mag-nitude of the subharmonic peak | ˜ m (1 / | and the average spin-wavedensity (cid:104) (cid:15) (cid:105) t in the ( α, λ ) plane. In the LMG limit ( λ = α = 0 ), no spin-wave excitationis generated and (cid:15) = 0 at all times. When departing fromsuch a limit, two scenarios are possible (Fig. 4a): (i) (cid:15) rapidlyreaches a plateau (cid:47) . (up to some small fluctuations), forwhich we consider the spin-wave approximation consistent,or (ii) (cid:15) rapidly grows to values (cid:39) , for which the spin-waveapproximation breaks down. Although the method is not exactand may fail to capture the very long-time physics, it suggeststhat (i) and (ii) correspond to prethermalization and thermal-ization, respectively [15, 35, 36].We observe that the higher-order DTCs are stable (at leastin a prethermal fashion) for sufficiently long-range interac-tions (that is, sufficiently small λ and α ), whereas thermal-ization quickly sets in for shorter-range interactions (Fig. 4a).The transition between these two dynamical phases is sharpand can be located comparing the spin-wave density time av- erage (cid:104) (cid:15) (cid:105) t with a threshold . (Fig. 4b). The stability of the n -DTC in the presence of competing power-law and nearest-neighbor interactions can be investigated in the ( α, λ ) planeplotting the amplitude of the subharmonic peak | ˜ m (1 /n ) | in Fig. 4c and the time-averaged spin-wave density (cid:104) (cid:15) (cid:105) t inFig. 4d. The n -DTC is stable for a whole region of the pa-rameter space surrounding the LMG point ( α = λ = 0 ),that is, if the interactions are sufficiently long-range. More-over, the DTC is also robust to arbitrary perturbations to theinitial state, as long as the initial spin-wave density (cid:15) (0) issufficiently small, as we have checked by injecting a smallamount of spin-wave excitations to the initial state. Finally,we note that the time-glassy character of the DTCs is main-tained throughout their whole stability region. In particular,we find small peaks at incommensurate frequencies that do(do not) vary when slightly perturbing the drive (initial condi-tions). Discussion and conclusion — Higher-order DTCs in cleanlong-range interacting systems are qualitatively distinct fromDTCs of MBL Floquet systems [32]. Indeed, the higher-orderDTCs require the establishment of order along directions dif-ferent from ± z . For instance, in the -DTC the spins are ap-proximately aligned along − y and + y at times t = 1 , , . . . and t = 3 , , . . . , respectively. In an MBL system, a disor-dered magnetic field or short-range interaction along z wouldimmediately scramble the system when the spins are far fromthe z axis, precluding the possibility of higher-order DTCs.Thus, our work establishes that translationally-invariant sys-tems with long-range interactions can circumvent these limi-tations [14, 32].In the LMG limit of all-to-all interactions, the model (1)has a low-dimensional semiclassical limit, which links the n -DTCs to the multifrequency mode locking of some nonlin-ear discrete maps, which is ubiquitous in the natural sciences[37–40]. On the one hand, our work establishes a connec-tion between this class of DCTs and dynamical system the-ory. On the other hand, it provide evidence for the stabilityof the higher-order DTCs in a whole region of the parameterspace surrounding the LMG limit, in the presence of compet-ing, mean-field breaking, long- and short-range interactions,that is in a genuinely quantum setting with no semiclassicalcounterpart.The choice of a continuous Floquet drive with constant-in-time interactions and monochromatic transverse magneticfield, together with the translational invariance, makes model(1) a prime candidate for experimental implementation. Forinstance, bosons in a double well [33] could be used to realizea truly all-to-all interacting, that is the LMG, model. In thiscase, the field pulses would be simply implemented loweringthe barrier between the two wells to allow particle tunnelling,and the fact that no time modulation for the particle-particleinteraction is necessary should provide a major simplification.Power-law interactions with tunable alpha ≤ α ≤ can in-stead be realized in trapped-ion experiments [17, 41, 42].In conclusion, we have discovered higher-order DTCs witha period that is not limited from above by the size of the lo-cal (or single-particle) Hilbert space. The dynamical phasespace fragments to host many higher-order n -DTCs with in-teger and even fractional n , at least in a prethermal fashion.Future work should attempt to gain further analytical under-standing regarding the role of long-range interactions in sta-bilizing the different higher-order DTCs. Most importantly, itshould be assessed what are the allowed fractions q/p that re-sult in a q/p -DTC? Further study should assess in more detailthe role of the Kac normalization and of the dimensionality onthe fate of the dynamical phases of matter presented here. Acknowledgements
It is a pleasure to thank F. Carollo, E. I. R. Chiacchio, F. M. Gam-betta, J. P. Garrahan and A. Lazarides for useful discussions. A. P. ac-knowledges support from the Royal Society. A. N. holds a Univer-sity Research Fellowship from the Royal Society and acknowledgesadditional support from the Winton Programme for the Physics ofSustainability.
Methods
Decorrelator.
We quantify the sensitivity to the initial conditions ofEq. (2) with a decorrelator d ( t ) [28, 43] d ( t ) = (cid:2) | ψ ↑ ( t ) | − | ψ (cid:48)↑ ( t ) | (cid:3) + (cid:2) | ψ ↓ ( t ) | − | ψ (cid:48)↓ ( t ) | (cid:3) , (4)measuring the distance in time between two initially very closecopies of the system evolving under Eq. (2). Specifically, we con-sider the following two close initial conditions ψ ↑ (0) = 1 , ψ ↓ (0) = 0 , (5) ψ (cid:48)↑ (0) = cos(∆ m ) e − i ∆ θ , ψ (cid:48)↓ (0) = sin(∆ m ) e + i ∆ θ , (6)with ∆ m = ∆ θ = 10 − . The decorrelator time-average (cid:104) d (cid:105) t isthen given by (cid:104) d (cid:105) t = 1 T + 1 T (cid:88) t =0 d ( t ) , (7)where T is the total simulation time, e.g. T = 10 in Fig. 3 andFig. 4. Spin-wave approximation.
Here, we briefly summarize the ideabehind the spin-wave approximation, which is thoroughly explainedin Refs. [35, 36] and the supplementary material therein, to whichwe redirect the reader for further details. In a DTC evolving froman initially z -polarized state | ψ (0) (cid:105) = |↑ , ↑ , . . . , ↑(cid:105) , the spins aremostly aligned at all times. Imperfections in the alignment can bedescribed within a Holstein-Primakoff transformation in terms ofbosonic spin-wave quasiparticles b † k . Crucially, the collective spin (cid:126)S = N (cid:80) Nj =1 (cid:126)σ j rotates in time, so that the Holstein-Primakofftransformation has to be performed in a rotating frame R (cid:48) =( X, Y, Z ) such that the, say, Z axis tracks the orientation of thecollective spin at all times. This tracking is encoded in the condi-tion (cid:104) S X (cid:105) = (cid:104) S Y (cid:105) = 0 , from which the dynamics of the rotatingframe is obtained self-consistently. The Holstein-Primakoff transfor-mation from spin degrees of freedom to bosonic degrees of freedomis then performed as σ Xj → b j + b † j , σ Yj → − i ( b j − b † j ) , and thespin-waves degrees of freedom are obtained after a Fourier trans-form ˜ b k = √ N (cid:80) Nj =1 e − ikj b j . On top of this, an approximation ismade in that the Hamiltonian is expanded to lowest non-trivial orderin the density of spin-wave excitations (cid:15) = N (cid:80) Nk (cid:54) =0 (cid:104) b † k b k (cid:105) , whichshould remain (cid:28) for the approximation to be consistent. This pro-cedure results in a set of N ordinary differential equations similar to Eqs. (26) and (29) in the Supplemental Material of Ref. [35], de-scribing the rotation of the new reference frame and the dynamics ofthe spin waves at the various momenta.Finally, we remark the main differences between the implementationof the spin-wave approximation in our work and in Ref. [35]. First,Ref. [35] considers a constant-in-time Hamiltonian, whereas the pa-rameters of the Hamiltonian in the present work are time-dependent.As a consequence, the parameters in the system of ordinary differ-ential equations become time-dependent. Second, Ref. [35] con-siders a nearest-neighbor interaction on top of an all-to-all one,whereas we consider the more general case of nearest-neighbor in-teraction on top of a power-law one. The cos k that appears in theequations of Ref. [35] is therefore substituted by a more generic ˜ J k = (cid:80) Nj =1 J r ,j e − ir ,j k in ours, where J r i,j contains both thenearest-neighbor and the power-law interactions. Supplementary Information for ”Higher-order and fractional discrete time crystals in cleanlong-range interacting systems”
Andrea Pizzi, Johannes Knolle and Andreas NunnenkampThe Supplementary Information are devoted to technical details of the derivations and complimentary results and is organizedas follows. In Section I we show the mapping of the Lipkin-Meshkov-Glick (LMG) model for fully-connected spins to a modelof bosons in a double well, and exploit it to obtain the Gross-Pitaevskii equation (GPE) in the limit of inifinite number of spins N . In Section II we present results from exact diagonalization, showing how the discrete time crystals (DTCs) emerge for anincreasing number of spins N . We show that higher-order n -DTCs emerge at much larger N as compared to the standard -DTC,which explains the difficulty to numerically observe them. In Section III we give an overview of the spin-wave approximation weemploy, whereas in Section IV we show that a binary Floquet protocol also leads to results similar to the ones of the continuousdriving discussed in the main text. I) GROSS-PITAEVSKII EQUATION
In this section, we derive the semiclassical GPE of motion for the LMG model of N fully-connected spins ( λ = α = 0 ). Tothis end, we first map the spin system into a model of bosons in a double well. We then derive Heisenberg equations for thebosonic operators, and we treat them semiclassically replacing the bosonic operators with complex numbers. Map to bosons in a double well
The Schwinger’s oscillator model of angular momentum connects the algebra of angular momentum and the algebra of twobosonic modes [44]. Here, we show this connection explicitly for the collective spin of a system of N spin- / .Given p = ( p , p , . . . , p N ) a permutation of the indexes , , . . . , N , we say P p the permutation operator acting as P p | s , s , . . . , s N (cid:105) = | s p , s p , . . . s p N (cid:105) , (S1)where s i ∈ {↑ , ↓} . The permutation operators are used to build the symmetrization operator S , defined as S = 1 √ N ! (cid:88) p P p . (S2)We say | n ↑ , n ↓ (cid:105) the symmetrized state with n ↑ spins up and n ↓ spins down, that is | n ↑ , n ↓ (cid:105) = 1 (cid:112) n ↑ ! n ↓ ! S| ↑ , ↑ , . . . , ↑ (cid:124) (cid:123)(cid:122) (cid:125) n ↑ , ↓ , ↓ , . . . , ↓ (cid:124) (cid:123)(cid:122) (cid:125) n ↓ (cid:105) . (S3)Since N is fixed and n ↓ = N − n ↑ , the notation | n ↑ , n ↓ (cid:105) is actually slightly redundant. In the following, we neverthelessprefer to stick with this notation for the sake of clarity. The states | n ↑ , n ↓ (cid:105) form a basis for the Hilbert subspace of completelysymmetrized states. Given an operator O commuting with the symmetrization operator S , the action of O on this subspace istherefore fully-characterized by its action on the states | n ↑ , n ↓ (cid:105) .In particular, we consider the ’collective’ operators (cid:16)(cid:80) Nj =1 σ αj (cid:17) with α = x, y, z , which indeed all commute with S .We have N (cid:88) j =1 σ xj | n ↑ , n ↓ (cid:105) = 1 (cid:112) n ↑ ! n ↓ ! S N (cid:88) j =1 σ xj | ↑ , . . . , ↑ (cid:124) (cid:123)(cid:122) (cid:125) n ↑ , ↓ , . . . , ↓ (cid:124) (cid:123)(cid:122) (cid:125) n ↓ (cid:105) (S4) = 1 (cid:112) n ↑ ! n ↓ ! n ↑ (cid:88) j =1 S| ↑ , . . . , ↑ (cid:124) (cid:123)(cid:122) (cid:125) n ↑ − , ↓ , . . . , ↓ (cid:124) (cid:123)(cid:122) (cid:125) n ↓ +1 (cid:105) + 1 (cid:112) n ↑ ! n ↓ ! N (cid:88) j = n ↑ +1 S| ↑ , . . . , ↑ (cid:124) (cid:123)(cid:122) (cid:125) n ↑ +1 , ↓ , . . . , ↓ (cid:124) (cid:123)(cid:122) (cid:125) n ↓ − (cid:105) (S5) = n ↑ (cid:112) ( n ↑ − n ↓ + 1)! (cid:112) n ↑ ! n ↓ ! | n ↑ − , n ↓ + 1 (cid:105) + n ↓ (cid:112) ( n ↑ + 1)!( n ↓ − (cid:112) n ↑ ! n ↓ ! | n ↑ + 1 , n ↓ − (cid:105) (S6) = (cid:113) n ↑ ( n ↓ + 1) | n ↑ − , n ↓ + 1 (cid:105) + (cid:113) ( n ↑ + 1) n ↓ | n ↑ + 1 , n ↓ − (cid:105) , (S7)and, similarly, N (cid:88) j =1 σ yj | n ↑ , n ↓ (cid:105) = i (cid:112) n ↑ ! n ↓ ! n ↑ (cid:88) j =1 S| ↑ , . . . , ↑ (cid:124) (cid:123)(cid:122) (cid:125) n ↑ − , ↓ , . . . , ↓ (cid:124) (cid:123)(cid:122) (cid:125) n ↓ +1 (cid:105) − i (cid:112) n ↑ ! n ↓ ! N (cid:88) j = n ↑ +1 S| ↑ , . . . , ↑ (cid:124) (cid:123)(cid:122) (cid:125) n ↑ +1 , ↓ , . . . , ↓ (cid:124) (cid:123)(cid:122) (cid:125) n ↓ − (cid:105) (S8) = i (cid:113) n ↑ ( n ↓ + 1) | n ↑ − , n ↓ + 1 (cid:105) − i (cid:113) ( n ↑ + 1) n ↓ | n ↑ + 1 , n ↓ − (cid:105) (S9)and, finally, N (cid:88) j =1 σ zj | n ↑ , n ↓ (cid:105) = ( n ↑ − n ↓ ) | n ↑ , n ↓ (cid:105) . (S10)Introducing standard bosonic operators a ↑ , a ↓ , a †↑ , a †↓ for the two bosonic modes labeled by ↑ and ↓ , we can thus write N (cid:88) j =1 σ xj = a †↑ a ↓ + a †↓ a ↑ , N (cid:88) j =1 σ yj = − i (cid:16) a †↑ a ↓ − a †↓ a ↑ (cid:17) , N (cid:88) j =1 σ zj = n ↑ − n ↓ (S11)with n ↑ = a †↑ a ↑ and n ↓ = a †↓ a ↓ .The Hamiltonian (1) in the LMG limit ( α = λ = 0 ) reads H = JN N (cid:88) i,j σ zi σ zj − πh [1 + sin(2 πt )] N (cid:88) j =1 σ xj , (S12)and is thus rewritten in terms of the bosonic operators as H = JN ( n ↑ − n ↓ ) − πh [1 + sin(2 πt )]( a †↑ a ↓ + a †↓ a ↑ ) . (S13)We elaborate on the first term of Eq. (S13) using m = n ↑ − n ↓ N and noting that (cid:40) n ↑ N = m n ↓ N = − m ⇒ n ↑ n ↓ N = 1 − m , (S14)from which m = (cid:18) n ↑ − n ↓ N (cid:19) = n ↑ + n ↓ − n ↑ n ↓ N = n ↑ + n ↓ N + m − , (S15)and, isolating m , m = 2 n ↑ + n ↓ N − n ↑ ( n ↑ −
1) + n ↓ ( n ↓ − N − N . (S16)Setting U = 4 J and τ ( t ) = πh [1 + sin 2 πt ] , up to irrelevant additional constant terms, the Hamiltonian thus reads H ( t ) = − τ ( t )( a †↑ a ↓ + a †↓ a ↑ ) + U N (cid:0) n ↑ ( n ↑ −
1) + n ↓ ( n ↓ − (cid:1) , (S17)where n ↑ = a †↑ a ↑ and n ↓ = a †↓ a ↓ . That is, the LMG model in the symmetric subspace is mapped to a model for N bosons in adouble well (the two wells being labeled by ↑ and ↓ ). Equations of motion
The bosonic representation of Eq. (S17) is particularly convenient to obtain dynamical equations. The Heisenberg equationsfor the bosonic operators read ( (cid:126) = 1 ) da ↑ d ( it ) = [ H ( t ) , a ↑ ] = τ ( t ) a ↓ − UN n ↑ a ↑ ,da ↓ d ( it ) = [ H ( t ) , a ↓ ] = τ ( t ) a ↑ − UN n ↓ a ↓ . (S18)In the limit N → ∞ , upon replacing a ↑ → √ N ψ ↑ and a ↓ → √ N ψ ↓ , with ψ ↓ and ψ ↓ complex fields, we finally derive thefollowing Gross-Pitaevskii equation dψ ↑ d ( it ) = τ ( t ) ψ ↓ − U | ψ ↑ | ψ ↑ ,dψ ↓ d ( it ) = τ ( t ) ψ ↑ − U | ψ ↓ | ψ ↓ . (S19)For an operator ˆ O = f ( a ↑ , a ↓ , a †↑ , a †↓ ) written as a function f of the bosonic operators, the beyond-mean-field dynamics of theexpectation value O ( t ) = (cid:104) ˆ O (cid:105) ( t ) can be generally computed within a Truncated Wigner approximation (TWA) as O ( t ) ≈ (cid:104) f ( ψ ↑ ( t ) , ψ ↓ ( t ) , ψ ∗↑ ( t ) , ψ ∗↓ ( t )) (cid:105) ψ ↑ (0) ,ψ ↓ (0) (S20)where (cid:104) . . . (cid:105) ψ ↓ (0) ,ψ ↑ (0) denotes the average over an ensemble of stochastic semiclassical initial conditions ψ ↑ (0) and ψ ↓ (0) thatare drawn according to the quantum initial condition, and then evolve in time with the GPE (S19).In particular, let us consider as initial condition the symmetrized state | ψ (cid:48) (0) (cid:105) with magnetization m (cid:48) = 1 − δ , with < δ (cid:28) and m (cid:48) N integer (which, since we assume N → ∞ , does not restrict the possible values for m (cid:48) ) | ψ (cid:48) (0) (cid:105) = (cid:12)(cid:12) n (cid:48)↑ , n (cid:48)↓ (cid:11) , n (cid:48)↑ = m (cid:48) N, n (cid:48)↓ = (1 − m (cid:48) ) N, (S21)for which the TWA is performed considering the following ensemble of initial conditions ψ (cid:48)↑ (0) = (cid:114) − δ e iθ ↑ (0) , ψ (cid:48)↓ (0) = (cid:114) δ e iθ ↓ (0) , (S22)with θ ↑ (0) and θ ↓ (0) independent uniform random numbers between and π . Thanks to a gauge transformation, we can alwayschange the initial conditions (S22) into ψ (cid:48)↑ (0) = (cid:114) − δ , ψ (cid:48)↓ (0) = (cid:114) δ e iθ , (S23)where θ is also a random number between and π . Consider now the limit of δ → , that is of | ψ (cid:48) (0) (cid:105) → | ψ (0) (cid:105) = |↑ , ↑ , . . . , ↑(cid:105) . In this limit, the ensemble of stochastic initial conditions (S23) shrinks in the phase space of complex coordinates ψ (cid:48)↑ (0) and ψ (cid:48)↓ (0) towards the point ψ ↑ (0) = 1 − ψ ↓ (0) = 1 . If the GPE is nonchaotic, the points of the shrinked ensemble followclose trajectories, so that the TWA average in Eq. (S20) can be actually replaced by the evaluation of f for a single trajectory ofthe ensemble, say the one starting in ψ ↑ (0) = 1 − ψ ↓ (0) = 1 . In contrast, if the GPE is chaotic, because of sensitivity to initialconditions, the ensemble quickly spreads, scrambling across the classical phase space. In this case the ensemble trajectories atlong-times in (S20) interfere destructively, washing out any time oscillation of O : the TWA results in thermalization.In the limit | ψ (cid:48) (0) (cid:105) → | ψ (0) (cid:105) = |↑ , ↑ , . . . , ↑(cid:105) , it is therefore convenient to consider the following ’single-shot GPE’ [28], tobe run just once O ( t ) → (cid:0) f ( ψ ↑ ( t ) , ψ ↓ ( t ) , ψ ∗↑ ( t ) , ψ ∗↓ ( t )) (cid:1) ψ ↑ (0)=1 ,ψ ↓ (0)=0 , (S24)which is then expected to be accurate when nonchaotic, and to signal quantum thermalization when chaotic, which motivates theuse of the symbol “ → ” , rather than “ = ” .In particular, considering the observable (cid:126)S = N (cid:80) Nj =1 (cid:126)σ j , from Eq. (S25) we thus write (cid:104) S x (cid:105) + i (cid:104) S y (cid:105) → ψ ∗↑ ψ ↓ , (cid:104) S x (cid:105) → | ψ ↑ | − | ψ ↓ | . (S25) II) EXACT DIAGONALIZATION
The LMG model is also particularly suitable for exact diagonalization techniques. In fact, the size of the symmetric subspace( N + 1 ) scales only linearly with the system size N , making exact diagonalization viable for systems which are much larger thanthe ones typically considered in less-symmetric 1D models. In Fig. S1 we show the Fourier transform ˜ m ( ν ) of the magnetization m = (cid:104) σ zj (cid:105) for an initially z -polarized state | ψ (0) (cid:105) = | N ↑ (cid:105) = |↑ , ↑ , . . . , ↑(cid:105) . For increasing N , we can observe the emergenceof the constant-frequency plateaus signaling the n -DTCs. First, this confirms that the GPE captures the correct dynamics forlarge N , as expected. Second, we notice that, for the -DTC and the -DTC, the plateau is clearly visible only for N (cid:39) and N (cid:39) , respectively. This is remarkable, as it suggests that higher-order DTCs may be numerically harder to observe asthey could appear at larger system sizes, generally beyond the reach of exact diagonalization techniques. This observation mightexplain why higher-order DTCs have so far remained elusive to numerical examples [32]. F r equen cy | (cid:1) | Magnetic field strength h0.47 N = 3 N = 10 N = 50 N = 1000.520.500.48 0.530.5 0.47 0.530.5 0.47 0.530.5 0.47 0.530.50.250.270.250.23 0.290.27 0.25 0.290.27 0.25 0.290.27 0.25 0.290.27 - D T C - D T C (cid:1) (cid:2) (GPE) | m ( (cid:1) ) | ∼ FIG. S1.
Exact diagonalization and system size scaling . For an initially z -polarized state | ψ (0) (cid:105) = |↑ , ↑ , . . . , ↑(cid:105) we plot the Fouriertransform ˜ m ( ν ) of the magnetization m = (cid:104) σ zj (cid:105) for various number of spins N and a fixed interaction strength J = 0 . . Results from theGPE in the limit N → ∞ are reported on the right as a reference. (top row) For h ≈ . , the -DTC plateau at frequency . emerges forincreasing N , being clearly visible already for N (cid:39) . (bottom row) For h ≈ . , the -DTC plateau at frequency . also emerges forincreasing N , but it does so only for considerably larger system sizes N (cid:39) , which might explain the difficulties in observing higher-orderDTCs. III) SPIN-WAVE APPROXIMATION
In this section, we briefly summarize the idea behind the spin-wave approximation, which is thoroughly explained in Refs. [35,36] and the supplementary material therein, to which we redirect the reader for further details.In a DTC evolving from an initially z -polarized state | ψ (0) (cid:105) = |↑ , ↑ , . . . , ↑(cid:105) , the spins are mostly aligned at all times. Imperfec-tions in the alignment can be described within a Holstein-Primakoff transformation in terms of bosonic spin-wave quasiparticles b † k . Crucially, the collective spin (cid:126)S = N (cid:80) Nj =1 (cid:126)σ j rotates as a function of time in the lab frame. Therefore, the Holstein-Primakoff transformation has to be performed in a rotating frame ( X, Y, Z ) such that the, say, Z axis tracks the orientation ofthe collective spin at all times. This tracking is encoded in the condition (cid:104) S X (cid:105) = (cid:104) S Y (cid:105) = 0 , from which the dynamics of therotating frame is obtained self-consistently. On top of this, an approximation is made in that the Hamiltonian is expanded tolowest non-trivial order in the density of spin-wave excitations (cid:15) = N (cid:80) Nk (cid:54) =0 (cid:104) b † k b k (cid:105) , which should remain (cid:28) for the approxi-mation to be consistent. This procedure results in a set of N ordinary differential equations similar to Eqs. (26) and (29) in theSupplemental Material of Ref. [35], describing the rotation of the new reference frame and the dynamics of the spin waves atthe various momenta.Finally, we remark the main differences between the implementation of the spin-wave approximation in our work and inRef. [35]. First, Ref. [35] considers a constant-in-time Hamiltonian, whereas the parameters of the Hamiltonian in the presentwork are time-dependent. As a consequence, the parameters in the system of ordinary differential equations become time-dependent. Second, Ref. [35] considers a nearest-neighbor interaction on top of a fully-connected one, whereas we consider themore general case of nearest-neighbor interaction on top of a power-law one. The cos k that appears in the equations of Ref. [35]is therefore substituted by a more generic ˜ J k = (cid:80) Nj =1 J r ,j e − ir ,j k in ours, where J r i,j contains both the nearest-neighbor andthe power-law interactions. IV) BINARY DRIVING
In this section, we compliment the results from the main paper showing that the fragmentation of the phase diagram to host amultitude of higher-order DTCs also occur for a binary Floquet protocol. In particular, we consider a periodic Hamiltonian with
Magnetic field strength h F r equen cy | (cid:1) | F 2 s-F (cid:1) )| ∼ FIG. S2.
Phase diagram fragmentation for a binary Floquet driving . This figure is in complete analogy with Fig. 1b of the main paper, butconsiders the binary Floquet protocol (S26). We plot the Fourier transform ˜ m ( ν ) of the magnetization m = (cid:104) σ zj (cid:105) in the plane of the transversefield strength h and of the frequency ν , for a fixed interaction strength J = 0 . and computed over periods. We observe that the spectrallines fragment in plateaus with constant frequency, signaling the n -DTCs, (dynamic) ferromagnet and stroboscopic-ferromagnet. In blue, weindicate the index n of some of the resolved DTCs. period alternating fully-connected interaction and transverse field H ( t ) = +2 JN N (cid:88) i,j =1 σ zi σ zj < t ≤ . − πh N (cid:88) j =1 σ xj . < t ≤ . (S26)From the Hamiltonian (S26), in the limit N → ∞ , we can derive a GPE in complete analogy with Eq. (S19). Solving itfor the initially z -polarized state | ψ (0) (cid:105) = |↑ , ↑ , . . . , ↑(cid:105) , we obtain the spectrum ˜ m ( ν ) of Fig. S2. Also for this model, it ispossible to check within a spin-wave approximation that the dynamical phases persist when deviating from the LMG limit offully-connected spins, as long as the interaction range is large enough, in complete analogy with the results of Fig. 4 of the mainpaper. [1] M. Heyl, A. Polkovnikov, and S. Kehrein, Phys. Rev. Lett , 135704 (2013).[2] M. Heyl, EPL Europhys. Lett. , 26001 (2019).[3] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papi´c, Nat. Phys. , 745 (2018).[4] K. Sacha, Phys. Rev. A , 033617 (2015).[5] D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. , 090402 (2016).[6] C. W. von Keyserlingk and S. L. Sondhi, Phys. Rev. B , 245146 (2016).[7] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, Physical Review B , 085112 (2016).[8] V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. , 250401 (2016).[9] K. Sacha and J. Zakrzewski, Rep. Prog. Phys. , 016401 (2017).[10] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, arXiv preprint arXiv:1905.13232 (2019).[11] V. Khemani, R. Moessner, and S. Sondhi, arXiv preprint arXiv:1910.10745 (2019).[12] A. Russomanno, F. Iemini, M. Dalmonte, and R. Fazio, Phys. Rev. B , 214307 (2017).[13] R. Moessner and S. L. Sondhi, Nat. Phys. , 424 (2017).[14] D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. X , 011026 (2017).[15] B. Zhu, J. Marino, N. Y. Yao, M. D. Lukin, and E. A. Demler, arXiv preprint arXiv:1904.01026 (2019).[16] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Phys. Rev. Lett. , 030401 (2017).[17] J. Zhang, P. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, et al. , Nature , 217 (2017).[18] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, et al. , Nature , 221 (2017).[19] J. Rovny, R. L. Blum, and S. E. Barrett, Phys. Rev. Lett. , 180603 (2018).[20] J. Smits, L. Liao, H. T. C. Stoof, and P. van der Straten, Phys. Rev. Lett. , 185301 (2018).[21] F. Gambetta, F. Carollo, A. Lazarides, I. Lesanovsky, and J. Garrahan, arXiv preprint arXiv:1904.01026 (2019).[22] A. Lazarides, S. Roy, F. Piazza, and R. Moessner, arXiv preprint arXiv:1904.04820 (2019).[23] F. Gambetta, F. Carollo, M. Marcuzzi, J. Garrahan, and I. Lesanovsky, Physical review letters , 015701 (2019). [24] G. J. Sreejith, A. Lazarides, and R. Moessner, Phys. Rev. B , 045127 (2016).[25] F. M. Surace, A. Russomanno, M. Dalmonte, A. Silva, R. Fazio, and F. Iemini, arXiv preprint arXiv:1811.12426 (2018).[26] K. Giergiel, A. Kosior, P. Hannaford, and K. Sacha, Phys. Rev. A , 013613 (2018).[27] P. Matus and K. Sacha, Phys. Rev. A , 033626 (2019).[28] A. Pizzi, J. Knolle, and A. Nunnenkamp, Phys. Rev. Lett. , 150601 (2019).[29] A. Lerose, B. ˇZunkoviˇc, A. Silva, and A. Gambassi, Phys. Rev. B , 121112 (2019).[30] F. Liu, R. Lundgren, P. Titum, G. Pagano, J. Zhang, C. Monroe, and A. V. Gorshkov, Phys. Rev. Lett , 150601 (2019).[31] M. C. Tran, A. Ehrenberg, A. Y. Guo, P. Titum, D. A. Abanin, and A. V. Gorshkov, arXiv preprint arXiv:1908.02773 (2019).[32] F. Machado, D. V. Else, G. D. Kahanamoku-Meyer, C. Nayak, and N. Y. Yao, arXiv preprint arXiv:1908.07530 (2019).[33] M. Albiez, R. Gati, J. F¨olling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Phys. Rev. Lett , 010402 (2005).[34] J. G. Cosme and O. Fialko, Phys. Rev. A , 053602 (2014).[35] A. Lerose, J. Marino, B. ˇZunkoviˇc, A. Gambassi, and A. Silva, Phys. Rev. Lett , 130603 (2018).[36] A. Lerose, B. ˇZunkoviˇc, J. Marino, A. Gambassi, and A. Silva, Phys. Rev. B , 045128 (2019).[37] M. H. Jensen, P. Bak, and T. Bohr, Phys. Rev. Lett. , 1637 (1983).[38] V. Belykh, N. F. Pedersen, and O. Soerensen, Phys. Rev. B , 4860 (1977).[39] M. R. Guevara and L. Glass, J. Math. Biol. , 1 (1982).[40] D. M. Bramble and D. R. Carrier, Science , 251 (1983).[41] R. Islam, C. Senko, W. Campbell, S. Korenblit, J. Smith, A. Lee, E. Edwards, C.-C. Wang, J. Freericks, and C. Monroe, Science ,583 (2013).[42] J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk, and J. J. Bollinger, Nature , 489(2012).[43] T. Bilitewski, S. Bhattacharjee, and R. Moessner, Phys. Rev. Lett.121