# Self-Organized Time Crystal in Driven-Dissipative Quantum System

Ya-Xin Xiang  Qun-Li Lei  National Laboratory of Solid State Microstructures
and School of Physics, Collaborative Innovation Center of Advanced
Microstructures, Nanjing University, Nanjing 210093, China  Zhengyang Bai
<EMAIL_ADDRESS>State Key Laboratory of Precision Spectroscopy, East
China Normal University, Shanghai 200062, China  Yu-Qiang Ma
<EMAIL_ADDRESS>National Laboratory of Solid State Microstructures and
School of Physics, Collaborative Innovation Center of Advanced
Microstructures, Nanjing University, Nanjing 210093, China

###### Abstract

Continuous time crystals (CTCs) are characterized by sustained oscillations
that break the time translation symmetry. Since the ruling out of equilibrium
CTCs by no-go theorems, the emergence of such dynamical phases has been
observed in various driven-dissipative quantum platforms. The current
understanding of CTCs is mainly based on mean-field (MF) theories, which fail
to address the problem of whether the long-range time crystalline order exists
in noisy, spatially extended systems without the protection of all-to-all
couplings. Here, we propose a new kind of CTC realized in a quantum contact
model through self-organized bistability (SOB). The exotic CTCs stem from the
interplay between collective dissipation induced by the first-order absorbing
phase transitions (APTs) and slow constant driving provided by an incoherent
pump. The stability of such oscillatory phases in finite dimensions under the
action of intrinsic quantum fluctuations is scrutinized by the functional
renormalization group method and numerical simulations. Occurring at the edge
of quantum synchronization, the CTC phase exhibits an inherent period and
amplitude with a coherence time diverging with system size, thus also
constituting a boundary time crystal (BTC). Our results serve as a solid route
towards self-protected CTCs in strongly interacting open systems.

††preprint: APS/123-QED

_Introduction._ —Time crystals are self-organized spatiotemporal structures,
first envisaged by Wilczek [1, 2], that spontaneously break the time-
translation symmetry imposed by the underlying Hamiltonians. Manifestations of
this idea include discrete time crystals [3, 4, 5, 6, 7] with subharmonic
responses under periodic driving, and CTCs [8, 9, 10, 11] with self-
oscillations out of continuous driving.

Since the advent of the no-go theorems stating that it is impossible to
observe a spontaneously oscillating ground-state (in thermal equilibrium) [12,
13], several efforts have been devoted to the study of CTCs in dissipative
quantum systems. By building up a limit cycle (LC), the arise of
synchronization in different types of physical platforms, such as
optomechanical oscillators [14], Rydberg gases [15, 16, 17] and hybrid atom-
cavity systems [18, 19], has been observed and related to the formation of
CTCs. Given that in open systems, the dissipation often associates with the
quantum Langevin noise, it is probable that the fluctuations would affect the
robustness of CTCs thereby destroying the crystalline order. Notwithstanding
the rapid advances in experimental studies, to what extend the CTCs predicted
by MF analysis remain intact under the action of intrinsic noise is an open
question worthy of theoretical endeavors.

Analogous to the famous notion of self-organized criticality (SOC), which is
related to self-organization to the critical point of a continuous APT [20,
21, 22, 23, 24, 25, 26, 27], the mechanism for SOB consists in a separation of
the time scale of the dynamics of the order parameter from that of the
corresponding control parameter. It triggers a LC phase of the hysteresis loop
of a first-order APT [28, 27]. In light of the common features shared by LCs
and CTCs, a new class of CTCs induced by SOB can be envisioned.

In this work, we theoretically investigate the formation and stability of CTCs
beyond MF approximation. Concretely, we consider a general contact model
characterized by the quantum and classical contact interactions between
quantum emitters [29]. In the classical regime, the system undergoes
continuous APTs. However, the transitions become discontinuous in the quantum
regime [29]. Being inherently dissipative, the system could support SOB if a
loading mechanism is added. The combined effects of the slow reloading and
fast collective dissipation bring about a non-stationary phase out of a
constant dynamical generator, where the average densities of total and active
quantum emitters change periodically, and in the same time, the system
undergoes repeated phase transitions [see Fig. 1(b)-(d)], and self-organizes
to a CTC phase.

We further investigate their robustness against the intrinsic quantum noise
and spatial inhomogeneities by theoretical analysis and numerical simulations,
and find that our CTCs exhibit strong dimensionality and system-size
dependence, reflecting the underlying first-order APTs and quantum
synchronization. Possible experimental implementations are briefly discussed
in the end.

Figure 1: (a) The effecitve four-level scheme. The quantum emitter in the
inactive state $|i\rangle$ (gray sphere) in proximity to emitters in the
active state $|a\rangle$ (large red sphere) can become active via the
processes of (in)coherent facilitated activation, and active emitters can
spontaneously decay into either inactive or removed states $|0\rangle$ (blue
sphere). Emitters in the $|p\rangle$ state (purple sphere) are incoherently
pumped to the inactive state $|p\rangle\to|i\rangle$. (b) Sketch of an SOB-
induced CTC, where the contact activation occurs on the facilitation shell
(large blue hollow sphere) of active emitters. The perpetual excursions into
the bistable region typically consist of subcritical (I) and supercritical
(II) states of low active densities, and supercritical highly active states
(III), corresponding to the purple, green, and orange open circles in panel
(c). (c) Time series of the average total $n$ and active $\rho$ densities,
where the dashed lines are equidistantly separated by half the period $T$. (d)
Snapshots of the active density field from simulations are sketched.
Parameters are $\kappa=0,\Omega=0.5,\lambda=3.2\times 10^{-3},L=64$.

_Model._ —We consider a general contact model adapted from Ref. [29]. As
depicted in Fig. 1(a), the model can be represented as an effective four-level
system, where quantum emitters in the active state $|a\rangle$ can
spontaneously decay into the inactive state $|i\rangle$ (with rate $\Gamma$),
and the inactive ones can be activated only in the vicinity of active ones
both incoherently and coherently (with rates $\kappa$ and $\Omega$,
respectively). In addition to the aforementioned contact processes, loss of
emitters due to the decaying of active emitters into the removed states
$|0\rangle$ (with rate $b\Gamma$), and an incoherent coupling
$|p\rangle\to|i\rangle$ that mimics injecting inactive emitters (with rate
$\lambda$) are added. We further restrict the contact processes to pairs of
emitters that are separated by a certain distance (the facilitation radius
$R_{\rm fac}$) from each other, and refer to the effective nearest neighbors
of an active emitter as emitters at the border of its facilitation sphere
[30]. Consequently, in free space, the propagation of activation occurs within
the facilitation shells [large blue hollow spheres in Fig. 1(b)].

Under the Markovian noise, the effective dynamics of this system permits a
microscopic description for the density operator $\hat{\rho}$ via a Lindblad
master equation
$\partial_{t}\hat{\rho}=-i[\hat{H},\hat{\rho}]+\sum_{\alpha}\mathcal{L}_{\alpha}\hat{\rho}$.
The coherent activation is described by the effective Hamiltonian
($\hbar\equiv 1$ henceforth)

$\hat{H}=\Omega\sum_{l}{\hat{C}}_{l}{\hat{\sigma}}_{l}^{x},$ (1)

where ${\hat{C}}_{l}=\sum_{k\in\partial
l}{{\hat{\sigma}}_{k}^{aa}},{\hat{\sigma}}_{l}^{\alpha\beta}\equiv|\alpha_{l}\rangle\langle\beta_{l}|$
($\alpha,\beta=a,i,p,0$), $l,k$ are indices for each emitter, and
$\sum_{k\in\partial l}$ denotes a summation of the effective nearest neighbors
of the $l$-th emitter. The operator
${\hat{\sigma}}_{l}^{x}=\hat{\sigma}^{-}_{l}+\hat{\sigma}^{+}_{l}$ flips the
quantum state with the ladder operators
$\hat{\sigma}^{+}_{l}\equiv\hat{\sigma}_{l}^{ai}$ and
$\hat{\sigma}^{-}_{l}\equiv\hat{\sigma}_{l}^{ia}$.

The dissipative dynamics is described by Lindblad terms
$\mathcal{L}_{\alpha}\hat{\rho}=\sum_{l}\left[{\hat{L}}_{\alpha,l}\hat{\rho}{\hat{L}}_{\alpha,l}^{\dagger}-\frac{1}{2}\left\\{{\hat{L}}_{\alpha,l}^{\dagger}{\hat{L}}_{\alpha,l},\hat{\rho}\right\\}\right]$.
The spontaneous inactivation of the active states is described by
${\hat{L}}_{d,l}=\sqrt{\Gamma}\hat{\sigma}^{-}_{l}$, and
${\hat{L}}_{p,l}=\sqrt{\gamma_{de}}{\hat{\sigma}}_{l}^{aa}$ represents
dephasing of quantum coherence with rate $\gamma_{de}$. Meanwhile, the loss
and reloading of inactive emitters are accounted for by
${\hat{L}}_{e,l}=\sqrt{b\Gamma}\hat{\sigma}_{l}^{0a}$ and
${\hat{L}}_{a,l}=\sqrt{\lambda}\hat{\sigma}_{l}^{ip}$, respectively. The
incoherent cotact processes are also included in Lindbladian, where the
respective jump operators for activation and inactivation of emitters are
given by ${\hat{L}}_{b,l}=\sqrt{\kappa}{\hat{C}}_{l}{\hat{\sigma}}_{l}^{+}$
and ${\hat{L}}_{c,l}=\sqrt{\kappa}{\hat{C}}_{l}{\hat{\sigma}}_{l}^{-}$.

The Heisenberg-Langevin equations of motion for the operators
$\hat{\sigma}^{x/aa}_{l}$,
$\hat{\sigma}^{y}_{l}=i{\hat{\sigma}}_{l}^{-}-i{\hat{\sigma}}_{l}^{+}$, and
$\hat{n}_{l}=\hat{\sigma}_{l}^{aa}+\hat{\sigma}_{l}^{ii}$ according to the
master equation read

$\displaystyle\partial_{t}\hat{\sigma}_{l}^{aa}=$
$\displaystyle-\Gamma\hat{\sigma}_{l}^{aa}+\Omega{\hat{C}}_{l}\hat{\sigma}_{l}^{y}+\kappa\hat{C}_{l}\left(\hat{n}_{l}-2\hat{\sigma}_{l}^{aa}\right)+\hat{\xi}_{l}^{aa}$
(2a) $\displaystyle\partial_{t}\hat{\sigma}_{l}^{x}=$
$\displaystyle-\frac{\kappa\hat{N}_{l}+\gamma}{2}\hat{\sigma}^{x}_{l}-\kappa\hat{C}_{l}\hat{\sigma}_{l}^{x}-\Omega{\hat{P}}_{l}\hat{\sigma}_{l}^{y}+\hat{\xi}_{l}^{x}$
(2b) $\displaystyle\partial_{t}\hat{\sigma}_{l}^{y}=$
$\displaystyle-\frac{\kappa{\hat{N}}_{l}+\gamma}{2}\hat{\sigma}^{y}_{l}-\kappa{\hat{C}_{l}}\hat{\sigma}_{l}^{y}+\Omega{\hat{P}}_{l}\hat{\sigma}_{l}^{x}$
(2c)
$\displaystyle+2\Omega\hat{C}_{l}\left(\hat{n}_{l}-2\hat{\sigma}_{l}^{aa}\right)+\hat{\xi}_{l}^{y}$
$\displaystyle\partial_{t}\hat{n}_{l}=$
$\displaystyle-b\Gamma\hat{\sigma}_{l}^{aa}+\lambda\hat{\sigma}_{l}^{pp}+\hat{\xi}^{n}_{l}$
(2d)

where ${\hat{N}}_{l}=\sum_{k\in\partial l}{{\hat{n}}_{k}}$,
${\hat{P}}_{l}=\sum_{k\in\partial l}{{\hat{\sigma}}_{k}^{x}}$, and
$\gamma=\Gamma+\gamma_{de}$. The Langevin noise operators
$\hat{\xi}_{l}^{x/y/aa/n}$ appear because the dissipation is attributed to the
coupling between the system and a large reservoir [31], and can by fixed in
accord with the jump operators via solving the Heisenberg equations under the
reservoir Hamiltonians in Born-Markov approximation [See Section I in
Supplemental Material (SM) for details]. From now on we set the time unit to
$\Gamma^{-1}=1$.

