# Swapping Floquet time crystal

Roberto Gargiulo<EMAIL_ADDRESS>Present affiliation:
Forschungszentrum Jülich GmbH, Peter Grünberg Institute, Quantum Control
(PGI-8), 52425 Jülich, Germany  Dipartimento di Fisica “E. Pancini”,
Università di Napoli “Federico II”, Monte S. Angelo, I-80126 Napoli, Italy
Gianluca Passarelli  Dipartimento di Fisica “E. Pancini”, Università di Napoli
“Federico II”, Monte S. Angelo, I-80126 Napoli, Italy  Procolo Lucignano
Dipartimento di Fisica “E. Pancini”, Università di Napoli “Federico II”, Monte
S. Angelo, I-80126 Napoli, Italy  Angelo Russomanno  Present affiliation:
Dipartimento di Fisica “E. Pancini”, Università di Napoli “Federico II”, Monte
S. Angelo, I-80126 Napoli, Italy  Scuola Superiore Meridionale, Università di
Napoli Federico II, Largo San Marcellino 10, I-80138 Napoli, Italy

###### Abstract

We propose a time-crystal model based on a disordered interacting long-range
spin chain where the periodic swapping of nearby spin couples is applied. This
protocol can be applied to systems with any local spin magnitude $s$ and in
principle also to systems with nonspin (fermionic or bosonic) local Hilbert
space. We explicitly consider the cases $s=1/2$ and $s=1$, using analytical
and numerical methods to show that the time-crystal behavior appears in a
range of parameters. In particular, we study the persistence of period-
doubling oscillations in time, the properties of the Floquet spectrum
($\pi$-spectral pairing and correlations of the Floquet states), and introduce
a quantity (the local imbalance) to assess what initial states give rise to a
period-doubling dynamics. We also use a probe of quantum
integrability/ergodicity to understand the interval of parameters where the
system does not thermalize, and a nontrivial persistent period-doubling
behavior is possible.

##  I Introduction

The experimental discovery [1, 2] of time-crystals few years after their
theoretical prediction [3, 4, 5] has been a real breakthrough. In analogy to
ordinary crystals, time crystals appear as a consequence of breaking time-
translation symmetry in the system [6, 7, 8, 9, 10]. Following earlier
attempts to identify systems able to display time-translation symmetry
breaking, a no-go theorem showed that this is not possible in the ground state
or in thermal equilibrium [11].

Among many possible non-equilibrium candidates, quantum periodically-driven
(Floquet) systems have proven to be the most promising realization. Here the
discrete time-translation symmetry is broken in the thermodynamic limit, by
the appearance of a response with a period multiple than the one of the
driving. Stimulated by the initial proposals [12, 13], a large body of
theoretical work has been performed [14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24, 25, 26, 27]. A common ingredient is the presence in the dynamics of a
sufficient number of constraints that introduce ergodicity-breaking, thus
preventing infinite-temperature thermalization, for instance disorder that
induces many-body localization (MBL); see [28] for a review.

In this context, the properties of the Floquet states [29, 30] – the
eigenstates of the periodically-driven dynamics – are crucial. In order for a
period-doubling Floquet time crystal to occur, all of the Floquet states must
be organized in doublets of cat states. Cat states are superpositions of two
macroscopically different classical spin configurations, and the time crystal
appears as a form of Rabi oscillation between these configurations, with
frequency provided by the splitting of the doublet. In the thermodynamic
limit, the oscillation frequencies of all the doublets become synchronized,
and this synchronization occurs at period double than the driving
($\pi$-spectral pairing), resulting in a time crystal. In the original
proposal all the Floquet spectrum was $\pi$-spectral paired, but there are
also works where only a fraction of the Floquet spectrum is paired, and a
period-doubling behavior for specific initial states can be observed [31, 32,
33].

Following these results, we propose a new periodic-driving protocol, which is
applied to a disordered spin chain with long-range interactions, and leads to
a Floquet time crystal. The drive consists in periodically swapping two
neighbouring spins, instead of flipping each individual spin (see Fig. 1). As
in the spin-flipping case, this protocol can be implemented by periodically
applying in sequence two different Hamiltonians, and can be applied to
interacting chains of particles with arbitrary value of the local spin
magnitude $s$, rather than only to spin-1/2 systems. Moreover, the swap
protocol has a straightforward physical meaning in systems with fermionic or
bosonic local Hilbert space, since it conserves particle number. Therefore it
has a potential range of applicability that goes beyond spin systems.

We find that a Floquet time crystal is realized in this case, although in a
slightly weaker sense than what has been found for the spin-flip driving.
While in Refs. [12, 13] all the Floquet spectrum is $\pi$-spectral paired,
here a fraction of the quasienergies does not have this property.
Nevertheless, this fraction tends exponentially to zero with the system size,
and we can observe persisting and robust period-doubling oscillations for
macroscopically many initial states.

The paper is organized as follows. In Sec. II we introduce our model
Hamiltonians – both for local spins with magnitude 1/2 and for generic
magnitudes – and study how the period doubling appears in the solvable case, a
special point in the parameter space where analytical computations can be
made.

In Sec. III we introduce the local imbalance, a quantity that allows
characterizing the initial classical spin configurations that give rise to a
period-doubling dynamics. We find that increasing the local spin magnitude,
the distribution of the imbalances on the spin configurations gets a larger
average and becomes narrower, marking that a period-doubling behavior occurs
more easily for a random initial spin configuration.

In Sec. IV we start with our numerical analysis, that we perform for the
models with local spin $s=1/2$ (spin-1/2 case) and $s=1$ (spin-1 case). The
numerical analysis is important in order to assess if the period-doubling
behavior is robust and persists also beyond the solvable point, which would
indicate a proper time-crystal phase supported by the interactions, rather
than an isolated point. We study the period doubling in the time domain,
looking at a collective quantity that is nonvanishing whenever period-doubling
oscillations of the local magnetizations are present.

Moving away from the solvable point, period doubling oscillations are still
there, but last a finite time that exponentially increases with the system
size. In this way period-doubling oscillations become persistent in the large-
size limit, as appropriate for a time crystal, for which the time-translation
symmetry breaking appears only in the thermodynamic limit, in analogy with
standard symmetry breaking [34]. In doing our analysis, we emphasize that for
some initial states only part of the onsite magnetizations can show period
doubling.

In order to see a nontrivial period-doubling behavior of the observables,
there must be no energy absorption and no infinite-temperature thermalization
(in this sense the dynamics must be regular or integrable-like). In order to
find the range of parameters where this occurs, we focus in Sec. IV.3 on the
average level-spacing ratio, a standard probe of the integrability/ergodicity
properties of the quantum dynamics [35, 36]. We find an interval of the
parameters where the dynamics is integrable-like, and see that this effect
becomes more marked with increasing system size. In this interval there is no
thermalization and we can see the period-doubling oscillations. We argue that
the physical effect leading to this breaking of ergodicity is disorder,
inducing MBL.

In Sec. V we study the $\pi$-spectral-pairing properties of the Floquet
states, and see that this phenomenon is strictly related to period doubling.
The $\pi$-spectral pairing can be found analytically at the critical point,
and we numerically find that it is robust in a well-defined range of
parameters that marks the time-crystal phase. This robustness is related to
the disordered long-range interactions that break the degeneracies in the
Floquet spectrum at the solvable point, and make the structure of Floquet
states and quasienergies robust under local perturbations (outside of the
solvable point, the $\pi$-spectral pairing appears rigorously only in the
thermodynamic limit. At finite size there is an error leading to the finite-
time duration of the period-doubling oscillations – see Appendix B).

In this section we also directly study the cat-state properties of the Floquet
states by using a global quantifier of the correlations. (These cat-state
properties are strictly related to $\pi$-spectral pairing.) The Floquet states
with the cat-state property are the correlated ones, and at the solvable point
we find that are those with a finite local imbalance (for $s=1/2$). Moving
away from the solvable point, we numerically show that the total amount of
correlations for an average Floquet state still increases with system size,
which indicates the presence of long-range correlations and further confirms
the robustness of the time crystal phase.

In Appendix A we show how the spin flipping for the spin-1/2 case occurs, in
Appendix B we discuss the deep relation between $\pi$-spectral pairing and
period doubling, in Appendix C we consider a different range of the
interactions, in Appendix D we discuss in detail the imbalance distributions
of the classical spin configurations for different values of $s$, and in
Appendix E we discuss the relation between local imbalance and correlations
for $s=1/2$ at the solvable point. In Sec. VI we draw our conclusions.

##  II Model Hamiltonians

###  II.1 Spin-1/2 model

####  II.1.1 Description of the model

We study a disordered spin-1/2 chain with power-law $ZZ$-interactions,
nearest-neighbour $XX$-interactions and a longitudinal field which is subject
to a periodic kick ($\hbar=1=T$). The periodic time-dependent Hamiltonian is

$\hat{H}(t)=\hat{H}_{\text{int}}+\hat{K}\sum_{n}\delta(t-n)\,,$ (1)

where the interaction Hamiltonian and the kicking Hamiltonian are respectively

$\displaystyle\hat{H}_{\text{int}}$
$\displaystyle\equiv-\sum_{k=1}^{L-1}J(\hat{\sigma}_{k}^{x}\hat{\sigma}_{k+1}^{x}+\hat{\sigma}_{k}^{y}\hat{\sigma}_{k+1}^{y})+\sum_{k=1}^{L-1}\sum_{q=k+1}^{L}V_{kq}\hat{\sigma}_{k}^{z}\hat{\sigma}_{q}^{z}+$
$\displaystyle+\sum_{k=1}^{L}h_{k}^{z}\hat{\sigma}_{k}^{z}\,,$
$\displaystyle\hat{K}$
$\displaystyle\equiv\left(\frac{\pi}{4}+\varepsilon\right)\sum_{k=1}^{L/2}\vec{\boldsymbol{\sigma}}_{2k-1}\cdot\vec{\boldsymbol{\sigma}}_{2k}\,,$
(2)

where $\epsilon\neq 0$ describes imperfect swaps. So one can construct the
time-evolution operator over one period (the so-called Floquet operator) as

$\hat{U}_{F}=e^{-\boldsymbol{i}\hat{K}}e^{-\boldsymbol{i}\hat{H}_{\text{int}}}\,.$
(3)

Here $\hat{\sigma}_{k}^{\beta=x,y,z}$ are the Pauli matrices; $V_{kq}$ follows
a power-law behavior plus disorder,

$V_{kq}=\frac{1}{\mathcal{N}_{L,\alpha}}\frac{\tilde{V}_{kq}}{|k-q|^{\alpha}},$
(4)

where $\mathcal{N}_{L,\alpha}$ is a suitable Kac normalization constant [37]
which makes the energy extensive defined as

$\mathcal{N}_{L,\alpha}=\begin{cases}1&\text{if $\alpha\geq 1$;}\\\ \ln
L&\text{if $\alpha=1$;}\\\ L^{1-\alpha}&\text{if $\alpha<1$.}\end{cases}$ (5)

$\tilde{V}_{kq},h_{k}$ are independent random variables taken from the uniform
box distributions $\tilde{V}_{kq}\in[V/2,3V/2]$, $h_{k}\in[-h,h]$.

Let us focus on the kicking Hamiltonian $\hat{K}$, with the associated unitary
transformation corresponding to one kick $\hat{P}=e^{-\boldsymbol{i}\hat{K}}$.
The operator $\hat{P}$ can be factored into a product of unitary operators
acting on neighboring sites (up to an irrelevant phase factor – see Appendix
A)

$e^{-\boldsymbol{i}\hat{K}}=\prod_{k=1}^{L/2}\frac{1+\vec{\boldsymbol{\sigma}}_{2k-1}\cdot\vec{\boldsymbol{\sigma}}_{2k}}{2}\,.$
(6)

The operator
$(1+\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j})/2$ is the
swap operator on two sites [38]. This acts on a two-site product state by
exchanging the states of the two sites (see Appendix A for a discussion)

$\frac{1+\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}{2}\ket{\psi_{i},\psi_{j}}=\ket{\psi_{j},\psi_{i}}\,.$
(7)

Therefore, the operator $\hat{P}$ exchanges the states of all odd sites with
their following even sites. Since this operation is well defined only for
pairs of spins, we restrict ourselves to the case of $L$ even. As we shall
see, this protocol will show several differences in terms of both dynamics and
spectrum with respect to the spin-flip protocol, such that it is not trivial
to relate the two models. Let us start considering a special point in the
parameter space where the dynamics is exactly solvable.

