Higher-order and fractional discrete time crystals in clean long-range interacting systems

Discrete time crystals are periodically driven systems characterized by a response with periodicity nT, with T the period of the drive and n > 1. Typically, n is an integer and bounded from above by the dimension of the local (or single particle) Hilbert space, the most prominent example being spin-1/2 systems with n restricted to 2. Here, we show that a clean spin-1/2 system in the presence of long-range interactions and transverse field can sustain a huge variety of different ‘higher-order’ discrete time crystals with integer and, surprisingly, even fractional n > 2. We characterize these (arguably prethermal) non-equilibrium phases of matter thoroughly using a combination of exact diagonalization, semiclassical methods, and spin-wave approximations, which enable us to establish their stability in the presence of competing long- and short-range interactions. Remarkably, these phases emerge in a model with continous driving and time-independent interactions, convenient for experimental implementations with ultracold atoms or trapped ions.


Supplementary Note 1: GROSS-PITAEVSKII EQUATION
In this Note, we derive the semiclassical GPE of motion for the LMG model of N fully-connected spins (λ = α = 0). To this end, we first map the spin system into a model of bosons in a double well. We then derive Heisenberg equations for the bosonic 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 two bosonic modes [? ]. Here, we show this connection explicitly for the collective spin of a system of N spin-1/2. Given p = (p 1 , p 2 , . . . , p N ) a permutation of the indexes 1, 2, . . . , N , we say P p the permutation operator acting as P p |s 1 , s 2 , . . . , s N = |s p1 , s p2 , . . . s p N , where s i ∈ {↑, ↓}. The permutation operators are used to build the symmetrization operator S, defined as We say |n ↑ , n ↓ the symmetrized state with n ↑ spins up and n ↓ spins down, that is Since N is fixed and n ↓ = N −n ↑ , the notation |n ↑ , n ↓ is actually slightly redundant. In the following, we nevertheless prefer to stick with this notation for the sake of clarity. The states |n ↑ , n ↓ form a basis for the Hilbert subspace of completely symmetrized states. Given an operator O commuting with the symmetrization operator S, the action of O on this subspace is therefore fully-characterized by its action on the states |n ↑ , n ↓ .
In particular, we consider the 'collective' operators N j=1 σ α j with α = x, y, z, which indeed all commute with S. We have = n ↑ (n ↓ + 1) |n ↑ − 1, n ↓ + 1 + (n ↑ + 1)n ↓ |n ↑ + 1, n ↓ − 1 , and, similarly, and, finally, Introducing standard bosonic operators a ↑ , a ↓ , a † ↑ , a † ↓ for the two bosonic modes labeled by ↑ and ↓, we can thus write with n ↑ = a † ↑ a ↑ and n ↓ = a † ↓ a ↓ . The Hamiltonian (1) in the LMG limit (α = λ = 0) reads and is thus rewritten in terms of the bosonic operators as We elaborate on the first term of Supplementary Eq. (13) using m = n ↑ −n ↓ N and noting that from which and, isolating m 2 , Setting U = 4J and τ (t) = πh[1 + sin 2πt], up to irrelevant additional constant terms, the Hamiltonian thus reads 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 a double well (the two wells being labeled by ↑ and ↓).

