# Spatiotemporally ordered patterns in a chain of coupled dissipative kicked
rotors

Angelo Russomanno  Scuola Superiore Meridionale, Università di Napoli Federico
II Largo San Marcellino 10, I-80138 Napoli, Italy and Dipartimento di Fisica,
Università di Napoli Federico II, Monte Sant’Angelo, I-80126 Napoli, Italy

###### Abstract

In this work we consider the dynamics of a chain of many coupled kicked rotors
with dissipation. We map a rich phase diagram with many dynamical regimes. We
focus mainly on a regime where the system shows period doubling, and forms
patterns that are persistent and depend on the stroboscopic time with period
double than that of the driving: The system shows a form of spatiotemporal
ordering analogous to quantum Floquet time crystals. Spatiotemporally ordered
patterns can be understood by means of a linear-stability analysis that
predicts an instability region that contains the spatiotemporally ordered
regime. The boundary of the instability region coincides with the lower
boundary of the spatiotemporally ordered regime, and the most unstable mode
has length scale double than the lattice spacing, a feature that we observe in
the spatiotemporally ordered patterns: Period doubling occurs both in time and
space. We propose an implementation of this model in an array of SQUID
Josephson junctions with a pulsed time-periodic flux.

##  I Introduction

In many-body dynamical systems out of equilibrium, ordered coherent patterns
in space and time naturally appear, from the convective Rayleigh-Bénard cells,
to the heart beats, to the Belusov-Zabotinskii reaction, to synchronization
[1, 2, 3]. In this context there are some universal properties. One of these
is the period-doubling bifurcation cascade [4, 5]. In simple systems
undergoing a periodic driving, like one-dimensional maps, one can see that
increasing a parameter, the system undergoes a sequence of period doublings.
At each of these transitions, a response of the system appears with period
double than the regime before. As a result, one has a response with period
$2^{q}$ times the driving, and for the parameter tending towards a finite
value, $q\to\infty$ and the system becomes chaotic. This phenomenon occurs
with universal scaling properties near the transition to chaos, independently
of the precise choice of the map [6, 7].

The period-doubling cascade has been observed in Nature in many contexts, from
the convection rolls in water and mercury [8, 9], to nonlinear electronic
circuits [10, 11], neurons [12], and infinite-range dissipative quantum spin
systems described by a mean-field theory in the thermodynamic limit [13]. All
these systems have in common the fact that they are nonlinear, are described
by few effective variables, and undergo a periodic driving.

When the periodic driving is applied to many-body systems, the situation
changes: As argued in [14], the generic response is the period doubling.
Responses at larger multiples of the period are possible, but only in case of
special symmetries. Any spatially ordered pattern with a period $m$ times the
driving with $m>2$ is doomed to be destroyed by the growth of bubbles. The
only possible stable patterns are the ones with period doubling, where at each
period the driving exchanges the inner and the outer of the bubble, that
therefore alternatively grows and shrinks.

This result is of great importance for the recent researches on discrete time
crystals in classical noisy periodically-driven systems [15, 16, 17], where a
persistent response at a multiple of the driving period appears only in the
thermodynamic limit [18], and in all the examples known at this time this
response occurs as a period doubling, in agreement with the results of [14]
described above. (See also Refs. [19, 20, 21] for general reviews on time
crystals.)

Motivated by this framework, we aim to understand how the period-doubling
cascade changes when a many-body context is considered. In order to do that,
we go beyond the more usual framework of coupled maps [22, 23, 24] and
consider an array of coupled kicked rotors with dissipation, a model easy to
numerically simulate, also at large sizes. Without dissipation, this model
reduces to a slight generalization [25] of the coupled kicked rotors, showing
Hamiltonian chaos [26, 27, 28, 29, 30, 31]. Without the coupling this model
reduces to the single dissipative kicked rotor, known also as Zaslavsky map
[32, 33]. This model shows a peculiar strange attractor and a very interesting
dynamics [32]. At the onset of chaos the single-rotor model shows a behavior
very similar to two parallel period doubling cascades, and we can study its
fate in the case of many coupled rotors.

We probe the system stroboscopically, that’s to say at discrete times, integer
multiple of the driving period, and we consider appropriate averages over
random initial conditions. Using some “period-doubling order parameters”
inspired by the literature on time crystals [34, 35] we see that the period-
doubling cascade is washed away and the model can show essentially only period
doubling in the regime of regular dynamics. In small parameter ranges there is
a response at a period 4 times the driving, but it is many order of magnitude
smaller than the one at period 2. So, we find that the findings of [14] are
essentially confirmed, with the small period 4 response due probably to the
fact that this model has continuous onsite variables.

Beyond the period $m$-tupling behavior, this model shows a very rich dynamical
behavior, and we summarize it in the phase diagram shown in Fig. 1.
Anticipating a little bit, on the axes there are two parameters, $J$ and $K$,
of the model (it is discussed in detail in Sec. II), while another parameter,
$\gamma$, is kept fixed at $\gamma=0.8$.

Let us first of all focus on the regime where period doubling occurs, that in
the phase diagram we term “Spatiotemporal ordering”. This name is due to the
fact that, whenever period doubling occurs, the system spontaneously organizes
in space, breaking the translation symmetry and giving rise to patterns. These
patterns are stable and persistent in time, change with a period twice the one
of the driving, and so give rise to the period doubling. They are an example
of an effect of nonlinear dynamics in spatially extended systems very common
in Nature [3], and the precise form of the patterns depend on the initial
state chosen.

It is important to remark that the spatiotemporally ordered regime can be
analytically understood by means of a linear stability analysis done according
to the method explained in [3]. We find an instability region that contains
the spatiotemporally ordered regime and shares with it a boundary. Moreover,
the eigenvalues in the unstable region are such that the unstable mode
increases changing sign at every period, consistent with a period-doubling
regime. (That sets on when nonlinear effects become strong enough to contrast
this increase.) Finally, the linear-stability analysis predicts that the most
unstable mode appears at momentum $k=\pi$, that’s to say at a length scale
double than the lattice spacing. Using the Fourier transform, we find that
also the fully formed patterns have the same typical length scale, and so
doubling of the time periodicity comes together with doubling of the space
periodicity. This analysis is very important: From one side provides
analytical predictions (of the lower boundary of the spatiotemporally ordered
regime, and of the typical length scale of the patterns); From the other shows
that spatiotemporally ordered patterns are true patterns in the sense of [3],
that’s to say instabilities with spatial organization around a uniform state,
that increase and are eventually stabilized by nonlinear effects.

Our spatiotemporal ordering has some analogy with time crystals, where time-
translation symmetry breaking comes together with the breaking of an internal
symmetry [34, 36, 37], but is a different phenomenon. Indeed, we see a
persistent period doubling response already at finite sizes, while a true time
crystal should break the space and time translation symmetry only in the
thermodynamic limit. We see here an effect of nonlinear dynamics, a mechanism
physically different to the quantum phase transition-like behavior involved in
quantum Floquet time crystals (see for instance the discussion in [35]).

There are other regimes in the phase diagram of Fig. 1. Let us start with the
“Pattern” one, where the system shows persistent patterns independent of the
stroboscopic time. In this regime the system breaks the space translation
symmetry, but not the discrete time translation symmetry. Unfortunately we are
not able to interpret these time-independent patterns in terms of a linear-
stability analysis. We see that the time-independent patterning exists only
for $K$ smaller than a threshold ($K<0.45$)and find a small region (surrounded
by yellow lines in Fig. 1), where “weak patterning” occurs: In a very jagged
and seemingly fractal way, there are points without patterning and points with
patterns with very small amplitude. The transition from time-independent
patterns to spatiotemporal ordering can be seen in the typical length scale of
the patterns, that suddenly drops at the boundary. This length scale in the
spatiotemporally ordered case approximately coincides with the one
analytically predicted using the linear-stability analysis.

In the region labeled as “Trivial”, the system relaxes to a uniform and time-
independent condition, where all the momenta are vanishing. In the chaotic
regime, in contrast, nearby trajectories in the phase space diverge
exponentially from each other, and the dynamics is given by aperiodic
oscillations in space and time. Here the largest Lyapunov exponent (LLE – the
measure of the rate of exponential divergence) is positive. The transition
from negative (regular dynamics) to positive (chaotic dynamics) is always
sharp, but along the segment marked in red, where the LLE is near to 0 and
intermittently becomes slightly positive.

Many of the transitions between the different regimes described above can be
seen in the behavior of the kinetic energy per site, a quantity often
considered in studies about kicked rotors [26, 27, 38, 25]. In contrast with
cases with Hamiltonian chaos, where this quantity increases steadily and
without a bound [26, 27, 38, 31], here the kinetic energy per site reaches an
asymptotic value, that does not scale with the system size.

This is an important information for experimental realizations. Indeed, we
propose an experimental realization of this model with an array of SQUID
Josephson junctions with a time-periodic pulsed magnetic flux. The kinetic
energy per site translates in the charging energy per site of the
superconducting system, and the fact that it is bounded provides the
possibility to tune the parameters so that this energy per site stays below
the superconducting gap. In this way the array of Josephson junctions can keep
superconductivity for long times, and be correctly described by our model.

