# High-order time crystal phases and their fractal nature

Guido Giachetti<EMAIL_ADDRESS>SISSA, via Bonomea 265, I-34136 Trieste,
Italy  INFN, Sezione di Trieste, I-34151 Trieste, Italy  Andrea Solfanelli
SISSA, via Bonomea 265, I-34136 Trieste, Italy  INFN, Sezione di Trieste,
I-34151 Trieste, Italy  Lorenzo Correale  SISSA, via Bonomea 265, I-34136
Trieste, Italy  INFN, Sezione di Trieste, I-34151 Trieste, Italy  Nicolò
Defenu  Institut für Theoretische Physik, ETH Zürich, Wolfgang-Pauli-Str. 27
Zürich, Switzerland

###### Abstract

Discrete Floquet time crystals (DFTC) are characterized by the spontaneous
breaking of the discrete time-translational invariance characteristic of
Floquet driven systems. In analogy with equilibrium critical points, also
time-crystalline phases display critical behaviour of different order, i.e.
oscillations whose period is a multiple $p>2$ of the Floquet driving period.
Here, we introduce a new order parameter which is able to unambiguously detect
crystalline phases regardless of the value of $p$ and, at the same time, is a
useful tool for chaos diagnostic. This new paradigm allows us to investigate
the phase diagram of the long-range (LR) kicked Ising model to an
unprecedented depth, unveiling a rich landscape characterized by self-similar
fractal boundaries. Our solid theoretical picture describes the emergence of
DFTCs phase both as a function of the strength and period of the Floquet
drive, capturing the emergent $\mathbf{Z}_{p}$ symmetry in the Floquet-Bloch
waves.

_Introduction:_ The efficacy of technological application of quantum
mechanics, such as reliable quantum communication [1], high-precision quantum
metrology [2], and fault-tolerant quantum computation [3], depends on the
capability of preserving systems out-of-equilibrium, evading the detrimental
effects of thermalization, which naturally leads to the loss of locally stored
quantum information. Accordingly, much theoretical and experimental effort has
been devoted to the study of a wide range out-of-equilibrium phenomena [4, 5,
6, 7] including, among the others, thermalization of isolated quantum many-
body systems [8, 9, 10, 11], dynamical phase transitions [12, 13, 14, 15] and,
finally, the celebrated Discrete Floquet Time Crystals (DFTC) [16, 17, 18, 19,
20, 21, 22].

The latter are systems where the discrete time-translation symmetry, encoded
in the driven Hamiltonian $H(t)=H(t+T)$, is spontaneously broken at the level
of the expectation values, which display oscillations with a period being an
integer multiple of $T$. More precisely, in such a driven system a DFTC phase
exists if, taken a class of states $\ket{\psi}$ with short-ranged connected
correlations [17], it always exists an observable $\hat{O}$ such that the
time-evolved expectation value in the thermodynamic limit $N\to\infty$,

$\displaystyle O(t)=\lim_{N\to\infty}\langle\psi(t)|O|\psi(t)\rangle,$ (1)

satisfies the following conditions [23]:

  1. 1.

Time-translation symmetry breaking: $O(t+T)\neq O(t)$, although $H(t)=H(t+T)$.
This is equivalent to have long-range (LR) correlated Floquet eigenstates of
the propagator $U_{F}=U(t+T,t)$ [17].

  2. 2.

Rigidity: $O(t)$ must display periodic oscillations, with some period $\tau$,
in a finite and connected region of the Hamiltonian parameters space.

  3. 3.

Persistence: in the large system size limit $N\to\infty$, the oscillations of
$O(t)$ must persist for infinitely long time.

We refer to $O(t)$ as the order parameter of the time-crystal.

Conditions 1-3 can not be satisfied by a generic many-body quantum system, in
which the presence of an external driving would lead to the relaxation on an
infinite-temperature state, ruling out long-lived oscillations. Protecting
ordering against relaxation necessitate a mechanism to keep the impact of
dynamically generated excitations under control. This is achievable, for
example, in systems with strong disorder, where the diffusion of such
excitations is suppressed by many-body localization (MBL) [17, 24, 25, 19, 5,
21].

The possibility of generating a DFTC in clean systems has been studied in the
context of LR interacting models, i.e. models in which the interaction between
different lattice sites $\mathbf{i}$, $\mathbf{j}$ decay as a power law
$J_{\mathbf{i},\mathbf{j}}\sim|\mathbf{i}-\mathbf{j}|^{-\alpha}$. LR systems
have sparked a lot of attention recently, due to the possibility of
experimental realizations in atomic, molecular and optical (AMO) systems [26,
27, 28, 29, 30, 31, 32, 33] and to the fact that they exhibit plenty of unique
features both for the classical [34] and the quantum [35] regime. Indeed,
while LR interactions are known to alter the universal behaviour of critical
systems at equilibrium [36, 37], they also generate a wide range of
unprecedented out-of-equilibrium phenomena. These include novel form of
dynamical phase transitions [38, 39], defect formation [40, 41, 42, 38],
anomalous thermalization [43], information spreading [44, 45, 46] and
metastable phases [47, 48], whose features have no-counterpart in local
dynamics.

In the context of Floquet driven systems, LR interactions are known to enhance
the robustness of collective oscillations [49] and the presence of DFTC phases
have been actually established in $\alpha<d$ (strong LR) regime [23, 25, 50],
while for $\alpha>d$ the oscillations are not persistent in the
$N\rightarrow\infty$ limit [21, 22, 51, 52]. For a long time, the presence of
a DFTC phase of order $p$ (that is, with period $pT$, for integer values of
$p$), was thought to be connected with the underlying $\mathbb{Z}_{p}$
symmetry of the model [23, 25], as the Floquet driving can be engineered in
such way that each spin approximately oscillates between the
$\mathbb{Z}_{p}$-connected states. Nevertheless, high-order DFTCs were
recently observed also in systems with only $\mathbb{Z}_{2}$ symmetry, where
the order parameter oscillations display a period $pT$ (with $p>2$) given a
fixed driving period $T$ [53]. This phenomenon was deeper investigated in Ref.
[50], where it was pointed out how higher-order DFTC phases are associated to
an emergent $\mathbb{Z}_{p}$ symmetry which is not apparent in the original
model. Despite these preliminary studies, the comprehension of high-order time
crystal phases has been fundamentally hindered by the lack of a unique order
parameter capable of detecting all the DFTC phases simultaneously.

In this Letter we solve these issues, providing a fully fledged
characterization of the dynamical phase diagram for the paradigmatic case of
the LR Ising model, as a function both the amplitude and the period of the
driving. By introducing a novel order parameter, we are able to characterize
the DFTC phases of any order $p$, both as a function of the drive’s strength
and period, and to diagnose the transition to chaos, only relying on geometric
features of the dynamics. Such tight connection between the geometry of the
dynamics and the shape of our order parameter shall allow its use in a wide
range of applications.

We exploit the generality of our approach to uncover the hidden fractal nature
of the time-crystalline phase diagram, which is rich of self-similar
structures, each of them representing a different high-order DFTC phase. In
spite of this complexity, we are able to develop a theoretical framework which
captures quantitatively the salient features of the phase diagram, both with
analytical and numerical tools. Our analysis reproduces the known features of
the time-crystals and displays the emergence of new universal properties.
Furthermore, by performing extensive numerical simulations, we find that the
emergent $\mathbb{Z}_{p}$ symmetry is present also at the level of the Floquet
eigenstates, which can be interpreted as Bloch superposition of $p$ localized
semi-classical states.

