Thermal melting of discrete time crystals: a dynamical phase transition induced by thermal fluctuations

Introduction – Spontaneous symmetry breakings (SSB) and universality classes are among the most fundamental concepts in modern physics: the former is used to characterize different phases of matter while the latter is to categorize the transitions between them. Compared to equilibrium systems, the physics in out-of-equilibrium systems is much richer but less known in general. The evolutionary processes give rise to increasing richness of the paradigm of spontaneous symmetry breaking and universality classes, which could take place not only in space but also in time, the latter opens up new opportunities to explore novel states of matter (the time crystal(TC)Wilczek (2012) for instance) and dynamical universality classes ( the critical dynamical phenomenaHohenberg and Halperin (1977); Taeuber (2014) or the Kardar-Parisi-Zhang (KPZ) universality classKardar et al. (1986) for instance) beyond the scope of equilibrium physics.

As a prototypical example of spontaneous symmetry breaking in time domain, the time crystal phaseWilczek (2012), in its different formsShapere and Wilczek (2012); Li et al. (2012); Wilczek (2013); Sacha (2015); Else et al. (2016); Khemani et al. (2016); Yao et al. (2017); Russomanno et al. (2017); Gong et al. (2018); Huang et al. (2018); Iemini et al. (2018); Das et al. (2018); Zhu et al. (2019); Liao et al. (2019); Kozin and Kyriienko (2019); Cai et al. (2020); Chinzei and Ikeda (2020); Lyu et al. (2020); Choudhury and Liu (2021), has attracted considerable interests in past decades. In spite of being proven absent in thermal equilibriumBruno (2013); Watanabe and Oshikawa (2015), such an intriguing phase has been observed in periodically driven non-equilibrium settingsChoi et al. (2017); Zhang et al. (2017); Autti et al. (2018); Träger et al. (2021); Keßler et al. (2021); Mi et al. (2021); Frey and Rachel (2021). These systems exhibit oscillations with period doubling with respect to that of the external driving, thus spontaneously break the discrete time translational symmetry(TTS) from symmetry group $\mathbb{Z}$ to $2\mathbb{Z}$. As a phase of matter, a profound question is its stability against perturbations, especially the thermal fluctuations inevitable in almost all the realistic experimental setups. From a theoretical point of view, this problem can be considered as a non-equilibrium analogue of thermal melting of the crystalline order, and is of fundamental interest due to its relevance to broader questions of the robustness of the temporal order against fluctuations, as well as the dynamical universality classes in non-equilibrium matter.

Recently, Yao et al explored this question by studying a classical Hamiltonian dynamics coupled to a finite-temperature bath, and found an activated discrete time crystal(DTC) whose crystalline order survives to long, but not infinite times at low temperature Yao et al. (2020). A question naturally raised is whether there exists DTC phases with true long-range crystalline order persisting forever Zhuang et al. (2021)? Is it possible for the thermal fluctuation to melt such a dynamical crystalline order, just as it does for conventional crystals? If so, how to characterize such a non-equilibrium phase transition and what’s the corresponding universality class?

In this paper, we attempt to answer these questions, focusing for simplicity on a three-dimensional (3D) classical periodically-driven interacting spin model, which could exhibit period doubling with respect to the driving (DTC phase) at zero temperature. The thermal fluctuations are introduced by coupling each spin to a heat bath, which provides both dissipation and noise. The problem is approached by solving a stochastic equation of motion (EOM), where the dynamics of each spin is governed by a Landau-Lifshitz-Gilbert (LLG) equationLakshmanan (2010) augmented by a random thermal forceBrown (1963). It has been shown that at high temperature, the DTC order parameter decays exponentially with time, indicating a finite life time of the DTC phase as in the 1D case (activated DTC). With decreasing temperature, the non-equilibrium system experiences a dynamical phase transition from an activated DTC to a DTC phase with true long-range crystalline order. It is interesting to show that despite the genuine non-equilibrium feature of such periodically driven systems, its critical properties are characterized by an Ising universality class with a dynamical exponent ($z=2$) identical with the value in the critical dynamics of undriven systems.