In the following, we consider the continuum limit, and after coarse-graining
transform the expectation values of the operators into classical fields. More
specifically, the fields of active and total densities are defined as
$\rho\left(\bm{r},t\right)\equiv\left\langle\text{Tr}\left\\{\hat{\sigma}^{aa}_{l}\hat{\rho}\right\\}\right\rangle_{\bm{r}}$
and
$n\left(\bm{r},t\right)\equiv\left\langle\text{Tr}\left\\{\hat{n}_{l}\hat{\rho}\right\\}\right\rangle_{\bm{r}}$,
respectively, where $\langle...\rangle_{\bm{r}}$ denotes an average over the
facilitation sphere centered at $\bm{r}$, and the $\sigma^{x/y}$ fields can be
defined similarly. The fields $\rho,n,\sigma^{x/y}$ evolve according to the
Langevin equations that follow from omitting the operator moments generated by
the two-body interaction in Eqs. (2a)-(2d). In the presence of dephasing and
classical fluctuations, the fields $\sigma^{x/y}$ relax faster than $\rho,n$
and can be regarded as stochastic background fields for the latter. We then
perturbatively eliminate the fast variables $\sigma^{x/y}$ from the coupled
Langevin equations (2a)-(2c) via the Janssen-De Dominicis-Martin-Siggia-Rose
procedure (see SM)[32, 33]. The resulting action for active density field
reads

$S[\rho,\tilde{\rho}]=\int{\tilde{\rho}\left[\left(\partial_{t}-D_{\rho}\nabla^{2}+u_{2}\right)\rho+u_{3}\rho^{2}+u_{4}\rho^{3}-\frac{\mu}{2}\tilde{\rho}\right]}$
(3)

where $\tilde{\rho}$ is the Martin-Siggia-Rose auxiliary field related to the
dynamic responses of $\rho$ to perturbations, and
$u_{2}=1-n\kappa-256n^{2}\Omega^{4}/\left(n\kappa+\gamma\right)^{7},u_{3}=2\left[\kappa-2n\Omega^{2}/(n\kappa+\gamma)\right],u_{4}=8\Omega^{2}/(n\kappa+\gamma)$,
$\mu=\left(1+n\kappa\right)\rho+4n\Omega^{2}\rho^{2}/\left(n\kappa+\gamma\right)^{2}$
are the coupling constants. The diffusion constant $D_{\rho}=D_{T}+n\kappa
R_{\text{fac}}^{2}/2$, where $D_{T}$ is the thermal diffusivity.

With the total density $n$ conserved ($b,\lambda=0$), the MF static phases are
determined by the uniform, noiseless saddle-point solutions to the action (3).
The corresponding phase boundaries conditional on $\kappa,\Omega$ are shown in
Fig. 2 (a). In the quantum regime ($\Omega\gg\kappa$), the bistable region
within the two boundaries indicates that the systems undergo discontinuous
APTs when the total density $n$ exceeds a critical value (see SM for MF static
phase diagrams), which is an element of SOB-induced CTCs to be discussed
later. Approaching the classical regime ($\kappa\gg\Omega$), the bistable
region shrinks and finally vanishes when the transition becomes continuous,
which has been related to SOC in driven-dissipative Rydberg gases [34, 30].

Figure 2: MF results. (a) The two surfaces are the phase boundaries.
Discontinuous APTs ensue when the two surfaces are separated from each other.
(b) Stable (solid) and unstable (dotted) fixed point of Eqs. (4a) and (4b) as
a function of the loading rate $\lambda$. (c) Dynamical phase diagram where
the surface represents the critical loading rate $\lambda_{c}$ (color-coded)
that separates a STA phase from an SOB-induced CTC phase, corresponding to
stable and unstable fixed points, respectively. (d) Phase-space trajectories
(upper) and periods (lower) as a function of $\lambda$ (color-coded).
Parameters for (b) and (d) are $\Omega=0.5,\kappa=0$.

_SOB-induced CTCs._ —In the regime for first-order APTs, sustained
oscillations can arise from the interplay between loss and reloading of
emitters ($b,\lambda\neq 0$). Our proposal for SOB-induced CTCs are encoded in
the Langevin equations for the density fields $\rho,n$ as follows

$\displaystyle\partial_{t}\rho=D_{\rho}\nabla^{2}\rho+\tau n-u_{2}\rho-
u_{3}\rho^{2}-u_{4}\rho^{3}+\eta$ (4a)
$\displaystyle\partial_{t}n=D_{T}\nabla^{2}n-b\rho+\lambda+\xi^{n}$ (4b)

where $\eta,\xi^{n}$ are Markovian white noise sources with vanishing mean and
respective variance $\mu+\tau n$ and $b\rho$. Without limiting the generality,
we set the $p$-state density to $n_{p}=1$. To prevent the system from trapping
in absorbing states where the dynamics stops, a small driving $\tau n$ is
added [27], which can be realized by the spontaneous activation process [30].
Previous studies [17, 16] of CTCs in atomic gases focused on the
$D_{T,\rho}\to\infty$ limit, where spatial inhomogeneities are suppressed, and
MF theories are accurate. In our case, we work with finite diffusivity.
Throughout the paper, we fix $\tau=10^{-7},b=0.01,D_{T}=1,$ and $\gamma=2$,
unless otherwise stated.

MF dynamical phases can be determined through a linear stability analysis of
the fixed point possessed by Eqs. (4a) and (4b) [see Fig. 2 (b)]. In the
presence of first-order APTs, one can identify a critical loading rate
$\lambda_{c}$ [surface plotted in Fig. 2 (c)] below which LC phases that break
time translation symmetry arise from the Hopf bifurcation [35], similar to
other instability-driven LC phases [36, 37]. The existence of such oscillatory
phases requires first-order APTs [shaded region in Fig. 2 (c)], because for
systems with continuous APTs, the fixed points are always stable, and
stationary (STA) states are expected. Results of numerical integration of the
MF equations are displayed in Fig. 2 (d), where the SOB-induced LCs alternate
between two phases with a period $\approx\lambda^{-1}$.

Figure 3: The effective potential $\Phi_{f}$ as a function of the active
density field $\rho$ for different total densities $n$ in (a) $d=1$, (b)
$d=2$, (c) $d=3$. The active/absorbing phases are determined by the local
minima (indicated by dots). (d) The corresponding phase diagram in comparison
with the MF predictions. The result for $d=3$ approaches the MF line.
Parameters for (a)-(d) are $\kappa=0,\Omega=0.5$.  Figure 4: Time series of
the total and active densities for (a) upper row: $\lambda=8\times 10^{-4}$
with $d=1,L=10^{4}$ (left) and $d=3,L=64$ (right); lower row:
$\lambda=1.2\times 10^{-3}$ with $d=2,L=256$ (left) and $d=3,L=64$ (right);
and (b) three-dimensional systems with $L=24$ (left) and $L=64$ (right). (c)
Probability distribution of the average total densities. (d) The
autocorrelation functions (upper panel) and their Fourier spectra (lower
panel) of the total density fields. (e) Rescaled period defined as the inverse
of the location of the highest peak of the Fourier spectra. (f) The coherence
time $\tau_{\text{CTC}}\equiv\pi/\Delta\omega$, with the half-width
$\Delta\omega$ determined by Lorentz fit of the heighest spectral peak, and
(g) the occurrence probability of king avalanches for different system sizes.
(h) The coherence time as a function of system size and loading rate.
Parameters for (a)-(h) are $\kappa=0,\Omega=0.5$. Other parameters are $d=3$
[for (b)-(h)] and $\lambda=3.2\times 10^{-3}$ [for (b)-(g)].

_Dimensionality dependence of first-order APTs_ —Let us discuss first how
dimensionality affects the stability of the SOB-induced CTCs. Intuitively, the
lower critical dimension of the CTC dictates that of the first-order APTs,
which in low dimensionalities can be rounded by fluctuations, and consequently
rendering CTC unstable. To understand the dimensionality dependence of APTs,
we adopt a functional renormalization group (fRG) approach via the Wetterich
equation [38, 39, 40] to obtain the flows of the phase structure with
decreasing infrared cutoff for various values of $n$ in different dimensions
(see SM). The resulting effective potential $\Phi_{f}(\rho)$ and the
corresponding phase diagram are shown in Fig. 3. In $d=1$, as $n$ increases,
the position of the local minimum shifts continuously from the origin to a
finite value, indicating a continuous transition [Fig. 3 (a)]. In $d\geq 2$,
however, increasing $n$ induces the appearance of a second local minimum at
the finite density, apart from the local minimum at the origin, with a barrier
in between, indicating a first-order transition [Fig. 3 (b) and (c)]. Besides,
the barrier is higher in $d=3$ than $d=2$, suggesting a weaker first-order
transition in lower dimensions. The phase diagram in accord is shown in Fig. 3
(d), where we can infer that discontinuous transitions are expected for $d\geq
2$, and in comparison with the MF results, the coexistence region becomes
narrower for a lower dimensionality. The fRG results indicate that CTCs are
possible in $d\geq 2$, with the three-dimensional CTCs being protected by a
higher barrier between the two coexisting phases and thus more stable compared
with the two-dimensional ones.

We then numerically simulate Eqs. (4a) and (4b) deploying the operator-
spliting scheme [41, 42], and record the time-series of the average total and
active densities. As shown in Fig. 4 (a), deep in the CTC regime, non-
stationary phases breaking time-translation invariance are observed for
$d=2,3$. The wax-and-wane dynamics signals the underlying spatiotemporal
correlations induced by many-body interactions. Whereas compared with those in
three-dimensional systems [right column in Fig. 4 (a)], time series in two
dimensions [lower panel of the left column in Fig. 4 (a)] are less regular,
implying the destruction of time crystalline order by fluctuations, in accord
with the much weaker first-order APTs therein, as revealed by the fRG
analysis. Given the enhanced stability of CTCs in higher dimensionalities, in
the following, we restrict our discussion to $d=3$. Results for $d=1,2$ can be
found in SM.

_Finite-size effects.—_ We can infer from Fig. 4 (b) that more ordered time
series and hence stable CTCs are expected in larger systems. Meanwhile, the
distribution of average total density remains broad as the system is enlarged
[Fig. 4 (c)], reflecting the neverending transitions from one phase to
another. A more precise characterization of the time crystalline order is
provided by the autocorrelation functions of the average density fields,
defined as $G_{x}(\Delta t)=\left\langle
x(t)\right\rangle_{t}^{-2}\left\langle x(t)x(t+\Delta t)\right\rangle_{t}$,
for $x=n,\rho$ [16]. The two-time correlation functions manifest constant
periodic oscillations for perfect time crystals and otherwise exhibit
modulations at times due to decoherence. Correspondingly, their Fourier
spectra $G_{x}(\omega)$ peak at the integer multiples of their respective
inherent frequencies $\omega_{m}$. Results for systems with various sizes are
displayed in Fig. 4 (d), where it is clear that, as the edge length $L$
increases, the amplitude of $G_{n}(\Delta t)$ varies more slowly, and the
Fourier spectra are more sharply peaked at $\omega/\omega_{m}=1,2,3…$, typical
of periodic structures in time. Besides, once the time crystalline order is
built, the period $T\equiv 2\pi/\omega_{m}$ remains invariant with diverging
$L$ [Fig. 4 (e)], and is thus inherent to CTCs. To quantify time crystalline
order, we follow [43] to estimate the coherence time
$\tau_{\text{CTC}}\equiv\pi/\Delta\omega$, where the half-width $\Delta\omega$
is determined by a Lorentz fitting of the Fourier spectra. As we can infer
from Fig. 4 (f), albeit finite, the coherence time increases with system size,
suggesting the occurrence of persistent oscillations in the $L\to\infty$
limit. The existence of sustained oscillations with an intrinsic amplitude and
frequency suggests that our CTCs are also a realization of BTCs [10].

