# Collective oscillations in a three-dimensional spin model with non-
reciprocal interactions

Laura Guislain and Eric Bertin  Univ. Grenoble Alpes, CNRS, LIPhy, F-38000
Grenoble, France

###### Abstract

We study the onset of collective oscillations at low temperature in a three-
dimensional spin model with non-reciprocal short-range interactions.
Performing numerical simulations of the model, the presence of a continuous
phase transition to global oscillations is confirmed by a finite-size scaling
analysis. By systematically varying the interaction range, we show that
collective oscillations in this spin model actually result from two successive
phase transitions: a mean-field phase transition over finite-size
neighborhoods, which leads to the emergence of local noisy oscillators, and a
synchronization transition of local noisy oscillators, which generates
coherent macroscopic oscillations. Using a Fokker-Planck equation under a
local mean-field approximation, we derive from the spin dynamics coupled
Langevin equations for the complex amplitudes describing noisy oscillations on
a mesoscopic scale. The phase diagram of these coupled equations is
qualitatively obtained from a fully-connected (mean-field) approximation. This
analytical approach allows us to clearly disentangle the onset of local and
global oscillations, and to identify the two main control parameters,
expressed as combinations of the microscopic parameters of the spin dynamics,
that control the phase diagram of the model.

##  1 Introduction

The onset of collective oscillations in far-from-equilibrium systems is a
ubiquitous phenomenon [1]. A first idea to model the emergence of collective
oscillations is to consider a large number of coupled oscillators, which reach
a synchronized state when interactions are strong enough. This is the basic
mechanism at play in the Kuramoto model [2, 3], which considers oscillators
with distributed frequencies. However, not all systems exhibiting collective
oscillations explicitly consist of a collection of coupled microscopic
oscillators. In many cases, individual degrees of freedom do not oscillate in
the absence of interactions. This is the case for different types of
experimental systems like biochemical clocks [4, 5, 6] or droplets in a fluid
binary mixture [7], as well as in stochastic models lacking detailed balance
like kinetic spin models [8, 9, 10, 11], socio-economic models [12, 13], or
models of population dynamics [14, 15].

For infinitely large systems, for which fluctuations can be neglected, the
onset of spontaneous oscillations can be accounted for in the frame of
dynamical system theory, and is described by a Hopf bifurcation [16]. In
contrast, oscillations occurring in mesoscopic systems are intrinsically noisy
(see, e.g., [17] in the context of biochemical clocks), notably leading to a
loss of phase coherence over a finite time [18, 19, 20, 21, 22]. Such noisy
oscillations may be phenomenologically described as a stochastic Hopf
bifurcation [23, 24]. A description rooted in the microscopic dynamics and
based on a large deviation approach has been proposed recently [25, 26] in
order to generalize the usual Landau theory of phase transitions to the case
of the onset of spontaneous oscillations. This Landau approach of the phase
transition complements earlier approaches based on nonequilibrium
thermodynamic concepts, like the presence of a singularity in the entropy
production [27, 28, 29, 30, 31, 32, 33, 34].

However, the Landau approach is restricted to mean-field systems, and a
natural question would be to try to characterize the onset of spontaneous
oscillations in finite-dimensional systems, as recently reported for instance
in experiments on cell monolayers [35, 36], in populations of biological cells
[37, 38], or in systems of active particles subjected to non-reciprocal
interactions [39, 40, 41]. This nonequilibrium phase transition has been
characterized using renormalization group methods in models of identical
coupled noisy oscillators [42, 43], confirming in particular the existence of
a continuous phase transition in three dimensions. In two dimensions, the
presence of defects has been recently shown to destroy the global oscillating
phase [44].

In this paper, we analyze a three-dimensional spin model with non-reciprocal
interactions exhibiting collective oscillations. We show using both numerical
and analytical methods that effective local oscillators emerge on a mesoscopic
scale from the microscopic dynamics through a finite-size, mean-field phase
transition. In a second stage, these local noisy oscillators may synchronize
on a global scale, leading to the onset of collective oscillations. We
describe in details the interplay of these two types of phase transitions. An
analytical approach allows us in particular to identify the few key
microscopic parameters that control the macroscopic phase diagram of the
model. The article is organized as follows. In Sec. 2, the model is introduced
and the phase transition to a globally oscillating state is characterized
numerically. The link with the synchronization transition of coupled
oscillators is explored in Sec. 3, where the mean-field phase diagram of
coupled noisy oscillators is also characterized. Finally, the influence of the
microscopic parameters of the three-dimensional spin models on the global
phase diagram obtained from coupled noisy oscillators is studied in Sec. 4.

##  2 Collective oscillations in a spin model with short-range interactions

###  2.1 Model and key observables

####  2.1.1 Definition of the model

We consider a model with $N=L^{d}$ spins $s_{i}$ and fields $h_{i}$ in
dimension $d=3$ placed on a cubic lattice, as illustrated in Fig. 1(a). Spins
and fields interact with neighboring sites. We consider spherical
neighborhoods of radius $a$, where $a$ can be non-integer [see Fig. 1(c)]. We
denote $n(a)$ the number of neighbors, and $d(i,j)$ the distance between two
sites $i$ and $j$. We consider a single spin $s_{i}$ or field $h_{i}$ flip
dynamic. The transition rates $W_{s}(i)$ and $W_{h}(i)$ to flip a spin
$s_{i}=\pm 1$ or a field $h_{i}=\pm 1$ are given by

$\displaystyle W_{s}(i)$ $\displaystyle=$ $\displaystyle[1+\exp(\Delta
E_{s}(i)/T)]^{-1},$ (1) $\displaystyle W_{h}(i)$ $\displaystyle=$
$\displaystyle[1+\exp(\Delta E_{h}(i)/T)]^{-1},$ (2)

where $T$ is the bath temperature; $\Delta E_{s}(i)$ and $\Delta E_{h}(i)$ are
the variations of interaction energies $E_{s}$ and $E_{h}$ when flipping spin
$s_{i}$ or field $h_{i}$ respectively. These interaction energies are defined
as [see Fig. 1(b)]:

$\displaystyle E_{s}$ $\displaystyle=$
$\displaystyle-\frac{1}{n(a)}\sum_{d(i,j)\leq
a}\left(\frac{J_{1}}{2}s_{i}s_{j}+s_{i}h_{j}+\frac{J_{2}}{2}h_{i}h_{j}\right),$
(3) $\displaystyle E_{h}$ $\displaystyle=$ $\displaystyle
E_{s}+\frac{\mu}{n(a)}\sum_{d(i,j)\leq a}s_{i}h_{j},$ (4)

where the sum $\sum_{d(i,j)\leq a}$ is performed over all pairs $(i,j)$
separated by a distance $d(i,j)$ less than or equal to $a$. The coupling
constants $J_{1}$ and $J_{2}$ control the ferromagnetic interactions between
the spins and between the fields. We take $J_{1}=0$ for the stochastic
simulations, and we briefly comment on the effect of $J_{1}$ in C. When
$\mu=0$, $E_{h}=E_{s}$ so that detailed balance is verified. Hence the
parameter $\mu$ controls the deviation from equilibrium, through the non-
reciprocity of interactions.

The limit where all spins and fields interact, corresponding to $n(a)=N$, is
studied in [25, 26]. In this limit, temporal oscillations appear for
$T<T_{0}$, with

$T_{0}=\frac{J_{1}+J_{2}}{2}$ (5)

provided that $\mu>\mu_{c}=1+(J_{1}+J_{2})^{2}/4$. At high $\mu$ and low $T$,
the system is then in an oscillating phase. At a fixed $\mu$, the transition
from the disordered to the oscillating phase corresponds to a Hopf bifurcation
in the deterministic limit. The phase diagram of the model is shown in Fig.
1(d). The low values of $\mu$ and $T$ correspond to a ferromagnetic phase with
two ferromagnetic points $(m_{0},h_{0}),(-m_{0},-h_{0})$.

Figure 1: (a) Schematic of the 3D model: $N$ spins and fields are placed on a
cubic lattice of size $L$. (b) Diagram of interactions between spins and
fields. (c) Representation for $a=2$ (2D section of the 3D lattice) of
interaction range, corresponding to a spherical neighborhood of radius $a$
around each site is considered in the interactions. (d) Phase diagram for
infinite range interactions ($n(a)=N\to\infty$) for $J_{1}=0$ and $J_{2}=1$,
from [25, 26]; (P) indicates the paramagnetic phase, (F) the ferromagnetic
phase and (O) the oscillating phase.

####  2.1.2 Order parameters characterizing the phase transition

The oscillating phase can be described using two order parameters. First, we
define the magnetization as

$m=N^{-1}\sum_{i=1}^{N}s_{i}$ (6)

and calculate $\langle m^{2}\rangle$ from stochastic simulations (see below).
A non-zero value of $\langle m^{2}\rangle$ corresponds to an ordered
(oscillating or ferromagnetic) phase, while $\langle m^{2}\rangle=0$ in the
limit $N\to\infty$ (or, more precisely, $\langle m^{2}\rangle\sim N^{-1}$ for
large $N$) corresponds to a paramagnetic phase. To distinguish a ferromagnetic
phase from an oscillating phase, we use the smoothed stochastic derivative of
the magnetization $\dot{m}$ introduced in [25] and defined as follows

$\dot{m}(\\{s_{i}\\},\\{h_{i}\\})=-\frac{2}{N}\sum_{i=1}^{N}s_{i}W_{s}(i).$
(7)

In the limit $N\to\infty$, the mean-square stochastic derivative of the
magnetization $\langle\dot{m}^{2}\rangle$ is non-zero only in the oscillating
phase. Unlike the magnetization, the stochastic derivative of the
magnetization is expensive to compute numerically, as all $W_{s}(i)$ must be
evaluated at each step. Thus, when extensive simulations need to be performed,
like to determine the critical temperature $T_{c}(\mu)$, it may be more
convenient to use the magnetization as order parameter, after checking from
the trajectories $m(t)$ or from $\langle\dot{m}^{2}\rangle$ that oscillations
are present.

####  2.1.3 Motivation for choosing to vary the range of interaction

First, we briefly consider the simplest version of the model where spins and
fields interact only with their nearest neighbors, corresponding to $a=1$ (and
therefore $n(a)=6$ in three dimensions) in the previously introduced model. We
have run simulations up to $L=14$ ($N=L^{3}=2744$) for different temperatures
$T$ and different $\mu$ (simulation results are presented in A.1 on Fig. 14)
for $J_{1}=0$ and $J_{2}=1$. For $\mu\geq 2$ and all temperatures, we observe
that $\langle m^{2}\rangle$ decreases in $N^{-1}$, indicating the absence of
an ordered phase. Thus, we observe that when considering nearest-neighbor
interactions ($a=1$) in this model and for the choice of $J_{1}=0$ and
$J_{2}=1$, the oscillating phase does not exist, whereas it does when
considering infinite-range interactions. The choice of $J_{1}$, $J_{2}$ is
arbitrary here, but simulations for lower $J_{2}$ values give similar results.
We note, however, that the absence of oscillations when considering nearest-
neighbor interactions does not seem to be a universal property of spin models
with non-reciprocal interactions [44], and it depends on the specific
definition of the model.

To sum up, the model exhibits collective oscillations when the interaction
range is infinite, but the oscillations disappear at very short range. We now
vary the range of interaction $a$ in order to understand how collective
oscillations appear. Intuitively, if the range of interaction is large enough,
noisy oscillators are expected locally, and then the onset of system-wide
oscillations is likely to be similar to a synchronization transition of noisy
oscillators.

