Local Dephasing Sustains Continuous Time Crystals in a Spin-1/2 Ensemble

Introduction

Continuous time crystals extend spontaneous ordering to persistent collective motion [1, 2, 3, 4, 5]. Unlike periodically driven discrete time crystals [6, 7, 8, 9, 10, 11, 12, 13], their oscillations arise under time-independent dynamics. Dissipation can support this order, extending its constructive role in quantum-state engineering [14, 15]. Collective dissipation couples constituents through shared environmental channels [16, 17, 18, 19]; with independent local dissipation, interactions can instead organize their motion. Such interaction-based oscillations have been explored in spin and atomic ensembles [20, 21, 22, 23, 24, 25], photonic systems [26, 27, 28, 29], and coupled ensembles [30, 31, 32]. Experimental platforms include atom-cavity systems [33, 34, 35], Rydberg gases [36, 37, 38, 39], rare-earth-doped solids [40], and spin ensembles [41, 42].

Many routes to dissipative collective oscillations employ population-changing processes, such as spontaneous decay or incoherent pumping, that directly redistribute local occupations. Pure dephasing presents a different possibility. It is ubiquitous in realistic spin ensembles and is generally regarded as unwanted decoherence: it suppresses coherence in its preferred basis without directly changing local populations or providing an obvious excitation channel. Although dephasing can strongly modify interacting nonequilibrium dynamics [43, 44, 45], its potential to sustain autonomous temporal order remains unclear. Can local pure dephasing itself become essential for continuous time-crystalline order?

Here we demonstrate how local pure dephasing sustains collective temporal order in a fully connected ferromagnetic spin-1/2 ensemble with independent local decay. Our results comprise four main contributions. (i) Dephasing-driven oscillations and their mechanism. An exact operator identity separates decay-driven relaxation of a shifted energy observable from a dephasing-induced source proportional to transverse pair correlations. Within homogeneous mean-field theory, we prove that periodic collective oscillations require pure dephasing when local decay is present. Energy charging alone, however, does not sustain motion: the finite energy-relaxation time generates delayed feedback on the transverse magnetization. We show that the cycle-integrated work of this feedback exactly compensates transverse damping along a periodic orbit. (ii) Quantum signatures of time-crystalline order. Beyond mean field, the same operator identity determines the relative amplitude and phase of energy and longitudinal magnetization in quantum oscillatory Liouvillian modes. Liouvillian spectra and stationary two-time correlations reveal decay rates decreasing algebraically with system size, frequencies approaching nonzero mean-field values, and appreciable per-spin collective oscillatory weights. These results connect the microscopic energy mechanism to collective quantum signatures of continuous time-crystalline order. (iii) Intertwined oscillations and metastability. Semiclassical stochastic trajectories reveal that the limit cycle coexists with two symmetry-related stable fixed points. Intrinsic fluctuations repeatedly interrupt collective oscillations with periods of residence near either fixed point and subsequently reactivate oscillatory motion. The oscillatory region also participates in switching between the two stationary regions, while occupation-time statistics characterize the resulting intermittency. (iv) Observability and robustness. Semiclassical photon-counting simulations show how the longitudinal oscillations can be resolved in detected spectra, and parameter-noise tests assess their robustness. Thus, a simple and widely used interacting-spin model with only time-independent local dissipative channels displays both sustained temporal order and a rich organization of dynamical attractors.

The paper is organized as follows. Section II introduces the model and mean-field dynamics. Sections III and IV establish the energy-feedback mechanism and its quantum signatures. Section V examines metastable intermittency, Sec. VI addresses observability and robustness, and Sec. VII discusses the implications and concludes.

Collective oscillations enabled by local dephasing

Model

We consider a homogeneous ensemble of $N$ spin-1/2 particles with all-to-all ferromagnetic interactions along $x$ and a uniform magnetic field along $z$ Figure 1(a). Each spin is subject to local decay and pure dephasing. We denote the Pauli operators for spin $i$ by $\sigma_i^\alpha$, with $\alpha=x,y,z$. Setting $\hbar=1$, the Hamiltonian reads

where $J$ is the ferromagnetic coupling strength and $h_z$ is the magnetic field strength. The factor $1/N$ ensures

that the interaction energy is extensive. In collective-spin notation, $S_\alpha=\frac{1}{2}\sum_i\sigma_i^\alpha$, the Hamiltonian is $H=h_zS_z-(2J/N)S_x^2+J/2$. It is thus a one-axis Lipkin–Meshkov–Glick-type model [46, 47]; the dissipative channels considered here act independently on individual spins. The combination of this interaction with local spin relaxation also underlies dissipative Ising criticality [48].