Our simulations reveal the finite-size effect, which can be readily explained
through phenomenological treatments like introducing to the MF equations a
noise source with covariance inversely proportional to the system volume [44,
45, 46]. The irregularity in the time series typically consists of a
significant increase in the total density followed by an abrupt decrease in
the active and total densities [see Fig. 4 (b), left panel]. For low loading
rates, such events have been attributed to the system falling into the
absorbing state, and the consequent overloading in turn brings about system-
spanning activation avalanches triggered by one active seed [47, 48]. We argue
that the same reasoning applies to the aperiodic oscillations in small systems
therein.

To test this idea, we count space-time activation avalanches by connecting
sites with active densities larger than a threshold ($\tau$) as neighbors in
the time-forward direction and grouping them into clusters [49]. The size of
the activation avalanches is defined to be the number of sites constituting
them. If two clusters collide, the collision points (interfaces) are chosen to
be the boundary between them. The occurrence probability of huge avalanches
(king avalanches, defined as those larger than half the system volume)
decreases as the system becomes larger [Fig. 4 (g)], in line with the longer
coherence time [Fig. 4 (f)]. A comparison of $\tau_{\text{CTC}}$ among various
loading rates and system sizes is displayed in Fig. 4 (h), the regime for CTCs
indicated by the significantly increased coherence time lies between that for
the aperiodic oscillations and the fluctuating uniform ones, and widens in
larger systems.

In small systems and for low loading rates, frequent huge avalanches lead to
coherent changes in active and total densities among a great many sites and
thus reflect the underlying quantum synchrony at its highest level. However, a
lack of synchrony results in stationary states with small fluctuations that
conserve time-translation invariance, as is the case for larger values of
$\lambda$. Sustained periodic oscillations reside in between the above two
scenarios, where discontinuous phase transitions spontaneously generate
finite-range correlations, which are enough to trigger coherence among local
sites and yet unable to support a global synchronization in infinite systems.
In other words, SOB-induced CTCs arise at the edge of quantum synchronization.

Our model can be implemented with coherent laser-driven Rydberg atoms in the
anti-blockade regime, where the electronic ground (Rydberg) states can be
mapped to the inactive (active) states, cf. [29, 34, 30]. The relative
importance of the coherent and incoherent activation processes is controlled
by dephasing. Because we have measured all frequencies with the decay rate and
the volume with that of the facilitation sphere, the assumption that the
density fields vary smoothly in space holds in the finite diffusivity regime
($D_{T}\geq n^{-2/d}\Gamma$) such that the Griffiths effects [58, 59] reported
in frozen gases ($D_{T}\approx 0$) are negligible, otherwise, addtional
geometric factor related to the thickness of the facilitation shell should be
included in $\kappa,\Omega$. Different dimensionalities can be realized by
changing the geometry of the optical trap.

_Conclusion and discussion._ —In this work, we propose a mechanism to realize
robust CTCs in far-from-equilibrium systems without either all-to-all coupling
or infinite thermal diffusivity. Our analysis is not restricted to APTs, and
also related to other systems with bi-/multi-stability, where collective
oscillatory phases can likewise be constructed. Our CTCs can be interpreted as
a BTC [10], where the reservoir plays the role analogous to the bulk
Hamiltonian therein, such that tracing out the corresponding degrees of
freedom yields the master equation (Lindbladian) for the (boundary) degrees of
freedom, where sustained oscillations with intrinsic amplitudes and periods
emerge in the thermodynamic limit.

Revealing how temporal organization arises from the SOB-induced bifurcation,
our study extends the dynamical phase diagram for both SOB and instability-
related LCs [50, 36, 27]. Moreover, the existence of CTCs in noisy,
diffusively coupled systems further enhances the potential of strongly
interacting, open quantum system to serve as controllable platforms for
studying self-organization in far-from-equilibrium settings. Furthermore, the
coarse-graining procedures to derive effective dynamics from microscopic ones
provide us with insight into the relation of other effective stochastic
equations to the microscopic underpinnings of real-world critical-like events
[51, 26, 52, 53, 47, 54], alongside the reaction-diffusion systems, where the
equations can be derived rigorously through the coherent-state path integral
formalism [55, 56, 57].

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
under Grants No. 11974175, No. 12247102 and No. 12274131. We are grateful to
the High Performance Computing Center (HPCC) of Nanjing University for
performing the numerical calculations in this paper on its blade cluster
system.

## References

  * Wilczek [2012] F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012).
  * Shapere and Wilczek [2012] A. Shapere and F. Wilczek, Phys. Rev. Lett. 109, 160402 (2012).
  * Choi  _et al._ [2017] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Nature 543, 221 (2017).
  * Smits  _et al._ [2018] J. Smits, L. Liao, H. T. C. Stoof, and P. van der Straten, Phys. Rev. Lett. 121, 185301 (2018).
  * Gambetta  _et al._ [2019] F. M. Gambetta, F. Carollo, M. Marcuzzi, J. P. Garrahan, and I. Lesanovsky, Phys. Rev. Lett. 122, 015701 (2019).
  * Keßler  _et al._ [2021] H. Keßler, P. Kongkhambut, C. Georges, L. Mathey, J. G. Cosme, and A. Hemmerich, Phys. Rev. Lett. 127, 043602 (2021).
  * Bluvstein  _et al._ [2021] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, S. Choi, M. Serbyn, M. Greiner, V. Vuletić, and M. D. Lukin, Science 371, 1355 (2021).
  * Piazza and Ritsch [2015] F. Piazza and H. Ritsch, Phys. Rev. Lett. 115, 163601 (2015).
  * Chan  _et al._ [2015] C.-K. Chan, T. E. Lee, and S. Gopalakrishnan, Phys. Rev. A 91, 051601 (2015).
  * Iemini  _et al._ [2018] F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte, and R. Fazio, Phys. Rev. Lett. 121, 035301 (2018).
  * Buča  _et al._ [2019] B. Buča, J. Tindall, and D. Jaksch, Nat. Commun. 10, 1 (2019).
  * Bruno [2013] P. Bruno, Phys. Rev. Lett. 111, 070402 (2013).
  * Watanabe and Oshikawa [2015] H. Watanabe and M. Oshikawa, Phys. Rev. Lett. 114, 251603 (2015).
  * Sheng  _et al._ [2020] J. Sheng, X. Wei, C. Yang, and H. Wu, Phys. Rev. Lett. 124, 053604 (2020).
  * Ding  _et al._ [2023] D.-S. Ding, Z. Bai, Z.-K. Liu, B.-S. Shi, G.-C. Guo, W. Li, and C. S. Adams, arXiv 2302.14145 (2023).
  * Wu  _et al._ [2023] X. Wu, Z. Wang, F. Yang, R. Gao, C. Liang, M. K. Tey, X. Li, T. Pohl, and L. You, arXiv:2305.20070 (2023).
  * Wadenpfuhl and Adams [2023] K. Wadenpfuhl and C. S. Adams, Phys. Rev. Lett. 131, 143002 (2023).
  * Keßler  _et al._ [2019] H. Keßler, J. G. Cosme, M. Hemmerling, L. Mathey, and A. Hemmerich, Phys. Rev. A 99, 053605 (2019).
  * Kongkhambut  _et al._ [2022] P. Kongkhambut, J. Skulte, L. Mathey, J. G. Cosme, A. Hemmerich, and H. Keßler, Science 377, 670 (2022).
  * Tang and Bak [1988] C. Tang and P. Bak, J. Stat. Phys. 51, 797 (1988).
  * Sornette [1992] D. Sornette, J. Phys. I 2, 2065 (1992).
  * Paczuski  _et al._ [1994] M. Paczuski, S. Maslov, and P. Bak, EPL 27, 97 (1994).
  * Grinstein [1995] G. Grinstein, Generic scale invariance and self-organized criticality, in  _Scale Invariance, Interfaces, and Non-Equilibrium Dynamics_, edited by A. McKane, M. Droz, J. Vannimenus, and D. Wolf (Springer US, Boston, MA, 1995) pp. 261–293.
  * Gil and Sornette [1996] L. Gil and D. Sornette, Phys. Rev. Lett. 76, 3991 (1996).
  * Dickman  _et al._ [1998] R. Dickman, A. Vespignani, and S. Zapperi, Phys. Rev. E 57, 5095 (1998).
  * Dickman  _et al._ [2000] R. Dickman, M. A. Muñoz, A. Vespignani, and S. Zapperi, Braz. J. Phys. 30, 27 (2000).
  * Buendía  _et al._ [2020] V. Buendía, S. di Santo, J. A. Bonachela, and M. A. Muñoz, Front. Phys. 8, 333 (2020).
  * di Santo  _et al._ [2016] S. di Santo, R. Burioni, A. Vezzani, and M. A. Muñoz, Phys. Rev. Lett. 116, 240601 (2016).
  * Marcuzzi  _et al._ [2016] M. Marcuzzi, M. Buchhold, S. Diehl, and I. Lesanovsky, Phys. Rev. Lett. 116, 245701 (2016).
  * Helmrich  _et al._ [2020] S. Helmrich, A. Arias, G. Lochead, T. M. Wintermantel, M. Buchhold, S. Diehl, and S. Whitlock, Nature 577, 481 (2020).
  * Pan  _et al._ [2020] L. Pan, X. Chen, Y. Chen, and H. Zhai, Nat. Phys. 16, 767 (2020).
  * Janssen [1976] H.-K. Janssen, Z. Physik B 23, 377 (1976).
  * Martin  _et al._ [1973] P. C. Martin, E. D. Siggia, and H. A. Rose, Phys. Rev. A 8, 423 (1973).
  * Klocke and Buchhold [2019] K. Klocke and M. Buchhold, Phys. Rev. A 99, 053616 (2019).
  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993).
  * Lefever and Nicolis [1971] R. Lefever and G. Nicolis, J. theor. Biol. 30, 267 (1971).
  * Lee  _et al._ [2011] T. E. Lee, H. Häffner, and M. C. Cross, Phys. Rev. A 84, 031402 (2011).
  * Wetterich [1993] C. Wetterich, Phys. Lett. B 301, 90 (1993).
  * Canet  _et al._ [2011] L. Canet, H. Chaté, and B. Delamotte, J. Phys. A: Math. Theor. 44, 495001 (2011).
  * Dupuis  _et al._ [2021] N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J. Pawlowski, M. Tissier, and N. Wschebor, Phys. Rep. 910, 1 (2021).
  * Pechenik and Levine [1999] L. Pechenik and H. Levine, Phys. Rev. E 59, 3893 (1999).
  * Dornic  _et al._ [2005] I. Dornic, H. Chaté, and M. A. Muñoz, Phys. Rev. Lett. 94, 100601 (2005).
  * Chen and Zhang [2023] Y. H. Chen and X. Zhang, Nat. Commun. 14, 6161 (2023).
  * Boland  _et al._ [2008] R. P. Boland, T. Galla, and A. J. McKane, J. Stat. Mech. 2008, P09001 (2008).
  * Cabot  _et al._ [2023] A. Cabot, L. S. Muhle, F. Carollo, and I. Lesanovsky, Phys. Rev. A 108, L041303 (2023).
  * Kampen [2007] N. V. Kampen,  _Stochastic Processes in Physics and Chemistry_ , 3rd ed., North-Holland personal library (Elsevier, Singapore, 2007).
  * Grassberger [2002] P. Grassberger, New J. Phys. 4, 17 (2002).
  * Kinouchi  _et al._ [2019] O. Kinouchi, L. Brochini, A. A. Costa, J. G. F. Campos, and M. Copelli, Sci. Rep. 9, 3874 (2019).
  * Ding  _et al._ [2020] D.-S. Ding, H. Busche, B.-S. Shi, G.-C. Guo, and C. S. Adams, Phys. Rev. X 10, 021023 (2020).
  * Epstein and Showalter [1996] I. R. Epstein and K. Showalter, J. Phys. Chem. 100, 13132 (1996).
  * Malamud  _et al._ [1998] B. D. Malamud, G. Morein, and D. L. Turcotte, Science 281, 1840 (1998).
  * Yoshioka [2003] N. Yoshioka, Earth Planet. Sp. 55, 283 (2003).
  * Beggs and Plenz [2003] J. M. Beggs and D. Plenz, J Neurosci 23, 11167 (2003).
  * Bonachela and Muñoz [2009] J. A. Bonachela and M. A. Muñoz, J. Stat. Mech. 2009, P09009 (2009).
  * Doi [1976] M. Doi, J. Phys. A: Math. Gen. 9, 1479 (1976).
  * Peliti [1985] L. Peliti, J. Physique 46, 1469 (1985).
  * Wiese [2016] K. J. Wiese, Phys. Rev. E 93, 042117 (2016).
  * Wintermantel  _et al._ [2021] T. M. Wintermantel, M. Buchhold, S. Shevate, M. Morgado, Y. Wang, G. Lochead, S. Diehl, and S. Whitlock, Nat. Commun. 12, 103 (2021).
  * Brady  _et al._ [2023] D. Brady, J. Bender, P. Mischke, T. Niederprüm, H. Ott, and M. Fleischhauer, arXiv 2302.14145 (2023).
  * Buchhold  _et al._ [2017] M. Buchhold, B. Everest, M. Marcuzzi, I. Lesanovsky, and S. Diehl, Phys. Rev. B 95, 014308 (2017).
  * Wehner and Wolfer [1983] M. F. Wehner and W. G. Wolfer, Phys. Rev. A 27, 2663 (1983).
  * Täuber [2014] U. C. Täuber,  _Critical Dynamics: A Field Theory Approach to Equilibrium and Non-Equilibrium Scaling Behavior_ (Cambridge University Press, New York, USA, 2014).