####  II.1.2 Analysis of the solvable case ($J=\varepsilon=0$)

To understand the behavior of the model at the level of the dynamics and
spectrum, it is instructive to look at the solvable case $J=\varepsilon=0$. In
this case, $\hat{H}_{\text{int}}$ commutes with all $\hat{\sigma}_{k}^{z}$ and
$e^{-\boldsymbol{i}\hat{K}}$ is a perfect swap. One can show that there is an
exact period-doubling dynamics, which can be seen at the level of the
operators themselves

$\hat{\sigma}_{k}^{z}(T)=\hat{U}_{F}^{\dagger}\hat{\sigma}_{k}^{z}\hat{U}_{F}=\hat{\sigma}_{\overline{k}}^{z},$
(8)

where $\overline{k}=k+1$ if $k$ is odd and $\overline{k}=k-1$ is $k$ is even.
Therefore, the magnetization serves as an order parameter, and period-doubling
will be observed for all initial states where, for some $k$,
$\braket{\hat{\sigma}_{k}^{z}}\neq\braket{\hat{\sigma}_{\overline{k}}^{z}}$,
since the local magnetization will oscillate between its initial value and
that of its paired neighbour. We graphically show this behavior in Fig. 1 and
compare it with the more known case of the spin flipping.

We note that the period-doubling will be a measurable effect in the
thermodynamic limit only if a finite fraction of sites satisfies the above
condition, i. e., there is an extensive quantity of not aligned spins.

As an example, one can imagine an initial state where the odd spins lie in the
$xy$ plane – such that $\braket{\hat{\sigma}_{2k-1}^{z}(t=0)}=0$ – and the
even spins lie in the $z$ direction – and so
$\braket{\hat{\sigma}_{2k}^{z}(t=0)}=\pm 1$. For the odd ones, the interaction
part results in a dephasing of the spins in the $xy$-plane [39] and the
$z$-magnetization remains zero. For the even ones, being $J=0$, the
interaction part of the Hamiltonian does not change the spins pointing along
the $z$ direction. So, due to the periodic swap, the magnetization at each
site will exhibit period-doubling by oscillating between $0$ and $\pm 1$.

Furthermore, certain configurations may display period-doubling restricted
only to a part of the chain, such as one which is part anti-ferromagnetic
($\\{+-,+-,\cdots\\}$) and part ferromagnetic ($\\{++,++,\cdots\\}$). In this
case the order parameter $\braket{\hat{\sigma}_{k}^{z}(t)}$ oscillates only in
a certain part of the chain. Nevertheless, if the anti-ferromagnetic part is
extensive, then we still consider such a state as breaking the time-
translation symmetry of the system.

We highlight that the number of spin configurations where no period-doubling
takes places increases with system size as $2^{L/2}$ and the number of
remaining configurations is $2^{L}-2^{L/2}$. As a consequence, the
configurations which display no period-doubling throughout the entire chain
must satisfy a very special condition, as their fraction vanishes with
increasing system size as $2^{-L/2}$ over the total number of spin
configurations.

Thus, we find that this swapping dynamics is in general quite different than
the prototypical model with spin-flip [40], where the magnetization always
(and only) changes sign, independently of neighbouring spins (see Fig. 1).

Figure 1: Comparison between swap and flip protocol. In the Flip case, all
spins with non-zero magnetization are affected. In the swap case, all pairs of
spins that are not aligned are affected.

###  II.2 Spin-1 model and beyond

The interesting aspect of spin swapping is that one can generalize the time-
crystal behavior of the spin-1/2 model to generic local spin magnitude. The
price to pay is that, for increasing spin magnitude $s$, the spin-swapping
operator becomes more and more complicated [38].

Let us focus here on the spin-1 case ($s=1$). The Hamiltonian for this case
has the same structure as Eq. (1). The interaction part has much the same form
as the spin-1/2 case, although the periodic swap also involves higher powers
of the spin operators

$\begin{split}\hat{H}_{\text{int}}&=-\sum_{k=1}^{L-1}J(\hat{S}_{k}^{x}\hat{S}_{k+1}^{x}+\hat{S}_{k}^{y}\hat{S}_{k+1}^{y})\\\
&\quad{}+\sum_{k=1}^{L-1}\sum_{q=k+1}^{L}\hat{S}_{k}^{z}\hat{S}_{q}^{z}+\sum_{k=1}^{L}h_{k}^{z}\hat{S}_{k}^{z},\\\
\hat{K}&=\left(\frac{\pi}{2}+\varepsilon\right)\sum_{k=1}^{L/2}\left[(\vec{\mathbf{S}}_{2k-1}\cdot\vec{\mathbf{S}}_{2k})^{2}+\vec{\mathbf{S}}_{2k-1}\cdot\vec{\mathbf{S}}_{2k}\right],\\\
\hat{U}_{F}&=e^{-\boldsymbol{i}\hat{K}}e^{-\boldsymbol{i}\hat{H}_{\text{int}}}.\end{split}$
(9)

and the operators $\hat{S}_{k}^{\beta=x,y,z}$ are now spin-1 operators.

$V_{kq},h_{k}$ have the same distributions as in the spin-1/2 case. Of course,
even though the Hamiltonian is similar in form, the spin-1 chain has a
different behavior with respect to the spin-1/2, due to the additional state
with $s_{k}=0$ [41].

The main difference with the spin-1/2 case is that the number of classical
configurations providing period doubling is given by $3^{L}-3^{L/2}$ [for
generic local spin magnitude $s$ it is $(2s+1)^{L}-(2s+1)^{L/2}$]. So, the
ratio of configurations that provide no period doubling is for generic spin
$(2s+1)^{-L/2}$. It is smaller than in the spin-1/2 case, and it decreases for
increasing value of the local spin $s$. Furthermore, for a random spin
configuration there will be a larger fraction of sites which participate to
the period-doubling (see see Appendix D).

It is important to define the spin basis $\ket{\\{s_{k}\\}}$. In the spin-1/2
case this basis is such that
$\hat{\sigma}_{k}^{z}\ket{\\{s_{k}\\}}=s_{k}\ket{\\{s_{k}\\}},\;s_{k}\in\\{+,-\\}$,
while in the generic case
$\hat{S}_{k}^{z}\ket{\\{s_{k}\\}}=s_{k}\ket{\\{s_{k}\\}},\;s_{k}\in\\{-s,-s+1,\ldots,s-1,s\\}$.
Another important property of our models is that the dynamics preserves the
total $z$ magnetization
($\hat{\mathcal{S}}_{z}=\frac{1}{2}\sum_{k=1}^{L}\hat{\sigma}_{k}^{z}$ for the
spin-1/2 case, $\hat{\mathcal{S}}_{z}=\sum_{k=1}^{L}\hat{S}_{k}^{z}$
otherwise), which will allow us to restrict our numerics to a subspace with
fixed $\mathcal{S}_{z}$. In the following are we going to study analytically
and numerically different aspects of the period-doubling behavior.

##  III Local imbalance

In order to easily identify which initial states lead to period-doubling we
define the following operator, which we call Local Imbalance
($\hat{\mathcal{I}}_{\rm LI}$). Our aim is to find an operator that, on
classical spin configurations, counts how many pairs of swapped spins have
different values of the magnetization
$\braket{\\{s_{k}\\}}{\hat{\mathcal{I}}_{\rm
LI}}{\\{s_{k}\\}}=\frac{1}{L/2}\sum_{k}(1-\delta_{s_{2k-1},s_{2k}})$. In the
spin-1/2 case, this operator is

$\hat{\mathcal{I}}_{\rm
LI}=\frac{1}{L/2}\sum_{k=1}^{L/2}\frac{(\hat{\sigma}_{2k}^{z}-\hat{\sigma}_{2k-1}^{z})^{2}}{4}\,,$
(10)

while, in the spin-1 case, it is

$\displaystyle\hat{\mathcal{I}}_{\rm LI}$
$\displaystyle=\frac{1}{L/2}\sum_{k=1}^{L/2}\big{[}(\hat{S}_{2k-1}^{z})^{2}+(\hat{S}_{2k}^{z})^{2}$
$\displaystyle-\frac{1}{2}\hat{S}_{2k-1}^{z}\hat{S}_{2k}^{z}(1+3\hat{S}_{2k-1}^{z}\hat{S}_{2k}^{z})\big{]}\,.$
(11)

This quantity – inspired by the charge-imbalance in fermionic chains [42, 43]
– measures the difference in alignment (spin analogue of the fermionic
occupation) in all the $L/2$ even-odd pairs, regardless of their sign (hence
the square).

Let us first focus on the solvable case, and consider the behavior of
$\hat{\mathcal{I}}_{\rm LI}$ on the spin basis $\ket{\\{s_{k}\\}}$. When
$\braket{\\{s_{k}\\}}{\hat{\mathcal{I}}_{\rm LI}}{\\{s_{k}\\}}=1$, pairs have
different spin values ($s_{2k}\neq s_{2k-1}$), the swap makes spins oscillate
between two values, and the dynamics is the same as for the spin-flip
protocol. When $\braket{\\{s_{k}\\}}{\hat{\mathcal{I}}_{\rm
LI}}{\\{s_{k}\\}}=0$, all pairs have the same spin value, and the state is
left unchanged by the periodic swap. In the general case
$0<\braket{\\{s_{k}\\}}{\hat{\mathcal{I}}_{\rm LI}}{\\{s_{k}\\}}<1$ only
certain pairs will give rise to period-doubling. This is in contrast with the
spin-flip case, where any classical configuration gives rise to period-
doubling oscillations at the solvable point.

Much of this also holds for states close to the spin basis. For instance,
focus on the spin-1/2 case and consider the product state where spins have
alignment close to the $z$ direction, i. e.,
$\bigotimes_{k=1}^{L}(\cos(\theta_{k})\ket{+}+\sin(\theta_{k})\ket{-})$ for
$\theta_{k}$ equal to $0$ or $\pi$ up to small corrections $\delta\theta\ll
1$. This state clearly displays period-doubling in every site and at leading
order in $\delta\theta$ we have $\mathcal{I}_{\rm LI}=1-C(\delta\theta)^{2}$
with $C$ a constant. Therefore, in this case, the high value of local
imbalance accurately predicts the presence of period-doubling.

Notice that at the solvable point $\mathcal{I}_{\rm LI}$ commutes with
$U_{F}$, meaning that it is provides also a good quantum number for the
eigenstates of $U_{F}$ and a useful signature of their properties, as we shall
see below. It is also important to remark that, for increasing local spin
magnitude $s$, the average local imbalance of the states in the spin basis
increases, and the variance decreases, as we show in detail in Appendix D.
This is interesting because initializing the system with a spin state with
nonvanishing imbalance leads to period-doubling oscillations lasting forever,
at least for the solvable case.

In the next section we move away from the solvable point and show that there
are still period-doubling oscillations that last for a time increasing with
the system size, so we have not a special isolated point but a robust time-
crystal phase.

##  IV Period-doubling oscillations

A time-crystal behavior is signaled by persistent period-doubling oscillations
of some observables. In order to speak about a phase and not an isolated
point, these oscillations must be robust to changes of the parameters. Let us
consider what happens when we move from the solvable point and take $J\neq 0$.
We are going to show that there is an interval of $J$ where period-doubling
oscillations appear, and last a time that exponentially increases with the
system size, for many possible initial states. We see that both in the case of
spin-1/2 model (Sec. IV.1), and spin-1 model (Sec. IV.2).

In the current and following sections we shall consider the disorder average
of various quantities, defined as
$\overline{(\ldots)}=\frac{1}{N_{d}}\sum_{i=1}^{N_{d}}(\ldots)$. The chosen
number of disorder realizations $N_{d}$ will typically get smaller with $L$,
since we take advantage of the self-averaging behaviour - specifically, we
take $20480\cdot 2^{1-L/2}$, unless otherwise mentioned, which allows for
sufficiently small statistical errors.

Furthermore, we fix $V=3$ and $h_{z}=16$ in all our computations and mostly
focus on the perfect swap case $\varepsilon=0$ \- except in certain cases,
which we explicitly mention. We shall consider two choices of interaction
range: $\alpha=0.5$ (long-range) and $\alpha=3$ (short-range). Since we have
found that the two values typically have only very slight quantitative
differences, we will mostly focus on the $\alpha=0.5$ case, and leave the
$\alpha=3$ case in Appendix C.