The density operator $\rho$ evolves under the Liouvillian superoperator $\mathcal{L}$ according to

where $\gamma>0$ and $\gamma_\phi\geq0$ are the local decay and pure-dephasing coefficients, respectively, and $\sigma_i^-=(\sigma_i^x-i\sigma_i^y)/2$ is the spin-lowering operator. The Lindblad dissipator associated with a jump operator $L$ is defined as $\mathcal{D}[L]\rho=L\rho L^\dagger-\{L^\dagger L,\rho\}/2$ [49]. All Hamiltonian parameters and dissipative coefficients are time independent. Local decay sets the longitudinal relaxation rate, $\Gamma_1=\gamma$, while both decay and pure dephasing contribute to the transverse relaxation rate, $\Gamma_2=\gamma/2+2\gamma_\phi$.

Mean-field dynamics

We first analyze the dynamics in the thermodynamic limit within the homogeneous mean-field approximation [50], in which we neglect connected correlations between distinct spins and assume identical local expectation values for all spins. Defining the Bloch-vector components as $m_\alpha=\langle M_\alpha\rangle$, where $M_\alpha=\frac{1}{N}\sum_{i=1}^{N}\sigma_i^\alpha$ are the normalized magnetization operators, we obtain the nonlinear Bloch equations

Figure 1 shows how local dephasing turns relaxation into sustained collective motion. Without pure dephasing, the longitudinal magnetization approaches a stationary value [Fig. 1(b)]. Raising $\gamma_\phi$ from zero to 1.5, with the Hamiltonian, $\Gamma_1$, and initial state unchanged, produces a persistent periodic response [Fig. 1(c)]. The added channel thus sustains oscillations despite enhancing transverse damping. This autonomous motion requires neither coherent driving nor incoherent pumping.

Besides the polarized fixed point $(0,0,-1)$, Eqs. (3) admit a symmetry-related pair $F_\pm$ when $2h_zJ>h_z^2+\Gamma_2^2$:

Linearizing the original Bloch equations at either member of this pair gives the characteristic equation

For $\Gamma_1>0$, all roots have negative real parts throughout the pair’s existence domain. The limit cycle therefore coexists with two locally stable fixed points rather than emerging from their Hopf destabilization. The asymptotic response depends on which basin of attraction contains the initial state, as in other time-crystal models with non-Hopf transitions or stationary–periodic coexistence [51, 52].

To characterize this initial-state dependence, we evolve 200 randomly sampled initial Bloch vectors at each parameter point and construct the dynamical response map in Figure 2(a), with $h_z=J$ and fixed $\gamma$. The hue indicates whether periodic or chaotic responses are more frequent, while the color intensity gives the fraction of all 200 trajectories assigned to that category. Nonstationary responses occupy a finite dephasing window, bounded by predominantly stationary regimes at weaker and stronger dephasing. Varying the ferromagnetic coupling shifts and reshapes this window, suggesting that sustained motion requires a suitable balance between interactions and dephasing. Near the boundaries between dynamical regimes, long-lived periodic and chaotic transients also occur. We include these transients in the corresponding periodic or chaotic response categories according to their behavior within the finite observation window. Representative fixed-point, periodic, and chaotic trajectories are shown in Figure 2(b)–Figure 2(d), respectively.

A readily prepared initial state is the fully $x$-polarized product state, $\boldsymbol{m}(0)=(1,0,0)$, obtainable by a global spin rotation of a polarized ensemble [53]. The inset in Figure 2(a) maps the collective response from this fixed initial state, providing a parameter guide for observing periodic oscillations with this preparation.

The Liouvillian $\mathcal{L}$ is invariant under the $\mathbb{Z}_2$ parity transformation $\rho\mapsto P\rho P$, where $P=\prod_{i=1}^{N}\sigma_i^z$. At the mean-field level, this transformation acts as $(m_x,m_y,m_z)\mapsto(-m_x,-m_y,m_z)$. For the symmetric periodic branch studied here, the parity transformation is equivalent to a half-period time shift along the orbit,

where $T$ is the full orbital period. The transverse and longitudinal signals have periods $T$ and $T/2$, corresponding to fundamental angular frequencies $\Omega=2\pi/T$ and $2\Omega$, respectively. To understand how local dephasing can sustain collective oscillations despite enhancing transverse damping, we next examine the energy balance of the interacting ensemble.

From local dephasing to collective energy feedback

The physical question is how a channel that damps transverse coherence can sustain collective motion. We first establish an exact energy balance for the quantum model, then use its homogeneous mean-field limit to prove the necessity of dephasing and identify the delayed feedback that compensates damping.

Exact interaction-energy balance