The paper is organized as follows. In Sec. II we introduce the model we study.
In Sec. III we perform the linear stability analysis around the trivial
regime, and show that there is an instability appearing exactly where the
spatiotemporally ordered regime sets in. In Sec. IV we discuss the patterns
and their properties – amplitude and typical length scale – and relate them to
the predictions of linear-stability analysis. In Sec. V we study chaos by
means of the largest Lyapunov exponent, and map the boundary line of the
“Chaos” region in Fig. 1. In Sec. VI we study the period $m$-tupling in the
single- and many- rotor cases, and show that in the many-rotor model only
period doubling survives. In Sec. VII we discuss the behavior of the kinetic
energy per site. In Sec. VIII we discuss how to realize our model in an array
of SQUID Josephson junctions with pulsed time-periodic magnetic flux. In Sec.
IX we draw our conclusions.

---
Figure 1: Sketch of the dynamical phase diagram for $\gamma=0.8$. (For other
values of $\gamma$ the situation is qualitatively similar.) We recognize the
trivial regime (where relaxation to a uniform asymptotic condition occurs),
the “Pattern” regime (where the system breaks space translation symmetry by
generating persistent time-independent patterns), the “Spatiotemporal ordered”
regime (where the persistent patterns depend on the stroboscopic time with a
period double the one of the driving, breaking thereby space and time-
translation symmetry), and the chaotic regime, where the dynamics is aperiodic
in space and time. Inside the yellow region there is the weak patterning
regime (see text for a description).

##  II Model

We add dissipation to the many-body generalization of the kicked rotor
considered in [25], a slight generalization to the paradigmatic model of many-
body Hamiltonian chaos theory studied in [26, 27, 28, 29, 30, 31], that can be
even realized experimentally with an array of bosonic Josephson junctions [39,
31]. The purely Hamiltonian model is given by the kicked Hamiltonian

$H(t)=\sum_{j=1}^{L}\left[\frac{p_{j}^{2}}{2}-\delta_{1}(t)\left(J\cos(\theta_{j}-\theta_{j+1})+K\cos(\theta_{j})\right)\right]\,,$
(1)

where $\theta_{j}$, $p_{j}$ are pairs of canonically conjugated variables. We
assume periodic boundary conditions ($\theta_{L+1}=\theta_{1}$), and we have
defined the periodic delta function
$\delta_{\tau}(t)\equiv\sum_{n\in\mathbb{N}}\delta(t-n\tau)$ as in [29]. The
fact that $\tau=1$ is not a limitation, because $\tau$ can always be rescaled
away in the definition of the couplings [38]) We see that the system is a
1-dimensional chain of rotors with kicked short-range interactions and we
focus on the stroboscopic dynamics, probing the system at discrete times
$t_{n}=n^{-}$, where the superscript “-” means that we look at the system just
before the $n$-the kick has been applied.

Writing the canonical equations of the dynamics

$\dot{\theta}_{j}=\partial_{p_{j}}H,\quad\dot{p}_{j}=-\partial_{\theta_{j}}H,$

and using a standard analysis (essentially the integration of the second
equation across the $\delta$ function – see for instance [38, 26, 27, 40]) we
see that the dynamics from $t_{n}$ to $t_{n+1}$ is described by a discrete map

$\displaystyle p_{j}^{(n+1)}$
$\displaystyle=p_{j}^{(n)}-J\left[\sin(\theta_{j}^{(n)}-\theta_{j+1}^{(n)})+\sin(\theta_{j}^{(n)}-\theta_{j-1}^{(n)})\right]$
$\displaystyle-K\sin(\theta_{j}^{(n)})$ $\displaystyle\theta_{j}^{(n+1)}$
$\displaystyle=\theta_{j}^{(n)}+p_{j}^{(n+1)}\,,$ (2)

for $j=1,\,\ldots,\,L$, where we have written for simplicity
$p_{j}(t_{n})=p_{j}(n)$, $\theta_{j}(n)=\theta_{j}(t_{n})$.

We add to this model a dissipation, so that between one kick and the next the
momenta are damped by a factor $0<\gamma<1$. The resulting map is

$\displaystyle p_{j}^{(n+1)}$ $\displaystyle=\gamma
p_{j}^{(n)}-J\left[\sin(\theta_{j}^{(n)}-\theta_{j+1}^{(n)})+\sin(\theta_{j}^{(n)}-\theta_{j-1}^{(n)})\right]$
$\displaystyle-K\sin(\theta_{j}^{(n)})$ $\displaystyle\theta_{j}^{(n+1)}$
$\displaystyle=\theta_{j}^{(n)}+p_{j}^{(n+1)}\,,$ (3)

for $j=1,\,\ldots,\,L$. This is the many-body generalization of the well-known
single kicked rotor with dissipation [33] – known also as Zaslavsky map [32]
–, to which our model reduces for $L=1$

$\displaystyle p^{(n+1)}$ $\displaystyle=\gamma p^{(n)}-K\sin(\theta^{(n)})$
$\displaystyle\theta^{(n+1)}$ $\displaystyle=\theta^{(n)}+p^{(n+1)}\,.$ (4)

Here we keep fixed $\gamma=0.8$ (other choices of $\gamma$ would give a
qualitatively similar phase diagram). Where not otherwise specified, we will
always consider appropriate averages over random initial conditions, taking in
each of them all the $p_{j}^{(0)}$ and the $\theta_{j}^{(0)}$ from a random
distribution uniform in the interval $[-1,1]$. Of this model we study chaotic
properties, period $m$-tupling properties, and patterning. For the first one
we use the largest Lyapunov exponent that we consider in Sec. V, for the
second we define an appropriate set of order parameters in Sec. VI, and
discuss patterning in Sec. IV. We can gain analytical insight in patterning
using a linear-stability analysis, as we are going to show in the next
section.

##  III Linear stability analysis

In order to understand when the patterning appears, we use the method
described in [3], and linearize the dynamical equations Eq. (II) around the
uniformly vanishing solution of Eq. (II), the one $p_{j}^{(n)}\equiv 0$,
$\theta_{j}^{(n)}\equiv 0$. This is the stable solution for $K$ and $J$ small
enough (“Trivial” regime), and corresponds to the trivial phase in Fig. 1. The
linearized equations are

$\displaystyle p_{j}^{(n+1)}$ $\displaystyle=\gamma
p_{j}^{(n)}-J\left[2\theta_{j}^{(n)}-\theta_{j-1}^{(n)}-\theta_{j+1}^{(n)}\right]-K\theta_{j}^{(n)}$
$\displaystyle\theta_{j}^{(n+1)}$
$\displaystyle=\theta_{j}^{(n)}+p_{j}^{(n+1)}\,.$ (5)

At this point let us apply the Fourier transform
$p_{j}^{(n)}=\sum_{k}p_{k}^{(n)}\operatorname{e}^{ikj}$,
$\theta_{j}^{(n)}=\sum_{k}\theta_{k}^{(n)}\operatorname{e}^{ikj}$, where $k$
can take the values consistent with periodic boundary conditions, that’s to
say $k=\frac{2\pi n}{L}$, with $n=0,\,1,\,\ldots,\,L-1$ integer. After the
Fourier transform we get

$\left(\begin{array}[]{c}p_{k}^{(n+1)}\\\
\theta_{k}^{(n+1)}\end{array}\right)=\left(\begin{array}[]{cc}\gamma&-\left[2J\left(1-\cos(k)\right)+K\right]\\\
\gamma&1-\left[2J\left(1-\cos(k)\right)+K\right]\end{array}\right)\left(\begin{array}[]{c}p_{k}^{(n)}\\\
\theta_{k}^{(n)}\end{array}\right)\,.$ (6)

To study stability, we need to evaluate the eigenvalues of this $2\times 2$
matrix. There is instability (and pattern generation) whenever the absolute
value of at least one of the eigenvalues is larger than 1. The eigenvalues
$\mu$ are the solutions of the equation

$\mu^{2}-\left\\{\gamma+1-2J\left(1-\cos(k)\right)-K\right\\}\mu+\gamma=0\,.$

These solutions are explicitly

$\displaystyle\mu_{\pm}^{(k)}=\frac{1}{2}\Big{[}\gamma+1-2J\left(1-\cos(k)\right)-K$
$\displaystyle\pm\sqrt{\left[\gamma+1-2J\left(1-\cos(k)\right)-K\right]^{2}-4\gamma}\Big{]}\,.$
(7)

In order to see when the uniformly vanishing solution becomes unstable, we
have to look at the maximum over $k$ of the absolute values of
$\mu_{\sigma}^{(k)}$. Let us call them
$\Lambda_{\sigma}=\max_{k}|\mu_{\sigma}^{(k)}|$, with $\sigma\in\\{-,\,+\\}$:
When at least one of these maxima becomes larger than 1 the instability sets
on and the system can generate patterns. We plot some example of
$\Lambda_{\pm}$ versus $K$ in Fig. 2(a).

