# Space-time symmetry breaking in nonequilibrium frustrated magnetism

Mingxi Yue  Wilczek Quantum Center and Key Laboratory of Artificial Structures
and Quantum Control, School of Physics and Astronomy, Shanghai Jiao Tong
University, Shanghai 200240, China  Zi Cai<EMAIL_ADDRESS>Wilczek Quantum
Center and Key Laboratory of Artificial Structures and Quantum Control, School
of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240,
China  Shanghai Research Center for Quantum Sciences, Shanghai 201315, China

###### Abstract

Spontaneous symmetry breaking is responsible for the rich phenomena in
equilibrium physics. Driving a system out-of-equilibrium can significantly
enrich the possibility of spontaneous symmetry breaking, which occurs not only
in space, but also in time domain. This study investigates a driven-
dissipative frustrated magnetic system. Results show that frustration in such
a far-from-equilibrium system could lead to a wealth of intriguing non-
equilibrium phases with intertwined space-time symmetry breaking, e.g., a
discrete time crystal phase accompanied by a time-dependent spatial order
oscillating between a long-range tripartite stripe and a short-range
ferromagnetic order.

Introduction – Frustration arises when interacting energies cannot be
simultaneously minimized for all bonds in a many-body system. It hosts
remarkable phenomena ranging from classical spin glassMezard et al. (1986) to
quantum spin liquidLacroix et al. (2011). Typically, frustration may lead to
macroscopic degeneracy in the classical ground-state manifold. However, this
degeneracy could be lifted by thermal or quantum fluctuations, which may
select particular configurations out of the degenerate manifoldVillain (1979);
Henley (1997); Chubukov (1992), or make a superposition among them to form
exotic quantum states without classical counterpartsAnderson (1973). For the
past decades, frustrated systems have been one of the central themes in
condensed matter physics. However, with a few exceptionsWan and Moessner
(2017, 2018); Bittner et al. (2020); Jin et al. (2022), most studies have
restricted their the scope within equilibrium physics, which is governed by
the paradigm of (free) energy minimization. The effect of frustration on
nonequilibrium systems is far from clear, because these systems can absorb
energy from external driving or environment, thus they are usually far from
the ground state and the energy minimization principle does not necessarily
apply.

Compared to equilibrium physics, non-equilibrium physics is much richer,
albeit less known. A prototypical example is the paradigm of spontaneous
symmetry breaking(SSB) that plays a crucial role in both equilibrium and
nonequilibrium physics. In contrast to SSB within the thermal equilibrium,
which is rooted in the variational principle of (free) energy minimization,
the richness of the spatiotemporal structures spontaneously emerging from far-
from-equilibrium systems can only be understood within a dynamical framework,
even for a nonequilibrium steady stateM.Cross and Greenside (2009). The
nonequilibrium phases of matter fundamentally differ from the equilibrium ones
in the sense that the time dimension plays an equally, if not more, important
role than the spatial dimensions in the classification of phases of matter.
For instance, incorporating the time direction significantly enriches the
possibility of SSB, giving rise to a wealth of intriguing nonequilibrium
phases, like the time crystal that spontaneously breaks the discrete or
continuous time translational symmetryWilczek (2012); Bruno (2013); Watanabe
and Oshikawa (2015); Sacha (2015); Else et al. (2016); Khemani et al. (2016);
Yao et al. (2017); Choi et al. (2017); Zhang et al. (2017); Autti et al.
(2018); Träger et al. (2021); Stehouwer et al. (2021); Mi et al. (2022); Frey
and Rachel (2022); Kongkhambut et al. (2022). Frustration is a source of the
exotic phase in equilibrium physics, thus, one may wonder whether it could
lead to novel phases of matter with an intriguing space-time structure in far-
from-equilibrium systems.

In this study, we attempt to address this question, by focusing on a driven-
dissipative interacting spin model. Periodic driving usually heats a generic
closed many-body system toward an infinite temperature state. To avoid this
indefinite energy absorption, we introduce dissipation by coupling each spin
to a heat bath, which will drive the spin system to thermal equilibrium in the
absence of periodic drivingSup . Considering the notorious difficulty of
dealing with quantum many-body dynamics, we focus on a classical system that
enables us to simulate high-dimensional systems up to a large system size. It
has been shown that the exotic nonequilibrium phases are not restricted to
quantum systems (e.g. the time crystal phases in classical many-body systems
have recently been investigatedLupo and Weber (2019); Yao et al. (2020);
Gambetta et al. (2019); Khasseh et al. (2019); Hurtado-Gutiérrez et al.
(2020); Ye et al. (2021); Pizzi et al. (2021)). Different from previous
studies about the prethermal dynamics of close systemsJin et al. (2022); Ye et
al. (2021); Pizzi et al. (2021), here we focus on the long-time asymptotic
behavior of the driven-dissipative system, and show that incorporating
frustration significantly enriches the categories of nonequilibrium phases of
matter, and leads to intriguing magnetic states with intertwined space-time
symmetry breaking.

Figure 1: (Color online)Sketch of a driven-dissipative magnetic system in
square lattice with next-nearest neighboring coupling (blue bonds).

Figure 2: (Color online)In the frustration-free case($\lambda=0$), (a)long-
time dynamics of AF order parameter $m_{A}(t)$ with a weak driving
($J^{\prime}=J$) and (b)long-time dynamics of the AF (red solid) and FM (blue
dash) order parameters $m_{A}(t)$ and $m_{F}(t)$ with a strong driving
($J^{\prime}=5J$). Parameters are chosen as $L=120$, $\eta=J$,
$\mathcal{D}=0.01J$, $h_{z}=1.5J$ and $\omega=2\pi J$.

Model and method – We start with a two-dimensional classical spin model. The
system Hamiltonian reads:

$H_{s}=H_{0}+\lambda H_{f},$ (1)

where $H_{0}$ is a frustration-free Hamiltonian (transverse Ising model)
defined in a $L\times L$ square lattice:

$H_{0}(t)=V(t)\sum_{\langle
ij\rangle}s_{i}^{x}s_{j}^{x}-\sum_{i}h_{z}s_{i}^{z}$ (2)

where the dynamical variable $\mathbf{s}_{i}$ is a three-dimensional classical
vector with a fixed length $|\mathbf{s}_{i}|=1$, and the summation of $\langle
ij\rangle$ is over the bonds connecting two adjacent sites in the square
lattice (green bonds, Fig.1). $h_{z}$ is the strength of the transverse field.
Periodic driving $V(t)=J+J^{\prime}\cos\omega t$ is imposed on the interaction
strength instead of the external field. $V(t)$ can be either positive or
negative, corresponds to antiferromagnetic(AF) or ferromagnetic(FM) coupling.
In our setup, alternate FM and AF couplings during the time evolution is
crucial for our discussion. The frustration is introduced via the Hamiltonian
$H_{f}$, whose strength is characterized by a dimensionless parameter
$\lambda$. The frustration interaction is defined on one diagonal of each
plaquette (blue bonds, Fig.1). The Hamiltonian $H_{f}$ reads:

$H_{f}(t)=V(t)\sum_{\langle\langle ij\rangle\rangle}s_{i}^{x}s_{j}^{x}$ (3)

Only one diagonal of each plaquette is included because in the undriven case
($J^{\prime}=0$) with $\lambda=1$, the Hamiltonian.(1) is reduced to an AF
model defined on a triangle lattice, a prototypical example of frustrated
magnetism. Throughout this paper, we assume $J>0$ and fix the driving
frequency as $\omega=2\pi$. In the absence of a thermal bath, the dynamics of
each spin can be described by the equation of motion (EOM):
$\dot{\mathbf{s}}_{i}=\mathbf{h}^{0}_{i}\times\mathbf{s}_{i}$, where the
effective magnetic field $\mathbf{h}^{0}_{i}=[-V(t)\bar{s}_{i}^{x},0,h_{z}]$
with $\bar{s}_{i}^{x}=\sum_{\langle
j\rangle}s_{j}^{x}+\lambda\sum_{\langle\langle j\rangle\rangle}s_{j}^{x}$.

The periodically-driven system is stabilized by introducing dissipation via
coupling of each spin to a thermal bath, which can be modeled using methods
familiar in the Brownian motion. The EOM for each spin is described by a
stochastic Landau-Lifshitz-Gilbert equationBrown (1963); Ma et al. (2010):

$\dot{\mathbf{s}}_{i}=\mathbf{h}_{i}\times\mathbf{s}_{i}-\eta\mathbf{s}_{i}\times(\mathbf{s}_{i}\times\mathbf{h}_{i})$
(4)

where $\eta$ is the dissipation strength fixed as $\eta=J$ and
$\mathbf{h}_{i}=\mathbf{h}^{0}_{i}+\bm{\xi}_{i}(t)$ is the effective magnetic
field, where $\bm{\xi}_{i}(t)$ is a three-dimensional zero-
mean($\langle\xi_{i}^{\alpha}(t)\rangle_{\bm{\xi}}=0$) stochastic magnetic
field representing a thermal noise. We further assume the local bath around
different sites is independent of each other, and the stochastic variables
satisfy:
$\langle\xi_{i}^{\alpha}(t)\xi_{j}^{\beta}(t^{\prime})\rangle_{\bm{\xi}}=\mathcal{D}^{2}\delta_{\alpha\beta}\delta_{ij}\delta(t-t^{\prime})$
where $\alpha,\beta$ are the index of three space dimensions and $\mathcal{D}$
is the noise strength. The ensemble average $\langle\rangle_{\bm{\xi}}$ is
over all the noise trajectories. For a thermal bath with temperature T, the
strengths of the dissipation and noise satisfy the fluctuation-dissipation
theorem $\mathcal{D}^{2}=2T\eta$. We fix $\mathcal{D}=0.01J$, which
corresponds to an extremely low temperature and does not plays a crucial role
here. The high-temperature case is discussed in the supplementary
material(SM)Sup .

We discretize the stochastic differential Eq.(5) by adopting Stratonovich’s
formula, and solve it by the standard Heun methodAment et al. (2016); Sup
with a time step of $\Delta t=10^{-3}J^{-1}$. In our simulation, we choose the
initial state for each spin as $[s_{i}^{x},0,s_{i}^{z}]$, where
$s_{i}^{x}\in[-1,1]$ is a random number, and $s_{i}^{z}$ is fixed accordingly
as $s_{i}^{z}=\sqrt{1-(s_{i}^{x})^{2}}$ (the initial state dependence has also
been checkedSup ). The system size in our simulation is up to $L=120$. Despite
the richness of the dynamical phase diagram of this driven-dissipative model,
below, we will only focus on the dynamical phases with SSB in both space and
time, and examine the roles of frustration and driving.

Figure 3: (Color online)In a fully frustrated case ($\lambda=1$), (a)
magnetization distribution $\\{s_{i}^{x}\\}$ of the long-time asymptotic state
without driving ($J^{\prime}=0$) in a real space (left panel) and its Fourier
transformation with the peaks at momenta $\mathbf{Q}_{0}$ and
$2\mathbf{Q}_{0}$ with $\mathbf{Q}_{0}=(\frac{2\pi}{3},\frac{2\pi}{3})$ (right
panel); (b)long-time dynamics of tripartite stripe order parameter $m_{3}(t)$
in the weakly driven case ($J^{\prime}=J$); and (c)long-time dynamics of the
tripartite stripe (purple solid) and the FM (blue dash) order parameters
$m_{3}(t)$ and $m_{F}(t)$ in the strongly driven case ($J^{\prime}=5J$).
$L=120$ except for (a) where $L=36$. Other parameters are chosen the same as
those in Fig.2.