###  2.2 Numerical results for $a>1$

####  2.2.1 Presence of collective oscillations

Figure 2: Evidence of collective oscillations. (a) Trajectory $m(t)$ and (b)
probability density $P_{N}(m,h)$ obtained from numerical simulations of the
three-dimensional spin model. Parameters: $J_{1}=0$, $J_{2}=1$, $\mu=2$,
$T=0.2$, $a=2$ and $L=21$. Figure 3: Evidence of a phase transition towards
collective oscillations. (a) $\langle m^{2}\rangle$ and (b)
$\langle\dot{m}^{2}\rangle$ for different temperatures $T$ for $L=21$. (c)
$\langle m^{2}\rangle$ and (d) $\langle\dot{m}^{2}\rangle$ for different
system sizes $N=L^{3}$; each color corresponds to different temperatures
[shown in (a) and (b)]. Parameters: $J_{1}=0$, $J_{2}=1$, $\mu=2$ et $a=2$.

We now consider a larger interaction range $a>1$. Before studying the
influence of varying the various parameters ($\mu$, $a$, $J_{1}$ and $J_{2}$),
we fix $J_{1}=0$, $J_{2}=1$, $\mu=2$ and $a=2$ (corresponding to a
neighborhood of $n(a)=32$ sites) and show that the model exhibits, for a
certain temperature range, collective temporal oscillations. A trajectory
$m(t)$ and the probability density $P_{N}(m,h)$ is shown in Fig. 2 for $T=0.2$
(low temperature) and a system size $L=21$ (corresponding to $N=21^{3}=9261$).
We observe spontaneous oscillations of the magnetization. Once the presence of
oscillations is confirmed by a nonzero measured value of
$\langle\dot{m}^{2}\rangle$ on a given system size, the amplitude of
oscillations is precisely characterized as a function of system size $L$ by
measuring $\langle m^{2}\rangle$, which is less costly numerically. In Figs.
3(a) and (b), we plot $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$
as a function of temperature for $L=21$. In order to verify that the observed
oscillations are not finite size effects, we study the evolution of the two
order parameters as a function of the system size $N=L^{3}$ for each
temperature [see Figs. 3(c) and (d)]. We observe that for high temperatures,
$\langle m^{2}\rangle\sim N^{-1}$ and $\langle\dot{m}^{2}\rangle\sim N^{-1}$
corresponding to a paramagnetic phase, whereas for low temperatures we have
$\langle m^{2}\rangle\neq 0$ and $\langle\dot{m}^{2}\rangle\neq 0$, which
corresponds to a phase with temporal oscillations of the magnetization.

Figure 4: (a), (b) Finite-size scaling analysis of the mean-squared
magnetization $\langle m^{2}\rangle$: (a) with the exponents of the 3D XY-
model and (b) with the exponents of the 3D Ising model. (c), (d) Finite-size
scaling analysis of the mean-squared derivative of the magnetization
$\langle\dot{m}^{2}\rangle$, with the 3D XY-model (c) or 3D Ising (d)
exponents. Parameters: $J_{1}=0$, $J_{2}=1$, $\mu=2$ et $a=2$.

To confirm the existence of a continuous phase transition and try to determine
its universality class, we use a finite-size scaling analysis of the
magnetization [45, 46]. We plot $\langle m^{2}\rangle L^{2\beta/\nu}$ as a
function of $(T-T_{c})L^{1/\nu}$ where $\beta$ and $\nu$ are the critical
exponents characterizing the magnetization and the correlation length
respectively. We then seek to obtain a collapse of the curves associated with
different sizes $N=L^{3}$. For comparison with the 3D XY-model, we use
$\beta=0.34$ and $\nu=0.66$ [47], while for the 3D Ising model we use
$\beta=0.31$ and $\nu=0.64$ [48]. We present the results in Fig. 3(a), (b).
Both sets of exponent values give a satisfactory data collapse, which confirms
the existence of a continuous phase transition to an oscillating state. We
also perform a similar finite-size scaling analysis on the derivative
$\dot{m}$ of the magnetization, by plotting $\langle\dot{m}^{2}\rangle
L^{2\beta/\nu}$ as a function of $(T-T_{c})L^{1/\nu}$ in Fig. 3(c), (d). Due
to the higher numerical cost of the evaluation of $\langle\dot{m}^{2}\rangle$
as compared to $\langle m^{2}\rangle$, less data are available, but the
resulting collapse is still reasonable. Given that exponent values of the
three-dimensional XY and Ising models are very close, our finite-size scaling
analysis does not allow us to distinguish between these two universality
classes. While in some situations the physics of the problem may help in
selecting the correct universality class, here both classes may be considered
as physically relevant: microscopically, the model is composed of Ising spins,
but at a mesoscopic scale oscillators similar to the continuous spins of the
XY model (driven at a given angular frequency) also appear, as discussed
below. Note that we find a similar behavior to previous studies of discrete
oscillators [44, 49, 50], in which simulations at larger system sizes were
performed, and where it was found that the critical exponents slowly approach
the XY universality class when increasing system size.

In summary, we observe that our model has no oscillatory phase when
considering only nearest-neighbor interactions ($a=1$), but by increasing the
interaction range $a$, an oscillatory phase appears, having critical
properties similar to those of noisy interacting oscillators [42, 43] and to
the oscillatory phase observed in similar discrete models [44, 49].

####  2.2.2 Influence of the range of interaction $a$

Figure 5: Influence of the interaction range $a$: $\langle m^{2}\rangle$ as a
function of temperature $T$ for different values of $a$. Parameters:
$J_{1}=0$, $J_{2}=1$, $\mu=2$, $L=31$.

We take $\mu=2$ and vary the interaction range $a$. While no oscillating phase
is found for $a=1$, for $a=2$ we find a transition to a spontaneously
oscillating phase at $T_{c}\approx 0.23$. In Fig. 5, we plot $\langle
m^{2}\rangle$ as a function of temperature for different values of the
interaction range $a$ for $L=31$ ($N=L^{3}=29791$). An oscillating phase
appears for $a=1.8$. The greater the interaction range, the greater the
amplitude of the oscillations and the more extended the oscillating phase
(i.e., the higher the critical temperature).

Figure 6: Evolution of the critical temperature $T_{c}$ with the different
parameters of the model ($T_{c}$ is set against the critical temperature
$T_{0}$ of the infinite range model): (a) with $n(a)$ for $\mu=2$, $J_{1}=0$
and $J_{2}=1$; (b) with $\mu$ for $a=3$, $J_{1}=0$ and $J_{2}=1$; (c) with
$J_{2}$ for $J_{1}=0$, $a=3$ and $\mu=2$. Each critical temperature is
obtained using the dependence with $N=L^{3}$ of $\langle m^{2}\rangle$ (for
sizes up to $L=41$).

To identify the dependence of the critical temperature, corresponding to the
onset of oscillations, on the interaction range, we study the dependence of
$\langle m^{2}\rangle$ with system size, for sizes up to $L=41$
($N=L^{3}=68921$). The critical temperature is displayed in Fig. 6(a) as a
function of the neighborhood size $n(a)$, and we find that

$\frac{T_{0}-T_{c}(a)}{T_{0}}\sim\frac{1}{n(a)},$ (8)

where $T_{0}=(J_{1}+J_{2})/2$ is the critical temperature of the model with
infinite interaction range.

####  2.2.3 Influence of the non-reciprocity parameter $\mu$

Figure 7: Influence of the non-reciprocity parameter $\mu$ on (a) $\langle
m^{2}\rangle$ and (b) $\langle\dot{m}^{2}\rangle$ as a function of temperature
for different values of $\mu$ indicated in the legend, for $L=14$. Same
colorcode on both panels. Parameters: $J_{1}=1$, $J_{2}=1$, $a=3$.

We now take $a=3$ and study the effect of the non-reciprocity parameter $\mu$.
In Fig. 7, we plot $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ as a
function of temperature $T$ for $L=14$ and for various high values of $\mu$
($\mu>2$). We comment on lower values of $\mu$ in A.2. We observe that
$\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ decrease with $\mu$ for
a fixed temperature and increase when $T$ decreases at a fixed $\mu$. In the
infinite interaction range limit, $m^{2}\sim 1/\mu^{2}$ while $\dot{m}^{2}\sim
1/\mu$, which explains why $\langle m^{2}\rangle$ decreases more strongly with
$\mu$ than $\langle\dot{m}^{2}\rangle$ in the three-dimensional model.

To distinguish between finite size effects ($\langle m^{2}\rangle\sim N^{-1}$)
and small values of $\langle m^{2}\rangle$ close to the transition (which
remain non-zero for $N\to\infty$), we calculate $\langle m^{2}\rangle$ for
different system sizes up to $L=31$ [see A.3] and plot the resulting phase
diagram on Fig. 6(b). As $\mu$ increases, the critical temperature decreases,
whereas in the limit of infinite interaction range, the critical temperature
is independent of $\mu$ [see the phase diagram in Fig. 1(d)]. The finite range
of interactions therefore reduces the extension of the oscillating phase in
the phase diagram $(T,\mu)$, and it is more difficult, very far from
equilibrium, to obtain collective oscillations than in the case of an infinite
interaction range.

####  2.2.4 Influence of the coupling constant $J_{2}$

Figure 8: Influence of the ferromagnetic coupling $J_{2}$ between fields:
$\langle m^{2}\rangle$ as a function of temperature for different values of
$J_{2}$. Parameters: $J_{1}=0$, $\mu=2$, $a=3$ and $L=31$.

We now briefly study the influence of the coupling constant $J_{2}$ between
field variables $h_{i}$. We take $J_{1}=0$ [there are no ferromagnetic
interactions between the spins]. We plot $\langle m^{2}\rangle$ as a function
of the reduced temperature $T/T_{0}$ (with $T_{0}=J_{2}/2$ when $J_{1}=0$) for
different values of $J_{2}$ in Fig. 8. We plot the critical temperature
$T_{c}$ for different values of $J_{2}$ in Fig. 6(c). We observe that the
stronger the ferromagnetic interactions between the fields (higher $J_{2}$),
the more the oscillating phase is extended (the higher the critical
temperature is).

####  2.2.5 Summary of numerical results

We carried out numerical simulations of the three-dimensional spin model
considering short-range interactions, varying the various parameters: $a$,
$\mu$, $J_{2}$. These simulations enabled us to conclude that:

  * •

As the interaction range increases, the oscillating phase becomes more
extended (the critical temperature $T_{c}$ approaches $T_{0}$), and the
oscillating phase appears to have a critical behavior similar to that of noisy
interacting oscillators. The oscillating phase seems to disappear when the
interaction range is very small.

  * •

The greater the non-reciprocity parameter $\mu$ (also controlling the distance
to equilibrium), the less extended the oscillating phase. The oscillating
phase seems to disappear very far from equilibrium.

  * •

The lower the ferromagnetic interactions between the fields (controlled by
$J_{2}$), the less extended the oscillating phase, and the smaller the
oscillation amplitude. The oscillating phase seems to disappear when $J_{2}$
is low.

We now seek to understand these behaviors using analytical tools, by linking
them to a model of interacting noisy oscillators.

##  3 Link with the synchronization transition of noisy oscillators

###  3.1 Effective model of interacting noisy oscillators