Define the transverse pair-correlation operator and a shifted, rescaled energy operator by

where $I$ is the identity operator. The local dissipators give $\mathcal{L}_{\mathrm{diss}}^\dagger(M_z+I)=-\Gamma_1(M_z+I)$ and $\mathcal{L}_{\mathrm{diss}}^\dagger U_N=-2\Gamma_2U_N$. Since $[H,Q_N]=0$, the coherent contribution to the energy derivative vanishes. Using $2\Gamma_2=\Gamma_1+4\gamma_\varphi$, we obtain the exact identity

The first form separates relaxation of $Q_N$ from an interaction-energy source carried by transverse pair correlations. The second expresses the same response in terms of longitudinal magnetization and will determine the quantum-mode relations in Sec. IV B.

This identity has a simple microscopic interpretation. Parallel transverse spins, whether aligned along $+x$ or $-x$, lower the ferromagnetic interaction energy. Dephasing in the $z$ basis suppresses the local coherence supporting that alignment and reduces the transverse pair correlations. For $e_N=\langle H\rangle/N$, its direct contribution is

Thus, for positive transverse ferromagnetic correlations, dephasing raises the interaction energy even though it does not directly change the local $z$ populations. The source is state dependent at finite $N$; the exact identity does not assert positive charging for arbitrary pair correlations. Its role becomes especially transparent in homogeneous mean field, where $\langle U_N\rangle\to m_x^2$.

Mean-field energy balance and the necessity of dephasing

The Hamiltonian energy per spin provides the natural starting point. Under homogeneous factorization, $\langle U_N\rangle\to m_x^2$ and $\langle Q_N\rangle\to Q$. The Hamiltonian in Eq. (1) gives

The factor $1/2$ in the interaction term follows from the $N(N-1)/2$ distinct pairs together with the $1/N$ Kac normalization. To remove the constant offset in the local-decay contribution to the energy balance, we introduce the shifted and rescaled energy variable

With this choice, local decay contributes the simple relaxation term $-\Gamma_1Q$. Substituting Eqs. (3) into the time derivative of Eq. (11) gives

This is the mean-field limit of the exact identity (10). Local decay relaxes $Q$ toward zero, whereas dephasing supplies a positive source proportional to $m_x^2$. On the periodic response considered below, $Q\geq0$, so Eq. (12) directly separates interaction-energy charging from energy loss.

The same variable feeds back on the transverse motion. Using

to eliminate $m_y$ and $m_z$ from the Bloch equations gives

Equation (13) describes a damped nonlinear oscillator coupled to the relaxing energy variable $Q$. Thus $m_x^2$ drives the energy variable through Eq. (12), while the energy variable acts back through the force $-2JQm_x$ in Eq. (13).

This form provides a direct proof that pure dephasing is necessary for periodic motion in the homogeneous mean-field dynamics of this model when $\Gamma_1>0$. At $\gamma_\varphi=0$, Eq. (12) gives $Q(t)=Q(0)e^{-\Gamma_1t}$, so periodicity requires $Q\equiv0$. When Eq. (13) is multiplied by $\dot{m}_x$ and integrated over a period $T$, all conservative terms reduce to endpoint differences and vanish. Hence

Because $\Gamma_2=\Gamma_1/2>0$, $m_x$ must be constant; the relations above then imply that $m_y$ and $m_z$ are constant as well. With nonzero local decay, pure dephasing is therefore necessary for any periodic orbit of the homogeneous mean-field dynamics.

Delayed feedback and sustained oscillations

Charging and loss can balance at a fixed point, so their presence alone does not explain sustained motion. What distinguishes the periodic response is that the energy follows the transverse order with a finite relaxation time and therefore acts back with a delay. Solving Eq. (12) with an integrating factor from an arbitrary initial time $t_0$ gives

For the bounded long-time periodic response, the initial term has decayed. Setting $s = t-t'$ gives

The kernel weights the recent history of $m_x^2$ over the memory time $\Gamma_1^{-1}$. It also makes $Q \ge0$ on the periodic orbit. Averaging Eq. (16) over a period $T_Q$ shared by $Q$ and $m_x^2$ gives the global charging–loss balance

The half-period symmetry in Eq. (8) makes $Q$ and $m_x^2$ repeat after $T_Q=T/2$.

This balance fixes the cycle-averaged energy, while the relative timing of $Q$ and $m_x^2$ within the cycle is determined by the memory kernel in Eq. (16). To quantify this temporal response, we expand both quantities in Fourier harmonics and evaluate the memory integral for each component with time dependence $e^{i\nu t}$, obtaining

At $\nu=0$, this relation reproduces the cycle-averaged balance in Eq. (17).