Frustration-free case: an antiferromagnetic discrete time crystal (AF-DTC) –
We start with a frustration-free case ($\lambda=0$), where we monitor the
magnetization dynamics $s_{\mathbf{i}}^{x}(t)$ based on Eq.(5). Without
driving($J^{\prime}=0$), the system will relax to an equilibrium AF state with
a nonzero order parameter
$m_{A}=\frac{1}{L^{2}}\sum_{\mathbf{i}}(-1)^{i_{x}+i_{y}}s_{\mathbf{i}}^{x}$.
Fig.2 indicates that such an AF order also persists in the presence of
driving. Fig.2 (a) shows that in the weak driving case ($J^{\prime}=J$) where
the coupling is always AF ($V(t)\geq 0$), $m_{A}(t)$ finally oscillates around
a finite value with a frequency synchronizing with driving. At strong driving
($J^{\prime}=5J$), one can also observe a time-dependent AF order (Fig.2 b),
however, such an AF state fundamentally differs from its equilibrium
counterpart in two aspects. First, the long-range AF order is present at any
time, even at those durations with FM coupling ($V(t)<0$). By contrast, the FM
order parameter
$m_{F}(t)=\frac{1}{L^{2}}\sum_{\mathbf{i}}s_{\mathbf{i}}^{x}(t)$ (the blue
dash line, Fig.2 b) vanishes during the whole evolution. Furthermore,
different from the weakly driven case, $m_{A}(t)$ oscillates with a period
doubling with respect to that of driving, thereby spontaneously breaking the
discrete time translational symmetry from the symmetry group $\mathbb{Z}$ to
$2\mathbb{Z}$. Consequently, such a state simultaneously breaks the space and
time translational symmetry, thus it is an AF-DTC.

Figure 4: (Color online)In a fully frustrated case ($\lambda=1$) with
intermediate driving ($J^{\prime}=3J$), (a) long-time dynamics of the
tripartite stripe (purple solid) and FM (blue dash) order parameters
$m_{3}(t)$ and $m_{F}(t)$ (inset magnifies $m_{F}(t)$), (b)typical
magnetization distribution $\\{s_{i}^{x}\\}$ in a real space at an time slice
with maximal FM coupling; and (c) comparison of the dynamics of excess energy
with respect to the perfect FM state at the time slices with maximal FM
coupling for the cases with intermediate ($J^{\prime}=3J$) and strong
($J^{\prime}=5J$) driving. Other parameters are chosen the same as in Fig.2.

Fully frustrated case: A zoo of nonequilibrium phases of matter – Now we
consider a fully frustrated case with $\lambda=1$, where the lattice is
equivalent to a triangle lattice (a study of intermediate frustration with
$0<\lambda<1$ can be found in SMSup ). In equilibrium magnetism, frustration
works against AF order. One may wonder whether it plays a similar role of
suppressing the forementioned AF-DTC order in this non-equilibrium setup. If
so, which order does frustration suppress? Is it the spatial AF or temporal
DTC, or both? What kinds of space-time structures will emerge once the AF-DTC
is destroyed? In the absence of driving ($J^{\prime}=0$), the system will
relax toward an equilibrium state close to the ground state of the
Hamiltonian.(1). The transverse field distinguishes our model from the pure
Ising model in the triangle lattice with extensive ground state degeneracy.
Fig.3 (a) depicts the steady state magnetization with a tripartite structure
and a stripe order along the diagonal direction. This magnetic order is
characterized by the order parameter
$m_{3}(t)=\frac{1}{L^{2}}\sum_{\mathbf{i}}\sin[\mathbf{Q}_{0}\cdot\mathbf{i}]s_{\mathbf{i}}^{x}(t)$
with $\mathbf{Q}_{0}=(\frac{2\pi}{3},\frac{2\pi}{3})$.

The stripe order parameter $m_{3}(t)$ starts to oscillate once the periodic
driving is switched on. With weak driving ($J^{\prime}=J$), Fig.3 (b) shows
that $m_{3}(t)$ oscillates around a finite value with a period the same as
driving, thus the tripartite stripe order still persists in this
nonequilibrium case. For strong driving ($J^{\prime}=5J$), $V(t)$ alternates
between FM and AF coupling during the evolution, thereby significantly
changing the dynamics. Fig.3 (c) illustrates that at the stripe/FM long-range
order is built during AM/FM coupling. However, this does not mean that the
system adiabatically follows the instantaneous ground state of
Hamiltonian.(1). Instead, it is in a genuine nonequilibrium state because both
$m_{3}(t)$ and $m_{F}(t)$ develop a DTC order in the time domain, indicating
that a temporal correlation is dynamically built. In other words, the
asymptotic state with strong driving is a space-time crystal that
simultaneously breaks not only the space and time translational symmetry, but
also the $Z_{2}$ symmetry of the Hamiltonian.(1) in spin space. This space-
time crystal is significantly different from the AF-DTC phase observed in the
frustration-free case, where the FM order is absent even during FM coupling.

The situation is more interesting with intermediate coupling (e.g.
$J^{\prime}=3J$). Fig.4 (a) illustrates that the tripartite stripe and
associated DTC orders still persist. However, different from that in the
strongly driven case, the long-range FM order is not built during the whole
period. Its order parameter $m_{F}(t)$ (dash blue line,Fig.4 a) stochastically
oscillates with a small amplitude that decreases with system sizeSup . A
typical magnetization configuration $\\{s_{i}^{x}\\}$ at a time slice with a
maximum FM coupling is plotted in Fig.4(b), which exhibits a wealth of FM
domain walls(DW). To measure the density of the DWs, we define an excess
energy density with respect to the perfect FM state as $\delta
e(t)=\frac{1}{L^{2}}\langle H^{I}_{s}(t)\rangle-e_{0}$, where $\langle
H^{I}_{s}(t)\rangle$ is the instantaneous interacting energy at time t (we
only focus on the time slices with maximal FM coupling) and $e_{0}$ is the
energy density of a perfect FM state along x-direction. In Fig.4 (c), for
$J^{\prime}=3J$, $\delta e(t)/e_{0}$ ultimately saturates toward a large
value, indicating that it is far from a perfect FM state and the density of
DWs does not decay in time. By contrast, for the strong driving case with
$J^{\prime}=5J$, excess energy is very small, indicating that the system could
reach an almost perfect FM state at the maximal FM coupling. The results show
that at intermediate coupling, even though the system builds a short-range FM
correlation, it has no time to develop long-range FM order before the coupling
turns back to AF within a driving circle. The difference between the
intermediate and strong coupling cases is because the FM coupling duration in
the former is shorter than that in the latter. The FM duration is not long
enough for the system to build up long-range FM correlation.

Discussion – We present several tunable parameters in our model; hence, one
may wonder whether the forementioned space-time crystals require fine-tuned
parameters in the massive parameter space. We will answer this question by
discussing the role of different model parameters. As for the bath, the
dissipation parameter $\eta$ is introduced to stabilize the system, thereby
controlling the transient relaxation dynamics. Small thermal noises
($\mathcal{D}\ll J$) are introduced to avoid the system from being trapped in
the metastable states. It will not change the nature of the space-time
crystals (Note, however, that a strong thermal fluctuation could melt the
space-time crystalsSup ; Yue et al. (2022)). Including a transverse field
($h_{z}$) makes the spin dynamics nontrivial, otherwise it is only a
precession along the x-direction.

The driving amplitude $J^{\prime}$, frequency $\omega$ and the frustration
strength $\gamma$ are crucial in determining the different space-time patterns
in long-time asymptotic states. As previously discussed, the DTC order can
only emerge when the coupling $V(t)$ alternates between AF and FM during the
evolution. With this, a sufficiently large $J^{\prime}$ is required. As for
$\omega$, in the limit $\omega\gg J$, a fast periodic driving only manifests
itself via its mean value thus doesn’t play a key role here. In the opposite
limit where the driving period is much longer than the relaxation time
$\omega\ll\eta$, the system is always close to the thermodynamical equilibrium
states with FM or stripe orders depending on the sign of $V(t)$. However,
these FM/stripe order parameters cannot spontaneously organize themselves into
a DTC in the time domain because in such an adiabatic limit, the time slices
with maximal AF and FM coupling are separated so far that a temporal
correlation cannot be built up between themSup . Therefore, to realize
nontrivial space-time patterns, periodic driving must be slow enough for the
system to react accordingly and develop different magnetic orders in space,
but not too slow so these magnetic order parameters could build a long-range
(DTC) correlation in time domain. In short, the realization of the space-time
crystals requires a strong driving amplitude and an intermediate driving
frequency. It is a robust phase in the phase diagram.

The frustration in our model suppressed the AF order and induced other
magnetic states, similar as it does in equilibrium magnetism. However, in our
nonequilibrium setup where $V(t)$ keeps changing its sign, frustration may
facilitate an FM order which is absent in the frustration-free case. This is
because in during FM coupling ($V(t)<0$), the Hamiltonian.(3) no longer leads
to “frustration”, instead it increases the effective FM coupling, while during
AF coupling, frustration still suppress the AF order. As a consequence,
frustration makes the AF order observed in the frustration-free case give way
to a FM order. This can be seen more clearly in an intermediate frustration
regime ($\lambda\in[0.68,0.9]$), where we find a DTC phase with only FM but no
stripe orderSup .

Conclusion and outlook – In this work, we introduce frustration into a driven-
dissipative magnetic system, and show that it gives rise to a wealth of
nonequilibrium phases with intriguing SSB in both space and time.