We see that $\Lambda_{-}$ increases with $K$, and there is a critical value
$K^{*}$ beyond which $\Lambda_{-}$ becomes definitively larger than 1. $K^{*}$
therefore sets the onset of the instability. We plot $K^{*}$ versus $J$ in
Fig. 2(b), and see that perfectly coincides with the lower boundary line of
the period-doubling spatiotemporally ordered regime (see the phase diagram in
Fig. 1 and the analysis in Sec. VI.3, to see how this boundary line is
numerically obtained). So, the linear-stability analysis is able to
quantitatively reproduce the setting on of the spatiotemporal ordering, but
does not describe the regime of time-independent patterning.

(a)
---
(b)
Figure 2: (a) Maximum absolute value of the eigenvalues $\Lambda_{\sigma}$
versus $K$, for $\sigma\in\\{-,\,+\\}$ and two values of $J$. There is a
$K^{*}$ such that, for $K>K^{*}$, $\Lambda_{-}$ becomes definitively larger
than 1 (red horizontal line), marking the setting on of the instability.
$\Lambda_{+}$ versus $K$ is independent of $J$ and equals 1 when $K=0$. (b)
Comparison between the instability boundary line provided by the linear
stability analysis and the lower boundary line of the spatiotemporally ordered
regime: They coincide perfectly and are very good approximated by the
analytical formula in Eq. (8).

When there is instability, the maximum of $|\mu_{-}^{(k)}|$ occurs always for
$k_{\rm max}^{+}=\pi$, as we have numerically checked. So, $k=\pi$ is the
momentum of the most unstable mode. This implies that the instability provides
a pattern with a typical length scale double than the lattice step. So, we see
that doubling of the time periodicity comes together with doubling of the
space periodicity of the model, a fact that is confirmed by the numerical
analysis of the fully nonlinear model (see Sec. IV.3).

Looking at $\mu_{\pm}^{(k_{\rm max}^{\pm})}$ we find that when the instability
sets on, the eigenvalues are both fully real and negative. This implies that,
in the long-time limit, the most unstable mode increases changing sign at
every period, consistently with the presence of period doubling. (The period
doubling fully develops after a transient, when nonlinear effects set in, and
balance the unbounded exponential growth).

This information allows us to provide an analytical approximation for the
boundary line. We have the most unstable mode for $k=\pi$ and it appears when
both the eigenvalues become real. So, we can substitute $k=\pi$ inside Eq.
(III) and impose that the argument of the root is vanishing. In this way we
find two solutions, one of them is

$K^{*}=-4J+2\sqrt{\gamma}+\gamma+1\,.$ (8)

We plot it in Fig. 2(b) and we see that it provides a very good approximation
to the instability boundary line.

On the opposite, $\Lambda_{+}$ – the maximum of $|\mu_{+}^{(k)}|$ – is
independent of $J$, as we can see in Fig. 2(a). In the interval where
$\Lambda_{+}$ is constant, the $|\mu_{+}^{(k)}|$ are independent of $k$,
otherwise we have checked that the maximum occurs for $k_{\rm max}^{-}=0$.
This happens for instance near $K=0$. In this region, substituting $k=0$ in
Eq. (III), we can analytically check that $\Lambda_{+}=1$ for $K=0$, whenever
$0<\gamma<1$ (see also Fig. 2). For $K=0$, therefore, there is marginally
stable mode with $|\lambda_{+}^{(0)}|=1$, that is fragile to the tiniest
perturbation (see for instance [41]). It is therefore not surprising that in
the limit $J\gg K$ the system can only show chaos (see the phase diagram in
Fig. 1).

##  IV Pattern formation

###  IV.1 Examples of persistent patterns

When the system is in the regime “Trivial” in Fig. 1, all the values of
$p_{j}$ relax to 0. Otherwise, in the regimes “Pattern” and “Spatiotemporal
order”, after a transient is died away, the system spontaneously forms
patterns of $p_{j}$ that are persistent in time. We find that different
initial conditions give rise to different asymptotic patterns, a phenomenon
common in nonlinear dynamics.

In absence of period doubling (the regime “Pattern”) the patterns are
independent of the stroboscopic time. We show an example thereof in Fig. 3(a).
We initialize with one random initial condition and wait that the initial
transient dies out. We see that the pattern is constant in the stroboscopic
time, so it lies unchanged whatever the value of $n$. Because the pattern
comes back to itself after each cycle, the rotors behave in a perfectly
synchronized way.

In the regime “Spatiotemporal order”, instead, the patterns are associated
with period doubling. In this regime, the persistence in time of the pattern
means that it changes at each stroboscopic time, and comes back to itself
after two cycles (so also here the rotors are perfectly synchronized, with a
period double than the driving). We show an example thereof in Fig. 3(b). The
period doubling appears in the fact that the pattern has a constant form for
$n$ even ($n=3\cdot 10^{4}$, $n=6\cdot 10^{4}$) and a different equally
constant form for $n$ odd ($n=3\cdot 10^{4}+1$, $n=6\cdot 10^{4}+1$).

Notice the staggered spatial ordering, with oscillations at a period double
than the one of the lattice, separated by defects. This suggests that the most
unstable mode with momentum $k=\pi$ provided by linear-stability analysis (see
Sec. III) keeps its spatial periodicity double than the lattice also when it
develops a stable pattern, as we are going to quantitatively confirm in the
subsequent analysis.

---
Figure 3: Examples of asymptotic patterns. (a) Example for the regime
“Pattern” in Fig. 1. Here the pattern is persistent and independent of the
stroboscopic time $n$ (Numerical parameters $J=0.1,\,K=2.71,\,L=512$, one
given random initial condition). (b) Example for the regime “Spatiotemporal
ordering” in Fig. 1. Here the pattern is persistent and has a periodicity
double than the driving, so the patterns for $n$ even coincide with each
other, and the same with the patterns for $n$ odd. Notice the regions with
space periodicity double than the lattice separated by defects. (Numerical
parameters $J=0.1,\,K=4.04,\,L=512$, one given random initial condition. For
clarity we have plotted only part of the pattern.)

###  IV.2 Analysis of the pattern amplitude

In order to study the existence and the properties of the patterns, it is
important to quantify them. For that purpose we introduce two quantities. The
first one is the pattern amplitude. To evaluate it, we fix $n\gg 1$, consider
the variance of $p_{j}^{(n)}$ over space, average it over random initial state
realizations, and evaluate the square root, namely

$\delta
p_{n}=\left[\;\overline{\frac{1}{L}\sum_{j}(p_{j}^{(n)})^{2}-\left(\frac{1}{L}\sum_{j}p_{j}^{(n)}\right)^{2}}\;\right]^{1/2}\,.$
(9)

This quantity is very important, because marks the existence of the patterns
(it vanishes in the trivial state). We show some examples of $\delta p$ versus
$K$ in Fig. 4. In all these figures we consider two values of $n$ ($n=3\cdot
10^{4}$ and $n=6\cdot 10^{4}$) to show that $\delta p_{n}$ has converged in
time, and consider $L=512$, large enough that all finite-size effects have
disappeared. Let us first focus our attention on the case $J=0.1$ [Fig. 4(a)].
We see first of all that $\delta p_{n}$ is independent of $n$ for $n\geq
3\cdot 10^{4}$. We see many features, let us discuss them moving from the
right to the left.

At the onset of chaos (red vertical line on the extreme right) we notice that
$\delta p_{n}$ starts abruptly to increase with a discontinuous derivative. In
the chaotic regime, patterns depend on time in an aperiodic fashion, but the
average over initial-state realizations provides a $\delta p_{n}$ that does
not depend on $n$. At the onset of the period doubling (green vertical line,
second from the right) we see that $\delta p_{n}$ shows no special feature.

For small $J$ the period-doubling regime is fully contained in a region where
there is patterning ($\delta p_{n}>0$) and the disappearance of period
doubling gives no discontinuity, neither on $\delta p_{n}$ nor on its
derivative. The situation is so for $J<0.45$; in contrast for $J\geq 0.45$ the
threshold for patterning coincides with the one for period doubling. [We can
see an example of that also in the plot in Fig. 4(c)]. It is important to
emphasize that period doubling appears always in association with the
appearance of patterns, that’s why we define “spatiotemporally ordered” the
range of parameters where period doubling appears.

Between the two yellow vertical lines in Fig. 4(a) there is the weak
patterning regime. It is characterized by a jagged profile where very small
values of $\delta p_{n}$ (order $10^{-3}$) alternate with vanishing values
(and then no pattern at all). We show some magnification of this jagged
profile in Fig. 5. This regime disappears for $J\simeq 0.14$, but also for
larger $J$ we can see a marked local minimum of $\delta p_{n}$ for $J$ just
below the onset of period doubling [see Fig. 4(b)].

For $J=0.5$ we are well inside the regime where the onset of patterning and of
period doubling coincide, and patterns in the regular regime show always a
time dependence of period 2 (spatiotemporal order) [Fig. 4(c)]. At the
transition to chaos one can see the same features as in the other two cases.

