Introduction

Due to its foundational and technological relevance, the study of condensed matter systems out of equilibrium has attracted growing interest in recent years, accounting among others for the discovery of dynamical phase transitions^1,2, quantum scars³, and, particularly, discrete time crystals (DTCs)^4–12. A DTC is a nonequilibrium phase of matter breaking the discrete time translational symmetry of a periodic (i.e., Floquet) drive. In the thermodynamic limit, the defining feature of an n-DTC is a subharmonic response at 1/nth of the drive frequency (n > 1), which is robust to perturbations and which persists up to infinite time¹². Following the first seminal proposals^4–8, DTCs have been widely investigated both theoretically and experimentally^12–23.

In this context, most work has focused on spin-1/2 systems, which have largely been shown to exhibit a 2-DTC where at every Floquet period each spin (approximately) oscillates between the states $$ ∣\uparrow ⟩$$ and $$ ∣\downarrow ⟩$$ leading to period doubling (n = 2). This fact naturally emerges from the dimension 2 of the local Hilbert space of the spins⁶, and can be generalized to n-DTCs in models of n-dimensional clocks^24,25. Another well-studied setting is that of bosons in a gravitational field bouncing on an oscillating mirror⁹, where the single-particle Hilbert space dimension is infinite (as the particle’s position is continuous) and where n-DTCs with arbitrary integer²⁶ and fractional²⁷n have been shown.

In these systems, heating to a featureless “infinite temperature” state is typically avoided by introducing disorder, which leads to a (Floquet) many-body-localized (MBL) phase^5,8. Alternatively, in clean (i.e., non-disordered) systems, heating can be escaped with all-to-all interactions^12,25,28, or significantly slowed down with long-range interactions^29–32. Very recently, ref. ³² has provided the theoretical framework 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-DTCs with n larger than the size of the local (or single-particle) Hilbert space, their concrete examples are limited to n = 2. From our analysis below, we see that part of the difficulty in numerically observing what we call “higher-order” DTCs may lie in their emergence at system sizes that are typically beyond the reach of exact diagonalization.

Here, we overcome this limitation by considering a system amenable to a set of complimentary methods, which enable us to discover an unusually rich dynamical phase diagram hosting a zoo of novel, exotic (and arguably prethermal) nonequilibrium phases of matter. More specifically, we show that a clean spin-1/2 chain in the presence of long-range interactions (Fig. 1a) can sustain robust higher-order n-DTCs with integer and, remarkably, even fractional n > 2 (e.g., n = 3, 4, 8/3, and beyond). These novel dynamical phases give rise to a peculiar fragmentation of the magnetization spectrum, which is intriguingly reminiscent of the plateau structure of the fractional quantum Hall effect.

Fig. 1

Higher-order and fractional discrete time crystals.

a A spin-1/2 chain with long-range interactions and initial z-polarization is driven with a monochromatic transverse magnetic field of strength h, inducing a spin precession around x. b The time crystallinity is probed by the Fourier transform $$ \mid \tilde{m}\left(\nu \right)\mid $$ of the magnetization along z. The spectrum fragments in a multitude of plateaus with constant frequency 1/n for a magnetic field strength h in a finite range ≈1/n, each of which signals a higher-order n-DTC robust to perturbations of the drive (n is indicated in blue font for some of the resolved DTCs). Especially remarkable are fractional n-DTCs, with n = p/q and p and q some coprime integers. This spectrum refers to the LMG limit (α = 0, λ = 0), at fixed interaction J = 0.5, restricting to the first frequency Brillouin zone −0.5 ≤ ν ≤ 0.5, for 500 and 2000 drive periods in the top and bottom panels, respectively.

In the following, we present a rather general model of long-range interacting spins, thoroughly study its semiclassical (i.e., mean-field) limit, and finally show that the physics observed extends far beyond the fine-tuned limit. We note that our work is distinct from traditional MBL DTCs, which do not have such a semiclassical limit. On a conceptual level, our analysis is closer to that of equilibrium statistical physics, where, for example, the ferromagnetic phase in the Ising model is best understood in a mean-field description, which is exact in the limit of all-to-all interactions. In our out-of-equilibrium and clean setting, the existence of a conceptually simple mean-field limit is particularly valuable, and highlights profound differences between the clean DTCs considered here and the pioneering works on MBL DTCs.

Results

We consider a one-dimensional chain of N spins in the thermodynamic limit (N → ∞), driven according to the following time-periodic Hamiltonian