Future developments will include the generalization of these results into
models with different lattice and spin symmetries. For instance, in other
geometric frustrated lattices (e.g. kagome or pyrochlore), one may expect
nonequilibrium phases with other magnetic patterns (e.g. nematic or spin ice)
in space and temporal orders (e.g. algebraic temporal correlationProkofev and
Svistunov (2018)) in time. A more exciting possibility is the realization of
nonequilibrium states with intertwined space-time symmetries that cannot be
decomposed into a direct product of spatial and temporal symmetriesXu and Wu
(2018); Gao and Niu (2021). As for the spin symmetry, the Hamiltonian.(1)
preserves the Ising symmetry, generalizing the spin symmetry to continuous
ones (e.g. $U(1)$) may lead to intriguing non-equilibrium phenomena (e.g. a
Berezinskii-Kosterlitz-Thouless-like phase transition in such a driven-
dissipative system, where the traditional binding-unbinding picture of
vortexKosterlitz and Thouless (1973) may be modified in the context of
nonequilibrium physicsAltman et al. (2015).

A more fundamental question here is to understand the “excitations” in such
space-time crystals, e.g. its response to external stimulation. In equilibrium
physics, elementary excitation is closely related with SSB (e.g. Goldstone
theorem), one may wonder whether similar relation holds for these
nonequilibrium phases with intriguing SSBYang and Cai (2021). If so, how these
excitations (e.g. nonequilibrium counterparts of magnon, soliton or phonon)
are characterized and realized? This problem is not only of fundamental
interest, but also of immense practical importance considering its relevance
to the experimental observable effect of these nonequilibrium phases.

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

## References

  * Mezard et al. (1986) M. Mezard, G. Parisi, and M. Virasoro, _Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications_ ( World Scientific, 1986).
  * Lacroix et al. (2011) C. Lacroix, P. Mendels, and F. Mila, _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_ (Springer, 2011).
  * Villain (1979) J. Villain, Z. Phys. B 33, 31 (1979).
  * Henley (1997) C. L. Henley, J. Appl. Phys. 61, 3962 (1997).
  * Chubukov (1992) A. Chubukov, Phys. Rev. Lett. 69, 832 (1992).
  * Anderson (1973) P. W. Anderson, Materials Research Bulletin 8, 153 (1973).
  * Wan and Moessner (2017) Y. Wan and R. Moessner, Phys. Rev. Lett. 119, 167203 (2017).
  * Wan and Moessner (2018) Y. Wan and R. Moessner, Phys. Rev. B 98, 184432 (2018).
  * Bittner et al. (2020) N. Bittner, D. Golež, M. Eckstein, and P. Werner, Phys. Rev. B 102, 235169 (2020).
  * Jin et al. (2022) H.-K. Jin, A. Pizzi, and J. Knolle, arXiv e-prints arXiv:2204.01761 (2022), eprint 2204.01761.
  * M.Cross and Greenside (2009) M.Cross and H. Greenside, _Pattern Formation and Dynamics in Nonequilibrium Systems_ (Cambridge University Press, Cambridge, 2009).
  * Wilczek (2012) F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012).
  * 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).
  * Sacha (2015) K. Sacha, Phys. Rev. A 91, 033617 (2015).
  * Else et al. (2016) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016).
  * Khemani et al. (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
  * Yao et al. (2017) N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Phys. Rev. Lett. 118, 030401 (2017).
  * Choi et al. (2017) S. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, et al., Nature 543, 221 (2017).
  * Zhang et al. (2017) J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, et al., Nature 543, 217 (2017).
  * Autti et al. (2018) S. Autti, V. B. Eltsov, and G. E. Volovik, Phys. Rev. Lett. 120, 215301 (2018).
  * Träger et al. (2021) N. Träger, P. Gruszecki, F. Lisiecki, F. Groß, J. Förster, M. Weigand, H. Głowiński, P. Kuświk, J. Dubowik, G. Schütz, et al., Phys. Rev. Lett. 126, 057201 (2021).
  * Stehouwer et al. (2021) J. N. Stehouwer, H. T. C. Stoof, J. Smits, and P. van der Straten, Phys. Rev. A 104, 043324 (2021).
  * Mi et al. (2022) X. Mi, M. Ippoliti, C. Quintana, A. Greene, Z. Chen, J. Gross, F. Arute, K. Arya, J. Atalaya, R. Babbush, et al., Nature 601, 531 (2022).
  * Frey and Rachel (2022) P. Frey and S. Rachel, Sci.Adv. 8, 7652 (2022).
  * Kongkhambut et al. (2022) P. Kongkhambut, J. Skulte, L. Mathey, J. G. Cosme, A. Hemmerich, and H. Kessler, Science 337, 670 (2022).
  * (27) See the supplementary material for the details of the Heun algorithm and a benchmark of the numerical results, and a numerical check the dependence of our results on the finite discrete time step $\Delta t$, the system size, noise trajectories, and initial states. The role of the driving frequency, frustration and thermal fluctuation are also discussed.
  * Lupo and Weber (2019) C. Lupo and C. Weber, Phys. Rev. B 100, 195431 (2019).
  * Yao et al. (2020) N. Y. Yao, C. Nayak, L. Balents, and M. P. Zaletel, Nat. Phys. 16, 438 (2020).
  * Gambetta et al. (2019) F. M. Gambetta, F. Carollo, A. Lazarides, I. Lesanovsky, and J. P. Garrahan, Phys. Rev. E 100, 060105 (2019).
  * Khasseh et al. (2019) R. Khasseh, R. Fazio, S. Ruffo, and A. Russomanno, Phys. Rev. Lett. 123, 184301 (2019).
  * Hurtado-Gutiérrez et al. (2020) R. Hurtado-Gutiérrez, F. Carollo, C. Pérez-Espigares, and P. I. Hurtado, Phys. Rev. Lett. 125, 160601 (2020).
  * Ye et al. (2021) B. Ye, F. Machado, and N. Y. Yao, Phys. Rev. Lett. 127, 140603 (2021).
  * Pizzi et al. (2021) A. Pizzi, A. Nunnenkamp, and J. Knolle, Phys. Rev. Lett. 127, 140602 (2021).
  * Brown (1963) W. F. Brown, Phys. Rev. 130, 1677 (1963).
  * Ma et al. (2010) P.-W. Ma, S. L. Dudarev, A. A. Semenov, and C. H. Woo, Phys. Rev. E 82, 031111 (2010).
  * Ament et al. (2016) S. Ament, N. Rangarajan, A. Parthasarathy, and S. Rakheja, arXiv e-prints arXiv:1607.04596 (2016), eprint 1607.04596.
  * Yue et al. (2022) M. Yue, X. Yang, and Z. Cai, Phys. Rev. B 105, L100303 (2022).
  * Prokofev and Svistunov (2018) N. Prokofev and B. Svistunov, J. Exp. Theor. Phys. 127, 860 (2018).
  * Xu and Wu (2018) S. Xu and C. Wu, Phys. Rev. Lett. 120, 096401 (2018).
  * Gao and Niu (2021) Q. Gao and Q. Niu, Phys. Rev. Lett. 127, 036401 (2021).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Altman et al. (2015) E. Altman, L. M. Sieberer, L. Chen, S. Diehl, and J. Toner, Phys. Rev. X 5, 011017 (2015).
  * Yang and Cai (2021) X. Yang and Z. Cai, Phys. Rev. Lett. 126, 020602 (2021).

Supplementary material

In this supplementary material, we first provide some details of the Heun
algorithm as well as a benchmark of the numerical results. Then we numerically
check the convergence of our results with the finite discrete time step
$\Delta t$ and the dependence of our results on the system size, noise
trajectories, and initial states. The role of the driving frequency,
frustration and thermal fluctuation is also discussed.

##  I Details about the numerical simulation

###  I.1 Heun algorithm

Here, we first derive the discrete formalism of stochastic Landau-Lifshitz-
Gilbert (LLG) equation based on Stratonovich’s formula, then formulate the
Heun algorithm to solve this Stochastic Differential equations (SDE). A
stochastic LLG equation reads:

$\dot{\mathbf{s}}_{i}=\frac{1}{\eta^{2}+1}[\mathbf{h}_{i}\times\mathbf{s}_{i}-\eta\mathbf{s}_{i}\times(\mathbf{s}_{i}\times\mathbf{h}_{i})]$
(5)

where $\mathbf{s}_{i}$ is a unit vector.
$\mathbf{h}_{i}(t)=\mathbf{h}^{0}_{i}(t)+\bm{h}^{T}_{i}(t)$ is the effective
magnetic field. $\mathbf{h}^{0}_{i}(t)$ comes from the interaction between
spin i and its neighbors. $\bm{h}^{T}_{i}(t)$ is a three-dimensional random
magnetic field representing the thermal noise, which satisfies:

$\displaystyle\langle h_{i}^{T\alpha}(t)\rangle_{\bm{\xi}}$ $\displaystyle=$
$\displaystyle 0$ (6) $\displaystyle\langle
h_{i}^{T\alpha}(t)h_{j}^{T\beta}(t^{\prime})\rangle_{\bm{\xi}}$
$\displaystyle=$
$\displaystyle\mathcal{D}^{2}\delta_{\alpha\beta}\delta_{ij}\delta(t-t^{\prime})$
(7)

where $\alpha,\beta$ are the index of three spatial dimensions and
$\mathcal{D}$ is the strength of the noise. $\langle\rangle_{\bm{\xi}}$ is the
ensemble average over all the trajectories of noises. According to the
Fluctuation-dissipation theorem, the strength of the thermal noise and the
dissipation satisfies the relation:

$\mathcal{D}^{2}=2T\eta.$ (8)

To solve this SDE numerically, we first discretize the time with the time step
of $\Delta t$. Let the spin configuration in the $m$-th time step
($t_{m}=m\Delta t$) be $\\{\mathbf{s}_{i}^{m}\\}$, the calculation of
$\\{\mathbf{s}_{i}^{m+1}\\}$ can be divided into two steps in the Heun
algorithm.

In the first step, we derive an intermediate spin configuration
$\\{\tilde{\mathbf{s}}_{i}^{m+1}\\}$:

$\tilde{\mathbf{s}}_{i}^{m+1}=\mathbf{s}_{i}^{m}+\frac{1}{\eta^{2}+1}[\mathbf{h}^{m}_{i}\times\mathbf{s}_{i}^{m}-\lambda\mathbf{s}_{i}^{m}\times(\mathbf{s}_{i}^{m}\times\mathbf{h}^{m}_{i})]\Delta
t$ (9)

with $\mathbf{h}_{i}^{m}=\mathbf{h}^{0}_{i,m}+\bm{\tilde{h}}^{T}_{i,m}$, where
$\mathbf{h}^{0}_{i,m}=\mathbf{h}^{0}_{i}(t_{m})$ and
$\bm{\tilde{h}}^{T}_{i,m}$ is a stochastic magnetic field satisfying:

$\tilde{h}^{T\alpha}_{i,m}=\frac{\mathcal{D}}{\sqrt{\Delta
t}}\xi_{i,m}^{\alpha}$ (10)

where $\xi_{i}^{\alpha}$ is a random number satisfying the Gaussian
distribution with zero mean and unit variance:
$\langle\xi_{i}^{\alpha}\rangle_{\bm{\xi}}=0$,
$\langle(\xi_{i}^{\alpha})^{2}\rangle_{\bm{\xi}}=1$.

In the Heun algorithm, $\mathbf{s}_{i}$ at the $m+1$-th time step can be
expressed as:

$\displaystyle\mathbf{s}_{i}^{m+1}$
$\displaystyle=\mathbf{s}_{i}^{m}+\frac{\Delta
t}{2}\\{\mathbf{h}^{m}_{i}\times\mathbf{s}_{i}^{m}-\lambda\mathbf{s}_{i}^{m}\times(\mathbf{s}_{i}^{m}\times\mathbf{h}^{m}_{i})$
(11) $\displaystyle+$
$\displaystyle\tilde{\mathbf{h}}^{m+1}_{i}\times\tilde{\mathbf{s}}_{i}^{m+1}-\lambda\tilde{\mathbf{s}}_{i}^{m+1}\times(\tilde{\mathbf{s}}_{i}^{m+1}\times\tilde{\mathbf{h}}^{m+1}_{i})\\}$

where $\tilde{\mathbf{s}}_{i}^{m+1}$ has been defined in Eq.(9), and
$\tilde{\mathbf{h}}_{i}^{m+1}=\mathbf{h}^{0}_{i,m+1}+\bm{\tilde{h}}^{T}_{i,m}$.

Figure 5: (Color online) Relaxation dynamics predicted by Eq.(5) in the
absence of driving. (a) the dynamics of the energy for a single spin model
predicted by the EOM.(5), whose long-time asymptotic value approaches that
predicted by the statistical ensemble average; (b) the distribution of energy
predicted by the EOM.(5), which agrees very well with the Boltzmann
distribution. (c) The dynamics of the energy for a two-spin model predicted by
the EOM.(5) and the value predicted by statistical ensemble average (red
dashed line). The parameters are chosen as $h_{x}=J$ , $\eta=J$ and the
temperature $T=J$.

Figure 6: (Color online) (a) Comparison between the dynamics of FM order
parameter $m_{F}(t)$ with different $\Delta t$ with $J^{\prime}=5J$, $L=30$.
Comparison between the dynamics of FM order parameter $m_{F}(t)$ with
different system size $L$ (b) in the strong driving case ($J^{\prime}=5J$) and
(c) in the intermediate driving case ($J^{\prime}=3J$). Other parameters are
chosen as $\eta=J$, $\mathcal{D}=0.01J$,$\lambda=J$,
$h_{z}=1.5J$,$\omega=2\pi$.

###  I.2 Benchmark: spin models without driving

It is known that once a system couples to a heat bath, it will finally relax
to a thermodynamical equilibrium state at the same temperature of the bath,
irrespective of its initial state. To verify the validity of the EOM.(5), we
consider two simple spin models as benchmarks, which show that the long-time
asymptotic state derived by EOM.(5) is indeed the thermodynamic equilibrium
states at a temperature determined by Eq.(8).

The first model we consider is a single spin model with the Hamiltonian:

$H^{s}_{1}=h_{z}s^{z}$ (12)

By solving EOM.(5) using the Heun algorithm, we can find the energy of the
system $E(t)=\langle H(t)\rangle$ quickly relaxes to a value with small
statistical fluctuation. According to the statistical physics, for a
thermodynamical equilibrium state, the long-time average of the system energy
is supposed to be the same with the value predicted by the statistical
ensemble, which is

$E_{s}=\frac{1}{Z}\int_{0}^{\pi}\sin\theta d\theta[h_{z}\cos\theta]e^{-\beta
h_{z}\cos\theta}$ (13)

where $\theta$ is the angle between the spin vector and z-axis, and
$Z=\int_{0}^{\pi}\sin\theta d\theta e^{-\beta h_{z}\cos\theta}$ is the
partition function. As shown in Fig.5 (a), the time average value of $\langle
H(t)\rangle$ agrees very well with the ensemble averaged value $E_{s}$ within
the statistical error. In addition, one can study the statistical distribution
of $E(t)$ during the long-time dynamics, $P(E)$ agrees very well with the
Boltzmann distribution, as shown in Fig.5 (b).

We also check a two-spin model with the Hamiltonian:

$H_{s}^{2}=Js_{1}^{x}s_{2}^{x}+h_{z}(s_{1}^{z}+s_{2}^{z})$ (14)

We study the evolution of the system, and focus on its energy. As shown in
Fig.5 (c), in the long-time limit, the system energy will approach the value
predicted by canonical ensemble accompanied by small statistical fluctuations.

Figure 7: (Color online) Comparison between the dynamics of FM order parameter
$m_{F}(t)$ with (a) different noise trajectories and (b) different initial
states with $J^{\prime}=5J$ and $L=30$. Other parameters are chosen the same
as in Fig.6.

###  I.3  Convergence of numerical results

Finite $\Delta$t – Throughout the maintext, we choose a discrete time step of
$\Delta t=10^{-3}J^{-1}$. For a stochastic differential equation, the choice
of $\Delta t$ is more subtle than that in the deterministic equation since the
random variable depends on $\Delta t$ as shown in Eq.(10). To examine the
convergence of our result with $\Delta t$, we choose different $\Delta
tJ=2\times 10^{-3},10^{-3}$ and $5\times 10^{-4}$, and compare their results.
As shown in Fig.6 (a), the results with different $\Delta t$ agree with each
other very well, which indicates that the $\Delta t$ chosen in our simulation
is sufficiently small, thus enables us to ignore the errors induced by the
discretization of time.

Finite L – In the maintext, the system we simulated is up to a system size
with $L=120$. One needs to check the system size dependence of our results.
For an ordered phase (e.g., the DTC phase with strong driving
($J^{\prime}=5J$)), as shown in Fig.6 (b), the deviation between the results
with $L=30$, $60$ and $120$ is pretty small. However, in the intermediate
driving regime without true long-range ferromagnetic (FM) order, the FM order
parameter $m_{F}(t)$ strongly depends on the system size. As shown in Fig.6
(c), the amplitude of $m_{F}(t)$ significantly decreases with $L$. This is due
to the fact that in the presence of intermediate driving, the long-range FM
order has not been built up during the FM coupling. Instead, the system is
spontaneously separated into different FM domains, and the overall FM order
parameter is a summation of the magnetization of them. Within each FM domain,
the magnetization oscillates as a DTC, but the phases of these DTCs are not
coherent, thus the magnetization in different domains at a given time cancels
with each others. As a consequence, the overall FM order parameter decreases
with the system size.

The different finite size dependence of the FM order parameter between the
strong and intermediate driving can be considered as the signature of the
different dynamical phases with long-range and short-range FM orders
respectively. In addition, the finite size effect is supposed to be important
near the dynamical critical points, which is an important issue but not
addressed in this work.

Noise trajectories – In principle, one needs to simulate over different noise
trajectories and perform the ensemble average over them. However, in our
simulation, we only randomly choose one noise trajectory for each set of
parameters. This is because we are only interested in the situation with weak
noise($\mathcal{D}=0.01J$), where a small thermal fluctuation does not change
the nature of the phases with discrete symmetry breaking (see the comparison
between the dynamics over two different noise trajectories in Fig.7 a).
However, for some special initial states, it is possible that the system could
be trapped into a metastable state if there is no thermal fluctuations. The
role of noise in our simulation is to thermally activate the system and make
it escape from the metastable state after sufficiently long time and enter the
genuine asymptotic long-time states discussed in the maintext.

Initial states – In our simulation, we start from a spatially inhomogeneous
random initial state: for each site, we choose its initial state as
$\bm{s}_{i}^{0}=[s_{i}^{x},0,s_{i}^{z}]$, where $s_{i}^{x}$ is an random
number different from site by site and uniformly distributed within [-1,1],
the z-component of the spin is chosen correspondingly as
$s_{i}^{z}=\sqrt{1-[s_{i}^{x}]^{2}}$. Since we focus on the non-equilibrium
phases with spontaneously symmetry breaking, it is well known that in this
case, the final state is supposed to be highly sensitive to the initial state.
This statement does not only work for equilibrium phases (FM or AF order), but
also for non-equilibrium phases. For instance, for a DTC phase with
spontaneously $Z_{2}$ time translational symmetry breaking, as show in Fig.7
(b), starting from different initial states, the system could finally fall
into either one of the $Z_{2}$ breaking phases, each of which is related with
the other by a half-period shift in time domain.

Figure 8: (Color online) The dynamics of the 3-period stripe order parameter
$m_{3}(t)$ in the presence of (a) fast driving with a frequency $\omega=8\pi
J$ and (b) slow driving with a frequency $\omega=0.02\pi$. The x-axis are
plotted in the unit of the corresponding driving periods with
$T_{1}=0.25J^{-1}$ and $T_{2}=100J^{-1}$. The parameters are chosen as $L=30$,
$J^{\prime}=5J$. Other parameters except $\omega$ are chosen the same as in
Fig.6.

##  II The role of different parameters in the model

In the maintext, we only studied several representative non-equilibrium phases
by focusing on special points in the parameter space. In the following, we
will systematically examine the role of different parameters of the models in
determining the space-time patterns of dynamical phases. The dynamical phase
transitions between them have also been studied.

Figure 9: (Color online) (a) The time-independent AF ($m_{A}$), FM ($m_{F}$)
and tripartite stripe ($m_{3}$) order parameters as a function of frustration
strength $\lambda$ with $J^{\prime}=5J$ and $\mathcal{D}=0.01J$. (b) In the
strongly driven case ($J^{\prime}=5J$), the DTC order parameter corresponding
to the FM order as a function of noise strength $\mathcal{D}$. (c)In the
intermediate driven case ($J^{\prime}=3J$), the DTC order parameter
corresponding to the 3-period stripe order as a function of noise strength
$\mathcal{D}$. $\lambda=J$ for (b) and (c). Other parameters are chosen as
$L=30$, $\eta=J$, $h_{z}=1.5J$, $\omega=2\pi$.

###  II.1 Drving frequency $\omega$

In the maintext, we changed the driving amplitude but fixed the driving
frequency as $\omega=2\pi J$. Here, we will exam the role of $\omega$ in
determining the space-time patterns of our model.

In the fast driving limit where $\omega\gg J$, the periodic driving oscillates
too fast to be followed by the system. In such a high frequency limit, one can
derive an effective time-independent Hamiltonian to describe the stroboscopic
dynamics of this periodically driven system, similar as the Floquett analysis
in quantum systems, where the effective time-dependent Hamiltonian can be
expressed as a expansion in terms of $1/\omega$. In the high frequency limit,
the dominant term in the expansion is an average of the Hamiltonian over one
period, where $\frac{1}{T_{0}}\int_{0}^{T_{0}}dtV(t)=J$, therefore, the
dynamics in this case is similar with the relaxation dynamics in the undriven
case ($J^{\prime}=0$), where the steady state is a stripe phase with a
nonvanishing order parameter
$m_{3}=\frac{1}{L^{2}}\sum_{\mathbf{i}}\sin[\mathbf{Q}_{0}\cdot\mathbf{i}]s_{\mathbf{i}}^{x}(t)$
with $\mathbf{Q}_{0}=(\frac{2\pi}{3},\frac{2\pi}{3})$. For a large but finite
frequency, the order parameter of the stripe phase will oscillates around its
equilibrium value with a frequency the same as the driving, as shown in Fig.8
(a) where the frequency is chosen as $\omega=8\pi J$.

In the opposite limit of slow driving, where the period of the driving is much
longer than the relaxation time $\omega\ll\eta$ ($\eta$ is the dissipation
strength), the system have sufficient time to relax, thus at any given time,
the system is always close to an equilibrium state. As a consequence, both the
FM and the stripe spatial order can be developed depending on the sign of
$V(t)$ in the instantaneous Hamiltonian. However, a thermalization of a system
means that it has no memory of the information of its initial state, or the
previous states far away from it. As a consequence, for magnetic states with
SSB, the system will randomly choose one state among the degenerate manifold
with SSB, while these symmetry breaking phases at different time slices barely
correlate with each other, thus cannot form a long-range order (DTC) in time
domain, as shown in Fig.8 (b) where the frequency is chosen as
$\omega=0.02\pi$.

###  II.2 Frustration $\lambda$

In the maintext, we only studied the unfrustrated ($\lambda=0$) and full-
frustrated ($\lambda=1$) cases, which exhibits non-equilibrium phases with
different SSB. Here, we will systematically study the role of frustration by
continuously tuning the frustration strength $\lambda$. We focus on strongly
driving case ($J^{\prime}=5J$). We find that the frustration does not change
the DTC nature of the phase, but is crucial in determining the spatial order
of the non-equilibrium phases. In general, the magnetic order parameters in
these dynamics phases keep oscillating in time, therefore, to characterize
their strength, we need to derive a time-independent order parameter. For
instance, for the FM order parameter $m_{F}(t)$, we choose those time slices
with $t=t_{n}^{F}$ when $m_{F}(t)$ reach its n-$th$ maximum, and perform the
average over them to derive a time-independent order parameter $m_{F}=\langle
m_{3}(t_{n}^{F})\rangle$, and use it to characterize the strength of the FM
order parameter in these dynamic phases.

We plot the time-independent AF, FM and tripartite stripe order parameters
$m_{A}$, $m_{F}$ and $m_{3}$ as a function of $\lambda$ in Fig.9 (a), from
which we can find that for a small frustration, the AF-DTC order persist until
$\lambda=0.68J$, where the AF order gives way to a FM order via a first order
(1st) phase transition. For a frustration within the regime
$\lambda\in[0.68J,0.9J]$, one can find an intermediate phase where the long-
range FM order has been built up in the duration of FM coupling, while there
is no long-range tripartite stripe order during the AF coupling. In another
word, frustration suppresses the AF order even in this non-equilibrium driven
case, while it facilitates the FM orders, since in the duration with FM
coupling ($V(t)<0$), the next-nearest-neighboring coupling no longer leads to
“frustration”, instead it increases the effective FM coupling. When the
frustration further increases, the system experience another 1st order
dynamical phase transition at $\lambda=0.9J$ characterized by the sudden onset
of the tripartite stripe order, and the system enters a dynamics phase
oscillating between the states with long-range FM and tripartite stripe
orders, which has been discussed in maintext of the fully-frustrated strongly
driving case.

###  II.3 Thermal fluctuation $\mathcal{D}$

In the maintext, we focus on the case with weak thermal fluctuation
($\mathcal{D}=0.01J$), which does not change the nature of the phases with
discrete symmetry breaking. However, it is known that the thermal fluctuation
works against spontaneous symmetry breaking, and a strong thermal fluctuation
could melt the ordered phase and restore the symmetries. Since the space-time
crystals phase discussed here simultaneously breaks different symmetries
(space and time translational symmetries and $Z_{2}$ spin symmetry), one may
wonder what’s the effect of the thermal fluctuations on these dynamical
orders.

As an example, we focus on the DTC order with spontaneous $Z_{2}$ symmetry
breaking in time domain, and check whether it is possible for thermal
fluctuation to restore this symmetry. Since at low temperature, both
$m_{F}(t)$ and $m_{3}(t)$ exhibit DTC order in time domain, one needs to
distinguish their corresponding DTC order via different order parameters as:

$DTC_{a}=\frac{2}{t_{0}}\int_{\frac{t_{0}}{2}}^{t_{0}}dte^{i\pi t}m_{a}(t)$
(15)

with $a=F$ or $3$ indicates the $FM$ or tripartite stripe order parameter
respectively. $t_{0}=2000J^{-1}$ is our simulation time and the Fourier
transformation is performed over the second half of the full simulation time,
during which the system has reached the long-time asymptotic state.

We focus on both cases with intermediate ($J^{\prime}=3J$) and strong
($J^{\prime}=5J$) driving and study the corresponding DTC order parameters
$DTC_{3}$ and $DTC_{F}$ as a function of $\mathcal{D}$. The ensemble average
is performed over $\mathcal{N}$ noise trajectories, with $\mathcal{N}=10$. For
the case with strong driving, the DTC order parameter corresponding to the FM
order persist until $\mathcal{D}\approx 0.15J$, above which the discrete time
translational symmetry has been restored (Fig. 9 b), however, the DTC order
corresponding to the tripartite stripe order is much more fragile, it vanishes
for $\mathcal{D}>0.03J$ (not shown here). In summary, in the case with strong
driving, there exists an intermediate noise regime where the FM-DTC order
survives but the stripe-DTC order does not, similar with what happened in the
intermediate frustration regime discussed above. In the case with intermediate
driving, the long-range FM order has not been built up even in the presence of
weak noise and only tripartite stripe order exists. However, different from
its counterpart in the strongly driven case, the DTC order corresponding to
this tripartite stripe phase in the intermediately driven case is quite robust
against thermal fluctuation, as shown in Fig. 9 (c). The thermal fluctuation-
induced transitions discussed in these two cases seem to be continuous.
However, a precise determination of the position of the phase transition point
and the critical properties calls for a finite-size scaling analysis, which
will be left in the future work.


# Confined modes of single-particle trajectories induced by stochastic
resetting

Aleksander A. Stanislavsky*  —1<EMAIL_ADDRESS>Aleksander
Weron*  —2<EMAIL_ADDRESS>1Institute of Radio Astronomy, 4
Mystetstv St., 61002 Kharkiv, Ukraine
2Faculty of Pure and Applied Mathematics, Hugo Steinhaus Center, Wrocław
University of Science and Technology, Wyb. Wyspiańskiego 27, 50-370 Wrocław,
Poland

###### Abstract

Random trajectories of single particles in living cells contain information
about the interaction between particles, as well as, with the cellular
environment. However, precise consideration of the underlying stochastic
properties, beyond normal diffusion, remains a challenge as applied to each
particle trajectory separately. In this paper, we show how positions of
confined particles in living cells can obey not only the Laplace distribution,
but the Linnik one. This feature is detected in experimental data for the
motion of G proteins and coupled receptors in cells, and its origin is
explained in terms of stochastic resetting. This resetting process generates
power-law waiting times, giving rise to the Linnik statistics in confined
motion, and also includes exponentially distributed times as a limit case
leading to the Laplace one. The stochastic process, which is affected by the
resetting, can be Brownian motion commonly found in cells. Other possible
models producing similar effects are discussed.

##  I Introduction

Combination of the live-cell single-molecule imaging with single-particle
tracking (SPT) methods has allowed a revolutionary breakthrough for the
quantitative analysis of dynamic processes in living cells [1, 2, 3, 4]. This
approach gives the visualization of the movement of individual “particles”,
the latter being single molecules, macromolecular complexes, as well as,
viruses or organelles in physiological conditions [5]. As such, it has been
crucial to study the mechanisms of intracellular transport, cell membrane
dynamics, and viral infection [6, 7, 8]. Stochastic processes associated with
the movement of particles are directly affected by interactions that occur
with other cellular structures or components [9]. Therefore, single-particle
dynamics often deviates from Brownian motion and exhibits heterogeneous
behavior characterized by changes in diffusion, transient confinement,
immobilization or anomalous diffusion [10, 11, 12]. The development of
theoretical frameworks for the robust analysis of random trajectories
implemented in biological scenarios is thus of fundamental importance to
understand molecular mechanisms of interaction [13].

As an example, detection of the transient confinement with high precision
requires the knowledge of position or displacement statistics [14]. Restrained
trajectories have shown to obey the Gaussian statistics, as well as, the
Laplace one (see, for example, [15, 16]). The Gaussian confinement is often
described by the well-known Ornstein-Uhlenbeck (OU) model [17]. In this case,
the stationary probability distribution function (PDF) is Gaussian, and the
random trajectories follow Gaussian statistics with an obvious physical
interpretation. Unfortunately, the OU model, based on the ordinary Langevin
equation with a harmonic potential, does not provide a description for
confined trajectories with the Laplace statistics.

An alternative scenario, leading to confinement of single-particle
trajectories with the Laplace statistics, can be based on the stochastic
resetting methodology. The resetting of a stochastic process describes
evolution of a stochastic system that is returned repeatedly to a steady (or
equilibrium) state [18], as it occurs in target search with home returns [21,
22, 23]. If the resetting process is Poissonian and independent of the random
motion that undergoes resetting, then the latter has a steady PDF. For
trajectories undergoing Brownian motion with Poissonian resetting, the steady
PDF has the Laplace form. It is noteworthy that subordinated Brownian
processes, leading to subdiffusion, under Poissonian resetting also produces a
stationary state with the Laplace distribution, but with another scale
parameter [24]. In contrast, the engagement of Lévy motion with the stochastic
resetting produces confinement with the Linnik statistics, in which the
Laplace one is a particular case (see [25]). Experimental works have provided
examples of the occurrence of confinement with Laplace statistics [15, 14,
16], but to our knowledge no examples of confined trajectories with Linnik
PDFs have been reported.

In this article, considering the frequent occurrence of stochastic resetting
in biological systems, we propose an analysis pipeline to robustly determine
the testing and apply it to several experimental datasets. Our study shows
that the Linnik statistics does occur in single-particle trajectories of both
G proteins and coupled receptors in living cells. After providing a brief
mathematical background of the stochastic resetting, we describe the
analytical framework and its application to the experiments that led to the
quantification of confined trajectories and their segments following the
Laplace and Linnik statistics. We discuss the results and further propose a
possible explanation for their occurrence. In our consideration, the transient
behavior means that the set of trajectories contains segments in different
diffusive modes: free, confined, and others. Each segment obeys one of the
modes, however trajectories not segmented can be also found. The case when
either segmental or non-segmental trajectories have different statistics
indicates the inhomogeneity of the medium and the mode contribution ratio
determines how heterogeneous the medium is. If all the trajectories follow the
Brownian motion, the medium is homogeneous.

##  II Models of confined modes

There are at least three ways to get a stationary PDF for stochastic
processes, namely

  * •

subordination of random processes;

  * •

stochastic differential equations (SDEs) with an attractive potential;

  * •

stochastic resetting.

Each of them may be considered as a mode leading to confinement. Sometimes,
they have similarities in limit cases, but, generally, the limits are
different. Let us consider the above approaches and their features below.

###  II.1 Model 1: subordination of random processes

The subordination consists of time randomization of a stochastic process with
the help of another independent random process [26]. Confined trajectories
obeying the Laplace PDF can be produced by the use of a specific subordinator,
closely related to the one providing tempered subdiffusion. The conjugate
property of Bernstein functions connects the tempered subdiffusion with the
confinement [24]. Interpretation of anomalous diffusion tending to the
confinement is that diffusive motion, accompanied by multiple-trapping events
with infinite mean sojourn time, is transformed into pure jumps, restricted in
a confined environment. If the Laplace exponent of a tempered $\alpha$-stable
process is $(s+\chi)^{\alpha}-\chi^{\alpha}$, where $\chi$ is a positive
constant and $0<\alpha<1$, then its conjugate partner from the set of
Bernstein functions takes the form $s/((s+\chi)^{\alpha}-\chi^{\alpha})$ [27].
The PDF of the operational time is easy to present in the Laplace transform
with respect to $t$, i.e.,

$\tilde{f}(\tau,s)=\frac{1}{(s+\chi)^{\alpha}-\chi^{\alpha}}\,e^{-\tau
s/((s+\chi)^{\alpha}-\chi^{\alpha})}\,.$ (1)

The propagator of such a subordinated process can be written in the integral
form

$p(x,t)=\int^{\infty}_{0}h(x,\tau)\,f(\tau,t)\,d\tau\,,$ (2)

where $h(x,\tau)$ is the PDF of a parent process, whereas $f(\tau,t)$ is the
PDF of a directing one. If the ordinary Brownian motion is a parent process,
and the directing process is described by Eq. (1), it is not difficult to find
a stationary distribution in the Laplace form explicitly

$p_{B}(x,\infty)=\frac{1}{\sqrt{2D\alpha\chi^{\alpha-1}}}\,\exp\left(-\frac{2|x|}{\sqrt{2D\alpha\chi^{\alpha-1}}}\right)\,,$
(3)

where $D$ is the diffusivity constant for the Brownian motion. Any
subordinated non-Brownian motion, in which the subordinator is defined by the
Laplace exponent conjugate to a tempered $\alpha$-stable process, has a
confined probability distribution itself [14, 16]. Their forms are simpler to
present through the Fourier transform, giving characteristic functions. This
procedure covers a wide class of geometrically infinitely divisible
distributions as a confined case of the non-Brownian motion subordinated by a
special subordinator responsible for the confinement.

###  II.2 Model 2: coupled SDEs with a potential

The coupled Langevin equations for position $x_{t}$ and diffusivity $D_{t}$
opens new possibilities in description of the confinement. The system of SDEs
reads

$\left\\{\begin{array}[]{ll}dx_{t}=-\frac{1}{\tau_{x}}(x_{t}-\bar{x})dt+\sqrt{D_{t}}dW_{t}^{(1)}\\\
dD_{t}=-\frac{1}{\tau_{D}}(D_{t}-\bar{D})dt+\sigma\sqrt{D_{t}}dW_{t}^{(2)}\end{array}\right.\,,$
(4)

where $\tau_{x}$ and $\tau_{D}$ are the correlation times for $x$ and $D$,
$\bar{x}$ and $\bar{D}$ denote the average position and diffusivity, whereas
$\sigma$ is the “speed” of fluctuations for the diffusion coefficient [28].
Here, $W_{t}^{(1)}$ and $W_{t}^{(2)}$ are independent Wiener processes. Note
that the second equation of the system (4) is independent of the first,
whereas the latter is “driven” by the former. Moreover, the former corresponds
to a Cox-Ingersoll-Ross (CIR) process [29]. The stationary solution of the CIR
SDE is presented, for example, in [30]. This PDF is expressed in terms of the
gamma distribution. Next, the stationary state for the first SDE in the system
(4) leads to the Gaussian PDF, if its standard deviation became constant as in
the ordinary OU model. To find the stationary PDF $F_{\lambda}(x)$ of the
coupled equations (4), we integrate the Gaussian and the gamma PDFs over $D$.
Consequently, the stationary PDF takes the following form:

$\displaystyle F_{\lambda}(x)$ $\displaystyle=$
$\displaystyle\frac{2\eta^{-\lambda-1/2}}{\Gamma(\lambda)\sqrt{\pi}}|x-\bar{x}|^{\lambda-1/2}$
(5) $\displaystyle\times$ $\displaystyle
K_{\lambda-1/2}\left(\frac{2|x-\bar{x}|}{\eta}\right)\,,$

where $K_{\nu}(z)$ is the modified Bessel function [31] (here
$\nu=\lambda-1/2$), $\eta=\sqrt{\tau_{x}\tau_{D}\sigma^{2}}$, and
$\lambda=\frac{2\bar{D}}{\tau_{D}\sigma^{2}}$. The index $\lambda-1/2$ has the
meaning of the shape parameter [32]. It is also interesting to mention the
characteristic function for this PDF. It reads

$\int_{-\infty}^{\infty}F_{\lambda}(x)\,\cos(qx)\,dx=\frac{e^{iq\bar{x}}}{\left(1+q^{2}\eta^{2}/4\right)^{\lambda}}\,,$
(6)

showing similarity with both the Laplace and the Linnik distributions. Note,
if $\lambda=1$, then Eq.(5) takes the form of the Laplace PDF

$F_{1}(x)=\frac{1}{\eta}\exp\left(-\frac{2|x-\bar{x}|}{\eta}\right)\,$ (7)

where $\bar{x}$ is the location parameter and $\eta>0$ is the scale parameter
(sometimes referred to as the diversity).

###  II.3 Model 3: stochastic resetting

The Brownian motion under Poissonian resetting can be described by a
starightforward mathematical framework [18, 19, 20]. The Green function in the
absence of resetting (or in other terms the propagator) is given by

$G_{1}(x,t|x_{0})=\frac{1}{\sqrt{4\pi
Dt}}\exp\left[-\frac{(x-x_{0})^{2}}{4Dt}\right]\,,$ (8)

where $x_{0}$ is the initial position, and $D$ is the diffusion constant. The
initial condition takes the Dirac delta-function, namely
$G_{1}(x,0|x_{0})=\delta(x-x_{0})$. Due to the Poissonian resetting with the
rate $r$, the PDF $p_{1}(x,t|x_{0})$ satisfies a renewal equation [21]. Its
form is expressed by a sum of two terms, described as

$\displaystyle p_{1}(x,t|x_{0})$ $\displaystyle=$ $\displaystyle
e^{-rt}G_{1}(x,t|x_{0})$ (9) $\displaystyle+$ $\displaystyle
r\int_{0}^{t}e^{-r\tau}\,G_{1}(x,\tau|X_{r})\,d\tau\,,$

where $X_{r}$ is the position to which the particle returns after resetting.
It should be noticed that the renewal equation Eq. (9) holds for more general
stochastic processes, with propagators having a more general form than
$G_{1}(x,t|x_{0})$ [21, 24]. In the stationary state with $t\to\infty$, the
first term of Eq. (9) may be neglected, and the second term can be exactly
calculated [18], providing the typical Laplace PDF [25]

$p_{1}(x,\infty|x_{0})=\frac{c_{1}}{2}e^{-c_{1}|x-X_{r}|}\,,$ (10)

where the scale parameter $c_{1}=\sqrt{r/D}$ depends on the rate $r$. Hence,
particles undergoing Brownian diffusion with resetting behave as performing
confined motion with the Laplace statistics. Poissonian resetting also yields
a similar behavior for subdiffusion [24], i.e. random walks characterized by a
nonlinear mean squared displacement $t^{\alpha}$ with $\alpha\in(0,1)$. The
only difference in this case, is that the scale parameter is given by
$c_{\alpha}=\sqrt{r^{\alpha}/D}$. A similar result can also be obtained in a
more general case. In fact, for a subordinator with an inverse infinitely-
divisible distribution, the subordinated Brownian motion under Poissonian
resetting tends to a stationary PDF in the Laplace form with the scale
parameter $c_{\Psi}=\sqrt{\bar{\Psi}(r)/D}$, where the Laplace exponent
$\bar{\Psi}(r)$ is expressed in terms of Bernstein functions [33]. Therefore,
Poissonian resetting can force many (but not any) nonstationary stochastic
processes to a steady state with the Laplace PDF, and this general capability
may explain why this distribution is often detected in confined trajectories
of single particles.

If, instead of Brownian motion, we consider the $\beta$-stable Lévy motion
with $\beta\in$ (0, 2) as a parent process [34], then the Poissonian resetting
of such random motion leads to a stationary characteristic function [24]

$\hat{p}(k,\infty|X_{r})=\frac{e^{ikX_{r}}}{1+D^{*}|k|^{\beta}/r}\,,$ (11)

corresponding to the symmetric Linnik PDF [35]. In Eq. (11) the term
$e^{ikX_{r}}$ shows that the PDF maximum is located at $X_{r}$, and $D^{*}$ is
constant. When $\beta=2$, the Linnik density reduces to the Laplace case. Both
Laplace and Linnik PDFs describe confinement with jumps (unlike the Brownian
confinement with continuous trajectories) [24], but the difference between
them is that the Linnik PDF has a heavier tail than the Laplace one [25].
Therefore, the Linnik confinement is characterized by longer jumps.

###  II.4 Advantages and disadvantages of models

As was shown above, there are different mathematical approaches leading to
confinement. Although each model may be implemented in the dependence of
physical conditions, they have their pros and cons, which are useful to list:

  * •

subordination of random processes leads to the Laplace PDF in many cases of
non-Gaussian processes, but the subordinator is too specific;

  * •

SDEs with an attractive potential are good for Gaussian confinement, but the
Laplace PDF exists for the only value of the parameter describing the
stationary PDF;

  * •

stochastic resetting has great potential for explaining non-equilibrium states
in physics, chemistry and biology, and this model provides more grounds for
understanding the diversity of confined trajectories.

The statistics of confined trajectories can play the key role in finding
physical processes responsible for the occurrence of confinement. Therefore,
discriminating tests of the random trajectories on possible PDFs are very
important for the study of confined modes. It should be noticed that
stochastic resetting is not so anomalous in living cells as it may seem. In
particular, the diffusive backward motion of paused RNA polymerases is a
diffusion process with stochastic resetting [36].

Figure 1:  Examples of random one-segment trajectories, focusing on confined
modes, in the motion of G proteins (a) and coupled receptors (b) observed in
living cells. The type of segments is identified by the DC–MSS algorithm and
statistical tests (see more details in Section III).

##  III Analysis of experimental data

To choose the most suitable model of confinement, we analyze stochastic
trajectories detected in the SPT experiment [2].

###  III.1 Classification of diffusion modes

Truthful classification of random trajectories in cells is of great importance
[37]. At the current level of development of science and technology, it builds
bridges between physics, biology, biochemistry, and biophysics necessary for
understanding how living cells function on the microscopic basis [1]. We used
the data of motion for G proteins and coupled receptors on the surface of
living cells (Figure 1). Their primary analysis can be started with
classification of the trajectories according to the standardized maximal
distance [10]. This method is useful because it allows us to evaluate quickly
the contribution of confined trajectories [14, 24]. However, it is rather
rough and does not take into account the segmentation of trajectories.
Segments can be different: immobile, confined, free (Brownian motion), and
directed (diffusion with drift). Moreover, multi-segment trajectories show
various modes in segments. Therefore, we used an accurate and computationally
efficient transient motion analysis algorithm, termed “divide-and-conquer
moment scaling spectrum” (DC–MSS) [38]. This approach includes three stages:
initial track segmentation, initial segment classification, and final
validation. The first stage calculates the maximum pairwise distance (MPD)
between particle positions within the window. The value of MPD reveals the
switches between different diffusive modes: MPD (immobile) $<$ MPD (confined)
$<$ MPD (free) $<$ MPD (directed). In the next stage the track segments
identified in the first stage are classified, using the MSS analysis of
molecule displacements. The final stage compensates for initial track
oversegmentation by testing the merger of adjacent segments. Using DC–MSS, we
have carried out the segment analysis of trajectories for G proteins and
coupled receptors in living cells. Their classification shows that there is a
significant number of confined segments (see Table 1). The difference in the
number of segments between G proteins and coupled receptors is due to a purely
experimental situation: G proteins were presented with fewer trajectories than
coupled receptors. Confined segments may have various statistics. Therefore,
the next step of our data analysis is to discriminate them in statistics.

Table 1: Classification of experimental data with random-walk segments in the
trajectories of G proteins and coupled receptors under basal conditions along
the coordinates $x$ and $y$ with the cutoff length of trajectories more and
equal to 50.

Types | G protein | Coupled
---|---|---
|  | Receptor
Coordinates | x | y | x | y
Confined segm. | 6513 | 6513 | 9637 | 9637
Normal PDF | 3373 | 3371 | 6684 | 6746
Laplace PDF | 2245 | 2285 | 2353 | 2232
Linnik PDF | 895 | 857 | 600 | 659

###  III.2 Test of confined motion

To group confined tracks of G proteins and coupled receptors by statistics, we
use a simple statistical test, based on the logarithm of the ratio of
maximized likelihoods between normal and Laplace distributions [39]. The ratio
$Q$ reads

$Q=\frac{n}{2}\ln(2)-\frac{n}{2}\ln(\pi)+n\ln(\hat{\theta})-n\ln(\hat{\sigma})+\frac{n}{2}\,,$
(12)

where the terms are dependent on the following values

$\displaystyle\hat{\theta}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}|Y_{i}-\hat{\eta}|\,,$
$\displaystyle\hat{\eta}$ $\displaystyle=$ $\displaystyle{\rm
median}\\{Y_{1},Y_{2},\dots,Y_{n}\\}\,,$ (13) $\displaystyle\hat{\sigma}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-\hat{\mu})^{2}\,,\quad\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}Y_{i}\,.$