At fixed $J$, $\Gamma_1$, and $\nu>0$, the response gain is proportional to $\gamma_\phi$, whereas the finite energy-relaxation time produces the phase lag $\delta_\nu$.

Figure 3(a) resolves this phase lag in the time domain, while Figure 3(b) displays the associated loop in $(m_x^2,Q)$. At a fixed value of $m_x^2$, $Q$ is smaller on one branch of the cycle and larger on the other. Since the transverse feedback force is $-2JQm_x$, this asymmetry produces a nonzero net contribution over a closed cycle.

Multiplying Eq. (13) by $\dot m_x$ and integrating over one feedback period of the symmetric orbit gives the global balance

The positive sign follows directly by setting $u=m_x^2$ and using Eq. (12). Integration by parts over a cycle gives

where $\langle Q\dot Q\rangle_{T_Q}=0$ by periodicity. Thus the delayed response supplies strictly positive feedback work on a periodic orbit. The left-hand side is the transverse damping accumulated over the cycle, while the right-hand side is the generalized work of the internal feedback force. Equation (19) is therefore a cycle-integrated description of sustained motion: the oriented loop in Figure 3(b) exactly compensates the damping of the reduced transverse dynamics. The physical energy charging and loss themselves remain those in Eq. (12).

If the energy followed $m_x^2$ instantaneously, $Q=(4J\gamma_\phi/\Gamma_1)m_x^2$, the path in the $(m_x^2,Q)$ plane would retrace a line and the contour integral in Eq. (19) would vanish. Dephasing-induced charging is therefore necessary, but finite response time is what converts that charging into the delayed feedback capable of sustaining the cycle.

This mechanism characterizes the periodic orbit once it exists; its emergence, coexistence with stationary states, and basin structure are governed by the full nonlinear dynamics of Eqs. (3).

Together, the necessity proof and the cycle balance explain how local dephasing, acting through ferromagnetic interactions, can sustain collective oscillations despite enhancing transverse damping. The exact identity (10) carries this energy-response structure into the quantum analysis below.

This interaction-energy mechanism complements symmetry-based routes to nonstationarity and synchronization [54, 55, 56, 57, 58] and dephasing-assisted temporal order under periodic driving [59].

Finite-size quantum dynamics and time-crystalline order

Continuous time-crystalline order requires nondecaying collective oscillations at a nonzero frequency in the thermodynamic limit, reflecting spontaneous breaking of continuous time-translation symmetry. To characterize this order beyond mean field, we identify the oscillatory modes of the Liouvillian $\mathcal{L}$ that contribute to collective two-time correlations and analyze three aspects of their system-size dependence: decay rates, oscillation frequencies, and collective oscillatory weights. We then investigate the energy response associated with these modes, connecting the finite-size quantum dynamics to the feedback mechanism developed in Sec. III.

Finite-size scaling and temporal order

Exploiting permutational invariance, we construct a Dicke-block representation of the Liouvillian $\mathcal{L}$ defined in Eq. (2) and analyze its eigenmodes [19]. A nonstationary right eigenmatrix $R_k$ satisfies

where $\kappa_k$ and $\omega_k$ denote the decay rate and angular frequency, respectively. Trace preservation implies $\operatorname{Tr}R_k=0$ for every nonzero eigenvalue. An oscillatory mode ($\omega_k\ne0$) evolves with time dependence $e^{-\kappa_k t+\mathrm{i}\omega_k t}$ and has a lifetime $\kappa_k^{-1}$ for $\kappa_k>0$.

The parity symmetry introduced in § 1.1 allows the modes to be classified into odd and even sectors. The symmetric mean-field orbit in Eq. (8) provides the reference frequencies $\Omega$ and $2\Omega$ for the transverse and longitudinal oscillations, respectively. To assess continuous time-crystalline order, we examine how the corresponding oscillatory branches evolve with system size. Over the studied range $N = 80$–$300$, both branches become longer lived, with the power-law fits in Figure 4(a) yielding $\kappa_{\mathrm{odd}} \propto N^{-0.323}$ and $\kappa_{\mathrm{even}} \propto N^{-0.458}$. Their frequencies move toward the mean-field values $\Omega$ and $2\Omega$, respectively, as shown in Figure 4(b).

To establish how these modes contribute to collective fluctuations, we calculate stationary connected two-time correlations using the quantum regression theorem (QRT) [60, 61]. For observables $A$ and $B$ and delay $\tau\ge0$, we define