---
Figure 4: $\delta p_{n}$ versus $K$ for different values of $J$ and of $n$
(chosen such that convergence has been attained). The vertical lines mark the
boundaries of the different regimes listed in Fig. 1: the red line bounds from
below the chaotic regime, the green one the period-doubling (spatiotemporal
ordering) regime, the purple one the patterning regime and between the yellow
lines lies the weak patterning regime. The weak patterning regime exists only
in panel (a), and in panel (c) the onset of patterning coincides with the
onset of period doubling. Numerical parameters: $N_{\rm r}=10^{3},\,L=512$;
(a) $J=0.1$, (b) $J=0.2$; (c) $J=0.5$.
---
Figure 5: Examples of $\delta p_{n}$ versus $K$ in the weak patterning regime.
Numerical parameters: $L=160$, $N_{\rm r}=10^{3}$.

###  IV.3 Analysis of the typical length scale of the pattern

Another important property is the shape of the pattern. We can see in Fig. 3
that the patterns oscillate in space, and to these oscillations we can
associate a wavelength, providing us the typical length scale of the pattern.
We can estimate this length scale in the following way. Exploiting that the
number of space oscillations is independent of $n$ (we can see some examples
thereof in Fig. 3), we choose $n\gg 1$ and numerically perform the space
Fourier transform of $p_{j}^{(n)}$ as

$f_{n}(k)=\frac{1}{L}\sum_{j=1}^{L}\operatorname{e}^{ikj}p_{j}^{(n)}\,,$ (10)

where $k=2\pi\ell/L$, with $\ell\in\\{1,\,\ldots,\,L\\}$ integer. Then we
evaluate the power spectrum $|f_{n}(k)|^{2}$ and choose the value $k_{\rm
max}^{(n)}$ where this power spectrum shows a maximum. The corresponding
wavelength is $\lambda_{\rm max}^{(n)}=2\pi/k_{\rm max}^{(n)}$. We perform the
logarithmic average of this quantity over random realizations of the initial
state and we get an estimate of the typical length scale of the patterns for a
given set of parameters

$\lambda_{n}\equiv\exp\left(\overline{\log\lambda_{\rm max}^{(n)}}\right)\,.$
(11)

(We choose the logarithmic average because the distributions of the
$\lambda_{\rm max}$ are broad).

---
Figure 6: Logarithm of the typical length scale of the patterns versus
$\gamma$ for different system sizes and values of $J$ at time $n=10^{5}$. The
errorbars are evaluated as $\delta\log_{10}\lambda^{*}=\frac{1}{\sqrt{N_{\rm
r}}}\left[\overline{(\log_{10}\lambda)^{2}}-(\overline{\log_{10}\lambda})^{2}\right]^{1/2}$.
The blue horizontal line marks the value $\lambda_{n}=2$ predicted by the
linear-stability analysis in presence of instability for the most unstable
mode. There is a good agreement with $\lambda_{n}$ in the spatiotemporally
ordered phase, near its lower phase boundary. The parameters for each panel
are written in the corresponding headings. $N_{\rm r}=10^{3}$.

We plot some examples of $\log_{10}\lambda_{n}$ versus $K$, for $L=512$ and
different values of $J$, in Fig. 6. For $J=0.1$ [Fig. 6(a)] we see a sudden
drop inside the weak-patterning regime (in part of which $\lambda_{n}$ it is
not even defined), and at the onset of the spatiotemporal ordered regime (dark
green vertical line). So, although the amplitude of the patterns was not able
to see the transition to the spatiotemporal ordering, it can be clearly seen
in the typical length scale of the pattern.

The blue horizontal lines mark the value $\lambda_{n}=2\pi/\pi=2$ predicted by
the linear-stability analysis for the most unstable mode. We find a good
agreement with $\lambda_{n}$ inside the spatiotemporally ordered regime, when
we consider $K$ near the lower phase boundary [the same occurs also for other
values of $J$ – see Fig. 6(b,c)]. So, the linear-stability analysis predicts
correctly not only the position of the lower phase boundary of the
spatiotemporally ordered regime, but also the typical length scale of the
pattern provided by the instability. At the onset of chaos (red vertical line)
the typical length scale shows ha huge peak, increasing of one order of
magnitude and then suddenly dropping. In some sense at this transition the
range of the correlations of the system increases.

For $J=0.2$ [Fig. 6(b)] the weak patterning regime has disappeared, but
$\log_{10}\lambda_{n}$ shows still a sudden drop at the onset of the
spatiotemporal ordered regime and the peak at the onset of chaos.

For $J=0.5$ [Fig. 6(c)] there are no time-independent patterns and the
threshold for period doubling coincides with that of patterning, and there is
no peak of the typical length scale at the onset of chaos.

##  V Chaos

To quantify chaos, that’s to say exponential increase in time of the distance
of two nearby trajectories, we use the largest Lyapunov exponent (LLE – see
for instance [42, 43]). It is defined in the following way. Considering two
dynamics ${\bf
X}^{(n)}=(p_{1}^{(n)},\,\ldots,\,p_{L}^{(n)},\,\theta_{1}^{(n)},\,\ldots,\,\theta_{L}^{(n)})$,
${{\bf
X}^{(n)}}^{\prime}=({p_{1}^{(n)}}^{\prime},\,\ldots,\,{p_{L}^{(n)}}^{\prime},\,{\theta_{1}^{(n)}}^{\prime},\,\ldots,\,{\theta_{L}^{(n)}}^{\prime})$
with different initial conditions ${\bf X}^{(0)}$, ${{\bf X}^{(0)}}^{\prime}$
such that $||{\bf X}^{(0)}-{{\bf X}^{(0)}}^{\prime}||=\epsilon>0$, the LLE is
defined as

$\text{LLE}=\lim_{\epsilon\to
0}\lim_{n\to\infty}\frac{1}{n}\log\left(\frac{||{\bf X}^{(n)}-{{\bf
X}^{(n)}}^{\prime}||}{\epsilon}\right)\,,$ (12)

where $||\ldots||$ is the Euclidean norm. In a chaotic dynamics, this quantity
evaluates the rate at which nearby trajectories exponentially separate from
each other. So, when the dynamics is chaotic this quantity is positive, in
absence of chaos it is vanishing or negative. To numerically evaluate the LLE
we use the method explained in [44].

The method goes as follows. One considers as above two initial conditions
distant $\epsilon$. After the first cycle one evaluates the distance
$d_{1}=||{\bf X}^{(1)}-{{\bf X}^{(1)}}^{\prime}||$. Then one redefines ${{\bf
X}^{(1)}}^{\prime}$ as

${{\bf X}^{(1)}}^{\prime\prime}={{\bf X}^{(1)}}+\frac{\epsilon}{d_{1}}({{\bf
X}^{(1)}}^{\prime}-{\bf X}^{(1)})\,,$ (13)

so that the distance between ${\bf X}^{(1)}$ and ${{\bf
X}^{(1)}}^{\prime\prime}$ becomes $\epsilon$ again. With these initial
conditions one performs another stroboscopic-evolution step getting some ${\bf
X}^{(2)}$ and ${{\bf X}^{(2)}}^{\prime}$. So one gets another value of the
distance $d_{2}$, and performs a redefinition of ${{\bf X}^{(2)}}^{\prime}$ as
in Eq. (13). This cycle is repeated many times, do that one gets a sequence of
distances $d_{1},\,d_{2},\,d_{3},\,\ldots,\,d_{n}$ and the Lyapunov exponent
is given by

$\text{LLE}=\frac{1}{\mathcal{T}}\sum_{k=1}^{\mathcal{T}}\log\frac{d_{k}}{\epsilon}\,,$
(14)

where $\mathcal{T}$ is large enough and $\epsilon$ small enough so that
convergence has been attained.

We use precisely this formula to get the LLE. We choose ${\bf X}^{(0)}$ taking
all the $p_{j}^{(0)}$ and the $\theta_{j}^{(0)}$ from a random distribution
uniform in the interval $[-1,1]$. We take ${{\bf X}^{(0)}}^{\prime}$ equal to
${\bf X}^{(0)}$ everywhere but on the coordinate ${p_{1}^{(0)}}^{\prime}$ that
we take ${p_{j}^{(0)}}^{\prime}=p_{j}^{(0)}+\epsilon$. To make convergence
faster, we average Eq. (14) over $N_{\rm r}$ realizations of the random ${\bf
X}^{(0)}$.

Fixing $J$, we find a quite sudden transition in $K$ from regular behavior
($\text{LLE}<0$) to chaotic behavior ($\text{LLE}>0$), provided the system
size is large enough. This allows to map the boundary line of the chaotic
region shown in Fig. 1. The only region where the mapping of this line is
problematic is around $J=0.04$. Here, there is not a sharp transition from a
negative LLE to a positive one, but a range (the segment in red in Fig. 7)
where the Lyapunov exponent lies near 0 and often becomes slightly positive
(see inset of Fig. 7). This value of $J$ marks an abrupt change in the
behavior of the boundary line of the chaotic region, that for $J<0.04$ keeps a
value similar to the single-particle case ($J=0$) and for $J>0.04$ starts
going down as a straight line (see Fig. 1).

In the limit of $J\gg K$ there is only chaos, consistently with the numerical
stability analysis described above, that for $K=0$ provides a marginally
stable (and therefore very fragile) mode.

