††thanks: These two authors contributed equally to this work
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>††thanks: These two authors contributed equally to this work
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
# Symmetries and conserved quantities of boundary time crystals
in generalized spin models

Giulia Piccitto  SISSA, Via Bonomea 265, I-34135 Trieste, Italy  Matteo
Wauters  Niels Bohr International Academy and Center for Quantum Devices,
Copenhagen University, Universitetsparken 5, 2100 Copenhagen, Denmark  Franco
Nori  Theoretical Quantum Physics Laboratory, RIKEN Cluster for Pioneering
Research, Wakoshi, Saitama 351-0198, Japan  RIKEN Center for Quantum Computing
(RQC), Wakoshi, Saitama 351-0198, Japan  Physics Department, University of
Michigan, Ann Arbor, Michigan, 48109-1040, USA  Nathan Shammah  Unitary Fund
Theoretical Quantum Physics Laboratory, RIKEN Cluster for Pioneering Research,
Wakoshi, Saitama 351-0198, Japan  Quantum Technology Lab, Dipartimento di
Fisica, Università degli Studi di Milano, 20133 Milano, Italy

###### Abstract

We investigate how symmetries and conserved quantities relate to the
occurrence of the boundary time crystal (BTC) phase in a generalized spin
model with Lindblad dissipation. BTCs are a non-equilibrium phase of matter in
which the system, coupled to an external environment, breaks the continuous
time translational invariance. We perform a detailed mean-field study aided by
a finite-size analysis of the quantum model of a $p$,$q$-spin-interaction
system, a generalized $p$-spin-interaction system, which can be implemented in
fully-connected spin-$\frac{1}{2}$ ensembles. We find the following conditions
for the observation of the BTC phase: First, the BTC appears when the discrete
symmetry held by the Hamiltonian, ${\mathbb{Z}}_{2}$ in the considered models,
is explicitly broken by the Lindblad jump operators. Second, the system must
be coupled uniformly to the same bath in order to preserve the total angular
momentum during the time evolution. If these conditions are not satisfied, any
oscillatory behavior appears only as a transient in the dynamics and a time-
independent stationary state is eventually reached. Our results suggest that
these two elements may be general requirements for the observation of a stable
BTC phase relating symmetries and conserved quantities in arbitrary spin
models.

##  I Introduction