In the entire section we will focus on stroboscopic dynamics, i. e., the time
$t=nT$ will be always an integer multiple ($n\in\mathbb{N}$) of the period
$T=1$.

###  IV.1 Spin-1/2 case

In this section we consider the problem of period-doubling dynamics by looking
at two initial states in the spin configuration with finite $\mathcal{I}_{\rm
LI}$, such that
$\braket{\hat{\sigma}_{k}^{z}}\neq\braket{\hat{\sigma}_{\overline{k}}^{z}}$
over a finite fraction of the total sites. Specifically, we consider an
initial Néel state $\ket{+-+-\cdots}$ (fully anti-ferromagnetic) and a “half-
Néel” state $\ket{+-+-\cdots++++\cdots}$ (anti-ferromagnetic on the length
$l_{\mathrm{AF}}=2\lfloor L/4\rfloor$, and the rest ferromagnetic), as
pictured in Fig. 2.

Figure 2: Initial Néel (upper panel) and Half-Néel (lower panel), with their
respective dynamics, which are exact only in the solvable case.

With the first state one gets the same kind of dynamics as the one observed in
the spin-flip protocol, since all spins get inverted after each kick. Indeed,
it has maximal local imbalance $\mathcal{I}_{\rm LI}=1$. With the second
initial state, “half-Néel”, instead one gets period-doubling dynamics only in
the first half of the chain and $\mathcal{I}_{\rm LI}=\lfloor
L/4\rfloor/(L/2)$.

In order to monitor the period-doubling dynamics obtained with these initial
states, we consider the following quantity, analogous to the one used in [40]

$\begin{split}Z(t)=\frac{1}{L/2}\sum_{k=1}^{L/2}\operatorname{sign}(\braket{\hat{\sigma}_{2k-1}^{z}(0)}-\braket{\hat{\sigma}_{2k}^{z}(0)})\times\\\
\times(\braket{\hat{\sigma}_{2k-1}^{z}(t)}-\braket{\hat{\sigma}_{2k}^{z}(t)})\,.\end{split}$
(12)

$Z(t)$ is defined in such a way that it is always non-negative at $t=0$ and at
the solvable point ($J=\varepsilon=0$) oscillates between a positive and
negative value every period $Z(t)=(-1)^{t}Z(0)$, displaying in this way period
doubling. It has zero contributions from regions where
$\braket{\hat{\sigma}_{2k-1}^{z}}=\braket{\hat{\sigma}_{2k}^{z}}$ and due to
the normalization factor $L/2$, only states with an extensive quantity of
unaligned spin pairs have a finite $Z(0)$ in the thermodynamic limit.

Let us average over disorder and focus on the quantity
$(-1)^{t}\overline{Z}(t)$, that in presence of period doubling is constantly
positive and does not change sign at every period. We show an example of
$(-1)^{t}\overline{Z}(t)$ versus $t$ with the initial Néel state for many
system system sizes is shown in Fig. 3(a). We can qualitatively see that
$(-1)^{t}\overline{Z}(t)$ decays over a time scale that increases with system
size.

\begin{overpic}[width=227.62204pt]{Z_avg_0-5}\put(0.0,95.0){(a)}\put(0.0,55.0){(b)}\end{overpic}
Figure 3: Spin-1/2 model. (Panel a) Dynamics of $(-1)^{t}\overline{Z}(t)$ for
initial Néel state and different system sizes. (Panel b) Dynamics of
$(-1)^{t}\overline{Z}(t)$ for initial Half-Néel state (solid line) and Néel
state (dashed line) and different system sizes. Initial values have all been
normalized to $+1$ to make the comparison clearer. Numerical parameters:
$J=0.1$, $\varepsilon=0.01$.

Let us study this behavior of the decay time more quantitatively. For a given
disorder realization, we define the decay time $\tau$ as the first time where
$Z(t)$ no longer manifests period-doubling behaviour, i. e.,
$Z(\tau-1)Z(\tau)>0$ (if multiple periods are skipped the condition becomes
$(-1)^{n}Z(\tau-n)Z(\tau)<0$). In Fig. 4 we show $\overline{\tau}$ versus $L$
for different system sizes. We notice that the decay time exponentially
increases with system size for a range of parameters. This is a signature of
stable period-doubling in the thermodynamic limit [44, 45], meaning that there
is time crystal behaviour in this range of parameters.

Figure 4: Spin-1/2 model. Decay times $\overline{\tau}$ of $(-1)^{t}Z(t)$ for
an initial Néel state. Here we have considered time evolutions while skipping
over $n\geq 1$ periods, such that the time step stays below $\sim 1\%$ of the
standard deviation of $\tau$.

We show $\overline{Z}(t)$ versus $t$ initializing with the half-Néel state in
Fig. 3(b). We notice that the curves do not seem to depend on the total length
of the chain $L$, but rather on the length of the anti-ferromagnetic part
$l_{\mathrm{AF}}=2\lfloor L/4\rfloor$. Each curve for the Néel state closely
follows the one for the half Néel state if the antiferromagnetic parts are of
the same size.

It is also worth noticing that even though $Z(t)$ decays to zero at long
times, the local magnetizations $\braket{\hat{\sigma}_{k}^{z}(t)}$ themselves
do not: The ones in the ferromagnetic part are frozen. This is an effect of
many-body localization (see Sec. IV.3) that hinders energy absorption and
infinite-temperature thermalization and induces space localization of
excitations.

Figure 5: Spin-1/2 model. Comparison of
$\braket{\overline{\hat{\sigma}_{k}^{z}}(t)}$ at $k=L-\lfloor L/4\rfloor$
(dashed line) and $(-1)^{t}\braket{\overline{\hat{\sigma}_{k}^{z}}(t)}$ at
$k=\lfloor L/4\rfloor$ (solid line), up to $t=10^{4}$, for an initial Half-
Néel state. Parameters: $J=0.1$

###  IV.2 Spin-1 case

Let us now look at the dynamics of the spin-1 chain, where once again we
require that in the initial state
$\braket{\hat{S}_{k}^{z}}\neq\braket{\hat{S}_{\overline{k}}^{z}}$ for an
extensive number of sites. Specifically, we consider the initial state “Up-
Zero” given by the alternating sequence $\ket{+1,0,+1,0,\cdots}$. This state
gives an example of a “non-flip” dynamics already within the spin
configuration basis, in the sense that each site oscillates between
$\braket{\hat{S}_{k}^{z}(t)}=0$ and $\braket{\hat{S}_{k}^{z}(t)}=+1$ (solvable
case) and does not reduce to a spin-flip. Furthermore, since the local
magnetization changes value at each site, its local imbalance [Eq. (III)] has
maximal value $1$. Once again, in order to study the dynamics of the system we
use the observable $Z(t)$ defined in Eq. (12), with $\hat{S}_{k}^{z}$ instead
of $\hat{\sigma}_{k}^{z}$. The behavior of its disorder average
$\overline{Z}(t)$ versus $t$ for different system sizes is shown in Fig. 6(a).

\begin{overpic}[width=227.62204pt]{sigmaz_Dynamics_two_panel_1-2_0-5}\put(0.0,95.0){(a)}\put(0.0,55.0){(b)}\end{overpic}
Figure 6: Spin-1 model. (Panel a) Dynamics of $(-1)^{t}\overline{Z}(t)$ for an
initial Up-Zero state. (Lower panel) Dynamics of
$(-1)^{t}(\braket{\overline{\hat{S}_{k}^{z}(t)}}-0.5)+0.5$ at $k=2\lfloor
L/4\rfloor$ (solid line) and $k=2\lfloor L/4\rfloor+1$ (dashed line) for an
initial Up-Zero state. Parameters: $J=0.1$.

We find again a range of parameters where $\overline{Z}(t)$ oscillates between
a positive and negative value and decays to zero in a time exponentially
increasing with the system size, as we can see in the behavior of
$\overline{\tau}$ versus $L$ shown in Fig. 7. Again, this means that there is
time crystal behaviour in this range of parameters.

Figure 7: Spin-1 model. Decay times $\overline{\tau}$ of $(-1)^{t}Z(t)$,
averaged over disorder for an initial Up-Zero state.

The period-doubling oscillations and the fact that they last for a time
exponential in the system size can be seen also at the level of the
expectations of the onsite magnetizations [see Fig. 6(b)]. Specifically, by
taking the middle of the chain as representative, we see that both sites
(starting at either $\braket{S_{k}^{z}}=0$ or $1$) oscillate around
$\braket{S_{k}^{z}}=0.5$ and decay towards this middle value at exponentially
long times.

###  IV.3 Regular dynamics

In order to see persistent period doubling, there must be no thermalization of
local observables. (Local thermalization is a common behavior of many-body
quantum systems [46]). The dynamics of the system must be regular, that’s to
say similar to an integrable quantum system. Here this effect can be provided
by disorder, that is known to break ergodicity and induce many-body
localization [47, 28, 48].

In order to check the range of $J$ where integrable-like behavior occurs, we
diagonalize the Floquet operator defined in Eq. (3) as
$\hat{U}_{F}\ket{\psi_{\beta}}=e^{-\boldsymbol{i}\mu_{\beta}}\ket{\psi_{\beta}}$.
The eigenvectors $\ket{\psi_{\beta}}$ are the so-called Floquet states and the
phases of the eigenvalues $\mu_{\beta}$ are the corresponding quasienergies
[30, 29]. Let us evaluate the average level spacing ratio $\braket{r}$ [35]
defined as

$\braket{r}=\frac{1}{\mathcal{N}}\sum_{\beta=1}^{\mathcal{N}_{M}-2}\overline{\left[\frac{\min(\mu_{\beta+2}-\mu_{\beta+1},\mu_{\beta+1}-\mu_{\beta})}{\max(\mu_{\beta+2}-\mu_{\beta+1},\mu_{\beta+1}-\mu_{\beta})}\right]}$
(13)

where quasienergies are restricted to the first quasienergy Brillouin zone
[49, 50] $[-\pi,\pi]$ (they are periodic of period $2\pi$) and taken in
increasing order. $\mathcal{N}$ is the dimension of the sector of the Hilbert
space with vanishing total magnetization, $s_{z}=0$. We consider in Eq. (13)
only the quasienergies $\mu_{\beta}$ corresponding to Floquet states in this
sector. In this way we are restricting to an irreducible invariant subspace of
the Hamiltonian, a condition required in order for the distribution of the
level spacings $\delta_{\beta}=\mu_{\beta+1}-\mu_{\beta}$ (and the related
ratio $\braket{r}$) to be a meaningful ergodicity indicator [51].

When the driven system is ergodic, i. e., locally thermalizing with infinite
temperature [52], the Floquet operator $\hat{U}_{F}$ belongs to the circular
orthogonal ensemble (COE) of symmetric unitary matrices (because of the time-
reversal invariance) [53, 54, 55, 36]. In this case, the distribution of the
level spacings $\delta_{\beta}=\mu_{\beta+1}-\mu_{\beta}$ is of the COE type
and the average level spacing ratio acquires the value $\braket{r}_{\rm
COE}\simeq 0.5269$.

A level spacing distribution of the Poisson type, instead, corresponds to an
integrable dynamics [56] and is observed for instance for many-body-localized
systems, both autonomous [35] and periodically-driven [57]. Like the present
ones, these are systems with disorder, and they never thermalize due to the
existence of a superextensive number of localized integrals of motion [58, 59,
60] that deeply affect dynamics, and forbid thermalization of local
observables. In the periodically-driven case this hinders energy absorption
and avoids infinite-temperature thermalization. Integrability and Poisson
level-spacing distribution correspond to an average level spacing ratio
$\braket{r}_{P}\simeq 0.386$, and due to the presence of disorder we expect to
see them here in some range of parameters.

To inquire this point, we plot $\braket{r}$ versus $J$ for different system
sizes in the spin-1/2 model [Fig. 8(a)] and the spin-1 model [Fig. 8(b)]. When
the size $L$ is large enough, we see that there is an interval where
$\braket{r}$ is stuck to the Poisson value. So there is a parameter region –
in all the considered cases $J\lesssim 10^{-1}$ – where the dynamics is
regular. This property persists as the system size is increased (at least for
the system sizes we can numerically attain). Then we can talk about an
integrable-like phase where there is no thermalization and a nontrivial
dynamics like the one leading to persistent period-doubling oscillations can
appear.