Model and method – The time crystal phases in classical many-body systems have attracted considerable interest recentlyLupo and Weber (2019); Yao et al. (2020); Gambetta et al. (2019); Khasseh et al. (2019); Hurtado-Gutiérrez et al. (2020); Ye et al. (2021); Pizzi et al. (2021). We start with a 3D classical spin model (transverse Ising model) with a time-dependent coupling strength. The system Hamiltonian reads:

where the dynamical variable $\mathbf{s}_{i}$ is a classical vector with fixed length $|\mathbf{s}_{i}|=1$. The summation $\langle ij\rangle$ is over the nearest-neighboring (NN) sites of site $i$ in the cubic lattice with length $L$. The strength of the interaction $V(t)$ periodically oscillates as $V(t)=V_{0}\cos\omega_{0}t$. $h_{z}$ is the strength of a uniform magnetic field along z-direction. Throughout this paper we fix $h_{z}=J$. In the absence of thermal bath, the dynamics of each spin can be described by the EOM: $\dot{\mathbf{s}}_{i}=\mathbf{h}^{0}_{i}\times\mathbf{s}_{i}$, where the effective magnetic field $\mathbf{h}^{0}_{i}=[V(t)\bar{s}_{i}^{x},0,h_{z}]$ with $\bar{s}_{i}^{x}=\frac{1}{2}\sum_{\langle j\rangle}s_{j}^{x}$ (the summation is over the NN sites of site i).

Even though our driven spin system is a genuine non-equilibrium system whose temperature is ill-defined, usually the degrees freedom of bath are much larger than those of the system thus the back action of the system to the bath can be neglected, which allows us to consider a thermal bath with a well-defined temperature (T) that doesn’t change during the dynamical process. The effect of the thermal bath can be modeled by methods familiar in the theory of the Brownian motion and other stochastic processes, where the EOM of each spin is described by a stochastic LLG equationLakshmanan (2010) as:

where $\lambda$ is the strength of the friction provided by the local thermal bath, which is fixed as $\lambda=J$ throughout this paper. $\mathbf{h}_{i}=\mathbf{h}^{0}_{i}+\bm{\xi}_{i}(t)$ is the effective magnetic field, where $\bm{\xi}_{i}(t)$ is a three-dimensional(3D) stochastic magnetic field representing the thermal noise, which can be approximated as white noise if assuming that the typical time scale of the bath is much shorter than that of the system. We further assume the local bath around each spin is independent with each other thus the stochastic variables satisfy:

where $\alpha,\beta=x,y,z$, $\mathcal{D}$ is the strength of the noise, and the average $\langle\rangle_{\bm{\xi}}$ is over all the noise trajectories. If the bath is in thermal equilibrium with temperature T, the fluctuation-dissipation theorem (FDT) indicates that the strengths of the friction and noise satisfy

Numerically, we adopt Stratonovich’s formula of the stochastic differential Eq.(7), and solve it using the standard Heun methodAment et al. (2016) with the time step of $\Delta t=10^{-3}$, the convergence of which has been checked numerically (see the Supplmentary Material(SM)Sup ). The system size in our simulation ranges from $L=8$ to $L=28$, which enable us to systematically analyze the finite-size effect. The ensemble average over the noise trajectories can be performed by directly sampling over $\mathcal{N}$ sets of noise realizations. ($\mathcal{N}$ ranges from $4\times 10^{3}$ to $10^{5}$ depending on the simulated system size). In our simulation, we calculate the evolution of the average magnetization along x-direction $M(t)=\langle\frac{1}{L^{3}}\sum_{i}s_{i}^{x}(t)\rangle_{\bm{\xi}}$ to characterize various dynamical behaviors. We choose the initial state as the ground state of Hamiltonian.(1) with $t=0$, since we start from a spatially uniform initial state, $M(t)$ is proportional to the auto-correlation function $C(t)=\langle\frac{1}{L^{2}}\sum_{i}s_{i}^{x}(0)s_{i}^{x}(t)\rangle_{\bm{\xi}}$, which characterizes the memory effect of the initial state information. However, we find that the long-time behavior doesn’t change even if we start from a non-uniform random initial state, that different spins will finally synchronize with each otherSup .