_The model:_ As a paradigmatic example of the phenomenology we want to study,
let us consider the case of a chain of $N$ spin $1/2$ particles, interacting
through the Floquet-driven LR Ising Hamiltonian:

$H=-\frac{J}{2N_{\alpha}}\sum_{i>j}\frac{\hat{\sigma}^{i}_{x}\hat{\sigma}^{j}_{x}}{|i-j|^{\alpha}}+h(t)\sum_{i}\hat{\sigma}^{i}_{z}\
,$ (2)

where $\alpha<d$,
$\hat{\sigma}^{i}_{x}$,$\hat{\sigma}^{i}_{y}$,$\hat{\sigma}^{i}_{z}$ are the
Pauli operator relative to the lattice site $i$, $h(t)$ is a driving periodic
with period $T$ and $N_{\alpha}=\sum_{j\neq 0}|j|^{-\alpha}$ is the Kac
scaling factor needed in order to have an extensive energy [34]. At $t=0$, the
system is initialized in the ground state of the $h(t)=0$ Hamiltionian,
$\ket{\Psi_{0}}=\ket{\rightarrow\dots\rightarrow}$. Without loss of generality
we fix the energy scale such that $J=1$.

As shown in Refs. [47, 48] the out-of-equilibrium dynamics of the system for
any $\alpha<d$ is expected to converge to the $\alpha=0$ case in the
thermodynamic limit, so that in the following we are going to set $\alpha=0$
and $N_{\alpha}=N$. The dynamics under consideration is the one of the kicked
Ising model, in which the driving is given by

$h(t)=\psi\sum^{\infty}_{n=1}\delta(t-nT)\ .$ (3)

We choose to investigate the evolution of the observables
$m_{a}(t)=\frac{1}{N}\sum_{j}\left\langle\hat{\sigma}^{j}_{a}\right\rangle$
(where $a=x,y,z$) which correspond to the components of the magnetization of
the system. In the $N\to\infty$ limit, it is possible to write a closed
evolution equation for $\mathbf{m}(t)$ [54], which leads, in terms of the
stroboscopic magnetization $\mathbf{m}_{n}\equiv\mathbf{m}(nT)$, to the map

$\mathbf{m}_{n+1}=f(\mathbf{m}_{n})\equiv
R_{z}(2\psi)R_{x}(-m_{x,n}T)\mathbf{m}_{n}\ ,$ (4)

with $\mathbf{m}_{0}=(1,0,0)$ and where $R_{x,y,z}(\xi)$ is the rotation
matrix of an angle $\xi$ around the corresponding axis (see Supp. Mat. [55]
and Ref. [50]). The constraint $\mathbf{m}^{2}=1$ is preserved, as consequence
of the global angular momentum conservation, so that the motion takes place on
a Bloch sphere [54]. The $\mathbb{Z}_{2}$ symmetry of the model is instead
encoded into the dynamical symmetry $\psi\rightarrow\psi+\pi/2$,
$\mathbf{m}_{n}\rightarrow R_{z}(\pi n)\mathbf{m}_{n}$ of Eq. (4).

Figure 1: Left panel: Color plot of the order parameter $\zeta$ as a function
of the amplitude $\psi$ and the period $T$ of the driving, saturated at the
value $\zeta=\sqrt{2/3}$, with $n_{\rm max}=300$, $\delta\psi=1,6\cdot
10^{-3}$. Top-right panel: Normalized occurrence frequency $p(\zeta)$ of
$\zeta$. The gap between the DFTC phase, localized around $\zeta=0$, and the
quasi-periodic one, localized around $\zeta\approx 0.36$ is apparent. On the
right, the profile of the Gaussian distribution around $\zeta=\sqrt{2/3}$,
characteristic of the chaotic phase. Bottom-right panel: Detail of the $T<3$
region of the phase diagram. The order of the principal DFTC phases is
indicated, along with the theoretical prediction for small $T$ of the
boundaries between the phases (blue, red, green solid lines respectively for
the $p=1,2$, $p=4$, $p=3,6$ islands). The prolongation of the boundaries of
the $p=3$, $p=6$ islands (dashed red line) gives a good estimate of the onset
of caos, which disrupts the time-crystal phases at large $T$.

_Results:_ Depending on the values of $T$ and $\psi$ we observe different
dynamical phases. In particular we can individuate three different phases:
$(i)$ a phase in which $\mathbf{m}_{n}$ shows a chaotic, non-periodic
behavior, with a a strong sensitivity to initial conditions, $(ii)$ a quasi-
periodic phase, in which $\mathbf{m}_{n}$ has a non-integer period, sensitive
to small changes in $\psi$, $(iii)$ the periodic one in which $\mathbf{m}_{n}$
exhibits an integer period $p$, robust with respect to small changes of
$\psi$, with an additional small quasi-periodic modulation (referred as
micromotion in Ref. [53]). Within the latter scenario we can distinguish the
trivial $p=1$ case (in which the discrete time translation symmetry is not
broken) the period-doubling DFTC for $p=2$, and the higher-order DFTC phases
for $p>2$. Due to the dynamical $\mathbb{Z}_{2}$ symmetry of the system, the
presence of an order $p$ DFTC for some $\psi$ automatically implies the
presence of DFTC in correspondence of $\psi-\pi/2$, of order $p$ if $p$ is
even, $2p$ if $p$ is odd.

In order to distinguish between these phases, and in particular to pinpoint
the periodic phase, we need an order parameter which does not depend on the
order $p$ of the DFTC. To this extent we introduce the quantity

$\zeta^{2}=\frac{1}{n_{\rm max}}\sum^{n_{\rm
max}}_{n=0}\left(m_{x,n}(\psi+\delta\psi)-m_{x,n}(\psi)\right)^{2}\ ,$ (5)

in the limit $\delta\psi\rightarrow 0$, $n_{\rm max}\rightarrow\infty$ with
$n_{\rm max}\ \delta\psi=O(1)$ fixed. Intuitively this choice can be
understood by noticing that in the DFTC phase $\zeta$ only depends on the
fluctuations induced by the micro-motion and it is expected to be small and
possibly vanishing, while this is no longer true in the quasi-periodic region,
in which the period depends on $\psi$. As shown in the Supp. Mat. [55],
$\zeta$ actually exhibits a finite jump as we pass from the DFTC to the quasi-
periodic or chaotic phases. In the chaotic phase, on the other hand, the two
series $\mathbf{m}_{n}(\psi)$, $\mathbf{m}_{n}(\psi+\delta\psi)$ become
unrelated, so that $\zeta^{2}$ is expected to be distributed as a Gaussian
around $\zeta^{2}=2/3$ with a width of the order $O(n_{\rm max}^{-1/2})$ (see
Supp. Mat. [55]). This picture is indeed confirmed by the Fig. (1) (top-right
panel), in which the occurrence frequency of different values of $\zeta$ is
plotted, showing a gap between the DFTC phase ($\zeta\lesssim 0.25$) and the
quasi-periodic one ($\zeta\gtrsim 0.36$), with an additional Gaussian peak
around $\zeta=\sqrt{2/3}\approx 0.816$, corresponding to the chaotic phase.