If the statistical test shows $Q>0$, the sample satisfies the normal
distribution, but $Q<0$ suggests the Laplace distribution or a similar one to
it. This procedure took into account only segments with a length greater than
and equal to 50. Unfortunately, this test cannot refer segments with the
Linnik statistics to the Laplace case due to their close relationship.

###  III.3 Detection of Linnik confinement

To recognize segments with the Linnik PDF, we apply the test suggested by
Anderson and Arnold [40]. It allows estimating the parameter $\beta$ of the
Linnik PDF for chosen segments. This test uses the minimization of the
following objective function

$\bar{I}_{L}=\int_{0}^{\infty}\left(\bar{\eta}(z)-(1+|\xi
z|^{\beta})^{-1}\right)^{2}\,\exp(-z^{2})\,dz\,,$ (14)

where $\bar{\eta}(z)=n^{-1}\sum_{j=1}^{n}\cos(z\,y_{j})$, and
$y_{1},y_{2},\dots,y_{n}$ is a data sampling. The expression is minimized with
respect to two parameters, $\beta$ and $\xi$ (scale parameter). The presence
of the weight function $\exp(-z^{2})$ provides fast convergence of the
integral. Thus, Eq. (14) can be calculated numerically at no extra cost. The
results of applying this test for G proteins and coupled receptors are shown
in Table 1. Note that the diffusive motion of these particles differs in
statistics along $x$ and $y$. Basically, the testing supports our above
statistical analysis of confined segments. Really, the most segments with the
jump-like confinement has the index $\beta$ equal or close to 2, typical for
the Laplace statistics. Nevertheless, there is also a significant part of
segments with $\beta\approx 1.203\pm 0.002$ corresponding to the Linnik case.
They are detected in both cases, for G proteins and coupled receptors. This
shows that the effect is not an anomaly or exotic.