In order to make the link between the microscopic model of spins and fields,
and a model where oscillators would appear at the mesoscopic level, we
introduce a local averaging procedure. To do this, we divide space into
distinct boxes of size $n_{b}=(2\lfloor{a}\rfloor+1)^{3}$. We denote as
$N_{G}=N/n_{b}$ the total number of boxes. We relabel the spins and the fields
as follows: the boxes are indexed with an index $k$ and within each box, the
$n_{b}$ sites are indexed by a position $i$ with $i=0$ the central position.
We then denote $\sigma_{k}^{1}(i)$ the spins and $\sigma_{k}^{2}(i)$ the
fields. The division into boxes can be done in different ways, depending on
the choice of central spins, and we consider an arbitrary choice of central
spins. We denote $M_{k}(i)$ the average magnetization of spins separated by a
distance less than $a$ from the spin of the box $k$ at position $i$:

$M_{k}(i)=\frac{1}{n(a)}\sum_{d(({k^{\prime}},{i^{\prime}}),(k,i))\leq
a}\sigma_{k^{\prime}}^{1}(i^{\prime}),$ (9)

where $\sum_{d(({k^{\prime}},{i^{\prime}}),(k,i))\leq a}$ corresponds to a sum
over $k^{\prime}$ and $i^{\prime}$ such that the site of box $k^{\prime}$ at
position $i^{\prime}$ is located at a distance less than $a$ from the site of
box $k$ at position $i$. We denote as

$M_{k}=M_{k}(0)$ (10)

the central spin of box $k$. We also introduce $H_{k}(i)$ as

$H_{k}(i)=\frac{1}{n(a)}\sum_{d(({k^{\prime}},{i^{\prime}}),(k,i))\leq
a}\sigma_{k^{\prime}}^{2}(i^{\prime}),$ (11)

and $H_{k}=H_{k}(0)$ the local averages of the fields.

We now seek to describe the system only on the $2^{2N_{G}}$ configurations
$(\vec{M},\vec{H})=(M_{1},...,M_{N_{G}},H_{1},....,H_{N_{G}})$ instead of the
$2^{2N}$ configurations of spins and fields
$(\vec{s},\vec{h})=(s_{1},...,s_{N},h_{1},...,h_{N})$. We therefore introduce
the probability $P_{N_{G}}(\vec{M},\vec{H})$ defined as

$P_{N_{G}}(\vec{M},\vec{H})=\sum_{(\vec{s},\vec{h})\in
S(\vec{M},\vec{H})}P(\vec{s},\vec{h})$ (12)

where $S(\vec{M},\vec{H})$ represents the set of configurations with local
magnetization $M_{k}$ and local field $H_{k}$ for all groups $k\in[1,N_{G}]$.
Calculations are detailed in B.1, and we now explain the main steps of the
reasoning. To obtain an equation on $P_{N_{G}}(\vec{M},\vec{H})$, we inject
Eq. (12) into the microscopic master equation governing the evolution of the
probability of configurations $(\vec{s},\vec{h})$. We make the assumption that
$n(a)\gg 1$ (the equations obtained will therefore not be valid for small $a$)
in order to perform a perturbative expansion in $1/n(a)$ and keep only terms
up to second order in this expansion. Next, we assume that in the same box
$k$, the values of the local magnetization (and field) are all comparable,
namely $M_{k}(i)\approx M_{k}$ and $H_{k}(i)\approx H_{k}$. For the zeroth-
order term in $M_{k}(i)-M_{k}$ and $H_{k}(i)-H_{k}$, all spins (and fields) in
box $k$ have identical contributions, as the transition rates depend only on
$M_{k}(i)$ and $H_{k}(i)$. The zeroth-order term therefore gives for each
individual box $k$ a contribution similar to that of the mean-field spin
model. The sums of the first-order terms in $M_{k}(i)-M_{k}$ and in
$H_{k}(i)-H_{k}$ can be approximated by Laplacian terms, proportional to an
interaction constant $C$. Finally, as we wish to study the onset of
oscillations, we also assume that $M_{k}$ and $H_{k}$ are small, in order to
simplify the resulting expressions. We then obtain the following equation on
$P_{N_{G}}(\vec{M},\vec{H})$:

$\begin{split}\partial_{t}&P_{N_{G}}(\vec{M},\vec{H})=\\\
\sum_{k}\Bigg{[}&\partial_{M_{k}}\left(M_{k}-\tanh\left[\frac{J_{1}M_{k}+H_{k}}{T}\right]-\frac{C}{T}\Delta(J_{1}M_{k}+H_{k})\right)P_{N_{G}}(\vec{M},\vec{H})\\\
+&\partial_{H_{k}}\left(H_{k}-\tanh\left[\frac{J_{2}H_{k}+(1-\mu)M_{k}}{T}\right]-\frac{C}{T}\Delta(J_{2}H_{k}+(1-\mu)M_{k})\right)P_{N_{G}}(\vec{M},\vec{H})\\\
+&\frac{1}{n(a)}\partial_{M_{k}}^{2}P_{N_{G}}(\vec{M},\vec{H})+\frac{1}{n(a)}\partial_{H_{k}}^{2}P(\vec{M},\vec{H})\Bigg{]}.\end{split}$
(13)

where $\Delta$ is the Laplacian. Eq. (13) corresponds to a Fokker-Planck
equation on $P_{N_{G}}(\vec{M},\vec{H})$, from which we can write the
associated Langevin equations:

$\begin{split}\frac{dM_{k}}{dt}&=-M_{k}+\tanh[(J_{1}M_{k}+H_{k})/T]+\frac{C}{T}\Delta(J_{1}M_{k}+H_{k})+\sqrt{\frac{2}{n(a)}}\xi_{k}^{1},\\\
\frac{dH_{k}}{dt}&=-H_{k}+\tanh[(J_{2}H_{k}+(1-\mu)M_{k})/T]+\frac{C}{T}\Delta[J_{2}H_{k}+(1-\mu)M_{k}]+\sqrt{\frac{2}{n(a)}}\xi^{2}_{k},\end{split}$
(14)

with $\xi^{i}_{k}$ independent Gaussian white noises with correlations

$\langle\xi^{l}_{k}(t)\xi^{l^{\prime}}_{k^{\prime}}(t^{\prime})\rangle=\delta_{k,k^{\prime}}\delta_{l,l^{\prime}}\delta(t-t^{\prime}).$
(15)

These equations are similar to those describing interacting noisy oscillators
[42]. Without interaction ($C=0$), each oscillator corresponds to an
oscillator obtained from a mean-field model of $n(a)$ spins and fields, which
therefore exhibits oscillations for $T<T_{0}$ and $\mu>1+(J_{1}-J_{2})^{2}/4$.
Each oscillator is subject to a noise whose amplitude is proportional to
$1/\sqrt{n(a)}$. For $C\neq 0$, oscillators interact with their nearest
neighbors.

To understand the influence of each parameter, we derive the normal form of
these equations, obtained for $M_{k},H_{k}\ll 1$. We introduce

$Z_{k}=(-1+J_{1}/T)M_{k}+H_{k}/T+i\sqrt{v_{0}}M_{k},$ (16)

with

$v_{0}=\frac{1}{T}\sqrt{\mu-1-\frac{(J_{1}-J_{2})^{2}}{4}}.$ (17)

We change the angular reference frame to suppress the average frequency
dependence, by defining $\tilde{Z_{k}}(t)=e^{-i\sqrt{v_{0}}t}Z_{k}(t)$. We
neglect the imaginary part of the coefficients of the terms in
$|\tilde{Z_{k}}|^{2}\tilde{Z_{k}}$ and in $\Delta\tilde{Z_{k}}$ as they
contribute to a change of frequency, which is small compared to the mean
frequency $\sqrt{v_{0}}$. We find the following Langevin equations on
$\tilde{Z_{k}}$:

$\frac{d\tilde{Z}_{k}}{dt}=(T_{0}/T-1)\tilde{Z}_{k}-U|\tilde{Z}_{k}|^{2}\tilde{Z}_{k}+\frac{T_{0}C}{T}\Delta\tilde{Z}_{k}+\sqrt{\frac{\alpha}{n(a)}}\left[\eta_{k}^{R}+i\eta_{k}^{I}\right]$
(18)

with $T_{0}$ defined in Eq. (5), and

$\displaystyle\alpha=\frac{1}{T^{2}}\left[(J_{1}-T)^{2}+1+(\mu-1)-\frac{(J_{1}-J_{2})^{2}}{4}\right],$
(19) $\displaystyle
U(T=T_{0})=\frac{(J_{1}+(\mu-1)J_{2})(\mu-1+J_{1}J_{2})}{(J_{1}+J_{2})[4(\mu-1)-(J_{1}-J_{2})^{2}]}.$
(20)

In Eq. (18), $\eta_{k}^{R}$ and $\eta_{k}^{I}$ are independent Gaussian white
noises, satisfying
$\langle\eta_{k}^{R}(t)\eta_{k^{\prime}}^{R}(t^{\prime})\rangle=\langle\eta_{k}^{I}(t)\eta_{k^{\prime}}^{I}(t^{\prime})\rangle=\delta_{k,k^{\prime}}\delta(t-t^{\prime})$
and $\langle\eta_{k}^{R}\eta_{k}^{I}\rangle=0$. To reduce the number of
parameters involved in these equations, we reabsorb the parameter $U$ into a
new variable $z_{k}$ defined as

$z_{k}=\sqrt{U}\tilde{Z}_{k}.$ (21)

We further introduce the two parameters:

$D=\frac{U\alpha}{n(a)}$ (22)

and

$\gamma=T_{0}/T,$ (23)

so that Eq. (18) becomes

$\frac{dz_{k}}{dt}=(\gamma-1)z_{k}-|z_{k}|^{2}z_{k}+\gamma C\Delta
z_{k}+\sqrt{D}\left[\eta_{k}^{R}+i\eta_{k}^{I}\right].$ (24)

The parameter $C$ can be estimated from the microscopic model [see B.1], and
is found to depend only weakly on the value of interaction range $a$, at odds
with the parameters $D$ and $\gamma$.

We now discuss the properties of Eq. (24) for generic values of the
coefficients, temporarily ignoring the relation to the microscopic model. The
link between this equation and the microscopic model will then be discussed in
section 4. Eq. (24) corresponds to a stochastic Ginzburg-Landau equation [43,
51], which describes interacting noisy oscillators. The parameter $D$ controls
the noise amplitude. The parameter $\gamma$ controls the amplitude of the
oscillations. Indeed, when the oscillators are neither noisy ($D=0$) nor
interacting ($C=0$), then oscillations of amplitude $z_{k}=\sqrt{\gamma-1}$
appear for $\gamma\geq 1$ on each oscillator. Finally, the parameter $C$
controls the amplitude of the interactions. For $C=0$, the oscillators do not
interact and are thus unsynchronized. In models of noisy interacting
oscillators, the noise is often neglected on the amplitude of the
oscillations, allowing the equations to be reduced to a phase equation for
constant-amplitude oscillators [2]. However, in the model presented here, the
oscillation amplitude may be small compared to the noise amplitude
($\gamma-1\ll D$), so that the assumption of constant-amplitude oscillators is
no longer valid.

###  3.2 Phase diagram of the synchronization of noisy oscillators in the mean
field limit

####  3.2.1 Mean field limit