---
Figure 7: LLE versus $K$ for different values of $J$ (In the inset $J=0.04$ is
plotted). For each curve we choose the value of $\mathcal{T}$, $N_{\rm r}$ and
$L$ so that convergence is attained (for $J=0$, $\mathcal{T}=10^{5},\,N_{\rm
r}=100,\,L=1$; for $J=0.04$, $\mathcal{T}=3\cdot 10^{5},\,N_{\rm
r}=300,\,L=160$; for $J=0.1$, $\mathcal{T}=10^{5},\,N_{\rm r}=100,\,L=40$; for
$J=0.5$, $\mathcal{T}=5\cdot 10^{5},\,N_{\rm r}=500,\,L=40$.

##  VI Period $m$-tupling

###  VI.1 Bifurcation diagram of the single rotor

In order to discuss period $m$-tupling, let us start with the single
dissipative rotor model Eq. (II). This model displays a period-doubling
cascade similar to the standard one seen in the logistic map. Before
discussing it more quantitatively, it is useful to show it by means of a
bifurcation diagram. In the plots of Fig. 8 we put $K$ on the horizontal
coordinate, and – for each value of $K$ – we write on the vertical coordinate
the last $10^{3}$ stroboscopic values of $p^{(n)}$, for an evolution lasting
$\mathcal{T}=10^{5}$ periods. If for that $K$ the system relaxes to an
asymptotic stroboscopic value, we see a single value on the vertical
coordinate; if there is a period doubling we will see two values, if there
period $m$-tupling for generic $m$ we will see $m$ values.

---
Figure 8: (Panel a) Bifurcation diagram. (Panel b) A magnification thereof
near the upper branch just below the onset of chaos. Notice the two main
parallel bifurcation cascades (for each value of $K$ there are points on only
one of the two).

What we see is a period-doubling cascade, that’s to say a sequence of points
of $K$ where pitchfork bifurcations occur and the number of points doubles.
(See [4], [2], [5] for more details on pitchfork bifurcations and period-
doubling cascade). So one moves to period doubling, to period 8-tupling, to
period 16-pling (all the powers of 2). The period doubling cascade ends at the
onset of chaos, that here starts at $K\simeq 5.978$ as one can see in the
bifurcation plot as the appearance of a region with randomly scattered points
(one can confirm this chaos threshold with the LLE analysis explained in Sec.
V). This period doubling cascade is similar to the logistic map, but there is
more, because there are two parallel bifurcation cascades, as one can see in
Fig. 8(b), and beyond that also some small ones near $K=5.97$. Quite
remarkably, for each value of $K$, the system chooses just one of the
bifurcation cascades, in an apparently random way.

###  VI.2 Definition of the period $m$-tupling order parameters

In order to better quantify these phenomena in a way that can be generalized
to the many-body case, we take inspiration from the literature on discrete
time crystals [34], and define the following period $m$-tupling onsite order
parameters

$\mathcal{O}_{n}^{m}(j)=\left[{p_{j}^{(n)}}\cos\left(\frac{2\pi
n}{m}\right)\right]\,,$ (15)

then we average them over time and sites, and

$\mathcal{O}^{m}=\frac{1}{n_{0}}\sum_{n=\mathcal{T}-n_{0}}^{\mathcal{T}}\frac{1}{L}\sum_{j=1}^{L}\mathcal{O}_{n}^{m}(j)\,,$
(16)

where we choose $n_{0}$ and $\mathcal{T}$ so that the initial transient is
vanished and the sum over $n$ is converged. Our order parameter is the average
over random initial-state realizations [this average is represented by the
symbol $\overline{(\ldots)}$] of the absolute value of this quantity

$\mathcal{O}^{(m)}=\overline{|\mathcal{O}^{m}|}\,.$ (17)

It is not difficult to convince oneself, that – for $m>1$ – if
$\mathcal{O}_{n}^{m}$ shows a response at frequency $m$, then the average over
an infinite time of $\mathcal{O}_{n}^{m}$ is nonvanishing. Notice that the
system could show the linear superposition of responses with different
periods, so a nonvanishing $\mathcal{O}_{n}^{m}$ is a necessary condition for
a period $m$-tupling, but not sufficient. For instance, there can be also a
response with period $2m$ – thereby a period $2m$-tupling – and we could still
find a nonvanishing $\mathcal{O}_{n}^{m}$, as we are going to see below.

It is important to stress that these order parameters are of some usefulness
only in the regular regime where $\text{LLE}<0$. Where there is chaos
($\text{LLE}>0$) the dynamics of the $p_{j}^{(n)}$ shows aperiodic
oscillations, there is a response at all the frequencies (see for instance
[45]), and so all the $\mathcal{O}^{(m)}$, providing no information on the
existence of a period doubling, but just trivially witnessing the chaos in the
dynamics. So we will focus our period $m$-tupling analysis on the regime of
regular dynamics.

We define also the quadratic average

$\mathcal{O}^{(2,m)}=\overline{\left[\frac{1}{n_{0}}\sum_{n=\mathcal{T}-n_{0}}^{\mathcal{T}}\frac{1}{L}\sum_{j=1}^{L}\mathcal{O}_{n}^{m}(j)\right]^{2}}\,,$
(18)

and obtain the uncertainty on $\mathcal{O}^{(m)}$ as

$\delta\mathcal{O}^{(m)}=\frac{1}{\sqrt{N_{r}}}\sqrt{\mathcal{O}^{(2,m)}-(\mathcal{O}^{(m)})^{2}}\,,$
(19)

where $N_{r}$ is the total number of randomness/noise realizations. In all the
numerical analyses that follow, we choose $\mathcal{T}$ finite large enough so
that the limits in Eqs. (17) and (18) are converged.

###  VI.3 Results of the order-parameter analysis

Let us start with the case of the single rotor. We show $\mathcal{O}^{(m)}$
versus $K$ for $m=2,\,4,\,8$ in Fig. 9. We see that the behavior of the order
parameters closely mirrors the one of the bifurcation diagram in Fig. 8.
$\mathcal{O}^{(2)}$ is nonvanishing whenever there is at least period doubling
[Fig. 8(a)], $\mathcal{O}^{(4)}$ is nonvanishing whenever there is at least
period 4-tupling, and $\mathcal{O}^{(8)}$ is nonvanishing whenever there is at
least period 8-tupling [Fig. 8(b)]. We see also that the response at $m=4$ is
two orders of magnitude smaller than the one at $m=2$, and the one at $m=8$
even smaller. Also this property closely mirrors the fact that at each
bifurcation the outcoming branches are much nearer than the ones at the
bifurcation before (see Fig. 8 and the self-similarity analysis in [6]).

Of course the analysis is meaningful whenever the dynamics is regular. After
the onset of chaos (marked by the vertical line in Fig. 8) all the order
parameters are nonvanishing and of the same order of magnitude, because the
chaotic dynamics has contributions at all the frequencies.

(a)
---
(b)
Figure 9: Single rotor. (Panel a) Period-doubling order parameters versus $K$.
(Panel b) Period $4$-tupling and period $8$-pling order parameter versus $K$.
The vertical line marks the onset of chaos, as found with the LLE analysis.
Numerical parameters $\mathcal{T}=4\cdot 10^{4},\,n_{0}=10^{3},\,N_{\rm
r}=10^{3}$.

If we add the interactions we see that the period-doubling cascade is washed
away, as we see in Fig. 10. Focusing on the regular-dynamics regime, we see
that the period-doubling response is dominant. There are still some small
regions with period 4-tupling, but the response at $m=4$ is three orders of
magnitude smaller that the one at $m=2$, making it quite negligible. Again in
the chaotic regime there is a response at all frequencies due to chaos and we
do not consider it. In the figure we choose a specific value of $J$, but we
have checked that the thing is general.

So, we essentially find confirmation of the results of [14]: whenever the
dynamics is regular, and then makes sense speaking about period $m$-tupling,
the dominant response is period doubling, with some negligible contribution at
larger frequency. The fact that this contribution exists is probably due to
the fact that our system has continuous local variables, in contrast with the
discrete local variables of cellular automata [14], and the known classical
dissipative time crystals [15, 16, 17].

(a)
---
(b)
Figure 10: Many rotors with $J=0.1$. (Panel a) Period-doubling order
parameters versus $K$ for different values of $L$. (Panel b) Period
$4$-tupling order parameter versus $K$ for different values of $L$. The
vertical line marks the onset of chaos, as found with the LLE analysis.
Numerical parameters $\mathcal{T}=4\cdot 10^{4},\,n_{0}=10^{3},N_{\rm
r}=10^{3}$.

We see in Fig. 10(a) that $\mathcal{O}^{(2)}$ is vanishing for $K$ up to a
threshold and then abruptly moves from 0. The threshold where this happens
marks the line bounding from below the “Spatiotemporal ordering” regime in
Fig. 1.

##  VII Behavior of the kinetic energy

Some properties of the patterns can be read in the kinetic energy per site
defined as

$E(n)=\frac{1}{2L}\sum_{j}[p_{j}^{(n)}]^{2}\,.$ (20)

We consider its average over initial-state realizations and time

$E=\frac{1}{n_{0}}\sum_{n=\mathcal{T}-n_{0}}^{\mathcal{T}}\overline{E(n)}\,,$
(21)

where $\mathcal{T}$ and $n_{0}$ are chosen so that convergence is reached. We
see first of all that this quantity is finite, also in the chaotic regime (see
Fig. 11). This is an important difference compared with the Hamiltonian case,
where chaos leads to an unbounded increase in time of energy. Moreover we see
in this quantity the same features that we see in $\delta p$ (compare with
Fig. 4). In particular it vanishes in the trivial regime, and shows a
discontinuity in the derivative and a sudden increase at the onset of chaos.
Therefore, also the kinetic energy can be used to find if there is pattern
formation.

---
Figure 11: Onsite kinetic energy averaged over time and initial-state
realizations [Eq. (21)] versus $K$, for different values of $J$ and $L$.
Numerical parameters: $N_{\rm r}=10^{3},\,n_{0}=10^{3},\,\mathcal{T}=6\cdot
10^{4}$.

In Fig. 11 we see also that the energy per site is bounded and does not scale
with the system size (it is actually size-independent for the values of $L$ we
are considering). This is a key point to make possible an experimental
realization with an array of Josephson junctions, as we are going to show in
the next section.

##  VIII Proposal for experimental realization

The Hamiltonian in Eq. (1) can be realized by means of an array of Josephson
junctions (see for instance [46, 47, 48] for details). We show it in Fig. 12.
With the crosses we mark SQUIDs, that’s to say Josephson junctions whose
energy can be tuned with the application of a magnetic flux. So, we have
$E_{J}=E_{J}(\Phi)=E_{J}^{(0)}\cos(2\pi\Phi/\Phi_{0})$ for the junctions in
the upper row and $E_{K}=E_{K}(\Phi)=E_{K}^{(0)}\cos(2\pi\Phi/\Phi_{0})$ for
the junctions in the lower row, where $\Phi$ is the magnetic flux and
$\Phi_{0}=hc/(2e)$ the Cooper-pair flux quantum. For the junctions in the
upper row we neglect the capacitance, for each of those in the lower row the
capacitance is $C$, so that the corresponding charging energy is
$E_{C}=e^{2}/C$. For each of the junction of the lower row the corresponding
capacity is connected in parallel, and each of these parallel circuits is
connected to the ground by means of a resistance $R$. As in [49, 50] we assume
we can neglect the linear inductances that are always present in the wires.

On each site (the green balls) we have two canonically conjugate dynamical
variables, the gauge invariant superconducting phase $\theta_{j}$ and the
charge $q_{j}$ (expressed in units of the electron charge $e$). If we assume
that the superconducting pieces just above the resistances have all the same
phase $\varphi$, up to a shift of the $\theta_{j}$, we can set the phases to 0
($\varphi=0$). When there is no resistance, the circuit is described by the
Hamiltonian

$H=\sum_{j=1}^{L}\frac{1}{2}E_{C}q_{j}^{2}-E_{J}(\Phi)\cos(\theta_{j}-\theta_{j+1})-E_{K}(\Phi)\cos(\theta_{j})\,.$
(22)

---
Figure 12: Josephson junction array for the experimental realization of the
model. Here we show open boundary conditions. In order to get periodic
boundary conditions, one more junction $E_{J}$ is needed connecting the first
and the last sites.

In Appendix A we show that the dynamics with no resistance is described by
this Hamiltonian. In the figure open boundary conditions are represented, but
we can impose periodic boundary conditions by adding a junction $E_{J}$
connecting the first and the last sites. We apply the following time-periodic
protocol:

  1. 1.

For a time $T_{1}$ the flux $\Phi$ is kept equal to $\Phi=\Phi_{0}/4$ (all the
SQUIDs have vanishing Josephson energy and they are closed);

  2. 2.

For a time $T_{2}$ the flux $\Phi$ is kept equal to $\Phi=0$ (all the SQUIDs
have maximum Josephson energy and they are open);

Physically, this corresponds to a periodic protocol with period $T_{1}+T_{2}$
where a kick lasting $T_{2}$ is applied, because during the interval $T_{1}$
the Josephson energies of the SQUIDs vanish and there is no Josephson energy
term acting.

We consider the stroboscopic evolution, looking at what happens from just
before one kick to just before the next, that’s to say from time
$t_{n}=n(T_{1}+T_{2})+T_{1}$ to $t_{n+1}=(n+1)(T_{1}+T_{2})+T_{1}$. The
dynamics is provided by the canonical equations

$\displaystyle\frac{\hbar}{{2e}}\dot{q}_{j}$
$\displaystyle=-\partial_{\theta_{j}}H\,,$
$\displaystyle\frac{\hbar}{{2e}}\dot{\theta}_{j}$
$\displaystyle=\phantom{-}\partial_{q_{j}}H\,.$ (23)

(The coefficient in front of the derivative is immaterial and can be
eliminated by appropriately rescaling the Hamiltonian or the canonical
coordinates.) We can solve these equations and, if the time of the kick is
much smaller than
$T_{2}\ll\min\left(\frac{\sqrt{E_{J}^{(0)}E_{C}}}{\hbar},\,\frac{\sqrt{E_{J}^{(0)}E_{C}}}{\hbar}\right)$,
we can approximate the evolution from time $t_{n}=n(T_{1}+T_{2})+T_{1}$ to
$t_{n+1}=(n+1)(T_{1}+T_{2})+T_{1}$ as a discrete map

$\displaystyle\frac{\hbar}{2e}q_{j}(t_{n+1})$
$\displaystyle=\frac{\hbar}{2e}q_{j}(t_{n})-T_{2}E_{J}^{(0)}(\sin(\theta_{j}-\theta_{j+1})+\sin(\theta_{j}-\theta_{j-1}))$
$\displaystyle-T_{2}E_{K}^{(0)}\sin(\theta_{j})$
$\displaystyle\frac{\hbar}{2e}\theta_{j}(t_{n+1})$
$\displaystyle=\frac{\hbar}{2e}\theta_{j}(t_{n})+T_{2}E_{C}q_{j}(t_{n+1})\,.$
(24)

Applying the change of variables $\theta_{j}^{(n)}=\theta_{j}(t_{n})$,
$p_{j}^{(n)}=2e\frac{T_{2}E_{C}}{\hbar}q_{j}(t_{n})$, and defining $J\equiv
4e^{2}T_{2}^{2}\frac{E_{C}E_{J}^{(0)}}{\hbar^{2}}$, $K\equiv
4e^{2}T_{2}^{2}\frac{E_{C}E_{K}^{(0)}}{\hbar^{2}}$ we get back the Hamiltonian
mapping Eq. (II).

In order to get the model with dissipation, we consider a nonvanishing
resistance $R$. About the presence of the resistances, we emphasize that they
are quite realistic, because they are naturally present and the difficult
thing is removing them, rather than adding. With the resistance, between one
kick and the next, the junctions are switched off and the charges damp as in
the $RC$ circuit, so that the map Eq. (VIII) becomes

$\displaystyle\frac{\hbar}{2e}q_{j}(t_{n+1})$
$\displaystyle=\operatorname{e}^{-RCT_{1}}\big{[}\hbar
q_{j}(t_{n})-T_{2}E_{J}^{(0)}(\sin(\theta_{j}-\theta_{j+1})$
$\displaystyle+\sin(\theta_{j}-\theta_{j-1}))-T_{2}E_{K}^{(0)}\sin(\theta_{j})\big{]}$
$\displaystyle\frac{\hbar}{2e}\theta_{j}(t_{n+1})$
$\displaystyle=\hbar\theta_{j}(t_{n})+E_{C}\left(\frac{1-\operatorname{e}^{-RCT_{1}}}{RC}\right)q_{j}(t_{n+1})\,,$
(25)

provided that $T_{2}\ll RC$. Applying the change of variables
$\theta_{j}^{(n)}=\theta_{j}(t_{n})$,
$p_{j}^{(n)}=2e\frac{E_{C}}{\hbar}\left(\frac{1-\operatorname{e}^{-RCT_{1}}}{RC}\right)q_{j}(t_{n})$,
and defining

$\displaystyle\gamma$ $\displaystyle\equiv\operatorname{e}^{-RCT_{1}}$
$\displaystyle J$ $\displaystyle\equiv
4e^{2}T_{2}\frac{E_{C}E_{J}^{(0)}}{\hbar^{2}}\operatorname{e}^{-RCT_{1}}\left(\frac{1-\operatorname{e}^{-RCT_{1}}}{RC}\right)$
$\displaystyle K$ $\displaystyle\equiv
4e^{2}T_{2}\frac{E_{C}E_{K}^{(0)}}{\hbar^{2}}\operatorname{e}^{-RCT_{1}}\left(\frac{1-\operatorname{e}^{-RCT_{1}}}{RC}\right)$
(26)

we get back the dissipative mapping Eq. (II).

There is an important consideration about the charging energy per site, given
by $E_{q}(n)=\frac{E_{C}}{2L}\sum_{j=1}^{L}[q_{j}(t_{n+1})]^{2}$. This
quantity is proportional to the kinetic energy per site we have discussed in
Sec. VII, and we have seen there that in the dissipative model it attains an
asymptotic value order 1. By appropriately tuning the parameters, there is the
possibility to keep this asymptotic value smaller than the superconducting gap
of the system. In this way, the model obeys for long times the dynamics we
have described here. In contrast with that, in case of Hamiltonian evolution,
the energy at some point starts increasing in an unbounded way [31, 26, 27,
31] and at some point it goes beyond the gap and superconductivity is lost.

##  IX Conclusion

In conclusion we have studied a model of coupled kicked rotors with
dissipation. In the case of a single rotor, the model reduces to the Zaslavsky
map and shows a behavior strictly similar to the period-doubling route to
chaos, that is so widespread in Nature. (Actually we see two parallel period-
doubling cascades – see Fig. 8.

We focus on the momenta probed at discrete stroboscopic times (at each period
of the driving) and consider the appropriate averages over many random initial
conditions. Using some period $m$-tupling order parameters inspired to the
time-crystal literature, we see that in the many-rotor case the period-
doubling cascade disappears and one can essentially see only period doubling,
with here and there much smaller contributions at period four times the
driving. So, we find essentially confirmation of the findings of Ref. [14],
with the presence of the small period-doubling contributions probably related
to the fact that this model has local continuous variables, in contrast with
the discrete ones of cellular automata and classical time crystals.

The dynamics of this model is very rich. First of all we see that the period
doubling always appears in association with the spontaneous formation of
patterns, that are stable and persistent in time, and depend on the
stroboscopic time with periodicity double than the one of the driving.
Therefore the system breaks at the same time the discrete time translation
symmetry and the space translation symmetry, therefore we talk about
“spatiotemporal ordering”. A similar phenomenon occurs in quantum Floquet time
crystals, but the physics is different: Here an effect of the classical
nonlinear dynamics, there a sort of quantum phase transition. Indeed, here we
see a persistent spatiotemporal ordering already at finite sizes, while in
time crystals it persists only in the thermodynamic limit.

Most remarkably, the spatiotemporal order can be predicted analytically, by
means of a linear stability analysis around the state with all vanishing
momenta and coordinates. This analysis predicts an instability region, whose
boundary coincides precisely with the lower boundary of the spatiotemporally
ordered phase. Moreover the eigenvalues of the stroboscopic dynamics are real
and negative, so that the most unstable mode grows changing sign at any
period, until this growth stops due to nonlinearities, and a period-doubled
steady pattern develops. Another interesting finding is that the momentum of
the most unstable mode is $k=\pi$, giving rise to a typical length scale
double than the lattice spacing. We find that this prediction of the linear-
stability analysis is valid for the fully-developed patterns in the
spatiotemporally ordered regime, that have precisely this typical length
scale. Therefore period doubling in space is associated to period doubling in
time.

There is also a phase where the system breaks only the space-translation
symmetry and the patterns are stable and independent from the stroboscopic
time. Inside this region of the phase space there is a smaller region that we
call “Weak patterning” and where points with no patterns alternate with point
with small-amplitude patterns, in a jagged and apparently fractal way.

Beyond that, there is the trivial regime, where the model relaxes to a uniform
condition with vanishing momentum, and the chaotic regime, where the largest
Lyapunov exponent is larger than 0, nearby trajectories in phase space diverge
exponentially, and the dynamics is aperiodic and irregular in space and time.

The transition from different regimes can be seen in the properties of the
patterns. For instance, the pattern typical length has a sudden drop moving
from time-independent patterning to spatiotemporal ordering, and – in a range
of parameters – shows a peak at the onset of chaos where it increases of one
order of magnitude. Inside the spatiotemporally ordered phase, near the lower
phase boundary, we find also a typical length scale near the value $\lambda=2$
predicted by linear stability analysis.

Many transitions in the dynamics can be seen also in the behavior of the
pattern amplitude, that anyway misses the threshold between time-independent
patterning and spatiotemporal ordering, because changes in a continuous and
regular way at this threshold. The same information given by the amplitude can
be understood from the behavior of the asymptotic average kinetic energy per
site. This quantity is finite, also in the chaotic regime, in contrast with
the unbounded increase of energy in case of Hamiltonian chaos. Moreover, in
our case, this quantity does not scale with the system size.

This is an important information for experimental realization. Indeed, we
propose to realize this model in an array of SQUID Josephson junctions, and
the fact that the system does not heat up above a threshold gives the
opportunity to choose the parameters in such a way that the system stays
superconducting, and our model is a good description thereof for long times.

About future developments, it is first of all compelling to understand why
linear stability analysis does not predict the time-independent patterns.
Another possibility is studying the patterns in different geometries (for
instance in a 2-dimensional lattice) and see if similar dynamical behaviors
appear if one couples other systems showing the period doubling cascade (for
instance the nonlinear electric circuits of [10, 11]). One might think also to
use the methods of [33] to quantize the model of coupled dissipative rotors
and, of course, realize the experimental proposal of Sec. VIII.

###### Acknowledgements.

I thank A. Delmonte, R. Fazio, D. Mukamel, G. Passarelli, S. Ruffo, and G. E.
Santoro for interesting discussions, and V. Russomanno for having drawn my
attention on Ref. [4] and the problems of chaos and period-doubling cascade. I
acknowledge P. Lucignano for the access to the qmat machine where most of the
numerics for this project was performed, and thank the ICTP for the warm
hospitality received (under ERC Project 101053159 – RAVE) during the
completion of this work. I acknowledge financial support from PNRR MUR Project
PE0000023-NQSTI.

##  Appendix A Validiy of the Hamiltonian Eq. (22) in the case without
resistance

In order to show the validity of the Hamiltonian Eq. (22), let us write the
Kirchhoff’s equations for the circuit in Fig. 12 without resistance. In order
to do that, let us call $I_{j}$ the current moving from the island $j$ to the
island $j+1$, $I_{j-1}$ the current from the island $j-1$ to the island $j$,
$I_{S\,j}$ the current from the superconducting island $j$ to the ground (see
Fig. 13).

---
Figure 13: Currents at the node where the island $j$ is.

Using the first Kirchhoff law at the node of the island we get

$I_{j-1}=I_{j}+I_{S\,j}\,.$

Using the definition of Josephson current we get
$I_{j-1}=\frac{2e}{\hbar}E_{J}(\Phi)\sin(\theta_{j}-\theta_{j-1})$,
$I_{j}=\frac{2e}{\hbar}E_{J}(\Phi)\sin(\theta_{j+1}-\theta_{j})$. Using the
Kirchhoff’s law at the bifurcation to the parallel of the Josephson junction
below and the capacitor, we get

$I_{S\,j}=-\dot{q}_{j}-\frac{2e}{\hbar}E_{K}(\Phi)\sin(\theta_{j})\,.$

Putting all these formulae together we get

$\dot{q}_{j}=\frac{2e}{\hbar}\left[E_{J}\sin(\theta_{j+1}-\theta_{j})-E_{J}\sin(\theta_{j}-\theta_{j+1})-E_{K}\sin\theta_{j}\right]\,.$
(27)

It is easy to check that this formula is provided by the canonical equation
$\dot{q}_{j}=-\frac{2e}{\hbar}\partial_{\theta_{j}}H$ [see Eq. (VIII)], with
$H$ given by Eq. (22).

We can also easily write the equations for the derivative of the phase.
Applying the known formula [47] at the junction below, we get

$\hbar\dot{\theta}_{j}=2e\Delta V_{j}=2e\frac{q_{j}}{C}\,,$

where $\Delta V_{j}$ is the voltage difference between the island $j$ and the
ground. The point is that $\theta_{j}$ is the phase difference across the
Josephson junction below, and the voltage difference across this junction is
provided by the voltage difference across the capacitor, $\Delta
V_{j}=\frac{q_{j}}{C}$. It is easy to show that this equation for the
derivative of $\theta_{j}$ is provided by the canonical equation
$\dot{\theta}_{j}=\frac{2e}{\hbar}\partial_{q_{j}}H$ [see Eqs. (VIII)], with
$H$ provided by Eq. (22).

## References

  * Ball [1998] P. Ball,  _The Self-Made Tapestry: Pattern Formation in Nature_ (Oxford Univ Press, Oxford, 1998).
  * Strogatz [2015] S. Strogatz,  _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering_ , $2^{\rm nd}$ ed. (Taylor & Francis, London, 2015).
  * Cross and Greenside [2009] M. Cross and H. Greenside,  _Pattern Formation and Dynamics in Nonequilibrium Systems_ (Cambridge University Press, Cambridge, 2009).
  * Nagashima and Baba [1999] H. Nagashima and Y. Baba,  _Introduction to Chaos: Physics and Mathematics of Chaotic Phenomena_ (Institute of Physics Publishing, Bristol and Philadelphia, 1999).
  * Goldstein  _et al._ [2014] H. Goldstein, C. P. Poole, and J. L. Safko,  _Classical Mechanics_ , new international edition ed. (Pearson, Edinburgh Gate, 2014).
  * Feigenbaum [1976] M. Feigenbaum, Universality in complex discrete dynamics (1975-1976), Los Alamos Theoretical Division Annual Report.
  * Briggs [1997] K. Briggs, Feigenbaum scaling in discrete dynamical systems (1997), PhD Thesis, University of Melbourne.
  * Giglio  _et al._ [1981] M. Giglio, S. Musazzi, and U. Perini, Transition to chaotic behavior via a reproducible sequence of period-doubling bifurcations, Phys. Rev. Lett. 47, 243 (1981).
  * Libchaber  _et al._ [1982] A. Libchaber, C. Laroche, and S. Fauve, Period doubling cascade in mercury, a quantitative measurement, J. Physique Lett. 43, 211 (1982).
  * Linsay [1981] P. S. Linsay, Period doubling and chaotic behavior in a driven anharmonic oscillator, Phys. Rev. Lett. 47, 1349 (1981).
  * Testa  _et al._ [1982] J. Testa, J. Pérez, and C. Jeffries, Evidence for universal chaotic behavior of a driven nonlinear oscillator, Phys. Rev. Lett. 48, 714 (1982).
  * Jia  _et al._ [2011] B. Jia, H. Gu, L. Li, and X. Zhao, Dynamics of period-doubling bifurcation to chaos in the spontaneous neural firing patterns, Cogn. Neurodyn. 6, 89 (2011).
  * Hartmann  _et al._ [2017] M. Hartmann, D. Poletti, M. Ivanchenko, S. Denisov, and P. Hänggi, Asymptotic floquet states of open quantum systems: the role of interaction, New Journal of Physics 19, 083011 (2017).
  * Bennett  _et al._ [1990] C. H. Bennett, G. Grinstein, Y. He, C. Jayaprakash, and D. Mukamel, Stability of temporally periodic states of classical many-body systems, Phys. Rev. A 41, 1932 (1990).
  * Zhuang  _et al._ [2021] Q. Zhuang, F. Machado, N. Y. Yao, and M. P. Zaletel, An absolutely stable open time crystal (2021), arXiv:2110.00585 [quant-ph] .
  * Pizzi  _et al._ [2021] A. Pizzi, A. Nunnenkamp, and J. Knolle, Bistability and time crystals in long-ranged directed percolation, Nature Communications 12, 10.1038/s41467-021-21259-4 (2021).
  * Gambetta  _et al._ [2019] F. M. Gambetta, F. Carollo, A. Lazarides, I. Lesanovsky, and J. P. Garrahan, Classical stochastic discrete time crystals, Phys. Rev. E 100, 060105(R) (2019).
  * [18] See also [51] for a long-lasting metastable period doubling in a system with continuous onsite variables.
  * Sacha and Zakrzewski [2017] K. Sacha and J. Zakrzewski, Time crystals: a review, Reports on Progress in Physics 81, 016401 (2017).
  * Else  _et al._ [2020] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Discrete time crystals, Annual Review of Condensed Matter Physics 11, 467 (2020), https://doi.org/10.1146/annurev-conmatphys-031119-050658 .
  * Zaletel  _et al._ [2023] M. P. Zaletel, M. Lukin, C. Monroe, C. Nayak, F. Wilczek, and N. Y. Yao, Colloquium: Quantum and classical discrete time crystals (2023), arXiv:2305.08904 [quant-ph] .
  * Kaneko [1984] K. Kaneko, Prog. Theor. Phys. 72, 480 (1984).
  * Kapral [1985] R. Kapral, Pattern formation in two-dimensional arrays of coupled, discrete-time oscillators, Phys. Rev. A 31, 3868 (1985).
  * Kaneko [1993] K. Kaneko,  _Theory and Applications of Coupled Map Lattices_ (Wiley, Chichester, 1993)).
  * Notarnicola  _et al._ [2018] S. Notarnicola, F. Iemini, D. Rossini, R. Fazio, A. Silva, and A. Russomanno, From localization to anomalous diffusion in the dynamics of coupled kicked rotors, Phys. Rev. E 97, 022202 (2018).
  * Kaneko and Konishi [1989] K. Kaneko and T. Konishi, Diffusion in hamiltonian dynamical systems with many degrees of freedom, Phys. Rev. A 40, 6130 (1989).
  * Konishi and Kaneko [1990] T. Konishi and K. Kaneko, Diffusion in hamiltonian chaos and its size dependence, Journal of Physics A: Mathematical and General 23, L715 (1990).
  * Falcioni  _et al._ [1991] M. Falcioni, U. M. B. Marconi, and A. Vulpiani, Ergodic properties of high-dimensional symplectic maps, Phys. Rev. A 44, 2263 (1991).
  * Chirikov and Vecheslavov [1993] B. V. Chirikov and V. V. Vecheslavov, Theory of fast arnold diffusion in many-frequency systems, J. Stat. Phys. 71, 243 (1993).
  * Chirikov and Vecheslavov [1997] B. V. Chirikov and V. V. Vecheslavov, Arnold diffusion in large systems, J. Exp. Theor. Phys. 85, 616 (1997).
  * Rajak  _et al._ [2019] A. Rajak, I. Dana, and E. G. Dalla Torre, Characterizations of prethermal states in periodically driven many-body systems with unbounded chaotic diffusion, Phys. Rev. B 100, 100302(R) (2019).
  * Zaslavsky [1978] G. M. Zaslavsky, The simplest case of a strange attractor, Physics Letters A 69, 145 (1978).
  * Dittrich and Graham [1986] T. Dittrich and R. Graham, Quantization of the kicked rotator with dissipation, Z. Phys. B 62, 515 (1986).
  * Else  _et al._ [2016] D. V. Else, B. Bauer, and C. Nayak, Floquet time crystals, Physical Review Letters 117, 090402 (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).
  * von Keyserlingk  _et al._ [2016] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, Absolute stability and spatiotemporal long-range order in floquet systems, Physical Review B 94, 085112 (2016).
  * Khemani  _et al._ [2016] V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phase structure of driven quantum systems, Physical Review Letters 116, 250401 (2016).
  * Chirikov [1991] B. V. Chirikov, Time-dependent quantum systems, in  _Chaos and Quantum Mechanics_ , Les Houches Lecture Series, Vol. 52, edited by M.-J. Giannoni, A. Voros, and J. Zinn-Justin (Elsevier, Amsterdam, 1991) pp. 443–545.
  * Cataliotti  _et al._ [2001] F. S. Cataliotti, S. Burger, C. Fort, P. Maddaloni, F. Minardi, A. Trombettoni, A. Smerzi, and M. Inguscio, Josephson junction arrays with bose-einstein condensates, Science 293, 843 (2001), https://www.science.org/doi/pdf/10.1126/science.1062612 .
  * Stöckmann [2007] H.-J. Stöckmann,  _QUantum Chaos: An Introduction_ (Cambridge University Press, Cambridge, UK, 2007).
  * Berry [1978] M. V. Berry, Regular and Irregular Motion, in  _Topics in Nonlinear Mechanics_, Vol. 46, edited by S. Jorna (Am.Inst.Ph., 1978) pp. 16–120.
  * Pikovsky and Politi [2016] A. Pikovsky and A. Politi,  _Lyapunov Exponents: A Tool to Explore Complex Dynamics_ (Cambridge University Press, Cambridge, 2016).
  * Ott [2002] E. Ott,  _Chaos in Dynamical Systems_ , $2^{\rm nd}$ ed. (Cambridge University Press, Cambridge, 2002).
  * Benettin  _et al._ [1976] G. Benettin, L. Galgani, and J.-M. Strelcyn, Kolmogorov entropy and numerical experiments, Phys. Rev. A 14, 2338 (1976).
  * Berge  _et al._ [1987] P. Berge, C. Vidal, and Y. Pomeau,  _Order within Chaos: Towards a Deterministic Approach to Turbulence_ (John Wiley and Sons, New York, 1987).
  * Leggett [1987] A. J. Leggett, Quantum Mechanics at the Macroscopic Level, in  _Le hasard et la matière/Chance and matter_ , Les Houches, Session XLVI, 1986, edited by J. Souletie and J. Vannimenus and R. Stora (Elsevier Science, 1987).
  * Tinkham [1996] M. Tinkham,  _Introduction to Superconductivity_ , $2^{\text{nd}}$ ed. (Mc Graw-Hill, 1996).
  * Devoret [1997] M. H. Devoret, Quantum Fluctuations in Electrical Circuits, in  _Fluctuations Quantiques/Quantum Fluctuations_, Les Houches, Session LXIII, 1995, edited by S. Reynaud ND E. Giacobino and J. Zinn-Justin (Elsevier Science, 1997).
  * Pino  _et al._ [2016] M. Pino, L. B. Ioffe, and B. L. Altshuler, Nonergodic metallic and insulating phases of josephson junction chains, Proceedings of the National Academy of Sciences 113, 536 (2016), https://www.pnas.org/content/113/3/536.full.pdf .
  * Mithun  _et al._ [2019] T. Mithun, C. Danieli, Y. Kati, and S. Flach, Dynamical glass and ergodization times in classical josephson junction chains, Phys. Rev. Lett. 122, 054102 (2019).
  * Yao  _et al._ [2020] N. Y. Yao, C. Nayak, L. Balents, and M. P. Zaletel, Classical discrete time crystals, Nature Physics 16, 438 (2020).