\begin{overpic}[width=227.62204pt]{Gap_Ratio_1-2_0-5}\put(0.0,48.0){(a)}\end{overpic}
---
\begin{overpic}[width=227.62204pt]{Gap_Ratio_1_0-5}\put(0.0,55.0){(b)}\end{overpic}
Figure 8: Average level spacing ratio $\braket{r}$ versus $J$ in the spin-1/2
case (panel a) and spin-1 case (panel b) for different $L$ and $\alpha=0.5$.

##  V Properties of the Floquet spectrum

In order to understand the behavior of the model throughout the entire Hilbert
space, it is useful to look at the Floquet spectrum. Thus, we discuss the
properties of the quasienergies $\\{\mu_{\beta}\\}$ (Sec. V.1) and the Floquet
states $\\{\ket{\psi_{\beta}}\\}$ (Sec. V.2) in relation with the period
doubling.

###  V.1 $\pi$-spectral pairing of the quasienergies

####  V.1.1 Theoretical discussion

A necessary condition for the presence of period doubling is the so called
$\pi$-spectral pairing [40, 61, 62, 63, 64], i. e., each quasienergy must have
a partner shifted by $\pi$. (See [22] and Appendix B for an analysis of how
$\pi$-spectral pairing gives rise to period doubling). At the solvable point
our spin-1/2 model shows $\pi$-spectral pairing, as can be easily checked.

Let us define $\ket{\\{s_{\overline{k}}\\}}$ as the transformed configuration
of $\ket{\\{s_{k}\\}}$ under spin swap, and consider that the energy of the
classical configuration is

$E(\\{s_{k}\\})=\sum_{k=1}^{L-1}\sum_{q=k+1}^{L}V_{kq}s_{k}^{z}s_{q}^{z}+\sum_{k=1}^{L}h_{k}^{z}s_{k}^{z}\,.$

We can write the action of the Floquet operator Eq. (3) on the spin
configurations mentioned above

$\displaystyle\hat{U}_{F}\ket{\\{s_{k}\\}}=e^{-\boldsymbol{i}E(\\{s_{k}\\})}\ket{\\{s_{\overline{k}}\\}}\,,$
$\displaystyle\hat{U}_{F}\ket{\\{s_{\overline{k}}\\}}=e^{-\boldsymbol{i}E(\\{s_{\overline{k}}\\})}\ket{\\{s_{k}\\}}\,.$
(14)

From these relations one can diagonalize $\hat{U}_{F}$ and find the Floquet
states as

$\ket{\psi_{\beta}^{\pm}}=\frac{1}{\sqrt{2}}[e^{\boldsymbol{i}E(\\{s_{k}\\})/2}\ket{\\{s_{k}\\}}\pm
e^{\boldsymbol{i}E(\\{s_{\overline{k}}\\})/2}\ket{\\{s_{\overline{k}}\\}}]\,,$
(15)

and the corresponding quasienergies respectively as

$\displaystyle\mu_{\beta}^{+}=[E(\\{s_{k}\\})+E(\\{s_{\overline{k}}\\})]/2\,,$
$\displaystyle\mu_{\beta}^{-}=\mu_{\beta}^{+}+\pi\,,\quad{\rm
whenever}\quad\\{s_{\overline{k}}\\}\neq\\{s_{k}\\}\,.$ (16)

Therefore a part of the Floquet states show $\pi$-spectral pairing, i. e.,
$\mu_{\beta}^{-}-\mu_{\beta}^{+}=\pi$. This is valid only for the
configurations such that $\\{s_{\overline{k}}\\}\neq\\{s_{k}\\}$, that are
$(2s+1)^{L}-(2s+1)^{L/2}$ in number (see Sec. II.1.2), that is, the
overwhelming majority. If the configurations $\\{s_{k}\\}$ and
$\\{s_{\overline{k}}\\}$ differ over an extensive number of sites, the Floquet
states $\ket{\psi_{\beta}^{\pm}}$ are cat states and are globally different
from each other.

We can interpret the oscillations of the magnetization as Rabi oscillations of
a spin configuration $\ket{\\{s_{k}\\}}$ with its swapped partner
$\ket{\\{s_{\overline{k}}\\}}$. Indeed these configurations are given by
linear combinations of the two Floquet states
$\ket{\psi_{\beta,+}}\pm\ket{\psi_{\beta,-}}$, and the relative phase $\pm
1=\mp e^{\boldsymbol{i}\pi}$ gets inverted every time the Floquet operator
acts on the state. Specifically, if all eigenstates are $\pi$-paired, this
guarantees that there will be period-doubled oscillations for any generic
initial state.

We note here that the presence of disorder and power-law interactions is
fundamental, as it removes all degeneracies between distinct quasienergies at
the solvable point. Due to this gap, the application of a small perturbation
is expected to preserve the structure of the Floquet states and the spectral
pairing.

More precisely, states that differ only locally (i. e are coupled only by
operators localized on a specific small range of sites [40]) will have
different quasienergies. This property is called local spectral gap, and in
the presence of this property we expect that for all eigenstates local
perturbations perturb locally, in a way that is analogous to the ground states
of gapped systems [65]: These perturbations couple globally different Floquet
states at order $\mathcal{O}(L)$ in perturbation theory, and so do not mix
them at any perturbative order in the limit $L\to\infty$ [66, 67]. This
implies that there exists a local unitary transformation which relates the
original eigenstates to the perturbed eigenstates [68, 69] and in particular
the properties of the Floquet spectrum are conserved under local perturbations
[40, 66]. As a consequence, we expect that even in the presence of
perturbations ($J,\varepsilon\neq 0$) there will still be $\pi$-spectral
pairing and thus a time crystal phase over a finite region of parameters. This
expectation is fully confirmed by the numerical results we are going to
discuss.

####  V.1.2 Numerical analysis

To study the persistence of the spectral pairing outside the solvable case, we
consider the following quantities [70, 21]

$\begin{split}\Delta_{\beta}^{0}&=\mu_{\beta+1}-\mu_{\beta}\,;\\\
\Delta_{\beta}^{\pi}&=\min_{\gamma}|\mu_{\gamma}-(\mu_{\beta}+\pi)_{1}|\,.\end{split}$
(17)

All quasienergies are taken in the first Brillouin zone
$\mu_{\beta},(\mu_{\beta}+\pi)_{1}\in(-\pi,\pi)$ and the differences between
quasi-energies are also taken $\text{mod }2\pi$, i. e., we evaluate
$\Delta_{\beta}^{\pi}$ by minimizing the $\pi$-distance between quasienergies
on the unit circle, as shown in Fig. 9.

In case the system features $\pi$-spectral pairing in the thermodynamic limit,
the $\pi$-shifted spectral gaps $\Delta_{\beta}^{\pi}$ should scale to zero
with system size faster than the spectral gaps $\Delta_{\beta}^{0}$. In order
to quantitatively probe this property, we numerically diagonalize
$\hat{U}_{F}$ in the $\mathcal{S}_{z}=0$ subspace and compute the logarithmic
averages [70, 71]

$\displaystyle\braket{\log_{10}\Delta^{\pi}}\equiv\frac{1}{\mathcal{N}}\sum_{\beta}\overline{\log_{10}\Delta_{\beta}^{\pi}}\,,\quad{\rm
and}$
$\displaystyle\braket{\log_{10}\Delta^{0}}\equiv\frac{1}{\mathcal{N}}\sum_{\beta}\overline{\log_{10}\Delta_{\beta}^{0}}\,,$
(18)

where the sum over $\beta$ is restricted to the quasienergies of Floquet
states in the sector of the Hilbert space with $\mathcal{S}_{z}=0$.

Figure 9: Choice of $\pi$-paired eigenvalue $\mu_{\gamma^{*}}$ corresponding
to a given $\mu_{\beta}$ and its $\pi$-shifted spectral gap
$\Delta_{\beta}^{\pi}$.

In order to understand if there is $\pi$ spectral pairing in the thermodynamic
limit, we consider the spectral-pairing parameter

$\ell_{\Delta}\equiv\braket{\log\Delta^{\pi}}-\braket{\log\Delta^{0}}\,.$ (19)

In the case of spectral pairing this difference should decrease with system
size $L$. We plot $\ell$ versus $J$ for different values of $L$ in Fig. 10(a)
for the spin-1/2 and Fig. 10(b) for the spin-1 case.

\begin{overpic}[width=227.62204pt]{Log_Delta_1-2_0-5}\put(0.0,48.0){(a)}\end{overpic}
---
\begin{overpic}[width=227.62204pt]{Log_Delta_1_0-5}\put(0.0,50.0){(b)}\end{overpic}
Figure 10: Spectral-paring parameter $\ell_{\Delta}$ versus $J$ in the
spin-1/2 case (panel a) and spin-1 case (panel b) for different $L$ and
$\alpha=0.5$.

We see that in all the considered cases there is a threshold value for $J\sim
10^{-1}$ where all the curves appear to cross. Below this threshold,
$\ell_{\Delta}$ decreases with the system size, marking the presence of
$\pi$-spectral pairing (the results at large $L$ and very small $J$ are flawed
by the finite representation of floating-point numbers in our numerics). In
contrast, above this threshold $\ell_{\Delta}$ no longer decreases with $L$
(either saturates or gets closer to $0$), and so there is no $\pi$-spectral
pairing.

In order to further understand the stability of the time crystal phase, we
also measure the spectral pairing for different non-zero values of $J$ and
$\varepsilon$ in the spin-1/2 case, as seen in Fig. 11. We find again that
even with both $J$ and $\varepsilon$ non-zero $\ell_{\Delta}$ decreases
exponentially in system size, albeit this is possible for fairly small values
of the perturbations.

So, we have $\pi$-spectral pairing not only at the solvable point $J=0$, but
also in a full range around this value. Therefore, there is a phase where
time-translation symmetry breaking appears, as we theoretically expected.
Quite nicely, the range of parameters of $J$ – at $\varepsilon=0$ – where
$\pi$-spectral pairing appears coincides with the range where the system has a
Poisson level statistics (see Fig. 8), and so the dynamics is non-
thermalizing, and persistent period-doubling oscillations can appear.

Figure 11: Spectral-paring parameter $\ell_{\Delta}$ versus $L$ in the
spin-1/2 case for different values of $J$, $\varepsilon$ and $\alpha=0.5$.
Numerical parameters: $N_{d}$ varies between $\sim 10^{2}$ and $\sim 10^{4}$.

###  V.2 Correlations in the Floquet states

Let us discuss the structure of the Floquet states. As we have seen above, the
majority of these states are the superposition of two different classical spin
configurations (cat states) similar to Eq. (15). The same occurs in the case
of the usual spin-flip time crystal [45, 40]. In that case all the Floquet
states are cat states, and all them show long-range correlations. Defining
$C_{\beta}(i,j)\equiv\braket{\psi_{\beta}}{\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}}{\psi_{\beta}}-\braket{\psi_{\beta}}{\hat{\sigma}_{i}^{z}}{\psi_{\beta}}\braket{\psi_{\beta}}{\hat{\sigma}_{j}^{z}}{\psi_{\beta}}$,
one has $C_{\beta}(1,L)\neq 0$, whatever the state $\ket{\psi_{\beta}}$ and
the length $L$ chosen.

In our case the situation is different. Here there can be spin pairs left
unchanged by the swap, and the cat-state features might not be fully captured
by looking at the correlations between first and last site. For instance,
consider the state which is ferromagnetic in the first half of the chain and
antiferromagnetic in the second half, and construct the Floquet states
according to Eq. (15)
$\ket{\psi_{\beta}}\propto\ket{++,++,\cdots,+-,+-}+e^{\boldsymbol{i}\phi}\ket{\text{swap}}$.
If we measure correlations only in the first half, or between the first and
second half, we would incorrectly classify this as a short-range correlated
state, whereas within the second half of the chain there are indeed long-range
correlations. Therefore, to properly measure correlations in such state it is
necessary to take into account a global correlation quantifier. We use the sum
of all two-point connected correlation function, divided by the number of
spins

$\Sigma_{\beta}\equiv\frac{1}{L}\sum_{i<j}|C_{\beta}(i,j)|\,,$ (20)