Equations of motion
The bosonic representation of Supplementary Eq. (17) is particularly convenient to obtain dynamical equations. The Heisenberg equations for the bosonic operators read ( = 1) In the limit N → ∞, upon replacing a ↑ → √ N ψ ↑ and a ↓ → √ N ψ ↓ , with ψ ↓ and ψ ↓ complex fields, we finally derive the following Gross-Pitaevskii equation For an operatorÔ = f (a ↑ , a ↓ , a † ↑ , a † ↓ ) written as a function f of the bosonic operators, the beyond-mean-field dynamics of the expectation value O(t) = Ô (t) can be generally computed within a Truncated Wigner approximation (TWA) as where . . . ψ ↓ (0),ψ ↑ (0) denotes the average over an ensemble of stochastic semiclassical initial conditions ψ ↑ (0) and ψ ↓ (0) that are drawn according to the quantum initial condition, and then evolve in time with the GPE (19).
In particular, let us consider as initial condition the symmetrized state |ψ (0) with magnetization m = 1 − δ, with 0 < δ 1 and m N integer (which, since we assume N → ∞, does not restrict the possible values for m ) for which the TWA is performed considering the following ensemble of initial conditions with θ ↑ (0) and θ ↓ (0) independent uniform random numbers between 0 and 2π. Thanks to a gauge transformation, we can always change the initial conditions (22) into where θ 0 is also a random number between 0 and 2π. Consider now the limit of δ → 0, that is of |ψ (0) → |ψ(0) = |↑, ↑, . . . , ↑ . In this limit, the ensemble of stochastic initial conditions (23) shrinks in the phase space of complex coordinates ψ ↑ (0) and ψ ↓ (0) towards the point ψ ↑ (0) = 1 − ψ ↓ (0) = 1. If the GPE is nonchaotic, the points of the shrinked ensemble follow close trajectories, so that the TWA average in Supplementary Eq. (20) can be actually replaced by the evaluation of f for a single trajectory of the ensemble, say the one starting in ψ ↑ (0) = 1 − ψ ↓ (0) = 1. In contrast, if the GPE is chaotic, because of sensitivity to initial conditions, the ensemble quickly spreads, scrambling across the classical phase space. In this case the ensemble trajectories at long-times in (20) interfere destructively, washing out any time oscillation of O: the TWA results in thermalization.
In the limit |ψ (0) → |ψ(0) = |↑, ↑, . . . , ↑ , it is therefore convenient to consider the following 'single-shot GPE' [? ], to be run just once which is then expected to be accurate when nonchaotic, and to signal quantum thermalization when chaotic, which motivates the use of the symbol " → ", rather than " = ". In particular, considering the observable S = 1 N N j=1 σ j , from Supplementary Eq. (11) we thus write

Supplementary Note 2: 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 than the ones typically considered in less-symmetric 1D models. In Fig. 1 we show the Fourier transform m(ν) of the magnetization m = σ z j for an initially z-polarized state |ψ(0) = |N ↑ = |↑, ↑, . . . , ↑ . For increasing N , we can observe the emergence of the constant-frequency plateaus signaling the n-DTCs. First, this confirms that the GPE captures the correct dynamics for large N , as expected. Second, we notice that, for the 2-DTC and the 4-DTC, the plateau is clearly visible only for N 10 and N 100, respectively. This is remarkable, as it suggests that higher-order DTCs may be numerically harder to observe as they could appear at larger system sizes, generally beyond the reach of exact diagonalization techniques. This observation might explain why higher-order DTCs have so far remained elusive to numerical examples [? ]. For h ≈ 0.25, the 4-DTC plateau at frequency 0.25 also emerges for increasing N , but it does so only for considerably larger system sizes N 100, which might explain the difficulties in observing higher-order DTCs.

Supplementary Note 3: SPIN-WAVE APPROXIMATION
In this Note, we briefly summarize the idea behind the spin-wave approximation, which is thoroughly explained in Refs. [? ? ] 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) = |↑, ↑, . . . , ↑ , the spins are mostly aligned at all times. Imperfections in the alignment can be described within a Holstein-Primakoff transformation in terms of bosonic spin-wave quasiparticles b † k . Crucially, the collective spin S = 1 N N j=1 σ 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 of the collective spin at all times. This tracking is encoded in the condition S X = S Y = 0, from which the dynamics of the rotating frame is obtained self-consistently. On top of this, an approximation is made in that the Hamiltonian is expanded to lowest non-trivial order in the density of spin-wave excitations = 2 N N k =0 b † k b k , which should remain 1 for the approximation to be consistent. This procedure results in a set of 2N ordinary differential equations similar to Eqs. (26) and (29) in the Supplemental Material of Ref. [? ], describing the rotation of the new reference frame and the dynamics of the spin waves at the various momenta.
Finally, we remark the main differences between the implementation of the spin-wave approximation in our work and in Ref. [? ]. First, Ref.
[? ] considers a constant-in-time Hamiltonian, whereas the parameters of the Hamiltonian in the present work are time-dependent. As a consequence, the parameters in the system of ordinary differential equations become time-dependent. Second, Ref.
[? ] considers a nearest-neighbor interaction on top of a fullyconnected one, whereas we consider the more general case of nearest-neighbor interaction on top of a power-law one. The cos k that appears in the equations of Ref.
[? ] is therefore substituted by a more genericJ k = N j=1 J r1,j e −ir1,j k in ours, where J ri,j contains both the nearest-neighbor and the power-law interactions.

