Classical Stochastic Discrete Time Crystals

TL;DR

This paper introduces a classical thermal paradigm for discrete time crystals using a 2D Ising model, demonstrating stable subharmonic dynamics and phase transitions akin to equilibrium Ising behavior, but as a genuine non-equilibrium state.

Contribution

It proposes a simple classical thermal framework for DTCs in a 2D Ising model, showing stable subharmonic phases and their relation to equilibrium phase transitions.

Findings

Stable DTC phases emerge over a wide parameter range.

DTC phase transition resembles 2D Ising class in stroboscopic observation.

DTCs are confirmed as genuine non-equilibrium states.

Abstract

We describe a possible general and simple paradigm in a classical thermal setting for discrete time crystals (DTCs), systems with stable dynamics which is subharmonic to the driving frequency thus breaking discrete time-translational invariance. We consider specifically 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 a thermal dynamics and periodically kicked. We show that for a wide parameter range a stable DTC emerges. The phase transition to the DTC state appears to be in the equilibrium 2D Ising class when dynamics is observed stroboscopically. However, we show that the DTC is a genuine non-equilibrium state. More generally, we speculate that systems with thermal phase transitions to multiple competing phases can give rise to DTCs when appropriately driven.

Figures

Equations

$$E(\mathbf{\sigma })\equiv -J\sum _{⟨i,j⟩}{\sigma }_i{\sigma }_j,$$

$${\mathcal{F}}_{\tau }\equiv \left[\prod _{i=1}^{\mathcal{N}}\left(\epsilon {\mathbb{I}}_i+\mu {\hat{\sigma }}_i^x\right)\right]e^{\mathcal{L}\tau },$$

$$\mathcal{L}\equiv \sum _{i=1}^{\mathcal{N}}({\hat{\sigma }}_i^x-{\mathbb{I}}_i){\Gamma }_i(\mathbf{\sigma })$$

$$m_{\mathrm{st}}(n)\equiv (-1{)}^{n}m(n\tau ),\text{with }n=0,1,\dots$$

$$m(t)\equiv \frac{1}{\mathcal{N}}\sum _{i=1}^{\mathcal{N}}{\sigma }_i(t).$$

$$⟨m_{\mathrm{st}}⟩=\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=1}^N{m}_{\mathrm{st}}(n).$$

$$Q_{2p}=\frac{⟨m_{\mathrm{st}}^{2p}⟩}{{⟨|m_{\mathrm{st}}{|}^p⟩}^2},$$

$$C_{i,j}(t)\equiv ⟨{\sigma }_i(t){\sigma }_j(t)⟩-⟨{\sigma }_i(t)⟩⟨{\sigma }_j(t)⟩,$$

$$\chi \equiv L^2\left(⟨m_{\mathrm{st}}^2⟩-{⟨m_{\mathrm{st}}⟩}^2\right).$$

$$⟨m_{\mathrm{st}}⟩\sim |\mu -{\mu }_c{|}^{\beta },\xi \sim |\mu -{\mu }_c{|}^{-\nu },\chi \sim |\mu -{\mu }_c{|}^{-\gamma },$$

$${\mu }_{\max }\sim L^{-1/\nu },{\chi }_{\max }\sim L^{\gamma /\nu },⟨m_{\mathrm{st}}({\mu }_c)⟩\sim L^{-\beta /\nu }.$$

Full text