where $\rho_{\mathrm{ss}}$ is the stationary density operator, $I$ is the identity operator, $\delta A = A - \langle A\rangle_{\mathrm{ss}} I$, and $\langle A\rangle_{\mathrm{ss}} = \operatorname{Tr}(A\rho_{\mathrm{ss}})$; $\delta B$ is defined analogously. For $A = B = M_\alpha$, we write $G_{\alpha\alpha}$, with the per-spin normalization already included in $M_\alpha$. These correlations probe collective fluctuations even when the stationary mean magnetization is time independent.

For a diagonalizable Liouvillian with a complete eigenbasis, the magnetization autocorrelations can be expressed as a sum of Liouvillian modes,

where $w_{\alpha k}$ denotes the generally complex spectral weight of mode $k$ in the $\alpha$ channel. For a parity-invariant stationary state, only odd-parity modes contribute to $G_{xx}$, whereas only even-parity modes contribute to $G_{zz}$.

Slow Liouvillian modes describe relaxation and metastable dynamics [62, 63, 64, 65]. To isolate the collective oscillations, we separate each correlation into a slowly decaying background and an oscillatory contribution,

The backgrounds $B_x(\tau)$ and $B_z(\tau)$ are obtained from fits to the long-delay correlations. Their assignments to nonoscillatory transverse and low-frequency longitudinal modes, respectively, are supported by independently computed Liouvillian spectra. After background subtraction, the transverse and longitudinal correlations in Figure 4(c,Figure 4d) exhibit damped oscillations whose decay rates and frequencies track the odd- and even-parity branches identified above. This correspondence connects the spectral modes to oscillations of the collective magnetization. As $N$ increases, these correlations persist over longer delays, while the spectral peaks shown in the insets narrow and shift toward the mean-field frequencies.

To quantify the collective oscillatory weights on the absolute per-spin scale, we fit the extracted correlations to

Here $A_\alpha(0)$ is the fitted oscillatory envelope amplitude extrapolated to zero delay; $\kappa_\alpha$, $\omega_\alpha$, and $\varphi_\alpha$ are the decay rate, angular frequency, and phase offset of the oscillatory contribution.

Figure 5 shows these amplitudes for $N = 80$–$250$. The transverse central estimate changes from approximately $5.4 \times10^{-2}$ to $4.1 \times10^{-2}$. The longitudinal estimate changes from $3.9 \times10^{-3}$ to $5.5 \times10^{-3}$, reaching $6.8 \times10^{-3}$ at $N = 150$. Within each channel, the oscillatory weight remains of comparable absolute magnitude across the studied sizes. The error bars span the amplitude estimates obtained using different fitting windows, with the larger longitudinal spans reflecting sensitivity to the background separation and oscillatory fitting interval.

Continuous time-crystalline order requires an oscillatory contribution whose decay rate vanishes in the thermodynamic limit, while its frequency and weight in the per-spin collective correlation remain nonzero [3, 4, 51]. Here the thermodynamic limit is taken before the long-delay limit. The algebraic decrease of the fitted decay rates, convergence toward nonzero mean-field frequencies, and comparable per-spin oscillatory amplitudes jointly support dephasing-induced continuous time-crystalline order.

Energy response of collective oscillations

The exact identity (10) now connects the collective oscillatory modes to the interaction-energy feedback Sec. III. Because the identity holds at finite $N$, this connection does not require factorization of the pair correlations.

To resolve this energy response within the oscillatory modes identified in Sec. IV A, we define

Since $Q_N$, $U_N$, and $M_z$ are parity even, we evaluate these projections for the even-parity oscillatory branch. For nonzero $u_k$ and $z_k$, projecting Eq. (10) onto $R_k$ and using $\operatorname{Tr} R_k=0$ gives the relative energy responses

These ratios are independent of the normalization of $R_k$. The first extends the mean-field harmonic response in Eq. (18) to a quantum oscillatory mode, with $i\nu$ replaced by $\lambda_k=-\kappa_k+i\omega_k$. The finite-size decay rate therefore enters the response denominator together with the oscillation frequency, determining both the relative energy amplitude and its phase. The second ratio expresses the same energy response relative to the longitudinal magnetization.

Figure 6(a,b) shows the relative energy amplitude and phase obtained from the even-parity eigenmodes, in agreement with Eq. (26). Across the studied sizes, the relative energy amplitude increases, accompanied by an overall reduction in the magnitude of the phase lag. For $\omega_k>0$ and $\gamma_\phi>0$, the energy phase relative to the longitudinal magnetization is

The phases below $-90^\circ$ in Figure 6(b) arise from $\kappa_k>2\Gamma_2$, which makes the real part of the response denominator negative. The phase shift thus reflects the energy response within a decaying collective oscillation. As $\kappa_k$ approaches zero, Eq. (26) reduces to the purely oscillatory response, connecting the finite-size modes to the harmonic description in Sec. III.