with the two-site correlation $C_{\beta}(i,j)$ defined above. Such a quantity
measures the amount of correlations per site. Heuristically, we expect two
distinct behaviors for short-range and long-range correlated states: for
short-range “physical” states, a given site will be correlated only with a
small (non-extensive) amount of neighbors, meaning $\Sigma_{\beta}\sim O(1)$;
for long-range “unphysical” states, a given site will be potentially
correlated with a finite fraction of the entire chain, meaning that
$\Sigma_{\beta}\sim O(L)$.

Let us first focus on the solvable case. Looking at the definition of the
local imbalance Eq. (10) we see that it is kept unchanged under swap, so it
has a well-defined value over the Floquet states (that are superposition of
two classical configuration related by the swapping operation). With a bit of
work (see Appendix E.1) we can find for $s=1/2$ the simple relation

$\Sigma_{\beta}={\mathcal{I}_{\rm LI}}_{,\beta}({\mathcal{I}_{\rm
LI}}_{,\beta}\cdot L-1)/2\,.$ (21)

From it one can see that the sum of correlations increases (quadratically)
with $L$ if and only if the local imbalance is nonvanishing. This shows why
the local imbalance is a natural measure of the time-translation symmetry
breaking properties of a Floquet state: It quantifies the amount of long-range
correlations related to its cat-state nature. In particular, notice that when
$\mathcal{I}_{\rm LI}=0$, the eigenstates are trivial spin configurations
$\ket{\psi_{\beta}}=\ket{\\{s_{k}\\}}$, where there are no long-range
correlations and no spectral pairing. Eq. (21) is valid only for $s=1/2$. For
generic $s$ a similar relation is valid if correlations are evaluated using
the mutual information (see Appendix E.2) instead of the connected correlator.

Let us move to the numerical results in the presence of perturbations.
Specifically, we look at the average $\braket{\Sigma}$, defined as

$\braket{\Sigma}=\frac{1}{\mathcal{N}}\overline{\sum_{\beta}\Sigma_{\beta}}\,,$

with $\beta$ restricted the $\mathcal{S}_{z}=0$ subspace. We plot
$\braket{\Sigma}$ versus $L$ in Fig. 12(a). We see two distinct behaviors of
$\Sigma_{\beta}$. At small $J$, close to the solvable point, it increases
linearly with $L$, as a signature of long-range correlations, as we have
mentioned above. The range of $J$ where this fact occurs is in agreement with
the range where spectral pairing occurs [see Fig. 10(a)] and the dynamics is
regular [see Fig. 8(a)]. At large $J$ instead $\braket{\Sigma}$ gets smaller
and possibly saturates with the system size, and accordingly there is no
spectral pairing and no regular dynamics.

\begin{overpic}[width=227.62204pt]{CORR_1-2_0-5}\put(0.0,48.0){(a)}\end{overpic}
---
\begin{overpic}[width=227.62204pt]{CORR_1_0-5}\put(0.0,65.0){(b)}\end{overpic}
Figure 12: $\braket{\Sigma}$ vs $L$ for different values of $J$ in the
spin-1/2 (panel a) and spin-1 (panel b) cases. Numerical parameters:
$\alpha=0.5$, $N_{d}=5120\cdot 2^{1-L/2}$.

A similar analysis can be performed in the spin-1 case replacing
$\hat{\sigma}_{i}^{z}$ with $\hat{S}_{i}^{z}$ in the definition of
$C_{\beta}(i,j)$ and computing the resulting $\Sigma_{\beta}$ over the Floquet
states. At the solvable point no straightforward relation may be found between
this quantity and local imbalance. Nevertheless, when we apply the
perturbation $J\neq 0$, we still find that $\braket{\Sigma}$ linearly
increases with $L$ in a range of $J$ [see Fig. 12(b)], marking long-range
correlations. In this range there is also spectral pairing [see Fig. 10(b)]
and regular dynamics [see Fig. 8(b)]. For larger $J$ there is none of these
properties and no time crystal. [72]

##  VI Conclusions

In conclusion we have described a new mechanism to realize time crystals in
Floquet systems, based on a spin Hamiltonian with disorder and long-range
interactions, to which a time-periodic swapping of nearby sites is applied. In
contrast with the existing examples, based on the flipping of the local spins,
this protocol is not restricted to spin-1/2 Hamiltonian, but can be applied to
models with local spins of any magnitude.

After that, we have introduced the model Hamiltonians, and have analytically
discussed the case of the solvable point, where the period-doubling
oscillations of the local magnetizations persist forever, also at finite
system size. If one takes as initial state a classical spin configuration, one
can see period-doubling oscillations for the large majority of these
configurations, with the exception of a fraction exponentially decaying to
zero with system size. In order to characterize which configurations give rise
to a period-doubling dynamics, we have introduced a quantity called local
imbalance ($\mathcal{I}_{\rm LI}$). We have seen that, for increasing local
spin magnitude $s$, the spin configurations have in average a larger
$\mathcal{I}_{\rm LI}$, with a smaller variance, and that implies an increased
tendency to provide a period-doubling behavior.

Then we have moved away from the solvable point, to see if the period-doubling
oscillations are robust to such a perturbation. This is an important point to
understand if we are dealing with a time crystal, i. e. a robust phase of the
matter stabilized by the interactions over a finite parameter range. For that
purpose we have numerically studied the model with $s=1/2$ (spin-1/2 case) and
the one with $s=1$ (spin-1 case). For perturbations small enough but finite we
see that the period-doubling oscillations persist for a time that
exponentially increases with the system size: Only in the large-size limit
they persist forever. In this way the time-translation symmetry breaking
occurs only in the thermodynamic limit, in analogy with the standard symmetry
breaking of quantum phase transitions.

We have then focused on the properties of Floquet quasienergies and states. We
have considered a probe of quantum regularity/ergodicity, the average level
spacing ratio of Floquet quasienergies, to understand in which range of
parameters the quantum dynamics is regular, thermalization of local
observables is absent, and then period-doubling oscillations are possible. We
have found a range of parameters where this ergodicity breaking occurs, and
the effect becomes more marked when the system size is increased, so that we
deal with a regular-dynamics phase. In this range of parameters a nontrivial
period-doubling dynamics is possible, and the ergodicity breaking is an effect
of disorder that induces many-body localization.

We have moved to period-doubling related properties of Floquet quasienergies.
We have analytically shown that at the solvable point, up to a fraction
exponentially vanishing with the system size, quasienergies show
$\pi$-spectral pairing, a property strictly associated to period doubling.
Moving away from the solvable point, we have numerically shown that the
$\pi$-spectral pairing is robust in a well-defined range of parameters, that
defines the time-crystal phase. This robustness is related to disorder and
long-range interactions, that break degeneracies of the quasienergies at the
solvable point, making properties of quasienergies and Floquet states robust
under local perturbations, and the range of parameters where it occurs lies
inside the range where the dynamics is regular and there is no thermalization.

Looking at the properties of the Floquet states, we have defined a global
correlation quantifier. At the solvable point we have found that it is
strictly related to the $\mathcal{I}_{\rm LI}$ of the Floquet state: a finite
$\mathcal{I}_{\rm LI}$ corresponds to a Floquet state with the cat-state
structure related to $\pi$-spectral pairing. Moving away from the solvable
point, we have found that there is a finite range of parameters where the
correlations increase linearly with system size, which indicates cat-state
behaviour and symmetry-breaking.

The ranges of parameters where there is regular dynamics, $\pi$-spectral
pairing and long-range correlations agree with each other. Together with the
results about the dynamics, this suggests the existence of a time-crystal
phase.

Perspectives of future work include the application of this protocol to models
with different local Hilbert space (for instance bosons or rotors), and the
use of different physical effects to induce ergodicity-breaking (for instance
the Stark many-body localization [73, 74]).

### Acknowledgements

We thank M. Fava and R. Fazio for insightful comments on the manuscript. We
acknowledge the Peter-Grünberg Institute, Quantum Control (PGI-8),
Forschungszentrum Jülich for providing computational resources that
facilitated the completion of this work. G. P. and A. R. acknowledge financial
support from PNRR MUR project PE0000023-NQSTI. P.L. acknowledges computational
resources from MUR, PON “Ricerca e Innovazione 2014-2020”, under Grant No.
PIR01_00011 - (I.Bi.S.Co.). P.L. acknowledges financial support from PNRR MUR
project CN_00000013-ICSC, as well as from the project QuantERA II Programme
STAQS project that has received funding from the European Union’s Horizon 2020
research and innovation programme under Grant Agreement No 101017733. R. G.
acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – Cluster of
Excellence Matter and Light for Quantum Computing (ML4Q) EXC 2004/1 –
390534769.

##  Appendix A Spin swapping

Here we show that the operator in Eq. (7) realizes swaps product states. It
suffices to show this in the $\hat{\sigma}_{z}$ basis $\ket{s_{i},s_{j}}$, and
then extends the result linearly to all product states. So we wish to show

$\frac{1+\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}{2}\ket{s_{i},s_{j}}=\ket{s_{j},s_{i}}.$
(22)

We use the defining identities of the Pauli operators
$\hat{\sigma}_{z}\ket{s}=s\ket{s}$, $\hat{\sigma}_{x}\ket{s}=\ket{-s}$ and
$\hat{\sigma}_{y}=\boldsymbol{i}\hat{\sigma}_{x}\hat{\sigma}_{z}$, hence:

$\begin{split}\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}&=\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}+\hat{\sigma}_{i}^{x}\hat{\sigma}_{j}^{x}+\hat{\sigma}_{i}^{y}\hat{\sigma}_{j}^{y}=\\\
&=\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}+\hat{\sigma}_{i}^{x}\hat{\sigma}_{j}^{x}(1-\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}),\\\
\frac{1+\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}{2}&=\frac{1+\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}}{2}+\hat{\sigma}_{i}^{x}\hat{\sigma}_{j}^{x}\frac{1-\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}}{2}.\end{split}$
(23)

Since the $\hat{\sigma}_{k}^{z}$ operators act diagonally on the spin basis,
they serve purely as a filter to either the condition
$\delta_{s_{i},s_{j}}=\frac{1+s_{i}s_{j}}{2}$ (first term, aligned spins) or
the condition
$\delta_{s_{i},-s_{j}}=1-\delta_{s_{i},s_{j}}=\frac{1-s_{i}s_{j}}{2}$ (second
term, unaligned spins):

$\frac{1+\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}{2}\ket{s_{i},s_{j}}=\begin{cases}\ket{s_{i},s_{j}}&s_{i}=s_{j},\\\
\ket{-s_{i},-s_{j}}&s_{i}=-s_{j},\end{cases}=\ket{s_{j},s_{i}},$ (24)

as expected. From this relation one can also see how in the spin-1 case the
swap operator has to involve higher power of the spin operators, since already
the ”filter” itself $\delta_{s_{i},s_{j}}$ (and its complement) involve terms
such as $s_{i}^{2}+s_{j}^{2}$, as seen in Eq. (III).

We now aim to demonstrate Eq. (6). This is equivalent to demonstrate the
following equality

$e^{-\boldsymbol{i}(\pi/4)\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}=e^{\boldsymbol{i}\phi}\frac{1+\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}{2}\,,$
(25)

with $i$ and $j$ arbitrary sites and some phase $\phi$. To do this, we use the
fact that the swap operator is an involution and thus squares to the identity:

$\left(\frac{1+\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}{2}\right)^{2}=1.$
(26)

This allows us to expand the exponential in a power series, which corresponds
to Euler’s identity

$e^{-\boldsymbol{i}\frac{\pi}{2}\frac{1+\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}{2}}=-\boldsymbol{i}\frac{1+\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}{2}$
(27)

and leads to the result we were looking for

$e^{-\boldsymbol{i}\frac{\pi}{4}\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}=e^{-\boldsymbol{i}\pi/4}\frac{1+\vec{\boldsymbol{\sigma}}_{i}\cdot\vec{\boldsymbol{\sigma}}_{j}}{2}.$
(28)

A similar derivation can also be done for the general spin-$s$ case, by
starting with the swap operator $\hat{P}$, which satisfies $\hat{P}^{2}=1$:

$e^{-\boldsymbol{i}\frac{\pi}{2}(\hat{P}-1)}=\hat{P}.$ (29)

In the case of spin-1 we arrive at the desired result:

$e^{-\boldsymbol{i}\frac{\pi}{2}((\vec{\mathbf{S}}_{i}\cdot\vec{\mathbf{S}}_{j})^{2}+\vec{\mathbf{S}}_{i}\cdot\vec{\mathbf{S}}_{j})}=-((\vec{\mathbf{S}}_{i}\cdot\vec{\mathbf{S}}_{j})^{2}+\vec{\mathbf{S}}_{i}\cdot\vec{\mathbf{S}}_{j}-1).$
(30)

##  Appendix B $\pi$ spectral pairing and period doubling

Let us show the relation between $\pi$-spectral pairing and period doubling
(see also [49, 22]). Let us assume that at least a part of the Floquet
spectrum is $\pi$-spectral paired, that is to say, it is organized in pairs
$\mu_{\beta}^{-}$, $\mu_{\beta}^{+}$ such that
$\mu_{\beta}^{-}-\mu_{\beta}^{+}=\pi$. We initialize the system in the state
$\ket{\psi(0)}$ and study the stroboscopic dynamics of some observable
$\mathcal{O}(t)=\braket{\psi(t)}{\hat{\mathcal{O}}}{\psi(t)}$, with $t$ an
integer multiple of the driving period $T=1$. We can expand the stroboscopic-
time-evolved state in the Floquet basis as
$\ket{\psi(t)}=\sum_{\beta}R_{\beta}e^{-i\mu_{\beta}t}\ket{\psi_{\beta}}$,
where $R_{\beta}\equiv\braket{\psi_{\beta}}{\psi(0)}$. Expanding the
expression for the expectation of the observable we get

$\displaystyle\mathcal{O}(t)$
$\displaystyle=\underbrace{\sum_{\beta}|R_{\beta}|^{2}\braket{\psi_{\beta}}{\hat{\mathcal{O}}}{\psi_{\beta}}}_{\text{diagonal}}$
$\displaystyle+\underbrace{\sum_{\beta}(R_{\beta}^{+})^{*}R_{\beta}^{-}\braket{\psi_{\beta}^{-}}{\hat{\mathcal{O}}}{\psi_{\beta}^{+}}e^{i(\mu_{\beta}^{-}-\mu_{\beta}^{+})t}}_{\text{period-
doubling}}$
$\displaystyle+\underbrace{\sum_{\beta\neq\gamma,\,\nu\neq\rho}(R_{\beta}^{\nu})^{*}(R_{\gamma}^{\rho})\braket{\psi_{\beta}^{\nu}}{\hat{\mathcal{O}}}{\psi_{\gamma}^{\rho}}e^{i(\mu_{\beta}^{\nu}-\mu_{\gamma}^{\rho})t}}_{\text{off-
diagonal}}\,.$ (31)

The diagonal term provides a constant contribution, while the off-diagonal one
(where $\mu_{\beta}^{\nu}$ and $\mu_{\gamma}^{\rho}$ are different and
unpaired) gives a superposition of incoherent oscillations that, after
disorder-average, vanish after a short transient. The interesting term is the
one denoted “period-doubling”, where one can factorize a term $e^{i\pi t}$
that provides the period-doubling oscillations, due to the fact that
$\mu_{\beta}^{-}-\mu_{\beta}^{+}=\pi$.

We have seen in Sec. V.1.2 that when we are not at the solvable point but $J$
is small enough, the $\pi$-spectral pairing appears only in the thermodynamic
limit. At any finite size one has
$\mu_{\beta}^{-}-\mu_{\beta}^{+}=\pi+\delta_{\beta}$, where $\delta_{\beta}$
is a small correction decreasing with the system size. Due to this correction,
the period-doubling term averaged over disorder provides period-doubling
oscillations that last for a finite time, provided roughly by the inverse of
the typical value of $\delta_{\beta}$ (call it $\delta$). Because $\delta$
decreases with the system size (see Sec. V.1.2), the time of the period-
doubling oscillations increases exponentially with it, as we have seen in
Secs. IV.1 and IV.2.

We emphasize that the period-doubling term is nonvanishing only for an
appropriate choice of operator and initial state. For instance, in Fig. 5, the
period-doubling term vanishes for the onsite magnetizations at the sites
corresponding to the half of the chain where the initial state is uniform.

##  Appendix C Short-range results ($\alpha=3$)

Here we show the results for the short-range case, $\alpha=3$. As mentioned in
the main body of the paper, we only found minor quantitative differences
between the two. In this section it is taken for granted that $\alpha=3$ in
every figure. Let us briefly list the figures:

  * •

In Fig. 13(a) we show $(-1)^{t}\overline{Z}(t)$ versus $t$ for different
system sizes and initial states in the spin-1/2 model. This is the equivalent
of Fig. 3.

  * •

In Fig. 14(a) we show the decay time $\overline{\tau}$ versus $L$ for
different parameters in the spin-1/2 model. This is the equivalent of Fig. 7.
In Fig. 14(b) we show the same the spin-1 model (equivalent of Fig. 7).

  * •

In Fig. 15 we show $\braket{r}$ versus $J$ for both models and different
system sizes. This is the equivalent of Fig. 8.

  * •

In Fig. 16 we show the spectral-pairing parameter defined in Eq. (19) versus
$J$ for both models and different $L$. This is the equivalent of Fig. 10.

  * •

In Fig. 17 we show $\braket{\Sigma}$ vs $L$ for different values of $J$ and
both models. This is the equivalent of Fig. 12.

Comparing Figs. 15, 16, and 17 we see that also for $\alpha=3$ regular
dynamics, spectral pairing and long.range correlations occur in the same range
of $J$.

\begin{overpic}[width=227.62204pt]{Z_avg_1-2_3}\put(0.0,54.0){(b)}\put(0.0,92.0){(a)}\end{overpic}
Figure 13: Spin-1/2 model. (Panel a) $(-1)^{t}\overline{Z}(t)$ versus $t$ for
initial Néel state and different system sizes. (Panel b)
$(-1)^{t}\overline{Z}(t)$ versus $t$ for initial Half-Néel state (solid line)
and Néel state (dashed line) for different system sizes. Initial values have
all been normalized to $+1$ to make the comparison clearer. Parameters:
$J=0.1$, $\varepsilon=0.01$.

\begin{overpic}[width=227.62204pt]{Decay_Times_1-2_3}\put(0.0,52.0){(a)}\end{overpic}

\begin{overpic}[width=227.62204pt]{Decay_Times_1-2_3}\put(0.0,52.0){(b)}\end{overpic}

Figure 14: (Panel a) Spin-1/2 model. Decay time $\overline{\tau}$ versus $L$
for initial Néel state and different values of $J$. Here we have considered
time evolutions while skipping over $n\geq 1$ periods, such that the time step
stays below $\sim 1\%$ of the standard deviation of $\tau$. (Panel b) Spin-1
model. Decay time $\overline{\tau}$ versus $L$ for initial Up-Zero state.
\begin{overpic}[width=227.62204pt]{Gap_Ratio_1-2_3}\put(0.0,52.0){(a)}\end{overpic}
---
\begin{overpic}[width=227.62204pt]{Gap_Ratio_1_3.pdf}\put(0.0,52.0){(b)}\end{overpic}
Figure 15: Average level spacing ratio $\braket{r}$ versus $J$ in the spin-1/2
case (panel a) and spin-1 case (panel b) for different $L$ and $\alpha=3$.
\begin{overpic}[width=227.62204pt]{Log_Delta_1-2_3}\put(0.0,52.0){(a)}\end{overpic}
---
\begin{overpic}[width=227.62204pt]{Log_Delta_1_3}\put(0.0,52.0){(b)}\end{overpic}
Figure 16: Spectral-pairing parameter $\ell_{\Delta}$ versus $J$ in the
spin-1/2 case (panel a) and spin-1 case (panel b) for different $L$ and
$\alpha=3$.
\begin{overpic}[width=227.62204pt]{CORR_1-2_3}\put(0.0,52.0){(a)}\end{overpic}
---
\begin{overpic}[width=227.62204pt]{CORR_1_3}\put(0.0,52.0){(b)}\end{overpic}
Figure 17: $\braket{\Sigma}$ vs $L$ for different values of $J$ in the
spin-1/2 (panel a) and spin-1 (panel b) case. Numerical parameters:
$\alpha=3$, $N_{d}=5120\cdot 2^{1-L/2}$.

##  Appendix D Local imbalance for spin states

As discussed in section II.2, the spin-1 case displays many differences with
the spin-1/2 case. From the point of view of the spectrum, at the solvable
point, we can already notice some variations. Apart from the fact that a
larger fraction of the spectrum participates to period-doubling at any $L$, a
generic spin configuration will also tend to have a larger local imbalance.
Quantitatively, this may be seen by looking at the asymptotic behavior
$L\to\infty$ of the dimensions of the eigenspaces at fixed local imbalance
(see Figure 18), normalized in such a way that the area under the graph is
$1$.

This may be seen more generally in the case of a chain of length $L$ where at
each site there corresponds a local Hilbert space of dimension $d=2s+1$, where
$s$ is the value of the local spin. It is easy to see that the
$\mathcal{I}_{\rm LI}$ can only take values

$\mathcal{I}_{\rm LI}=\frac{N}{L/2}\,,$

where $N\leq L/2$ is the number of pairs whose spins differ from each other.
There are many spin configurations associated to a single value of
$\mathcal{I}_{\rm LI}$. This degeneracy, i. e. the dimension
$D_{L}(\mathcal{I}_{\rm LI})$ of an eigenspace corresponding to a given value
of local imbalance $\mathcal{I}_{\rm LI}$ is

$D_{L}\left(\mathcal{I}_{\rm
LI}=\frac{N}{L/2}\right)=(d(d-1))^{N}d^{L/2-N}\binom{L/2}{N}\,.$ (32)

One can find this result considering that each pair where spins differ can
appear in $d(d-1)$ different ways, while each pair where spins are equal can
appear in $d$ different ways. So the probability distribution of the
imbalances

$p(\mathcal{I}_{\rm LI})=\frac{1}{d^{L}}D_{L}\left(\mathcal{I}_{\rm
LI}=\frac{N}{L/2}\right)$

is a binomial. Therefore, the asymptotic distribution $p(\mathcal{I}_{\rm
LI})$ is a normal distribution with mean

$\braket{\mathcal{I}_{\rm LI}}=\frac{(d-1)}{d}$ (33)

and variance

$\sigma^{2}(\mathcal{I}_{\rm LI})=p(1-p)/(L/2)=\frac{2(d-1)}{Ld^{2}}$ (34)

vanishing for $L\to\infty$ and $s\to\infty$. This proves more concretely that
increasing $s$ improves the behavior of the period-doubling caused by the
swap, because the typical value of the local imbalance increases and the
fluctuations around it decrease. Although no such direct analytical result can
be found restricting to the subspace $\mathcal{S}_{z}=0$, the same qualitative
behavior is observed numerically at large $L$ (not shown). This suggests once
again that the subspace $\mathcal{S}_{z}=0$ we numerically consider here
faithfully represents the entire Hilbert space.

Figure 18: Probability Density Function of the Local Imbalance at different
local Hilbert space dimensions $d$ and $L=120$ (solid line) and its
approximation as a normal distribution (dashed line). Spin-1/2 and spin-1
correspond to $d=2$ and $d=3$, respectively.

##  Appendix E Correlations in the exact Floquet states

In this Appendix we find the relation between the local imbalance and the
global correlation quantifier Eq. (20). In Sec. E.1 we consider the case where
the global quantifier is constructed using the connected correlations and the
relation is valid only for local spin magnitude $s=1/2$. In Sec. E.2 we
evaluate correlations using the mutual information and get a relation valid
for generic $s$.

###  E.1 Connected correlations

Here we derive Eq. (21) relating the local imbalance and the global
correlation quantifier Eq. (20). Eq. (21) is valid in the solvable case
($\varepsilon=0=J$), for the spin-1/2 chain. In this case the Floquet states
are given by superpositions of a classical spin configuration and its swapped
partner:

$\ket{\psi_{\beta}}=c_{+}\ket{\\{s_{k}\\}}+c_{-}\ket{\\{s_{\overline{k}}\\}},$
(35)