Our goal is now to study the mean-field phase diagram of the coupled noisy
oscillators described by Eq. (24). While the parameters $D$ and $\gamma$ have
significant variations as a function of microscopic parameters of the three-
dimensional spin model (like the interaction range $a$), the parameter $C$
remains almost constant with a value $C\approx 0.025$ [see B.1]. We thus take
$C$ as a constant in the following, and focus on the influence of the
parameters $D$ and $\gamma$ on the behavior of Eq. (24). As a mean-field
approximation to the three-dimensional model described in Eq. (24), we
consider the solvable case where all the individual noisy oscillators interact
with each other in a mean-field way:

$\Delta z_{k}\to 6(\overline{z}-z_{k})$ (25)

where $6$ corresponds to the number of nearest neighbors in a cubic lattice in
dimension $d=3$ and

$\overline{z}=\frac{1}{N_{G}}\sum_{k=1}^{N_{G}}z_{k}$ (26)

is the mean value of $z_{k}$ over all boxes. To lighten notations, we
introduce $c=6C\approx 0.15$. In polar coordinates, we write:

$\overline{z}=Re^{i\Theta},\quad z=re^{i\theta}.$ (27)

We take the limit $N_{G}\to\infty$, and since we are considering the limit
where all boxes interact, we introduce $P_{1}(x,y)$, the distribution for a
single box, such that we assume

$P_{N_{G}}(\\{r_{1},\theta_{1},...,r_{N_{G}},\theta_{N_{G}}\\})=\prod_{i=1}^{N}P_{1}(r_{i},\theta_{i}).$
(28)

We can then obtain $P_{1}(r,\theta)$ as [see B.2 for detailed calculations]:

$P_{1}(r,\theta)=p_{0}\exp{\Bigg{[}-\frac{2}{D}\left(\frac{1}{4}r^{4}-\frac{(\gamma-1)}{2}r^{2}+\frac{\gamma
c}{4}\left(R^{2}+r^{2}-2Rr\cos(\Theta-\theta)\right)\right)\Bigg{]}},$ (29)

where $p_{0}$ is a normalization constant such that

$\int drrd\theta P_{1}(r,\theta)=1.$ (30)

The values of $R$ and $\Theta$ are determined from the following self-
consistent equation:

$Re^{i\Theta}=\frac{\int_{0}^{\infty}dr\int_{-\pi}^{\pi}d\theta
r^{2}e^{i\theta}\exp\left[{-\frac{2}{D}\left(\frac{1}{4}r^{4}-\frac{\gamma-1}{2}r^{2}+\frac{\gamma
c}{4}\left(r^{2}+R^{2}-2Rr\cos(\Theta-\theta)\right)\right)}\right]}{\int_{0}^{\infty}dr\int_{-\pi}^{\pi}d\theta
r\exp\left[{-\frac{2}{D}\left(\frac{1}{4}r^{4}-\frac{\gamma-1}{2}r^{2}+\frac{\gamma
c}{4}\left(r^{2}+R^{2}-2Rr\cos(\Theta-\theta)\right)\right)}\right]}$ (31)

which comes from the definition: $\overline{z}=Re^{i\Theta}=\int rdrd\theta
re^{i\theta}P_{1}(r,\theta)$. A change of variable
$\theta^{\prime}=\theta-\Theta$ in the integrals shows that Eq. (31) is
independent of $\Theta$. The global phase $\Theta$ is therefore arbitrary, and
only the relative phase $\theta-\Theta$ indicates the presence or absence of
synchronization.

Figure 9: Phase diagram of the noisy coupled oscillator model in the mean-
field limit. Amplitudes $R$ (left panel), $r_{max}$ (middle panel) and
$\sqrt{\langle r_{k}^{2}\rangle}$ (right panel) in the $(\gamma,D)$ plane
obtained numerically from Eqs. (29) and (31).  Figure 10: Top: Phase diagram
of the noisy coupled oscillator model in the mean-field limit. The blue curve
corresponds to $\gamma_{c}(D)$ [see Eq. (32)], the red curve corresponds to
$\gamma_{loc}$ [see Eq. (34)]. Three different phases are observed: a phase
with global oscillations where $r_{max}\neq 0$ and $R\neq 0$, a phase with
unsynchronized local oscillations where $r_{max}\neq 0$ and $R=0$, and a phase
with noise $r_{max}=0$ and $R=0$. Bottom: Examples of probability
distributions $P_{1}(r,\theta)$, plotted with coordinates
$x=r\cos(\theta-\Theta)$ and $y=r\sin(\theta-\Theta)$, for each region of the
phase diagram, marked by (A), (B) and (C) in the top panel.

####  3.2.2 Presence of collective oscillations

Since $R$ corresponds to the amplitude of $\overline{z}$, there is a
transition to (partially) synchronized oscillations (with frequency
$\sqrt{v_{0}}$) when Eq. (31) admits a solution with $R\neq 0$. In Fig. 9, we
show with a colorcode in the $(\gamma,D)$-plane the value of $R$ ($0\leq R\leq
1$), obtained from Eq. (31). We obtain that for high values of noise $D$ and
low values of $\gamma$, there are no macroscopic oscillations, $R=0$, while
such oscillations are present for low values of noise. Expanding Eq. (31) for
$R\ll 1$, we find that a non-zero $R$ solution appears when the pair
$(\gamma,D)$ verifies the following equation:

$\begin{split}1=\frac{\gamma
c}{D}\Bigg{[}\gamma(1-c)-1+\frac{\sqrt{\frac{2D}{\pi}}\exp\left[{-\frac{(\gamma(1-c)-1)^{2}}{2D}}\right]}{\text{erfc}\left(\frac{-\gamma(1-c)+1}{\sqrt{2D}}\right)}\Bigg{]}\end{split}$
(32)

where $\text{erfc}(x)=1-\text{erf}(x)$, and $\text{erf}(x)$ is the error
function. This equation can be solved numerically to obtain the critical
parameters $\gamma_{c}(D)$ (or $D_{c}(\gamma)$) where a solution $R\neq 0$
appears. The results of which are shown in Fig. 10. Taking the limit $D\ll 1$
and $\gamma\sim 1$, we obtain

$\gamma_{c}=1+\frac{2D}{c}.$ (33)

In the low-noise limit ($D\to 0$), we obtain $\gamma_{c}=1$ when $c\neq 0$: as
soon as oscillations appear locally ($\gamma=1$), they are synchronized.

Models of oscillator synchronization transitions often start with oscillators
that oscillate at a given frequency and amplitude, usually with a small noise
compared to the oscillation amplitude. These models can then be used to study
whether the oscillators are synchronized or not [2] by neglecting the noise on
the oscillator amplitude. In the model introduced here, the oscillation
amplitude can become negligible compared to the noise, which adds a third
phase to the phase diagram, with neither collective nor local oscillations. We
now turn to the condition for oscillations to occur locally, but to be
desynchronized on a global scale.

####  3.2.3 Presence of local oscillators

We try here to determine the parameter values for which oscillators exist
locally, but are unsynchronized on a global scale. First, we need to choose
the criterion for the presence of local oscillations. Oscillations are
considered to exist when the maximum of $P_{1}(r,\theta)$ occurs at a non-zero
value of $r$, denoted $r_{max}$. In Fig. 9, we show in color the value of
$r_{max}$ in the plane $(\gamma,D)$ obtained numerically from equations (29)
and (31). And from Eq. (29), we obtain that local oscillations appear for
$\gamma=\gamma_{loc}$, with

$\gamma_{loc}=(1-c)^{-1}.$ (34)

For low values of $\gamma$ (i.e., $\gamma<\gamma_{loc}$), there are no local
oscillations. The interaction amplitude $c$ is therefore the quantity that
controls the extent to which oscillators can exhibit local oscillations. As
$c$ is fixed in the model, we cannot control the onset of local oscillations.
The $\gamma_{loc}$ curve, which is independent of $D$, is plotted in Fig. 10
in the ($\gamma,D$) plane.

####  3.2.4 Summary of the phase diagram

To sum up, we have obtained a phase diagram in the space $(\gamma,D)$ shown in
Fig. 10. We observe a phase where there are global oscillations ($R\neq 0$)
for small values of $D$ and give an example of the probability density
$P_{1}(x,y)$ in Fig. 10(C). The average phase of each oscillator is
$\theta_{k}=\Theta$, and there is little noise in both amplitude and phase.
There is a second phase, for low values of $\gamma$ and high values of $D$,
where there are no oscillators locally, and the $z_{k}$ variables are only
subject to noise. An example of $P_{1}(x,y)$ in this phase is given in Fig.
10(A). Finally, there is a phase, for high values of noise $D$ and high values
of $\gamma$, where effective local oscillators emerge (i.e., $P_{1}(x,y)$ is
maximal for a non-zero value of $r=\sqrt{x^{2}+y^{2}}$) but are not
synchronized. In this case, the phase $\theta$ is uniformly distributed and
therefore $R=0$. An example of $P_{1}(x,y)$ in this phase is shown in Fig.
10(B).

####  3.2.5 Influence of noise on local oscillation detection

It is sometimes interesting to look only at moments of the probability density
$P_{1}(r,\theta)$, such as the second moment $\langle r_{k}^{2}\rangle$. In
this case, the condition for the existence of local oscillations is different,
as noise can hide the presence of very low-amplitude oscillations. Indeed,
when $R=0$, we obtain from Eq. (29),

$\langle
r_{k}^{2}\rangle=\gamma(1-c)-1+\sqrt{\frac{2D}{\pi}}\frac{\exp\left[{-\frac{(\gamma(1-c)-1)^{2}}{2D}}\right]}{1+\text{erf}\left(\frac{\gamma(1-c)-1}{\sqrt{2D}}\right)}.$
(35)

We therefore find that if

$\gamma(1-c)-1\ll D^{1/2},$ (36)

fluctuations dominate and even if $r_{max}\neq 0$, it is difficult to
distinguish this case from the one where $P_{1}(r,\theta)$ is maximal in
$r=0$. We calculate $\langle r_{k}^{2}\rangle$ from Eq. (29) and we represent
$\sqrt{\langle r_{k}^{2}\rangle}$ in the $(\gamma,D)$-plane in Fig. 9. The
second moment $\langle r_{k}^{2}\rangle$ is not sufficient to determine
whether local oscillations are present or not.

###  3.3 Comparison with 3D lattice simulations of interacting noisy
oscillators

Figure 11: Numerical study of the phase diagram of the three-dimensional noisy
coupled oscillator model defined by Eq. (24), simulated on a cubic lattice of
size $N=10^{3}$. Amplitudes $R$ and $\sqrt{\langle r_{k}^{2}\rangle}$ are
plotted in the ($\gamma,D)$ plane, to detect the presence of global (left
panel) or local (right panel) oscillations.

In Sec. 3.2, we considered a mean-field limit where all boxes interact with
each other. We now compare these mean-field results with the behavior of the
three-dimensional model of noisy oscillators given in Eq. (24), which was
derived from the three-dimensional spin model. We therefore run simulations of
the Langevin equation Eq. (24) on a 3D lattice to compare with the mean-field
results. We take $N=1000$ oscillators in the numerical simulations. The
amplitude $R$ of $\overline{z}$ and the second moment $\sqrt{\langle
r_{k}^{2}\rangle}$ are plotted in the $(\gamma,D)$ plane in Fig. 11, and we
observe a qualitative agreement with results obtained in the mean-field limit
[Fig. 9]. The phase without global oscillations corresponds to high $D$ and
low $\gamma$, the synchronized phase corresponds to low $D$ and high $\gamma$,
and there is a phase with local but unsynchronized oscillators for high
$\gamma$ and high $D$ (which can be confirmed by looking at the probability
density $P_{1}(x,y)$ not shown here, insofar as $R$ and $\sqrt{\langle
r_{k}^{2}\rangle}$ are not enough to identify this phase).