Supplemental Material for

Self-Organized Time Crystal in Driven-Dissipative Quantum System

Ya-Xin Xiang, Qun-Li Lei, Zhengyang Bai and Yu-Qiang Ma

This supplemental material contains details on the derivation of the noise
operators (Sec. I), the the coarse-graining procedure (Sec. II), the Janssen-
De Dominicis-Martin-Siggia-Rose procedure for the construction of dynamic path
integrals (Sec. III), the functional renormalization group analysis of the
first-order APTs (Sec. IV), the numerical integration scheme for
multiplicative noise (Sec. V), and simulation results in one and two
dimensions (Sec. VI).

## I. The microscopic model

Figure S5: The facilitation activation occurs within the facilitation shell
(with a radius of $R_{\text{fac}}$ and thickness of $r_{\text{fac}}$) of
active emitters.

###  .1 Effective nearest neighbors

To begin with, we introduce the notion of the effective nearest neighbor. As
illustrated in Fig. S5, an active emitter facilitates the activation of
inactive emitters if the relative distance between them lies within the range
$\left[R_{\rm fac}-r_{\rm fac},R_{\rm fac}+r_{\rm fac}\right]$ (cf. [30]) with
a radius of $R_{\text{fac}}$ and thickness of $r_{\text{fac}}$. The inactive
emitter is in this sense considered to be the effective nearest neighbor
(referred to as “nn” in what follows). To count the number of “nn” in free
space, we introduce for later convenience the function
$f\left(a;R,r\right)\equiv\theta\left(a-R-r\right)-\theta\left(a-R+r\right)$,
with $\theta\left(x\right)$ being the Heaviside step function, one can write

$\displaystyle{\hat{C}}_{l}=\sum_{k}{\hat{\sigma}}_{k}^{aa}f\left(\left|\bm{r}_{lk}\right|;R_{\text{fac}},r_{\text{fac}}\right)$
(S5a)
$\displaystyle{\hat{N}}_{l}=\sum_{k}{\hat{n}}_{k}f\left(\left|\bm{r}_{lk}\right|;R_{\text{fac}},r_{\text{fac}}\right)$
(S5b)
$\displaystyle{\hat{P}}_{l}=\sum_{k}{\hat{\sigma}}_{k}^{x}f\left(\left|\bm{r}_{lk}\right|;R_{\text{fac}},r_{\text{fac}}\right)$
(S5c)

where $\hat{\sigma}_{l}^{bc}\equiv|b_{l}\rangle\langle c_{l}|$, $\bm{r}_{lk}$
is the relative distance between the $k$-th and $l$-th emitters, and
$\hat{\sigma}_{l}^{x}=\hat{\sigma}^{+}_{l}+\hat{\sigma}^{-}_{l}$ with the
ladder operators defined as
$\hat{\sigma}^{+}_{l}\equiv\sigma_{l}^{ai},\hat{\sigma}^{-}_{l}\equiv\hat{\sigma}_{l}^{ia}$

###  .2 Reservoir Hamiltonian

We follow [60, 29] to write the three local reservoir Hamiltonian, where
$\hat{H}_{d}$ represents the spontaneous decay, and $\hat{H}_{b}$ and
$\hat{H}_{c}$ account for the incoherent contact activation and the inverse
process, respectively.

$\displaystyle\hat{H}_{d}=$
$\displaystyle\sum_{q}\lambda_{q}^{d}\left(\hat{\sigma}^{+}\hat{d}_{q}+\hat{d}_{q}^{\dagger}\hat{\sigma}^{-}\right)+\sum_{q}\omega_{q}^{d}\hat{d}_{q}^{\dagger}\hat{d}_{q}$
(S6a) $\displaystyle\hat{H}_{b}=$
$\displaystyle\sum_{q,\text{nn}}\lambda_{q}^{b}\left(\hat{\sigma}^{aa}_{\text{nn}}\hat{b}_{q}^{\dagger}\hat{\sigma}^{+}+\hat{\sigma}^{aa}\hat{b}_{q,\text{nn}}^{\dagger}\hat{\sigma}^{+}_{\text{nn}}+\text{H.c.}\right)$
(S6b)
$\displaystyle+\sum_{q}\omega_{q}^{c}\left(\hat{b}_{q}^{\dagger}\hat{b}_{q}+\sum_{\text{nn}}\hat{b}_{q,\text{nn}}^{\dagger}\hat{b}_{q,\text{nn}}\right)$
$\displaystyle\hat{H}_{c}=$
$\displaystyle\sum_{q,\text{nn}}\lambda_{q}^{b}\left(\hat{\sigma}^{aa}_{\text{nn}}\hat{\sigma}^{+}\hat{c}_{q}+\hat{\sigma}^{aa}\hat{\sigma}^{+}_{\text{nn}}\hat{c}_{q,\text{nn}}+\text{H.c.}\right)$
(S6c)
$\displaystyle+\sum_{q}\omega_{q}^{c}\left(\hat{c}_{q}^{\dagger}\hat{c}_{q}+\sum_{\text{nn}}\hat{c}_{q,\text{nn}}^{\dagger}\hat{c}_{q,\text{nn}}\right)$

From now on, the indices for emitters within these Hamiltonians are dropped
because the noise is Markovian, and the operators
$\hat{d}_{q},\hat{b}_{q},\hat{c}_{q}$ are bosonic bath modes with their
respective frequencies $\omega_{q}^{d/b/c}$, and the coupling strength
$\lambda_{q}^{d/b/c}$. Besides, we discard the possible noise correlations
between a pair of neighboring sites, which arises from the Hamiltonians
$\hat{H}_{b/c}$ acting on the neighbors of a given site and vice versa [60].
The noise operators in accord with the jump operators in the main text are
fixed via solving the Heisenberg equations under the above Hamiltonians in
Born-Markov approximation.

###  .3 Spontaneous decay

For starter, let us consider the spontaneous decay (inactivation). The
Heisenberg equations under the action of $\hat{H}_{d}$ read

$\displaystyle\partial_{t}\hat{\sigma}^{aa}$
$\displaystyle=i\left[\hat{H}_{d},\hat{\sigma}^{aa}\right]=i\sum_{q}\lambda_{q}^{d}\left(\hat{d}_{q}^{\dagger}\hat{\sigma}^{-}-\hat{\sigma}^{+}\hat{d}_{q}\right)$
(S7a) $\displaystyle\partial_{t}\hat{\sigma}^{+}$
$\displaystyle=i\left[\hat{H}_{d},\hat{\sigma}^{+}\right]=-i\sum_{q}\lambda_{q}^{d}\hat{d}_{q}^{\dagger}\hat{\sigma}^{z}$
(S7b) $\displaystyle\partial_{t}\hat{d}_{q}$
$\displaystyle=i\left[\hat{H}_{d},\hat{d}_{q}\right]=-i\lambda_{q}^{d}\hat{\sigma}^{-}-i\omega_{q}^{d}\hat{d}_{q}$
(S7c)

Formally integration of the last equation (S7c) leads to

$\hat{d}_{q}(t)=\hat{d}_{q}(0)e^{-i\omega_{q}^{d}t}-i\lambda^{d}_{q}\int_{0}^{t}dt^{\prime}\hat{\sigma}^{-}(t^{\prime})e^{-i\omega_{q}^{d}(t-t^{\prime})}$
(S8)

Then the rest two equations become

$\displaystyle\partial_{t}\hat{\sigma}^{aa}=$ $\displaystyle
i\sum_{q}\lambda_{q}^{d}\left[\hat{d}_{q}^{\dagger}(0)\hat{\sigma}^{-}(t)e^{i\omega_{q}^{d}t}-\text{H.c.}\right]$
(S9a)
$\displaystyle-\sum_{q}(\lambda_{q}^{d})^{2}\int_{0}^{t}dt^{\prime}\left[\hat{\sigma}^{+}\left(t^{\prime}\right)\hat{\sigma}^{-}(t)e^{i\omega_{q}^{d}(t-t^{\prime})}+\text{H.c.}\right]$
$\displaystyle\partial_{t}\hat{\sigma}^{+}=$
$\displaystyle-i\sum_{q}\lambda_{q}^{d}\hat{d}_{q}^{\dagger}(0)\hat{\sigma}^{z}(t)e^{i\omega_{q}^{d}t}$
(S9b)
$\displaystyle+\sum_{q}(\lambda_{q}^{d})^{2}\int_{0}^{t}dt^{\prime}\hat{\sigma}^{+}(t^{\prime})\hat{\sigma}^{z}(t)e^{i\omega_{q}^{d}(t-t^{\prime})}$

The first addends on the right hand side of Eqs. (S9a) and (S9b) are the noise
operators corresponding to the spontaneous decay process,

$\displaystyle\hat{\xi}^{aa}_{d}(t)$
$\displaystyle=i\sum_{q}\lambda_{q}^{d}\left[\hat{d}_{q}^{\dagger}(0)\hat{\sigma}^{-}(t)e^{i\omega_{q}^{d}t}-\text{H.c.}\right]$
(S10a) $\displaystyle\hat{\xi}^{+}(t)$
$\displaystyle=-i\sum_{q}\lambda_{q}^{d}\hat{d}_{q}^{\dagger}(0)\hat{\sigma}^{z}(t)e^{i\omega_{q}^{d}t}$
(S10b)

and the second addends represent the dissipation.

In the Born-Markov approximation, the summation over all frequencies
contributes to a delta function $\delta(t-t^{\prime})$, then the dissipation
rate are related to the effective coupling constant
$\lambda^{d}(0)=\sum_{q}\delta(\omega_{q}^{d})\lambda_{q}^{d}$ and the bath
density of states at zero frequency $D(0)=\sum_{q}\delta(\omega_{q}^{d})$, via
$2\pi D(0)[\lambda^{d}(0)]^{2}=1$. It follows from
$\langle...\rangle=\text{Tr}\left\\{(...)\hat{\rho}^{0}_{d}\right\\}$ and
$\langle\hat{d}_{q}^{\dagger}\hat{d}_{q}\rangle=n_{q}=0$, that
$\langle\hat{\xi}^{aa/+}_{d}\rangle=0$, and the non-vanishing covariance reads

$\displaystyle\langle\hat{\xi}^{aa}_{d}(t)\hat{\xi}^{aa}_{d}(t^{\prime})\rangle=\hat{\sigma}^{aa}\delta(t-t^{\prime})$
(S11a)
$\displaystyle\langle\hat{\xi}^{aa}_{d}(t)\hat{\xi}^{+}_{d}(t^{\prime})\rangle=\hat{\sigma}^{+}\delta(t-t^{\prime})$
(S11b)
$\displaystyle\langle\hat{\xi}^{-}_{d}(t)\hat{\xi}^{+}_{d}(t^{\prime})\rangle=\hat{n}\delta(t-t^{\prime})$
(S11c)