Introduction.– Discrete time crystals (DTCs) Sacha (2015); Khemani et al. (2016); Else et al. (2016); von Keyserlingk and Sondhi (2016); Khemani et al. (2017); Yao et al. (2017); Huang et al. (2018); Choi et al. (2017); Zhang et al. (2017); Ho et al. (2017); Else et al. (2016, 2017); Russomanno et al. (2017); Rovny et al. (2018a, b); Pal et al. (2018); Yu et al. (2019); Barfknecht et al. (2019); Huang et al. (2018); Flicker (2018); Yao et al. (2018); Heugel et al. (2019) have recently emerged as a novel form of non-equilibrium quantum matter. A time crystal is a system in which the time translation invariance of the dynamics is spontaneously broken asymptotically Wilczek (2012); Shapere and Wilczek (2012). While time crystallinity is impossible in thermal equilibrium states of systems with a local time-independent Hamiltonian Bruno (2013); Nozières (2013); Watanabe and Oshikawa (2015); Kozin and Kyriienko (2019), it can emerge in excited states Syrwid et al. (2017) or through spontaneous breaking of an underlying discrete time symmetry Sacha (2015); Else et al. (2016); Khemani et al. (2016). Concrete examples of the latter mechanism have been found in closed quantum systems with time-periodic Hamiltonians (“Floquet” systems Shirley (1965); Sambe (1973); Grifoni and Hänggi (1998)) such as the $\pi$-spin glass and related models Else et al. (2016); Khemani et al. (2016); von Keyserlingk and Sondhi (2016); Khemani et al. (2017); Yao et al. (2017); Choi et al. (2017); Zhang et al. (2017); Ho et al. (2017). The usual setting in such unitary spin systems is that the dynamics is split between evolution with an interacting disordered Hamiltonian followed by a rotation of the spins. The usual Floquet heating towards an infinite temperature state D’Alessio and Rigol (2014); Lazarides et al. (2014); Ponte et al. (2015a) is avoided by exploiting localization Lazarides et al. (2015); Ponte et al. (2015b) so that spatio-temporal order can be established without fine-tuning. For reviews see, e.g., Refs. Moessner and Sondhi (2017); Sacha and Zakrzewski (2018); Else et al. (2019).

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 Sacha (2015); Else et al. (2017); Russomanno et al. (2017); Yu et al. (2019); Huang et al. (2018); Barfknecht et al. (2019); Giergiel et al. (2018); Schäfer et al. (2019) and whether DTCs can survive in the presence of dissipation Lazarides and Moessner (2017); Iemini et al. (2018); Wang et al. (2018); Gong et al. (2018); Gambetta et al. (2019); Pal et al. (2018); O’Sullivan et al. (2018); Tucker et al. (2018); Zhu et al. (2019); Lazarides et al. (2019); Droenner et al. (2019); Buča et al. (2019); Calvanese Strinati et al. (2019); Carlo et al. (2019); Bello et al. (2019). 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 Lazarides and Moessner (2017); Zhang et al. (2017); Choi et al. (2017); Lazarides et al. (2019). 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 Wang et al. (2018); Gong et al. (2018); O’Sullivan et al. (2018); Droenner et al. (2019); Cosme et al. (2019); Keßler et al. (2019); Lledó et al. (2019), with some candidate open quantum systems argued to do the same away from mean-field Gambetta et al. (2019); Lazarides et al. (2019); Zhu et al. (2019). An important open question concerning the emergence of DTCs in open settings is their survival in the presence of thermal noise Yao et al. (2018); Heugel et al. (2019); Oberreiter et al. (2019); Lazarides et al. (2019). 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. Yao et al. (2018); Heugel et al. (2019). Interestingly, in Ref. Yao et al. (2018) 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. Heugel et al. (2019), the stabilization of truly many-body and robust period-doubled states have been demonstrated in a network of coupled non-linear 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 two-dimensional (2D) Ising model as an obvious example - which is periodically driven. A period of the dynamics of duration $\tau$ consists of stochastic evolution under conditions such that asymptotically (i.e., if $\tau$ was to diverge) the symmetry-broken state would be stable, followed by a sudden random inversion of a fraction $\mu$ of the spins, see Fig. 1. We show that there exists a wide range of values of $(\tau,\mu)$ in which a stable DTC emerges. This driven dynamics, when observed stroboscopically, leads to a non-equilibrium stationary state (NESS), in terms of which we construct the associated phase diagram as function of $(\tau,\mu)$. 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 neighbour couplings and no magnetic field. The energy function is

