Bistability and time crystals in long-ranged directed percolation

Stochastic processes govern the time evolution of a huge variety of realistic systems throughout the sciences. A minimal description of noisy many-particle systems within a Markovian picture and with a notion of spatial dimension is given by probabilistic cellular automata, which typically feature time-independent and short-ranged update rules. Here, we propose a simple cellular automaton with power-law interactions that gives rise to a bistable phase of long-ranged directed percolation whose long-time behaviour is not only dictated by the system dynamics, but also by the initial conditions. In the presence of a periodic modulation of the update rules, we find that the system responds with a period larger than that of the modulation for an exponentially (in system size) long time. This breaking of discrete time translation symmetry of the underlying dynamics is enabled by a self-correcting mechanism of the long-ranged interactions which compensates noise-induced imperfections. Our work thus provides a firm example of a classical discrete time crystal phase of matter and paves the way for the study of novel non-equilibrium phases in the unexplored field of driven probabilistic cellular automata.

P ercolation theory describes the connectivity of networks, with applications pervading virtually any branch of science 1 , including economics 2 , engineering 3 , neurosciences 4 , social sciences 5 , geoscience 6 , food science 7 and, most prominently, epidemiology 8 . Among the multitude of phenomena described by percolation, of predominant importance are spreading processes, in which time plays a crucial role and that can be studied within models of directed percolation (DP) 9 . Characterized by universal scalings in time 10 , in their discretized versions these models are probabilistic cellular automata (PCA), that is, dynamical systems with a state evolving in discrete time according to a set of stochastic and generally short-ranged update rules. To account for certain realistic situations, e.g. of long-distance travels in epidemic spreading, DP has been extended to long-ranged updates 11,12 leading to a change of the universal scaling exponents 13 .
Despite their wide applicability, PCAs have surprisingly remained an outlier in a branch of non-equilibrium physics that has recently experienced a tremendous amount of excitementthat of discrete time crystals (DTCs) [14][15][16][17][18][19][20] . In essence, DTCs are systems that, under the action of a time-periodic modulation with period T, exhibit a periodic response at a different period T 0 ≠T, thus breaking the discrete time-translational symmetry of the drive and of the equations of motion. DTCs thus extend the fundamental idea of symmetry breaking 21 to non-equilibrium phases of matter. Following the pioneering proposals in the context of many-body-localized (MBL) systems 17,18 , DTCs have been observed experimentally 22,23 , and their notion has been extended beyond MBL [24][25][26][27] .
More recently, Yao and collaborators have fleshed out the essential ingredients of a classical DTC phase of matter 28 . Namely, in a classical DTC, many-body interactions should allow for an infinite autocorrelation time, which should be stable in the presence of a noisy environment at finite temperature, a subtle requirement that rules out the vast class of long-known deterministic dynamical systems. Despite various efforts [28][29][30][31] , an example of such a classical DTC has mostly remained elusive, and proving an infinite autocorrelation time robust to noise and perturbations for this phase of matter is an outstanding problem. The general expectation is in fact that PCAs and other minimal models for noisy systems in one spatial dimension can only show a transient subharmonic response because noise-induced imperfections generically nucleate and spread, destroying true infiniterange symmetry breaking in time 28,32 .
Here we overcome these difficulties by introducing a simple and natural generalization of DP in which the dynamical rules are governed by power-law correlations. This leads to qualitative changes of the system behaviour and, crucially, the emergence of a bistable phase of long-ranged DP, enabled by the ability of longrange interactions to counteract the dynamic proliferation of defects. By adding a periodic modulation to the update rules, we then study a version of periodically driven DP and show that the underlying bistable phase intimately connects to a stable DTC. In this non-equilibrium phase, the system is able to self-correct noise-induced errors and the autocorrelation time grows exponentially with the system size, thus becoming infinite in the thermodynamic limit. In analogy to the one-dimensional Ising model for which, at equilibrium, long-range interactions enable a normally forbidden finite-temperature magnetic phase 33,34 , in our model, out of equilibrium, the long-range interactions lead to a classical time-crystalline phase. Crucially, our results appear naturally in a minimal model of long-ranged DP but are expected to find applications in many different contexts of dynamical many-body systems.
Basic understanding of new concepts has historically been built around the study of minimal models, such as the Ising model for magnetism at equilibrium 33,34 , the kicked transverse field Ising chain for DTCs 17,18 , or the prototypical Domany-Kinzel (DK) PCA for DP 35 . In this paper, we start our discussion with a brief review of the DK model and then generalize it to include power-law interactions. We characterize its phase diagram and show that its long-range nature is the key ingredient for the emergence of a bistable phase. Finally, we include a periodic drive for the long-ranged DP process and show with a careful scaling analysis that the autocorrelation time of the subharmonic response is exponential in system size. In the thermodynamic limit, our model provides therefore the first example of a PCA behaving as a classical DTC, which is persistent and stable to the continuous presence of noise. Lastly, we conclude with a summary of our findings and an outlook for future research.