Together with
$\hat{\xi}_{d}^{x}=\hat{\xi}_{d}^{+}+\hat{\xi}_{d}^{-},\hat{\xi}_{d}^{y}=-i\hat{\xi}_{d}^{+}+i\hat{\xi}_{d}^{-}$,
one arrives at

$\displaystyle\langle\hat{\xi}^{aa}_{d}(t)\hat{\xi}^{x}_{d}(t^{\prime})\rangle=\hat{\sigma}^{+}\delta(t-t^{\prime})$
(S12a)
$\displaystyle\langle\hat{\xi}^{aa}_{d}(t)\hat{\xi}^{y}_{d}(t^{\prime})\rangle=-i\hat{\sigma}^{+}\delta(t-t^{\prime})$
(S12b)
$\displaystyle\langle\hat{\xi}^{x}_{d}(t)\hat{\xi}^{y}_{d}(t^{\prime})\rangle=-i\hat{n}\delta(t-t^{\prime})$
(S12c)
$\displaystyle\langle\hat{\xi}^{x}_{d}(t)\hat{\xi}^{x}_{d}(t^{\prime})\rangle=\langle\hat{\xi}^{y}_{d}(t)\hat{\xi}^{y}_{d}(t^{\prime})\rangle=\hat{n}\delta(t-t^{\prime})$
(S12d)

###  .4 Incoherent contact processes

Hereafter we introduce
$\lambda_{q}=\lambda_{q}^{b/c},\omega_{q}=\omega_{q}^{b/c}$, since the
frequencies and the coupling of the two baths are the same. In the same spirit
as the Sec. IC, the Hamitonian $\hat{H}_{b}$ yields the following equations

$\displaystyle\partial_{t}\hat{\sigma}^{+}=-i\sum_{q,\text{nn}}\lambda_{q}\left[\hat{\sigma}_{\text{nn}}^{aa}\hat{\sigma}^{z}\hat{b}_{q}-\left(\hat{\sigma}^{+}\hat{\sigma}_{\text{nn}}^{-}\hat{b}_{q,\text{nn}}+\text{H.c.}\right)\right]$
(S13a)
$\displaystyle\partial_{t}\hat{\sigma}^{aa}=i\sum_{q,\text{nn}}\lambda_{q}\left(\hat{\sigma}^{aa}_{\text{nn}}\hat{\sigma}^{-}\hat{b}_{q}-\text{H.c.}\right)$
(S13b)
$\displaystyle\partial_{t}\hat{b}_{q}=-i\sum_{\text{nn}}\lambda_{q}\hat{\sigma}_{\text{nn}}^{aa}\hat{\sigma}^{+}-i\omega_{q}\hat{b}_{q}$
(S13c)
$\displaystyle\partial_{t}\hat{b}_{q,\text{nn}}=-i\lambda_{q}\hat{\sigma}^{aa}\hat{\sigma}^{+}_{\text{nn}}-i\omega_{q}\hat{b}_{q,\text{nn}}$
(S13d)

The equations under the Hamiltonian $\hat{H}_{c}$ read

$\displaystyle\partial_{t}\hat{\sigma}^{+}=-i\sum_{q,\text{nn}}\lambda_{q}\left[\hat{\sigma}_{\text{nn}}^{aa}\hat{c}_{q}^{\dagger}\hat{\sigma}^{z}-\left(\hat{\sigma}^{+}\hat{\sigma}_{\text{nn}}^{+}\hat{c}_{q,\text{nn}}+\text{H.c.}\right)\right]$
(S14a)
$\displaystyle\partial_{t}\hat{\sigma}^{aa}=i\sum_{q,\text{nn}}\lambda_{q}\left(\hat{\sigma}_{\text{nn}}^{aa}\hat{c}_{q}^{\dagger}\hat{\sigma}^{-}-\text{H.c.}\right)$
(S14b)
$\displaystyle\partial_{t}\hat{c}_{q}=-i\sum_{\text{nn}}\lambda_{q}\hat{\sigma}_{\text{nn}}^{aa}\hat{\sigma}^{-}-i\omega_{q}\hat{c}_{q}$
(S14c)
$\displaystyle\partial_{t}\hat{c}_{q,\text{nn}}=-i\lambda_{q}\hat{\sigma}^{aa}\hat{\sigma}_{\text{nn}}^{-}-i\omega_{q}\hat{c}_{q,\text{nn}}$
(S14d)

Similar to the procedure in Sec. IC, by solving the Eqs. (S13) and (S14), the
spin operators $\hat{\sigma}^{aa/+}$ yielding the following noise operators,

$\displaystyle\hat{\xi}^{+}_{b}(t)$
$\displaystyle=-i\sum_{q}\lambda_{q}\hat{C}(t)\hat{\sigma}^{z}(t)\hat{b}_{q}(0)e^{-i\omega_{q}t}$
(S15a)
$\displaystyle+i\sum_{q,\text{nn}}\lambda_{q}\left[\hat{\sigma}^{+}(t)\hat{\sigma}_{\text{nn}}^{-}(t)\hat{b}_{q,\text{nn}}(0)e^{-i\omega_{q}t}+\text{H.c.}\right]$
$\displaystyle\hat{\xi}^{aa}_{b}(t)$
$\displaystyle=i\sum_{q}\lambda_{q}\hat{C}(t)\left[\hat{\sigma}^{-}(t)\hat{b}_{q}(0)e^{-i\omega_{q}t}-\text{H.c.}\right]$
(S15b) $\displaystyle\hat{\xi}^{+}_{c}(t)$
$\displaystyle=-i\sum_{q}\lambda_{q}\hat{C}(t)\hat{c}_{q}^{\dagger}(0)\hat{\sigma}^{z}(t)e^{i\omega_{q}t}$
(S15c)
$\displaystyle+i\sum_{q,\text{nn}}\lambda_{q}\left[\hat{\sigma}^{+}(t)\hat{\sigma}^{+}_{\text{nn}}(t)\hat{c}_{q,\text{nn}}(0)e^{-i\omega_{q}t}+\text{H.c.}\right]$
$\displaystyle\hat{\xi}^{aa}_{c}(t)$
$\displaystyle=i\sum_{q}\lambda_{q}\hat{C}(t)\left[\hat{c}_{q}^{\dagger}(0)\hat{\sigma}^{-}(t)e^{i\omega_{q}t}-\text{H.c.}\right]$
(S15d)

One can compute the noise covariance in the Born-Markov approximation $2\pi
D(0)[\lambda^{b/c}(0)]^{2}=\kappa/2$, of which the non-vanishing ones read

$\displaystyle\langle\hat{\xi}_{b}^{+}(t)\hat{\xi}_{b}^{-}(t^{\prime})\rangle=$
$\displaystyle\frac{\kappa}{2}\left[\hat{C}(t)\hat{\sigma}^{ii}(t)+\hat{N}(t)\hat{\sigma}^{aa}(t)\right]\delta(t-t^{\prime})$
(S16a)
$\displaystyle\langle\hat{\xi}_{b}^{-}(t)\hat{\xi}_{b}^{+}(t^{\prime})\rangle=$
$\displaystyle\frac{\kappa}{2}\left[\hat{N}(t)-\hat{C}(t)\right]\hat{\sigma}^{aa}(t)\delta(t-t^{\prime})$
(S16b)
$\displaystyle\langle\hat{\xi}_{b}^{+}(t)\hat{\xi}_{b}^{aa}(t^{\prime})\rangle=$
$\displaystyle-\frac{\kappa}{2}\hat{C}(t)\hat{\sigma}^{+}(t)\delta(t-t^{\prime})$
(S16c)
$\displaystyle\langle\hat{\xi}_{b}^{aa}(t)\hat{\xi}_{b}^{-}(t^{\prime})\rangle=$
$\displaystyle-\frac{\kappa}{2}\hat{C}(t)\hat{\sigma}^{-}(t)\delta(t-t^{\prime})$
(S16d)
$\displaystyle\langle\hat{\xi}_{b}^{aa}(t)\hat{\xi}_{b}^{aa}(t^{\prime})\rangle=$
$\displaystyle\frac{\kappa}{2}\hat{C}(t)\hat{\sigma}^{ii}(t)\delta(t-t^{\prime})$
(S16e)

for facilitated activation, and

$\displaystyle\langle\hat{\xi}_{c}^{+}(t)\hat{\xi}_{c}^{-}(t^{\prime})\rangle=$
$\displaystyle\frac{\kappa}{2}\hat{C}(t)\hat{\sigma}^{aa}(t)\delta(t-t^{\prime})$
(S17a)
$\displaystyle\langle\hat{\xi}_{c}^{-}(t)\hat{\xi}_{c}^{+}(t^{\prime})\rangle=$
$\displaystyle\frac{\kappa}{2}\hat{C}(t)\left[\hat{\sigma}^{aa}(t)+\hat{n}(t)\right]\delta(t-t^{\prime})$
$\displaystyle\langle\hat{\xi}_{c}^{-}(t)\hat{\xi}_{c}^{aa}(t^{\prime})\rangle=$
$\displaystyle\frac{\kappa}{2}\hat{C}(t)\hat{\sigma}^{-}(t)\delta(t-t^{\prime})$
(S17b)
$\displaystyle\langle\hat{\xi}_{c}^{aa}(t)\hat{\xi}_{c}^{+}(t^{\prime})\rangle=$
$\displaystyle\frac{\kappa}{2}\hat{C}(t)\hat{\sigma}^{+}(t)\delta(t-t^{\prime})$
(S17c)
$\displaystyle\langle\hat{\xi}_{c}^{aa}(t)\hat{\xi}_{c}^{aa}(t^{\prime})\rangle=$
$\displaystyle\frac{\kappa}{2}\hat{C}(t)\hat{\sigma}^{aa}(t)\delta(t-t^{\prime})$
(S17d)

for the inverse process.

Then the covariance of the noise operators $\hat{\xi}^{x/y/aa}_{f}$ arising
from the above two processes can be obtained through the combinations
according to

$\displaystyle\hat{\xi}^{+/aa}_{f}=\hat{\xi}^{+/aa}_{b}+\hat{\xi}^{+/aa}_{c}$
(S18a) $\displaystyle\hat{\xi}^{x}_{f}=\hat{\xi}^{+}_{f}+\hat{\xi}^{-}_{f}$
(S18b) $\displaystyle\hat{\xi}^{y}_{f}=-i\hat{\xi}^{+}_{f}+i\hat{\xi}^{-}_{f}$
(S18c)

Then the non-vanishing covariance is

$\displaystyle\langle\hat{\xi}_{f}^{x}(t)\hat{\xi}_{f}^{x}(t^{\prime})\rangle=\kappa\left[\hat{C}(t)\hat{n}(t)+\hat{N}(t)\hat{\sigma}^{aa}(t)\right]\delta(t-t^{\prime})$
(S19a)
$\displaystyle\langle\hat{\xi}_{f}^{y}(t)\hat{\xi}_{f}^{y}(t^{\prime})\rangle=\kappa\left[\hat{C}(t)\hat{n}(t)+\hat{N}(t)\hat{\sigma}^{aa}(t)\right]\delta(t-t^{\prime})$
(S19b)
$\displaystyle\langle\hat{\xi}_{f}^{aa}(t)\hat{\xi}_{f}^{aa}(t^{\prime})\rangle=\frac{\kappa}{2}\hat{C}(t)\hat{n}(t)\delta(t-t^{\prime})$
(S19c)
$\displaystyle\langle\hat{\xi}_{f}^{aa}(t)\hat{\xi}_{f}^{x}(t^{\prime})\rangle=i\frac{\kappa}{2}\hat{C}(t)\hat{\sigma}^{y}(t)\delta(t-t^{\prime})$
(S19d)
$\displaystyle\langle\hat{\xi}_{f}^{aa}(t)\hat{\xi}_{f}^{y}(t^{\prime})\rangle=-i\frac{\kappa}{2}\hat{C}(t)\hat{\sigma}^{x}(t)\delta(t-t^{\prime})$
(S19e)