Non-equilibrium quantum many-body systems are one of the modern paradigms of
research in quantum physics. Indeed, they represent one of the key elements to
control and manipulate mesoscopic quantum systems and, therefore, they lie at
the root of the growing field of quantum technologies [1, 2, 3, 4]. Moreover
they display a wide range of non-trivial dynamics [5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17] and host a plethora of interesting phenomena with no
analogue in equilibrium states [18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. Of
particular interest are many-body driven-dissipative quantum systems [28, 29],
where dissipation and decoherence can radically change the critical properties
of phase transitions [30, 31] and induce an extremely rich steady-state phase
diagram.

In recent years, considerable attention has been drawn by the so-called time
crystals [32, 33, 34], a phase of matter characterized by a spontaneous
breaking (in the thermodynamic limit) of the time-translational symmetry, that
in short-ranged systems is impossible to realize at equilibrium [35, 36]. In
their most common form, they are a sub-harmonic response to a periodic
driving, thus breaking a discrete time-translational symmetry [37, 38, 39].

More interestingly, a time crystal might appear also when there is no explicit
time dependence in the Hamiltonian but the system is subjected to incoherent
processes [40, 41, 42]. In this particular case, known as Boundary Time
Crystals (BTC), an open system self-organizes in a time-periodic pattern, with
a frequency that depends on the ratio between the model’s relevant energy
scales and the strength of the coupling with the environment. From an
experimental point of view, this phenomenon is very appealing because it
relies on the presence of dissipation, while, usually, time-order is destroyed
by the inevitable decoherence effects [43] present in experimental setups.

BTCs have been first introduced in the context of fully connected spin models,
where permutational invariance [44, 45, 46] allows the exact numerics up to
rather large systems [47, 48, 49, 50, 51, 52, 53, 54, 55] and makes the mean-
field description particularly reliable [44, 56, 57, 58, 59, 60, 61]. In
contrast to discrete time crystals, however, little is still understood about
the conditions needed for the emergence of the BTC phase.

An interesting topic is that of the relation between BTCs and symmetries. The
role of the symmetries in the open-system dynamics has been already under
discussion [62, 63], especially in the case of open systems with decoherence-
free subspaces [64, 65, 66, 67, 68, 62, 69, 70]. Although BTCs share some
phenomenology with them, they have different physical roots. The BTC phase,
indeed, is well defined only in the thermodynamic limit, being associated with
the spontaneous symmetry breaking of the time translational invariance. On the
other hand, if the dynamics has decoherence-free subspaces, a non-stationary
behavior may emerge after an initial relaxation at any system size. In
particular, symmetry-preserving dissipation concurs to split the non-decaying
part of the Hilbert space into disjoint sectors, leading to permanent
oscillations with a frequency directly inherited from the Hamiltonian [71]. In
BTCs, instead, the frequency depends on a nontrivial interplay between the
parent closed system dynamics and the incoherent dissipation or driving.

In Ref. [40], it has been shown that both a free spin in a static magnetic
field and the fully connected Ising model exhibit a transition from a normal
phase to a time-ordered one, driven by the system-environment coupling
strength, when all spins coupled to the same bath through a collective [72]
operator. However, an apparent contradiction emerged in Ref. [73], where no
evidence of time-ordering was found in the closely related $p$-spin model,
where $p$ refers to the order of the highest interaction in the Hamiltonian,
$p\geq 2$.

In this work we solve this seeming contrast and clarify the mechanism behind
the onset of persistent oscillations in the expectation values of physical
observables. To this purpose, we consider a generalized $p$-spin Hamiltonian,
coupled to the environment via Lindblad operators – the model can be
implemented in a fully-connected spin-$\frac{1}{2}$ ensemble with collective
dissipation, analogously to Refs. [40, 73]. The dissipative dynamics is
studied by combining a mean-field approach and an exact finite-size analysis,
where we solve directly the Lindblad master equation [74, 75] for the density
matrix $\dot{\hat{\rho}}=\mathcal{L}[\hat{\rho}]$, where the evolution is
determined by the Liouvillian superoperator $\mathcal{L}$. In this setting,
the BTC phase is manifested by the spontaneous breaking of the time-
translational invariance ${\rm
e}^{\mathcal{L}t}\hat{\rho}_{0}=\hat{\rho}_{0}$. In other words, in the
thermodynamic limit the system does not reach a stationary state, as it does,
instead, for finite sizes.

We find that BTCs survive when two elements are simultaneously present. First,
the Hamiltonian part of the model possesses a discrete symmetry,
${\mathbb{Z}}_{2}$ for the collective spin-$\frac{1}{2}$ system, which must be
explicitly broken by the operators coupling the system with the environment.
The model used in Refs. [73], contrarily to Ref. [40], does not satisfy this
symmetry requirement, hence explaining the absence of BTCs.

Second, the system must have a strong symmetry, namely an observable that
commutes both with the Hamiltonian and the Lindblad operators. This
requirement is fulfilled by including in the open-system dynamics only
collective incoherent processes, as was guessed in Refs. [52, 41], meaning
that all sites need to be coupled uniformly to global baths. Physically, this
is related to the presence of a conserved quantity, which in the considered
spin models is the total angular momentum.

In the absence of either one of the two conditions, we find that the system
always relaxes towards a stationary state, time-translational invariance is
recovered in the long-time limit and all information on the initial condition
is lost. When both are satisfied, instead, the boundary time crystal phase is
present; at the mean-field level, it is identified by the presence of closed
periodic trajectories in the semiclassical phase space, while in finite sizes
it is revealed by oscillations with a damping rate diverging in the
thermodynamic limit. Supported by numerical evidence in a study including
dissipation to ever increasing bath sizes – ranging from local to collective–,
we argue that these elements may be necessary conditions for the existence of
BTCs in general spin systems, beyond homogeneous fully-connected spin models.

Interestingly, also the finite-size density matrix $\hat{\rho}_{0}$ of the
stationary state can identify the BTC phase: when time-order is present in the
thermodynamic limit, $\hat{\rho}_{0}$ cannot be approximated by a well-defined
mean-field ansatz, but rather it is the average over all the states on the
time-crystal trajectory [76].

The article is organized as follows: In Sec. II we introduce the model we
investigate and its mean-field description. In Sec. III we review some of the
main results known in the literature within the more general framework we
present in this paper. Our main results are presented in Sec. IV, where we
focus on the symmetry requirements to the BTC phase, and in Sec. V, where we
show the effects of non-collective Lindblad operators. In Sec.VII we discuss
possible experimental realizations of boundary time-crystals. Finally, we
summarize our findings in Sec. VIII.

##  II Theoretical framework

We consider a generalization of the well know $p$-interacting spin model [77,
78, 79, 80] described by the Hamiltonian ($\hbar=1$)

$\displaystyle\hat{H}$
$\displaystyle=-N\left(\omega_{z}\hat{J}_{z}^{p}+\omega_{x}\hat{J}_{x}^{q}\right),$
(1)

with $p,q\in\mathbb{N}$ and the spin algebra is given by
$[\hat{J}_{x},\hat{J}_{y}]=i2\hat{J}_{z}/N$,
$[\hat{J}_{+},\hat{J}_{-}]=2\hat{J}_{z}/N$, where, without loss of generality,
we consider these as collective spins from a system of $N$ all-to-all
interacting spin-$\frac{1}{2}$.
$\hat{J}_{\alpha}=\sum_{i}\hat{\sigma}^{\alpha}_{i}/N$ are the (collective)
magnetization operators and $\hat{\sigma}_{i}^{\alpha}$ are the Pauli matrices
acting on the $i$-th site. The dynamics we consider is given by the Lindblad
master equation [74, 81]

$\displaystyle\frac{{\rm d}}{{\rm d}t}{\hat{\rho}}=-$ $\displaystyle
i[\hat{H},\hat{\rho}]+N\Gamma^{\uparrow}\left(\hat{J}^{+}\hat{\rho}\hat{J}^{-}-\frac{1}{2}\\{\hat{J}^{-}\hat{J}^{+},\hat{\rho}\\}\right)$
(2) $\displaystyle+$ $\displaystyle
N\Gamma^{\downarrow}\left(\hat{J}^{-}\hat{\rho}\hat{J}^{+}-\frac{1}{2}\\{\hat{J}^{+}\hat{J}^{-},\hat{\rho}\\}\right)\
,$

where $\hat{J}^{\pm}=\hat{J}_{x}\pm i\hat{J}_{y}$ describe two collective
incoherent processes with associated rates $\Gamma^{\uparrow(\downarrow)}$.
Note that Eq. (2) in general couples the spin ensemble to two separate baths,
as the $\Gamma^{\uparrow(\downarrow)}$ rates can be independently chosen. This
includes the case of a single bath at detailed balance, in which
$\Gamma^{\uparrow}=\Gamma_{0}\left(1+n_{T}\right)$, and
$\Gamma^{\downarrow}=\Gamma_{0}n_{T}$ due to the favored spin alignment in the
upward direction of Eq. (1), with $n_{T}$ the thermal occupation number, which
goes to zero at $T=0$, and $\Gamma_{0}$ a coefficient fixed for the considered
model. Depending on the largest between $\Gamma^{\downarrow}$ and
$\Gamma^{\uparrow}$, the Lindblad operators favor spin alignment along the
upward ($\hat{J}^{+}$) or the downward ($\hat{J}^{-}$) $z$ direction,
explicitly breaking the ${\mathbb{Z}}_{2}$ symmetry the Hamiltonian displays
for $p$ even.

Notice that, since
$[\hat{H},\hat{J^{2}}]=[\hat{J}_{-},\hat{J^{2}}]=[\hat{J}_{+},\hat{J^{2}}]=0$,
the total angular momentum $\hat{J^{2}}$ is a _strong symmetry_ [65, 66] of
the set of generalized spin models considered and Eq. (2) preserves the total
angular momentum $\hat{J^{2}}$111We want to remark that Noether theorem can
not be generalized at the operator level in the case of open quantum systems.
For this reason there is no one-to-one correspondence between symmetries and
conserved quantities. In particular, it is possible to find conserved
quantities that are not associated to any symmetry, but any strong symmetry in
the model has an associated conserved quantity[66].. Moreover, besides for the
treatment of local baths performed in Sec. V, hereafter we can generalize the
models beyond their derivation in terms of collections of spin-$\frac{1}{2}$,
and consider them genuine $(N+1)$-spins, further opening up quantum simulation
experimental possibilities, as detailed in Sec. VII.

The mean-field equations of motion for the expectation values of the
magnetization operators
$X,Y,Z=\braket{\hat{J}_{x}},\braket{\hat{J}_{y}},\braket{\hat{J}_{z}}$ are
obtained with a Gutzwiller approximation,
$\braket{\hat{J}_{\alpha}\hat{J}_{\beta}}=\braket{\hat{J}_{\alpha}}\braket{\hat{J}_{\beta}}$,
where $\alpha,\beta=x,y,z$ and $\braket{\cdot}=\text{Tr}[\hat{\rho}\ \cdot]$
is the expectation value on the state $\hat{\rho}$. With this approximation,
the equations for the macroscopic variables are

$\left\\{\begin{aligned} &\dot{X}=2p\omega_{z}Z^{p-1}Y-2\delta\Gamma ZX,\\\
&\dot{Y}=2XZ\left(q\omega_{x}X^{q-2}-p\omega_{z}Z^{p-2}\right)-2\delta\Gamma
ZY,\\\
&\dot{Z}=-2q\omega_{x}YX^{q-1}+2\delta\Gamma(1-Z^{2}),\end{aligned}\right.$
(3)

with $\delta\Gamma=\Gamma^{\uparrow}-\Gamma^{\downarrow}$ (the full derivation
can be found in App. A).

We introduce in notation $\dot{\bf R}=(\dot{X},\dot{Y},\dot{Z})$ and study the
solutions of this system for $\dot{\bf R}=0$, i.e. ${\bf
R}=(X_{\text{st}},Y_{\text{st}},Z_{\text{st}})$ are the stationary states of
the dynamics. In what follows, we refer to a _ferromagnetic_ stationary state
if $|Z_{\text{st}}|>0$ and to a _paramagnetic_ one if $Z_{\text{st}}=0$. The
real part of the eigenvalues of the associated Jacobian, i.e. the real part of
the Lyapunov exponents, provides a classification of the stationary states: if
all the eigenvalues have positive (negative) real part, the trajectories are
attracted toward (repelled from) the relative state. More interestingly, if
the real part is zero and the Lyapunov exponents are purely imaginary 222The
presence of purely imaginary Lyapunov exponents is necessary for the existence
of the BTCs but it is not sufficient. For instance, it is possible to have
purely imaginary Lyapunov exponents associated to slow spiraling trajectories
in which the damping of the oscillations is a non-linear effect of the
dynamics. we are in presence of a marginal fixed point generating periodic
orbits, and the stationary state is associated with a BTC trajectory 333Notice
that in the BTC phase every initial condition is associated to a different
periodic orbit that keeps memory of the evolution of the initial state. This a
substantial difference with the dynamics in presence of limit cycles in which
all the trajectories eventually reach the same periodic orbit independently on
the initial conditions..

When studying a finite size system, a fundamental role is played by the
structure of the Liouvillian superoperator [81, 85, 86, 87, 88, 89]
$\mathcal{L}$, defined through Eq. (2) as
$\dot{\hat{\rho}}=\mathcal{L}[\rho]$. For a finite system size $N$, in
general, the spectrum of $\mathcal{L}$ is gapped and the dynamic relaxes
towards a stationary state $\hat{\rho}_{0}$, defined as the right eigenvector
of $\mathcal{L}$ with eigenvalue $\lambda_{0}=0$. The other eigenvectors,
instead, are associated with eigenvalues $\lambda_{i\neq 0}$ with
$\mathbb{R}\mathrm{e}\left[\lambda_{i\neq 0}\right]<0$ and give information on
the transient dynamics. For instance, the real part of the first nonzero
eigenvalue $|\mathbb{R}\mathrm{e}\left[\lambda_{1}\right]|$, the Liouvillian
gap, describes the relaxation rate.

The BTC phase, instead, is characterized by the presence of at least one
complex eigenvalue (and its complex conjugate) whose real part vanishes in the
thermodynamic limit, while the imaginary part saturates to a constant value
[90, 76]. In this case, the real and the imaginary parts of the eigenvalues
are associated with the damping of the oscillations in finite systems and to
their frequency, respectively.

##  III Boundary time crystals in $p,q$-interacting spin models: An overview

To exemplify our findings, we start with reviewing the results known in the
literature: the Ising Hamiltonian with $z$–interactions ($p=2$ and $q=1$), and
the Ising Hamiltonian with $x$–interactions ($p=1$ and $q=2$), both in the
presence of dissipation, which have been studied in Refs. [40] and [73],
respectively. Here we want to summarize the main features that will be
essential for what comes next.

###  III.1 Free spins ($p=1,q=1$)

Before proceeding let us briefly comment on the case of free spins in a
magnetic field lying in the $x$-$z$-plane. The Hamiltonian reads ($p=1,q=1$ in
Eq. (1))

$\hat{H}^{\text{free}}=-N\left(\omega_{z}\hat{J}_{z}+\omega_{x}\hat{J}_{x}\right).$
(4)

Without loss of generality, we take $\omega_{x},\ \omega_{z}\geq 0$ throughout
the whole paper.

The limit $\omega_{x}=0$ is trivial: we have two possible ferromagnetic
stationary solutions ${\bf R}=\left(0,\ 0,\ \pm 1\right)$, obtained by solving
$\dot{\bf R}=0$; one of them is stable and the other is unstable, depending on
the sign of $\delta\Gamma$. When $\omega_{z}=0$, as discussed in [40], the
system has two different phases divided by the critical value of the
dissipation $\delta\Gamma_{c}=\omega_{x}$: for $\delta\Gamma<\delta\Gamma_{c}$
the system is in a BTCs phase, while for $\delta\Gamma>\delta\Gamma_{c}$ the
system again has two stationary states ${\bf R}=\left(0,\ 0,\ \pm 1\right)$.
As soon as $\omega_{z}>0$, the system falls inevitably in the ferromagnetic
phase, where the time crystal order is destroyed and there are only trivial
solutions with magnetization
$Z=\omega_{z}/\sqrt{\omega_{x}^{2}+\omega_{z}^{2}-\delta\Gamma^{2}}$.

As a side remark, we would like to comment on the analogies between the mean-
field limit in Eq. (3), for $p=q=1$, and the Lorenz equations [91], which are
one of the best known examples of deterministic chaos in classical systems.
Despite this similarity, our model we cannot have chaotic behavior because of
the conservation of angular momentum $X^{2}+Y^{2}+Z^{2}=1$, due to the
collective jump operators, that reduces the effective dimension of the phase
space. Chaotic dynamics, however, might be recovered by introducing an
explicit time dependence in the equations of motion.

###  III.2 $z$–interactions ($p=2,q=1$)

Let us now move to the fully connected Ising chain in transverse field, with
interactions along the $z$ direction described by the Hamiltonian ($p=2,q=1$
in Eq. (1))

$\displaystyle\hat{H}^{(z)}$
$\displaystyle=-N\left(\omega_{z}\hat{J}_{z}^{2}+\omega_{x}\hat{J}_{x}\right).$
(5)

Figure 1: Phase diagram of the system in Eq. (5), with $z$-interactions and
$x$-transverse field in the Hamiltonian and collective dissipation greater
than collective incoherent pumping, $\delta\Gamma>0$ . We plot the class of
solutions as a function of the transverse field $\omega_{x}$ and the
dissipation rate $\delta\Gamma$ both in units of $\omega_{z}$. F means that we
have ferromagnetic solutions only, BTCs that we have marginal solutions only,
F + BTCs means that both phases coexist. The red square marks the parameters
used in Fig. 2.

Without loss of generality, let us assume $\omega_{x}>0$. The mean-field
equations in Eq. (3) have four fixed points, two ferromagnetic solutions, and
two paramagnetic ones.

The ferromagnetic solutions are given by

${\bf
F_{\pm}}=\left(\frac{2\omega_{x}\omega_{z}}{4\omega_{z}^{2}+\delta\Gamma^{2}},\
\frac{\delta\Gamma\omega_{x}}{4\omega_{z}^{2}+\delta\Gamma^{2}},\
\pm\sqrt{1-\frac{\omega_{x}^{2}}{4\omega_{z}^{2}+\delta\Gamma^{2}}}\right),$
(6)

one of which is stable and the other unstable, in fact the system is attracted
toward the positive or negative magnetization state depending on the sign of
$\delta\Gamma$. From the constraint $X^{2}+Y^{2}+Z^{2}=1$ we deduce the
existence condition for these solutions
$\omega_{x}<\omega_{x_{c}}=\sqrt{4\omega_{z}^{2}+\delta\Gamma^{2}}$.

The two paramagnetic solutions are

${\bf P_{\pm}}=\left(\pm\sqrt{1-\frac{\delta\Gamma^{2}}{\omega_{x}^{2}}},\
\frac{\delta\Gamma}{\omega_{x}},\ 0\right),$ (7)

which exist only when the dissipation rate is smaller than the transverse
field $\delta\Gamma\leq\omega_{x}$. Among the two, the most interesting is the
negative one, since it is the marginal point that acts as the generator of the
periodic orbits characterizing the BTC phase.

Having gathered this information, we can draw the phase diagram in Fig. 1, in
which we show the phases of the system as a function of $\omega_{x}$ and
$\delta\Gamma$, both in unit of $\omega_{z}$: the system has a ferromagnetic
phase (F) for $\delta\Gamma>\omega_{x}$, a BTCs phase for
$\omega_{x}>\sqrt{4\omega_{z}^{2}+\delta\Gamma}$ and a region where they
coexist. Note that values of $\delta\Gamma>\omega_{z}$ are attainable in
driven systems whose interaction-picture Hamiltonian takes the form of Eq.
(3), even if dissipation mechanisms are perturbative effects [92]. In the
coexistence phase the dynamics depends on the initial conditions: the closer
the parameters $\omega_{x}$ and $\delta\Gamma$ are to the purely ferromagnetic
phase, the larger is the basin of attraction for the corresponding fixed
point. In the opposite situation, close to the boundary with the BTC phase,
most of the initial conditions lead to periodic trajectories while a small
number of them relaxes towards the ferromagnetic point.

Figure 2: Phase portrait of the dynamical system in Eq. (3) for $p=2$ and
$q=1$. The parameters $\omega_{z}=\omega_{x}=1$, $\delta\Gamma=0.2\omega_{x}$
are chosen to be in the coexistence region. The black lines are the
trajectories associated to different initial conditions, while the red dots
are the stationary states of the system, two ferromagnetic solutions and two
paramagnetic ones. Of the two ferromagnetic solutions, the positive one (solid
circle) is stable and attracts the trajectories, while the negative one
(dashed circle) is unstable and therefore repels them. The negative
paramagnetic solutions (solid square) generate periodic orbits associated to
the BTC behavior, while the positive paramagnetic one (dashed square)
separates the ferromagnetic basin of attraction from the BTC one.

A more intuitive way to visualize this information is by looking at the phase
portrait obtained by studying the dynamics in polar coordinates, which
correspond to the parametrization

${\bf R}=\left(\sin\theta\cos\varphi,\ \sin\theta\sin\varphi,\
\cos\theta\right),$ (8)

see appendix B for more details on the derivation.

In Fig. 2 we show the phase portrait for the phase coexistence region, with
$\omega_{x}=\omega_{z}$ and $\delta\Gamma=0.2\omega_{z}$ as a function of the
polar angle $\varphi$ and its conjugate variable $\cos\theta$, with $\theta$
the azimuthal angle. The red dots are the solutions of the system in Eq. (3),
with the left-hand side equal to zero, while the grey lines are trajectories
associated with different initial conditions. Each fixed point clearly
influences the trajectories in a different way: the ferromagnetic state with
positive (solid circle) magnetization ($\cos\theta>0$) is stable and therefore
attracts the nearby trajectories, while the negative(dashed circle) one repels
them. The marginal paramagnetic (solid square) solution at
$(\varphi=\pi,\cos\theta=0)$, instead, generates the periodic orbits
associated with the time crystal behavior. Finally, there is a second
paramagnetic “saddle” (dashed square) fixed point at
$(\varphi=0,\cos\theta=0)$ which separates the region of influence of all
other stationary solutions.

###  III.3 $x$–interactions ($p=1,q=2$)

Let us now switch the direction of the interaction and the transverse field in
the Hamiltonian, i.e. we consider a fully connected Ising chain in transverse
field with interactions along the $\hat{\bf x}$ direction, described by
($p=1,q=2$ in Eq. (1)):

$\displaystyle\hat{H}^{(x)}$
$\displaystyle=-N\left(\omega_{x}\hat{J}_{x}^{2}+\omega_{z}\hat{J}_{z}\right)\
.$ (9)

Although the closed systems described by Eqs. (5) and (9) are identical, the
specific form of the Lindblad operators leads to very different dynamics when
coupling to the environment. In this case, there is a single line in the phase
diagram with time crystal order, corresponding to $\omega_{z}=0$: as soon as
$\omega_{z}>0$ BTCs are destroyed and the system relaxes toward a stable fixed
point, that can be either x-paramagnetic ($X=0$), or x-ferromagnetic
($|X|>0$). In Fig. 3 we show the phase portrait for $\omega_{z}=\omega_{x}=1$
and $\delta\Gamma=0.2\omega_{x}$. In order to avoid the singularities at the
poles of the unit sphere, here we use a different parametrization with respect
to standard polar coordinates

${\bf R}=\left(\cos\theta^{\prime},\ \sin\theta^{\prime}\cos\varphi^{\prime},\
\sin\theta\sin\varphi^{\prime}\right)\,$ (10)

where $\varphi^{\prime}$ and $\theta^{\prime}$ are the polar and azimuthal
angle defined with respect to the $x$-axis. Looking at the phase portrait in
panel (a) of Fig. 3 we identify two attractive x-ferromagnetic states (solid
circles), reflecting the ${\mathbb{Z}}_{2}$ symmetry of the Hamiltonian, and
one repelling x-paramagnetic state (dashed circle) with $X=0$ and $Z=-1$.

Although, at first sight, it might appear that this point is surrounded by
closed orbits, actually the trajectories are escaping from the fixed point,
describing a _spiral_. This is more evident by looking at the behavior on the
Bloch sphere in panel (b) of Fig. 3 of the corresponding blue trajectory in
the phase space. Thus, in contrast with the previous situation, there are no
periodic orbits associated with a time-crystalline phase. This holds for any
nonzero value of $\omega_{z}$: the stronger the transverse field, the faster
the trajectories will collapse. Finally, the last x-paramagnetic solution
(dashed square) separates the basin of attraction of the positive
x-ferromagnetic solution to that of the negative x-ferroamagnetic one.

Figure 3: Panel (a): Phase portrait of the dynamical system in Eq. (3) for
$p=1$ and $q=2$. The parameters $\omega_{z}=\omega_{x}$,
$\delta\Gamma=0.2\omega_{x}$ correspond to a ferromagnetic phase. The black
lines are the trajectories associated to different initial conditions, while
the red dots are the stationary states of the system, two x-ferromagnetic
($|X|>0$) solutions (solid circles) and two x-paramagnetic ($X=0$) ones. We
notice that the x-ferromagnetic solutions are both stable, reflecting the
${\mathbb{Z}}_{2}$ symmetry of the Hamiltonian. The negative x-paramagnetic
solution (dashed circle) is a repeller, while the other one (dashed square)
separates the positive x-ferromagnetic basin of attraction from the negative
x-ferromagnetic one. The Bloch sphere in panel (b) shows a particular
trajectory which well represent the behavior close to attractive and repulsive
fixed points. The same trajectory is depicted with a continuous blue line also
in the main panel.

##  IV BTCs existence condition 1: Symmetry

In this section, we present the mean-field and exact finite-$N$ dynamics of
the generalized $p$-spin model. From our results, it emerges that to observe
BTCs the Hamiltonian needs to be ${\mathbb{Z}}_{2}$ invariant and the Lindblad
operators must break explicitly this symmetry. This implies that the BTC phase
exists only for even values of $p$, while $q$ can be either even or odd.

Although in closed system there is a perfect correspondence when exchanging
$p\longleftrightarrow q$ and $\omega_{z}\longleftrightarrow\omega_{x}$, this
is no longer true in the presence of dissipation. In the case of even $p$, the
${\mathbb{Z}}_{2}$ symmetry is generated by
$G_{z}=\prod_{j}\hat{\sigma}^{x}_{j}$. This operator commutes with the
Hamiltonian but acts non-trivially on the jump operators
$G_{z}\hat{J}^{\pm}G^{\dagger}_{z}=\hat{J}^{\mp}$. This action corresponds to
an effective switch of $\Gamma^{\uparrow}$ and $\Gamma^{\downarrow}$ and,
consequently, to an effective switch of the preferred alignment of the spins,
breaking the ${\mathbb{Z}}_{2}$ symmetry of the system. In the case of even
$q$, instead, the symmetry generator is $G_{x}=\prod_{j}\hat{\sigma}^{z}_{j}$
and its action on the jump operators is
$G_{x}\hat{J}^{\pm}G^{\dagger}_{x}=-\hat{J}^{\pm}$, hence $G_{x}$ does not
commute with $\hat{J}^{\pm}$. Despite this, we observe that, differently from
the previous case, the Lindblad equation remains unchanged under the effect of
the symmetry operator (the jump operators always come in pairs and the minus
signs cancel), leaving a weak ${\mathbb{Z}}_{2}$ symmetry in the open system.

Figure 4: Top panels: Phase portraits obtained by solving Eq. (3) with $q=1$,
$p=3$ (a) and $p=4$ (b), fixed $\omega_{z}=\omega_{x}$ and
$\delta\Gamma=0.2\omega_{z}$. The red trajectories correspond to the initial
condition $\varphi_{0}=0.10$ and $\theta_{0}=1.47$. Both diagrams present two
ferromagnetic solutions and two paramagnetic ones. The ferromagnetic solutions
are in both cases one attractive and the other repelling trajectories. The
nature of the paramagnetic solutions, instead, depends on the value of $p$. In
fact, only for $p=4$ (d), namely the ${\mathbb{Z}}_{2}$-symmetric case, we
have a marginal fixed point generating periodic orbits, while for $p=3$ (c) it
slowly attracts the trajectories. Bottom panels: typical trajectories depicted
on the Bloch sphere for $p=3$ (c) and $p=4$ (d); the same trajectories are
highlighted with the same color also in the phase portraits in (a) and (b).
The “thickness” of the red trajectory in (c) signals that it is slowly
spiralling towards the $X=1$ fixed point.

###  IV.1 Mean-field analysis

Here, we explore the mean-field dynamics of the generalized $p$-spin model in
Eq. (1) to discuss the symmetry conditions necessary to observe BTCs. As
explained in Sec. II, the Hamiltonian in Eq. (1) displays a ${\mathbb{Z}}_{2}$
symmetry on the $z$($x$) direction for even values of $p$($q$). The Lindblad
operators $\hat{J}^{+}$ and $\hat{J}^{-}$ favor the positive or negative
alignment of the spins along the $z$ direction depending on the sign of
$\delta\Gamma=\Gamma^{\uparrow}-\Gamma^{\downarrow}$. Our claim is that, for
such a dissipation, the BTC phase appears only in the presence of
${\mathbb{Z}}_{2}$ symmetry along the $z$ direction. Hence, we expect to
observe BTCs for even $p$, independently on the value of $q$.

To prove this, we present two different scenarios: First we consider no
symmetries along the $x$ directions ($q=1$) and we show that BTCs arise only
for even $p$; then we fix $p=2$ and we observe that BTCs exist for any value
of $q$.

The more intuitive way to do so is by looking at the phase portrait for
different values of $p$ and $q$. From now on we will parametrize the sphere as
in Eq. (8). All the results (unless specified) are obtained with
$\omega_{x}=\omega_{z}$ and $\delta\Gamma=0.2\omega_{z}$.

First, we fixed $q=1$. In the top panels of Fig. 4 we show the phase portrait
for $p=3$ (a) and $p=4$ (b), a non-symmetric and to a
${\mathbb{Z}}_{2}$-symmetric Hamiltonian, respectively. At first sight, the
two figures seem to share the same physics: the trajectories are either
attracted toward a positive ferromagnetic stationary state or stuck into
periodic orbits.

Figure 5: Top panels: phase portraits for $p=2$ and $q=1$ (a), $q=2$ (b),
$q=3$ (c), $q=4$ (d). We set $\omega_{x}=\omega_{z}$,
$\delta\Gamma=0.2\omega_{z}$. For all the values of $q$, we observe periodic
orbits, indicating BTCs. The blue and orange squares in panel (a) Fig. 5 mark
the parameters used in Fig. 7. The bottom panels show the behavior in time of
the longitudinal magnetization evaluated over the red trajectories in the
relative phase portrait (corresponding to the initial conditions
$\varphi_{0}=3.10$, $\theta_{0}=1.47$ for $q=1$ (e) and $\varphi_{0}=1.57$,
$\theta_{0}=1.47$ for $q=2$ (f), $q=3$ (g), $q=4$ (h)). We observe that the
oscillations persist in time independently on the value of $q$ confirming that
the BTC phase is not affected by the interactions along the $\hat{\bf x}$
direction.

Actually, this is not valid in the case $p=3$ because the trajectories
surrounding the paramagnetic solutions, as shown in panel (c), are very slowly
dampedtoward a time-independent value.

In the bottom panels of Fig. 4 we plot the trajectories of the longitudinal
magnetization, that in the phase portrait are highlited in blue and red on the
Bloch sphere. From these figures emerge that, while in the case $p=4$ (d) we
have a closed orbit, for $p=3$ (c) the trajectories are attracted toward the
ferromagnetic solution with damping (which can be very slow, depending on the
initial conditions) that eventually suppresses the oscillations, ruling out
the possibility of BTCs. An analysis of the oscillations amplitude $A(t)$
shows, indeed, that for $p>1$ it decreases in time with a $p$-dependent power
law. When $q=1$ the scaling is

$A(t)=B\ t^{-\frac{1}{p-1}}\ ,$ (11)

where the precise value of $B$ depends on the other system parameters.

Let us now consider the case of a ${\mathbb{Z}}_{2}$ preserving Hamiltonian by
fixing $p=2$. In the top panels of Fig. 5 we plot the phase portraits for
$q=1,2,3,4$. Although the complexity of the phase portrait increases by
increasing $q$, there are always some periodic orbits signaling the presence
of BTCs. The blue and the orange diamonds in panel (a) Fig. 5 mark the
parameters used in Fig. 7. The bottom panels show the trajectories in time of
the longitudinal magnetization evaluated over the red trajectories, confirming
the persistence of the oscillations and the emerging of a time-crystal order.

Although we do not have rigorous proof, our numerical analysis strongly
suggests our claim on the symmetry requirement: if the dissipation does not
break the discrete symmetry of the Hamiltonian, closed periodic orbits do not
appear.

To conclude this section, we remark that the presence of the BTC phase,
regardless of the value of $q$, has an important implication on the robustness
of the phase itself. When $q=1$, the $x$-term in the Hamiltonian can be
interpreted as a transverse field or coherent driving. In this case, any
fluctuations on its direction could break the ${\mathbb{Z}}_{2}$ symmetry and
inevitably destroy the BTC phase. If $q=p=2$, instead, the Hamiltonian
contains only interaction terms and no external field or driving. In this case
the interactions concur to stabilize the BTCs phase since the symmetry
${\mathbb{Z}}_{2}\times{\mathbb{Z}}_{2}$ cannot be broken by noise in the
Hamiltonian parameters.

###  IV.2 Finite-$N$ analysis

To corroborate the mean-field analysis done so far, we solve the exact quantum
dynamics in finite systems and present three indicators that identify the BTC
phase: (i) the oscillations of the expectation value of the magnetization and
the link between the damping rate and the system size $N$, (ii) the part of
the Liouvillian spectrum closest to the origin of the complex plane, (iii) the
structure of the stationary state density matrix $\hat{\rho}_{0}$.

####  IV.2.1 Oscillations of the magnetization

Since the Lindblad equation conserves the total spin, we can restrict the
analysis to the sector of the Hilbert space with total spin
$S^{2}=\frac{N}{2}\left(\frac{N}{2}+1\right)$. This allows us to access
numerically systems of the order $N\sim 10^{2}-10^{3}$, depending on whether
we are interested in the Liouvillian spectrum or only in the time evolution.
All numerical results have been obtained using the QuTiP [93, 94] Python
package. In what follows, we mainly focus our attention on the Hamiltonian in
Eq. (5), corresponding to the choice $p=2$, $q=1$.

First, let us consider the finite-$N$ dynamics in the BTC phase. Fig. 6(a)
shows the main features of a typical trajectory of the magnetization in the
BTCs phase. Differently from the mean-field case, the magnetization presents
oscillations decaying with a damping strength that decreases for increasing
system size $N$, which finally vanishes in the thermodynamic limit. By
performing a data collapse of the oscillation amplitudes, we found that the
decay rate follows a power-law compatible with $N^{-0.4}$, as can be observed
in Fig. 6(b). This exponent, however, is non-universal and it depends on $p$
and $q$. For instance, in the free spin case in Eq. (4) with $\omega_{z}=0$
studied in detail in Ref. [40], the decay rate decreases as $N^{-1}$.

The period of the oscillations, instead, is practically size-independent and,
in fact, the Fourier spectrum is always peaked around the same frequency with
a small broadening that reduces by increasing $N$. In the BTC phase, this
behavior is independent of the initial condition and, qualitatively, it holds
$\forall p$ (even).

Figure 6: (a): Magnetization $\langle\hat{J}_{z}\rangle$ versus time $t$ in
the BTC phase, for several values of the system size $N$ and interaction ranks
$p=2$, $q=1$. The dashed line corresponds to the trajectory in the
thermodynamic limit, obtained through the mean-field equations (3). The other
parameters of the model are $\omega_{x}=3\omega_{z}$ and
$\delta\Gamma=0.2\omega_{z}$. (b): Oscillations amplitude as a function of the
rescaled time $\tilde{t}_{N}=tN^{-0.4}$. The nice collapse of the different
curves suggests that the decay time has a power-law divergence $N^{0.4}$.

In the coexistence phase (BTC+F), periodic orbits and relaxation dynamics
coexist in the thermodynamic limit. When one focuses on finite-size effects,
this behavior persists for a finite-time window, as shown in Fig. 7. Depending
on the initial conditions, the trajectory might display both $N$-dependent
damped oscillations or a size-independent relaxation toward the ferromagnetic
stationary state. Eventually, this apparent bistability breaks after a time
scale $\tau_{N}\sim N$, when the trajectory deviates from the time-crystal
density matrix to slowly approach the unique steady state (inset of Fig. 7).
The relative Liuvillian spectrum is qualitatively similar to that of the BTC
phase with a stationary state that is ferromagnetic instead of paramagnetic.

Figure 7: Comparison of the time dependence of the average magnetization for
two different initial conditions in the coexistence phase: all spins aligned
along the direction $(\theta_{0},\phi_{0})=(-\frac{\pi}{6},0)$ (solid blue
curve) and $(\theta_{0},\phi_{0})=(\frac{\pi}{2},0)$ (dashed orange curve).
The system size is $N=100$ and the other Hamiltonian parameters are $p=2$,
$q=1$, $\omega_{x}=\omega_{z}$. The dissipation is
$\delta\Gamma=0.1\omega_{z}$. The inset shows the same date on a longer time
scale, where one can see the slow relaxation of the trajectory corresponding
to the boundary time crystal towards the ferromagnetic stationary state.

####  IV.2.2 Liouvillian spectrum

The properties depicted above are reflected in the eigenspectrum of the
Liouvillian. In Fig. 8 we show the behavior of the real and the imaginary part
(top and bottom panels, respectively) of the eigenvalues as a function of the
system size $N$ for the BTCs phase, in (a) and (c), and the ferromagnetic
phase in (b) and (d).

In the BTC phase, see Fig. 8(a), the real part of the spectrum roughly divides
into a set of real eigenvalues that decreases as
$\mathbb{R}\mathrm{e}\left[\lambda^{(N)}\right]\sim N^{-1}$, and a set of
complex eigenvalues associated with the BTCs (blue circles) with a much slower
decay. Again, the only exception is the free spin case in which, as mentioned,
we find a power-law behavior. Because of the computational effort needed to
diagonalize the Liouvillian superoperator, it is hard to extract a precise
scaling; however, we found that its behavior is compatible with the power-law
decay $N^{-0.4}$, in agreement with the damped oscillations shown in Fig. 6.

The imaginary part in Fig. 8(c), instead, clearly saturates to values
independent from $N$. Surprisingly, the oscillation frequency of the
magnetization is not given by the “quantization” of the imaginary part of the
spectrum, in contrast to to what happens in the free case. In absence of
interactions, i.e. $p=0$, $q=1$, the frequency is indeed independent on the
initial condition and can be extracted directly from the Liouvillian [40]. In
the present case, instead, the $\hat{J}_{z}^{p}$ interaction introduces a
dependence of the frequency on the initial conditions, which might also be due
to the presence of Liouvillian eigenstates with imaginary parts that are not
perfectly commensurate.

In (b) and (d), we plot the Liouvillian eigenvalues for the ferromagnetic
phase. In this case, the behavior of the real part of the eigenvalues is
markedly different from the previous one. There is still at least one real
eigenvalue decreasing as $N^{-1}$, suggesting the possibility of multiple
steady-states in the thermodynamic limit. Despite this, the Louvillian gap is
much larger than the corresponding one in the BTC phase, as one can appreciate
from the different scales of the two vertical axes. More importantly, the real
part of the lowest complex eigenvalues (blue circles) do not display any
dependence on $N$, which means that any oscillatory behavior is doomed to
decay also in the thermodynamic limit.

Figure 8: Real and imaginary parts vs $1/N$ of the 21 eigenvalues of the
Liouvillian with $\mathbb{R}\mathrm{e}\left[\lambda^{(N)}_{i}\right]$ closest
to zero. The large blue circles in the upper panels highlight eigenvalues with
a nonzero imaginary part. The dashed blue line in (a) and highlights the $1/N$
scaling of the Liouvillian gap (green squares); the solid red line shows the
possible $N^{-0.4}$ scaling of the eigenvalues with nonzero imaginary parts.
(a) and (c): BTC phase with $p=2$, $q=1$, $\omega_{x}=3\omega_{z}$, and
$\delta\Gamma=0.2\omega_{z}$. (b) and (d): ferromagnetic phase with $p=2$,
$q=1$, $\omega_{x}=0.25\omega_{z}$, and $\delta\Gamma=0.5\omega_{z}$.

####  IV.2.3 Steady-state density matrix

Finally, we would like to briefly comment on the stationary state density
matrix $\hat{\rho}_{0}$. When the system is in the BTC phase, the density
matrix $\hat{\rho}_{0}$ of the stationary state for finite size $N$ is
independent on the choice of individual parameters and acquires the
characteristic structure shown in Fig. 9(a). Besides some corrections of the
order $O(N^{-2})$, $\hat{\rho}_{0}$ is proportional to the identity matrix,
where all eigenstates of $\hat{S}_{z}$ are populated uniformly. This suggests
that there is an emergent ${\mathbb{Z}}_{N+1}$ symmetry, generated by a
unitary ladder operator $\hat{Q}$ that acts on the eigenstates $\ket{m}$ of
$\hat{S}_{z}$ as $\hat{Q}\ket{m}=\ket{m+1}$.

At finite size $N$ this symmetry is not exact, but it gets closer to an actual
${\mathbb{Z}}_{N+1}$ invariance for increasing values of $N$, until it becomes
a U$(1)$ symmetry in the thermodynamic limit. A similar mechanism to simulate
a U(1) lattice gauge theory has been proposed in Ref. [95]. In the present
case, we argue that this U$(1)$ symmetry is associated to the closed orbits
arising in the mean-field description of the thermodynamic limit, where the
dynamics keeps the memory of the initial state, modulo a rotation. This result
is further validated by the correspondence between the emergence of time
crystals and dissipative phase transitions as shown in Ref. [69].

As a consequence of its structure, the BTC steady-state density matrix cannot
be described with a mean-field approximation, meaning that it cannot be
written as a direct product of single spin density matrices, since the
fluctuations of the magnetization components diverge in the thermodynamic
limit. By a simple calculation, see App. C, it is possible to show that it
does not exist a factorized density matrix $\hat{\rho}_{0}$ giving
$\mathrm{Tr}\\!\left[\hat{\rho}_{0}\hat{S}_{x}\right]\simeq\mathrm{Tr}\\!\left[\hat{\rho}_{0}\hat{S}_{y}\right]\simeq\mathrm{Tr}\\!\left[\hat{\rho}_{0}\hat{S}_{z}\right]\simeq
0$, while preserving $\hat{S}^{2}$. It is more appropriate to think of the
finite size steady state as the average of the magnetization vector over a BTC
trajectory in the thermodynamic limit. This picture holds for any parameter
choice where the BTC phase exists and it is in perfect agreement with the
analysis of the semiclassical trajectories in Ref. [76].

However, it is important to stress that this structure emerges independently
in each eigenspace of the total spin $\hat{S}^{2}$. Since each of these
sectors of the Hilbert space has a different dimension, the full density
matrix is

$\hat{\rho}_{0}\simeq\sum_{n=0}^{N/2}\left[\frac{g(n)}{2n+1}\sum_{m=-n}^{n}\ket{n,m}\bra{n,m}+O(n^{-2})\right]\
,$ (12)

where $n$ and $m$ are the quantum numbers associated with $\hat{S}^{2}$ and
$\hat{S}_{z}$ respectively, and we assumed the number of spin variables $N$ to
be even. Above, $g(n)$ is the degeneracy of the corresponding subspace. Note
also that the corrections to the identity-like structure becomes more and more
relevant the smaller the total magnetization $n$ is. Hence, $\rho_{0}$ is not
at all close to the identity in the full Hilbert space but only if we look at
a fixed (large) value of of $n$, as we did in our analysis with $n=N/2$,
corresponding to the maximally polarized subsector.

In the ferromagnetic and mixed phases, the stationary-state density matrix
resembles, instead, a coherent spin state, with a degree of spin squeezing
[96] which depends on the values of the Hamiltonian parameters, $p$ and
$\delta\Gamma$ in particular. As an example, in Fig. 9(b) we show the absolute
value of $\hat{\rho}_{0}$ for $p=3$. It is evident that, differently from the
BTC steady state, the steady state also displays appreciable coherences. This
feature provides further interesting elements on the generality of the model
of Eq. (1) also with regard to the steady-state properties of the driven-
dissipative system, e.g., for subradiant and subradiance, coherence and spin
squeezing [97, 62].

To compare quantitatively the two cases, in Fig. 9(c) we plot the purity of
the density matrix $\mathrm{Tr}\\!\left[\hat{\rho}_{0}^{2}\right]$ as a
function of the system size $N$. In the BTC phase, corresponding to $p$ even,
the purity is practically the lowest possible $(N+1)^{-1}$, the same as the
identity matrix. $\hat{\rho}_{0}$ is indeed a mixed state without coherences
which satisfies
$\mathrm{Tr}\\!\left[\hat{\rho}_{0}\hat{S}_{x}\right]\simeq\mathrm{Tr}\\!\left[\hat{\rho}_{0}\hat{S}_{y}\right]\simeq\mathrm{Tr}\\!\left[\hat{\rho}_{0}\hat{S}_{z}\right]\simeq
0$. For odd values of $p$ the purity is much higher, showing that the steady-
state density matrix is quantitatively and qualitatively different in the two
phases.

Figure 9: (a) and (b): absolute value of the stationary-state density matrix
for $N=30$, $q=1$, $\omega_{x}=3\omega_{z}$ and $\delta\Gamma=0.2\omega_{z}$.
In (a) $p=2$, while in (b) $p=3$. The density matrix elements are computed in
the basis of the $z$-magnetization $\hat{S}_{z}=\frac{N}{2}\hat{J}_{z}$. (c):
purity of the stationary-state density matrix
$\mathrm{Tr}\\!\left[\rho_{0}^{2}\right]$ versus $1/(N+1)$ for several values
of $p$ and $q$. When the BTC phase is present, the purity is very close to the
lowest possible in the subspace of the Hilbert space where the dynamics is
constrained, namely $(N+1)^{-1}$. When the system relaxes towards a state with
finite $z$-magnetization, instead, the purity is much larger.

##  V BTCs existence condition 2: Collective dissipation

In this section, we show that only in the presence of collective dissipation
processes only, the dynamics of Eq. (3) displays time crystallinity. To this
purpose, we follow the dynamics of the simplest time crystal Hamiltonian
described in Eq. (5), i.e. when $p=2$ and $q=1$, coupled to a bath through a
string of jump operators $J^{\pm}_{s}=\sum_{i}^{N_{s}}\sigma_{i}^{\pm}/N_{s}$,
with $N_{s}\leq N$. Note that these Lindblad operators do not conserve the
total angular momentum $J^{2}$. The mean-field equations of motion read

$\left\\{\begin{aligned} &\dot{X}=2p\omega_{z}Z^{p-1}Y-2\left(\delta\Gamma
Z_{s}+\frac{\overline{\Gamma}}{N_{s}}\right)X_{s},\\\
&\dot{Y}=2XZ\left(q\omega_{x}X^{q-2}-p\omega_{z}Z^{p-2}\right)-2\left(\delta\Gamma
Z_{s}+\frac{\overline{\Gamma}}{N_{s}}\right)Y_{s},\\\
&\dot{Z}=-2q\omega_{x}YX^{q-1}+2\delta\Gamma\left(1-Z^{2}+\frac{1}{N_{s}}\right)-\frac{2\overline{\Gamma}}{N_{s}}Z_{s},\end{aligned}\right.$
(13)

where we have defined
$\overline{\Gamma}=\Gamma_{\uparrow}+\Gamma_{\downarrow}$ and
$\alpha_{s}=\braket{J_{\alpha}^{s}}$. The results in Sec. III are recovered in
the thermodynamic limit when $N_{s}=N$. In (a) of Fig. 10 we show the
trajectories of the magnetization $Z$ as a function of time for three
different values of $N_{s}=10,20,50$; with $\delta\Gamma=0.2\omega_{z}$ and
$\omega_{x}=1.1\omega_{z}$. By comparing with the collective dissipation
trajectory (gray dashed line) it emerges a twofold effect, both on the phase
and on the amplitude of the oscillations, due to the finite range of the
dissipation operators. The trajectories for different values of $N_{s}$ have
slightly different periods, leading to a progressive dephasing which, however,
has no relevant effect on the time-order. The more relevant consequence of the
finite range of the bath operators is the damping that reduces the amplitude
of the oscillations, eventually recovering a time-translational invariant
state.

To analyze the $N_{s}$ dependence of the time scale over which time order is
destroyed, in (b) of Fig. 10 we plot the amplitude of the oscillations of $Z$
as a function of time, for different $N_{s}$.

Figure 10: (a): Trajectories of the longitudinal magnetization $Z$ as a
function of time for $N_{s}=10,100,1000$, fixed $\delta\Gamma=0.2\omega_{z}$
and $\omega_{x}=1.1\omega_{z}$. The dashed line is the benchmark case with
global dissipation processes $N_{s}=N$. Reducing the value of $N_{s}$ the
trajectories are affected by dephasing and damping effects that are as
stronger as smaller $N_{s}$. (b): Oscillation amplitude of the $Z$ component
of the magnetization (full squares) as a function of time, for different
values of $N_{s}$, and $\delta\Gamma=0.2\omega_{z}$ and
$\omega_{x}=1.1\omega_{z}$. The data are fitted with functions parametrized as
$f(t)=\exp\left(-\beta t/N_{s}\right)$ (dashed lines), with the fit parameter
being $\beta=0.11$.

The dashed lines, obtained by fitting the data, behave as

$f(t)=\exp\left(-\beta t/N_{s}\right)\ ,$ (14)

with $\beta=0.11$. This means that the oscillations are exponentially damped
in time, with a rate which decreases as a power law of the length of the
string of the dissipation operators. Hence, only in the limit of global
dissipation processes, $N_{s}\to N$, the oscillations persist for arbitrarily
long time. This result confirms the intuition by Riera-Campeny and
collaborators [98] that global dissipation processes are a key ingredient to
observe the boundary time crystal. Our intuition is that this is an
application of a more general condition that requires at least a conserved
quantity in the dynamics. If so, in fact, the system cannot loose the
information on the initial state during the evolution and is prevented to
relax toward a time-independent state.

##  VI Physical interpretation

In our investigation, we relied only on numerical evidences to support our
claims, while rigorous proofs are left for future studies. However, we
conjecture some physical interpretations and connections with discrete time
crystals to help reaching a more complete understanding of the phenomenon.

Of the two conditions, the claim that the Hamiltonian should be
${\mathbb{Z}}_{2}$ invariant and the dissipation should break this symmetry is
the harder to physically justify. The role of the symmetry in the Hamiltonian
might be similar to that played in Floquet time crystals. There, the breaking
of discrete time translation usually comes from a periodic driving that at
each period connects different symmetry-broken sectors of the system’s
$\mathbb{Z}_{n}$ symmetric Hamiltonian [33]. Similarly, we observe that BTC
trajectories oscillate between two finite values of the magnetization, i.e.
between two states that break the Hamiltonian ${\mathbb{Z}}_{2}$ symmetry. If
there is indeed a connection with Floquet time crystals, this suggests that
the BTC phase can be observed when the Hamiltonian has a general discrete
symmetry.

A characteristic behavior for even $p$, as explained in Sec. IV, is that the
operators $\hat{J}^{+}$ and $\hat{J}^{-}$ play alternatively the role of an
incoherent dissipation and driving, depending on the sign of the
magnetization. This continuous switching of the roles of the Lindbald
operators concurs to stabilize the periodic orbits, which otherwise collapse
to a stationary state. However, since the presence of an additional symmetry
along the $x$ direction does not preclude the emergence of BTCs, our first
condition can probably be generalized as the condition for the Hamiltonian to
have a discrete symmetry group which is reduced to a smaller one when the
system is coupled with the environment.

From a more mathematical perspective, it is known that dissipation can lead to
nontrivial asymptotic states in the presence of degeneracies in the
Hamiltonian spectrum [99, 100, 101].

The second condition we address is that we need the jump operators to act
collectively on the whole system, i.e. to have the form
$J^{\pm}=\sum_{i}^{N}\sigma_{i}^{\pm}/N$. In this way, the Lindblad operators
commute with the total angular momentum $\hat{J}^{2}$, which is a strong
symmetry of the complete model and, consequently, a conserved quantity [66].
Thus, in the present paper, the condition of collective pumping/dissipation
and the conservation of $\hat{J}^{2}$ are used interchangeably.

In a more general framework, a conservation law reduces the effective
dimension of the space in which the dynamics takes place. In the BTC phase,
this reduction prevents the system to loose memory on its initial conditions
and to relax toward a stationary state, allowing the emergence of periodic
trajectories. We believe that the key element here is the presence of a
(quasi)conserved quantity during the open dynamics, independently of it being
also a strong symmetry of the Lindblad equation, thus leading to a weaker
constraint. This can be observed in spin models with power-law decaying
interactions, where $\hat{J}^{2}$ does not commute with the Hamiltonian but
the system might be still rigid enough to observe the emergence of BTCs.

##  VII Experimental realizations

In this section, we discuss how to experimentally implement the proposed
model. We discuss two main aspects: relevant quantum simulation platforms and
the engineering of higher-order spin interactions.

We have shown that the spin algebra can be implemented in ensembles of
spin-$\frac{1}{2}$ with all-to-all connectivity. Rydberg atoms [102], trapped
ions [103], artificial qubits in superconducting circuits [104], and color
defects in diamond [105] are the most prominent platforms in which such all-
to-all connectivity has been demonstrated. Let us point out that hybrid
quantum systems [2, 106], such as color vacancies in diamond coupled to a
single superconducting resonator offer the possibility to smooth out a long-
range photon-mediated interaction on several emitters [105].

With regards to tuning the order of interactions, $p$, in the Hamiltonian of
Eq. (1), the task can be divided in the standard cavity QED interactions,
$p,q\leq 2$ and in the more exotic many-body interactions, $p,q>2$. In
general, up to dipole approximation, the light-matter interaction that
characterizes cavity QED systems involves only two particles per time, i.e.
annihilation/creation of bosonic/fermionic particles. The interaction mediated
by bosons can in several cases be traced out, leading to effective spin-spin
interactions ($p=2$). Spin-spin interactions are also the outcome of several
microscopic processes, such as genuine magnetic interaction, internal-external
degrees of freedom coupling in ions, Van der Waals forces in neutral atoms and
so on. Using the quantum simulation paradigm [107], quantum phase transitions
and dissipative phase transitions with long-range interactions have been
engineered in trapped-ion and Rydberg-atom simulators [103, 108, 109, 110],
although a power-law decay in the interaction length is common. In artificial
atoms such as superconducting qubits, effective spin-spin interactions have
been first proposed [111] and more recently implemented [104, 112],
reproducing the physics of paradigmatic spin models, the Ising and Lipkin-
Meshkov-Glick models.

The couplings for $p,q>2$ generally represent an indirect process and have
been theoretically proposed in the context of spin glasses [113], which have
become recently amenable to quantum simulation [114, 115, 116]. Direct
interactions show up in the non-perturbative coupling between light and matter
in the so-called ultrastrong and deep-strong coupling regimes [117], where it
is not possible to arbitrarily tune higher-order coupling strengths, and where
they occur in the boson-fermion framework. Other promising candidates to
realize three-body interactions are hard-core bosons[118] or polar molecules
in optical lattices [119, 120, 121], and trapped Rydberg ions [122, 123],
which represent the state-of-the-art for controlling complex many-body
systems.

Recently, the multiple excited states of natural and artificial atoms have
been exploited to implement quantum systems beyond qubit systems, but with
qudits or in general bosonic degrees of freedom, e.g., molecules or transmons
in superconducting devices. It is important to point out that spin models can
be mapped to bosons (and viceversa), e.g., with a Holstein-Primakoff
approximation, but with the limitation of being in the diluted regime, i.e.
where nonlinear effects are negligible [52].

With this regard, higher order interactions in bosonic systems based on
superconducting circuit devices have recently been implemented, providing more
freedom to tune higher-order interactions and suppress lower-order ones,
exploiting nonlinearities in Josephson junctions [124, 125, 126, 127], which
can also be used to simulate spins with angular momentum greater than
$\frac{1}{2}$ [128, 129].

To conclude this section, it is worth to mention that the all-to-all
connectivity condition can probably be relaxed. Restricting ourselves to
$p,q\leq 2$, it is a well established fact [130, 24] that a long-range Ising
model, with interactions between sites $i$ and $j$ decaying with distance as
$J_{i,j}\propto|i-j|^{-\alpha}$, belongs to the same universality class of the
fully-connected model considered in this paper when $\alpha<1$. We believe
that the “rigidity” of the collective spin dynamics arising from the long
range interactions is sufficient to observe boundary time crystals in such
systems, at least for sufficiently small values of $\alpha$. It is hard to
anticipate, though, how the long time behavior changes if $\alpha\geq 1$ when
fluctuations between different magnetization sectors may destroy the BTC
phase.

##  VIII Discussion and conclusion

In this work, we discussed two conditions for the existence of boundary time-
crystals in a large class of generalized $p$-interaction spin models with
dissipation, showing that the emergence of this non-equilibrium phase is
strongly related to the symmetries and the conserved quantities of the
Hamiltonian and how they are affected by the dissipation.

The first condition is that BTCs arise only if the spin Hamiltonian is
${\mathbb{Z}}_{2}$-symmetric and the Lindblad operators explicitly break this
symmetry. If so, the mean-field trajectories of the magnetization display
persistent oscillations in time. This result resolves the apparent
contradiction between Ref. [40] and Ref. [73]: in the latter the symmetry
condition is not satisfied, hence the absence of BTCs.

This reflects an emerging U$(1)$ symmetry of the equation of motion: closed
orbits are invariant under rotations along the trajectory itself, meaning that
the system preserves only partial information on the initial condition since
all points on the same orbit lead to the same time-crystal dynamics. At
finite-size $N$, instead, the oscillations are affected by a damping rate that
decreases by increasing $N$. The U$(1)$ symmetry of the thermodynamic limit
now seems to appear as an approximated discrete ${\mathbb{Z}}_{N+1}$ symmetry,
similar to that argued in Ref. [95], as suggested by the characteristic
profile of the stationary-state density matrix $\hat{\rho}_{0}$ associated to
the BTC phase. This condition provides a nice parallelism between boundary and
discrete time crystals, a phase in which the subharmonic response is due to
the exploration of different subsectors of some symmetric Hamiltonian or
evolution operator [39, 37, 38, 43].

The second condition regards the operators coupling the system with the bath:
the time-crystal order is destroyed by dissipation processes that do not
conserve the total angular momentum. In particular, we showed that if only a
portion $N_{s}<N$ of the system is coupled to the same external bath, the
amplitude of oscillations decreases exponentially with $N_{s}$.

This result suggests that, more generally, BTCs arise only when there is at
least one strong symmetry in the dynamics that prevents the system to lose
information on the initial conditions and to attain a time-independent steady
state. This draws an interesting parallelism with results that have shown how
the symmetry sector to which an initial state belongs to can determine the
current flows and quantum transport properties of the steady state [131].
However, we do not exclude that other kind of long range couplings between the
system and the environment may still allow for a stable BTC phase. This might
be an interesting route to pursue to better understand if the robustness of
the BTC phase is determined by the correlation length of the jump operators or
strictly requires a homogeneous collective bath coupling.

Local dissipation induces also a dephasing in the mean-field trajectories, on
top of the usual damping. This dephasing is probably due to the exploration of
sectors with different total angular momentum; in the Dicke representation of
permutational-invariant systems [52], this is equivalent to jumps between
different Dicke ladders. Since different sectors have different Liouvillian
spectra, the difference between the imaginary parts of the eigenvalues might
indeed be responsible for the progressive dephasing of the oscillations.
Understanding this mechanism could be an important step toward the full uptake
of the BTC phenomenon. Moreover, it could give some hints on possible
configurations in which the competition between global and local dissipation
processes could give rise to interesting non-trivial phases of matters.

It is important to observe that the two conditions must be met in order to
enter a BTC phase in the thermodynamic limit. A different, but similar
phenomenology is that of decoherence-free subspaces [64, 65, 66, 70], where
time-dependent oscillations for long times are present at any system size, as
investigated in the context of artificial giant atoms [64, 67, 71]. An example
for the considered model is obtained for $p=0$, hence
$\hat{H}\propto\hat{J}_{x}^{q}$ and a Lindblad jump operator equal to
$\hat{J}_{x}$ (which could be interpreted as a stochastic jump or global
dephasing), or another jump operator commuting with the Hamiltonian, such as
spin-squeezing dissipation, $\hat{J}_{x}^{2}$, or a collective depolarizing
channel.

These models preserve the ${\mathbb{Z}}_{2}$ symmetry but are characterized by
persistent oscillations at any system size $N$. Indeed, the Liouvillian
spectrum at finite $N$ is qualitatively independent from the coupling strength
with the bath, indicating that collective dephasing does not induce different
phases in our model and acts somehow “trivially” on the system.

All the results of this work have been derived for a general class of
dissipative spin models, which can be implemented in fully connected
spin-$\frac{1}{2}$ ensembles, with a specific choice of the Lindblad
dissipators, which can be engineered in quantum simulators, e.g., in
superconducting-circuit-based quantum devices. However, the fact that BTCs are
clearly associated with symmetry properties and not with the specific form of
the Hamiltonian suggest that these criteria might apply for a wider class of
systems and Lindblad operators. We anticipate that one of the most interesting
directions for future investigations would be to explore the universality of
our results, aiming for a general formulation of the existence conditions of
BTCs in terms of symmetries and conserved quantities that may apply to a
generic quantum many-body system.

###### Acknowledgements.

The authors acknowledge useful discussions with M. Burrello, M. Dalmonte, R.
Fazio, F. Minganti, D. Rossini, G.E. Santoro, A. Silva and F.M. Surace. G.P.
thanks Giovanni Piccitto for useful comments. F.N. is supported in part by:
Nippon Telegraph and Telephone Corporation (NTT) Research, the Japan Science
and Technology Agency (JST) [via the Quantum Leap Flagship Program (Q-LEAP),
the Moonshot R&D Grant Number JPMJMS2061, and the Centers of Research
Excellence in Science and Technology (CREST) Grant No. JPMJCR1676], the Japan
Society for the Promotion of Science (JSPS) [via the Grants-in-Aid for
Scientific Research (KAKENHI) Grant No. JP20H00134 and the JSPS–RFBR Grant No.
JPJSBP120194828], the Army Research Office (ARO) (Grant No. W911NF-18-1-0358),
the Asian Office of Aerospace Research and Development (AOARD) (via Grant No.
FA2386-20-1-4069), and the Foundational Questions Institute Fund (FQXi) via
Grant No. FQXi-IAF19-06. M.W. is supported by the Villum Foundation (Research
Grant No. 25310). This project has received funding from the European Union’s
Horizon 2020 research and innovation program under the Marie Sklodowska-Curie
grant agreement No. 847523 ‘INTERACTIONS’

##  Appendix A Derivation of the equations of motion

In this section we derive the mean-field equations of motion of Eq. (13). Let
us assume to have a dynamics described by the master equation

$\displaystyle\dot{\rho}=-$ $\displaystyle
i[\hat{H},\rho]+N\Gamma^{\uparrow}\left(J^{+}_{s}\rho
J^{-}_{s}-\frac{1}{2}\\{J^{-}_{s}J^{+}_{s},\rho\\}\right)$ (15)
$\displaystyle+$ $\displaystyle N\Gamma^{\downarrow}\left(J^{-}_{s}\rho
J^{+}_{s}-\frac{1}{2}\\{J^{+}_{s}J^{-}_{s},\rho\\}\right)$ $\displaystyle=$
$\displaystyle\dot{\rho}^{c}+\dot{\rho}^{d}.$

This evolution has two different contributions, a term that accounts for the
coherent evolution, $\dot{\rho}^{c}$, and a term that accounts for the
incoherent one, $\dot{\rho}^{d}$. The coherent part is simply obtained by
evaluating the commutator with the Hamiltonian leading to

$\displaystyle\dot{J}_{x}^{c}=2\sum_{n=0}^{p-1}J_{z}^{p-1-n}J_{y}J_{z}^{n},$
(16)
$\displaystyle\dot{J}_{y}^{c}=-2\sum_{n=0}^{p-1}J_{z}^{p-1-n}J_{x}J_{z}^{n}+2gJ_{z},$
$\displaystyle\dot{J}_{z}^{c}=-2gJ_{y}.$

We want to evaluate the incoherent contribution to the equations of motion for
the expectation values of the spin operators. This is given by
$\braket{\dot{J^{\alpha}}}=\text{Tr}\left(\dot{\rho}^{d}J^{\alpha}\right)$
leading to

$\displaystyle\braket{\dot{J_{\alpha}}}=-$ $\displaystyle
N\Gamma^{\uparrow}\text{Tr}\left(J^{+}_{s}\rho
J^{-}_{s}J_{\alpha}-\frac{1}{2}\\{J^{-}_{s}J^{+}_{s},\rho\\}J_{\alpha}\right)$
(17) $\displaystyle+$ $\displaystyle
N\Gamma^{\downarrow}\text{Tr}\left(J^{-}_{s}\rho
J^{+}_{s}J_{\alpha}-\frac{1}{2}\\{J^{+}_{s}J^{-}_{s},\rho\\}J_{\alpha}\right).$

By manipulating the equations we obtain

$\displaystyle\braket{\dot{J_{\alpha}}}=$
$\displaystyle\frac{Ni\delta\Gamma}{2}\braket{\left\\{J^{x}_{s},\left[J_{\alpha},J_{s}^{y}\right]\right\\}-\left\\{J^{y}_{s},\left[J_{\alpha},J_{s}^{x}\right]\right\\}}$
(18) $\displaystyle+$
$\displaystyle\frac{N\overline{\Gamma}}{2}\braket{\left[J^{x}_{s},\left[J_{\alpha},J_{s}^{x}\right]\right]+\left[J^{y}_{s},\left[J_{\alpha},J_{s}^{y}\right]\right]}$
$\displaystyle=$
$\displaystyle\frac{i\delta\Gamma}{2N_{s}^{2}}\sum_{i,j}^{N_{s}}\sum_{l}^{N}\braket{\left\\{\sigma^{x}_{i},\left[\sigma^{\alpha}_{l},\sigma_{j}^{y}\right]\right\\}-\left\\{\sigma^{y}_{i},\left[\sigma^{\alpha}_{l},\sigma_{j}^{x}\right]\right\\}}$
$\displaystyle+$
$\displaystyle\frac{\overline{\Gamma}}{2N_{s}^{2}}\sum_{i,j}^{N_{s}}\sum_{l}^{N}\braket{\left[\sigma^{x}_{i},\left[\sigma^{\alpha}_{l},\sigma_{j}^{x}\right]\right]+\left[\sigma^{y}_{i},\left[\sigma^{\alpha}_{l},\sigma_{j}^{y}\right]\right]}$
$\displaystyle=-$
$\displaystyle\frac{\delta\Gamma}{2N_{s}^{2}}\sum_{i,j}^{N_{s}}\braket{\left\\{\sigma^{x}_{i},\epsilon_{\alpha
y\beta}\sigma^{\beta}_{j}\right\\}-\left\\{\sigma^{y}_{i},\epsilon_{\alpha
x\beta}\sigma^{\beta}_{j}\right\\}}$ $\displaystyle-$
$\displaystyle\frac{\overline{\Gamma}}{2N_{s}^{2}}\sum_{i,j}^{N_{s}}\braket{\left[\sigma^{x}_{i},\epsilon_{\alpha
x\beta}\sigma^{\beta}_{j}\right]+\left[\sigma^{y}_{i},\epsilon_{\alpha
y\beta}\sigma^{\beta}_{j}\right]},$

where $\delta\Gamma=\Gamma^{\uparrow}-\Gamma^{\downarrow}$ and
$\overline{\Gamma}=\Gamma^{\uparrow}+\Gamma^{\downarrow}$.

##  Appendix B Derivation of the phase portrait

In this appendix we provide some details on the derivation of the phase
portrait. In general, the spin operators in the semi-classical approximation
can be parametrized on the Bloch sphere by giving the radius $r$ of the
sphere, the polar angle $\varphi$ and the azimuthal one $\theta$. This mapping
is singular in $\theta=0,\pi$, therefore some care is needed in choosing the
axis defining $\varphi$ and $\theta$. The results of the paper have been
derived by parametrizing the spin as

$\left(X=r\sin\theta\cos\varphi,\ Y=r\sin\theta\sin\varphi,\
Z=r\cos\theta\right).$ (19)

The only exceptions are the results in Subsec. III.3, in which we assumed the
polar angle to span the $yz$ plane (namely, we fixed $X=\cos\theta$). We will
omit the derivation of this case since it can be obtained by simply extending
the following calculations.

Eq. (19) can be inverted leading to

$\displaystyle r=\sqrt{X^{2}+Y^{2}+Z^{2}},$ (20)
$\displaystyle\varphi=\arctan(Y/X),$
$\displaystyle\theta=\arctan((X^{2}+Y^{2})/Z).$

By deriving and substituting the expressions of $X,Y,Z$ in terms of
$r,\varphi,\theta$ we obtain, for the collective decay, the following
equations of motion (since classically we have that $\varphi$ and $\cos\theta$
are conjugate variables, it is better to derive the equations for
$\cos(\theta)$ instead of $\theta$)

$\displaystyle\dot{r}=2r\delta\Gamma\cos\theta\left(1-\ r^{2}\right),$ (21)
$\displaystyle\dot{\varphi}=-2\omega_{z}p\cos^{p-1}\theta+2\omega_{x}q\cos\theta\sin^{q-2}\theta\cos^{q}\varphi,$
$\displaystyle\dot{\cos\theta}=-\left(2\omega_{x}\sin^{q}\theta\cos^{q-1}\varphi\sin\varphi+2\delta\Gamma\left(1-\cos^{2}\theta\right)\right)/r.$

We notice that $r=1$ is a fixed point, hence the dynamics is constrained on
the surface of the Bloch sphere.

##  Appendix C Mean-field description of the totally symmetric subsector

Let us consider a generic single spin density matrix, parametrized as

$\hat{\rho}=\left[\begin{array}[]{cc}a&b{\rm e}^{-i\phi}\\\ b{\rm
e}^{i\phi}&1-a\end{array}\right]$ (22)

where $a$ and $b$ are both real numbers. Its purity is defined as
$\mathcal{P}=\mathrm{Tr}\\!\left[\hat{\rho}^{2}\right]=\frac{1}{2}\left[(2a-1)^{2}+1+4b^{2}\right]\leq
1$. The reason why write it in this way will be clearer in the following.

We can construct a mean-field ansatz for the system by taking the tensor
product of $N$ identical spins described by density matrices as Eq. (22)

$\hat{\rho}_{\rm mf}=\bigotimes_{j=1}^{N}\hat{\rho}_{j}\ .$ (23)

We can compute the expectation value of the total spin
$\hat{S}^{2}=\hat{S}_{x}^{2}+\hat{S}_{y}^{2}+\hat{S}_{z}^{2}$ on this state.
First, let us rewrite it in terms of single-spin Pauli matrices

$\hat{S}^{2}=\frac{3}{4}N+\frac{1}{4}\sum_{k,j\neq
k}\hat{\sigma}^{x}_{j}\hat{\sigma}^{x}_{k}+\hat{\sigma}^{y}_{j}\hat{\sigma}^{y}_{k}+\hat{\sigma}^{z}_{j}\hat{\sigma}^{z}_{k}\
.$ (24)

Then, simple algebra leads to

$\mathrm{Tr}\\!\left[\hat{\rho}_{\rm
mf}\hat{S}^{2}\right]=\frac{3}{4}N+\frac{1}{4}N(N-1)\left[(2a-1)^{2}+4b^{2}\right]\
,$ (25)

which can conveniently be rewritten in terms of the purity of the single-spin
density matrix

$\mathrm{Tr}\\!\left[\hat{\rho}_{\rm
mf}\hat{S}^{2}\right]=\frac{3}{4}N+\frac{1}{4}N(N-1)(2\mathcal{P}-1)\ .$ (26)

If $\hat{S}^{2}$ is conserved and we consider in the maximally-polarized
subsector, the total spin has to be equal to
$\frac{N}{2}\left(\frac{N}{2}+1\right)$, which is compatible with Eq. (26)
only if $\mathcal{P}=1$. This means that a mean field ansatz for the density
matrix cannot describe a mixed state in the maximally polarized subspace, but
only pure ones. That is why to recover a mean-field approximation of mixed
states such as that shown in Fig. 9, it is necessary to average over all
magnetization vectors belonging to the same time-crystal trajectory.

## References

  * MacFarlane  _et al._ [2003] A. G. J. MacFarlane, J. P. Dowling, and G. J. Milburn, Quantum technology: the second quantum revolution, Philos. T. R. Soc. A 361, 1655 (2003).
  * Xiang  _et al._ [2013] Z.-L. Xiang, S. Ashhab, J. Q. You, and F. Nori, Hybrid quantum circuits: Superconducting circuits interacting with other quantum systems, Rev. Mod. Phys. 85, 623 (2013).
  * Kurizki  _et al._ [2015] G. Kurizki, P. Bertet, Y. Kubo, K. Mølmer, D. Petrosyan, P. Rabl, and J. Schmiedmayer, Quantum technologies with hybrid systems, Proc. Natl. Acad. Sci. USA 112, 3866 (2015).
  * Acín  _et al._ [2018] A. Acín, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, D. Esteve, N. Gisin, S. J. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. F. Riedel, P. O. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. K. Wilhelm, The quantum technologies roadmap: a European community view, New J. Phys. 20, 080201 (2018).
  * Babadi  _et al._ [2015] M. Babadi, E. Demler, and M. Knap, Far-from-equilibrium field theory of many-body quantum spin systems: Prethermalization and relaxation of spin spiral states in three dimensions, Phys. Rev. X 5 (2015).
  * Kollar  _et al._ [2011] M. Kollar, F. A. Wolf, and M. Eckstein, Generalized Gibbs ensemble prediction of prethermalization plateaus and their relation to nonthermal steady states in integrable systems, Phys. Rev. B 84 (2011).
  * Stark and Kollar [2013] M. Stark and M. Kollar, Kinetic description of thermalization dynamics in weakly interacting quantum systems (2013), arXiv:1308.1610 [cond-mat.str-el] .
  * Fratus and Srednicki [2016] K. R. Fratus and M. Srednicki, Eigenstate thermalization and spontaneous symmetry breaking in the one-dimensional transverse-field Ising model with power-law interactions (2016), arXiv:1611.03992 [cond-mat.stat-mech] .
  * Safavi-Naini  _et al._ [2019] A. Safavi-Naini, M. L. Wall, O. L. Acevedo, A. M. Rey, and R. M. Nandkishore, Quantum dynamics of disordered spin chains with power-law interactions, Phys. Rev. A 99 (2019).
  * Lerose  _et al._ [2019] A. Lerose, B. Žunkovič, A. Silva, and A. Gambassi, Quasilocalized excitations induced by long-range interactions in translationally invariant quantum spin chains, Phys. Rev. B 99 (2019).
  * Pappalardi  _et al._ [2018] S. Pappalardi, A. Russomanno, B. Žunkovič, F. Iemini, A. Silva, and R. Fazio, Scrambling and entanglement spreading in long-range spin chains, Phys. Rev. B 98 (2018).
  * Koffel  _et al._ [2012] T. Koffel, M. Lewenstein, and L. Tagliacozzo, Entanglement entropy for the long-range Ising chain in a transverse field, Phys. Rev. Lett. 109 (2012).
  * Gong  _et al._ [2017] Z.-X. Gong, M. Foss-Feig, F. G. Brandão, and A. V. Gorshkov, Entanglement area laws for long-range interacting systems, Phys. Rev. Lett. 119 (2017).
  * Piccitto and Silva [2019a] G. Piccitto and A. Silva, Dynamical phase transition in the transverse field Ising chain characterized by the transverse magnetization spectral function, Phys. Rev. B 100 (2019a).
  * Kastner [2017] M. Kastner, N-scaling of timescales in long-rangeN-body quantum systems, J. Stat. Mech. 2017, 014003 (2017).
  * Hauke and Tagliacozzo [2013] P. Hauke and L. Tagliacozzo, Spread of correlations in long-range interacting quantum systems, Phys. Rev. Lett. 111, 207202 (2013).
  * Cevolani  _et al._ [2015] L. Cevolani, G. Carleo, and L. Sanchez-Palencia, Protected quasilocality in quantum systems with long-range interactions, Phys. Rev. A 92, 041603 (2015).
  * Lang  _et al._ [2018] J. Lang, B. Frank, and J. C. Halimeh, Dynamical quantum phase transitions: A geometric picture, Phys. Rev. Lett. 121 (2018).
  * Weidinger  _et al._ [2017] S. A. Weidinger, M. Heyl, A. Silva, and M. Knap, Dynamical quantum phase transitions in systems with continuous symmetry breaking, Phys. Rev. B 96 (2017).
  * Sandri  _et al._ [2012] M. Sandri, M. Schiró, and M. Fabrizio, Linear ramps of interaction in the fermionic Hubbard model, Phys. Rev. B 86 (2012).
  * Gurarie [2009] V. Gurarie, Nonequilibrium dynamics of weakly and strongly paired superconductors, Phys. Rev. Lett. 103 (2009).
  * Sciolla and Biroli [2010] B. Sciolla and G. Biroli, Quantum quenches and off-equilibrium dynamical transition in the infinite-dimensional Bose-Hubbard model, Phys. Rev. Lett. 105 (2010).
  * Piccitto  _et al._ [2019] G. Piccitto, B. Žunkovič, and A. Silva, Dynamical phase diagram of a quantum Ising chain with long-range interactions, Phys. Rev. B 100 (2019).
  * Piccitto and Silva [2019b] G. Piccitto and A. Silva, Crossover from fast to slow dynamics in a long range interacting Ising chain, J. Stat. Mech. 2019, 094017 (2019b).
  * Langen  _et al._ [2013] T. Langen, R. Geiger, M. Kuhnert, B. Rauer, and J. Schmiedmayer, Local emergence of thermal correlations in an isolated quantum many-body system, Nature Phys. 9, 640–643 (2013).
  * Rosch  _et al._ [2008] A. Rosch, D. Rasch, B. Binz, and M. Vojta, Metastable superfluidity of repulsive fermionic atoms in optical lattices, Phys. Rev. Lett. 101 (2008).
  * Berges  _et al._ [2004] J. Berges, S. Borsanyi, and C. Wetterich, Prethermalization, Phys. Rev. Lett. 93, 142002 (2004).
  * Sieberer  _et al._ [2013] L. M. Sieberer, S. D. Huber, E. Altman, and S. Diehl, Dynamical critical phenomena in driven-dissipative systems, Phys. Rev. Lett. 110, 195301 (2013).
  * Sieberer  _et al._ [2016] L. M. Sieberer, M. Buchhold, and S. Diehl, Keldysh field theory for driven open quantum systems, Rep. Prog. Phys. 79, 096001 (2016).
  * Altman  _et al._ [2015] E. Altman, L. M. Sieberer, L. Chen, S. Diehl, and J. Toner, Two-dimensional superfluidity of exciton polaritons requires strong anisotropy, Phys. Rev. X 5, 011017 (2015).
  * Maghrebi and Gorshkov [2016] M. F. Maghrebi and A. V. Gorshkov, Nonequilibrium many-body steady states via Keldysh formalism, Phys. Rev. B 93, 014307 (2016).
  * Wilczek [2012] F. Wilczek, Quantum time crystals, Physical Review Letters 109 (2012).
  * 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, Rep. Prog. Phys. 81, 016401 (2017).
  * Watanabe and Oshikawa [2015] H. Watanabe and M. Oshikawa, Absence of quantum time crystals, Phys. Rev. Lett. 114 (2015).
  * Kozin and Kyriienko [2019] V. K. Kozin and O. Kyriienko, Quantum time crystals from Hamiltonians with long-range interactions, Phys. Rev. Lett. 123, 210602 (2019).
  * Russomanno  _et al._ [2017a] A. Russomanno, F. Iemini, M. Dalmonte, and R. Fazio, Floquet time crystal in the Lipkin-Meshkov-Glick model, Phys, Rev. B 95 (2017a).
  * Else  _et al._ [2016] D. V. Else, B. Bauer, and C. Nayak, Floquet time crystals, Phys. Rev. Lett. 117 (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).
  * Iemini  _et al._ [2018] F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte, and R. Fazio, Boundary time crystals, Phys. Rev. Lett. 121, 035301 (2018).
  * Riera-Campeny  _et al._ [2020a] A. Riera-Campeny, M. Moreno-Cardoner, and A. Sanpera, Time crystallinity in open quantum systems, Quantum 4, 270 (2020a).
  * Hurtado-Gutiérrez  _et al._ [2020] R. Hurtado-Gutiérrez, F. Carollo, C. Pérez-Espigares, and P. I. Hurtado, Building continuous time crystals from rare events, Phys. Rev. Lett. 125, 160601 (2020).
  * Lazarides  _et al._ [2020] A. Lazarides, S. Roy, F. Piazza, and R. Moessner, Time crystallinity in dissipative Floquet systems, Phys. Rev. Research 2 (2020).
  * Sciolla and Biroli [2013] B. Sciolla and G. Biroli, Quantum quenches, dynamical transitions, and off-equilibrium quantum criticality, Phys. Rev. B 88 (2013).
  * Lerose and Pappalardi [2020a] A. Lerose and S. Pappalardi, Origin of the slow growth of entanglement entropy in long-range interacting spin systems, Phys. Rev. Research 2 (2020a).
  * Lerose and Pappalardi [2020b] A. Lerose and S. Pappalardi, Bridging entanglement dynamics and chaos in semiclassical systems (2020b), arXiv:2005.03670 [quant-ph] .
  * Lee  _et al._ [2013] T. E. Lee, S. Gopalakrishnan, and M. D. Lukin, Unconventional magnetism via optical pumping of interacting spin systems, Phys. Rev. Lett. 110, 257204 (2013).
  * Biella  _et al._ [2015] A. Biella, L. Mazza, I. Carusotto, D. Rossini, and R. Fazio, Photon transport in a dissipative chain of nonlinear cavities, Phys. Rev. A 91, 053815 (2015).
  * Biella  _et al._ [2017] A. Biella, F. Storme, J. Lebreuilly, D. Rossini, R. Fazio, I. Carusotto, and C. Ciuti, Phase diagram of incoherently driven strongly correlated photonic lattices, Phys. Rev. A 96, 023839 (2017).
  * Biella  _et al._ [2018] A. Biella, J. Jin, O. Viyuela, C. Ciuti, R. Fazio, and D. Rossini, Linked cluster expansions for open quantum systems on a lattice, Phys. Rev. B 97, 035103 (2018).
  * Minganti  _et al._ [2016] F. Minganti, N. Bartolo, J. Lolli, W. Casteels, and C. Ciuti, Exact results for schrödinger cats in driven-dissipative systems and their feedback control, Sci. Rep. 6, 26987 EP (2016).
  * Shammah  _et al._ [2018] N. Shammah, S. Ahmed, N. Lambert, S. De Liberato, and F. Nori, Open quantum systems with local and collective incoherent processes: Efficient numerical simulations using permutational invariance, Phys. Rev. A 98, 063815 (2018).
  * Huybrechts and Wouters [2019] D. Huybrechts and M. Wouters, Cluster methods for the description of a driven-dissipative spin model, Phys. Rev. A 99, 043841 (2019).
  * Huybrechts  _et al._ [2020] D. Huybrechts, F. Minganti, F. Nori, M. Wouters, and N. Shammah, Validity of mean-field theory in a dissipative critical system: Liouvillian gap, $\mathbb{PT}$-symmetric antigap, and permutational symmetry in the $\mathit{XYZ}$ model, Phys. Rev. B 101, 214302 (2020).
  * Wauters  _et al._ [2020] M. M. Wauters, G. B. Mbeng, and G. E. Santoro, Polynomial scaling of the quantum approximate optimization algorithm for ground-state preparation of the fully connected $p$-spin ferromagnet in a transverse field, Phys. Rev. A 102, 062404 (2020).
  * Brandes [2005] T. Brandes, Coherent and collective quantum optical effects in mesoscopic systems, Phys. Rep. 408, 315–474 (2005).
  * Dylewsky  _et al._ [2016] D. Dylewsky, J. K. Freericks, M. L. Wall, A. M. Rey, and M. Foss-Feig, Nonperturbative calculation of phonon effects on spin squeezing, Phys. Rev. A 93 (2016).
  * Matsunaga  _et al._ [2014] R. Matsunaga, N. Tsuji, H. Fujita, A. Sugioka, K. Makise, Y. Uzawa, H. Tera, Z. Wang, H. Aoki, and R. Shimano, Light-induced collective pseudospin precession resonating with Higgs mode in a superconductor, Science 345 (2014).
  * Bapst and Semerjian [2012a] V. Bapst and G. Semerjian, On quantum mean-field models and their quantum annealing, J. Stat. Mech. 2012, P06007 (2012a).
  * Vodola  _et al._ [2015] D. Vodola, L. Lepori, E. Ercolessi, and G. Pupillo, Long-range Ising and Kitaev models: phases, correlations and edge modes, New J. Phys. 18, 015001 (2015).
  * Carollo and Lesanovsky [2020] F. Carollo and I. Lesanovsky, Proof of a nonequilibrium pattern-recognition phase transition in open quantum multimode Dicke models (2020), arXiv:2009.13932 [cond-mat.stat-mech] .
  * Sánchez Muñoz  _et al._ [2019] C. Sánchez Muñoz, B. Buča, J. Tindall, A. González-Tudela, D. Jaksch, and D. Porras, Symmetries and conservation laws in quantum trajectories: Dissipative freezing, Phys. Rev. A 100, 042113 (2019).
  * Tindall  _et al._ [2020] J. Tindall, C. Sánchez Muñoz, B. Buča, and D. Jaksch, Quantum synchronisation enabled by dynamical symmetries and dissipation, New J. Phys. 22, 013026 (2020).
  * Baumgartner and Narnhofer [2008] B. Baumgartner and H. Narnhofer, Analysis of quantum semigroups with GKS–lindblad generators: II. general, J. Phys. A 41, 395303 (2008).
  * Buča and Prosen [2012] B. Buča and T. Prosen, A note on symmetry reductions of the Lindblad equation: transport in constrained open spin chains, New J. Phys. 14, 073007 (2012).
  * Albert and Jiang [2014] V. V. Albert and L. Jiang, Symmetries and conserved quantities in lindblad master equations, Phys. Rev. A 89, 022118 (2014).
  * Kockum  _et al._ [2018] A. F. Kockum, G. Johansson, and F. Nori, Decoherence-free interaction between giant atoms in waveguide quantum electrodynamics, Phys. Rev. Lett. 120 (2018).
  * Manzano and Hurtado [2018] D. Manzano and P. Hurtado, Harnessing symmetry to control quantum transport, Advances in Physics 67, 1 (2018).
  * Minganti  _et al._ [2020] F. Minganti, I. I. Arkhipov, A. Miranowicz, and F. Nori, Correspondence between dissipative phase transitions of light and time crystals (2020), arXiv:2008.08075 [quant-ph] .
  * Lieu  _et al._ [2020] S. Lieu, R. Belyansky, J. T. Young, R. Lundgren, V. V. Albert, and A. V. Gorshkov, Symmetry breaking and error correction in open quantum systems, Physical Review Letters 125, 10.1103/physrevlett.125.240405 (2020).
  * Buča  _et al._ [2019] B. Buča, J. Tindall, and D. Jaksch, Non-stationary coherent quantum many-body dynamics through dissipation, Nat. Commun. 10, 1730 (2019).
  * [72] Throughout the paper, the word “collective” identifies an operator which is the sum over all sites of single spin objects, $\hat{O}=\sum_{j=1}^{N}\hat{o}_{j}$. Notice that $\hat{O}$ is invariant under permutation of the spin indices.
  * Wang and Fazio [2020] P. Wang and R. Fazio, Dissipative phase transitions in the fully-connected Ising model with $p$-spin interaction (2020), arXiv:2008.10045 [cond-mat.quant-gas] .
  * Manzano [2020] D. Manzano, A short introduction to the lindblad master equation, AIP Advances 10, 025106 (2020).
  * Am-Shallem  _et al._ [2015] M. Am-Shallem, A. Levy, I. Schaefer, and R. Kosloff, Three approaches for representing lindblad dynamics by a matrix-vector notation (2015), arXiv:1510.08634 [quant-ph] .
  * Seibold  _et al._ [2020] K. Seibold, R. Rota, and V. Savona, Dissipative time crystal in an asymmetric nonlinear photonic dimer, Phys. Rev. A 101 (2020).
  * Gardner [1985] E. Gardner, Spin glasses with p-spin interactions, Nucl. Phys. B 257, 747 (1985).
  * Filippone  _et al._ [2011] M. Filippone, S. Dusuel, and J. Vidal, Quantum phase transitions in fully connected spin models: An entanglement perspective, Phys. Rev. A 83, 022327 (2011).
  * Bapst and Semerjian [2012b] V. Bapst and G. Semerjian, On quantum mean-field models and their quantum annealing, J. Stat. Mech. 2012, P06007 (2012b).
  * Wauters  _et al._ [2017] M. M. Wauters, R. Fazio, H. Nishimori, and G. E. Santoro, Direct comparison of quantum and simulated annealing on a fully connected Ising ferromagnet, Phys. Rev. A 96, 022326 (2017).
  * Englert and Morigi [2002] B.-G. Englert and G. Morigi, Five lectures on dissipative master equations, Lecture Notes in Physics , 55–106 (2002).
  * Note [1] We want to remark that Noether theorem can not be generalized at the operator level in the case of open quantum systems. For this reason there is no one-to-one correspondence between symmetries and conserved quantities. In particular, it is possible to find conserved quantities that are not associated to any symmetry, but any strong symmetry in the model has an associated conserved quantity[66]..
  * Note [2] The presence of purely imaginary Lyapunov exponents is necessary for the existence of the BTCs but it is not sufficient. For instance, it is possible to have purely imaginary Lyapunov exponents associated to slow spiraling trajectories in which the damping of the oscillations is a non-linear effect of the dynamics.
  * Note [3] Notice that in the BTC phase every initial condition is associated to a different periodic orbit that keeps memory of the evolution of the initial state. This a substantial difference with the dynamics in presence of limit cycles in which all the trajectories eventually reach the same periodic orbit independently on the initial conditions.
  * Patra  _et al._ [2020] A. Patra, B. L. Altshuler, and E. A. Yuzbashyan, Driven-dissipative dynamics of atomic ensembles in a resonant cavity: Quasiperiodic route to chaos and chaotic synchronization, Ann. Phys. 417, 168106 (2020).
  * Gegg and Richter [2016] M. Gegg and M. Richter, Efficient and exact numerical approach for many multi-level systems in open system cqed, New J. Phys. 18, 043037 (2016).
  * Daley [2014] A. J. Daley, Quantum trajectories and open many-body quantum systems, Adv. Phys. 63, 77–149 (2014).
  * Morrison and Parkins [2008] S. Morrison and A. S. Parkins, Dissipation-driven quantum phase transitions in collective spin systems, J. Phys. B 41, 195502 (2008).
  * Lee  _et al._ [2014] T. E. Lee, C.-K. Chan, and S. F. Yelin, Dissipative phase transitions: Independent versus collective decay and spin squeezing, Phys. Rev. A 90 (2014).
  * Russomanno  _et al._ [2017b] A. Russomanno, F. Iemini, M. Dalmonte, and R. Fazio, Floquet time crystal in the Lipkin-Meshkov-Glick model, Phys. Rev. B 95, 214307 (2017b).
  * Lorenz [1963] E. N. Lorenz, Deterministic nonperiodic flow, J. Atmos. Sci. 20, 130 (01 Mar. 1963).
  * Dimer  _et al._ [2007] F. Dimer, B. Estienne, A. S. Parkins, and H. J. Carmichael, Proposed realization of the Dicke-model quantum phase transition in an optical cavity QED system, Phys. Rev. A 75, 013804 (2007).
  * Johansson  _et al._ [2012] J. Johansson, P. Nation, and F. Nori, QuTiP: An open-source python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 183, 1760 (2012).
  * Johansson  _et al._ [2013] J. Johansson, P. Nation, and F. Nori, QuTiP 2: A python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 184, 1234 (2013).
  * Zohar  _et al._ [2013] E. Zohar, J. I. Cirac, and B. Reznik, Quantum simulations of gauge theories with ultracold atoms: Local gauge invariance from angular-momentum conservation, Phys. Rev. A 88, 023617 (2013).
  * Ma  _et al._ [2011] J. Ma, X. Wang, C. Sun, and F. Nori, Quantum spin squeezing, Phys. Rep. 509, 89–165 (2011).
  * Dalla Torre  _et al._ [2016] E. G. Dalla Torre, Y. Shchadilova, E. Y. Wilner, M. D. Lukin, and E. Demler, Dicke phase transition without total spin conservation, Phys. Rev. A 94, 061802 (2016).
  * Riera-Campeny  _et al._ [2020b] A. Riera-Campeny, M. Moreno-Cardoner, and A. Sanpera, Time crystallinity in open quantum systems, Quantum 4, 270 (2020b).
  * Harbola  _et al._ [2006] U. Harbola, M. Esposito, and S. Mukamel, Quantum master equation for electron transport through quantum dots and single molecules, Phys. Rev. B 74, 235309 (2006).
  * Cattaneo  _et al._ [2020] M. Cattaneo, G. L. Giorgi, S. Maniscalco, and R. Zambrini, Symmetry and block structure of the liouvillian superoperator in partial secular approximation, Phys. Rev. A 101, 042108 (2020).
  * Dorn  _et al._ [2021] G. Dorn, E. Arrigoni, and W. von der Linden, Efficient energy resolved quantum master equation for transport calculations in large strongly correlated systems, Journal of Physics A: Mathematical and Theoretical 54, 075301 (2021).
  * Henriet  _et al._ [2020] L. Henriet, L. Beguin, A. Signoles, T. Lahaye, A. Browaeys, G.-O. Reymond, and C. Jurczak, Quantum computing with neutral atoms, Quantum 4, 327 (2020).
  * Zhang  _et al._ [2017] J. Zhang, G. Pagano, P. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. Gorshkov, Z. Gong, and C. Monroe, Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator., Nature 551, 601 (2017).
  * Puri  _et al._ [2017] S. Puri, C. K. Andersen, A. L. Grimsmo, and A. Blais, Quantum annealing with all-to-all connected nonlinear oscillators, Nat. commun. 8, 1 (2017).
  * Angerer  _et al._ [2018] A. Angerer, K. Streltsov, T. Astner, S. Putz, H. Sumiya, S. Onoda, J. Isoya, W. J. Munro, K. Nemoto, J. Schmiedmayer, _et al._ , Superradiant emission from colour centres in diamond, Nat. Phys. 14, 1168 (2018).
  * Clerk  _et al._ [2020] A. Clerk, K. Lehnert, P. Bertet, J. Petta, and Y. Nakamura, Hybrid quantum systems with circuit quantum electrodynamics, Nat. Phys. , 1 (2020).
  * Georgescu  _et al._ [2014] I. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. Mod. Phys. 86, 153–185 (2014).
  * Keesling  _et al._ [2019] A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, _et al._ , Quantum Kibble–Zurek mechanism and critical dynamics on a programmable Rydberg simulator, Nature 568, 207 (2019).
  * Scholl  _et al._ [2020] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Programmable quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms (2020), arXiv:2012.12268 [quant-ph] .
  * Ebadi  _et al._ [2020] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletic, and M. D. Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator (2020), arXiv:2012.12281 [quant-ph] .
  * Tsomokos  _et al._ [2008] D. I. Tsomokos, S. Ashhab, and F. Nori, Fully connected network of superconducting qubits in a cavity, New J. Phys. 10, 113020 (2008).
  * Xu  _et al._ [2020] K. Xu, Z.-H. Sun, W. Liu, Y.-R. Zhang, H. Li, H. Dong, W. Ren, P. Zhang, F. Nori, D. Zheng, and et al., Probing dynamical phase transitions with a superconducting quantum simulator, Science Advances 6, eaba4935 (2020).
  * Kirkpatrick and Thirumalai [1987] T. R. Kirkpatrick and D. Thirumalai, p-spin-interaction spin-glass models: Connections with the structural glass problem, Phys. Rev. B 36, 5388 (1987).
  * Roy  _et al._ [2019] N. Roy, A. Sharma, and R. Mukherjee, Quantum simulation of long-range $xy$ quantum spin glass with strong area-law violation using trapped ions, Phys. Rev. A 99, 052342 (2019).
  * Harris  _et al._ [2018] R. Harris, Y. Sato, A. Berkley, M. Reis, F. Altomare, M. Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, _et al._ , Phase transitions in a programmable quantum spin glass simulator, Science 361, 162 (2018).
  * Rotondo  _et al._ [2015] P. Rotondo, M. Cosentino Lagomarsino, and G. Viola, Dicke simulators with emergent collective quantum computational abilities, Physical Review Letters 114, 10.1103/physrevlett.114.143601 (2015).
  * Frisk Kockum  _et al._ [2019] A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Ultrastrong coupling between light and matter, Nat. Rev. Phys. 1, 19–40 (2019).
  * Pachos and Plenio [2004] J. K. Pachos and M. B. Plenio, Three-spin interactions in optical lattices and criticality in cluster hamiltonians, Phys. Rev. Lett. 93, 056402 (2004).
  * Büchler  _et al._ [2007] H. P. Büchler, A. Micheli, and P. Zoller, Three-body interactions with cold polar molecules, Nature Physics 3, 726 (2007).
  * Capogrosso-Sansone  _et al._ [2009] B. Capogrosso-Sansone, S. Wessel, H. P. Büchler, P. Zoller, and G. Pupillo, Phase diagram of one-dimensional hard-core bosons with three-body interactions, Phys. Rev. B 79, 020503 (2009).
  * Bonnes  _et al._ [2010] L. Bonnes, H. Büchler, and S. Wessel, Polar molecules with three-body interactions on the honeycomb lattice, New Journal of Physics 12, 053027 (2010).
  * Gambetta  _et al._ [2020a] F. M. Gambetta, W. Li, F. Schmidt-Kaler, and I. Lesanovsky, Engineering nonbinary Rydberg interactions via phonons in an optical lattice, Phys. Rev. Lett. 124, 043402 (2020a).
  * Gambetta  _et al._ [2020b] F. M. Gambetta, C. Zhang, M. Hennrich, I. Lesanovsky, and W. Li, Long-range multibody interactions and three-body antiblockade in a trapped Rydberg ion chain, Phys. Rev. Lett. 125, 133602 (2020b).
  * Kockum and Nori [2019] A. F. Kockum and F. Nori, Quantum bits with josephson junctions, in  _Fundamentals and Frontiers of the Josephson Effect_ (Springer, 2019) pp. 703–741.
  * Touzard  _et al._ [2018] S. Touzard, A. Grimm, Z. Leghtas, S. O. Mundhada, P. Reinhold, C. Axline, M. Reagor, K. Chou, J. Blumoff, K. M. Sliwa, S. Shankar, L. Frunzio, R. J. Schoelkopf, M. Mirrahimi, and M. H. Devoret, Coherent oscillations inside a quantum manifold stabilized by dissipation, Phys. Rev. X 8, 021005 (2018).
  * Lescanne  _et al._ [2020] R. Lescanne, M. Villiers, T. Peronnin, A. Sarlette, M. Delbecq, B. Huard, T. Kontos, M. Mirrahimi, and Z. Leghtas, Exponential suppression of bit-flips in a qubit encoded in an oscillator, Nat. Phys. 16, 509–513 (2020).
  * Vrajitoarea  _et al._ [2020] A. Vrajitoarea, Z. Huang, P. Groszkowski, J. Koch, and A. A. Houck, Quantum control of an oscillator using a stimulated josephson nonlinearity, Nat. Phys. 16, 211 (2020).
  * Neeley  _et al._ [2009] M. Neeley, M. Ansmann, R. C. Bialczak, M. Hofheinz, E. Lucero, A. D. O’Connell, D. Sank, H. Wang, J. Wenner, A. N. Cleland, _et al._ , Emulation of a quantum spin with a superconducting phase qudit, Science 325, 722 (2009).
  * Nori [2009] F. Nori, Quantum football, Science 325, 689 (2009).
  * Dutta and Bhattacharjee [2001] A. Dutta and J. K. Bhattacharjee, Phase transitions in the quantum Ising and rotor models with a long-range interaction, Phys. Rev. B 64, 184106 (2001).
  * Manzano and Hurtado [2014] D. Manzano and P. I. Hurtado, Symmetry and the thermodynamics of currents in open quantum systems, Phys. Rev. B 90, 125138 (2014).