The dynamics from an initially z-polarized state $$ ∣\psi \left(0\right)⟩=∣\uparrow ,\uparrow ,\dots ,\uparrow ⟩$$ is integrable in the noninteracting limit J = λ = 0, for which the magnetization $$ m\left(t\right)=⟨{\sigma }_j^z⟩\left(t\right)$$ at stroboscopic times t = 0, 1, 2, … reads $$ m\left(t\right)=\cos \left(2\pi ht\right)$$, that is h is the system’s characteristic frequency. The essential question to diagnose a n-DTC is whether, upon switching on the interactions, there exists a finite range of h for which the system’s characteristic frequency ν remains instead locked to a constant value 1/n < 1, signaling the stability of the DTC to perturbations of the drive. In the following, we answer this question affirmatively not only for the well-known n = 2 case, but, if the interactions are sufficiently long range, also for integer and even fractional n > 2, corresponding to the higher-order DTCs. Note, the $$ {\mathbb{Z}}_2$$ symmetry of the Hamiltonian in Eq. (1) raises, for a 2-DTC, the conceptual issue whether the subharmonic response stems from the time-symmetry breaking itself or it rather “piggybacks” on an underlying breaking of the $$ {\mathbb{Z}}_2$$ symmetry^7,14,33. This issue, however, disappears for the higher-order n-DTCs that, because the Hamiltonian lacks any $$ {\mathbb{Z}}_n$$ symmetry (n > 2), must indeed be a “genuine” manifestation of time-symmetry breaking.

For the sake of clarity, we first focus on the LMG limit of all-to-all interactions (α = λ = 0), which allows for a conceptually simple semiclassical interpretation of the various dynamical phases. The dynamics of the system is in this case described by a semiclassical Gross–Pitaevskii equation (GPE) for the complex fields ψ_↑ and ψ_↓ (details in Supplementary Note 1)

The dynamics of the magnetization m is obtained integrating the GPE (2) from an initially z-polarized state (ψ_↑(0) = 1, ψ_↓(0) = 0), and the corresponding Fourier transform $$ \mid \tilde{m}\left(\nu \right)\mid $$ versus the magnetic field strength h is plotted in Fig. 1b. As it is well known¹², the 2-DTC results in the system characteristic frequency ν being locked to 1/2 for h ≈ 1/2. Surprisingly, the same locking occurs at frequencies 1/n with integer and fractional n > 2 (e.g., n = 3, 4, 8/3), giving rise to a fragmentation of the spectral line of $$ \tilde{m}\left(\nu \right)$$ in a sequence of plateaus of constant frequency for a finite range of h ≈ 1/n. Each of these plateaus signals a higher-order (possibly fractional) DTC, the width of the plateau being a signature of the DTC’s robustness to drive perturbations. The “halos” surrounding the plateaus in Fig. 1 correspond to incommensurate (non-subharmonic) frequencies adding a time-glassy aspect to the DTCs. The magnitude of these secondary peaks is in the order of a few percent compared to the dominant subharmonic peak, resulting in weak aperiodic fine features on top of the subharmonic response.

The plateau at ν = 0 for h ≈ 0 signals the tendency of the spins to remain aligned along z in a dynamical ferromagnetic phase (F). This corresponds to macroscopic quantum self-trapping of weakly driven bosons in a double well³⁴, which can in fact be exactly mapped to the LMG limit (details in Supplementary Note 1). For h ≈ 1, 2, 3, …, the spins complete approximately 1, 2, 3, … revolutions around the Bloch sphere at each drive period, respectively, and yet maintain a preferential alignment along z at stroboscopic times, in what may be called a stroboscopic-ferromagnetic phase (sF).

Our results are confirmed by exact diagonalization studies. Owing to the all-to-all coupling of the LMG limit, the dynamics is in fact confined to the symmetric sector, whose size grows only linearly with the number of spins N. This allows a scaling analysis extended up to large system sizes, showing a progressive emergence of the spectral line plateaus for an increasing number of spins N. For the standard 2-DTC, the plateau is clearly visible already for N ⪆ 10, whereas, crucially, for the 4-DTC it appears only for N ⪆ 100 (see details in Supplementary Note 2). This observation strongly suggests that signatures of the higher-order n-DTCs arise for larger system sizes as compared to the standard 2-DTC, making them generally elusive to exact diagonalization techniques. This fact might explain the difficulties in observing higher-order DTCs in the past and motivates the choice of model (1) in the first place.

