Classical stochastic discrete time crystals

We describe a general and simple paradigm for discrete time crystals (DTCs), systems with a stable subharmonic response to an external driving ﬁeld, in a classical thermal setting. We consider, speciﬁcally, an Ising model in two dimensions, as a prototypical system with a phase transition into stable phases distinguished by a local order parameter, driven by thermal dynamics and periodically kicked with a noisy protocol. By means of extensive numerical simulations for large sizes—allowed by the classical nature of our model—we show that the system features a true disorder-DTC order phase transition as a function of the noise strength, with a robust DTC phase extending over a wide parameter range. We demonstrate that, when the dynamics is observed stroboscopically, the phase transition to the DTC state appears to be in the equilibrium two-dimensional Ising universality class. However, we explicitly show that the DTC is a genuine nonequilibrium state. More generally, we speculate that systems with thermal phase transitions to multiple competing phases can give rise to DTCs when appropriately driven.

The discovery of DTCs in unitary disordered Floquet systems raises numerous questions. Two important ones are the possibility of realizing DTCs in clean quantum systems [1,7,11,12,16,17,39,40] and whether DTCs can survive in the presence of dissipation [15,[41][42][43][44][45][46][47][48][49][50][51][52][53][54]. Studying these two issues is part of the more general search for an understanding of the range of mechanisms through which time crystalline order can be stabilized. In the case of quantum systems coupled to an external environment, the problem one faces is that of the natural tendency of dissipation to destabilize order [8,9,41,49]. In this respect, several mean-field or fully connected model systems have been shown to display DTC behavior with an appropriate engineering of the dissipative processes [43,44,46,50,[55][56][57], with some candidate open quantum systems argued to do the same away from mean field [45,48,49]. An important open question concerning the emergence of DTCs in open settings is their survival in the presence of thermal noise [19,20,49,58]. In this context, the emergence of DTCs in classical systems coupled to an external environment, inducing both damping and stochastic forces, have been recently investigated in Refs. [19,20]. Interestingly, in Ref. [19] an activated first-order dynamical DTC phase transition has been found in a one-dimensional driven Frenkel-Kontorova model at finite temperature while, in Ref. [20], the stabilization of truly many-body and robust period-doubled states have been demonstrated in a network of coupled nonlinear oscillators.
In this work, we shed light on the interplay between driving, noise, and interactions and the associated disorder-DTC order phase transition by considering a prototypical setting for DTCs in fully classical and thermal many-body systems. The generic scheme we propose is extremely simple, but to our knowledge has not been presented elsewhere before. The setup is that of a classical system with an equilibrium symmetry-breaking transition-below we consider a twodimensional (2D) Ising model as an obvious example-which is periodically driven. A period of the dynamics of duration τ consists of stochastic evolution under conditions such that asymptotically (i.e., if τ was to diverge) the symmetry-broken state would be stable, followed by a sudden random inversion of a fraction μ of the spins (see Fig. 1). We show that there exists a wide range of values of (τ, μ) in which a stable DTC emerges. This driven dynamics, when observed stroboscopically, leads to a nonequilibrium stationary state (NESS), in terms of which we construct the associated phase diagram as a function of (τ, μ). We then show that, despite the out-of-equilibrium nature of the dynamics, the transition from disorder to DTC order appears to be in the 2D Ising universality class.
Model and dynamical protocol. We study a classical Ising model on a 2D square lattice of linear size L and periodic boundary conditions, with nearest neighbor couplings and no magnetic field. The energy function is where σ i = ±1, sites are labeled by i, j = {1, . . . , L 2 }, i, j denotes nearest neighbors, σ = {σ i , . . . , σ N } with N = L 2 indicates a whole configuration of spins, and the coupling is ferromagnetic, i.e., J > 0. Below we set J = 1 for simplicity. The dynamics is periodic [59] according to the following two-step protocol within each period τ (see Fig. 1): (1) The first step evolves the system for time τ with a single spin-flip thermal dynamics which obeys detailed balance with respect to the equilibrium Boltzmann distribution associated to the energy function of Eq. (1) at temperature T . As we are interested in dynamics where a symmetry-broken state is stable, we will restrict to T = 0 below.
(2) The second step instantaneously inverts each spin with probability μ. When μ < 1 not all spins are inverted, and we denote with ε = 1 − μ the "error" introduced in the inversion step. This driven dynamics is described, within a period, by the evolution operator withσ x i and I i the x Pauli and the identity matrices acting on the ith spin, respectively. The rightmost factor implements the generator of zero-temperature Glauber dynamics [60][61][62].
where σ i indicates the configuration obtained from σ by flipping the ith spin, and θ (z) the Heaviside function. Therefore, only moves that do not increase the energy are allowed. The leftmost factor in Eq. (2) implements step (2) above. After a transient, the system reaches a time-periodic NESS. The dynamics depends on two parameters: the period τ and the rotated fraction of spins μ. To investigate the DTC phase transition it is convenient to consider the stroboscopic-time staggered magnetization, where m(t ) denotes the total magnetization per spin at time t, For a large enough number of cycles, we might expect that a stroboscopic NESS will be established. To characterize it, a natural order parameter is the average of the stroboscopic-time staggered magnetization of Eq. (4), Note that m st corresponds to the usual order parameter used to identify DTC order, given by the Fourier component of the magnetization at half the driving frequency [36][37][38]. DTC state and phase transition. We numerically investigate the dynamics given in Eq. (2) by using a continuous-time Monte Carlo algorithm [61,63,64] for the thermal part of the protocol. In Figs. 2(a) and 2(b) we show the (stroboscopic) behavior of the magnetization at times which are multiples of the period, m(nτ ), with n = 0, 1, 2, . . . , for τ = 1 and two different values of μ. Two phases can be clearly identified: Fig. 2(a) shows m(nτ ) for μ = 0.95 > μ c . Here, the system displays stable DTC oscillations with twice the period of the driving (ordered DTC phase). For the value of τ considered in Fig. 2, our estimate for the transition point between the ordered DTC and the disordered driven paramagnet is μ c ≈ 0.925 (see below for details). For a smaller value of μ = 0.85 < μ c (i.e., for a larger number of mistakes in the inversion) the oscillations in the magnetization decay quickly with time [see Fig. 2(b)], corresponding to the disordered (paramagnetic) phase (PM).
In terms of the staggered magnetization of Eq. (4), the disordered phase has m st = 0 while in the DTC phase m st = 0. This suggests that a phase transition takes place at the transition point μ c . To determine its nature and to precisely locate the transition point, we then inspect the behavior of Binder cumulants of the staggered magnetization [61,64,65]  with p 2. As shown in Fig. 2(c), the curves for Q 2p as a function of μ for different system sizes L and fixed τ cross at a single point. The size independence of Q 2p at this point implies the emergence of a diverging correlation length ξ (see below for details) associated with a continuous phase transition. Thus, the crossing point allows for a precise location of the critical point μ c for a given value of τ . The corresponding phase diagram of the disorder-DTC order phase transition in the (τ, μ) plane can be obtained by repeating the procedure described in Fig. 2(c) for different values of τ and it is shown in Fig. 2(d). Here, for a given τ , there is a μ c (τ ) such that μ < μ c (τ ) corresponds to the disordered phase, while μ > μ c (τ ) to the DTC phase. μ c (τ ) decreases with increasing τ as the longer annealing time allows for a larger error density before the DTC order becomes unstable. Figure 3 shows stroboscopically sampled configurations along representative trajectories in the various regimes of the dynamics. The top row corresponds to the DTC, with μ > μ c . Here, the magnetization at even and odd number of cycles alternates, displaying subharmonic behavior. The bottom row corresponds to the paramagnetic stroboscopic NESS at μ < μ c . Beyond fluctuations, there is no distinction between even and odd cycles, and discrete time symmetry remains unbroken. The middle row corresponds to conditions near criticality, μ ≈ μ c . Here, the presence of a continuous phase transition should result in scale invariance [66]. Indeed, this can be seen in the inset to the rightmost panel of the middle row, in which we consider two typical configurations of the system for two different sizes, L = 64 and L = 128, respectively. The latter spin configuration is re-scaled through a majority voting coarse graining on 2 × 2 plaquettes. Clearly almost perfect self-similarity, thus confirming scale invariance at criticality.
Critical properties of the DTC transition. Results from previous sections are indicative of a continuous phase transition with an emerging divergent correlation length at criticality. The latter describes the behavior of the (equal-time) connected correlation function, at large (spatial) distances. In this limit, we expect that C i, j ∼ e −r i j /ξ , with ξ being the correlation length and r i j ≡ |r i − r j | the distance between spins i and j. In a continuous phase transition a power-law behavior emerges in various thermodynamic quantities as the critical point is approached [66]. Examples for such quantities include ξ , the order parameter m st , and its susceptibility χ , defined as From the theory of continuous phase transitions, the following critical scaling relations are expected to hold in the stroboscopic NESS near the critical point with β, ν, and γ being the critical exponents corresponding to the various thermodynamic quantities. Despite that Eqs. (10) strictly hold in the thermodynamic limit, critical exponents can be efficiently extracted through a finite-size scaling analysis [61,64,65]. In a finite system the correlation length ξ is bounded by the system size L and all the various thermodynamic quantities (e.g., χ ) saturate when ξ ∼ L. Close to the critical point, by using that |μ − μ c | ∼ L −1/ν [see the  Fig. 2(c). The slope β of the fitting curve (black, dashed) has been fixed at the value obtained from the main panel. (c) Maximum value of the susceptibility, χ max ∼ L γ /ν , as a function of L. The fit leads to γ /ν = 1.69 ± 0.08. Using the result for ν from panel (a) we obtain γ = 1.8 ± 0.6. In panels (b) and (c), μ c is extracted as explained in Fig. 2. equation in the middle of Eqs. (10)], the following scaling relations can be obtained [64,65]: For a given L, μ max corresponds to the value of μ in which χ takes its maximum value, χ (μ max ) = χ max , and m st (μ c ) is the value of the staggered magnetization at the critical point. From Eqs. (11), by inspecting different system sizes, it is possible to extract the values of ν, β/ν, and γ /ν. The behavior of μ max , χ max , and m st (μ c ) as a function of L is shown in Fig. 4. By fitting with the power laws in Eqs. (11), we obtain the following values for the critical exponents of the disorder-DTC phase diagram: ν = 0.9 ± 0.3, β = 0.12 ± 0.04, and γ = 1.8 ± 0.6. These values are compatible with those of the equilibrium 2D Ising model [66], (ν 2D , β 2D , γ 2D ) = (1, 1/8, 7/4), even though in our case they correspond to a nonequilibrium phase transition, as shown in the next section. To benchmark the validity of the fitting procedure, in the inset of Fig. 4(b) we plot m st as a function of μ c = (μ − μ c )/μ c > 0 for different sizes and checked that, for large L, it follows the behavior predicted in Eqs. (10), i.e., m st ∼ |μ − μ c | β . Nonequilibrium nature of the DTC phase. To assess the nonequilibrium nature of the DTC phase transition, we now analyze the properties of the Floquet evolution operator F τ . At stationarity we have F τ |ρ NESS = |ρ NESS , with |ρ NESS the system's stationary state. Here, the NESS probability of finding the system in a given spin configuration σ is ρ σ = σ|ρ NESS , while the transition rate from the configuration σ to σ is T σ→σ = σ |F τ |σ . In these discrete time dynamics, one can define the average elementary currents between two configurations as J σ→σ ≡ ρ σ T σ→σ − ρ σ T σ →σ . In equilibrium settings, detailed balance requires ρ σ T σ→σ = ρ σ T σ →σ and, thus, on average no currents can be observed between any pairs of spin configurations σ and σ , i.e., J σ→σ = 0 ∀σ, σ . As such, in order to prove the nonequilibrium character of the phase transition that we are investigating, it would be sufficient to show the presence of at least one nonzero current [67]. For small sizes, the NESS |ρ NESS can be found through exact diagonalization of F τ , while the probabilities ρ σ and the transition rates T σ→σ can be directly evaluated by listing all the possible spin configurations. One can then construct the matrix current, whose entries are the elementary currents between all the possible spin configurations J σ→σ , and evaluate its L 1 norm (defined as σ,σ |J σ→σ |/2). Figure 5(a) gives the L 1 norm of the current matrix for a system with L = 3 and for a range of values of (τ, μ) and shows that, except for the limits τ = 0 and μ = 1, there exists at least one nonvanishing average current, resulting in a violation of detailed balance and, therefore, demonstrating that the disorder-DTC order is a nonequilibrium phase transition.
Such a nonequilibrium nature of the dynamics is also clearly apparent in continuous-time resolved dynamical realizations, as displayed in Fig. 5(b) for the energy and the magnetization, which are manifestly not time reversible.
Having shown that the asymptotic state is a nonequilibrium one, we can ask whether the analogy with the 2D Ising model extends beyond the stroboscopic NESS. A simple test is to consider the dynamics starting from an initial state with zero magnetization (such as a quench from a random configuration). Figure 5(c) shows that under DTC conditions the stroboscopic dynamics exhibits coarsening, displaying progressive growth of domains of both magnetized states (before eventual collapse to one of the two in the last panel, which is a finite-size effect). As the relaxation step of the protocol is standard Glauber dynamics at zero temperature we could have expected coarsening only within a period, but the fact that it survives for long times despite the periodic driving is nontrivial, and constitutes another indication that the DTC state is robust.
Conclusions. We have shown that a DTC phase can be obtained in a fully classical and thermal setting-and in the absence of disorder or any form of classical localization-as a symmetry-breaking transition of a driven 2D Ising model. The classical nature of the problem we consider allowed us, by means of numerical simulations of large system sizes and finite-size scaling (something often outside the range of possibility in both closed and open quantum systems), to demonstrate convincingly that the transition to a DTC is indeed a phase transition. Our results suggest that when observed stroboscopically, this phase transition in the NESS is in the 2D Ising universality class. The mechanism for DTCs described here is simple and easily generalizable. For example, we expect that driven Potts models will lead to DTCs with periods larger than two [68,69]. The simplicity of the scheme also suggests that it should be possible to observe these classical DTCs in experiments such as magnetic solid state systems and other platforms with order-disorder transitions with a scalar order parameter. In particular, the dynamics we have studied in this work can be implemented in a setup of dissipative Rydberg atoms [70][71][72] in the limit of strong dephasing [73][74][75] and driven by a periodic sequence of slightly imperfect π -pulse rotations, which simulates the random inversion of a fraction of spins.