Zero temperature phase diagram: – It is shown that even at zero temperature ($\mathcal{D}=0$), the classical EOM.(7) could exhibit rich long-time dynamical behaviors, which can be found in the $V_{0}-\omega_{0}$ phase diagram as shown in Fig.1. For sufficiently large $V_{0}$ and small $\omega_{0}$, one can find a $DTC$-2 phase, where the ferromagnetic order parameter $M(t)$ oscillate with a period twice of that of the driving ($\frac{2\pi}{\omega_{0}}$). Besides that, one can also find other DTC phases whose periods are other integer multiples of $\frac{2\pi}{\omega_{0}}$ ($DTC$-3 and $DTC$-4 phases for instance). If both $V_{0}$ and $\omega_{0}$ are large enough, one can find a synchronization between $M(t)$ ad $V(t)$, both of which oscillates in the same period(synchronized phase). Besides these periodic oscillations, one can also find other dynamical phases that $M(t)$ could decay to zero (paramagnetic phase), or oscillate chaotically (chaos).

The fate of DTC phase at various temperatures: – In the following, we will focus on the $DTC$-2 phase by fixing $V_{0}=5J$, $\omega_{0}=2\pi J$ and change the temperature. The noises, in general, will suppress the long-range order. As shown in Fig.2 (a), $M(t)$ exhibits a damped oscillation whose amplitude decays exponentially with time ($\sim e^{-\Delta t}$ as shown in the inset) at a relatively high temperature (e.g $\mathcal{D}=0.29J$ corresponding to $T=0.042J$) . Such an exponential decay can be understood as a consequence of the “kinks” activated by thermal fluctuations: the noise can induce tunnelings from one of the degenerate DTC phases to the other by a $\pi$-phase shiftingYang and Cai (2021), as shown in Fig.2 (b) for the time evolution of M(t) under a single noise trajectory. Such topological “excitations”, no matter how rare they are, will inevitably destroy the long-range crystalline order along the temporal direction via a similar mechanism responsible for the absence of spatial long-range order in 1D thermal equilibrium systems. At low temperature (e.g $\mathcal{D}=0.26J$ corresponding to $T=0.038J$), the thermal fluctuation is not strong enough to induce a $\pi$-phase shifting, but can only modulate the amplitude of the oscillation, thus we find a persistent oscillation of $M(t)$ as shown in Fig.2 (a), which indicates a DTC phase with true long-range order in time domain.

The critical behavior of the non-equilibrium phase transition– As the temperature is lowered, the exponent of the exponential decay $\Delta$ decreases and finally vanishes when the temperature drops below a critical temperature $T_{c}=\frac{\mathcal{D}_{c}^{2}}{2J}$ with $D_{c}=0.2782J$, which indicates an absence of the exponential decay in the thermodynamic limit. For a sufficiently large system, it is shown in Fig.2 (c) that at the critical temperature ($T=T_{c}$), the amplitude of $M(t)$ decays algebraically with time ($\sim t^{-\alpha_{c}}$ with $\alpha_{c}=0.33(1)$), qualitatively different from either the exponential decay at the high temperature or the persistent oscillation at low temperature. Such a critical exponent $\alpha_{c}=0.33$ is different from the value observed recently ($\alpha_{c}=0.5$) in a mean-field analysis of the critical properties of a dissipative DTC phaseChinzei and Ikeda (2021).