##  4 Comparison between the spin model and the interacting noisy oscillator
model

We have obtained the phase diagram for interacting noisy oscillators, where
the two main control parameters are the noise amplitude $D$ and the parameter
$\gamma$ controlling the amplitude of oscillations in the absence of
interactions [see Eqs. (22) and (23) for the definition of $D$ and $\gamma$].
We now come back to the three-dimensional spin model, and seek to explain
qualitatively the results observed in the simulations by comparing numerical
results with the analytical results of the mean-field model of noisy
oscillators, taking into account the expression of the parameters $D$ and
$\gamma$ in terms of the microscopic parameters of the three-dimensional spin
model.

###  4.1 Influence of interaction range $a$

The interaction range $a$ enters only the definition of the parameter $D$,
with $D\propto n(a)^{-1}$ where $n(a)$ is the number of sites with which a
spin or a field interacts. We obtained that $\gamma_{c}-1\sim 2D/c$ in the
mean-field limit of interacting noisy oscillators [see Eq. (33)]. Since
$\gamma=T_{0}/T$, the critical temperature depends on $a$ as

$\frac{T_{0}-T_{c}}{T_{0}}\sim\frac{1}{n(a)}.$ (37)

This is the dependence obtained for stochastic simulations of the three-
dimensional spin model [see Fig. 6(a)]. Qualitatively, when the neighborhood
size is large, low-noise oscillators appear locally and synchronize, whereas
for small neighborhoods, the noise is too high and oscillations are destroyed.

###  4.2 Influence of the non-reciprocity parameter $\mu$

Figure 12: Influence of $\mu$ on parameters (a) $U$ controlling oscillation
amplitude and (b) $\alpha U$ controlling noise amplitude ($D\propto\alpha U$).
(c) Comparison of critical temperatures $T_{c}$ normalized by $T_{0}$ as a
function of $\mu$ obtained simulations of the three-dimensional spin model
[blue dots, from simulations reported in Fig. (7)], with the noisy coupled
oscillator model in the mean-field approximation [red curve, solution of Eq.
(33)]. Parameters: $J_{1}=0$, $J_{2}=1$ [for (a), (b), (c)] and $a=3$ [for
(c)].

Before comparing critical temperatures, we note a strong dependence of $m$ and
$\dot{m}$ with $\mu$. The equation that controls the behavior of interacting
noisy oscillators is Eq. (24). To arrive at this equation, we started with the
variables $M_{k}$ and $H_{k}$, corresponding to the local magnetization and
the local mean field, and made the following change of variable:

$(-1+J_{1}/T)M_{k}+H_{k}/T+i\sqrt{v_{0}}M_{k}=e^{i\sqrt{v_{0}}t}z_{k}/\sqrt{U},$
(38)

with $v_{0}=\sqrt{\mu-\mu_{0}}/T$, where
$\mu_{0}(T_{0})=1+(J_{1}-J_{2})^{2}/4$. In terms of the variable $z_{k}$ (or
$\overline{z}$) we have

$M_{k}\sim z_{k}/\sqrt{v_{0}U},\quad m\sim\overline{z}/\sqrt{v_{0}U}.$ (39)

The variable $\dot{m}$ defined in Eq. (7), and corresponding to the stochastic
derivative of the magnetization, can be rewritten as

$\dot{m}=\frac{1}{N_{G}}\sum_{k=1}^{N_{G}}\sum_{i=1}^{n_{b}}\frac{-2\sigma_{k}^{l}(i)}{n_{b}}\frac{1}{1+\exp\left[-2\sigma_{k}^{l}(i)[J_{1}M_{k}(i)+H_{k}(i)]/T\right]}.$
(40)

Keeping only the lowest order in $M_{k}(i)-M_{k}$ and $H_{k}(i)-H_{k}$, we
then obtain the following approximation for $\dot{m}$:

$\dot{m}=\frac{1}{N_{G}}\sum_{k=1}^{N_{G}}\left[-M_{k}+\tanh\left(\frac{J_{1}M_{k}+H_{k}}{T}\right)\right].$
(41)

At lowest order in $M_{k}$ and $H_{k}$, we get
$\dot{m}\approx(-1+J_{1}/T)m+h/T$, which leads to the following scaling law
[see Eq. (38)]

$\dot{m}\sim\overline{z}/\sqrt{U}.$ (42)

The dependence of $U$ on $\mu$ is given in Eq. (20)(a). We plot $U$ as a
function of $\mu$ in Fig. 12(a) for $J_{1}=0$ and $J_{2}=1$, and obtain that
$U$ increases with $\mu$. As $m\sim z/\sqrt{v_{0}U}$ we therefore obtain that
$m\sim z/\mu^{2}$ (neglecting $\mu_{0}$ with respect to $\mu$), and
$\dot{m}\sim z/\mu$. We can therefore explain the trend observed in Fig. 7:
$\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ decrease with $\mu$,
but the decrease in $\langle m^{2}\rangle$ is much stronger. This explains
why, for high values of $\mu$, $m$ and $\dot{m}$ are very small, making it
difficult to observe oscillations numerically for finite system sizes.

For values of $\mu$ up to $\mu\approx 4$, we were still able to measure
critical temperatures [see Fig. 7]. We obtain that $T_{c}$ decreases when $U$
increases. To understand this behavior, we determine the dependence of the
coefficient $D$ on $\mu$. We have

$D\propto\alpha U,$ (43)

where $\alpha$ and $U$ are given in Eqs. (19) and (20). The dependence of
$\alpha U$ is shown in Fig. 12(b) for $J_{1}=0$ and $J_{2}=1$. Except for low
values of $\mu$ (where the oscillating phase no longer exists), $\alpha U$ and
therefore the noise parameter $D$ increase with $\mu$. According to the phase
diagram shown in Fig. 10, as the noise parameter $D$ increases, the transition
takes place at higher $\gamma=T_{0}/T$ and therefore the critical temperature
$T_{c}$ decreases. In Fig. 12(c), we plot the evolution of $T_{c}$ for
different values of $\mu$. The points correspond to stochastic simulations of
the three-dimensional spin model [results shown in Fig. 7] and the curve
corresponds to $T_{c}(\mu)$ obtained for the noisy oscillator model with the
mean-field couplings [see Eq. (33)]. A quantitative difference between the two
models is observed, but the overall trend is similar.

###  4.3 Influence de $J_{2}$

Figure 13: Influence of $J_{2}$ on parameters (a) $U$ controlling oscillation
amplitude and (b) $\alpha U$ controlling noise amplitude ($D\propto\alpha U)$.
(c) Comparison of critical temperatures $T_{c}$ normalized by $T_{0}$ as a
function of $J_{2}$ obtained from simulations of the three-dimensional spin
model [blue dots, from simulations reported in Fig. (8)] with the noisy
coupled oscillator model in the mean-field approximation [red curve, solution
of Eq. (33)]. Parameters: $\mu=2$, $J_{1}=0$ [panels (a), (b) and (c)] and
$a=3$ [panel (c)].

We finally study the influence of the amplitude of the ferromagnetic
interactions between the fields $J_{2}$ (for $J_{1}=0$). For the spin model
simulations, we observed two main features: $\langle m^{2}\rangle$ increases
with $J_{2}$ and the transition temperature also increases with $J_{2}$. For
the model of noisy oscillators interacting with mean-field couplings, the
coefficient $U$ depends only weakly on $J_{2}$ [see Fig. 13(a)], so that the
value of $\dot{m}$ changes only slightly if the oscillators are completely
synchronized. However, $m\sim(J_{1}+J_{2})\overline{z}/\sqrt{U}$ as
$v_{0}\propto 1/(J_{1}+J_{2})^{2}$ and therefore $m$ depends strongly on the
value of $J_{2}$. Furthermore, we find that the coefficient $D\propto\alpha U$
increases as $J_{2}$ increases [see Fig. 13(b)], which implies that the
transition temperature increases as $J_{2}$ increases. In Fig. 13(c), we plot
the scaled critical temperature $T_{c}/T_{0}$ obtained from Eq. (32) for
different values of $J_{2}$ (with $J_{1}=0$, $a=3$ and $\mu=2$) and compare
with the critical temperatures obtained from simulations of the three-
dimensional spin model with the same parameters. A qualitative agreement
between the two models is observed.

The influence of the ferromagnetic interaction $J_{1}$ between the spins is
briefly discussed in C. We obtain that the coupling constant $J_{1}$ has a
similar role as the coupling constant $J_{2}$.

##  5 Conclusion

We have considered a stochastic model consisting of spins and fields arranged
on a three-dimensional cubic lattice, with interactions on a spherical
neighborhood. We studied the influence of the interaction range using
numerical simulations and analytical tools. We have shown that decreasing the
range of ferromagnetic couplings, or strongly increasing the non-reciprocity
of on-site interactions between spins and fields, hinders the onset of
collective oscillations. To understand the phase diagram of the model in terms
of microscopic parameters, we linked the three-dimensional spin model, in
which individual degrees of freedom do not oscillate in the absence of
interactions, to an effective model of coupled noisy oscillators on a cubic
lattice, which we analyzed numerically as well as analytically within a mean-
field approximation. This mapping allowed us to propose the following
qualitative interpretation of the transition scenario in three dimensions. At
the local level, noisy oscillations emerge (similar to the mean-field case,
but for a finite-size system), whose amplitude and fluctuations are related to
the interaction range. Subsequently, these local oscillators may synchronize
at system level. The coherence of the emerging picture a posteriori justifies
the approximations made in the analytical derivations.

As for future work, it would be useful to study such finite-dimensional
systems more in depth, notably as a function of space dimension. One may
expect that the onset of local noisy oscillations combined with a large-scale
synchronization transition may be a generic two-stage scenario for the onset
of collective oscillations in finite-dimensional spin system, that may be
worth confirming in different models and in different dimensions. However,
this scenario is likely to be modified in low dimensions. In particular, for
two-dimensional systems, the oscillating phase has been shown to be destroyed
by the propagation of defects [44]. In addition, it would also be of interest
to develop a renormalization group approach starting from the spin level of
description, thereby extending the results of [42, 43] based on the
synchronization of oscillators.

##  Appendix A Simulation methods

###  A.1 Simulations with nearest-neighbor interactions

Figure 14: Order parameters (a) $\langle m^{2}\rangle$ and (b)
$\langle\dot{m}^{2}\rangle$ obtained from numerical simulations of the spin
model on a cubic lattice in three dimensions with nearest-neighbor
interactions ($a=1$) for different $\mu$ (different colors indicated in the
legend) and different temperature $T$ (different transparencies,
$T\in[0.01,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55]$ from lightest
to darkest). Parameters: $J_{1}=0$ and $J_{2}=1$.

Here we present results of simulations of the three-dimensional spin model on
a cubic lattice, with nearest-neighbor interactions ($a=1$), to illustrate the
results presented in Sec. 2.1.3, indicating that for $a=1$ we observe no phase
transition. In Fig. 14, we plot $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ as a function of system size $N$ for different
values of the parameters $T$ and $\mu$. As a reminder, in the presence of
oscillations, $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ go to a
nonzero limit when $N\to\infty$. We plot them for different values of $\mu$
(different colors) and for different temperatures $T$ (different
transparency). For all parameters explored, $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ decrease in $1/N$, indicating that there is no
oscillating phase.