The same energy-response relation also governs the stationary two-time correlations. Using Eq. (22), we denote the energy–magnetization and pair–magnetization correlations by $G_{Qz}$ and $G_{Uz}$, corresponding to $(A,B)=(Q_N,M_z)$ and $(U_N,M_z)$, respectively. Equation (10) gives

Its solution separates the relaxation of the initial covariance from the response to longitudinal correlations,

The convolution makes explicit how the energy correlation responds to the history of the longitudinal collective fluctuations. The pair–magnetization correlation follows from $G_{Uz}=(h_zG_{zz}-G_{Qz})/J$. Figure 6(c) shows these energy and pair responses at $N=200$, reconstructed from the full longitudinal correlation and the stationary initial covariance. Together with the mode projections, these time-domain responses show how the energy-response structure developed in Sec. III is encoded in the collective quantum oscillations.

Intertwined time-crystalline dynamics and metastability

The mean-field limit cycle coexists with two stable fixed points. We now examine how intrinsic fluctuations connect these attractors, first through semiclassical trajectories and then through fixed-point occupation statistics. Figure 7 resolves repeated capture near fixed points and reactivation of collective oscillations within this three-region dynamics.

Noise-induced switching between stationary and oscillatory states

The mean-field periodic attractor identified in § 1.1 coexists with two symmetry-related stable fixed points, $F_+$ and $F_-$. The labels follow the sign of $m_x^F$ in Eq. (6): $m_x^F>0$ at $F_+$ and $m_x^F<0$ at $F_-$. To resolve how fluctuations connect these coexisting attractors, we describe the magnetization using a semiclassical diffusion approximation that combines the mean-field drift with Gaussian noise from the local dissipative channels [66, 67]:

Equation (30) is interpreted in the Itô sense. Here $\boldsymbol{m} = (m_x,m_y,m_z)$ denotes the random magnetization in a single semiclassical realization, $\boldsymbol{F}$ is the vector field in Eqs. (3), and $\boldsymbol{W}_t$ contains three independent Wiener processes satisfying $\mathbb{E}[dW_\alpha dW_\beta] = \delta_{\alpha\beta}\,dt$. The expectation $\mathbb{E}[\cdot]$ averages over stochastic realizations. The factor $N^{-1/2}$ reflects the $1/N$ scaling of the noise variance with system size, resulting from independent local noise contributions to the normalized collective magnetization. This collective diffusion model resolves fluctuation-driven excursions around the mean-field dynamics; it is not an exact quantum-jump unraveling of Eq. (2). Related phase-space approaches provide complementary tools for large dissipative spin ensembles [68, 69].

We obtain the diffusion matrix by evaluating the local dissipators in the symmetrized second-moment equations and subtracting the drift contributions [66, 70]:

For the parameters in Figure 7, we initialize the dynamics in the mean-field limit-cycle basin and follow trajectories at $N = 10^5$. Here $T$ and $\Omega= 2\pi/T$ refer to the noiseless cycle. Figure 7(a,Figure 7b) show repeated alternation between oscillatory motion and residence near either fixed point. Capture is temporary: intrinsic fluctuations can reactivate oscillations, as the complete $4000T$ records for three initial conditions confirm [Fig. 7(c)]. We also observe reactivation from an initial state placed directly at a stable mean-field fixed point.

The limit-cycle region has a dual role. It supports many rapid oscillations as a metastable dynamical region, and it provides an observed intermediate channel between the stationary regions: escape from one fixed-point neighborhood can be followed by an extended oscillatory episode and capture near the other. These paths intertwine fast temporal order with slower transitions between regions. This is distinct from phase diffusion, which can erase phase memory without leaving the oscillatory region; reactivation likewise need not restore the earlier phase. The observed switching channel involves the cycle’s neighborhood, without implying that every noisy path must cross the deterministic orbit.

Coexistence and switching are established phenomena, but the relevant precedents involve different dynamical structures. Stationary-state switching has been studied in Rydberg, photonic, and Dicke systems [66, 70, 71, 72, 73]. In particular, Stitely et al. have resolved switching among a normal state and two symmetry-related superradiant states, with oscillatory dynamics studied in a separate regime [74]. Mattes et al. have established stationary–limit-cycle coexistence and hysteresis [52]; quantum-trajectory intermittency at a time-crystal transition [75] and phase-diffusive versus activated relaxation in distinct oscillatory and stationary regimes [51] have also been characterized. Here the distinguishing structure is the recurrent participation of a dephasing-sustained limit cycle in switching between two stationary regions at fixed parameters. Temporal order is therefore part of an intermediate metastable episode, not merely a separate phase or a transient preceding stationary relaxation.