##  IV Resetting at power-law times

Let us mention that from the three models: subordination, system with
potential and stochastic resetting only the last one can lead to the emergence
of Linnik confinement from Brownian motion, represented certainly in many
trajectories. Moreover, the stochastic resetting has the ability to be
generalized from an exponential distribution of resetting times to a power-law
one by using the Mittag-Leffler function. The Brownian motion with Poissonian
resetting, leading to exponentially distributed times, has been studied in
great details (see [18, 21] and references therein), while the power-law
models have been explored less. For the first time stochastic resetting
following the power-law distributed times, was considered in [41, 42]. Below,
we study such a model. It includes both the resetting with power-law times and
the exponential case as a limit.

A more general resetting protocol is to take the sequence of resetting times
to be generated by a probability density function $\psi(\tau)$. This is the
so-called non-Poissonian resetting [21]. It follows that
$\Psi(\tau)=1-\int_{0}^{\tau}\psi(s)\,ds$ is the survival probability, i.e.,
the probability that no resetting has occurred up to time $\tau$. In
particular, for Poissonian resetting, the functions read
$\psi(\tau)=re^{-r\tau}$ and $\Psi(\tau)=e^{-r\tau}$. The first renewal
equation (9) can be easily generalized

$\displaystyle p_{\psi}(x,t|x_{0})$ $\displaystyle=$
$\displaystyle\Psi(t)\,G_{1}(x,t|x_{0})$ (15) $\displaystyle+$
$\displaystyle\int_{0}^{t}\psi(\tau)\,G_{1}(x,\tau|X_{r})\,d\tau\,.$