###  A.2 Comment on the low-temperature phase close to equilibrium

In this article, we do not study lower values of $\mu$ for which a
ferromagnetic phase would exist in the long-range interaction limit
($a\to\infty$). According to [44], the presence or absence of the finite-
dimensional ferromagnetic phase in a model of non-equilibrium spins with
ferromagnetic interactions is not clearly identified. At finite sizes, the
system may remain locked in an ordered phase. Simulations of very large system
sizes would be needed to see that this phase is not stable. Our simulations do
not provide a clear-cut answer to this question, so that we do not present
numerical results for low $\mu$ values of the three-dimensional model.

###  A.3 Identifying the transition temperature as $\mu$ varies

Figure 15: Numerical simulations of the spin model on a three-dimensional
lattice with interaction range $a=3$ for different values of $\mu$ and
temperatures $T$. (a) $\langle m^{2}\rangle$ as a function of temperature for
different values of $\mu$ (different colors shown in legend) for different
system sizes $N=L^{3}$, $L\in[7,10,14,21,31]$ from most transparent to least
transparent. (b)-(f) $\langle m^{2}\rangle$ as a function of system size $N$
for different temperatures
$T\in[0.01,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55]$ from
yellow/light to violet/dark. Each panel corresponds to a different $\mu$, as
indicated. Parameters: $J_{1}=0$ and $J_{2}=1$.

We explain here how the transition temperature $T_{c}(\mu)$ was obtained. In
Fig. 7, we presented $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ as
a function of temperature for $L=14$. The magnetization is less numerically
expensive than the stochastic derivative of the magnetization, hence to
determine the critical temperature $T_{c}(\mu)$ we use only $\langle
m^{2}\rangle$. In Fig. 15(a), we plot $\langle m^{2}\rangle$ as a function of
temperature $T$ for different values of $\mu$ (different colors) and different
system sizes $N$ (different transparency). To obtain the critical temperature,
we plot $\langle m^{2}\rangle$ as a function of system size $N$ for different
temperatures $T$ at fixed $\mu$ [see Figs. 15(b)-(f)]. When $\langle
m^{2}\rangle$ is constant with $N$, oscillations are presence (since we have
seen, for lower system sizes, that $\langle\dot{m}^{2}\rangle\neq 0$) whereas
$\langle m^{2}\rangle\propto 1/N$ corresponds to the paramagnetic phase. The
critical temperature $T_{c}$ is the temperature separating these two regimes.

##  Appendix B Analytical methods

###  B.1 Details of the derivation of a Fokker-Planck equation

In this appendix, we detail the reasoning and calculations that led to Eq.
(13), a Fokker-Planck equation on the probability density of local, coarse-
grained magnetizations $M_{k}$ and fields $H_{k}$. The reader is referred to
Sec. 3.1 for notation definitions. We change the notation of the transition
rates to take into account the position in the box, $W_{s}(i)$ becomes

$w_{1}(k,i,\sigma_{k}^{1}(i))=\left[1+e^{2\sigma_{k}^{1}(i)[J_{1}M_{k}(i)+H_{k}(i)]/T}\right]^{-1}$
(44)

and $W_{h}(i)$ turns into

$w_{2}(k,i,\sigma^{2}_{k}(i))=\left[1+e^{2\sigma_{k}^{2}(i)[J_{2}H_{k}(i)+(1-\mu)M_{k}(i)]/T}\right]^{-1}.$
(45)

We aim to describe the system only in the $2^{2N_{G}}$ configurations of
$M_{k}$ and $H_{k}$ instead of the $2^{2N}$ configurations of the spins and
fields $\sigma_{k}^{l}$. We write

$\vec{\Sigma}=(M_{1},...,M_{N_{G}},H_{1},...,H_{N_{G}})$ (46)

and $\vec{\sigma}$ is a configuration of all spins $\sigma_{k}^{1}(i)$ and
fields $\sigma_{k}^{2}(i)$. We introduce

$P_{N_{G}}(\vec{\Sigma})=\sum_{\vec{\sigma}\in\mathcal{S}(\vec{\Sigma})}P(\vec{\sigma})$
(47)

where $\mathcal{S}(\vec{\Sigma})$ represents the set of configurations with
local magnetization $M_{k}$ and local field $H_{k}$ for each group
$k\in[1,N_{G}]$. The master equation for microscopic $\vec{\sigma}$
configurations reads as follows

$d_{t}P(\vec{\sigma})=\sum_{\vec{\sigma}^{\prime}(\neq\vec{\sigma})}[W(\vec{\sigma}|\vec{\sigma}^{\prime})P(\vec{\sigma}^{\prime})-W(\vec{\sigma}^{\prime}|\vec{\sigma})P(\vec{\sigma})].$
(48)

Using the new notations (44) and (45) for the transition rates, Eq. (48) can
be rewritten as

$d_{t}P(\vec{\sigma})=\sum_{l=1,2}\sum_{k=1}^{N_{G}}\sum_{i=0}^{n_{b}-1}\left[w_{l}\left(k,i,-\sigma^{l}_{k}(i)\right)P(\vec{\sigma}^{l}_{k}(i))-w_{l}(k,i,\sigma^{l}_{k}(i))P(\vec{\sigma})\right],$
(49)

where $\vec{\sigma}^{l}_{k}(i)$ is the configuration obtained from
$\vec{\sigma}$ after flipping the spin ($l=1$) or field ($l=2$) at position
$i$ in box $k$. We introduce a sum over $\sigma=\pm 1$ to replace the spin
$\sigma_{k}^{l}(i)$ in the transition rates, also introducing the
corresponding Kronecker deltas:

$d_{t}P(\vec{\sigma})=\sum_{l=1,2}\sum_{k=1}^{N_{G}}\sum_{i=0}^{n_{b}-1}\sum_{\sigma=\pm
1}w_{l}(k,i,\sigma)\left[\delta(\sigma+\sigma^{l}_{k}(i))P(\vec{\sigma}^{l}_{k}(i))-\delta(\sigma-\sigma^{l}_{k}(i))P(\vec{\sigma})\right].$
(50)

Injecting Eq. (47) into the latter equation, we find

$\begin{split}d_{t}P_{N_{G}}(\vec{\Sigma})=\sum_{l=1,2}\sum_{k=1}^{N_{G}}\sum_{\vec{\sigma}\in\mathcal{S}(\vec{\Sigma})}\sum_{i=0}^{n_{b}-1}\sum_{\sigma=\pm
1}w_{l}(k,i,\sigma)\delta(\sigma+\sigma^{l}_{k}(i))P(\vec{\sigma}^{l}_{k}(i))\\\
-\sum_{l=1,2}\sum_{k=1}^{N_{G}}\sum_{\vec{\sigma}\in\mathcal{S}(\vec{\Sigma})}\sum_{i=0}^{n_{b}-1}\sum_{\sigma=\pm
1}w_{l}(k,i,\sigma)\delta(\sigma-\sigma^{l}_{k}(i))P(\vec{\sigma}).\end{split}$
(51)

Spins that are in box $k$ but do not contribute to $M_{k}$ because they are
outside the vicinity of the central spin give a zero contribution due to the
sum on $\vec{\sigma}$. We can therefore transform $\sum_{i=1}^{n_{b}-1}$ into
$\sum_{i=1}^{n(a)}$ where we recall that we labeled from $i=0$ to $i=n(a)$ the
spins (and fields) belonging to the neighborhood of the central spin (and
field) $i=0$.

We introduce

$A_{k}^{l}(\vec{\Sigma})=\sum_{i=0}^{n(a)}\sum_{\vec{\sigma}\in\mathcal{S}(\vec{\Sigma})}\delta(\sigma^{l}_{k}(i)-\sigma)w_{l}(k,i,\sigma)P(\vec{\sigma})$
(52)

corresponding to the last term of the previous equation without the sum over
$l$, $k$ and $\sigma$. We seek to write the first term of Eq. (51) as a
function of $A_{k}^{l}$. We denote $\Delta_{k}^{l}$ the operator calculating
the difference when flipping a spin ($l=1$) or a field ($l=2$) $\sigma$ in box
$k$. We have $\Delta_{k}^{l}\vec{\Sigma}=-\sigma\vec{d}_{k}^{l}/n(a)$ with

$\displaystyle\vec{d}^{1}_{k}=$
$\displaystyle(0,\dots,\underbrace{2}_{k\text{th
position}},\dots,0)\in\mathbb{R}^{2N_{G}},\vec{d}^{2}_{k}=$
$\displaystyle(0,\dots,\underbrace{2}_{\text{position
}N_{G}+k},\dots,0)\in\mathbb{R}^{2N_{G}}.$ (53)

We transform the sum on the configurations
$\vec{\sigma}\in\mathcal{S}(\vec{\Sigma})$ with $\sigma_{k}^{l}(i)=-\sigma$
into a sum on configurations
$\vec{\sigma}^{\prime}\in\mathcal{S}(\vec{\Sigma}+\sigma\vec{d}_{k}^{l}/n(a))$
with a spin (or field) $\sigma_{k}^{l}(i)=+\sigma$ of a configuration
$\vec{\sigma}^{\prime}$. Eq. (51) thus becomes

$d_{t}P_{N_{G}}(\vec{\Sigma})=\sum_{l=1,2}\sum_{k=1}^{N_{G}}\sum_{\sigma=\pm
1}\left[A_{k}^{l}\left(\vec{\Sigma}+\frac{\sigma\vec{d}_{l,k}}{n(a)}\right)-A_{k}^{l}(\vec{\Sigma})\right].$
(54)

We now consider a limit of large $n(a)$, which allows us to do a perturbative
expansion in $1/n(a)$ of
$A_{k}^{l}\left(\vec{\Sigma}+\frac{\sigma\vec{d}_{l,k}}{n(a)}\right)$. We keep
only the terms up to order $1/n(a)$, leading to the following equation:

$\begin{split}d_{t}P_{N_{G}}(\vec{\Sigma})=\sum_{k=1}^{N_{G}}\left[\underset{E_{1}}{\underbrace{\sum_{\sigma=\pm
1}\frac{2\sigma}{n(a)}\partial_{M_{k}}A_{k}^{1}(\vec{\Sigma})}}+\underset{F_{1}}{\underbrace{\sum_{\sigma=\pm
1}\frac{2}{n(a)^{2}}\partial^{2}_{M_{k}}A_{k}^{1}(\vec{\Sigma})}}\right.\\\
+\left.\underset{E_{2}}{\underbrace{\sum_{\sigma=\pm
1}\frac{2\sigma}{n(a)}\partial_{H_{k}}A_{k}^{2}(\vec{\Sigma})}}+\underset{F_{2}}{\underbrace{\sum_{\sigma=\pm
1}\frac{2}{n(a)^{2}}\partial^{2}_{H_{k}}A_{k}^{2}(\vec{\Sigma})}}\right].\end{split}$
(55)

In the previous equation, we introduced four different terms
$E_{1},E_{2},F_{1}$ and $F_{2}$, which we now seek to simplify. We assume that
$M_{k}(i)\approx M_{k}$ and $H_{k}(i)\approx H_{k}$ and keep only the first
order in $M_{k}(i)-M_{k}$ and in $H_{k}(i)-H_{k}$. We thus have