Supplementary Note 4: BINARY DRIVING
In this Note, we compliment the results from the main paper showing that the fragmentation of the phase diagram to host a multitude of higher-order DTCs also occur for a binary Floquet protocol. In particular, we consider a  (26). We plot the Fourier transformm(ν) of the magnetization m = σ z j in the plane of the transverse field strength h and of the frequency ν, for a fixed interaction strength J = 0.5 and computed over 500 periods. We observe that the spectral lines fragment in plateaus with constant frequency, signaling the n-DTCs, (dynamic) ferromagnet and stroboscopic-ferromagnet. In blue, we indicate the index n of some of the resolved DTCs.
periodic Hamiltonian with period 1 alternating fully-connected interaction and transverse field From the Hamiltonian (26), in the limit N → ∞, we can derive a GPE in complete analogy with Supplementary Eq. (19). Solving it for the initially z-polarized state |ψ(0) = |↑, ↑, . . . , ↑ , we obtain the spectrumm(ν) of Fig. 2. Also for this model, it is possible to check within a spin-wave approximation that the dynamical phases persist when deviating from the LMG limit of fully-connected spins, as long as the interaction range is large enough, in complete analogy with the results of Fig. 4 of the main paper.

Supplementary Note 5: STABILITY AGAINST Z2 SYMMETRY-BREAKING TERMS
The Hamiltonian H(t) of Eq. 1 in the main text is Z 2 symmetric, since its commutator with a global spin flip N j=1 σ x j vanishes. For a 2-DTC, a question raises whether the subharmonic response truly stems from the breaking of the time-translational symmetry, or it rather "piggybacks" on an underlying breaking of the Z 2 symmetry[? ? ? ]. This question motivates an analysis to assess whether the 2-DTC is stable in the presence of a term breaking the underlying Z 2 symmetry explicitly, such as a longitudinal field. For higher-order and fractional n-DTCs, instead, no such an issue exist, because the Hamiltonian lacks any Z n symmetry (n > 2), and the DTCs must be a "genuine" manifestation of time-symmetry breaking alone, breaking an emergent Z n symmetry. Nonetheless, it is still interesting to investigate the effects that the explicit breaking of the Z 2 symmetry has on the n-DTCs. For completeness, in this Note, we therefore perform a stability analysis upon introducing a longitudinal field.
Consider the following Hamiltonian in which, on top of the Hamiltonian H(t) from Eq. (1) in the main text, we add a longitudinal field of strength πh z (the factor π comes from the analogy with the definitions for the transverse field) In Fig. 3, we investigate the stability of the 4-DTC in the presence of the longitudinal field h z breaking the Z 2 symmetry of the Hamiltonian. In the LMG limit (α = λ = 0), the introduction of a small field h z has qualitatively no significant effect. In the Poincaré map, the 4 islands characterizing the 4-DTC are in fact just slightly deformed for a small h z , suggesting that the subharmonic response lasts up to infinite time, similarly to h z = 0. Away from the LMG limit, instead, the dynamics within the spin-wave approximation shows a proliferation of the spin waves after a long prethermal regime. Although it is hard to assess the accuracy of this approximate method, the results suggest that the expected prethermal nature of the DTC beyond the LMG limit is only changed quantitatively, with the heating time scaling exponentially as t th ∼ e 1/hz . For h z = 0 we recover the results for which the 4-DTC appears infinitively long-lived within the spin-wave approximation.