The stroboscopic dynamics generated by the GPE (2) can be conveniently described with Poincaré maps, popular tools in dynamical systems theory that here provide an immediate interpretation of much of the underlying physics, which to some extent characterizes the dynamical phases also when deviating from the LMG limit. In Fig. 2, the trajectory starting in the z-polarized state (green asterisk) is highlighted with red markers. For a weak drive h ≈ 0, the spins tend to remain aligned along z in a dynamical ferromagnetic phase (a), giving rise to a Poincaré map which closely resembles the phase portrait of undriven bosons in a double well³⁴. For h ≈ 1, the micromotion consists of approximately an entire revolution of the spins around the Bloch sphere per period, with a preferential z alignment restored at stroboscopic times despite the detuning in the magnetic field strength (b). For h ≈ 1/n and n = 2, 4, 8 in (c), (d), and (e), respectively, the n-DTC results in the presence of n “islands” in the phase space, which the system visits sequentially jumping from one to the next at each drive period. In n drive periods, the system visits all the n islands once, and the magnetization m completes one oscillation.

Fig. 2

Phase space structure of the dynamical phases.

Poincare maps of the semiclassical dynamics (2) for various magnetic field strengths h and a fixed interaction J = 0.5. Red markers highlight the trajectory starting in the z-polarized state (m = 1, θ = 0, green asterisk). a Dynamical ferromagnet (F): the magnetization m remains ≈1 at all times; b stroboscopic ferromagnet (sF): the magnetization m changes sign during the micromotion and yet it remains positive at stroboscopic times; c 2-DTC: the system alternatively visits two islands of the phase space—one with m ≈ 1 (numbered as 0) at even times, and the other with m ≈ − 1 (numbered as 1) at odd times; d, e higher-order n-DTCs with integer n = 4, 8, respectively: the system visits cyclically n islands of the phase space (accordingly numbered in red), with one tour 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/3: it takes p revolutions of the spins for the system to tour q islands of the phase space, resulting in a sharp magnetization oscillation frequency ν = p/q. The insets on the right zoom on the island visited at times t = 8k + 5, k = 0, 1, 2, … for the 8-DTC (top) and the 8/3-DTC (bottom).

Furthermore, for h ≈ 3/8 the system behaves as a n-DTC with fractional n = q/p = 8/3 (f). In this case, the system cyclically visits q islands of the phase space in q drive periods. Differently from a q-DTC, however, during this time the magnetization m completes p oscillations, resulting in a characteristic frequency p/q. Finally, for larger interactions J the Poincaré maps become chaotic (as in ref. ¹²), signaling thermalization³⁵.

Since the island-to-island hopping that underpins the subharmonic response holds for any point of any island, the islands themselves can be interpreted as stability regions of the DTCs with respect to perturbations of the initial state. The Poincaré maps also provide an interpretation of the subdominant frequencies visible as halos around the plateaus in Fig. 1, which are in fact associated with the revolution period of the intra-island orbits (e.g., those shown in the insets of Fig. 2). This also explains why these frequencies are sensitive to perturbations of the drive, which deform the shape of the islands and thus the orbits’ revolution periods, but are only weakly sensitive to perturbations of the initial state.

It is well established for the standard 2-DTC that the robust subharmonic response hinges on the interaction being sufficiently strong. The fact that interactions are necessary for the robustness of DTCs is critical, as it underpins the many-body nature of the DTCs and it justifies their classification as nonequilibrium phases of matter¹⁰. It becomes thus of primary importance to assess the role of the interactions also for the higher-order DTCs. To this end, as a concrete example, in Fig. 3 we investigate the effects of the interaction strength J on the 4-DTC. If the interaction is weak, a slightly mistaken magnetic field strength h = 1/4 + ϵ, with ϵ ≪ 1, originates in envelopes (i.e., beatings) with period ~1/ϵ, resulting in the Fourier transform $$ \tilde{m}$$ being peaked at ν ≈ h and in trivial dynamics (a). Crucially, stronger interactions can compensate the mistake in the flipping field (b): the envelopes in m(t) disappear, the peak in $$ \tilde{m}$$ is set back to the subharmonic frequency ν = 1/n, and the discrete time symmetry is broken. The time-glassy character of the DTC is observed in a small aperiodic modulation of the magnetization on top of the subharmonic response.

Fig. 3

Many-body nature of the higher-order discrete time crystals.