[TABLE]

where $\sigma_{i}=\pm 1$, sites are labelled by $i,j=\{1,...,L^{2}\}$, $\left<i,j\right>$ denotes nearest neighbours, ${\bm{\sigma}}=\{\sigma_{i},\ldots,\sigma_{\mathcal{N}}\}$ with $\mathcal{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 Note (1) according to the following two-step protocol within each period $\tau$, see Fig. 1:

The first step evolves the system for time $\tau$ 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. 2.

The second step instantaneously inverts each spin with probability $\mu$. When $\mu<1$ not all spins are inverted, and we denote with $\varepsilon=1-\mu$ the “error” introduced in the inversion step.

This driven dynamics is described, within a period, by the evolution operator

[TABLE]

with $\hat{\sigma}^{x}_{i}$ and $\mathbb{I}_{i}$ the $x$ Pauli and the identity matrices acting on the $i$-th spin, respectively. The rightmost factor implements step 1 above, with

[TABLE]

the generator of zero-temperature Glauber dynamics Glauber (1963); Binder and Heermann (2010); Krapivsky et al. (2010). Here, $\Gamma_{i}({\bm{\sigma}})=\theta[E({\bm{\sigma}})-E({\bm{\sigma}}_{i})]$, where ${\bm{\sigma}}_{i}$ indicates the configuration obtained from ${\bm{\sigma}}$ by flipping the $i$-th spin, and $\theta(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 $\tau$ and the rotated fraction of spins $\mu$. To investigate the DTC phase transition it is convenient to consider the stroboscopic-time staggered magnetization,

[TABLE]

where $m(t)$ denotes the total magnetization per spin at time $t$,

[TABLE]

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),

[TABLE]

Note that $\left<m_{\rm st}\right>$ corresponds to the usual order parameter used to identify DTC order, given by the Fourier component of the magnetization at half the driving frequency Moessner and Sondhi (2017); Sacha and Zakrzewski (2018); Else et al. (2019).

DTC state and phase transition.– We numerically investigate the dynamics given in Eq. (2) by using a continuous-time Monte Carlo algorithm Bortz et al. (1975); Newman and Barkema (1999); Binder and Heermann (2010) for the thermal part of the protocol. In Figs. 2(a,b) we show the (stroboscopic) behavior of the magnetization at times which are multiples of the period, $m(n\tau)$, with $n=0,1,2,\ldots$, for $\tau=1$ and two different values of $\mu$. Two phases can be clearly identified: Fig. 2(a) shows $m(n\tau)$ for $\mu=0.95>\mu_{c}$. Here, the system displays stable DTC oscillations with twice the period of the driving (ordered DTC phase). For the value of $\tau$ considered in Fig. 2, our estimate for the transition point between the ordered DTC and the disordered driven paramagnet is $\mu_{c}\approx 0.925$ (see below for details). For a smaller value of $\mu=0.85<\mu_{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 $\left<m_{\rm st}\right>=0$ while in the DTC phase $\left<m_{\rm st}\right>\neq 0$. This suggests that a phase transition takes place at the transition point $\mu_{c}$. To determine its nature and to precisely locate the transition point, we then inspect the behavior of Binder cumulants of the staggered magnetization Newman and Barkema (1999); Binder and Heermann (2010); Sandvik (2010)

[TABLE]

with $p\geq 2$. As shown in Fig. 2(c), the curves for $Q_{2p}$ as a function of $\mu$ for different system sizes $L$ and fixed $\tau$ cross at a single point. The size-independence of $Q_{2p}$ at this point implies the emergence of a diverging correlation length $\xi$ (see below for details) associated with a continuous phase transition. Thus, the crossing point allows for a precise location of the critical point $\mu_{c}$ for a given value of $\tau$. The corresponding phase diagram of the disorder-DTC order phase transition in the $(\tau,\mu)$ plane can be obtained by repeating the procedure described in Fig. 2(c) for different values of $\tau$ and it is shown in Fig. 2(d). Here, for a given $\tau$, there is a $\mu_{c}(\tau)$ such that $\mu<\mu_{c}(\tau)$ corresponds to the disordered phase, while $\mu>\mu_{c}(\tau)$ to the DTC phase. $\mu_{c}(\tau)$ decreases with increasing $\tau$ 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 $\mu>\mu_{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 $\mu<\mu_{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, $\mu\approx\mu_{c}$. Here, the presence of a continuous phase transition should result in scale invariance Chaikin and Lubensky (1995). 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\times 2$ plaquettes. Clearly, the $L=64$ and the re-scaled $L=128$ configurations show an 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,

[TABLE]

at large (spatial) distances. In this limit, we expect that $C_{i,j}\sim e^{-r_{ij}/\xi}$, with $\xi$ being the correlation length and $r_{ij}\equiv|\bm{r}_{i}-\bm{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 Chaikin and Lubensky (1995). Examples for such quantities include $\xi$, the order parameter $\left<m_{\mathrm{st}}\right>$ and its susceptibility $\chi$, defined as

[TABLE]

From the theory of continuous phase transitions, the following critical scaling relations are expected to hold in the stroboscopic NESS near the critical point

[TABLE]

with $\beta$, $\nu$, and $\gamma$ being the critical exponents corresponding to the various thermodynamic quantities. Despite Eq. (10) strictly holds in the thermodynamic limit, critical exponents can be efficiently extracted through a finite-size scaling analysis Newman and Barkema (1999); Binder and Heermann (2010); Sandvik (2010). In a finite system the correlation length $\xi$ is bounded by the system size $L$ and all the various thermodynamic quantities (e.g., $\chi$) saturate when $\xi\sim L$. Close to the critical point, by using that $|\mu-\mu_{c}|\sim L^{-1/\nu}$ [see the equation in the middle of Eq. (10)], the following scaling relations can be obtained Newman and Barkema (1999); Sandvik (2010)

[TABLE]

For a given $L$, $\mu_{\mathrm{max}}$ corresponds to the value of $\mu$ in which $\chi$ takes its maximum value, $\chi(\mu_{\mathrm{max}})=\chi_{\mathrm{max}}$, and $\left<m_{\mathrm{st}}(\mu_{c})\right>$ is the value of the staggered magnetization at the critical point. From Eq. (11), by inspecting different system sizes, it is possible to extract the values of $\nu$, $\beta/\nu$, and $\gamma/\nu$. The behavior of $\mu_{\mathrm{max}}$, $\chi_{\mathrm{max}}$, and $\left<m_{\mathrm{st}}(\mu_{c})\right>$ as a function of $L$ is shown in Fig. 4. By fitting with the power laws in Eq. (11), we obtain the following values for the critical exponents of the disorder-DTC phase diagram: $\nu=0.9\pm 0.3$, $\beta=0.12\pm 0.04$ and $\gamma=1.8\pm 0.6$. These values are compatible with the ones of the equilibrium 2D Ising model Chaikin and Lubensky (1995), $(\nu_{2D},\beta_{2D},\gamma_{2D})=(1,1/8,7/4)$, even though in our case they correspond to a non-equilibrium phase transition, as shown in next section. To benchmark the validity of the fitting procedure, in the Inset of Fig. 4(b) we plot $\left<m_{\mathrm{st}}\right>$ as a function of $\mu_{c}=(\mu-\mu_{c})/\mu_{c}>0$ for different sizes and checked that, for large $L$, it follows the behavior predicted in Eq. (10), i.e., $\left<m_{\mathrm{st}}\right>\sim|\mu-\mu_{c}|^{\beta}$.

Non-equilibrium nature of the DTC phase.— To assess the non-equilibrium nature of the DTC phase transition, we now analyze the properties of the Floquet evolution operator $\mathcal{F}_{\tau}$. At stationarity we have $\mathcal{F}_{\tau}\left|\rho_{\mathrm{NESS}}\right>=\left|\rho_{\mathrm{NESS}}\right>$, with $\left|\rho_{\mathrm{NESS}}\right>$ the system’s stationary state. Here, the NESS probability of finding the system in a given spin configuration $\bm{\sigma}$ is $\rho_{\bm{\sigma}}=\left<\bm{\sigma}|\rho_{\mathrm{NESS}}\right>$, while the transition rate from the configuration $\bm{\sigma}$ to $\bm{\sigma^{\prime}}$ is $\mathcal{T}_{\bm{\sigma}\rightarrow\bm{\sigma}^{\prime}}=\left<\bm{\sigma}^{\prime}\right|\mathcal{F}_{\tau}\left|\bm{\sigma}\right>$. In these discrete time dynamics, one can define the average elementary currents between two configurations as $J_{{\bm{\sigma}}\to{\bm{\sigma}}^{\prime}}\equiv\rho_{{\bm{\sigma}}}\mathcal{T}_{{\bm{\sigma}}\to{\bm{\sigma}}^{\prime}}-\rho_{{\bm{\sigma}}^{\prime}}\mathcal{T}_{{\bm{\sigma}}^{\prime}\to{\bm{\sigma}}}$. In equilibrium settings, detailed balance requires $\rho_{{\bm{\sigma}}}\mathcal{T}_{{\bm{\sigma}}\to{\bm{\sigma}}^{\prime}}=\rho_{{\bm{\sigma}}^{\prime}}\mathcal{T}_{{\bm{\sigma}}^{\prime}\to{\bm{\sigma}}}$ and, thus, on average no currents can be observed between any pairs of spin configurations ${\bm{\sigma}}$ and ${\bm{\sigma}}^{\prime}$, i.e., $J_{{\bm{\sigma}}\to{\bm{\sigma}}^{\prime}}=0\ \forall\bm{\sigma},\bm{\sigma^{\prime}}$. As such, in order to prove the non-equilibrium character of the phase transition that we are investigating, it would be sufficient to show the presence of at least one non-zero current Seifert (2012). For small sizes, the NESS $\left|\rho_{\mathrm{NESS}}\right>$ can be found through exact diagonalization of $\mathcal{F}_{\tau}$, while he probabilities $\rho_{\bm{\sigma}}$ and the transition rates $\mathcal{T}_{{\bm{\sigma}}\to{\bm{\sigma}}^{\prime}}$ 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_{{\bm{\sigma}}\to{\bm{\sigma}}^{\prime}}$, and evaluate its $L_{1}$-norm (defined as $\sum_{{\bm{\sigma}},{\bm{\sigma}}^{\prime}}|J_{{\bm{\sigma}}\to{\bm{\sigma}}^{\prime}}|/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 $(\tau,\mu)$ and shows that, except for the limits $\tau=0$ and $\mu=1$, there exists at least one non-vanishing average current, resulting in a violation of detailed balance and, therefore, demonstrating that the disorder-DTC order is a non-equilibrium phase transition.

Such a non-equilibrium 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 non-equilibrium 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 non-trivial, 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 Surace et al. (2019); Pizzi et al. (2019). 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 Löw et al. (2012); Lee et al. (2012); Marcuzzi et al. (2014a) in the limit of strong dephasing Lesanovsky and Garrahan (2013); Marcuzzi et al. (2014b); Levi et al. (2016) and driven by a periodic sequence of slightly imperfect $\pi$-pulse rotations, which simulates the random inversion of a fraction of spins.

11footnotetext: Ising models driven via a periodically oscillating magnetic field has been studied in Refs. Chakrabarti and Acharyya (1999); Korniss et al. (2000); Park and Pleimling (2013). However, they do not lead to DTC behavior. Note that the inversion step of Eq. (2) in our protocol cannot be realized by a time-dependent change of the energy function (such as a with an oscillatory magnetic field).