Introducing a vectorial representation
$\bm{\hat{\xi}}(\bm{r},t)=\begin{pmatrix}\hat{\xi}_{l}^{x}(\bm{r},t)&\hat{\xi}_{l}^{y}(\bm{r},t)&\hat{\xi}_{l}^{aa}(\bm{r},t)\end{pmatrix}$
, one can put the the noise correlations covariance in a matrix form

$\langle\bm{\hat{\xi}}(\bm{r},t)\bm{\hat{\xi}}(\bm{r}^{\prime},t^{\prime})\rangle=\delta\left(\bm{r}-\bm{r}^{\prime}\right)\delta\left(t-t^{\prime}\right)\bm{\hat{M}}$
(S20)

where $\bm{\hat{M}}$ is of the following form

$\begin{pmatrix}\kappa\left(\hat{C}\hat{n}+\hat{N}\hat{\sigma}^{aa}\right)+\hat{n}&-i\hat{n}&\hat{\sigma}^{-}-i\frac{\kappa}{2}\hat{C}\hat{\sigma}^{y}\\\
i\hat{n}&\kappa\left(\hat{C}\hat{n}+\hat{N}\hat{\sigma}^{aa}\right)+\hat{n}&i\hat{\sigma}^{-}+i\frac{\kappa}{2}\hat{C}\hat{\sigma}^{x}\\\
\hat{\sigma}^{+}+i\frac{\kappa}{2}\hat{C}\hat{\sigma}^{y}&-i\hat{\sigma}^{+}-i\frac{\kappa}{2}\hat{C}\hat{\sigma}^{x}&\frac{\kappa}{2}\hat{C}\hat{n}+\hat{\sigma}^{aa}\end{pmatrix}$
(S21)

## II. Coarse-grained dynamics

We coarse-grain the dynamics over the facilitation radius $R_{\text{fac}}$ to
obtain a set of equations for the density fields,

$\displaystyle\rho\left(\bm{r},t\right)=\mathcal{V}_{\text{fac}}^{-1}\sum_{l}\theta\left(R_{\text{fac}}-\left|\bm{r}_{l}-\bm{r}\right|\right)\text{Tr}\left\\{\hat{\sigma}^{rr}_{l}\hat{\rho}\right\\}$
(S22a) $\displaystyle
n\left(\bm{r},t\right)=\mathcal{V}_{\text{fac}}^{-1}\sum_{l}\theta\left(R_{\text{fac}}-\left|\bm{r}_{l}-\bm{r}\right|\right)\text{Tr}\left\\{\hat{n}_{l}\hat{\rho}\right\\}$
(S22b)
$\displaystyle\sigma^{x}\left(\bm{r},t\right)=\mathcal{V}_{\text{fac}}^{-1}\sum_{l}\theta\left(R_{\text{fac}}-\left|\bm{r}_{l}-\bm{r}\right|\right)\text{Tr}\left\\{\hat{\sigma}_{l}^{x}\hat{\rho}\right\\}$
(S22c)
$\displaystyle\sigma^{y}\left(\bm{r},t\right)=\mathcal{V}_{\text{fac}}^{-1}\sum_{l}\theta\left(R_{\text{fac}}-\left|\bm{r}_{l}-\bm{r}\right|\right)\text{Tr}\left\\{\hat{\sigma}_{l}^{y}\hat{\rho}\right\\}$
(S22d)

where $\mathcal{V}_{\text{fac}}$ is the volume of the facilitation sphere. In
the long-wavelength regime, we adopt the substitution as follows

$\sum_{\text{nn}}O_{l,\text{nn}}\rightarrow\mathcal{V}_{\text{shell}}\left(1+\frac{R_{\text{fac}}^{2}}{2}\nabla^{2}\right)O_{l}$
(S23)

where $\mathcal{V}_{\text{shell}}=8\pi R_{\text{fac}}^{2}r_{\text{fac}}$ is
the volume of the facilitation shell. As a final step, we rescale the rates
$\kappa$, $\Omega$ as
$\kappa\rightarrow\kappa\mathcal{V}_{\text{shell}},~{}\Omega\rightarrow\Omega\mathcal{V}_{\text{shell}}$.

The coarse-grained Langevin equations read

$\displaystyle\partial_{t}\rho=$
$\displaystyle-\rho+\Omega\sigma^{y}\left(1+D_{\rho}\nabla^{2}\right)\rho$
(S24a)
$\displaystyle+\kappa\left(n-2\rho\right)\left(1+D_{\rho}\nabla^{2}\right)\rho+\xi^{\rho}$
$\displaystyle\partial_{t}\sigma^{x}=$
$\displaystyle-\frac{\sigma^{x}}{2}\kappa\left(1+D_{\rho}\nabla^{2}\right)n-\kappa\sigma^{x}\left(1+D_{\rho}\nabla^{2}\right)\rho$
$\displaystyle-\frac{\gamma}{2}\sigma^{x}-\Omega\sigma^{y}\left(1+D_{\rho}\nabla^{2}\right)\sigma^{x}+\xi^{x}$
(S24b) $\displaystyle\partial_{t}\sigma^{y}=$
$\displaystyle-\frac{\sigma^{y}}{2}\kappa\left(1+D_{\rho}\nabla^{2}\right)n-\kappa\sigma^{y}\left(1+D_{\rho}\nabla^{2}\right)$
$\displaystyle-\frac{\gamma}{2}\sigma^{y}+\Omega\sigma^{x}\left(1+D_{\rho}\nabla^{2}\right)\sigma^{x}\rho$
(S24c)
$\displaystyle+2\Omega\left(n-2\rho\right)\left(1+D_{\rho}\nabla^{2}\right)\rho+\xi^{y}$
$\displaystyle\partial_{t}n=$ $\displaystyle D_{T}\nabla^{2}n-b\rho+\xi^{n}$
(S24d)

Here the diffusion coefficient

$D_{\rho}=\kappa R_{\text{fac}}^{2}/2+D_{T}$ (S25)

and $D_{T}$ is the thermal diffusion coefficient of atoms.

The classical noise covariance $M_{lm}$ is the Hermitian parts of the quantum
noise covariance (S21), i.e.,
$M_{lm}=\langle\xi^{l}\xi^{m}\rangle=\left\langle\hat{M}_{ml}+\hat{M}_{lm}\right\rangle/2$,
where $l,m\in\\{x,y,\rho\\}$. Since the noise sources are Markovain, we
discard the contributions of the spatial gradients to the noise covariance,
then the covariance matrix $\bm{M}$ reads

$\bm{M}=\begin{pmatrix}\left(1+2\kappa\rho\right)n&0&\frac{\sigma^{x}}{2}\\\
0&\left(1+2\kappa\rho\right)n&\frac{\sigma^{y}}{2}\\\
\frac{\sigma^{x}}{2}&\frac{\sigma^{y}}{2}&\left(\kappa
n+1\right)\rho\end{pmatrix}$ (S26)

With respect to the total density field $n$, there is

$\langle\xi^{n}\left(\bm{r},t\right)\xi^{n}\left(\bm{r}^{\prime},t^{\prime}\right)\rangle=b\rho(\bm{r},t)\delta\left(\bm{r}-\bm{r}^{\prime}\right)\delta\left(t-t^{\prime}\right)$
(S27)

Given that the number of emitters in the $p$-state is large enough to be
regarded constant throughout the dynamics of other density fields,
contribution of the loading process to (S27) is ignored.

## III. Dynamic path integral

###  .5 Janssen-De Dominicis-Martin-Siggia-Rose procedure

From the Langevin equations (S24a)-(S24c), one can construct a functional
integration representation in terms of the Janssen-de Dominicis formalism
[32], upon introducing an auxiliary Martin-Siggia-Rose response field [33].

Dynamics of the fields $\sigma^{x/y}$ described by Eqs. (S24b) and (S24c)
relaxes over the time scale of $\left(\kappa n+\gamma\right)^{-1}$, which is
typically shorter than those of the two density fields. Treating
$\sigma^{x/y}$ as uniform fields and the slow variables $\rho,n$ as external
parameters, a conditional transitional probability corresponding to Eqs.
(S24b) and (S24c) in Itō’s discretization can be written following Refs. [61,
39]. The resulting transition probability can be recast into a Gaussian
integral, upon introducing the imaginary, response fields
$\tilde{\sigma}^{x},\tilde{\sigma}^{y}$ conjugate to $\sigma^{x},\sigma^{y}$,
repectively, one obtains

$P\left(\sigma^{x}_{j+1},\sigma^{y}_{j+1}|\sigma^{x}_{j},\sigma^{y}_{j}\right)=\prod_{l=x,y}\int_{-i\infty}^{i\infty}\frac{d\tilde{\sigma}^{l}_{j+1}}{2\pi
i}\exp\left[-\tilde{\sigma}^{l}_{j+1}\left(\sigma^{l}_{j+1}-\sigma^{l}_{j}-dt\mathcal{F}_{\sigma^{l},j}\right)+\frac{1}{2}dtM_{ll,j}\left(\tilde{\sigma}^{l}_{j+1}\right)^{2}\right]$
(S28)

where the subscripts $j$ and $j+1$ indicate the discrete time steps,
$dt=t_{j+1}-t_{j}$ is the time step, $M_{ll,j}=\text{Var}(\xi^{l}_{j})$ is the
noise covariance, and $\mathcal{F}_{\sigma^{l},j}$ is the MF part of the
Langevin Eqs. (S24b) and (S24c).

The transition probability during a finite time follows from summation of all
trajectories with fixed initial and ending points, which can be omitted for
nonequilibrium steady states [39]. By summing up all spatiotemporal
configurations instead of trajectories, we obtain a path integral in the
$dt\to 0$ limit

$\mathcal{Z}=\int\mathcal{D}\left[\sigma^{x},\tilde{\sigma}^{x},\sigma^{y},\tilde{\sigma}^{y}\right]\mathcal{J}\left[\sigma^{x},\sigma^{y}\right]e^{-S_{\sigma}}$
(S29)

in which the Jacobian $\mathcal{J}\left[\sigma^{x},\sigma^{y}\right]=1$ upon a
precise ordering of times such that all the response fields appear at a time
larger than or equal to the times of their conjugate fields [62, 39] (i.e. in
Itō’s discretization), and the action $S_{\sigma}$ takes the form

$-\int\left[\tilde{\sigma}^{x}\left(\mathcal{F}_{\sigma^{x}}+\frac{M_{xx}}{2}\tilde{\sigma}^{x}\right)+\tilde{\sigma}^{y}\left(\mathcal{F}_{\sigma^{y}}+\frac{M_{yy}}{2}\tilde{\sigma}^{y}\right)\right]$
(S30)

where the time derivatives $\partial_{t}\sigma^{x/y}$ vanish in the
nonequilibrium steady states. After dropping all the spatial derivatives in
Eqs. (S24b) and (S24c), the generalized forces $\mathcal{F}_{\sigma^{x/y}}$
are given by

$\displaystyle\mathcal{F}_{\sigma^{x}}=$
$\displaystyle-\frac{n\kappa+\gamma}{2}\sigma^{x}-\kappa\sigma^{x}\rho-\Omega\sigma^{y}\sigma^{x}$
(S31a) $\displaystyle\mathcal{F}_{\sigma^{y}}=$
$\displaystyle-\frac{n\kappa+\gamma}{2}\sigma^{y}+\Omega\left(\sigma^{x}\right)^{2}$
(S31b)
$\displaystyle+\left[2\Omega\left(n-2\rho\right)-\kappa\sigma^{y}\right]\rho$

Since the density fields $\rho,n$ serve as external parameters, we exclude the
effects of their fluctuations on the elements of the noise covariance matrix
$\bm{M}$. Those contributing to the action (S30) are $M_{xx},M_{yy}$. From now
on, we use the shorthand notation $\int=\int dtd\bm{r}$. We note that the
action $S_{\sigma}=S_{\sigma,0}+S_{\sigma,\text{int}}$, in which the former
represents respective relaxation of the $\sigma^{x/y}$ with