The robustness of the higher-order time crystals is induced by the interactions, justifying their classification as nonequilibrium phases of matter. For concreteness, we show this for the 4-DTC in the LMG limit. a, b Magnetization m(t) at stroboscopic times (left) and respective Fourier transform $$ \mid \tilde{m}\left(\nu \right)\mid $$ (right) for a slightly detuned magnetic field strength h = 1/4 + 0.01. For a weak interaction J = 0.1 in (a), the system trivially oscillates at frequencies ν ≈±0.26, whereas a larger interaction J = 0.5 in (b) re-establishes a robust subharmonic response at frequency ν = 1/4. c Amplitude of the subharmonic peak $$ \mid \tilde{m}\left(1/4\right)\mid $$ in the (J, h) plane. The 4-DTC phase opens up from the integrable point J = 0, h = 0.25, which is the interaction that makes the 4-DTC robust. d A thermal region of the phase space is characterized by a finite value ~1 of the decorrelator time average $$ {⟨d^2⟩}_t$$, corresponding to semiclassical chaos. Both (c) and (d) are computed over 10³ drive periods.

The subharmonic peak magnitude $$ \mid \tilde{m}\left(1/4\right)\mid $$ can be used to trace out the 4-DTC phase in the (J, h) plane (c). The 4-DTC phase opens up from the integrable point J = 0, h = 1/4 for increasing interactions, in analogy with the opening of the standard 2-DTC from J = 0, h = 1/2¹⁶. This opening, which in dynamical system theory would be referred to as Arnold’s tongue, confirms that larger interactions J allow the higher-order DTCs to bear larger detunings in the field h. However, at even larger J ⪆ 0.8 semiclassical chaos sets in and the time crystalline order is broken irrespectively of h. To see this, we introduce a decorrelator $$ {⟨d^2\left(t\right)⟩}_t$$ (see “Methods” and Supplementary Note 3), measuring the average distance between two initially very close copies of the system evolving under Eq. (2). $$ {⟨d^2⟩}_t~1$$ corresponds to sensitivity to the initial conditions, that is, to classical chaos, which in turn signals quantum thermalization³⁵.

As shown, the DTCs rely on the interactions being sufficiently (but not too) strong. Crucially, in contrast to the standard 2-DTC, higher-order DTCs also necessitate the interactions to be sufficiently long range. We now probe the robustness of the higher-order DTCs along yet a different direction in the drive space, exploring the effects of non-all-to-all interaction on higher-order DTCs, particularly assessing their stability upon breaking the mean-field solvability of the dynamics with power-law (α > 0) and nearest-neighbor (λ > 0) interactions. In this case, the system is no longer described as a collective spin, and spin-wave excitations are rather generated. To account for them, we adopt a spin-wave approximation (see “Methods”), in which the central dynamical variable is the density of spin-wave excitations ϵ(t)

In the LMG limit (λ = α = 0), no spin-wave excitation is generated and ϵ = 0 at all times. When departing from such a limit, two scenarios are possible (Fig. 4a): (i) ϵ rapidly reaches a plateau ⪅ 0.1 (up to some small fluctuations), for which we consider the spin-wave approximation consistent, or (ii) ϵ rapidly grows to values ⪆ 1, for which the spin-wave approximation breaks down. Although the method is not exact and may fail to capture the very long-time physics, it suggests that (i) and (ii) correspond to prethermalization and thermalization, respectively^15,36,37.

Fig. 4

Stability and prethermalization with power-law and nearest-neighbor interactions.

The higher-order and fractional DTCs survive, most likely in a prethermal fashion, when deviating from the LMG limit. For concreteness, we focus on the 4-DTC at h = 0.27 and J = 0.5, and consider the effects of power-law (α > 0) and nearest-neighbor (λ > 0) interactions. a If the interactions are sufficiently long-range (that is α is small enough, here for a fixed λ = 0.03), the density of spin-wave excitations ϵ remains small throughout several time decades. Conversely, shorter-range interactions lead to the proliferation of spin-wave excitations that makes the system quickly thermalize destroying any time crystalline order¹⁵. b A sharp transition between these two regimes is highlighted by the time average 〈ϵ〉_t over 10³ periods versus α (at a fixed λ = 0.03). The critical α_c at which 〈ϵ〉_t crosses the threshold 0.1 (inset), grows and possibly saturates with the system size N, suggesting the stability of the 4-DTC in the thermodynamic limit N → ∞. c, d The stability of the 4-DTC for a whole region of the parameter space surrounding the LMG point α = λ = 0 is highlighted plotting the magnitude of the subharmonic peak $$ \mid \tilde{m}\left(1/4\right)\mid $$ and the average spin-wave density 〈ϵ〉_t in the (α, λ) plane.