For a finite system, the dynamics of $M(t)$ at the critical temperature can be separated into three regimes: a transient regime that depends on the initial state, a long-time exponential regime dominated by the finite size effect, and an intermediate universal algebraic regime sandwiched between them, as shown in Fig.2 (d). The finite size effect leads to an exponential decay after sufficiently long time. However, as the system size increases, a finite size scaling of $\Delta_{c}(L)$ (the inset of Fig.2 (d)) suggests that $\Delta_{c}(L)$ vanishes in the thermodynamic limit as $\Delta_{c}(L)\sim L^{-z}$ with $z=2$. This result indicates a divergence of the relaxation time with the system size as $\tau_{c}\sim L^{z}$, where z is the dynamical critical exponent characterizing the dynamical universality classHohenberg and Halperin (1977). It is known that in a relaxation dynamics of a critical kinetic Ising model without periodical driving, $z$ could be either $2$Glauber (1963) or $3$Huse (1986) depending on the conservation law during the relaxation dynamics. Here, we find the dynamical critical exponent $z$ of our model agrees well with the value in relaxation dynamics of undriven kinetic 2D Ising model without conservation law (Glauber modelGlauber (1963)). In the thermodynamic limit, $\Delta_{c}(L)$ vanishes and the exponential regime gives ways to the algebraic regime, which persists forever.

To study other critical exponents of this dynamical phase transition, we focus on the divergence of the relaxation time $\tau$ (the inverse of $\Delta$) close to the critical point. We first define a renormalized relaxation time $\tilde{\tau}=\tau/L^{z}$ with $z=2$, and plot its dependence on $\mathcal{D}$ for systems with various system sizes $L$ around the critical point. As shown in Fig.2 (e), one can find a scaling invariant point at $\mathcal{D}=\mathcal{D}_{c}$, a signature of a continuous phase transition. Near the critical point, the relaxation time diverges as $\tau\sim|T-T_{c}|^{-z\nu}$. To obtain the critical exponent $\nu$, we perform the data collapses as shown in the insets of Fig.2 (e), from which we can extract that the critical exponents $\nu=0.64$, which agrees with the value ($\nu=0.642$) in the 3D Ising universality class .

Another important quantity to characterize the critical behavior is the correlation length $l$, which can be extracted from the equal-time correlation function $S(\mathbf{r},t)=\langle s_{1}^{x}(t)s_{\mathbf{r}}^{x}(t)\rangle_{\bm{\xi}}-\langle s_{1}^{x}(t)\rangle_{\bm{\xi}}\langle s_{\mathbf{r}}^{x}(t)\rangle_{\bm{\xi}}$ asSandvik (2010):

where $S(\mathbf{Q},t)=\frac{1}{L^{3}}\sum_{\mathbf{r}}e^{i\mathbf{Q}\cdot\mathbf{r}}S(\mathbf{r},t)$ is the structure factor of the spin-spin correlation. $\mathbf{q}_{x}=[\frac{2\pi}{L},0,0]$ and $l_{x}(t)$ is the correlation length along x-direction at time t. As shown in Fig.2 (f), at the critical point, the correlation length diverges algebraically in time as $l(t)\sim t^{\beta}$ with $\beta=0.55(3)$, which agrees with the dynamical critical exponent $\beta=1/z$ in the Glauber model.

Discussion – Despite the presence of the periodical driving of our model, we find that its dynamical critical behaviors are similar with the relaxation dynamics of the undriven kinetic Ising model near the critical point (the model A in Hohenberg-Halperin’s alphabetic classification of dynamical critical phenomenaHohenberg and Halperin (1977)). This result can be understood by a time-independent description of stroboscopic dynamics under periodic driving. Such a time-independent Hamiltonian is similar with the Floquet Hamiltonian in the periodically-driven quantum systems, and can be derived using an approximation scheme similar with the Magnus expansion, a series expansion in terms of the driving periodMagnus (1954). Even though higher order terms (the four-site or multi-site interactions for instance) may appear in the effective time-independent Hamiltonian, they could be irrelevant for the long-time dynamics near the critical point, thus the stroboscopic dynamics of our periodically driven system is similar with the relaxation dynamics of a time-independent kinetic Ising model near the critical point.