where $|c_{+}|=|c_{-}|=1/\sqrt{2}$ if $\mathcal{I}_{\rm LI}>0$ while
$c_{+}=1$, $c_{-}=0$ (or equivalently $c_{-}=1$, $c_{+}=0$) if
$\mathcal{I}_{\rm LI}=0$. Since we are dealing with expectation values of
$\hat{\sigma}^{z}$ operators and their products, the phases of $c_{+}$ and
$c_{-}$ are irrelevant and we ignore them for convenience of notation.

We first notice that the spins at sites $k$ where $s_{k}=s_{\overline{k}}$ are
independent from every other spin, since their value remains the same in the
superposition. This means that the state factorizes as a product state between
site $k$ and every other site:

$\begin{split}\ket{\psi_{\beta}}=c_{+}&\ket{s_{1}\cdots
s_{k}\cdots}+c_{-}\ket{s_{\overline{1}}\cdots
s_{k}\cdots}=\ket{\psi_{\beta}}_{\not{k}}\otimes\ket{s_{k}}\,.\end{split}$
(36)

Therefore, any connected correlation function between such sites $k$ and any
other site $q$ will be zero, $C_{\beta}(k,q)=0,\;\forall q\neq k$.

There only remains the case where both $k$ and $q$ satisfy
$s_{k}=-s_{\overline{k}}$, $s_{q}=-s_{\overline{q}}$. The evaluation of
$C_{\beta}(k,q)$ is now trivial if one notices that
$s_{k}s_{q}=s_{\overline{k}}s_{\overline{q}}$:

$\begin{split}\braket{\hat{\sigma}_{k}^{z}\hat{\sigma}_{q}^{z}}_{\beta}&=|c_{+}|^{2}\bra{\\{s_{i}\\}}\hat{\sigma}_{k}^{z}\hat{\sigma}_{q}^{z}\ket{\\{s_{i}\\}}+|c_{-}|^{2}\bra{\\{s_{\overline{i}}\\}}\hat{\sigma}_{k}^{z}\hat{\sigma}_{q}^{z}\ket{\\{s_{\overline{i}}\\}}\\\
&=s_{k}s_{q},\\\ \\\
\braket{\hat{\sigma}_{k}^{z}}_{\beta}&=|c_{+}|^{2}\bra{\\{s_{i}\\}}\hat{\sigma}_{k}^{z}\ket{\\{s_{i}\\}}+|c_{-}|^{2}\bra{\\{s_{\overline{i}}\\}}\hat{\sigma}_{k}^{z}\ket{\\{s_{\overline{i}}\\}}\\\
&=|c_{+}|^{2}(s_{k}-s_{\overline{k}})=0,\\\ \\\
C_{\beta}(k,q)&=\braket{\hat{\sigma}_{k}^{z}\hat{\sigma}_{q}^{z}}_{\beta}-\braket{\hat{\sigma}_{k}^{z}}_{\beta}\braket{\hat{\sigma}_{q}^{z}}_{\beta}=s_{k}s_{q}.\end{split}$
(37)

In the spin-1/2 case this can only take values $\pm 1$, meaning that its
absolute value is always $+1$. Finally, we have that $|C_{\beta}(k,q)|$
vanishes for every site $k$ aligned with its partner $\overline{k}$, while it
equals one for every choice of sites $k$ and $q$ that are both unaligned with
respect to their partners $\overline{k}$ and $\overline{q}$. By definition of
local imbalance, the number of such sites is equal to ${\mathcal{I}_{\rm
LI}}_{,\beta}\cdot L$ (double the number of unaligned pairs). Finally, the
computation of $\Sigma_{\beta}$ reduces to counting the number of distinct
pairs of sites within this region:

$\begin{split}\Sigma_{\beta}=\frac{1}{L}\sum_{i<j}|C_{\beta}(i,j)|=\frac{1}{L}\sum_{i|s_{i}\neq
s_{\overline{i}}}&\,\sum_{j>i|s_{j}\neq s_{\overline{j}}}+1=\\\
=\frac{{\mathcal{I}_{\rm LI}}_{,\beta}({\mathcal{I}_{\rm LI}}_{,\beta}\cdot
L-1)}{2},\end{split}$ (38)

which is the desired result.

###  E.2 Mutual information

A general, model-independent, alternative to connected correlation functions
which can measure correlations between two sites is given by the mutual
information $\mathcal{I}$:

$\mathcal{I}_{\beta}(i,j)=S_{i}^{(\beta)}+S_{j}^{(\beta)}-S_{i,j}^{(\beta)},$
(39)

where $S_{A}^{(\beta)}=-\mathrm{Tr}(\rho_{A}^{(\beta)}\ln\rho_{A}^{(\beta)})$
is the entanglement entropy of the Floquet state $\ket{\psi_{\beta}}$ over the
region $A$ (i.e. a set of sites in the chain). Calling $B$ the complement of
$A$, the reduced density matrix is defined as
$\rho_{A}^{(\beta)}=\mathrm{Tr}_{B}[\ket{\psi_{\beta}}\bra{\psi_{\beta}}]$.
The global correlation quantifier for the mutual information is defined as

$\Sigma_{\beta}(\mathcal{I})=\frac{1}{L}\sum_{i<j}\mathcal{I}_{\beta}(i,j)\,.$
(40)

We now show by using the mutual information instead of a connected correlation
function that we get the same scaling law in the general spin-$s$ case. Let us
first explicitly write the reduced density matrix of sites $i$ and $j$ for the
Floquet state
$\ket{\psi_{\beta}}=c_{+}\ket{\\{s_{k}\\}}+c_{-}\ket{\\{s_{\overline{k}}\\}}$:

$\begin{split}\rho_{i,j}^{(\beta)}&=\mathrm{Tr}_{\not{i}\not{j}}\ket{\psi_{\beta}}\bra{\psi_{\beta}}\\\
&=|c_{+}|^{2}\ket{s_{i},s_{j}}\bra{s_{i},s_{j}}+|c_{-}|^{2}\ket{s_{\overline{i}},s_{\overline{j}}}\bra{s_{\overline{i}},s_{\overline{j}}}+\\\
&+c_{+}c_{-}^{*}\delta_{\\{s_{k}\\}_{\not{i}\not{j}},\\{s_{\overline{k}}\\}_{\not{i}\not{j}}}\ket{s_{i},s_{j}}\bra{s_{\overline{i}},s_{\overline{j}}}+\mathrm{h.c.},\end{split}$
(41)

where each $s_{k}$ now takes $d=2s+1$ different possible values in
$\\{-s,-s+1,\cdots,s-1,s\\}$. If
$\\{s_{k}\\}_{\not{i}\not{j}}\neq\\{s_{\overline{k}}\\}_{\not{i}\not{j}}$ \-
which means that there is at least one site $q$ where $s_{q}\neq
s_{\overline{q}}$ \- then the reduced density matrix takes the simple form:

$\rho_{i,j}^{(\beta)}=|c_{+}|^{2}\ket{s_{i},s_{j}}\bra{s_{i},s_{j}}+|c_{-}|^{2}\ket{s_{\overline{i}},s_{\overline{j}}}\bra{s_{\overline{i}},s_{\overline{j}}}$
(42)

and similarly for the reduced density matrix of site $i$ and site $j$:

$\begin{split}\rho_{i}^{(\beta)}&=|c_{+}|^{2}\ket{s_{i}}\bra{s_{i}}+|c_{-}|^{2}\ket{s_{\overline{i}}}\bra{s_{\overline{i}}},\\\
\rho_{j}^{(\beta)}&=|c_{+}|^{2}\ket{s_{j}}\bra{s_{j}}+|c_{-}|^{2}\ket{s_{\overline{j}}}\bra{s_{\overline{j}}},\end{split}$
(43)

It is straightforward to see that when $|c_{+}|^{2}=|c_{-}|^{2}=1/2$ and
$s_{i}\neq s_{\overline{i}},\,s_{j}\neq s_{\overline{j}}$ the entropy of these
states is $\ln 2$ and consequently the mutual information is also $\ln 2$:

$\mathcal{I}_{\beta}(i,j)=S_{i}+S_{j}-S_{ij}=\ln 2.$ (44)

Outside of an irrelevant special case (at $\mathcal{I}_{\rm
LI}=\frac{1}{L/2}$), the scaling of the global correlation quantifier is still
respected, as only the pairs with unaligned spins will contribute

$\Sigma_{\beta}(\mathcal{I})=\frac{1}{L}\sum_{i<j}\mathcal{I}_{\beta}(i,j)=\frac{{\mathcal{I}_{\rm
LI}}_{,\beta}({\mathcal{I}_{\rm LI}}_{,\beta}\cdot L-1)}{2}\ln 2.$ (45)

We notice that for $L\gg 1$ a typical Floquet state $\ket{\psi_{\beta}}$ has
$(d-1)/d\cdot L/2$ pairs of unaligned neighbouring spins, corresponding to
$\mathcal{I}_{\rm LI}=(d-1)/d$. This typical value coincides with the average
Eq. (33), because the fluctuations Eq. (34) vanish for $L\to\infty$.
Substituting this typical value in Eq. (45) we find that the typical amount of
global correlations increases with local spin magnitude $s$

$\overline{\Sigma(\mathcal{I})}\sim\frac{L}{2}\frac{(d-1)^{2}}{d^{2}}\ln 2\,,$
(46)

and the fluctuations around this value become smaller, as we can see in Eq.
(34). This finding confirms that at higher values of spin the correlations are
stronger and thus the time crystalline behaviour becomes more robust with
increasing $s$.