Temporal statistics of metastable intermittency

To quantify this intermittency, we define the fixed-point occupation fraction as the fraction of an observation window spent in identified residence intervals near either $F_{+}$ or $F_{-}$:

where $\chi_{F}$ equals one during identified fixed-point residence intervals and zero otherwise. The remaining intervals predominantly correspond to limit-cycle-associated oscillatory motion, including noise-driven excursions away from the deterministic cycle; they also contain a small number of unresolved brief fixed-point captures.

At fixed $N = 10^{5}$, Figure 7(d) examines how the probability of exceeding a prescribed occupation threshold decreases with observation time. The thresholds of 17%, 20%, and 23% all exceed the sample mean fixed-point occupation fraction of approximately 12.8%. Over the fitting range $150 \leq L/T \leq500$, these upper-tail probabilities exhibit an approximately exponential decrease, consistent with a temporal large-deviation description [76, 77, 78, 79]. The effective decay rate increases with the threshold. Thus, sustaining an unusually large fixed-point occupation fraction over a longer observation time becomes increasingly unlikely. This is a statement about occupation-time fluctuations at fixed $N$, rather than a measurement of escape-rate scaling with system size.

The occupation statistic complements the trajectory geometry. The former quantifies how much of a record is interrupted by fixed-point residence, whereas the latter shows how oscillations resume and participate in switching. Together, they identify a recurrent oscillatory metastable region embedded in broader switching dynamics. Resolving residence times and transition paths could connect this phase-and-region description to effective metastable dynamics and its trajectory dependence [64, 80]. The long oscillatory intervals at $N = 10^{5}$ provide windows for continuous observation, which we examine next.

Observability and robustness

The preceding sections distinguish the mechanism sustaining temporal order, its quantum spectral signatures, and its fluctuating realization among coexisting attractors. We now ask how this dynamics appears in an experimental record. Using the semiclassical trajectories of § 1.4, we examine photon-counting readout of the longitudinal oscillations, then test the response to fluctuations of the interaction and local dissipative rates.

Continuous photon readout

We examine how clearly the longitudinal collective oscillations can be resolved in continuous photon-counting records. Continuous optical monitoring has revealed time-crystalline motion in atom–cavity experiments [33, 35]. In the present model, radiative local spin decay provides a direct probe of the excited-state population fraction, $p_{\uparrow}(t) = [1 + m_{z}(t)]/2$. For photons from independent local decay channels counted with a common efficiency $\eta$, the mean detection rate along a given semiclassical trajectory is [81]

Within the semiclassical photodetection model, the detected photon number $K_{b}$ in each time bin $[t_{b}, t_{b} + \Delta t]$ is Poisson distributed conditional on the population trajectory [82]:

The corresponding estimate of the bin-averaged population is

Writing the population averaged over the bin as $\bar{p}_{\uparrow,b}$, the conditional counting variance is

Thus, for a fixed population fraction, detection efficiency, and bin width, the standard deviation of the counting error decreases as $N^{-1/2}$. Here the conditional Poisson sampling is a semiclassical readout model applied to the magnetization record, rather than a simulation of the exact measurement-conditioned density operator.

Figure 8 compares the population trajectories and their photon-counting readout at $N=10^4$, $10^5$, and $10^6$, using the same initial preparation, detection efficiency $\eta=0.3$, and bin width $\Delta t=0.08T$. To distinguish the spectral structure of the underlying dynamics from the additional counting background, we evaluate the power spectra of the bin-averaged population and its count estimate, denoted $S_{\mathrm{traj}}$ and $S_{\mathrm{count}}$. Both spectra use the same $128T$ record following an $8T$ preparation interval and are computed with the Welch method [83], using $32T$ Hann segments with 50% overlap. Their common normalization in absolute population units allows the oscillatory peaks and backgrounds to be compared directly.

The particle-number dependence reflects the combined effects of intrinsic finite-size fluctuations and photon-counting noise. At $N=10^4$, intrinsic fluctuations strongly disrupt the regularity of the population dynamics, leaving no well-resolved oscillation peak even in the trajectory spectrum. The absence of a clear spectral signature therefore originates in the underlying dynamics, rather than solely in the detection process.

At $N=10^5$, the weaker intrinsic fluctuations allow a more regular oscillatory signal to develop, but photon shot noise still largely obscures individual cycles in the time-domain record. Spectral analysis nevertheless resolves a peak near $2\Omega$ approximately an order of magnitude above the surrounding background. This distinction arises because the oscillatory power is concentrated near its characteristic frequency, whereas counting noise contributes a broadband background [81]. The collective oscillation can therefore be identified spectrally even when its individual cycles are difficult to distinguish directly.