This similarity reminds us of the dynamical phase transition observed in a kinetic Ising model in a periodically oscillating magnetic fieldKorniss et al. (2000), which also experiences a non-equilibrium phase transition from a symmetry-restoring oscillations to symmetry-breaking oscillations for the time-dependent magnetization, where the dynamical critical behavior is also found to be the same as those of the undriven casesFujisaka et al. (2001). In spite of the similarity of this modelKorniss et al. (2000) and our model, there is an important difference that in our model, instead of external field, the periodic driving is imposed on the strength of interaction, which give rise to the consequence of the spontaneous breaking of discrete TTS. While in model with oscillating fieldKorniss et al. (2000), the period of the magnetization and external field is always the same, thus the discrete TTS is not broken, corresponding to “synchronized phase” in our language.

Another important issue has not been discussed so far is the quantum effect: what will happen if thermal noises are replaced by quantum fluctuations? In another word, if the Hamiltonian.(1) is a quantum transverse Ising model, can we observe DTC behavior in such a periodical-driven quantum system? The effect of quantum fluctuations on systems with continuous TTS breaking has been studied in the single-particle level. For instance, for a Van der Pol oscillator with a classical trajectory of limiting circle, quantum fluctuation induces a phase diffusion in the limiting circleNavarrete-Benlloch et al. (2017), which recovers the continuous TTS symmetry thus destroys the long-range TC order. On the contrary, our model breaks the discrete TTS, which is expected to be more stable compared to the continuous ones. The competition between the quantum fluctuations and time crystalline order may give rise to rich dynamical phenomena, thus is worthwhile to be explored in the future work, even though a numerical simulation of such a non-equilibrium 3D (open) quantum many-body system is a formidable, if not impossible, task.

Conclusion and outlook – In conclusion, we study the stability of DTC phase in the presence of thermal bath, and find that despite its genuine non-equilibrium feature, this DTC phase shares a lot of common properties with equilibrium ordered phases with spatially $Z_{2}$ symmetry breaking, especially its finite-size effect and critical behavior. Future developments will include an analytic explanation of the critical behavior, which calls for a coarse-grained effective description and a non-equilibrium field theory analysis of our system in the Keldysh formalismMartin et al. (1973); Kamenev (2011); Altland and Simons (2010). Recently, Natsheh et al analytically studied the critical behavior of a periodically driven quantum O(N) model with DTC phaseNatsheh et al. (2021), a generalization of this work may shed light on our problem. Another important question involves the generality of our conclusion: whether similar behavior can be observed in other DTC models (e.g. the coupled driven-dissipative nonlinear pendulumsYao et al. (2020))? Last but not the least, whether it is possible to find the melting transitions with the dynamical universality that go beyond the undriven cases? After all, the Magnus expansion could generate long-range interactions in the effective Hamiltonian, especially in the presence of low-frequency driving, and it is possible that such long-range interaction may qualitatively change the critical behavior of the systems and give rise to intriguing dynamical universality classes.

Acknowledgments.—This work is supported by the National Key Research and Development Program of China (Grant No.2020YFA0309000), Natural Science Foundation of China (Grant No.12174251), Natural Science Foundation of Shanghai (Grant No.22ZR142830), Shanghai Municipal Science and Technology Major Project (Grant No.2019SHZDZX01).

Supplemental material

In this supplementary material, we first provide some details of the Heun algorithm used in our simulation, then we numerically check the convergence of our results with the finite discrete time step $\Delta t$. Finally, we discuss the dependence of our results on the initial states as well as the finite size effect in our simulations.