The stochastic resetting with power-law times can be implemented with help of
the Mittag-Leffler distribution. Let us recall that a statistical distribution
in terms of the Mittag-Leffler function was defined by Pillai [43].
Differentiating it, the Mittag-Leffler PDF becomes
$y^{\alpha-1}E_{\alpha,\alpha}(-y^{\alpha})$ with $0<\alpha\leq 1$, where
$E_{\alpha,\beta}(y)=\sum_{k=0}^{\infty}y^{k}/\Gamma(\alpha k+\beta)$ is the
two-parameter Mittag-Leffler function [44]. If $\alpha\to 1$, the PDF tends to
exponential. Next, in a stationary state of Eq. (15) the PDF
$p_{\alpha}(x,t|x_{0})$ takes the following form

$\displaystyle p_{\alpha}(x,\infty|x_{0})$ $\displaystyle=$ $\displaystyle
r\int_{0}^{\infty}(r\tau)^{\alpha-1}E_{\alpha,\alpha}[-(r\tau)^{\alpha}]$ (16)
$\displaystyle\times$ $\displaystyle G_{1}(x,\tau|X_{r})\,d\tau\,.$

Note that the normalization condition
$\int_{-\infty}^{\infty}p_{\alpha}(x,\infty|x_{0})\,dx=1$ is connected with
the integral
$\int_{0}^{\infty}(r\tau)^{\alpha-1}E_{\alpha,\alpha}[-(r\tau)^{\alpha}]\,d\tau=1/r$.
Moreover, in this context the function $\Psi(t)=E_{\alpha,1}[-(rt)^{\alpha}]$
tends to $e^{-rt}$ for $\alpha\to 1$ and, as expected, the relation $r=0$
shows the absence of resetting in Eq. (15). If we pass from the PDF
$G_{1}(x,t|x_{0})$ to its characteristic function $e^{-k^{2}Dt}$ through the
Fourier transform, then it is not difficult to find the characteristic
function for $p_{\alpha}(x,\infty|x_{0})$, namely