We observe that the higher-order DTCs are stable (at least in a prethermal fashion) for sufficiently long-range interactions (i.e., sufficiently small λ and α), whereas thermalization quickly sets in for shorter-range interactions (Fig. 4a). The transition between these two dynamical phases is sharp and can be located comparing the spin-wave density time average 〈ϵ〉_t with a threshold 0.1 (Fig. 4b). The stability of the n-DTC in the presence of competing power-law and nearest-neighbor interactions can be investigated in the (α, λ) plane plotting the amplitude of the subharmonic peak $$ \mid \tilde{m}\left(1/n\right)\mid $$ in Fig. 4c and the time-averaged spin-wave density 〈ϵ〉_t in Fig. 4d. The n-DTC is stable for a whole region of the parameter space surrounding the LMG point (α = λ = 0), that is, if the interactions are sufficiently long range. The DTC is also robust to arbitrary perturbations to the initial state, as long as the initial spin-wave density ϵ(0) is sufficiently small, as we have checked by injecting a small amount of spin-wave excitations to the initial state. For completeness, in the Supplementary Note 5 we also investigate the stability of the 4-DTC against the introduction of a longitudinal field of strength h_z. We found in this case that the expected prethermal nature of the DTC is only changed quantitatively, with the heating time scaling exponentially with the ratio of the driving frequency (for us always 1) and the energy scale h_z. Finally, we note that the time-glassy character of the DTCs is maintained 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 conditions).

Discussion

Higher-order DTCs in clean long-range interacting systems are qualitatively distinct from DTCs of MBL Floquet systems³². Indeed, the higher-order DTCs require the establishment of order along directions different from ±z. For instance, in the 4-DTC the spins are approximately aligned along −y and +y at times t = 1, 5, … and t = 3, 7, …, respectively. In an MBL system, a disordered magnetic field or short-range interaction along z would immediately scramble the system when the spins are far from the z-axis, precluding the possibility of higher-order DTCs. Thus, our work establishes that translationally invariant systems with long-range interactions can circumvent these limitations^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 nonlinear discrete maps, which is ubiquitous in the natural sciences^38–41. On the one hand, our work establishes a connection between this class of DCTs and dynamical system theory. On the other hand, it provides evidence for the stability of the higher-order DTCs in a whole region of the parameter space surrounding the LMG limit, in the presence of competing, mean-field breaking, long- and short-range interactions, which is in a genuinely quantum setting with no semiclassical counterpart.

The choice of a continuous Floquet drive with constant-in-time interactions and monochromatic transverse magnetic field, together with the translational invariance, makes model (1) a prime candidate for experimental implementation. For instance, bosons in a double well³⁴ could be used to realize a truly all-to-all interaction, which is the LMG, model. In this case, the field pulses would be simply implemented lowering the barrier between the two wells to allow particle tunneling, and the fact that no time modulation for the particle–particle interaction is necessary should provide a major simplification. Power-law interactions with tunable alpha 0 ≤ α ≤ 3 can instead be realized in trapped-ion experiments^17,42,43. For completeness, we note, anyway, that a phenomenology similar to that presented here also emerges in the case of a binary drive (see Supplementary Note 4 for details).

In conclusion, we have discovered higher-order DTCs with a period that is not limited from above by the size of the local (or single-particle) Hilbert space. The dynamical phase space fragments to host many higher-order n-DTCs with integer and even fractional n, at least in a prethermal fashion. Future work should attempt to gain further analytical understanding regarding the role of long-range interactions in stabilizing the different higher-order DTCs. Most importantly, it should be assessed what are the allowed fractions q/p that result in a q/p-DTC? Further study should assess in more detail the role of the Kac normalization and of the dimensionality on the fate of the dynamical phases of matter presented here.

Methods

Decorrelator

We quantify the sensitivity to the initial conditions of Eq. (2) with a decorrelator d²(t)^28,44

measuring the distance in time between two initially very close copies of the system evolving under Eq. (2). Specifically, we consider the following two close initial conditions

with Δ_m = Δ_θ = 10⁻⁶. The decorrelator time average $$ {⟨d^2⟩}_t$$ is then given by

where T is the total simulation time, for example, T = 10³ in Figs. 3 and 4.