We conclude that $\zeta$ is not only a natural choice to define the DFTC and
measures its distance from the perfect crystalline case, but it turns out to
be an useful choice in chaos diagnostic as well. This considerations are
expected to be very general, depending only on the geometrical features of the
map .

The value of the order parameter $\zeta$ as a function of $\psi$ and $T$ is
shown in Fig. 1 (left-panel). The resulting structure is strikingly complex
and convoluted. The symmetry around the $\psi=\pi/4$ is a consequence of the
dynamical $\mathbb{Z}_{2}$ symmetry. For small values of $T$ the quasi-
periodic phase is prevalent, while small islands of the periodic phase appear
around some particular values of $\psi$ which corresponds to rational
multiples of $\pi$. Initially, the size of these islands grows with $T$ and,
as they get closer to each other, chaotic regions start to appear around their
boundaries. Finally, all the islands corresponding to a DFTC of order $p>2$
disappear gradually as they are swallowed by the chaotic phase, the last one
being the $p=4$ one. In correspondence of particular values of the driving
period we have a revival of the higher-order DFTC phases, especially visible
in correspondence of $p=4$.

The boundary between the chaotic and the DFTC phase is not smooth: rather it
presents plenty of self-similar patterns which are repeated at smaller and
smaller scales. Taking as an example the boundary of the $p=3$ island, the
numerical estimate of its Minkowski - Bouligand dimension [56], $d_{MB}$,
gives the value

$d_{MB}\approx 1.4(2)>1$ (6)

(see Supp. Mat. [55]) indicating a fractal nature. This is compatible with the
estimate of $d_{MB}$ in the entire phase diagram displayed on Fig. 1. The
appearance of fractal scaling for time crystalline phase boundaries draws a
direct analogy with similar phenomena in traditional critical systems,
especially percolation, self-avoiding random walks and Potts model [57, 58],
where a rigorous connection between conformal invariance and stochastic
evolution has been verified [59, 60].

As already noticed in Ref. [50], the formation DFTC islands can be understood
within the formalism of area-preserving maps [61], and in particular can be
linked to the existence of Arnold tongues [62]. By exploiting this picture, we
were able to obtain analytic predictions about the structure of the phase
diagram at leading order in $T$. As shown in the Supp. Mat. [55] indeed, the
formation of a DFTC phase of order $p$ can be linked to the properties of the
$p$-iterated map $f^{p}(\mathbf{m})$, close to the resonant points
$\psi=\psi_{r}\equiv r\pi$ with $r=q/p$, where $q$,$p$ are coprime integers.
It is possible to derive an effective Hamiltonian which generates the
$p$-iterated evolution in the DFTC phase. While the couplings of the effective
Hamiltonian get renormalized, its general structure is independent on
$\psi_{r}$, accounting thus for the self-similar structure of the phase
diagram.

In turn, this universal structure allows to derive the form of the boundary
between the lower-order DFTC regions and the quasi-periodic ones; which are
consistent with the ones obtained through the numerics, as shown in Fig. 1
(bottom-right panel). In particular we find that the boundaries corresponding
to $r=0,1$ are well described for small $T$ by the equation $T=4|\delta\psi|$
(blue solid lines); the ones of the $r=1/4$ island by the curve
$T^{2}=32|\delta\psi|$ (green solid lines); while the equation
$T^{2}=24\sqrt{3}(\pm\delta\psi)$ describes the right/left boundary
respectively of the $r=1/3$, $r=1/6$ islands (red solid lines). According to
the criterion introduced by Chirikov [63], the value of $T$ at which two of
these curves intersect can be taken as an estimate of the threshold $T_{*}$
beyond which the chaos takes over: this gives the estimate
$T_{*}=(12\pi^{2})^{1/4}\approx 3.299$ for onset of chaos in the $p=4$ island,
which is in excellent agreement with the numerics. The uncovering of such a
universal, scale-invariant, structure in the DFTC phase is both unprecedented
and unexpected and allows for a deeper understanding of the phenomenon.

Figure 2: Color plot of the overlap
$|\braket{\Omega_{\theta,\phi}}{\eta_{m}}|^{2}$ with the spin coherent state
$\ket{\Omega_{\theta,\phi}}$ for different Floquet eigenstates
$\ket{\eta_{m}}$, corresponding to different phases, for $N=800$,
$\psi=\pi/2+0.01$, $T=1$ (top) and $T=10$ (bottom). While in the chaotic phase
(bottom panel) the eigenstate has no structure, the eigenstate $(b)$ (top
panel), which correspond to the DFTC phase with $p=4$, clearly exhibits the
structure of a Bloch wave-function localized around the $\mathbb{Z}_{4}$
symmetric wells. The eigenstate $(b)$ (top panel) has maximum overlap with the
spin coherent state corresponding to the initial conditions $\cos\theta=0$,
$\phi=\pi/2$. Initial conditions localized around the eigenstates $(a)$ and
$(c)$ (top panel) instead correspond to a quasi-periodic phase.

_Numerical analysis:_ Our comprehensive description of the DFTC phases is
grounded on the solid foundation of chaos theory. This profound parallelism is
granted by the presence of LR interactions in the Hamiltonian in (2), which
suppress fluctuations and thermalization effects in the thermodynamic limit.
On the other hand, it was estabilished in Ref. [23] that the $p=2$ DFTC
(around $\psi=\pi/2$) survives to finite sizes $N$, where the Floquet
eigenstates have equal amplitudes on both the two $\mathbb{Z}_{2}$ symmetric
sectors of the Hilbert space. One may wonder if such a picture can be
generalized to higher order DFTC phases, to infer the resilience of our
thermodynamic description to finite $N$. To do so, we performed the exact
diagonalization of the $\alpha=0$ model, explicitly computing the structure of
the Floquet eigenstates $\ket{\eta_{m}}$.

The fact that the modulus of the total spin $\mathbf{S}$ of the system is
conserved allows us to restrict ourselves to the subspace with
$\mathbf{S}^{2}=s(s+1)$, with $s=N/2$, and thus to consider larger sizes
($N=800$) [64, 23]. In Fig. 2 we show the projection of the Floquet
eigenstates $\ket{\eta_{m}}$ onto the spin coherent states
$\ket{\Omega_{\theta,\phi}}$ with $s=N/2$, $\theta$, $\phi$ being the polar
coordinates relative to $z$-axis [65]. These eigenstates correspond to the
three different phases of the system: while no structure is present in the
chaotic phase (Fig. 2, bottom panel), in the quasi-periodic phase the
eigenstate is localized in a connected region of the $(\theta,\phi)$ space
(Fig. 2,top panel: curves (a) and (c)), while in the $p=4$ DFTC phase it
appears localized around four, $\mathbb{Z}_{4}$ symmetric, points (Fig. 2,top
panel: curve (b)).

This behavior can be explained within our picture. Indeed, close to a
resonance, the Floquet evolution can be interpreted semi-classically as a
tight-binding hopping between adjacent wells in the space $(\phi,\cos\theta)$,
each one localized around $\psi=k\psi_{r}$ (with $k=1,\cdots p-1$), so that
the corresponding $\ket{\eta_{m}}$ will take the form of a Bloch superposition