$\displaystyle\hat{p}_{\alpha}(k,\infty|x_{0})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}p_{\alpha}(x,\infty|x_{0})\,e^{ikx}\,dx$
(17) $\displaystyle=$
$\displaystyle\frac{e^{ikX_{r}}}{1+D^{\alpha}k^{2\alpha}/r^{\alpha}}\,.$

When $\alpha=1$, the expression $\hat{p}_{1}(k,\infty|x_{0})$ is the
characteristic function for Eq. (10), clearly describing the resetting with
exponentially distributed times. After the inverse Fourier transform, the
characteristic function (17) leads to the Linnik distribution generalizing the
Laplace distribution [25]. The evolution from the Brownian motion under non-
Poissonian resetting to the Linnik distribution is presented in Fig. 2. Thus,
the Linnik distribution as a stationary state in random processes under
stochastic resetting can occur not only, when the $\beta$-stable Lévy motion
for $0<\beta\leq 2$ is subject to the Poissonian resetting [24, 45], but also
when Brownian motion is under resetting at power-law times generated by the
Mittag-Leffler distribution in which $0<\alpha\leq 1$.

Figure 2:  Brownian motion under the Mittag-Leffler resetting ($r=3$, $D=2$,
$x_{0}=0.1$, $X_{r}=0.5$, $\alpha=0.75$).

The model of the Mittag-Leffler resetting is characterized by the same index
for small and large times. In Appendix A we present its generalization, using
the three-parameter Mittag-Leffler function, where the stochastic resetting
depends on two indices. This makes the model more flexible and more adequate
for the stochastic resetting behaving differently for small and large times.
Moreover, Appendix B considers the optimal search for the Mittag-Leffler
resetting.

##  V Conclusions

Single molecular imaging of G proteins and coupled receptors [46] as they
diffuse and interact on the surface of a living cell shows many diffusion
modes with frequent switching between them. Confined segments of their
trajectories are typical for the interaction of individual coupled receptors
and G proteins. Their behavior resembles the stochastic resetting. The coupled
receptors have to “find” and interact with G proteins on the cell membrane in
order to initiate and regulate intracellular processes [47]. The processes are
highly heterogeneous and complex. The interaction between G proteins and
coupled receptors occurs multiple times, and the time interval between
activation and deactivation turns out to be random. The stochastic resetting
makes this process optimal. If the stochastic resetting happens with the
exponentially distributed times, then the confined motion obeys the Laplace
distribution, but power-law distributed times generate the Linnik distribution
in confined modes. We analyzed the segmented experimental data with random
trajectories of coupled receptors and G proteins. Among the set of segments,
by using special tests, we found confined ones with the Laplace and Linnik
statistics. Our analysis to confined trajectories and fragments can serve as a
filter to identify the interaction of G proteins and coupled receptors at
power-law times.

Finally, we can summarize the following conclusions.

  * (i)

The binding of G proteins and coupled receptors [48] causes their confinement.

  * (ii)

The resetting corresponds to rebinding events. If this is the case, it means
that the confined segments actually correspond to binding-unbinding-rebinding
events where the unbinding step has not been detected.

  * (iii)

If the time spent unbound is short (the DC-MSS cannot detect segments shorter
than 20 steps), then assuming that between two consecutive binding events the
G proteins undergo Brownian motion, the Linnik/Laplace statistics might just
be the result of the mixing due to the suboptimal segmentation or the presence
of unbinding events shorter than the minimum detectable duration.

###### Acknowledgements.

A.W. thanks the support of Beethoven Grant No. DFG-NCN 2016/23/G/ST1/04083.
The authors are grateful to T. Sungkaworn and D. Calebiro for free access to
their experimental data analyzed in Section III as well as to C. Manzo for
useful discussions.

##  Appendix A Generalized Mittag-Leffler resetting

The three-parameter Mittag–Leffler function $E_{\alpha,\beta}^{\gamma}(y)$ is
a generalization of the two-parameter Mittag–Leffler function
$E_{\alpha,\beta}(y)$ mentioned in Section IV. It was first introduced by
Prabhakar [49] who defined it as

$E_{\alpha,\beta}^{\gamma}(y)=\sum_{k=0}^{\infty}\frac{(\gamma)_{k}\,y^{k}}{k!\Gamma(\alpha
k+\beta)}\,,$ (18)

where $k\in\mathbb{N}$, and $(\gamma)_{k}:=\Gamma(k+\gamma)/\Gamma(\gamma)$ is
the Pochhammer symbol. A random variable $Y$ has a generalized Mittag–Leffler
distribution [50, 51, 52], if its probability density function is

$p_{Y}(y)=\lambda^{\delta}y^{\delta\nu-1}E_{\nu,\delta\nu}^{\delta}(-\lambda
y^{\nu})$ (19)

with $y>0$, $\nu\in(0,1]$, shape $\delta\in\mathbb{R}$ and rate $\lambda>0$.
Then the first renewal equation (15) reads

$\displaystyle p_{\alpha,\delta}(x,t|x_{0})$ (20) $\displaystyle=$
$\displaystyle\left(1-(r\tau)^{\delta\alpha}E_{\alpha,\delta\alpha+1}^{\delta}[-(r\tau)^{\alpha}]\right)G_{1}(x,t|x_{0})$
$\displaystyle+$ $\displaystyle
r\int_{0}^{t}(r\tau)^{\delta\alpha-1}E_{\alpha,\delta\alpha}^{\delta}[-(r\tau)^{\alpha}]\,G_{1}(x,\tau|X_{r})\,d\tau.$

In a stationary state of Eq. (20), the PDF $p_{\alpha,\delta}(x,t|x_{0})$ is
easier to find through its characteristic function, namely

$\displaystyle\hat{p}_{\alpha,\delta}(k,\infty|x_{0})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}p_{\alpha,\delta}(x,\infty|x_{0})\,e^{ikx}\,dx$
(21) $\displaystyle=$
$\displaystyle\frac{e^{ikX_{r}}}{\left(1+D^{\alpha}k^{2\alpha}/r^{\alpha}\right)^{\delta}}\,.$

As a result, we obtain the characteristic function of the generalized Linnik
PDF [25]. When $\delta=1$, the PDF $p_{\alpha,\delta}(x,\infty|x_{0})$
simplifies to the ordinary Linnik form (see Section IV). If $\alpha=1$, then
the stationary PDF $p_{1,\delta}(x,\infty|x_{0})$ can be found explicitly in
the form (5) that is the same for the coupled Langevin equations mentioned in
Subsection II.2.

##  Appendix B Mean time to absorption

For the non-Poissonian resetting there is a convenient way to connect the
survival probability with resetting, denoted as $Q_{r}$, while one without
resetting as $Q_{1}$ [21]. It is based on the last renewal equation written as

$\displaystyle Q_{r}(x_{0},t)$ $\displaystyle=$
$\displaystyle\Psi(t)Q_{1}(x_{0},t)$ (22) $\displaystyle+$
$\displaystyle\int_{0}^{t}\psi(\tau)\,Q_{1}(X_{r},\tau)\,Q_{r}(x_{0},t-\tau)\,d\tau$

in the shorthand $Q_{r}(x_{0},t|X_{r})=Q_{r}(x_{0},t)$ and similarly for
$Q_{1}$. Taking the Laplace transform in time and setting $x_{0}=X_{r}$, we
have

$\bar{Q}_{r}(X_{r},s)=\frac{\int_{0}^{\infty}e^{-st}\,\Psi(t)\,Q_{1}(X_{r},t)\,dt}{1-\int_{0}^{\infty}e^{-st}\,\psi(t)\,Q_{1}(X_{r},t)\,dt}\,.$
(23)

Then the mean time to absorption is obtained by setting $s=0$ in Eq. (23), i.
e.

$\langle
T(X_{r})\rangle=\frac{\int_{0}^{\infty}\Psi(t)\,Q_{1}(X_{r},t)\,dt}{1-\int_{0}^{\infty}\psi(t)\,Q_{1}(X_{r},t)\,dt}\,.$
(24)

The survival probability of a diffusive particle starting from $x_{0}$ without
resetting is written as

$Q_{1}(x_{0},t)=\mathrm{erf}\left(\frac{x_{0}}{2(Dt)^{1/2}}\right)\,,$ (25)

where $\mathrm{erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-y^{2}}\,dz$ is the
Gauss error function [31]. In Eq. (24) only the functions $\Psi$ and $\psi$
depend on $r$. Therefore, by changing variables $rt\to t$ it is convenient to
collect all parameters ($r$, $X_{r}$, and $D$) in $Q_{1}(X_{r},t)$, i.e.,

$Q_{1}(X_{r},t)=\mathrm{erf}\left(\frac{\gamma}{2t^{1/2}}\right)\,,$ (26)

where $\gamma=X_{r}\sqrt{r/D}$, taking $\Psi(t,r=1)$ and $\psi(t,r=1)$ for the
fixed $r$. Consequently, Eq. (24) yields

$\langle
T(X_{r})\rangle=\frac{1}{r}\left(\frac{\int_{0}^{\infty}\Psi(t,r=1)\,Q_{1}(X_{r},t)\,dt}{1-\int_{0}^{\infty}\psi(t,r=1)\,Q_{1}(X_{r},t)\,dt}\right)\,.$
(27)

Since $\mathrm{erf}(z)\sim 2z/\sqrt{\pi}$ for $z$ tending to zero, the value
of $\langle T(X_{r})\rangle$ diverges as $r\to 0$ as $\langle
T(X_{r})\rangle\sim r^{-1/2}$. It is clear that in the absence of resetting,
the mean time is infinity. On the other hand, when $r\to\infty$, the survival
probability without resetting $Q_{1}(X_{r},t)$ is limited to one. Then
$\int_{0}^{\infty}\psi(t,r=1)\,dt\to 1$, and the asymptotic behavior of the
denominator in Eq. (27) is determined by $1-\mathrm{erf}(z)\sim
e^{-z^{2}}/(\sqrt{\pi}z)$ for $z\to\infty$. Therefore, $\langle
T(X_{r})\rangle$ also diverges as $r\to\infty$. This is not surprising. The
higher the reset rate, the shorter the time between resets to reach the
origin. These two divergences can surround one minimum of $\langle
T(X_{r})\rangle$.

Figure 3:  Mean time to absorption as a function of $r$ with the generalized
Mittag-Leffler resetting ($D=1$, $x_{0}=X_{r}=1$, $\alpha=0.85$ and
$\delta=0.5$).