$w_{1}(k,i,\sigma)\approx\frac{1}{1+e^{2\sigma(J_{1}M_{k}+H_{k})/T}}-\frac{2\sigma
e^{2\sigma(J_{1}M_{k}+H_{k})/T}[J_{1}(M_{k}(i)-M_{k})+H_{k}(i)-H_{k}]}{T(1+e^{2\sigma(J_{1}M_{k}+H_{k})/T})^{2}}.$
(56)

We first consider the lowest order term in $M_{k}(i)-M_{k}$ and
$H_{k}(i)-H_{k}$. This term corresponds to the transition rates for a mean-
field model inside the box $k$, since the transition rates no longer depend on
the position $i$ in the box. Its contribution to the $E_{1}$ term is therefore

$[M_{k}-\tanh(H_{k}/T)]P_{N_{G}}(\vec{\Sigma}).$ (57)

For the higher-order term in $M_{k}(i)-M_{k}$ and $H_{k}(i)-H_{k}$, it is
assumed that $M_{k}$ and $H_{k}$ are small, which is valid near the transition
temperature only, such that

$\frac{2\sigma
e^{2\sigma(J_{1}M_{k}+H_{k})/T}}{T(1+e^{2\sigma(J_{1}M_{k}+H_{k})/T})^{2}}\sim\frac{\sigma}{2T}.$
(58)

The contribution of the higher-order term in $M_{k}(i)-M_{k}$ and
$H_{k}(i)-H_{k}$ to the $E_{1}$ term thus becomes

$-\frac{1}{Tn(a)}\sum_{\sigma=\pm
1}\sum_{i=1}^{n(a)}\sum_{\vec{\sigma}\in\mathcal{S}(\vec{\Sigma})}\delta(\sigma^{1}_{k}(i)-\sigma)\sigma^{2}\left[J_{1}(M_{k}(i)-M_{k})+H_{k}(i)-H_{k}\right]P(\vec{\sigma})$
(59)

which, summed over $\sigma=\pm 1$, becomes

$-\frac{1}{Tn(a)}\sum_{i=1}^{n(a)}\sum_{\vec{\sigma}\in\mathcal{S}(\vec{\Sigma})}\left[J_{1}(M_{k}(i)-M_{k})+H_{k}(i)-H_{k}\right]P(\vec{\sigma}).$
(60)

We wish to write $M_{k}(i)-M_{k}$ and $H_{k}(i)-H_{k}$ solely as a function of
$M_{k}$ and $H_{k}$. We introduce a position variable $\mathbf{x}$ such that
two boxes are separated by a distance $|\Delta\mathbf{x}|=1$. Two neighboring
spins are then separated by the distance
$|\mathbf{x}(i)-\mathbf{x}(j)|=1/(2\lfloor a\rfloor+1)$. We write
$\mathbf{x}_{k}$ the position of the central spin ($i=0$) of box $k$. We write
$\mathbf{x}_{k}+\delta\mathbf{x}_{k}(i)$ the position of spin $i$ in box $k$.
We have :

$M_{k}(i)-M_{k}=\\!\\!\\!\sum_{m\in\\{1,2,3\\}}\\!\\!\\!\delta\mathbf{x}_{k}(i)\cdot\mathbf{u}_{m}\partial_{m}M(\mathbf{x}_{k})+\frac{1}{2}\\!\\!\sum_{m,n\in\\{1,2,3\\}}\\!\\!\\!\\!\delta\mathbf{x}_{k}(i)\cdot\mathbf{u}_{m}\partial_{m}\partial_{n}M(\mathbf{x}_{k})\delta\mathbf{x}_{k}(i)\cdot\mathbf{u}_{n},$
(61)

where $\mathbf{u}_{m}$ is a unit vector oriented along the $m\in\\{1,2,3\\}$
axis in $3$-dimensional space. For the linear and quadratic terms
corresponding to $m\neq n$, the sum over $i$ is zero by symmetry. So we have :

$\frac{1}{n(a)}\sum_{i=1}^{n(a)}(M_{k}(i)-M_{k})=\sum_{m\in\\{1,2,3\\}}\frac{1}{2n(a)}\sum_{i=1}^{n(a)}(\delta\mathbf{x}_{k}(i)\cdot\mathbf{u}_{m})^{2}\partial_{m}^{2}M(\mathbf{x}_{k}),$
(62)

and we introduce

$C=\frac{1}{2n(a)}\sum_{i=1}^{n(a)}(\delta\mathbf{x}_{k}(i)\cdot\mathbf{u}_{1})^{2}$
(63)

which, by symmetry, leads to the following equation

$\frac{1}{n(a)}\sum_{i=1}^{n(a)}(M_{k}(i)-M_{k})=C\nabla^{2}M(\mathbf{x}_{k}).$
(64)

Finally, we get

$E_{1}=\big{[}M_{k}-\tanh\left[(J_{1}M_{k}+H_{k})/T\right]-C\Delta(J_{1}M_{k}+H_{k})/T\big{]}P_{N_{G}}(\vec{\Sigma})$
(65)

and similarly, we can simplify $E_{2}$ into

$\begin{split}E_{2}=\left[H_{k}-\tanh[(J_{2}H_{k}+(1-\mu)M_{k})/T]-C\Delta(J_{2}H_{k}+(1-\mu)M_{k})/T\right]P_{N_{G}}(\vec{\Sigma}).\end{split}$
(66)

We come back to the quantity $C$ introduced in equation (63), which depends on
the lattice geometry. For a cubic lattice, it depends only on the range of the
interactions $a$, see Fig. 16. We find that $C$ depends little on $a$ and
$C\approx 0.025\pm 0.005$. As explained in the main text, we thus neglect the
dependence of $C$ on $a$ and take $C=0.025$.

Figure 16: Parameter $C$ introduced in Eq. (63) as a function of the
interaction range $a$ for a cubic lattice.

To simplify $F_{1}$ and $F_{2}$, we consider the lowest order in
$(M_{k},H_{k})$, leading to

$F_{1}=F_{2}=\frac{1}{n(a)}$ (67)

which simplifies Eq. (55) into the following equation

$\begin{split}\partial_{t}&P_{N_{G}}(\vec{\Sigma}=(M_{1},...,M_{N_{G}}H_{1},...,H_{N_{G}}))=\\\
\sum_{k=1}^{N_{G}}\Bigg{[}&\partial_{M_{k}}\left(M_{k}-\tanh\left[\frac{J_{1}M_{k}+H_{k}}{T}\right]-\frac{C}{T}\Delta(J_{1}M_{k}+H_{k})\right)P_{N_{G}}(\vec{\Sigma})\\\
+&\partial_{H_{k}}\left(H_{k}-\tanh\left[\frac{J_{2}H_{k}+(1-\mu)M_{k}}{T}\right]-\frac{C}{T}\Delta(J_{2}H_{k}+(1-\mu)M_{k})\right)P_{N_{G}}(\vec{\Sigma})\\\
+&\frac{1}{n(a)}\partial_{M_{k}}^{2}P_{N_{G}}(\vec{\Sigma})+\frac{1}{n(a)}\partial_{H_{k}}^{2}P(\vec{\Sigma})\Bigg{]}.\end{split}$
(68)

This equation corresponds to the Fokker-Planck equation on
$P_{N_{G}}(\vec{\Sigma})$ given in Eq. (13).

###  B.2 Mean-field coupling of noisy oscillators

In this appendix, we detail how the distribution $P_{1}(r,\theta)$ of local
variables in a given box, with mean-field couplings, was obtained in Eq. (29).
We start from the three-dimensional noisy equation on $z_{k}$ [Eq. (24)] :

$\frac{dz_{k}}{dt}=(\gamma-1)z_{k}-|z_{k}|^{2}z_{k}+\gamma C\Delta
z_{k}+\sqrt{D}\left[\eta_{k}^{R}+i\eta_{k}^{I}\right].$ (69)

Replacing the Laplacian by a mean-field coupling [see Eq. (25)], and
introducing $z_{k}=x_{k}+iy_{k}$, leads to

$\begin{split}\frac{dx_{k}}{dt}&=(\gamma-1)x_{k}-(x_{k}^{2}+y_{k}^{2})x_{k}+\gamma
c\left(\overline{x}-x_{k}\right)+\sqrt{D}\eta_{k}^{R}\\\
\frac{dy_{k}}{dt}&=(\gamma-1)y_{k}-(x_{k}^{2}+y_{k}^{2})y_{k}+\gamma
c\left(\overline{y}-y_{k}\right)+\sqrt{D}\eta_{k}^{I}.\end{split}$ (70)

These equations are Langevin equations associated with the following Fokker-
Planck equation:

$\begin{split}&\frac{\partial P_{N_{G}}(\vec{x},\vec{y})}{\partial t}\\\
&=-\sum_{k}\partial_{x_{k}}[(\gamma-1)x_{k}-(x_{k}^{2}+y_{k}^{2})x_{k}+\gamma
c\left(\overline{x}-x_{k}\right)]P_{N_{G}}(\vec{x},\vec{y})\\\
&-\sum_{k}\partial_{y_{k}}[(\gamma-1)y_{k}-(x_{k}^{2}+y_{k}^{2})y_{k}+\gamma
c\left(\overline{y}-y_{k}\right)]P_{N_{G}}(\vec{x},\vec{y})\\\
&+\frac{D}{2}\sum_{k}\left[\partial_{x_{k}}^{2}+\partial_{y_{k}}^{2}\right]P_{N_{G}}(\vec{x},\vec{y})\end{split}$
(71)

with the notations $\vec{x}=(x_{1},...,x_{N_{G}})$ and
$\vec{y}=(y_{1},...,y_{N_{G}})$. We take the limit $N_{G}\to\infty$, and since
we are considering the limit where all boxes interact, we introduce
$P_{1}(x,y)$, the distribution of variables in a single box, such that

$P_{N_{G}}(\vec{x},\vec{y})=\prod_{k=1}^{N_{G}}P_{1}(x_{k},y_{k}).$ (72)

By injecting the ansatz (72) into Eq. (71), we obtain an equation on
$P_{1}(x,y)$:

$\begin{split}\frac{\partial P_{1}(x,y)}{\partial
t}=&-\partial_{x}\left[(\gamma-1)x-(x^{2}+y^{2})x+\gamma c\int
dx^{\prime}dy^{\prime}P_{1}(x^{\prime},y^{\prime})(x^{\prime}-x)\right]P_{1}(x,y)\\\
&-\partial_{y}\left[(\gamma-1)y-(x^{2}+y^{2})y+\gamma c\int
dx^{\prime}dy^{\prime}P_{1}(x^{\prime},y^{\prime})(y^{\prime}-y)\right]P_{1}(x,y)\\\
&+\frac{D}{2}\left[\partial_{x}^{2}+\partial_{y}^{2}\right]P_{1}(x,y).\end{split}$
(73)

The stationary solution of this equation is:

$\begin{split}P_{1}(x,y)=p_{0}\exp\left[-\frac{2}{D}\left(\frac{1}{4}(x^{2}+y^{2})^{2}-\frac{\gamma-1}{2}(x^{2}+y^{2})+\frac{\gamma
c}{2}\left((\overline{x}-x)^{2}+(\overline{y}-y)^{2}\right)\\!\right)\\!\right]\end{split}$
(74)

where, by definition of $\overline{z}$,