$\braket{\Omega_{\theta,\phi}}{\eta_{m}}=\sum^{p-1}_{k=0}e^{2i\pi k/p}\
W_{m}(I,\phi-k\psi_{r})$ (7)

of the $p$ wavefunctions $W(I,\phi-k\phi_{r})$, each localized around
$\phi=k\psi_{r}$ and connected to the following by the Floquet propagator:
$U_{F}W_{m}(I,\phi-k\psi_{r})=e^{i\beta_{m}}W_{m}(I,\phi-(k+1)\psi_{r})$. As a
consequence, we expect this quantity to be localized around the fixed points
of the $p$-iterated evolution in the time-crystalline phase. The agreement
obtained between this semiclassical description and the numerical description
obtained by exact diagonalization proves the validity of our analysis also at
finite values of $N$ and $T$.

_Conclusions:_ In this work we have investigated the mechanism underlying the
formation of a higher-order Discrete Floquet Time-Crystal (DFTC) in clean
long-range (LR) systems by using the case of the LR kicked Ising model as a
paradigmatic example. The pressing need for this study was generated by the
recent depiction of high order DFTC phases [50, 53], which raised numerous
questions on the evolution and interplay of DFTC phases as a function of the
various drive parameters. Here, we presented a comprehensive theoretical
picture, substantiated by extensive numerical simulations, capable to answer
these questions. This was made possible by the introduction of a new, unique,
order parameter $\zeta$ which captures the features of the DFTC independently
of its order. This order parameter, moreover, turns out to be able to detect
quantum chaos and is a valuable tool in chaos diagnostic.

The resulting phase diagram of the model, as a function of the $\psi$, $T$
unveils a surprisingly rich and complex structure, which has not been noticed
before. It exhibits a self-similar structure and non integer fractal dimension
in the boundary of the chaotic phase, which we were able to numerically
estimate. DFTC phase islands develop as Arnold tongues [62] from the value
$T=0$, around the values $\psi_{r}=r\pi$, with $r\in\mathbb{Q}$, to be
swallowed by the onset of chaos as $T$ increases further. As discussed in Ref.
[50], these are intimately connected to the emergence of resonance in the
phase space, predicted by the Poincaré-Birkhoff theorem [66, 67] and its
quantum counterpart [68]. Through our theoretical analysis, however, we were
here able to draw quantitative predictions valid for the whole phase diagram,
and fully explain its emergent self-similar structure in terms of a an
effective Hamiltonian map with renormalized couplings.

Since our picture only relies on generic features of our system (such as the
LR nature of the interactions, the presence of an Hamiltonian structure, and
conservation laws) we expect it to apply to a wide class of models. This could
be verified explicitly in a large variety of systems, e.g. spin model with
continuous $O(n)$ symmetric LR interactions. Further investigation is needed
to check how our picture is modified in the prethermal regime $\alpha>d$.

Our newly introduced observable $\zeta$ is easily accessible in the
experiments, so that our predictions on the phase diagram may be tested e.g.
in NMR experiments on driven, ordered systems [22]. Our characterization of
the higher-order, stable DFTC phases could as well pave the way for a series
of advancement in the field of quantum technologies.

## Supplementary Material: Fractal nature of clean higher-order discrete time
crystals

###  .1 Derivation of Eq. (4)

Here we revise the derivation of the dynamic map of Eq. (4) of the main text,
i.e. the evolution equation of the magnetization
$m_{a}(t)=\sum_{j}\braket{\hat{\sigma}_{a}^{j}}/N$, at stroboscopic times
$t_{n}=nT$, in the thermodynamic limit, $N\rightarrow\infty$.

First we notice how, due to the impulsive nature of the magnetic field $h(t)$,
the Floquet propagator can be written as the product of two different
operators:

$U_{F}=e^{-2i\psi\hat{S}_{z}}e^{iJT\hat{S}_{x}^{2}/N}\ ,$ (8)

where we introduced the global spin operators

$\hat{S}_{a}=\frac{1}{2}\sum_{j}\hat{\sigma}_{a}^{j}\ ,$ (9)

with $a=x,y,z$. Being $\hat{S}_{x}$ the generator of $O(3)$ rotations around
the $x$ axis, the kick term $\exp(-2i\psi\hat{S}_{z})$ in Eq.(8) acts on the
observable $\mathbf{m}$ simply as a rotation around the $z$-axis. On the other
hand, the other term describes the evolution over one period $T$ of
$\mathbf{m}$ induced by the second term on the r.h.s of Eq. (8). The
Heisenberg equations of motion corresponding to such evolution for the
operators $\hat{S}_{a}$ are:

$\left\\{\begin{split}\frac{d}{dt}\hat{S}_{x}&=0\ ,\\\
\frac{d}{dt}\hat{S}_{y}&=\frac{J}{N}\left(\hat{S}_{x}\hat{S}_{z}+\hat{S}_{z}\hat{S}_{x}\right)\
,\\\
\frac{d}{dt}\hat{S}_{z}&=-\frac{J}{N}\left(\hat{S}_{x}\hat{S}_{y}+\hat{S}_{y}\hat{S}_{x}\right)\
.\end{split}\right.$ (10)

According to the general theory developed in Ref. [54], for $N\to\infty$ the
spin-spin correlation become negligible, that is
$\braket{\hat{S}_{a}\hat{S}_{b}}\simeq\braket{\hat{S}_{a}}\braket{\hat{S}_{b}}$,
so that Eqs. (10) become a closed set of equations for $m_{a}=2\left\langle
S_{a}\right\rangle/N$, namely:

$\left\\{\begin{split}\dot{m}_{x}&=0\ ,\\\ \dot{m}_{y}&=Jm_{x}m_{z}\ ,\\\
\dot{m}_{z}&=-Jm_{x}m_{y}\ .\end{split}\right.$ (11)

In turn this results after a time $T$ in a (clockwise) rotation around the
$x$-axis, of angle $JTm_{x}(t)$.
Posing $J=1$, the overall effect of the Eq. (8) on our observable $\mathbf{m}$
is the one of Eq.(4) of the main text.

###  .2 Hamiltonian chaos and the order parameter $\zeta$

In this Section we will derive in detail some properties of the order
parameter $\zeta$, in the context of the theory of symplectic maps.

Indeed, let us notice that the map of Eq. (4) inherits a Hamiltonian
structure, which is manifest in the fact that the area of any region on the
sphere $\mathbf{m}^{2}=1$ is preserved by the action of $f(\mathbf{m})$. In
terms of the usual polar coordinates along the $z$ axis,
$\mathbf{m}=(\sin\theta\cos\phi,\cos\theta\cos\phi,\cos\theta)$, the area
element can be written as $dS=d\cos\theta d\phi$, so that $\phi$ and
$I=\cos\theta$ are natural canonical conjugate variables for our system. This
is coherent with the picture developed in Ref.[54], and can be intuitively
understood by thinking of $I$ as the $z$ component of the angular momentum,
and $\phi$ the coordinate corresponding to rotation around the $z$ axis.

In particular, in the limit $T=0$, the map becomes

$\begin{split}I_{n+1}&=I_{n}\ ,\\\ \phi_{n+1}&=\phi_{n}+2\psi\ ,\end{split}$
(12)

with $I_{0}=0$, $\phi_{0}=\pi/2$. This correspond to the stroboscopic section
of an integrable dynamics, $\phi$ and $I$ playing the role of an angle-action
pair. In terms of the Floquet phases introduced in our paper, this implies a
quasi-periodic evolution of $\mathbf{m}_{n}$, with period $\pi/\psi$.

As a small $T$ is switched on, it can be treated as a perturbation to the map
(12). In this case, the fate of the system is described by the Kolmogorov-
Arnold-Moser theorem [69, 70, 71], according to which the torus $I=\rm const$
is only deformed as long as the corresponding frequency is not resonant. In
turns this means that the quasi-periodic phase survives as long as $\psi$ is
not close to a rational multiple of $\pi$. In the case of $\psi$ close to a
$q:p$ resonance, i.e. $\psi=\psi_{r}\equiv r\pi$, with $r=q/p$ and $p$ and $q$
coprime integers, according to the Poincare-Birkhoff theorem, pairs of
elliptic and unstable fixed points are expected to arise. In this case the
action of the $p$-iterated map $f^{p}(\mathbf{m})$ individuates different
regions in the phase space $(I,\phi)$, corresponding to different possible
behaviors of $\mathbf{m}_{n}$. In particular, if $(I_{0},\phi_{0})$ is far
from the fixed points, we have a rotation dynamics, with $\phi$ growing from
$0$ to $2\pi$, and we have a quasi-periodic behavior. If $(I_{0},\phi_{0})$ is
close to one of the centers instead, we have a libration dynamics, with $\phi$
oscillating around a finite value. As a consequence, $m_{n+p}$ remains close
to $m_{n}$, and we have a DFTC phase. At the boundaries between this two
regions a chaotic region is expected to arise which, as $T$ increases, grows
and possibly swallows the regular ones.

Let us analyze the consequence of this picture on the order parameter $\zeta$:
to consider different values of the amplitude, $\psi$, $\psi+\delta\psi$ it
means that we are considering to nearby initial conditions on the phase space.
Both in the DFTC phase and in the quasi-periodic one the evolution is not
chaotic, so that the two trajectories will diverge polynomially in time. As a
consequence

$\zeta^{2}=\frac{1}{n_{\rm max}}\sum^{n_{\rm
max}}_{n=0}\left(m_{x,n}(\psi+\delta\psi)-m_{x,n}(\psi)\right)^{2}\sim\frac{\ell}{n_{\rm
max}}\sum^{n_{\rm max}}_{n=0}\delta\psi^{2}n^{2}\sim\ell(\delta\psi
n_{\max})^{2}\ ,$ (13)

where $\ell$ depends on the the average distance between two randomly chosen
points of the two nearby trajectories. While the r.h.s. of Eq. (13) remains
finite and $O(1)$ as $n_{\rm max}\rightarrow\infty$, $\delta\psi\rightarrow
0$, $\ell$ jumps discontinuously as we pass from the libration regime
(corresponding to a DFTC phase) and the rotation one (corresponding to a
quasi-periodic phase). In particular, as we approach the the fixed point of
the iterate map, i.e. in the regime in which the micro-motion is becomes
negligible, $\zeta\rightarrow 0$, so that $\zeta$ is able to quantify how far
away is our system from the pure crystalline regime. Let us notice, however,
how the other values of $\zeta$ in this two phases are not universal, as they
depend on the choice of $\lim_{\delta\psi\rightarrow 0}n_{\rm max}\
\delta\psi$.

In the chaotic phase, instead, the trajectories diverge exponentially, so that
after a time-scale $n_{\rm max}\sim-\log(\delta\psi)$ the memory of the
initial condition is lost. In this case we can assume each $m_{x,n}(\psi)$ and
$m_{x,n}(\psi+\delta\psi)$ to be drawn from a set of equally distributed
random variables with zero mean. As a consequence, according to the central
limit theorem, $\zeta^{2}$ is distributed as a Gaussian around the value

$\left\langle\zeta^{2}\right\rangle=2\left\langle m_{x}^{2}\right\rangle$ (14)

with a variance $O(n_{\rm max}^{-1})$. Assuming furthermore the distribution
of the three components to be isotropic, and taking into account the
constraint $\mathbf{m}^{2}=1$, we have

$\left\langle
m_{x}^{2}\right\rangle=\frac{1}{3}\left\langle\mathbf{m^{2}}\right\rangle=\frac{1}{3}\
,$ (15)

so that $\left\langle\zeta^{2}\right\rangle=2/3$.

###  .3 Non-integer dimension of the chaotic phase boundary

Figure 3: Color plot of the order parameter $\zeta$ saturated at the value
$\zeta=\sqrt{2/3}$ with $n_{\rm max}=300$, $\delta\psi=10^{-4}$ (left) and
corresponding normalized occurrence frequency $p(\zeta)$ of $\zeta$ for this
region of the phase diagram.

The Minkowski-Bouligand, or box-counting, dimension is defined as follows: let
us cover our space with and evenly spaced square grid of side $\epsilon$. Said
$N(\epsilon)$ the number of boxes which lies on the boundary, then the
dimension is defined as

$d_{MB}\equiv-\lim_{\epsilon\rightarrow 0}\frac{\ln N(\epsilon)}{\ln\epsilon}\
.$ (16)

Since we are interested to the border between the chaotic and the DFTC phases,
we restricted ourselves to a region of the phase diagram in which only these
phases are present. Here we choose $\psi\in(0.94,1.01)$, $T\in(2.3,2.8)$,
which corresponds to the edge of the $p=3$ island. The region is shown in Fig.
3, left panel. To precisely define our boundary, we have here to define a
threshold $\zeta_{*}$ such that a point with $\zeta>\zeta_{*}$ is considered
to belong to the chaotic phase.

Since for any finite $n_{\rm max}$ the distribution of $\zeta$ in the chaotic
phase around $\zeta=0.816$ has a finite width, it is convenient, in this case,
to choose $\delta\psi$ such that, $n_{\rm max}\ \delta\psi<<1$. In this
regime, indeed, the distribution of $\zeta$ in the DFTC phase is sharply
peaked around $0$, and the separation of the two phases more pronounced. If,
however, we choose $\delta\psi$ to be too small, the diagnostic of the chaotic
region will not be accurate, as $\mathbf{m}_{n}(\psi)$ will not forget its
initial condition for $n\sim n_{\rm max}$. We checked that the choice $n_{\rm
max}=300$, $\delta\psi=10^{-4}$ is close to the optimal one. The separation
between the different phases, shown in Fig. 3, left panel, is clear.

The behavior of $d_{MB}$ as a function of the cutoff $\zeta^{*}$ Fig. 4, left
panel: as expected i, $d_{MB}$ shows a very weak dependence on $\zeta^{*}$
within a finite interval, signaling that the two phases are well distinct. The
corresponding value of the fractal dimension turns out to be $d_{MB}\approx
1.42$. Repeating the procedure with slightly different values of $n_{\rm max}$
and $\delta\psi$ gives an uncertainty of the order of order $10^{-2}$ on the
above result.

Figure 4: Behavior of the Minkowski-Bouligan dimension $d_{MB}$ of the
boundary between the chaotic DFTC phases as a function of the threshold
$\zeta_{*}$ chosen, for the case of the restricted area of Fig. (left) and for
the case of the whole phase diagram (right).

In Fig. 4 right panel, we reported the behavior of $d_{MB}$ as a function of
the threshold $\zeta_{*}$ for the whole phase-diagram of Fig. 1 of the main
text. To rule out the effect of the quasi-periodic phase, which is now
present, we have to restrict ourselves to the region $\zeta_{*}\gtrsim 0.36$.
As a consequence, the resulting estimate is far less accurate, and the
dependence of $d_{MB}$ on $\zeta_{*}$ more pronounced. We see, however, that
the our previous, local, estimate is fully compatible with this global
results.

###  .4 Small $T$ expansion

Let us now derive the results presented in the paper, which are valid for
small $T$ and $\psi$ close to $\psi_{r}=r\pi$, where $r=q/p$, $q$,$p$ being
coprime integers.

First, let us expand the map of Eq. (4) of the main text to the first order in
$T$. As in Section B, it is convenient to parametrize the magnetization as
$\mathbf{m}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$ and introduce
the canonical coordinates $I=\cos\theta$, $\phi$. We obtain

$\begin{split}I_{n+1}&=I_{n}-\frac{T}{2}(1-I^{2}_{n+1})\sin
2\phi_{n}+O(T^{2})\ ,\\\
\phi_{n+1}&=\phi_{n}+2\psi+\frac{T}{2}I_{n+1}\left(1+\cos
2\phi_{n}\right)+O(T^{2})\ .\end{split}$ (17)

Let us notice that the $O(T)$ terms on the r.h.s. of the above equations,
$I_{n+1}$ can be replaced by $I_{n}$ up to higher order corrections. This
choice guarantees that our approximation preserves the symplectic structure of
the original map. The above is invariant under the discrete translation
$\phi_{n}\rightarrow\phi_{n}+n\pi$; this is a consequence of the dynamical
$\mathbb{Z}_{2}$ symmetry of the original model.

Let us now build the equivalent of the map (17) for the generic
$\tilde{p}$-iterated map. For $T=0$ we have than that the action of the
$\tilde{p}$-iterated map is trivially
$\phi_{n+\tilde{p}}=\phi_{n}+2\tilde{p}\psi$, $I_{n+\tilde{p}}=I_{n}$. Then,
up to higher order in $T$ we can write:

$\begin{split}I_{n+\tilde{p}}&=I_{n}-\frac{\tilde{p}T}{2}(1-I^{2}_{n+p})U_{\tilde{p}}(\psi)\
\sin(2\phi_{n}+2(\tilde{p}-1)\psi)\ ,\\\
\phi_{n+\tilde{p}}&=\phi_{n}+2\tilde{p}\psi+\frac{\tilde{p}T}{2}I_{n+\tilde{p}}\left[1+U_{\tilde{p}}(\psi)\cos(2\phi_{n}+2(\tilde{p}-1)\psi)\right]\
,\end{split}$ (18)

where

$U_{\tilde{p}}(\psi)=\frac{\sin 2\tilde{p}\psi}{\tilde{p}\sin 2\psi}\ .$ (19)

Being $T$ small, in general the $2\tilde{p}\psi$ term in the r.h.s. of the
second equation od Eq (18) is going to dominate the evolution so that, we
expect to find ourselves in the quasi-periodic phase. However, as
$\tilde{p}\psi\sim k\pi$ for some integer $k$ we have that this terms becomes
small as well, signaling a the onset of the Poincaré-Birkhoff mechanism (for
$k$ odd, this is due to the $\mathbf{Z}_{2}$ symmetry, which allows us to
reabsorb the term by redefining $\phi_{n}\rightarrow\phi_{n}+n\pi$).

Let us put ourselves close to a resonance, i.e. let us consider the limit
$\psi=\psi_{r}+\delta\psi$ with $\delta\psi\ll\pi/p$. If $p$ is odd, the
smallest choice of $\tilde{p}$ which makes the term $\tilde{p}\psi$ small is
$\tilde{p}=p$; if $p$ is even, however, we have to choose $\tilde{p}=p/2$ (and
redefine $\phi_{n}\rightarrow\phi_{n}+n\pi$. In this limit the Eq. (18)
becomes

$\begin{split}I_{n+\tilde{p}}&=I_{n}-\frac{\tilde{p}T}{2}(1-I^{2}_{n+\tilde{p}})a_{r}\
\sin(2\phi_{n}-2\psi_{r})\ ,\\\
\phi_{n+\tilde{p}}&=\phi_{n}+2\tilde{p}\delta\psi+\frac{\tilde{p}T}{2}I_{n+\tilde{p}}\left(1+a_{r}\cos(2\phi_{n}-2\psi_{r})\right)\
,\end{split}$ (20)

where

$a_{r}=\begin{cases}1\ \ &\text{if}\ p=1,2\\\
2(-1)^{\tilde{p}-1}\csc(2\psi_{r})\delta\psi\ \ &\text{if}\ p\geq 3\
.\end{cases}$ (21)

Let us notice how in Eq. (20) the evolution of both $\phi$ and $I$ is now
slow, signaling that the $\tilde{p}$-iterated map can be approximated by a
continuous flow. In order to do so, we have to redefine the time scale
$\tilde{p}T\rightarrow T$, such that
$\phi_{n+\tilde{p}},I_{n+\tilde{p}}\rightarrow\phi_{n+1},I_{n+1}$ now by
introducing the time step $\Delta t=\tilde{p}T$, and expanding
$I_{n+1}=I_{n}+\dot{I}\ T+O(T^{2})$, $\phi_{n+1}=\phi_{n}+\dot{\phi}\
T+O(T^{2})$. We find then

$\begin{split}\dot{I}&=-\frac{a_{r}}{2}(1-I^{2})\sin(2\phi-2\psi_{r})\ ,\\\
\dot{\phi}&=2\tilde{p}\frac{\delta\psi}{T}+\frac{1}{2}I\left(1+a_{r}\cos(2\phi-2\psi_{r})\right)\
.\end{split}$ (22)

In turn, this can has the form of an Hamiltonian flow, generated by

$H(\phi,J)=2\frac{\delta\psi}{T}I-\frac{1}{4}(1-I^{2})\left(1+a_{r}\cos(2\phi-2\psi_{r})\right)\
.$ (23)

By taking into account our initial condition, namely $\phi(0)=\pi/2$, $J(0)=0$
we have that our $p$-iterated dynamics takes place along the curve
$H(I,\phi)=-\frac{1}{4}\left(1-a_{r}\cos(2\psi_{r})\right)$.

If this curve is bounded between two finite values of the angle $\phi$, the
dynamics will circle around a fixed point, signaling that we are within the
time crystalline phase. If, instead, the curves cover the all
$\phi\in[0,2\pi]$ interval, we have a quasi-periodic motion. Finally, close to
the separatrix between this two cases, chaos is expected to arise, so that the
boundary between these two regimes in terms of $T$, $\phi$ can be taken as an
estimate of the edge of the DFTC for small $T$.

For $p=1,2$, this criterion gives the condition $T=4|\delta\psi|$. For $p>2$,
instead, at the lowest order in $\delta\psi$ we find, for any
$\psi_{r}<\pi/2$, the condition

$T^{2}=\begin{cases}8\ \tilde{p}^{2}\tan\psi_{r}\ |\delta\psi|\ \ \
&\text{if}\ \delta\psi<0\ ,\\\ 8\ \tilde{p}^{2}\cot\psi_{r}\ |\delta\psi|\ \ \
&\text{if}\ \ \delta\psi>0\ .\\\ \end{cases}$ (24)

Let us notice how Eq.(24) captures both the non-analytic behavior of the
boundary of the DFTC and its lack of symmetry around the resonant value
$\psi_{r}$ for $r\neq 1/4$. However, in order for the result to be predictable
at the quantitative level we have to impose that the term $\delta\psi/T$ in
Hamiltonian (23) to be small, this implies that the steepest curve between the
two of Eq. (24) is not a good approximation of the boundary in the whole
region $|\delta\psi|\sim\pi/p$, and its coefficient is not reliable.

For $r=1/4$ ($\tilde{p}=2$) Eq. (24) gives

$T^{2}=32|\delta\psi|\ ,$ (25)

while for $r=1/6$ and $r=1/3$ ($\tilde{p}=3$) we find that the steepest edge
grows respectively as

$T^{2}=\pm 24\sqrt{3}\delta\psi\ .$ (26)

## References

  * Gisin  _et al._ [2002] N. Gisin, G. Ribordy, W. Tittel, and H. Zbinden, Quantum cryptography, Rev. Mod. Phys. 74, 145 (2002).
  * Giovannetti  _et al._ [2011] V. Giovannetti, S. Lloyd, and L. Maccone, Advances in quantum metrology, Nature Photonics 5, 222 (2011).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
  * Zhang  _et al._ [2017a] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator, Nature 551, 601 (2017a).
  * Jurcevic  _et al._ [2017] P. Jurcevic, H. Shen, P. Hauke, C. Maier, T. Brydges, C. Hempel, B. P. Lanyon, M. Heyl, R. Blatt, and C. F. Roos, Direct observation of dynamical quantum phase transitions in an interacting many-body system, Phys. Rev. Lett. 119, 080501 (2017).
  * Grifoni and Hänggi [1998] M. Grifoni and P. Hänggi, Driven quantum tunneling, Physics Reports 304, 229 (1998).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * Kaufman  _et al._ [2016] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Quantum thermalization through entanglement in an isolated many-body system, Science 353, 794 (2016), https://www.science.org/doi/pdf/10.1126/science.aaf6725 .
  * Mori  _et al._ [2018] T. Mori, T. N. Ikeda, E. Kaminishi, and M. Ueda, Thermalization and prethermalization in isolated quantum systems: a theoretical overview, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 112001 (2018).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Heyl  _et al._ [2013] M. Heyl, A. Polkovnikov, and S. Kehrein, Dynamical quantum phase transitions in the transverse-field ising model, Phys. Rev. Lett. 110, 135704 (2013).
  * Heyl [2018] M. Heyl, Dynamical quantum phase transitions: a review, Reports on Progress in Physics 81, 054001 (2018).
  * Žunkovič  _et al._ [2018] B. Žunkovič, M. Heyl, M. Knap, and A. Silva, Dynamical quantum phase transitions in spin chains with long-range interactions: Merging different concepts of nonequilibrium criticality, Phys. Rev. Lett. 120, 130601 (2018).
  * Karl  _et al._ [2017] M. Karl, H. Cakir, J. C. Halimeh, M. K. Oberthaler, M. Kastner, and T. Gasenzer, Universal equilibrium scaling functions at short times after a quench, Phys. Rev. E 96, 022110 (2017).
  * Wilczek [2012] F. Wilczek, Quantum time crystals, Phys. Rev. Lett. 109, 160401 (2012).
  * Else  _et al._ [2016] D. V. Else, B. Bauer, and C. Nayak, Floquet time crystals, Phys. Rev. Lett. 117, 090402 (2016).
  * Else  _et al._ [2020] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Discrete time crystals, Annual Review of Condensed Matter Physics 11, 467 (2020), https://doi.org/10.1146/annurev-conmatphys-031119-050658 .
  * Bordia  _et al._ [2017] P. Bordia, H. Lüschen, U. Schneider, M. Knap, and I. Bloch, Periodically driving a many-body localized quantum system, Nature Physics 13, 460 (2017).
  * Zhang  _et al._ [2017b] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, _et al._ , Observation of a discrete time crystal, Nature 543, 217 (2017b).
  * Choi  _et al._ [2017] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, _et al._ , Observation of discrete time-crystalline order in a disordered dipolar many-body system, Nature 543, 221 (2017).
  * Rovny  _et al._ [2018] J. Rovny, R. L. Blum, and S. E. Barrett, Observation of discrete-time-crystal signatures in an ordered dipolar many-body system, Physical review letters 120, 180603 (2018).
  * Russomanno  _et al._ [2017] A. Russomanno, F. Iemini, M. Dalmonte, and R. Fazio, Floquet time crystal in the lipkin-meshkov-glick model, Phys. Rev. B 95, 214307 (2017).
  * Khemani  _et al._ [2016] V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phase structure of driven quantum systems, Phys. Rev. Lett. 116, 250401 (2016).
  * Surace  _et al._ [2019] F. M. Surace, A. Russomanno, M. Dalmonte, A. Silva, R. Fazio, and F. Iemini, Floquet time crystals in clock models, Phys. Rev. B 99, 104303 (2019).
  * Haffner  _et al._ [2008] H. Haffner, C. Ross, and R. Blatt, Quantum computing with trapped ions, Phys. Rep. 469, 155 (2008).
  * Lahaye  _et al._ [2009] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, The physics of dipolar bosonic quantum gases, Rep. Prog. Phys. 72, 126401 (2009).
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with rydberg atoms, Rev. Mod. Phys. 82, 2313 (2010).
  * Ritsch  _et al._ [2013] H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger, Cold atoms in cavity-generated dynamical optical potentials, Rev. Mod. Phys. 85, 553 (2013).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 025001 (2021).
  * Mivehvar  _et al._ [2021] F. Mivehvar, F. Piazza, T. Donner, and H. Ritsch, Cavity qed with quantum gases: new paradigms in many-body physics, Adv. Phys. 70, 1 (2021).
  * Pagano  _et al._ [2018] G. Pagano, P. Hess, H. Kaplan, W. Tan, P. Richerme, P. Becker, A. Kyprianidis, J. Zhang, E. Birckelbaw, M. Hernandez, _et al._ , Cryogenic trapped-ion system for large scale quantum simulation, Quantum Science and Technology 4, 014004 (2018).
  * Campa  _et al._ [2014] A. Campa, T. Dauxois, D. Fanelli, and S. Ruffo,  _Physics of Long-Range Interacting Systems_ (Oxford Univ. Press, 2014).
  * Defenu  _et al._ [2021] N. Defenu, T. Donner, T. Macri, G. Pagano, S. Ruffo, and A. Trombettoni, Long-range interacting quantum systems (2021), arXiv:2109.01063 [cond-mat.quant-gas] .
  * Defenu  _et al._ [2017] N. Defenu, A. Trombettoni, and S. Ruffo, Criticality and phase diagram of quantum long-range o($n$) models, Phys. Rev. B 96, 104432 (2017).
  * Giachetti  _et al._ [2021] G. Giachetti, N. Defenu, S. Ruffo, and A. Trombettoni, Berezinskii-kosterlitz-thouless phase transitions with long-range couplings, Phys. Rev. Lett. 127, 156801 (2021).
  * Defenu  _et al._ [2019] N. Defenu, G. Morigi, L. Dell’Anna, and T. Enss, Universal dynamical scaling of long-range topological superconductors, Phys. Rev. B 100, 184306 (2019).
  * Halimeh  _et al._ [2020] J. C. Halimeh, M. Van Damme, V. Zauner-Stauber, and L. Vanderstraeten, Quasiparticle origin of dynamical quantum phase transitions, Phys. Rev. Research 2, 033111 (2020).
  * Acevedo  _et al._ [2014] O. L. Acevedo, L. Quiroga, F. J. Rodríguez, and N. F. Johnson, New dynamical scaling universality for quantum networks across adiabatic quantum phase transitions, Phys. Rev. Lett. 112, 030403 (2014).
  * Hwang  _et al._ [2015] M.-J. Hwang, R. Puebla, and M. B. Plenio, Quantum phase transition and universal dynamics in the rabi model, Phys. Rev. Lett. 115, 180404 (2015).
  * Defenu  _et al._ [2018] N. Defenu, T. Enss, M. Kastner, and G. Morigi, Dynamical critical scaling of long-range interacting quantum magnets, Phys. Rev. Lett. 121, 240403 (2018).
  * Van Regemortel  _et al._ [2016] M. Van Regemortel, D. Sels, and M. Wouters, Information propagation and equilibration in long-range kitaev chains, Phys. Rev. A 93, 032311 (2016).
  * Tran  _et al._ [2020] M. C. Tran, C.-F. Chen, A. Ehrenberg, A. Y. Guo, A. Deshpande, Y. Hong, Z.-X. Gong, A. V. Gorshkov, and A. Lucas, Hierarchy of linear light cones with long-range interactions, Phys. Rev. X 10, 031009 (2020).
  * Chen and Lucas [2019] C.-F. Chen and A. Lucas, Finite speed of quantum scrambling with long range interactions, Phys. Rev. Lett. 123, 250605 (2019).
  * Kuwahara and Saito [2020] T. Kuwahara and K. Saito, Strictly linear light cones in long-range interacting systems of arbitrary dimensions, Phys. Rev. X 10, 031010 (2020).
  * Defenu [2021] N. Defenu, Metastability and discrete spectrum of long-range systems, Proceedings of the National Academy of Sciences 118, e2101785118 (2021), https://www.pnas.org/doi/pdf/10.1073/pnas.2101785118 .
  * Giachetti and Defenu [2021] G. Giachetti and N. Defenu, Entanglement propagation and dynamics in non-additive quantum systems, arXiv preprint arXiv:2112.11488 (2021).
  * Lerose  _et al._ [2019] A. Lerose, J. Marino, A. Gambassi, and A. Silva, Prethermal quantum many-body kapitza phases of periodically driven spin systems, Phys. Rev. B 100, 104306 (2019).
  * Muñoz-Arias  _et al._ [2022] M. H. Muñoz-Arias, K. Chinni, and P. M. Poggi, Floquet time crystals in driven spin systems with all-to-all p-body interactions (2022), arXiv:2201.10692 [quant-ph] .
  * Machado  _et al._ [2020] F. Machado, D. V. Else, G. D. Kahanamoku-Meyer, C. Nayak, and N. Y. Yao, Long-range prethermal phases of nonequilibrium matter, Phys. Rev. X 10, 011043 (2020).
  * Collura  _et al._ [2021] M. Collura, A. D. Luca, D. Rossini, and A. Lerose, Discrete time-crystalline response stabilized by domain-wall confinement (2021), arXiv:2110.14705 [cond-mat.stat-mech] .
  * Pizzi  _et al._ [2021] A. Pizzi, J. Knolle, and A. Nunnenkamp, Higher-order and fractional discrete time crystals in clean long-range interacting systems, Nature Communications 12, 2341 (2021).
  * Sciolla and Biroli [2013] B. Sciolla and G. Biroli, Quantum quenches, dynamical transitions, and off-equilibrium quantum criticality, Physical Review B 88, 201110 (2013).
  * [55] See Supplemental Material.
  * Strogatz [2014] S. Strogatz,  _Nonlinear Dynamics and Chaos_ , Vol. 17 (Westview Press, 2014).
  * Hastings [2002] M. B. Hastings, Exact multifractal spectra for arbitrary laplacian random walks, Phys. Rev. Lett. 88, 055506 (2002).
  * Duplantier [2000] B. Duplantier, Conformally invariant fractals and potential theory, Physical Review Letters 84, 1363 (2000).
  * Kager and Nienhuis [2004] W. Kager and B. Nienhuis, A guide to stochastic löwner evolution and its applications, Journal of Statistical Physics 115, 1149 (2004).
  * Cardy [2005] J. Cardy, Sle for theoretical physicists, Annals of Physics 318, 81 (2005).
  * MacKay [1982] R. S. MacKay,  _Renormalisation in area preserving maps_ (Princeton University, 1982).
  * Vulpiani  _et al._ [2009] A. Vulpiani, F. Cecconi, and M. Cencini,  _Chaos: from simple models to complex systems_ , Vol. 17 (World Scientific, 2009).
  * Chirikov [1979] B. V. Chirikov, A universal instability of many-dimensional oscillator systems, Physics reports 52, 263 (1979).
  * Ribeiro  _et al._ [2008] P. Ribeiro, J. Vidal, and R. Mosseri, Exact spectrum of the lipkin-meshkov-glick model in the thermodynamic limit and finite-size corrections, Phys. Rev. E 78, 021106 (2008).
  * Auerbach [2012] A. Auerbach,  _Interacting electrons and quantum magnetism_ (Springer Science & Business Media, 2012).
  * Poincaré [1912] H. Poincaré, Sur un théoreme de géométrie, Rendiconti del Circolo Matematico di Palermo (1884-1940) 33, 375 (1912).
  * Birkhoff [1913] G. D. Birkhoff, Proof of poincaré’s geometric theorem, Transactions of the American Mathematical Society , 14 (1913).
  * Wisniacki  _et al._ [2011] D. A. Wisniacki, M. Saraceno, F. J. Arranz, R. M. Benito, and F. Borondo, Poincaré-birkhoff theorem in quantum mechanics, Phys. Rev. E 84, 026206 (2011).
  * Kolmogorov [1979] A. Kolmogorov, On the conservative of conditionally periodic motions for a small change in hamilton’s function (in russ.), dokl. acad. nauk sssr 98,(1954), 524–530; english transl, Lecture Notes in Physics 93, 51 (1979).
  * Arnold [2009] V. I. Arnold, Proof of a theorem of an kolmogorov on the invariance of quasi-periodic motions under small perturbations of the hamiltonian, Collected Works: Representations of Functions, Celestial Mechanics and KAM Theory, 1957–1965 , 267 (2009).
  * Möser [1962] J. Möser, On invariant curves of area-preserving mappings of an annulus, Nachr. Akad. Wiss. Göttingen, II , 1 (1962).