$\displaystyle S_{\sigma,0}$
$\displaystyle=\int\tilde{\sigma}^{x}\left(\frac{n\kappa+\gamma}{2}\sigma^{x}-\frac{n}{2}\tilde{\sigma}^{x}\right)$
(S32)
$\displaystyle+\int\tilde{\sigma}^{y}\left(\frac{n\kappa+\gamma}{2}\sigma^{y}-\frac{n}{2}\tilde{\sigma}^{y}\right)$

and the latter represents interactions between them.

Likewise, we introduce a response field $\tilde{\rho}$ to rewrite the Langevin
equation (S24a) into a path integral

$\int\mathcal{D}\left[\tilde{\rho},\rho\right]\exp{\left(-\int\tilde{\rho}\left[\left(\partial_{t}-D_{\rho}\nabla^{2}\right)\rho-\mathcal{F}_{\rho}-\frac{M_{\rho\rho}}{2}\tilde{\rho}\right]\right)}$
(S33)

in which the diffusivity $D_{\rho}$ is given in (S25) and

$\mathcal{F}_{\rho}=-\rho+\left[\Omega\sigma^{y}+\kappa\left(n-2\rho\right)\right]\rho$
(S34)

The contributions of $\sigma^{x/y},\rho$ to $D_{\rho}$ are neglected. The
action associated with Eq. (S33) consists of a bare,
$\sigma^{x/y}$-independent part $S_{\rho}$

$S_{\rho}=\int\tilde{\rho}\left[\left(\partial_{t}-D_{\rho}\nabla^{2}\right)\rho+\rho-\kappa\left(n-2\rho\right)\rho-\frac{1+\kappa
n}{2}\rho\right]$ (S35)

and a part $S_{\rho,\sigma}$ arising from the coupling between the slow and
fast variables,

$S_{\rho,\sigma}=-\int\Omega\tilde{\rho}\rho\sigma^{y}$ (S36)

The coarse-grained action $\left\langle S_{\rho,\sigma}\right\rangle$ follows
from

$\left\langle
S_{\rho,\sigma}\right\rangle=-\ln{\left(\frac{\int\prod_{i=x,y}\mathcal{D}\left[\sigma^{i},\tilde{\sigma}^{i}\right]e^{-S_{\rho,\sigma}-S_{\sigma}}}{\int\prod_{i=x,y}\mathcal{D}\left[\sigma^{i},\tilde{\sigma}^{i}\right]e^{-S_{\sigma,0}}}\right)}$
(S37)

In the presence of dephasing and decay, there is $(\kappa n+\gamma)\geq 2$,
the above integration can be evaluated perturbatively.

###  .6 Perturbative integration

In this subsection we present the details of perturbative integration of the
partition function over $\sigma^{x/y},\tilde{\sigma}^{x/y}$.

$\mathcal{Z}=\int\mathcal{D}\left[\sigma^{x},\tilde{\sigma}^{x},\sigma^{y},\tilde{\sigma}^{y}\right]e^{-S}$
(S38)

where the action $S=S_{\sigma}+S_{\rho,\sigma}$, follows from Eqs. (S30) and
(S36).

We start with integrating out $\sigma^{x},\tilde{\sigma}^{x}$, via averaging
$\exp{\left(-\int\frac{1}{2}\bm{\sigma}^{x}\bm{X}_{1}(\bm{\sigma}^{x})^{\text{T}}\right)}$
over
$\exp{\left(-\int\frac{1}{2}\bm{\sigma}^{x}\bm{X}_{0}(\bm{\sigma}^{x})^{\text{T}}\right)}$,
where the vectorial representation
$\bm{\sigma}^{x}=\begin{pmatrix}\sigma^{x}&\tilde{\sigma}^{x}\end{pmatrix}$,
and

$\displaystyle\bm{X}_{0}=\begin{pmatrix}0&\frac{n\kappa+\gamma}{2}\\\
\frac{n\kappa+\gamma}{2}&-n\end{pmatrix}$ (S39a)
$\displaystyle\bm{X}_{1}=\begin{pmatrix}-2\Omega\tilde{\sigma}^{y}&\kappa\rho+\Omega\sigma^{y}\\\
\kappa\rho+\Omega\sigma^{y}&-2\kappa\rho n\end{pmatrix}$ (S39b)

The action $S$ now reduces to

$\displaystyle\langle S\rangle_{\bm{X}_{0}}$
$\displaystyle=\int\left\\{\left[\frac{n\kappa+\gamma}{2}+\frac{16\Omega^{2}n}{\left(\kappa
n+\gamma\right)^{3}}+\kappa\rho\right]\tilde{\sigma}^{y}\sigma^{y}-\frac{1+2\kappa\rho}{2}n\left(\tilde{\sigma}^{y}\right)^{2}-\frac{2\Omega^{2}\left(\sigma^{y}\right)^{2}}{\left(n\kappa+\gamma\right)^{2}}\right\\}$
(S40) $\displaystyle-\int\left[2\Omega\rho\left(n-2\rho\right)+\frac{4\Omega
n}{\left(\kappa
n+\gamma\right)^{2}}\left(1+2\kappa\rho-\frac{4\kappa\rho}{n\kappa+\gamma}\right)\right]\tilde{\sigma}^{y}$
$\displaystyle-\int\left[\Omega\tilde{\rho}\rho-\frac{2\Omega}{\kappa
n+\gamma}\left(1-\frac{2\kappa\rho}{n\kappa+\gamma}\right)\right]\sigma^{y}$

We proceed to compute the average of Eq. (S40) over
$\sigma^{y},\tilde{\sigma}^{y}$ fields. This will lead to an action

$\langle S\rangle_{\bm{X}_{0},\bm{Y}_{0}}=-\ln{\left\langle
e^{\int\left[-\frac{1}{2}\bm{\sigma}^{y}\bm{Y}_{1}(\bm{\sigma}^{y})^{\text{T}}+\bm{b}(\bm{\sigma}^{y})^{\text{T}}\right]}\right\rangle_{\bm{X}_{0},\bm{Y}_{0}}}$
(S41)

where the vectorial notation
$\bm{\sigma}^{y}=\begin{pmatrix}\sigma^{y}&\tilde{\sigma}^{y}\end{pmatrix}$,
and the matrices $\bm{Y}_{0},\bm{Y}_{1}$, and $\bm{b}$ are given by

$\displaystyle\bm{Y}_{0}=\begin{pmatrix}0&\frac{n\kappa+\gamma}{2}\\\
\frac{n\kappa+\gamma}{2}&-n\end{pmatrix}$ (S42a)
$\displaystyle\bm{Y}_{1}=\begin{pmatrix}-\frac{4\Omega^{2}}{\left(n\kappa+\gamma\right)^{2}}&\kappa\rho+\frac{16\Omega^{2}n}{\left(n\kappa+\gamma\right)^{3}}\\\
\kappa\rho+\frac{16\Omega^{2}n}{\left(n\kappa+\gamma\right)^{3}}&-2\kappa
n\rho\end{pmatrix}$ (S42b)
$\displaystyle\bm{b}^{\text{T}}=\begin{pmatrix}\Omega\rho\tilde{\rho}-\frac{2\Omega}{n\kappa+\gamma}\left(1-\frac{2\kappa\rho}{n\kappa+\gamma}\right)\\\
2\Omega\rho\left(n-2\rho\right)+\frac{4n\Omega}{\left(n\kappa+\gamma\right)^{2}}\left(1+2\kappa\rho-\frac{4\kappa\rho}{n\kappa+\gamma}\right)\end{pmatrix}$
(S42c)

The leading and subleading order corrections $S_{1},S_{2}$ arise from
$\bm{b}\bm{Y}_{0}^{-1}\bm{b}$ and
$\bm{b}\bm{Y}_{0}^{-1}\bm{Y}_{1}\bm{Y}_{0}^{-1}\bm{b}$, where

$\displaystyle S_{1}=$
$\displaystyle-\frac{2n\Omega^{2}\rho^{2}\tilde{\rho}^{2}}{\left(n\kappa+\gamma\right)^{2}}-\frac{4\Omega^{2}\tilde{\rho}}{n\kappa+\gamma}\left(n\rho^{2}-2\rho^{3}\right)$
$\displaystyle-\frac{16n\Omega^{2}\kappa\tilde{\rho}\rho^{2}}{\left(n\kappa+\gamma\right)^{3}}\left(1-\frac{1}{n\kappa+\gamma}\right)$
(S43a) $\displaystyle S_{2}=$
$\displaystyle-\frac{256n^{2}\Omega^{4}\tilde{\rho}\rho}{\left(n\kappa+\gamma\right)^{7}}+\frac{96n^{2}\Omega^{4}\tilde{\rho}^{2}\rho^{2}}{\left(n\kappa+\gamma\right)^{6}}+\frac{8n\Omega^{2}\kappa\tilde{\rho}\rho^{3}}{\left(n\kappa+\gamma\right)^{2}}$
$\displaystyle+\frac{16n\Omega^{2}\kappa\tilde{\rho}\rho^{2}}{\left(n\kappa+\gamma\right)^{3}}+...$
(S43b)

where we retain the leading order corrections to each coupling constant. The
effective density field action in the main text is the sum of Eqs. (S35),
(S43a), (S43b).

###  .7 MF static phase diagram

Figure S6: MF results. (a) Phase boundaries. Static phase diagrams with
conserved ($b,\lambda=0$) total density for (b) $n=3$, (c) $\kappa=0.3$, and
(d) $\Omega=0.5$. The color represents the density of active emitters. The
lines for discontinuous (dashed) and continuous (solid) transitions end at the
tricritical points (white dots).

Equipped with the effective action for the active densities, we can determine
the MF static phases by the noiseless saddle-point solutions to the action
with the total density $n$ fixed. As shown in Fig. S6 (b), in the classical
regime ($\kappa\gg\Omega$), the systems undergo continuous APTs when the
incoherent activation rate $\kappa$ exceeds a critical value. In the quantum
regime ($\Omega\gg\kappa$), the order parameter changes discontinuously when
the coherent activation rate $\Omega$ is larger than the critical one. With
the activation rates serving as external parameters, one can determine the
phase boundaries $n(\kappa,\Omega)$. Discontinuous APTs are characterized by
the bistable region enclosed by the phase boundaries [S6 (a)], with the lower
surface corresponding to the dashed lines in the static phase diagrams [S6 (c)
and (d)]. Increasing the total density $n$ from below the lower surface, the
density of the ordered phase changes discontinuously from zero to a finite
value [Fig. S6 (c) and (d)]. In contrast, when there is a unique boundary
separating the active phase from the absorbing one, the order parameter
increases continuously with $n$.

## IV. functional renormalization group analysis

In this section, we provide details on the non-perturbative renormalization
group (NPRG) approach to discontinuous APTs based on the Wetterich equation
[38].

$\partial_{k}\Gamma_{k}=\frac{1}{2}\text{Tr}\left[\left(\Gamma^{(2)}_{k}+R_{k}\right)^{-1}\partial_{k}R_{k}\right]$
(S44)

Here, in out-of-equilibrium settings, $\Gamma_{k}$ is the nonequilibrium
equivalent of the Gibbs free energy in the presence of a mass-like regulator
$R_{k}$ [39, 40], defined as the Legendre transform of the generating
functional of connected functions
$\mathcal{W}_{k}=\ln{\mathcal{Z}_{k}\left[j,\tilde{j}\right]}$, viz.

$\Gamma_{k}[\rho,\tilde{\rho}]=\int{\left(j\rho+{\tilde{j}}{\tilde{\rho}}\right)}-\mathcal{W}_{k}\left[j,\tilde{j}\right]$
(S45)

where the scale-dependent partition function $\mathcal{Z}_{k}$ is built from
adding a momentum-dependent term $\Delta\mathcal{S}_{k}$ to the bare action
$\mathcal{S}$ evaluated with a uniform field:

$\mathcal{Z}_{k}\left[j,\tilde{j}\right]=\int\mathcal{D}\left[\rho,\tilde{\rho}\right]\exp{\left[-\mathcal{S}-\Delta\mathcal{S}_{k}+\int{\left(j\rho+{\tilde{j}}{\tilde{\rho}}\right)}\right]}$
(S46)

and

$\Delta\mathcal{S}_{k}=\frac{1}{2}\int{\bm{\rho}^{\text{T}}\bm{R}_{k}\bm{\rho}}$
(S47)