$\overline{z}=\overline{x}+i\overline{y}=\int dxdy(x+iy)P_{1}(x,y)$ (75)

which corresponds to a self-consistent equation. In polar coordinates, we
write $\overline{z}=Re^{i\Theta}$, $z=re^{i\theta}$, which leads to Eqs. (29)
and (31) of the main text.

##  Appendix C Comment on the influence of $J_{1}$

Figure 17: Influence of $J_{1}$ and $J_{2}$ on parameters (a) $U$ controlling
oscillation amplitude and (b) $\alpha U$ controlling noise amplitude, for
$\mu=2$.

Throughout this paper, we have considered $J_{1}=0$ in order to reduce the
space of parameters to be studied. On Fig. 17, we plot in color the amplitude
of $U$ and $\alpha U$ in the $(J_{1},J_{2})$-plane for $\mu=2$. For a given
value of $J_{2}$, increasing $J_{1}$ means increasing $\alpha U$ and therefore
lowering the transition temperature. As $U$ increases slightly, the amplitude
of the magnetization slightly decreases, when all boxes are synchronized. We
can see that the coupling constant $J_{1}$ plays a similar role to the
coupling constant $J_{2}$, so that considering $J_{1}=0$ is not a strong
restriction.

## References

  * [1] G. Nicolis.  Dissipative systems.  Rep. Prog. Phys. 49, 873 (1986).
  * [2] J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler.  The Kuramoto model: A simple paradigm for synchronization phenomena.  Rev. Mod. Phys. 77, 137 (2005).
  * [3] S. Gupta, A. Campa, and S. Ruffo.  Kuramoto model of synchronization: equilibrium and nonequilibrium aspects.  J. Stat. Mech.: Theor. Exp. page R08001 (2014).
  * [4] Y. Cao, H. Wang, Q. Ouyang, and Y. Tu.  The free-energy cost of accurate biochemical oscillations.  Nat. Phys. 11, 772 (2015).
  * [5] B. Nguyen, U. Seifert, and A. C. Barato.  Phase transition in thermodynamically consistent biochemical oscillators.  J. Chem. Phys. 149, 045101 (2018).
  * [6] L. Aufinger, J. Brenner, and F. C. Simmel.  Complex dynamics in a synchronized cell-free genetic clock.  Nat. Commun. 13, 2852 (2022).
  * [7] C. Devailly, C. Crauste-Thibierge, A. Petrosyan, and S. Ciliberto.  Phase-transition oscillations induced by a strongly focused laser beam.  Phys. Rev. E 92, 052312 (2015).
  * [8] F. Collet, M. Formentin, and D. Tovazzi.  Rhythmic behavior in a two-population mean-field Ising model.  Phys. Rev. E. 94, 042139 (2016).
  * [9] D. De Martino and A. C. Barato.  Oscillations in feedback driven systems: thermodynamics and noise.  Phys. Rev. E 100, 062123 (2019).
  * [10] P. Dai Pra, M. Formentin, and P. Guglielmo.  Oscillatory behavior in a model of non-Markovian mean-field interacting spins.  J. Stat. Phys. 179, 690 (2020).
  * [11] D Sinelshchikov, A. Poggialini, M.F. Abbate, and D. De Martino.  Emergence of collective self-oscillations in minimal lattice models with feedback.  Phys. Rev. E 108, 044204 (2023).
  * [12] S. Gualdi, J.-P. Bouchaud, G. Cencetti, M. Tarzia, and F. Zamponi.  Endogenous Crisis Waves: Stochastic Model with Synchronized Collective Behavior.  Phys. Rev. Lett. 114, 088701 (2015).
  * [13] S. D. Yi, S.-K. Baek, G. Chevereau, and E. Bertin.  Symmetry restoration by pricing in a duopoly of perishable goods.  J. Stat. Mech.: Theor. Exp. page P11001 (2015).
  * [14] B. Andrae, J. Cremer, T. Reichenbach, and E. Frey.  Entropy Production of Cyclic Population Dynamics.  Phys. Rev. Lett. 104, 218102 (2010).
  * [15] D. Duan, B. Niu, and J. Wei.  Hopf-Hopf bifurcation and chaotic attractors in a delayed diffusive predator-prey model with fear effect.  Chaos, Solitons and Fractals 123, 206 (2019).
  * [16] J. D. Crawford.  Introduction to bifurcation theory.  Rev. Mod. Phys. 63, 991 (1991).
  * [17] C. Fei, Y. Cao, Q. Ouyang, and Y. Tu.  Design principles for enhancing phase sensitivity and suppressing phase fluctuations simultaneously in biochemical oscillatory systems.  Nat. Commun. 9, 1434 (2018).
  * [18] P. Gaspard.  The correlation time of mesoscopic chemical clocks.  J. Chem. Phys. 117, 8905 (2002).
  * [19] A. C. Barato and U. Seifert.  Cost and precision of Brownian clocks.  Phys. Rev. X 6, 041053 (2016).
  * [20] A. C. Barato and U. Seifert.  Coherence of biochemical oscillations is bounded by driving force and network topology.  Phys. Rev. E 95, 062409 (2017).
  * [21] L. Oberreiter, U. Seifert, and A. C. Barato.  Universal minimal cost of coherent biochemical oscillations.  Phys. Rev. E 106, 014106 (2022).
  * [22] B. Remlein, V. Weissmann, and U. Seifert.  Coherence of oscillations in the weak-noise limit.  Phys. Rev. E 105, 064101 (2022).
  * [23] F. Sagués, J. M. Sancho, and J. García-Ojalvo.  Spatiotemporal order out of noise.  Rev. Mod. Phys. 79, 829 (2007).
  * [24] H.-Y. Xu, Y.-P. Luo, J.-W. Wu, and M.-C. Huang.  Langevin approach for intrinsic fluctuations of chemical reactions with Hopf bifurcation.  Physica D 411, 132612 (2020).
  * [25] L. Guislain and E. Bertin.  Non-equilibrium phase transition to temporal oscillations in mean-field spin models.  Phys. Rev. Lett. 130, 207102 (2023).
  * [26] L. Guislain and E. Bertin.  Discontinuous Phase Transition from Ferromagnetic to Oscillating States in a Nonequilibrium Mean-Field Spin Model.  Phys. Rev. E 109, 034131 (2024).
  * [27] L. Crochik and T. Tomé.  Entropy production in the majority-vote model.  Phys. Rev. E 72, 057103 (2005).
  * [28] T. J. Xiao, Z. Hou, and H. Xin.  Entropy production and fluctuation theorem along a stochastic limit cycle.  J. Chem. Phys. 129, 114506 (2008).
  * [29] T. Xiao, Z. Hou, and H. Xin.  Stochastic Thermodynamics in Mesoscopic Chemical Oscillation Systems.  J. Phys. Chem. B 113, 9316 (2009).
  * [30] A. C. Barato and H. Hinrichsen.  Entropy production of a bound nonequilibrium interface.  J. Phys. A: Math. Theor. 45, 115005 (2012).
  * [31] T. Tomé and M. J. de Oliveira.  Entropy production in nonequilibrium systems at stationary states.  Phys. Rev. Lett. 108, 020601 (2012).
  * [32] C. E. F. Noa, P. E. Harunari, M. J. de Oliveira, and C. E. Fiore.  Entropy production as tool for characterizing nonequilibrium phase transitions.  Phys. Rev. E 100, 012104 (2019).
  * [33] T. Martynec, S. H. L. Klapp, and S. A. M. Loos.  Entropy production at criticality in a nonequilibrium Potts model.  New J. Phys. 22, 093069 (2020).
  * [34] D. S. Seara, B. B. Machta, and M. P. Murrell.  Irreversibility in dynamical phases and transitions.  Nat. Commun. 12, 392 (2021).
  * [35] V. Petrolli, M. Le Goff, M. Tadrous, K. Martens, C. Allier, O. Mandula, L. Hervé, S. Henkes, R. Sknepnek, T. Boudou, G. Cappello, and M. Balland.  Confinement-Induced Transition between Wavelike Collective Cell Migration Modes.  Phys. Rev. Lett. 122, 168101 (2019).
  * [36] G. Peyret, R. Mueller, J. d’Alessandro, S. Begnaud, P. Marcq, R.-M. Mège, J. M. Yeomans, A. Doostmohammadi, and B. Ladoux.  Sustained Oscillations of Epithelial Cell Sheets.  Biophysical Journal 117, 464 (2019).
  * [37] K. Kamino, Y. Kondo, A. Nakajima, M. Honda-Kitahara, K. Kaneko, and S. Sawai.  Fold-change detection and scale invariance of cell–cell signaling in social amoeba.  Proc. Natl. Acad. Sci. USA 114, E4149 (2017).
  * [38] S.-W. Wang and L.-H. Tang.  Emergence of collective oscillations in adaptive cells.  Nat. Commun. 10, 5613 (2019).
  * [39] S. Saha, J. Agudo-Canalejo, and R. Golestanian.  Scalar active mixture: The nonreciprocal Cahn-Hilliard model.  Phys. Rev. X 10, 041009 (2020).
  * [40] Z. You, A. Baskaran, and M. C. Marchetti.  Nonreciprocity as a generic route to traveling states.  Proc. Natl. Acad. Sci. USA 117, 19767 (2020).
  * [41] D. Martin, D. Seara, Y. Avni, M. Fruchart, and V. Vitelli.  The transition to collective motion in nonreciprocal active matter: coarse graining agent-based models into fluctuating hydrodynamics.  (2024).  Preprint arXiv:2307.08251.
  * [42] T. Risler, J. Prost, and F. Jülicher.  Universal Critical Behavior of Noisy Coupled Oscillators.  Phys. Rev. Lett. 93, 175702 (2004).
  * [43] T. Risler, J. Prost, and F. Jülicher.  Universal Critical Behavior of Noisy Coupled Oscillators: A Renormalization Group Study.  Phys. Rev. E 72, 016130 (2005).
  * [44] Y. Avni, M. Fruchart, D. Martin, D. Seara, and V. Vitelli.  The Non-Reciprocal Ising Model.  (2023).  Preprint arxiv:2311.05471.
  * [45] V. Privman.  Finite Size Scaling and Numerical Simulation of Statistical Systems.  WORLD SCIENTIFIC 1990.
  * [46] M. E. Fisher and M. N. Barber.  Scaling Theory for Finite-Size Effects in the Critical Region.  Phys. Rev. Lett. 28, 1516 (1972).
  * [47] A. P. Gottlob and M. Hasenbusch.  Critical Behaviour of the 3D XY-model: A Monte Carlo Study.  Physica A 201, 593 (1993).
  * [48] K. Huang.  Statistical Mechanis, 2nd Ed. Wiley India Pvt. Limited 2008.
  * [49] K. Wood, C. Van den Broeck, R. Kawai, and K. Lindenberg.  Universality of Synchrony: Critical Behavior in a Discrete Model of Stochastic Phase-Coupled Oscillators.  Phys. Rev. Lett. 96, 145701 (2006).
  * [50] K. Wood, C. Van Den Broeck, R. Kawai, and K. Lindenberg.  Critical Behavior and Synchronization of Discrete Stochastic Phase-Coupled Oscillators.  Phys. Rev. E 74, 031113 (2006).
  * [51] I. S. Aranson and L. Kramer.  The world of the complex Ginzburg-Landau equation.  Rev. Mod. Phys. 74, 99 (2002).