At $N=10^6$, the relative counting fluctuations are further suppressed, and the detected record closely follows the underlying population modulation. The agreement between the two spectra around the dominant peak likewise shows that photon counting faithfully preserves the oscillatory signal. For the detection conditions simulated here, shot noise thus has little effect on resolving the collective oscillations in either domain.

Robustness to experimental parameter noise

In addition to the intrinsic fluctuations considered above, noise in experimental control parameters can affect the stability of collective oscillations. To assess their robustness against such noise, we introduce stochastic variations of the interaction strength $J$, local decay rate $\Gamma_1$, and pure-dephasing rate $\gamma_{\phi}$ into the mean-field Bloch equations. We vary one parameter at a time, keeping the remaining independent parameters fixed.

To keep the dissipative rates positive and preserve the ferromagnetic sign of the interaction, we use a common lognormal modulation,

where $a_0$ is the noiseless parameter value and $s_a$ is the logarithmic modulation amplitude. The stationary Ornstein–Uhlenbeck process $u_a(t)$ has zero mean, unit variance, and correlation $\mathbb{E}[u_a(t)u_a(t')]=e^{-\lvert t-t'\rvert/\tau_c}$ [84].

The exponential form guarantees $a(t) > 0$, while the offset $-s_a^2/2$ ensures $\mathbb{E}[a(t)] = a_0$. Thus, the modulation introduces parameter fluctuations while preserving the mean parameter values.

For correlation times satisfying $\omega\tau_c \ll1$, the parameter-noise spectrum is approximately flat over the relevant frequency range. We choose $\tau_c = 0.01T$, where $T$ is the full period of the noiseless mean-field cycle, so that the fluctuations approximate white noise over the frequency range relevant to the collective oscillations.

The noise amplitude is specified by the relative standard deviation,

and is reported as a percentage. At equal $r_a$ and $\tau_c$, the three channels have the same relative noise statistics, allowing their effects on the collective oscillations to be compared. The transverse relaxation rate follows $\Gamma_2(t) = \Gamma_1(t)/2 + 2\gamma_\phi(t)$, while time and frequency are expressed using the noiseless value of $\Gamma_2$. For each noise strength, $S_{zz}(\omega)$ in Figure 9 is the average of the longitudinal magnetization power spectra from 100 independent realizations.

The different noise responses in Figure 9 emphasize that robustness depends on which part of the feedback mechanism is perturbed. These tests concern fluctuations around an operating point with nonzero dephasing; they complement, rather than replace, the necessity result obtained by removing dephasing altogether.

Discussion and conclusions

Our results identify a route from a ubiquitous decoherence channel to autonomous temporal order in a familiar interacting spin model. Pure dephasing is energy neutral for isolated spins in the dephasing basis, but not for their transverse interaction energy: the loss of ferromagnetic correlations supplies energy that local decay subsequently removes. The key is not heating alone, but its delayed return as a collective feedback force.

Equation (18) makes this conversion explicit: the cycle-averaged feedback power is proportional to $\langle\dot{Q}^{2}\rangle$ and compensates damping, whereas instantaneous energy response produces no net work. The resulting design principle is to identify energy-bearing correlations, a dissipative channel that modifies them, and a relaxing internal variable that feeds back on collective motion. Unlike stabilization through reservoir memory [85], the delay here arises within dynamics generated by a time-independent Markovian master equation. This perspective suggests treating dephasing strength and basis as control parameters rather than simply minimizing decoherence. The exact amplitude–phase relation provides a diagnostic for pursuing this mechanism beyond fully connected ensembles, where spatial fluctuations can reshape temporal order [86]. Metastability adds a complementary control problem: an oscillatory attractor can be repeatedly reactivated and can participate in transitions between stationary regions. For sustained readout and time-crystal sensing [87, 88], controlling access to that attractor can therefore matter alongside stabilizing its oscillation frequency.

We have shown that local pure dephasing can sustain continuous time-crystalline order in a spin-1/2 ensemble with local decay. We have established an exact energy identity, proved the necessity of dephasing for periodic motion in the homogeneous mean-field model, and identified delayed energy feedback as the mechanism that balances damping. We have connected this mechanism to quantum Liouvillian modes and collective correlations, uncovered recurrent oscillatory episodes within metastable switching, and characterized photon-counting readout and robustness. These results have established a constructive role for local coherence loss and identified a route to autonomous temporal order through a dynamical loop of dephasing-induced energy charging, decay-driven relaxation, and delayed feedback on the collective magnetization.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grants No. 12375013 and No. 12547109) and the Guangdong Provincial Quantum Science Strategic Initiative (Grant No. GDZX2503008).