Results
Review of DP. We consider a triangular lattice in which one dimension can be interpreted as discrete space i and the other one as discrete time t = 1, 2, 3, … , see Fig. 1. To implicitly account for the triangular nature of the lattice, i runs over integers and halfintegers at odd and even times t, respectively. We denote L the spatial system size and are interested in the thermodynamic limit L → ∞. The site i at time t can be either occupied or empty, s i,t = 0,1. For a given time t, we call generation the collection of variables fs i;t g i specifying the system state. Initially, the sites are occupied with uniform probability p 1 > 0. A DP process is defined by a stochastic Markovian update rule with which, starting from the initial generation fs i;1 g i , all subsequent generations fs i;t g i are obtained one by one. The main observable we will focus on is the global density n(t) (henceforth just referred to as density for brevity) defined as where the inner and outer brackets denote average over the L sites and over R independent runs, respectively. Since n(1) = p 1 , we will often refer to p 1 as initial density. The simplest, and yet already remarkably rich, example of the above setting of DP is the DK model 35 . Here, we briefly review it adopting an unconventional notation that, making explicit use of a local density, will prove very convenient for a straightforward generalization to a model of long-ranged DP.
In the DK model, the probability of site i to be occupied at time t depends on the state of its neighbours i ± 1/2 at previous time t − 1. More specifically, as summarized in Fig. 1a, site i is: (i) empty if both its neighbours were empty, (ii) occupied with probability q 1 if one and just one of its neighbours was occupied, and (iii) occupied with probability q 2 if both its neighbours were occupied. To account for these possibilities in a compact fashion, we define a local density n i,t as and say that site i at time t is occupied with a probability p i,t given by In other words, the probability p i,t is a nonlinear function f q 1 ;q 2 ðn i;t Þ of the local density n i,t , with domain {0, 0.5, 1}. Since n i, t only involves the nearest neighbours of site i, the DK model of DP is obviously 'short-ranged'. In essence, s i,t is a Bernoullian random variable of parameter p i,t , which we compactly denote s i;t $ Bernoulli ðp i;t Þ. The complexity of this model arises from the fact that the value of the parameter p i,t is not known a priori, as it depends on the actual state of the system at previous time t − 1. Equipped with a random number generator, one can obtain all the generations one by one according to the above procedure, as schematically illustrated in the flowchart of Fig. 1b. Reiterating for several independent runs, one finally obtains the time series of the density n in Eq. (1). The DK model features two dynamical phases, shown in Fig. 1c, d. In the inactive phase, for small enough probabilities q 1 and q 2 , the system eventually reaches the completely unoccupied absorbing state, that is, no percolation occurs. In the active phase instead, for large enough probabilities q 1 and q 2 , a finite fraction of sites remains occupied up to infinite time, that is, the system percolates. For small initial probability p 1 ≪ 1, the critical line separating the two phases is characterized by a power-law growth of the density 36 , n~t θ , with exponent θ ≈ 0.31. As conjectured by Grassberger 37 , this exponent is universal for all systems in the DP universality class. Indeed, DP exemplifies how the unifying concept of universality pertaining to quantum and classical many-body systems 38 can be extended to non-equilibrium phenomena.
Important for our work is that, in the DK model, whether the system percolates or not depends on the parameters q 1 and q 2 but not on the initial density p 1 , at least as long as p 1 > 0. Indeed, the phase boundaries for initial densities p 1 = 0.01 and p 1 = 1 in Fig. 1c, d, respectively, coincide.
Long-ranged percolation and bistability. As the vast majority of PCA, the DK model features short-ranged update rules 9 . In realistic systems, however, it is often the case that the occupation of a site i is influenced not only by the neighbouring sites but also by farther sites j, with an effect decreasing with the distance r i,j between the sites. Building on an analogy with the DK model, we propose here a model for such a 'long-ranged' DP, whose protocol is explained in the flowchart of Fig. 2a. Specifically, we consider as a local density n i,t a power-law-weighted average of the previous generation fs j;tÀ1 g j centred around site i where the normalization factor N α;L ensures n i,t = 1 if all sites j are occupied and the adjective 'local' emphasizes the site dependence. The occupation probability p i,t then depends on the local density n i,t through some nonlinear function f μ that for concreteness we consider to be with μ ∈ (0, 1) a control parameter. The whole DP dynamics is determined via the occupations s i;t $ Bernoulli ðp i;t Þ and reiterating from one generation to the next. Note, our findings are not contingent on the specific choice of Eqs. (4) and (5) but are rather expected to hold generally for a broad class of long-ranged forms of the densities n i,t and of functions f μ -see 'Methods' section for details. We emphasize that Eqs. (4) and (5) 2) and (3) and Fig. 1b, respectively. Furthermore, whereas in the DK model the control parameters are the probabilities q 1 and q 2 , the control parameter is now μ. As an important difference, now the domain of f μ accounts for several (and α-dependent) values of n i,t , for which the piecewise definition of p i,t as in Eq. (3) would have been unpractical, and the compact form of Eq. (5) was necessary instead.
The introduction of a long-ranged local density n i,t in Eq. (4) has profound implications. Arguably, the most dramatic is the appearance of a bistable phase, in addition to the standard active and inactive ones. In the bistable phase, the ability of the system to percolate depends on the initial density p 1 , see the red lines in Fig. 2b, c. That is, the bistable phase features two basins of attraction, resulting into an asymptotically vanishing or finite n, respectively, and separated by some critical initial density p 1,c > 0. To characterize systematically the dynamical phases of our model, we plot in Fig. 2d, e the long-time density n(t = 10 3 ) as a suitable order parameter in the plane of the power-law exponent α and control parameter μ. Comparing the results obtained for a large and a small initial density p 1 , it is possible to sketch a phase diagram composed of three phases: (i) inactive-n decays to 0 at long times; (ii) active-n does not decay at long times; (iii) bistable-n either decays or not depending on p 1 being small or large. The existence of this bistable phase is in striking contrast with short-ranged models of DP such as the DK model and in fact appears only for α ⪅ 2, that is, when the local densities fn i;t g i are The probability p i,t of site i to be occupied at time t depends on the occupation of its nearest-neighbours i ± 1 2 at time t − 1 and can take discrete values 0, q 1 and q 2 . b Flowchart representation of the DK model. The initial occupation probability is uniform p i,t = 1 = p 1 . At time t, each site i is either occupied (s i,t = 1) or empty (s i,t = 0) with probability p i,t and 1 − p i,t , respectively. Time is advanced and local densities fn i;t g i are computed for each site i as averages of the nearest-neighbour occupations at previous time, and these densities determine the occupation probabilities for the next generation, see Eq. (3). The generations at all subsequent times are obtained by iteration. c, d The density n at late times can be used to discern the active and inactive phases, in which n(t = 10 3 ) >0 and ≈0, respectively. The dashed lines serve as a reference to locate the phase boundary and are the same for initial densities p 1 = 1 (c) and p 1 = 0.01 (d). The insets show representative single instances of the DP for the points in the (q 1 , q 2 ) plane marked with a cross. Here L = 100 and R = 10 3 . correlated over a sufficiently long range. To understand the origin of this rich phenomenology, we study the short-and infiniterange limits of our DP process.
In the short-range limit α → ∞, the local densities n i,t reduce to the averages of the nearest-neighbour occupations s iÀ 1 2 ;tÀ1 and s iþ 1 2 ;tÀ1 , that is, Eq. (4) recasts into Eq. (2) and the DK model is recovered. In the notation of Eq. (3), the DK parameters are q 1 = f μ (0.5) and q 2 = f μ (1). Therefore, we can move across the DK parameter space (q 1 ,q 2 ) varying μ, going from the inactive phase (μ < μ 1 c ) to the active one (μ > μ 1 c ), and no bistable phase is possible. We find that the transition happens at a critical μ 1 c ¼ 0:85ð7Þ. Note that, in the active phase, a completely empty state (p 1 = 0) remains trivially empty at all times. This behaviour is, however, unstable, because any p 1 > 0 leads to percolation (i.e. p 1,c = 0), and we therefore do not classify the active phase as bistable. At criticality, and for p 1 ≪ 1, the density grows as n~t θ with θ = 0.3(0), as expected for the DP universality class 9 . See Supplementary Fig. 2 for details.
In the infinite-range limit α → 0, and more generally for α ≤ 1, the factor N α;L in Eq. (4) diverges as L → ∞. Correspondingly, spatial stochastic fluctuations are suppressed, that is, all sites i share the same occupation probability p i,t+1 = p t and density n i,t = n(t) = p t . Therefore, in this limit the dynamics reduces to the deterministic 0-dimensional recurrence relation The system asymptotic behaviour can then be understood from the analysis of the fixed points (FPs) of the equation x = f μ (x), which is detailed in the 'Methods' section.
Driven percolation and time crystals. We have established that long-range correlated local densities fn i;t g i give rise to a bistable phase. We now show how, in a driven DP with periodically modulated update rules, this phase intimately relates to the emergence of a classical DTC. In this phase, as we shall see, the density n displays oscillations over a period larger than that of the drive and up to a time that, thanks to the long-range interactions and despite the presence of multiple sources of noise, is exponentially large in the system size, a feature that would generally be forbidden in short-ranged PCA 28 . In the thermodynamic limit L → ∞, these subharmonic oscillations are therefore persistent, that is, the system autocorrelation time diverges to infinity, breaking the time-translational symmetry and proving a classical DTC in a periodically driven PCA.
In the spirit of keeping the model as simple as possible, we consider a minimal drive in which, after every T iterations of the  4) and (5). b, c Time evolution of the density n for p 1 = 1 (b) and p 1 = 0.01 (c) for various representative values of the power-law exponent α and control parameter μ. Three dynamical phases can be distinguished: (i) inactive-the density n decays to 0 (blue); (ii) active-n does not decay to 0 (yellow); (iii) bistable-n either decays to 0 or not depending on whether the initial density p 1 is small or large (red). d, e Long-time density n(t = 10 3 ) in the plane of α and μ for p 1 = 1 (d) and p 1 = 0.01 (e). With the criterion used in b, c, we discern the three phases: inactive (light), active (dark), and bistable (light or dark depending on p 1 ). The dashed lines help locating the phases and coincide in d and e, and critical values μ 0 c and μ 1 c of μ in the limits α → 0 and α → ∞, respectively, are reported (the offset of μ 0 c from the dashed line, as well as the softening of the dashed line for α ≈ 1, are due to finite-size effects). Crucially, the bistable phase is present only for small enough α ≾ 2, that is, for a sufficiently longranged DP. Single instances of the DP for the three phases are shown in the insets, as obtained for the α and μ indicated with coloured dots, and corresponding to the parameters used in b, c. Here R = 10 4 and 10 2 in b, c and d, e, respectively, and L = 500. DP in Eqs. (4) and (5), empty sites are turned into occupied ones and vice versa, making the full equations of motion periodic with period T. As a further source of imperfections, adding to the underlying noisy DP, we also account for faulty swaps with probability p d . More explicitly, the periodic drive consists of the following transformation s i;1þkT ! 1 À s i;1þkT with probability 1 À p d s i;1þkT with probability p d : ( In Fig. 3a, b, we show the spatio-temporal pattern of single instances of the driven DP, alongside with the density n averaged over several independent runs. If the DP is shortranged enough, the spatio-temporal pattern at long times looks similar from one period to the next, that is, the density n synchronizes with the drive and eventually picks a periodicity T. On the contrary, for a long-ranged enough DP, the system keeps alternating at every period between a densely occupied regime and a sparsely occupied one, and n oscillates with period 2T, that is, the system breaks the discrete time-translation symmetry of the equations of motion. When using the tag 'classical DTC', special care should be reserved for showing the defining features of this phase, namely, its rigidity and persistence 28 . Our system is rigid in the sense that it does not rely on fine-tuned model parameters, e.g. μ, α or the initial density p 1 , and that noise, either in the form of the inherently stochastic underlying DP or of a small but non-zero drive defect density p d , does not qualitatively change the results. Moreover, in the limit L → ∞, our DTC is truly persistent. Indeed, one might expect that the accumulation of stochastic mistakes introduces phase slips and eventually leads to the (possibly slow but unavoidable) destruction of the subharmonic response. Although this expectation is generally correct for shortranged DP models, including our model at large α, it can fail for long-ranged DP models.
To show that, in the limit L → ∞, the lifetime of our DTC is infinite, we perform a scaling analysis comparing results for increasing system sizes L. First, we introduce an order parameter Φ(t), henceforth called subharmonicity, that is defined at stroboscopic times t = 1, 1 + T, 1 + 2T, … as If the density n oscillates with the same period T as the drive, then n(t) = n(t + T) and Φ(t) = 0. On the contrary, if n oscillates with a doubled period 2T, then n(t = 1 + kT) is positive and negative for even and odd k, respectively, and Φ(t) is finite and maintains a constant sign. Therefore, Φ(t) is a suitable diagnostics to track the degree of subharmonicity of n in time and to perform the scaling analysis.
In Fig. 3c, we show Φ(t) for various system sizes L. For both α = 1.4 and α = 1.8, the subharmonicity decays exponentially in time, ΦðtÞ $ expðÀ tÀ1 τT Þ. As shown in Fig. 3d, these two values of α are, however, crucially different in how the lifetime τT scales with the system size. In fact, τT is approximately independent of L for α = 1.8, whereas it scales exponentially as τ $ expðβLÞ for α = 1.4, for which the decay of the subharmonicity is therefore just a finite-size effect. The scaling coefficient β quantifies the time crystallinity of the system and can thus be used to obtain a full phase diagram as a function of the power-law exponent α, in Fig. 3e. We observe a phase transition between a DTC and a trivial phase at α ≈ 1.7. That is, if the DP is sufficiently longranged (α ⪅ 1.7), β is finite and in the thermodynamic limit L → ∞ the subharmonic response extends up to infinite time, as required for a true DTC. In contrast, for a shorter-range DP (α ⪆ 1.7), β ≈ 0 independently of L and the subharmonic response is always dynamically destroyed.

Discussion
We have shown that long-range DP and its periodically driven variant can give rise to a bistable phase and a DTC, respectively. At the core of our model in Eqs. (4) and (5) is the idea that the occupation of a given site depends on the state of all the other sites at the previous time. In this sense, our model is reminiscent of some SIR-type models of epidemic spreading in which not only a sick site can infect a susceptible site, but several infected sites can also cooperate to weaken a susceptible site and finally infect it 39,40 . This cooperation mechanism among an infinite number of parent sites, rather than a finite one as considered in previous works on long-ranged DP 13,41 , is the key feature allowing the emergence of the bistable phase that finds a transparent explanation in the infinite-range limit α → 0, where it corresponds to the equation x = f μ (x) having two stable FPs. Bistability also provides intuition on the origin of the DTC, to which it is deeply connected. Indeed, the drive in Eq. (7) switches the system from a densely occupied regime to a sparsely occupied one (and vice versa). If the underlying DP is bistable, these regimes fall each within different basins of attraction and can therefore be both stabilized by the contractive dynamics 25,29 . Ultimately, this double stabilization facilitates the establishment of the DTC with infinite autocorrelation time. Remarkably, this mechanism does not rely on the equations of motion being perfectly periodic, as required for DTCs in closed MBL systems 42 , and we expect that infinite autocorrelation times could be maintained even in the presence of aperiodic variations of the drive (although the nomenclature should be revised in this case, since the underlying discrete time symmetry would only be present on average but not for individual realizations). This is in contrast to DTCs in closed MBL systems 42 , in which the non-ergodic dynamics hinges on the peculiar mathematical structure of the Floquet operator, which, in turns, relies on the underlying equations being perfectly periodic.
The intimate connection between bistability and DTC is, however, not a strict duality, and the boundaries of the two phases, in the equilibrium and non-equilibrium phase diagrams, respectively, do not coincide. For instance, in our analysis we found that for μ = 0.9 the bistable phase extends up to α ≈ 1.6, whereas the DTC stretches slightly farther, up to α ≈ 1.7. The origins of this imperfect correspondence can be traced back to two competing effects. On the one hand, bistability may not be sufficient to stabilize a DTC. This can already be understood in the limit α → 0, in which the asymmetry of f μ and of its FPs does not guarantee the drive to switch the density n from one basin of attraction to the other, that is, across the critical probability p 1,c . This issue becomes even more relevant for larger α, for which the asymmetry is possibly accentuated and p 1,c can approach 0 (see for instance Supplementary Fig. 1). On the other hand, a perfect bistability may not even be necessary for a DTC to exist. In fact, for the stabilization of a DTC, it may be sufficient that, of the densely and sparsely occupied regimes of the underlying DP, only one is stable, and the other is just weakly unstable (that is, metastable), meaning that the time scales of the dynamics of the density n in the two regimes are very different. Loosely speaking, the stability of one regime might be able to compensate for the weaker instability of the other, resulting in an overall stable DTC. The asymmetry of the underlying DP and the mismatch between the bistable phase and the DTC highlight the purely dynamical nature of the latter, that cannot 'piggy-back' on any underlying symmetry.
While these considerations are model and parameters dependent, and it is ultimately up to numerics to find the bistable and the DTC phases, what is universal and far reaching here is the concept that long-ranged DP, and PCA more generally, can host novel dynamical phases, such as DTCs. As Yao and collaborators recently pointed out 28 , long autocorrelation times are in fact generally unexpected in 1 + 1-dimensional PCA, because imperfections and phase slips can nucleate, spread and destroy the order. Our work proves that this fate can be avoided, and time-crystalline order established, in long-ranged PCA. These systems enable in fact an error correction mechanism, in our case intimately related to the bistability, that would be impossible if correlations were limited to a finite radius. We may speculate that, in the physical picture of a Hamiltonian system coupled to a bath, this defect suppression would correspond to the cooling rate being larger than the heating rate.
In conclusion, we have studied the effects of long-range correlated update rules in a model of DP, which we built from an analogy with the prototypical (but short-ranged) DK PCA. First, we proved that, beyond the standard active and inactive phases, a new bistable phase emerges in which the system at long times is either empty or finitely occupied depending on whether it was initially sparsely or densely occupied. Second, in a driven DP with periodic modulation of the update rules, we showed that this bistable phase intimately connects with a DTC phase, in which the density oscillates with a period twice that of the drive. In this DTC phase, the autocorrelation time scales exponentially with the system size, and in the thermodynamic limit a robust and persistent breaking of the discrete time-translation symmetry is established.
As an outlook for future research, further work on the driven DP should better assess the nature of the transition between the DTC Single instances of the periodically driven DP, alongside with the density n averaged over multiple independent runs, for L = 500 sites. a For a power-law exponent α = 1.4, n oscillates subharmonically with a period that is twice that of the drive, whereas, for α = 1.8, n eventually picks the periodicity T enforced by the drive. c For finite system sizes L, the subharmonicity Φ(t) decays as ΦðtÞ $ expðÀ tÀ1 τT Þ due to the accumulation of phase slips, and, after a few time scales τT, the density n synchronizes with the drive and oscillates with period T. Exponential fits (dotted lines) can be used to extrapolate the lifetime τ of the subharmonic response, on which a scaling analysis is performed in d. For α = 1.4 (blue), the lifetime τ scales exponentially with the system size, τ~e βL , whereas no such a scaling is found for α = 1.8. The scaling coefficient is again found from an exponential fit (dotted line) and plotted in e versus the power-law exponent α. For small α, that is, long-ranged enough DP, the scaling coefficient β is finite, indicating that in the thermodynamic limit L → ∞ the subharmonic response is persistent and a DTC with infinite autocorrelation time emerges. On the contrary, β ≈ 0 for large α, indicating a trivial dynamical phase in which no stable subharmonic dynamics is established. Here we considered p 1 = 1, μ = 0.9, p d = 0.02, T = 20 and R = 2000. and the trivial phase, characterize more systematically the phase diagram in other directions of the parameter space, and, most interestingly, address the role of dimensionality. Indeed, it is well known that dimensionality can facilitate the establishment of ordered phases of matter at equilibrium, and the question whether this is the case also out of equilibrium remains open. A positive answer to this question is suggested by the fact that, in D + 1dimension with D ≥ 2, bistability can emerge even in short-ranged models of DP 40,43,44 . Another interesting question regards the fate of chaos and damage spreading in long-ranged DP 45 . Further research should then aim to gain analytical intuition into the problem. For instance, the critical α separating the various phases may be located using a field theoretical approach, which has been successful in similar contexts in the past 41 . Finally, on a broader perspective, our work paves the way towards the study of nonequilibrium phases of matter in the uncharted territory of driven PCA, with a potentially very broad range of applications throughout different branches of science. As a timely example, Floquet PCA may provide new insights into the understanding of seasonal epidemic spreading and periodic intervention efficacy.

Methods
Here we provide further technical details on our work. In Eq. (4), we considered as distance r i,j between sites i and j where the tangent accounts for periodic boundary conditions and makes the distance of the farthest sites with |i − j| = L/2 artificially diverge. This divergence is expected to reduce finite-size effects without changing the underlying physics, that is, in fact dominated by sites with |i − j| ≪ L, for which we get a natural r i,j ≈ |i − j|. Indeed, as we checked, similar results are obtained with r i;j ¼ minðji À jj; L À ji À jjÞ. The Kac-like normalization factor N α;L reads instead The phenomenology of the bistable phase can be understood from a graphical FP analysis of the equation f μ (x) = x illustrated in Fig. 4, which explains the dynamics for α < 1. Three scenarios are possible and interpreted in terms of the ways the graph of the function f μ intersects with the bisect. (i) Inactive-if μ < μ 0 c , the only FP is x 0 = 0, which is stable and corresponds to a completely empty state. The system moves towards this FP and p t t ! 1!0. (ii) Critical-if μ ¼ μ 0 c , a new semi-stable FP emerges at x c , which is attractive from its right and repulsive on its left. (iii) Bistable-if μ > μ 0 c , the semi-stable FP splits into an unstable FP x 1 > x 0 and a stable FP x 2 > x 1 . In this case, the system will reach either the unoccupied FP x 0 = 0 or the finitely occupied FP x 2 > 0 depending whether p 1 < x 1 or p 1 > x 1 , respectively. That is, the system is bistable, and the critical initial probability separating its two basins of attraction is p 1,c = x 1 (see also Supplementary Fig. 1). The critical value μ 0 c is obtained numerically solving for the condition of tangency between the graph of f μ and the bisect and gives μ 0 c ¼ 0:6550ð8Þ and x c = 0.5216 (9). For μ > μ 0 c , the FPs x 1 and x 2 are found solving for f μ (x) = x, and, for instance, we find x 1 = 0.3326(5) and x 2 = 0.7890(9) for μ = 0.8.
The FP analysis also clarifies the general features of f μ that allow for the emergence of bistability, that is, in fact not contingent on the choice of f μ made in Eq. (5). Indeed, the only requirement is that, for some parameter(s) μ, the equation f μ (x) = x has three FPs x 0 < x 1 < x 2 , of which x 0 and x 2 are stable, whereas x 1 is unstable. Put simply, f μ should be a nonlinear function with a graph looking qualitatively as that of Fig. 4c. This condition guarantees a bistable phase for α < 1, which can then possibly extend to α ≥ 1 and, in the presence of a periodic drive, facilitate the establishment of a DTC.
Finally, note that higher resolution and smaller fluctuations could be achieved in the figures throughout the paper if simulating larger system sizes L and/or considering a larger number of independent runs R. This could, for instance, allow a more accurate characterization of both the equilibrium and the non-equilibrium phase diagrams of our model, which could be explored in other directions of the parameter space for varying α, μ, p d and T. This would, however, require a formidable numerical effort and goes therefore beyond the scope of this work. As a reference, for instance, the generation of Fig. 3e for the parameters considered therein requires a computing time of approximately 4 × 10 3 h per 3 GHz core.

Data availability
No data sets were generated or analysed during the current study.

Code availability
The codes that support the findings of this study are available at https://figshare.com/ articles/software/Code/13468836. Received: 2 July 2020; Accepted: 19 January 2021;  a For a control parameter μ<μ 0 c ¼ 0:6550ð8Þ, the system is inactive, corresponding to a single FP x 0 = 0: at long times, the system ends up in the empty, absorbing state with state variables s i = 0 for all sites i. b At the critical point μ ¼ μ 0 c , a new semi-stable FP emerges at x c = 0.5216 (9), that is, unstable from his left and stable on his right. c Increasing μ above μ 0 c , the semi-stable FP splits into an unstable FP x 1 < x c and a stable FP x 2 > x c . Depending on whether the initial density p 1 is <x 1 or >x 1 , the system flows towards density n = x 0 = 0 or n = x 2 > 0, respectively, indicating bistability.