## References

  * Zhang  _et al._ [2017] C. Zhang, F. Pollmann, S. L. Sondhi, and R. Moessner, Density-matrix renormalization group study of many-body localization in Floquet eigenstates, Annalen der Physik 529, 1600294 (2017).
  * Choi  _et al._ [2017] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Observation of discrete time-crystalline order in a disordered dipolar many-body system, Nature 543, 221 (2017).
  * Wilczek [2012] F. Wilczek, Quantum time crystals, Phys. Rev. Lett. 109, 160401 (2012).
  * Shapere and Wilczek [2012] A. Shapere and F. Wilczek, Classical time crystals, Physical review letters 109, 160402 (2012).
  * Wilczek [2013] F. Wilczek, Superfluidity and space-time translation symmetry breaking, Phys. Rev. Lett. 111, 250402 (2013).
  * Khemani  _et al._ [2019] V. Khemani, R. Moessner, and S. L. Sondhi, A brief history of time crystals (2019), arXiv:1910.10745 [cond-mat.str-el] .
  * Sacha and Zakrzewski [2017] K. Sacha and J. Zakrzewski, Time crystals: a review, Reports on Progress in Physics 81, 016401 (2017).
  * 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 .
  * Zaletel  _et al._ [2023] M. P. Zaletel, M. Lukin, C. Monroe, C. Nayak, F. Wilczek, and N. Y. Yao, Colloquium: Quantum and classical discrete time crystals, Rev. Mod. Phys. 95, 031001 (2023).
  * Sacha [2020] K. Sacha,  _Time crystals_ (Springer, 2020).
  * Watanabe and Oshikawa [2015] H. Watanabe and M. Oshikawa, Absence of quantum time crystals, Phys. Rev. Lett. 114, 251603 (2015).
  * Else  _et al._ [2016a] D. V. Else, B. Bauer, and C. Nayak, Floquet time crystals, Phys. Rev. Lett. 117, 090402 (2016a).
  * Khemani  _et al._ [2016a] V. Khemani, F. Pollmann, and S. L. Sondhi, Obtaining highly excited eigenstates of many-body localized Hamiltonians by the density matrix renormalization group approach, Phys. Rev. Lett. 116, 247204 (2016a).
  * Khemani  _et al._ [2017] V. Khemani, C. W. von Keyserlingk, and S. L. Sondhi, Defining time crystals via representation theory, Phys. Rev. B 96, 115127 (2017).
  * Yao  _et al._ [2017a] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Discrete time crystals: Rigidity, criticality, and realizations, Phys. Rev. Lett. 118, 030401 (2017a).
  * Ho  _et al._ [2017] W. W. Ho, S. Choi, M. D. Lukin, and D. A. Abanin, Critical time crystals in dipolar systems, Phys. Rev. Lett. 119, 010602 (2017).
  * Pizzi  _et al._ [2020] A. Pizzi, D. Malz, G. De Tomasi, J. Knolle, and A. Nunnenkamp, Time crystallinity and finite-size effects in clean floquet systems, Physical Review B 102, 214207 (2020).
  * Smits  _et al._ [2018] J. Smits, L. Liao, H. Stoof, and P. van der Straten, Observation of a space-time crystal in a superfluid quantum gas, Phys. Rev. Lett. 121, 185301 (2018).
  * Pizzi  _et al._ [2019] A. Pizzi, J. Knolle, and A. Nunnenkamp, Period-n discrete time crystals and quasicrystals with ultracold bosons, Phys. Rev. Lett. 123, 150601 (2019).
  * Pizzi  _et al._ [2021a] A. Pizzi, J. Knolle, and A. Nunnenkamp, Higher-order and fractional discrete time crystals in clean long-range interacting systems, Nature communications 12, 1 (2021a).
  * 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).
  * Surace  _et al._ [2019] F. M. Surace, A. Russomanno, M. Dalmonte, A. Silva, R. Fazio, and F. Iemini, Floquet time crystals in clock models, Physical Review B 99, 10.1103/physrevb.99.104303 (2019).
  * Zhu  _et al._ [2019] B. Zhu, J. Marino, N. Y. Yao, M. D. Lukin, and E. A. Demler, Dicke time crystals in driven-dissipative quantum many-body systems, New Journal of Physics 21, 073028 (2019).
  * Lazarides and Moessner [2017] A. Lazarides and R. Moessner, Fate of a discrete time crystal in an open system, Physical Review B 95, 195135 (2017).
  * Pizzi  _et al._ [2021b] A. Pizzi, A. Nunnenkamp, and J. Knolle, Bistability and time crystals in long-ranged directed percolation, Nature communications 12, 1 (2021b).
  * Gong  _et al._ [2018] Z. Gong, R. Hamazaki, and M. Ueda, Discrete time-crystalline order in cavity and circuit qed systems, Phys. Rev. Lett. 120, 040404 (2018).
  * Else  _et al._ [2017] D. V. Else, B. Bauer, and C. Nayak, Prethermal phases of matter protected by time-translation symmetry, Physical Review X 7, 011026 (2017).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Sambe [1973] H. Sambe, Steady states and quasienergies of a quantum-mechanical system in an oscillating field, Physical Review A 7, 2203 (1973).
  * Shirley [1965] J. H. Shirley, Solution of the schrödinger equation with a Hamiltonian periodic in time, Physical Review 138, B979 (1965).
  * Huang [2023] B. Huang, Analytical theory of cat scars with discrete time crystalline dynamics in floquet systems (2023), arXiv:2302.08056 [cond-mat.str-el] .
  * Huang  _et al._ [2022] B. Huang, T.-H. Leung, D. M. Stamper-Kurn, and W. V. Liu, Discrete time crystals enforced by floquet-bloch scars, Phys. Rev. Lett. 129, 133001 (2022).
  * Huang  _et al._ [2018] B. Huang, Y.-H. Wu, and W. V. Liu, Clean floquet time crystals: Models and realizations in cold atoms, Phys. Rev. Lett. 120, 110603 (2018).
  * Sachdev and Keimer [2011] S. Sachdev and B. Keimer, Quantum criticality, arXiv preprint arXiv:1102.4628 (2011).
  * Pal and Huse [2010a] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010a).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Phys. Rev. Lett. 110, 084101 (2013).
  * Kac  _et al._ [1963] M. Kac, G. Uhlenbeck, and P. Hemmer, On the van der Waals Theory of the Vapor-Liquid Equilibrium. I. Discussion of a One-Dimensional Model, J. Math. Phys. 4, 216 (1963).
  * Zhang and Wang [2006] G.-M. Zhang and X. Wang, Spin swapping operator as an entanglement witness for quantum Heisenberg spin-s systems, Journal of Physics A: Mathematical and General 39, 8515 (2006).
  * Serbyn  _et al._ [2014] M. Serbyn, Z. Papić, and D. A. Abanin, Quantum quenches in the many-body localized phase, Phys. Rev. B 90, 174302 (2014).
  * Else  _et al._ [2016b] D. V. Else, B. Bauer, and C. Nayak, Floquet time crystals, Phys. Rev. Lett. 117, 090402 (2016b).
  * Richter  _et al._ [2019] J. Richter, N. Casper, W. Brenig, and R. Steinigeweg, Magnetization dynamics in clean and disordered spin-1 xxz chains, Phys. Rev. B 100, 144423 (2019).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science 349, 842 (2015), https://www.science.org/doi/pdf/10.1126/science.aaa7432 .
  * Kohlert  _et al._ [2019] T. Kohlert, S. Scherg, X. Li, H. P. Lüschen, S. D. Sarma, I. Bloch, and M. Aidelsburger, Observation of many-body localization in a one-dimensional system with a single-particle mobility edge, Physical Review Letters 122, 10.1103/physrevlett.122.170403 (2019).
  * Russomanno  _et al._ [2020] A. Russomanno, S. Notarnicola, F. M. Surace, R. Fazio, M. Dalmonte, and M. Heyl, Homogeneous floquet time crystal protected by gauge invariance, Phys. Rev. Res. 2, 012003 (2020).
  * Yao  _et al._ [2017b] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Discrete time crystals: Rigidity, criticality, and realizations, Phys. Rev. Lett. 118, 030401 (2017b).
  * [46] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalatorre, Nonequilibrium dynamics of closed interacting quantum systems, 2011 rev, Mod. Phys 83, 5.
  * Pal and Huse [2010b] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010b).
  * Tikhonov and Mirlin [2018] K. S. Tikhonov and A. D. Mirlin, Many-body localization transition with power-law interactions: Statistics of eigenstates, Phys. Rev. B 97, 214205 (2018).
  * Russomanno  _et al._ [2012] A. Russomanno, A. Silva, and G. E. Santoro, Periodic steady regime and interference in a periodically driven quantum system, Phys. Rev. Lett. 109, 257201 (2012).
  * Rudner  _et al._ [2012] M. Rudner, N. Lindner, E. Berg, and M. Levin, Anomalous edge states and the bulk-edge correspondence for periodically driven two-dimensional systems, Physical Review X 3, 031005 (2012).
  * Berry [1983] M. V. Berry, Semiclassical mechanics of regular and irregular motion, in  _Les Houches, Session XXXVI, 1981 — Chaotic Behaviour of Deterministic Systems_ , edited by R. S. G. Ioos, R. H. G. Helleman (North-Holland Publishing Company, 1983) pp. 174–271.
  * Lazarides  _et al._ [2014] A. Lazarides, A. Das, and R. Moessner, Equilibrium states of generic quantum systems subject to periodic driving, Phys. Rev. E 90, 012110 (2014).
  * D’Alessio and Rigol [2014] L. D’Alessio and M. Rigol, Long-time Behavior of Isolated Periodically Driven Interacting Lattice Systems, Phys. Rev. X 4, 041048 (2014).
  * Haake [2001] F. Haake,  _Quantum Signatures of Chaos ( $2^{\rm nd}$ ed)_ (Springer, 2001).
  * Eynard  _et al._ [2018] B. Eynard, T. Kimura, and S. Ribault, Random matrices (2018), arXiv:1510.04430 [math-ph] .
  * Berry and Tabor [1976] M. V. Berry and M. Tabor, Closed orbits and the regular bound spectrum, Proc. Roy. Soc. A349, 101 (1976).
  * Ponte  _et al._ [2015a] P. Ponte, Z. Papić, F. Huveneers, and D. A. Abanin, Many-body localization in periodically driven systems, Physical review letters 114, 140401 (2015a).
  * Serbyn  _et al._ [2013] M. Serbyn, Z. Papić, and D. A. Abanin, Local conservation laws and the structure of the many-body localized states, Phys. Rev. Lett. 111, 127201 (2013).
  * Ros  _et al._ [2015] V. Ros, M. Müller, and A. Scardicchio, Integrals of motion in the many-body localized phase, Nuclear Physics B 891, 420 (2015).
  * Imbrie [2016] J. Z. Imbrie, On many-body localization for quantum spin chains, Journal of Statistical Physics 163, 998 (2016).
  * Wolf  _et al._ [2008a] M. M. Wolf, F. Verstraete, M. B. Hastings, and J. I. Cirac, Area laws in quantum systems: Mutual information and correlations, Phys. Rev. Lett. 100, 070502 (2008a).
  * Jian  _et al._ [2015] C.-M. Jian, I. H. Kim, and X.-L. Qi, Long-range mutual information and topological uncertainty principle (2015), arXiv:1508.07006 [cond-mat.str-el] .
  * von Keyserlingk  _et al._ [2016a] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, Absolute stability and spatiotemporal long-range order in floquet systems, Phys. Rev. B 94, 085112 (2016a).
  * Khemani  _et al._ [2016b] V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phase structure of driven quantum systems, Phys. Rev. Lett. 116, 250401 (2016b).
  * Bachmann  _et al._ [2012] S. Bachmann, S. Michalakis, B. Nachtergaele, and R. Sims, Automorphic equivalence within gapped phases of quantum lattice systems, Communications in Mathematical Physics 309, 835 (2012).
  * Surace [2017] F. M. Surace, Floquet time crystals in clock models, https://etd.adm.unipi.it/theses/available/etd-07032017-183808/unrestricted/tesi.pdf (2017), Master Thesis.
  * Hastings [2010] M. B. Hastings, Quasi-adiabatic continuation for disordered systems: Applications to correlations, lieb-schultz-mattis, and hall conductance (2010), arXiv:1001.5280 [math-ph] .
  * Bauer and Nayak [2013] B. Bauer and C. Nayak, Area laws in a many-body localized state and its implications for topological order, Journal of Statistical Mechanics: Theory and Experiment 2013, P09005 (2013).
  * Ponte  _et al._ [2015b] P. Ponte, Z. Papić, F. m. c. Huveneers, and D. A. Abanin, Many-body localization in periodically driven systems, Phys. Rev. Lett. 114, 140401 (2015b).
  * von Keyserlingk  _et al._ [2016b] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, Absolute stability and spatiotemporal long-range order in floquet systems, Phys. Rev. B 94, 085112 (2016b).
  * not [a] The use of log-averages is here used to better represent the typical value of the spectral pairing, rather than directly looking at its average value, since, due to disorder, various orders of magnitude will compete with each other [75, 76].
  * not [b] We notice that the total amount of correlations in the spin-1 case, as measured by $\Sigma_{\beta}$, is smaller than the spin-1/2 case for the corresponding parameters. This may be explained by noticing that in this case there are other correlations which we have not taken into account, namely, the ones between $(\hat{S}_{k}^{z})^{2}$ at different sites and the ones between $\hat{S}_{k}^{z}$ and $(\hat{S}_{q}^{z})^{2}$. One could use a model-independent quantity to avoid this issue such as the two-site mutual information $I(i,j)=S_{i}+S_{j}-S_{ij}$, where $S_{A}$ is the entanglement entropy of region $A$. Such a quantity upper bounds all normalized connected correlation functions [77], thus including all possible such contributions. Furthermore, one can show (see Appendix E.2) that, even in the spin-$s$ case, at the solvable point, the average mutual information per site has a similar behavior as $\Sigma_{\beta}$ for the solvable point in the spin-1/2 case. Since on average the local imbalance increases with spin-$s$, this confirms that increasing $s$ results in more correlated Floquet states and a stronger time crystal phase.
  * Liu  _et al._ [2023] S. Liu, S.-X. Zhang, C.-Y. Hsieh, S. Zhang, and H. Yao, Discrete time crystal enabled by stark many-body localization, Phys. Rev. Lett. 130, 120403 (2023).
  * Schulz  _et al._ [2019] M. Schulz, C. Hooley, R. Moessner, and F. Pollmann, Stark many-body localization, Physical Review Letters 122, 040606 (2019).
  * Özgüler  _et al._ [2020] A. B. Özgüler, C. Xu, and M. G. Vavilov, Response of a quantum disordered spin system to a local periodic drive, Phys. Rev. B 101, 024204 (2020).
  * Binder and Young [1986] K. Binder and A. P. Young, Spin glasses: Experimental facts, theoretical concepts, and open questions, Rev. Mod. Phys. 58, 801 (1986).
  * Wolf  _et al._ [2008b] M. M. Wolf, F. Verstraete, M. B. Hastings, and J. I. Cirac, Area laws in quantum systems: Mutual information and correlations, Phys. Rev. Lett. 100, 070502 (2008b).