However, this will be the case if the numerator of Eq. (27) is finite. Note
that $\int_{0}^{\infty}\psi(t,r=1)\,Q_{1}(X_{r},t)\,dt<1$ for $r<\infty$.
Taking $\Psi(t,r=1)=E_{\alpha,1}(-t^{\alpha})$, the improper integral in the
numerator of Eq. (27) converges, when $\alpha>1/2$. Here the convergence
problem arises at the upper limit of the integral. For $t\to\infty$ the
integrand $\Psi(t,r=1)\,Q_{1}(X_{r},t)$ depends on
$E_{\alpha,1}(-t^{\alpha})\sim t^{-\alpha}/\Gamma(1-\alpha)$, having
$0<\alpha\leq 1$, and $Q_{1}(X_{r},t)\sim\gamma/\sqrt{\pi t}$. This leads to
the convergence of the numerator of Eq. (27) just as in the case of
$1/2<\alpha\leq 1$. When $\Psi(t,r=1)$ corresponds to the generalized Mittag-
Leffler resetting, we observe the same convergence condition of this
numerator, i. e. for $1/2<\alpha\leq 1$. If the numerator of Eq. (27) is
finite, between two divergences ($r\to 0$ and $r\to\infty$) there is a single
minimum. Such an example (related to the generalized Mittag-Leffler resetting)
is shown in Figure 3. In this case the minimal mean time is achieved when the
ratio of the distance $X_{r}$ (from the resetting position to the target) to
the typical diffusion length between resets is $\gamma=1.5124\dots$. Finally,
the Mittag-Leffler (generalized and ordinary) resetting with $0<\alpha\leq
1/2$ is characterized by divergence, and the mean time for a diffusive
particle to reach the origin becomes infinite for any $r$.

## References

  * [1] C. Manzo and M.F. Garcia-Parajo, A review of progress in single particle tracking, Rep. Prog. Phys. 78, 124601 (2015).
  * [2] T. Sungkaworn, M-L. Jobin, K. Burnecki, A. Weron, M. J. Lohse, and D. Calebiro, Single-molecule imaging reveals receptor-G protein interactions at cell surface hot spots, Nature (London) 550, 543 (2017).
  * [3] Y. Lanoiselée, J. Grimes, Z. Koszegi, and D. Calebiro, Detecting transient trapping from a single trajectory: A structural approach, Entropy 23, 1044 (2021).
  * [4] G. Mu$\tilde{\rm n}$oz-Gil, G. Volpe, M.A. Garcia-March, et al., Objective comparison of methods to decode anomalous diffusion, Nat. Commun. 12, 6253 (2021).
  * [5] N. Chenouard, I. Smal, F. de Chaumont, et al., Objective comparison of particle tracking methods, Nat. Methods 11, 281 (2014).
  * [6] M.J. Saxton and K. Jacobson, Single-particle tracking: applications to membrane dynamics, Annu. Rev. Biophys. Biomol. Struct. 26, 373 (1997).
  * [7] D. Wirtz, Particle-tracking microrheology of living cells: principles and applications, Annu. Rev. Biophys. 38, 301 (2009).
  * [8] B. Brandenburg and X. Zhuang, Virus trafficking – learning from single-virus tracking, Nat. Rev. Microbiol. 5, 197 (2007).
  * [9] W.F. Marshall, A. Straight, J.F. Marko, et al., Interphase chromosomes undergo constrained diffusional motion in living cells, Curr. Biol. 7, 930 (1997).
  * [10] V. Briane, C. Kervrann and M. Vimond, Statistical analysis of particle trajectories in living cells, Phys. Rev. E 97, 062121 (2018).
  * [11] A. Weron, J. Janczura, E. Boryczka, T. Sungkaworn, and D. Calebiro, Statistical testing approach for fractional anomalous diffusion classification, Phys. Rev. E 99, 042149 (2019).
  * [12] J. Janczura, P. Kowalek, H. Loch-Olszewska, J. Szwabiński, A. Weron, Classification of particle trajectories in living cells: Machine learning versus statistical testing hypothesis for fractional anomalous diffusion, Phys. Rev. E 102, 032402 (2020).
  * [13] K. Burnecki, E. Kepten, J. Janczura, I. Bronshtein, Y. Garini, and A. Weron, Universal algorithm for identification of fractional Brownian motion. A case of telomere subdiffusion, Biophys. J. 103, 1839 (2012).
  * [14] A. Stanislavsky and A. Weron, Look at tempered subdiffusion in a conjugate map: Desire for the confinement, Entropy 22, 1317 (2020).
  * [15] F.A. Espinoza, M.J. Wester, J.M. Oliver, B.S. Wilson, N.L. Andrews, D.S. Lidke, and S.L. Steinberg, Insights into cell membrane microdomain organization from live cell single particle tracking of the IgE high affinity receptor Fc$\epsilon$RI of mast cells, Bull. Math. Biol. 74, 1857 (2012).
  * [16] A. Stanislavsky and A. Weron, Confined random motion with Laplace and Linnik statistics, J. Phys. A: Math. Theor. 54, 055009 (2021).
  * [17] G.E. Uhlenbeck and L.S. Ornstein, On the theory of Brownian motion, Phys. Rev. 36 (5), 823 (1930).
  * [18] M.R. Evans and S.N. Majumdar, Diffusion with stochastic resetting, Phys. Rev. Lett. 106, 160601 (2011).
  * [19] S.M. Ross, Introduction to Probability Models, 11th Edition (Academic Press, New York, 2014).
  * [20] M. Magdziarz and K. Taźbirski, Stochastic representation of processes with resetting, Phys. Rev. E 106, 014147 (2022).
  * [21] M.R. Evans, S.N. Majumdar, G. Schehr, Stochastic resetting and applications, J. Phys. A: Math. Theor. 53, 193001 (2020).
  * [22] A. Pal, L. Kuśmierz, and S. Reuveni, Search with home returns provides advantage under high uncertainty, Phys. Rev. Research 2, 043174 (2020).
  * [23] A. Stanislavsky and A. Weron, Subdiffusive search with home returns via stochastic resetting: A subordination scheme approach, J. Phys. A: Math. Theor. 55, 074004 (2022).
  * [24] A. Stanislavsky and A. Weron, Optimal non-Gaussian search with stochastic resetting, Phys. Rev. E 104, 014125 (2021).
  * [25] S. Kotz, T. Kozubowski, and K. Podgórski, The Laplace Distribution and Generalizations: A Revisit with Applications to Communications, Economics, Engineering, and Finance (Birkhauser, Boston, 2001).
  * [26] S. Bochner, Diffusion equation and stochastic processes, Proc. Natl. Acad. Sci. U.S.A. 35, 368 (1949).
  * [27] A. Stanislavsky and A. Weron, Accelerating and retarding anomalous diffusion: A Bernstein function approach, Phys. Rev. E 101, 052119 (2020).
  * [28] Y. Lanoiselée, A. Stanislavsky, D. Calebiro, and A. Weron, Temperature and friction fluctuations inside a harmonic potential, Phys. Rev. E 106, 064127 (2022).
  * [29] J.C. Cox, J.E. Ingersoll, and S.A. Ross, A theory of the term structure of interest rates, Econometrica 53(2), 385 (1985).
  * [30] S.P. Blomberg, S.I. Rathnayake, and C.M. Moreau, Beyond Brownian motion and the Ornstein-Uhlenbeck process: stochastic diffusion models for the evolution of quantitative characters, The American Naturalist 195(2), 145-165 (2020).
  * [31] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables (Dover Publications, New York, 1965).
  * [32] O. Kropáč, Some properties and applications of probability distributions based on MacDonald function, Aplikace matematiky 27(4), 285–302 (1982).
  * [33] R. L. Schilling, R. Song, and Z. Vondra$\check{\rm c}$ek, Bernstein Functions: Theory and Applications (de Gruyter Studies, Berlin, 2010).
  * [34] W. Feller, Introduction to Probability Theory and Its Application, Vol. II (John Wiley & Sons, New York, 1967).
  * [35] Yu.V. Linnik, Linear forms and statistical tests, Ukrain. Math. Z. 5, pp. 207-243, English Translations in Mathematical statistics and Probability, American Mathematical Society, Providence, R.I., 3, pp. 1-40, pp. 41-90 (1962).
  * [36] É. Roldán, A. Lisica, D. Sánchez-Taltavull, and S. W. Grill, Stochastic resetting in backtrack recovery by RNA polymerases, Phys. Rev. E 93, 062411 (2016).
  * [37] K. Burnecki, G. Sikora, A. Weron, M.M. Tamkun, and D. Krapf, Identifying diffusive motions in single-particle trajectories on the plasma membrane via fractional time-series models, Phys. Rev E 99, 012101 (2019).
  * [38] A.R. Vega, S.A. Freeman, S. Grinstein, and Kh. Jaqaman, Multistep track segmentation and motion classification for transient mobility analysis, Biophys. J. 114, 1018 (2018).
  * [39] D. Kundu, Discriminating between normal and Laplace distributions, In: N. Balakrishnan, H.N. Nagaraja and N. Kannan (eds) Advances in Ranking and Selection, Multiple Comparisons, and Reliability. Statistics for Industry and Technology (Boston, Birkhäuser, 2005) pp. 65–79.
  * [40] D.N. Anderson and B.C. Arnold, Linnik distributions and processes, J. Appl. Probab. 30(2), 330 (1993).
  * [41] A. Nagar and S. Gupta, Diffusion with stochastic resetting at power-law times, Phys. Rev. E 93, 060102(R) (2016).
  * [42] A.S. Bodrova, A.V. Chechkin, and I.M. Sokolov, Nonrenewal resetting of scaled Brownian motion, Phys. Rev. E 100, 012119 (2019); Scaled Brownian motion with renewal resetting, Phys. Rev. E 100, 012120 (2019).
  * [43] R.N. Pillai, On Mittag-Leffler functions and related distributions, Ann. Inst. Stat. Math. 42, 157 (1990).
  * [44] A.M. Mathai and H.J. Haubold, Special Functions for Applied Scientists (Springer, Berlin, 2008).
  * [45] C. Di Bello, A.V. Chechkin, A.K. Hartmann, Z. Palmowski, R. Metzler, Time-dependent probability density function for partial resetting dynamics, New J. Phys. 25, 082002 (2023).
  * [46] D. Calebiro, Z. Koszegi, Y. Lanoiselée, T. Miljus, and S. O’Brien, G protein-coupled receptor-G protein interactions: a single-molecule perspective, Physiol. Rev. 101(3), 857 (2021).
  * [47] B. Trzaskowski, D. Latek, S. Yuan, U. Ghoshdastider, A. Debinski, and S. Filipek, Action of Molecular Switches in GPCRs – Theoretical and Experimental Studies, Curr. Med. Chem. 19(8), 1090 (2012).
  * [48] W.I. Weis and B.K. Kobilka, The molecular basis of G protein–coupled receptor activation, Annu. Rev. Biochem. 87, 897 (2018).
  * [49] T.R. Prabhakar, A singular integral equation with a generalized Mittag–Leffler function in the kernel, Yokohama Math. J. 19, 7 (1971).
  * [50] A. Jurlewicz, K. Weron, and M. Teuerle, Generalized Mittag-Leffler relaxation: Clastering-jump continous-time random walk approach, Phys. Rev. E 78, 011103 (2008).
  * [51] K.K. Jose, P. Uma, V. Seetha Lekshmi, and H.J. Haubold, Generalized Mittag-Leffler distributions and processes for applications in astrophysics and time series modeling, In: Proceedings of the Third UN/ESA/NASA Workshop on the International Heliophysical Year 2007 and Basic Space Science, Eds. H.J. Haubold and A.M. Mathai, (Springer, Berlin Heidelberg, 2010), pp. 79-92.
  * [52] D.O. Cahoy and F. Polito, Renewal processes based on generalized Mittag–Leffler waiting times, Commun. Nonlinear Sci. Numer. Simul. 18(3), 639 (2013).