where $\bm{\rho}^{\text{T}}=\begin{pmatrix}\rho&\tilde{\rho}\end{pmatrix}$,
and $\bm{R}_{k}$ is a $2\times 2$ matrix of quadratic functions that decouple
short-wavelength modes from long-wavelength ones. By suppressing the low-
energy fluctuations, the regulated action $\Gamma_{k}$ interpolates smoothly
between the microscopic action $\Gamma_{0}$ and the original MF one
$\Gamma_{\Lambda}=\mathcal{S}$ as the cutoff momentum $k$ varies within the
range $k\in\left[0,\Lambda\right]$. Rather than evaluating the exact equation
(S44) directly, which is numerically demanding, we resort to the well-
established local potential approximation (LPA) ansatz [40], according to
which the bare interaction part of the action $\Gamma_{k}$ is evaluated with a
uniform background field

$V^{-1}\ln{\Gamma_{0,k}}=U_{k}\left(\delta_{j}\ln{\mathcal{Z}_{0,k}}|_{j,\tilde{j}},\delta_{\tilde{j}}\ln{\mathcal{Z}_{0,k}}|_{j,\tilde{j}}\right)$
(S48)

where $V$ is the space volume and the running local action (analogous to the
local potential in the equilibrium cases)
$U_{k}\left(\rho,\tilde{\rho}\right)$ is given by

$U_{k}=\tilde{\rho}\left[u_{2,k}\rho+u_{3,k}\rho^{2}+u_{4,k}\rho^{3}\right]-\mu_{1,k}\tilde{\rho}\rho-\mu_{2,k}\tilde{\rho}\rho^{2}$
(S49)

The two (running) noiseless saddle-point solutions
{$\tilde{\rho}=0,\rho=\alpha_{k}$} and {$\tilde{\rho}=0,\rho=0$} to the action
$U_{k}$ correspond to the densities of the absorbing and active phases,
respectively, which are also the running local minima of the effective
potential

$\Phi_{k}(\rho)=\frac{u_{2,k}}{2}\rho^{2}+\frac{u_{3,k}}{3}\rho^{3}+\frac{u_{4,k}}{4}\rho^{4}$
(S50)

With a convenient choice of the theta regulator [40],

$\bm{R}_{k}(q)=\begin{pmatrix}0&R_{k}(q)\\\ R_{k}(q)&0\end{pmatrix}$ (S51)

where
$R_{k}(q)=k^{2}\left(1-q^{2}/k^{2}\right)\theta\left(1-q^{2}/k^{2}\right)$,
the flow equation for the dimensionless potential $U\to k^{-(d+2)}U$ in terms
of the dimensionless fields $\rho\to k^{-d}\rho,\tilde{\rho}\to\tilde{\rho}$
and the RG time $s=\ln{(k/\Lambda)}$ with vanishing anomalous dimension
follows from integrating over frequencies and momenta,

$\partial_{s}U=-\left(d+2\right)U+d\rho
U^{(1,0)}+\frac{V_{d}}{\sqrt{1-\frac{U^{(2,0)}U^{(0,2)}}{\left(U^{(1,1)}+1\right)^{2}}}}$
(S52)

where $V_{d}=\left[2^{d-1}d\pi^{d/2}\Gamma(d/2)\right]^{-1}$, and we have
introduced the notation
$U^{(n,m)}=\partial_{\rho}^{n}\partial_{\tilde{\rho}}^{m}U$ and scaled away
the $k$-independent diffusivity $D_{\rho}$. Hereafter, we omit the explicit
dependence on $s$ of the running parameters.

We evaluate the flow equation (S52) at the running uniform, noiseless saddle
point solution with finite active densities and arrive at

$\displaystyle\partial_{s}u_{2}=$
$\displaystyle-2u_{2}-\frac{V_{d}\left(c_{2}-\alpha
c_{3}\right)\mu_{1}}{\left(1+c_{1}\right)^{2}}-\frac{V_{d}c_{2}\alpha\left[\mu_{1}c_{2}-3c_{3}\alpha\left(\mu_{1}+\alpha\mu_{2}\right)\right]}{\left(1+c_{1}\right)^{3}}+\mathcal{O}\left[\left(1+c_{1}\right)^{-4}\right]$
(S53a) $\displaystyle\partial_{s}u_{3}=$
$\displaystyle-(2-d)u_{3}-\frac{V_{d}\left[c_{3}\left(\mu_{1}-\alpha\mu_{2}\right)+c_{2}\mu_{2}\right]}{(1+c_{1})^{2}}+\frac{2V_{d}c_{2}\left[\mu_{1}\left(c_{2}-6\alpha
c_{3}\right)-\alpha
c_{2}\mu_{2}\right]}{(1+c_{1})^{3}}+\mathcal{O}\left[\left(1+c_{1}\right)^{-4}\right]$
(S53b) $\displaystyle\partial_{s}u_{4}=$
$\displaystyle-2(1-d)u_{4}-\frac{V_{d}c_{3}\mu_{2}}{(1+c_{1})^{2}}+\frac{V_{d}\left\\{c_{3}^{2}\alpha\left(\mu_{1}+\alpha\mu_{2}\right)+c_{2}\left[3c_{3}\left(\mu_{1}+2\alpha\mu_{2}\right)+2c_{2}\mu_{2}\right]\right\\}}{(1+c_{1})^{3}}+\mathcal{O}\left[\left(1+c_{1}\right)^{-4}\right]$
(S53c) $\displaystyle\partial_{s}\mu_{1}=$
$\displaystyle-2\mu_{1}-\frac{2V_{d}\mu_{1}\mu_{2}}{(1+c_{1})^{2}}-\frac{4V_{d}\left\\{\alpha\mu_{1}\mu_{2}\left(c_{2}+10\alpha
c_{3}\right)+\left[\alpha^{2}\mu_{2}^{2}(9c_{2}+11\alpha
c_{3})+\mu_{1}^{2}(\alpha
c_{3}-c_{2})\right]\right\\}}{(1+c_{1})^{3}}+\mathcal{O}\left[\left(1+c_{1}\right)^{-4}\right]$
(S53d) $\displaystyle\partial_{s}\mu_{2}=$
$\displaystyle-(2-d)\mu_{2}-\frac{2V_{d}\mu_{2}^{2}}{(1+c_{1})^{2}}+\frac{2V_{d}\left[8c_{2}\mu_{2}\left(\mu_{1}+2\alpha\mu_{2}\right)+c_{3}\left(2\mu_{1}^{2}+13\alpha\mu_{1}\mu_{2}+13\alpha^{2}\mu_{2}^{2}\right)\right]}{(1+c_{1})^{3}}+\mathcal{O}\left[\left(1+c_{1}\right)^{-4}\right]$
(S53e)

where $\alpha$ is the running minimum related to the active phase, and we have
introduced
$c_{1}=u_{2}+2u_{3}\alpha+3u_{4}\alpha^{2},c_{2}=2u_{3}+6u_{4}\alpha,$ and
$c_{3}=6u_{4}$ to obtain a friendlier expression.

Starting with the MF values as the initial conditions, we numerically
integrate the flow equations until either the two minima become degenerate or
the values of the dimensionful counterparts of the running coefficients no
longer evolve with decreasing $s$. Plugging the renormalized coefficients into
(S50) and determining the phases accordingly yield the results shown in the
main text.

Figure S7: Simulation results in one and two dimensions. (a) Time series of
the total and active densities for upper: $\kappa=0.5,\Omega=0$, and lower:
$\kappa=0,\Omega=0.5$, with $\lambda=10^{-3}$. (b) Snapshots of the active
density fields taken at every half the period $T$ starting with time $t=t_{0}$
for $L=256$. The black and red dotted frames correspond to
$t-t_{0}=lT,(1+l)T/2,l=0,1,2...$, respectively, same as the dotted lines in
the associated time series show in (c). (d) The two-time correlation function
(right) and its Fourier spectra of (c). (e) Time series for $d=1,L=10^{4}$,
whose two-time correlation function and the Fourier spectra are shown in (f).
Parameters are $D_{\rho,T}=1,\gamma=2,\lambda=1.2\times 10^{-3}$.

## V. Numerical integration scheme

Numerical integration of the Langevin equations with multiplicative noise is
performed by an operator-splitting scheme [41, 42], which consists in
integrating the stochastic part first, by sampling the time-dependent solution
of the corresponding Fokker-Planck equation (FPE), and then using the
generated value to evolve the deterministic part of the equation by any
standard numerical integration method.

The Langevin equation for the active density field reads

$\partial_{t}\rho=\alpha+\beta\rho+f(\rho)+\sigma\sqrt{\rho}\xi$ (S54)

where

$\alpha=\tau
n+\frac{D_{\rho}}{\left(dx\right)^{2}}\sum_{i=1}^{2d}\rho\left(\bm{r}+\bm{e}_{i},t\right)$
(S55)

arising from discretizing the Laplacian $\nabla^{2}\rho$ of site $\bm{r}_{i}$
on a $d$-dimensional square lattice of mesh size $dx$, and

$\displaystyle\beta=-u_{2}-2d\frac{D_{\rho}}{\left(dx\right)^{2}}$ (S56a)
$\displaystyle f(\rho)=-u_{3}\rho^{2}-u_{4}\rho^{3}$ (S56b)
$\displaystyle\sigma^{2}=1+\kappa
n+\frac{4\Omega^{2}n\rho}{\left(n\kappa+\gamma\right)^{2}}$ (S56c)

The stochastic value $\rho^{\ast}$ is generated according to the conditional
transition probability density function
$P\left(\rho,t\right)=P\left(\rho(t)=\rho|\rho(0)=\rho_{0}\right)$, which
reads

$P\left(\rho,t\right)=\lambda e^{-\lambda\left(\rho_{0}e^{\beta
t}+\rho\right)}\left[\frac{\rho}{\rho_{0}e^{\beta
t}}\right]^{\mu/2}I_{\mu}\left(2\lambda\sqrt{\rho_{0}\rho e^{\beta t}}\right)$
(S57)

where $I_{\mu}$ is a Bessel function of order $\mu$, and

$\lambda=\frac{2\beta}{\sigma^{2}\left(e^{\beta
t}-1\right)},\,\mu=-1+\frac{2\alpha}{\sigma^{2}}$ (S58)

Using the Taylor-series expansion of the Bessel function, we can sample the
value $\rho^{\ast}$ according to the following mixture

$\rho^{\ast}\sim\text{Gamma}\left[\mu+1+\text{Poisson}\left[\lambda\rho_{0}e^{\beta
t}\right]\right]/\lambda$ (S59)

and then we use $\rho^{\ast}$ as the initial condition for the remaining part
of the equation

$\rho\left(\bm{r},t+dt\right)=\rho^{\ast}+f\left(\rho^{\ast}\right)dt$ (S60)

The equation for the total density is integrated via the Euler’s method. For
the simulations we fix $\gamma=2,D_{T}=D_{n}=1,\tau=10^{-7},b=0.01$.

## VI. Simulation results for low dimensionalities

Simulation results for $d=1,2$ are shown in Fig. S7. In two dimensions, self-
organization in systems undergoing continuous and discontinuous APTs gives
rise to stationary and oscillatory phases [Fig. S7 (a) upper and lower panel].
Within a system displaying a CTC phase, the active densities of various sites
oscillate periodically in synchrony [Fig. S7 (b)], leading to oscillatory
average active and total densities [Fig. S7 (c)]. Furthermore, the associated
two-time correlation function $G_{\rho}(\Delta t)$ maintains periodic for a
finite time, with its Fourier spectra $G_{\rho}(\omega)$ peaking at the
inherent frequency $\omega_{m}$ [Fig. S7 (d)]. Compared with CTCs in $d=3$,
the coherence time for two-dimensional CTCs is significantly shorter, and
there are fewer peaks in the Fourier spectra, indicating that the lower-
dimensional CTCs are less stable. The situation is very different for $d=1$,
where the time series [Fig. S7 (e)] show no sign of periodic oscillations and
the autocorrelation function as well as the spectra reveal a lack of time
crystalline order [Fig. S7 (f)]. This is consistent with the absence of the
first-order APTs in $d=1$.