# Discontinuous phase transition from ferromagnetic to oscillating states in a
nonequilibrium mean-field spin model

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

###### Abstract

We study a nonequilibrium ferromagnetic mean-field spin model exhibiting a
phase with spontaneous temporal oscillations of the magnetization, on top of
the usual paramagnetic and ferromagnetic phases. This behavior is obtained by
introducing dynamic field variables coupled to the spins through non-
reciprocal couplings. We determine a nonequilibrium generalization of the
Landau free energy in terms of the large deviation function of the
magnetization and of an appropriately defined smoothed stochastic time
derivative of the magnetization. While the transition between paramagnetic and
oscillating phase is continuous, the transition between ferromagnetic and
oscillating phases is found to be discontinuous, with a coexistence of both
phases, one being stable and the other one metastable. Depending on parameter
values, the ferromagnetic points may either be inside or outside the limit
cycle, leading to different transition scenarios. The stability of these
steady states is determined from the large deviation function. We also show
that in the coexistence region, the entropy production has a pronounced
maximum as a function of system size.

##  I Introduction

A number of systems driven far from equilibrium are known to exhibit
spontaneous collective oscillations. This is the case for instance for coupled
oscillators, like the Kuramoto model with distributed frequencies [1, 2], or
in models of identical coupled noisy oscillators [3, 4]. Interestingly,
spontaneous oscillations have also been reported in systems composed of a
large number of coupled units which individually do not oscillate in the
absence of interaction. Standard examples include different types of chemical
oscillators [5], and recent experimental and theoretical studies have also
reported spontaneous oscillations in populations of biological cells [6, 7],
assemblies of active particles with non-reciprocal interactions [8, 9],
biochemical clocks [10, 11, 12], droplets in a fluid binary mixture [13],
models of population dynamics [14, 15], socio-economic models [16, 17] or
nonequilibrium spin systems [18, 19, 20].

At the deterministic level of description, valid in the infinite system size
limit, the spontaneous emergence of oscillations can be described as a Hopf
bifurcation [21] in the framework of dynamical system theory. However, many
situations of experimental relevance involve mesoscopic systems for which
fluctuations cannot be neglected, as in the case of biochemical clocks for
instance [22]. An important consequence of the presence of fluctuations is
that the coherence time of oscillations becomes finite [23, 24, 25, 26, 27].
At a phenomenological level, the onset of oscillations in a fluctuating system
may be described as a stochastic Hopf bifurcation [28, 29]. Yet, a deeper
understanding would require to cast this phenomenon in the general framework
of nonequilibrium phase transitions, by explicitly connecting the collective
level of description to the microscopic dynamics. One may in particular
interpret the onset of oscillations in a large system of interacting degrees
of freedom by extending to far-from-equilibrium systems the thermodynamic
framework of phase transition introduced at equilibrium. Along this line, the
entropy production density has been shown to play the role of a generalized
thermodynamic potential, with a discontinuous derivative at the onset of
spontaneous oscillations [30, 31, 32, 14, 33, 34, 11, 35, 36, 37]. Another
approach to phase transitions consists in introducing order parameters
associated with spontaneously broken symmetries [38]. At a mean-field level of
description, one may then introduce a Landau free energy and determine its
expansion close to the phase transition. While this approach has been
originally designed for equilibrium systems, several recent works have
extended it to different types of nonequilibrium situations in the context of
spin models, to describe relaxation effects [39, 40], or the driving by an
oscillatory field or by multiple heat baths [41]. Based on a large deviation
theory approach, the spontaneous transition from a paramagnetic to an
oscillating state has been recently described in a nonequilibrium Landau
framework [42].

In this paper, we extend the results of Ref. [42] by considering within the
same nonequilibrium Landau framework the transition from a ferromagnetic state
to a state of spontaneous collective oscillations. We study a mean-field spin
model where spins are coupled to dynamic fields in a non-reciprocal way,
resulting in a breaking of detailed balance which allows for the onset of
oscillations in some parameter range. The presence of spontaneous oscillations
may be interpreted as an instance of a non-reciprocal phase transition [43,
44]. Both spin and field variables are also subjected to ferromagnetic
couplings, with different values. The phase transition is characterized by
determining a large deviation function of the magnetization and of a
stochastic variable playing the role of a smoothed time derivative of the
magnetization. The transition from the ferromagnetic state to the oscillating
state is found to be discontinuous, with a coexistence of the two phases in
the transition region. The large deviation function allows us to determine
which phase is stable or metastable. We also characterize finite size effects
in terms of entropy production.

The paper is organized as follows. The model is defined in Sec. II, the method
is presented in Sec. III and the main results of Ref. [42] on the continuous
phase transition from paramagnetic to oscillating states are summarized and
extended in Sec. IV. Then, Sec. V characterizes a first scenario of
discontinuous phase transition from ferromagnetic to oscillating states,
whereby the limit cycle surrounds the ferromagnetic points. A second scenario,
in which ferromagnetic points stand outside the limit cycle, is studied in
Sec. VI. Conclusions and perspectives are summarized in Sec. VII.

##  II Model description

###  II.1 Definition of the model

We consider a generalization of the kinetic mean-field Ising model with
ferromagnetic interactions introduced in [42], and sharing some similarities
with related models having two spin populations [45, 18] or subjected to a
feedback control [19, 46]. The model involves $2N$ microscopic variables: $N$
spins $s_{i}=\pm 1$ and $N$ fields $h_{i}=\pm 1$. We define the magnetization
$m$ and average field $h$ as

$m=\frac{1}{N}\sum_{i=1}^{N}s_{i}\,,\qquad
h=\frac{1}{N}\sum_{i=1}^{N}h_{i}\,.$ (1)

The stochastic dynamics consists in randomly flipping a single spin $s_{i}=\pm
1$ with rate $w_{s}^{\pm}$, or a single field $h_{i}=\pm 1$ with rate
$w_{h}^{\pm}$. In a mean-field spirit, the flipping rates $w_{s}^{\pm}$ and
$w_{h}^{\pm}$ are independent of $i$, and depend only on $m$ and $h$,

$w_{s}^{\pm}(m,h)=\frac{1}{1+e^{\beta\Delta E_{s}^{\pm}(m,h)}},\quad
w_{h}^{\pm}=\frac{1}{1+e^{\beta\Delta E_{h}^{\pm}(m,h)}},$ (2)

with $\beta=T^{-1}$ the inverse temperature and $\Delta E_{s,h}^{\pm}(m,h)$
the variation of $E_{s,h}(m,h)$ when flipping a spin $s_{i}=\pm 1$ or a field
$h_{i}=\pm 1$, where

$\displaystyle E_{s}(m,h)$
$\displaystyle=-N\left(\frac{J_{1}}{2}m^{2}+\frac{J_{2}}{2}h^{2}+mh\right),$
(3) $\displaystyle E_{h}(m,h)$ $\displaystyle=E_{s}(m,h)+\mu Nhm\,.$ (4)

When $\mu=0$, $E_{h}=E_{s}$ and the transition rates satisfy detailed balance
with respect to the equilibrium probability distribution $P_{\rm eq}\propto
e^{-\beta E_{s}}$. Detailed balance is broken as soon as $\mu\neq 0$, and
$\mu$ can thus be interpreted as a parameter controling the distance to
equilibrium. For fixed values of the interactions $J_{1}$ and $J_{2}$, the
temperature $T$ and the distance to equilibrium $\mu$ are the two control
parameters of the model.

We denote as $\mathcal{C}=(s_{1},\dots,s_{N},h_{1},\dots,h_{N})$ the
microscopic configuration of the system. Flips of the variables $s_{i}$ and
$h_{i}$ are encoded into formal transition rates
$W(\mathcal{C}^{\prime}|\mathcal{C})$ from a configuration $\mathcal{C}$ to a
configuration $\mathcal{C}^{\prime}$. The probability $P(\mathcal{C},t)$ of a
configuration $\mathcal{C}$ at time $t$ evolves according to the master
equation

$\partial_{t}P(\mathcal{C},t)=\sum_{\mathcal{C}^{\prime}(\neq\mathcal{C})}\big{[}W(\mathcal{C}|\mathcal{C}^{\prime})P(\mathcal{C}^{\prime},t)-W(\mathcal{C}^{\prime}|\mathcal{C})P(\mathcal{C},t)\big{]}.$
(5)

###  II.2 Phase diagram in the deterministic limit

We first study the bifurcation diagram of the system obtained in the
deterministic limit $N\to\infty$. We compute the time derivatives
$d_{t}\langle m\rangle$ and $d_{t}\langle h\rangle$ using the master equation
Eq. (5), where $\langle
x\rangle=\sum_{\mathcal{C}}m(\mathcal{C})P(\mathcal{C})$ for any observable
$x$. We assume that the law of large numbers holds in the limit $N\to\infty$
so that $\langle f(m,h)\rangle\to f(m,h)$ for any regular function $f$.
Finally we obtain a set of deterministic equations on $m(t)$ and $h(t)$ (see
Appendix A):

$\displaystyle\frac{dm}{dt}$ $\displaystyle=-m+\tanh[\beta(J_{1}m+h)],$ (6)
$\displaystyle\frac{dh}{dt}$ $\displaystyle=-h+\tanh[\beta(J_{2}h+(1-\mu)m)].$
(7)

We explore regimes where the magnetization $m(t)$ may exhibit oscillations. In
dynamical systems theory, a limit cycle may generally be described in the
plane defined by a variable and its time derivative, thus we introduce a new
variable $\dot{m}=dm/dt$. The set of deterministic equations become

$\displaystyle\frac{dm}{dt}=\dot{m}\,,\qquad\frac{d\dot{m}}{dt}=Y(m,\dot{m})$
(8)

where $Y(m,\dot{m})$ has a lengthy expression, given in Appendix B [see Eq.
(104)]. $Y(m,\dot{m})$ satisfies the symmetry $Y(-m,-\dot{m})=-Y(m,\dot{m})$.
To study the fixed points of Eq. (8) and their stability, we decompose
$Y(m,\dot{m})$ into a $\dot{m}$-independent contribution

$Y(m,0)=-V^{\prime}(m)$ (9)

[see Appendix B, Eq. (105) for its explicit expression] and a
$\dot{m}$-dependent contribution

$\dot{m}g(m,\dot{m})=Y(m,\dot{m})-Y(m,0)\,,$ (10)

which defines the function $g(m,\dot{m})$. From Eq. (8), the fixed points
$(m,\dot{m})=(m_{0},0)$ satisfy $Y(m_{0},0)=0$, and thus $V^{\prime}(m_{0})=0$
according to Eq. (9). One finds in particular that $(m,\dot{m})=(0,0)$ is
always a fixed point of the system, because $V^{\prime}(0)=0$ by symmetry.

Linearizing the dynamics around a fixed point $(m_{0},0)$, $m=m_{0}+\delta m$,
$\dot{m}=\delta\dot{m}$, one has from Eq. (8)

$\frac{d}{dt}\begin{pmatrix}\delta m\\\
\delta\dot{m}\end{pmatrix}=\mathbf{M}\begin{pmatrix}\delta m\\\
\delta\dot{m}\end{pmatrix}$ (11)

with

$\mathbf{M}=\begin{pmatrix}0&1\\\
-V^{\prime\prime}(m_{0})&g(m_{0},0)\end{pmatrix}.$ (12)

The linear stability of the fixed point $(m_{0},0)$ is determined by the sign
of the eigenvalues of the matrix $\mathbf{M}$,

$\lambda_{\pm}=\frac{1}{2}g(m_{0},0)\left(1\pm\sqrt{1-\frac{4V^{\prime\prime}(m_{0})}{g(m_{0},0)^{2}}}\right).$
(13)

The fixed point $(m_{0},0)$ is stable if both $\lambda_{+}$ and $\lambda_{-}$
are negative (or have a negative real part), implying that
$V^{\prime\prime}(m_{0})>0$ and $g(m_{0},0)<0$. We see in particular from Eq.
(13) that the fixed point $(m_{0},0)$ becomes unstable when $g(m_{0},0)$ is
positive. We define the critical temperature $T_{c}=(J_{1}+J_{2})/2$ and the
dimensionless deviation from $T_{c}$,

$\varepsilon=\frac{T_{c}-T}{T_{c}}\,.$ (14)

Using expression (104) of $Y(m,\dot{m})$, we get for $m_{0}=0$ and small
$\varepsilon$ that $g(0,0)=a_{0}\varepsilon$, with $a_{0}=2T/T_{c}$. We also
have $V^{\prime\prime}(0)=(\mu-\mu_{l}(T))/T^{2}$, where we define
$\mu_{l}(T)$ as

$\mu_{l}(T)=1-(J_{1}-T)(J_{2}-T).$ (15)

Hence, the fixed point $(m,\dot{m})=(0,0)$ is linearly stable for $T>T_{c}$
[$\varepsilon<0$, implying $g(0,0)<0$] provided that $\mu>\mu_{l}(T)$
[implying $V^{\prime\prime}(0)>0$], and unstable otherwise. We define
$\mu_{c}=\mu_{l}(T_{c})$.

Two examples of stability diagrams, obtained from the numerical evaluation of
the fixed points and their stability [given by the sign of the eigenvalues of
Eq. (13)], are shown in Fig. 1 for different values of $J_{1}$ and $J_{2}$.
Trajectories and existence of limit cycles are obtained from the numerical
integration of Eqs. (6) and (7). A stable paramagnetic fixed point [denoted as
P in Fig. 1(a,b)] is found at high enough temperature, while this point
becomes unstable at low temperature. For low values of $(T,\mu)$, two
symmetric ferromagnetic stable fixed points (F) are observed. At low $T$ and
high $\mu$, an oscillating behavior (O) is observed. The transition lines
between the three different behaviors meet at a tricritical point ($T_{c}$,
$\mu_{c}$), see Fig. 1(a,b).

Figure 1: Phase diagram of the deterministic dynamics Eq. (8) obtained
numerically for (a) $J_{1}=0.6$, $J_{2}=0.4$ and (b) $J_{1}=J_{2}=0.5$. Three
distinct behaviors are observed: a stable paramagnetic point (P), two stable
ferromagnetic points (F), and an oscillating phase (O). On the thick dashed
and dotted lines, the limit cycle and the ferromagnetic points coexist. On the
blue dotted line, the ferromagnetic points are inside the limit cycle (Type-I
coexistence) and on the red dashed line the ferromagnetic points are outside
the limit cycle (Type-II coexistence). When the two lines meet, the
bifurcation is heteroclinic. The orange line corresponds to $\mu_{l}(T)$ [Eq.
(15)]: in plain for $T>T_{c}$ where it represents the limit of stability of
the paramagnetic points and dotted for $T<T_{c}$ as indication. (c)-(f)
Examples of trajectories $m(t)$ and phase space ($m(t),\dot{m}(t))$: (c),(d)
Type-I coexistence for $J_{1}=0.6$, $J_{2}=0.4$, $T/T_{c}=0.9$ and
$(\mu-\mu_{c})/\mu_{c}=-3.5\times 10^{-3}$ (blue dot in (a)); (e),(f) Type-II
coexistence for $J_{1}=J_{2}=0.5$, $T/T_{c}=0.9$ and
$(\mu-\mu_{c})/\mu_{c}=8\times 10^{-3}$ (red dot in (b)).

Depending on the value of $\mu$, the bifurcation from the paramagnetic point
to a limit cycle at $T_{c}$ which occurs for $\mu>\mu_{c}$ can either be
continuous (supercritical Hopf bifurcation) or discontinuous (subcritical Hopf
bifurcation) [42]. It is generically continuous when the couplings $J_{1}$ and
$J_{2}$ are positive. The transition from the ferromagnetic points to a limit
cycle is discontinuous [except for the particular values $J_{2}=\pm 2+J_{1}$]
and we observe small regions of the parameter space where the ferromagnetic
points and the limit cycle coexist. In Fig. 1, they are represented by thick
dotted and dashed lines. Depending on the values of the parameters ($T$,
$J_{1}$ and $J_{2}$), the ferromagnetic points can be either inside or outside
the limit cycle, which leads to a topological classification of the transition
into two differents types. In the following, we call discontinuous transition
of Type I the case when the ferromagnetic points are inside the limit cycle,
and discontinuous transition of Type II the case when the ferromagnetic points
are outside the limit cycle. A discontinuous transition of Type I is typically
found close to $T_{c}$ for $J_{1}>J_{2}$ (under additional assumptions that
are specified below), as illustrated for $J_{1}=0.6$ and $J_{2}=0.4$ by the
dotted blue line in Fig. 1(a). An example of trajectory $m(t)$ and phase space
plot $(m(t),\dot{m}(t))$ is represented in Fig. 1(c),(d). A discontinuous
transition of Type II is found for $J_{1}\leq J_{2}$ (under additional
assumptions that are specified below) and $T<T_{c}$, see Fig. 1(b) for
$J_{1}=J_{2}=0.5$. The corresponding trajectory $m(t)$ and its phase space
representation $(m(t),\dot{m}(t))$ is plotted in Fig. 1(e),(f). The farther
from $T_{c}$, the closer the ferromagnetic points and the limit cycle are.

Note that in the case $J_{1}=0.6$ and $J_{2}=0.4$, one observes that for lower
temperatures, the ferromagnetic points sit outside the limit cycle [dashed red
line in Fig. 1(a)], similarly to the behavior displayed in Fig. 1(e), (f). At
the point where the dotted blue line meets the dashed red line, a limit cycle
joining the two ferromagnetic points and with an infinite period appears when
the ferromagnetic points loose stability, corresponding to a heteroclinic
bifurcation.

Note also that for $J_{2}=\pm 2+J_{1}$, the transition is neither of type I or
II, but is a continuous transition from the ferromagnetic points to the limit
cycle corresponding to a heteroclinic bifurcation. We do not study this
particular case in this paper, but a comment on the specificity of this case
is made in Sec. VI.5.

###  II.3 Close to the tricritical point

Figure 2: Sign of $v_{1}(T_{c},\mu_{c})$ from Eq. (18) in the plane ($J_{1}$,
$J_{2}$). In regions with $v_{1}>0$, the discontinuous transition between
ferromagnetic and oscillating phases in the $(T,\mu)$ phase diagram is of Type
I: a limit cycle coexists with ferromagnetic points located inside the cycle
(see Fig. 1(a) and Sec. V). In the case $v_{1}=0$ with $J_{1}=J_{2}$ (red
line), the discontinuous transition between ferromagnetic and oscillating
phases in the $(T,\mu)$ phase diagram is of Type II: a limit cycle coexists
with ferromagnetic points located outside the cycle (see Fig. 1(b) and Sec.
VI). The case $v_{1}=0$ with $J_{1}\neq J_{2}$ is not discussed in details in
this work.

In the following, we focus on the transition close to $T_{c}$, i.e., for small
$\varepsilon$ in order to use a perturbative theory. We consider that $m$ is
small such that only the first orders of the series expansion of $V(m)$ are
necessary. One finds at order $m^{4}$ for $V(m)$ and at order $m^{2}$ for
$g(m,0)$ that

$\displaystyle
V(m)=\frac{\mu-\mu_{l}(T)}{2T^{2}}m^{2}+\frac{v_{1}(T,\mu)}{4}m^{4}+V_{0},$
(16) $\displaystyle g(m,0)=a_{0}\varepsilon-a_{1}m^{2},$ (17)

where $V_{0}$ is at this stage an arbitrary constant, and where
$v_{1}(T,\mu)$, $a_{0}$ and $a_{1}$ are given in Appendix B. In particular one
has $a_{0},a_{1}>0$ and

$v_{1}(T_{c},\mu_{c})=\frac{(J_{1}-J_{2})[4-(J_{1}-J_{2})^{2}]}{12(J_{1}+J_{2})}\,.$
(18)

The sign of $v_{1}(T_{c},\mu_{c})$, which plays a key role in the behavior of
the model, thus depends on the relative values of $J_{1}$ and $J_{2}$ (see
Fig. 2).

When $v_{1}>0$, ferromagnetic points exist for $\mu<\mu_{l}(T)$, and are given
by

$m_{0}^{2}=\frac{\mu_{l}(T)-\mu}{T^{2}v_{1}}\,,$ (19)

i.e., nonzero solutions of the equation $V^{\prime}(m_{0})=0$. According to
Eq. (13), and given that $V^{\prime\prime}(m_{0})>0$, ferromagnetic points are
linearly stable when $g(m_{0},0)<0$, which corresponds to

$\mu<\mu_{F}(T)\equiv\mu_{l}(T)-\frac{\varepsilon a_{0}T^{2}v_{1}}{a_{1}}\,.$
(20)

Numerically, one observes that before ferromagnetic points become linearly
unstable upon increasing $\mu$, they coexist over a tiny parameter range with
a limit cycle that surrounds them. An example is given for $J_{1}=0.6$ and
$J_{1}=0.4$ in Fig. 3, which displays the ferromagnetic points and the
extension of the limit cycle as a function of $\mu$ at fixed $\varepsilon$
[Fig. 3(a)], as well as the corresponding trajectories $m(t)$ [Fig. 3(b)] and
the coexisting trajectories in the phase space $(m,\dot{m})$ [Fig. 3(c)].
Unlike for smaller values of $T$, we observe that close to the tricritical
point, the ferromagnetic points and the limit cycle are well separated. We
also observe in this regime that the two symmetries $m\,\mapsto\,-m$ and
$\dot{m}\,\mapsto\,-\dot{m}$ are separately valid to a good approximation,
while they were previously valid only under the simultaneous transformation
$(m,\dot{m})\,\mapsto\,(-m,-\dot{m})$.

When $v_{1}\leq 0$, higher order terms in the expansion of $V(m)$ are
necessary to obtain the ferromagnetic points and their stability. Numerically,
we observe two scenarios: the first one is a heteroclinic bifurcation, the
ferromagnetic points loose stability and a limit cycle with infinite period
arises. This is the case in particular for $J_{2}=\pm 2+J_{1}$, for which
$v_{1}(T_{c},\mu_{c})=0$. The second scenario is that before disappearing, the
ferromagnetic points coexist with a small elliptic limit cycle. This is the
case for $J_{1}=J_{2}=0.5$, see Fig. 4 where an example of the evolution with
$\mu$ (at fixed $\varepsilon$) of the ferromagnetic fixed points and of the
limit cycle are displayed, together with examples of trajectories $m(t)$. When
$v_{1}(T_{c},\mu_{c})<0$, one finds numerically that the ferromagnetic phase
and the paramagnetic phase coexist for $T\gtrsim T_{c},\mu\gtrsim\mu_{l}(T)$,
so that $(T_{c},\mu_{c})$ is no longer a tricritical point, whereas for
$v_{1}(T_{c},\mu_{c})=0$ which is verified for $J_{1}=J_{2}$ and for
$J_{2}=\pm 2+J_{1}$, one finds that the three transition lines meet at
$(T_{c},\mu_{c})$ (see Fig. 1 for $J_{1}=J_{2}=0.5$ for an example of
bifurcation diagram). In the following, we focus on the case where
$v_{1}(T_{c},\mu_{c})\geq 0$ where the three phases meet at the critical point
$(T_{c},\mu_{c})$, in order to perform a perturbative analysis close to the
tricritical point.

Figure 3: Type-I coexistence close to the tricritical point, corresponding to
ferromagnetic points inside the limit cycle, for
$\varepsilon=(T_{c}-T)/T_{c}=10^{-4}$ and $J_{1}=0.6$, $J_{2}=0.4$. (a) Values
of $m$ for the ferromagnetic points (orange lines) and for the limit cycle
(blue shaded area) along the transition. At $\mu=\mu_{F}(T)$ [Eq. (20)]
ferromagnetic points become linearly unstable. (b) Examples of trajectories
for $(\mu-\mu_{F})/\mu_{F}=-5\times 10^{-7}$. (c) Corresponding phase space
representation in the plane ($m$, $\dot{m}$). Figure 4: Type-II coexistence
close to the tricritical point, corresponding to ferromagnetic points outside
the limit cycle, for $\varepsilon=10^{-3}$ and $J_{1}=J_{2}=0.5$. (a) Values
of $m$ for the ferromagnetic points (orange lines) and for the limit cycle
(blue shaded area) along the transition. At $\mu=\mu_{F}(T)$ [here,
$(\mu_{F}-\mu_{c})/\mu_{c}\approx 1.5\times 10^{-5}$] the ferromagnetic points
become unstable. (b) Examples of trajectories for
$(\mu-\mu_{F})/\mu_{F}=-5\times 10^{-6}$. (c) Corresponding phase space
representation in the plane ($m$, $\dot{m}$).

##  III Generalized Landau theory

Figure 5: Trajectories $m(t)$: the dashed blue and orange lines correspond to
deterministic trajectories and the green line to a trajectory for a finite
system size $N=5000$ obtained from stochastic numerical simulations.
Parameters: $J_{1}=J_{2}=0.5$, $\varepsilon=5\times 10^{-2}$ and
$(\mu-\mu_{c})/\mu_{c}=3.88\times 10^{-3}$. Jumps between noisy oscillatory
states and ferromagnetic states are observed.

The deterministic limit provides knowledge on the different stable fixed
points or limit cycles that are present in the system. However it lacks
information on the behavior of the system at finite size $N$, such as
knowledge on the macroscopic fluctuations around the stable points or cycles.
But most importantly, in case of coexistence of solutions in the limit
$N\to\infty$, the deterministic approach fails to predict which solution is
the most stable one at finite but large size $N$. In addition, for moderate
size $N$, one observes jumps between noisy oscillatory states and
ferromagnetic states, as illustrated in Fig. 5. A statistical description of
such a situation where the ferromagnetic points and the limit cycle are both
linearly stable in the deterministic limit would thus be useful. We briefly
recall in this section the nonequilibrium generalization of the Landau theory
developed in [42], which allows for a description of phase transitions to
oscillating states.

###  III.1 Stochastic time derivative $\dot{m}$

We first introduce a new variable $\dot{m}$ that plays the role of a smoothed
time derivative of the magnetization for finite-size systems. Following [42],
we formally define the stochastic derivative $\dot{m}(\mathcal{C})$ of the
magnetization $m(\mathcal{C})$ as

$\dot{m}(\mathcal{C})=\sum_{\mathcal{C}^{\prime}(\neq\mathcal{C})}\left[m\left(\mathcal{C}^{\prime}\right)-m\left(\mathcal{C}\right)\right]W(\mathcal{C}^{\prime}|\mathcal{C}),$
(21)

such that on average $d\langle m\rangle/dt=\langle\dot{m}\rangle$. Eq. (21)
thus associates with each microscopic configuration $\mathcal{C}$ an
observable $\dot{m}(\mathcal{C})$, which is a smoothed time derivative of $m$
because it is averaged over all possible transitions
$\mathcal{C}\to\mathcal{C}^{\prime}$, for a fixed configuration $\mathcal{C}$.
The advantage of this definition is that fluctuations of $\dot{m}$ are
typically on the same scale as that of $m$, which is a key property for the
large deviation approach described below. Taking instead the time derivative
of $m\big{(}\mathcal{C}(t)\big{)}$ would lead to diverging, white-noise-like
fluctuations which are not appropriate to develop a generalized Landau theory.

Under the mean-field assumption, the formal transition rate
$W(\mathcal{C}^{\prime}|\mathcal{C})$ can be reexpressed in terms of the
flipping rates $w_{s}^{\pm}(m,h)$ to flip a spin $s_{i}=\pm 1$ defined in Eq.
(2). When flipping a spin $s_{i}=\pm 1$, the magnetization change is given by
$m(\mathcal{C}^{\prime})-m(\mathcal{C})=\mp 2/N$. Since there are $N(1\pm
m)/2$ possibilities to choose a spin $s_{i}=\pm 1$, one finds:

$\dot{m}=(1-m)w_{s}^{-}(m,h)-(1+m)w_{s}^{+}(m,h).$ (22)

Using the expression (2) of the flipping rates $w_{s}^{\pm}(m,h)$, Eq. (22)
becomes:

$\dot{m}=-m+\tanh[\beta(J_{1}m+h)].$ (23)

Note that the functional relation $\dot{m}(m,h)$ turns out to be identical to
the functional relation (6) obtained in the deterministic limit $N\to\infty$.
However, Eq. (23) is valid for any finite $N$, and the variables $m$ and $h$
are here stochastic variables.

###  III.2 Large deviation function

At finite size $N$, the dynamics of the system is determined by the master
equation (5). Instead of considering $P(\mathcal{C})$ which involves $2^{2N}$
configurations $\mathcal{C}=\\{s_{1},\dots,s_{N},h_{1},\dots,h_{N}\\}$, we
consider $P_{N}(m,\dot{m})$ the joint stationary probability density of the
global observables $m$ and $\dot{m}$. The variations of $m$ and $\dot{m}$
during a transition scale as $1/N$. We introduce $\mathbf{d}_{k}$ such that
$(\Delta m,\Delta\dot{m})=\pm\mathbf{d}_{k}/N$ with $k=1$ when flipping a spin
$s_{i}=\pm 1$ and $k=2$ when flipping a field $h_{i}=\pm 1$, so that we have:

$\displaystyle\textbf{d}_{1}$ $\displaystyle=\left(-2,2-2\beta J_{1}+2\beta
J_{1}(m+\dot{m})^{2}\right),$ (24) $\displaystyle\textbf{d}_{2}$
$\displaystyle=-\left(0,-2\beta+2\beta(m+\dot{m})^{2}\right).$ (25)

We note $NW_{k}^{\pm}$ the coarse-grained transition rates from a
configuration $(m,h)$ to $(m^{\prime},h^{\prime})$, with $m^{\prime}=m\mp 2/N$
and $h^{\prime}=h$ if $k=1$, and $m^{\prime}=m$ and $h^{\prime}=h\mp 2/N$ if
$k=2$. One has:

$\displaystyle W_{1}^{\pm}$ $\displaystyle=\frac{(1\pm m)/2}{1+\exp[\pm
2\beta(J_{1}m+h)]},$ (26) $\displaystyle W_{2}^{\pm}$
$\displaystyle=\frac{(1\pm h)/2}{1+\exp[\pm 2\beta(J_{2}h+(1-\mu)m)]}.$

The coarse-grained master equation governing the evolution of
$P_{N}(m,\dot{m})$ reads [42]:

$\displaystyle\partial_{t}$ $\displaystyle
P_{N}(m,\dot{m})=N\sum_{k,\sigma}\bigg{[}-W_{k}^{\sigma}(m,\dot{m})P_{N}(m,\dot{m})$
(27)
$\displaystyle+W_{k}^{\sigma}\left((m,\dot{m})-\frac{\sigma\textbf{d}_{k}}{N}\right)P_{N}\left((m,\dot{m})-\frac{\sigma\textbf{d}_{k}}{N}\right)\bigg{]}$

where $k=1,2$ and $\sigma=\pm$. For large $N$, the stationary joint
distribution $P_{N}(m,\dot{m})$ takes a large deviation form [47]

$P_{N}(m,\dot{m})\underset{N\to\infty}{\sim}\exp[-N\phi(m,\dot{m})]\,,$ (28)

a property justified by the theory of Markov jump processes with vanishing
jump size [48]. Beside providing information on the fluctuations at finite
system size, the large deviation function (or rate function) $\phi(m,\dot{m})$
determines the macroscopic phase of the system. Linearly stable solutions of
the deterministic equations correspond to local minima of the large deviation
function. When two or more linearly stable solutions are present in the
deterministic equations, the global minima of $\phi$ gives the macroscopic
phase of the system (i.e., the most stable one).

Injecting the large deviation form (28) into Eq. (1) gives to order
$(\nabla\phi)^{2}$,

$\displaystyle\dot{m}\partial_{m}\phi+Y(m,\dot{m})\partial_{\dot{m}}\phi+D_{11}(\partial_{m}\phi)^{2}$
(29)
$\displaystyle+2D_{12}\partial_{m}\phi\partial_{\dot{m}}\phi+D_{22}(\partial_{\dot{m}}\phi)^{2}=0\,,$

with

$(\dot{m},Y(m,\dot{m}))=\sum_{k,\sigma}-\sigma\mathbf{d}_{k}W_{k}^{\sigma}\,.$
(30)

The function $Y(m,\dot{m})$ is found to be the same function as the one
introduced in the deterministic limit in Eq. (8). We introduce
$\mathbf{D}=\\{D_{ij}(m,\dot{m})\\}$ as

$\mathbf{D}\equiv\sum_{k}\mathbf{d}_{k}\cdot\mathbf{d}_{k}^{T}W_{k}^{\sigma},$
(31)

whose explicit expression is given in Appendix B.

We use the decomposition introduced in Eq. (10),

$Y(m,\dot{m})=-V^{\prime}(m)+\dot{m}g(m,\dot{m})\,$ (32)

and we focus, in this paper, on obtaining the large deviation function in
regions where a fixed point ($m_{0},0$) looses stability in the deterministic
limit, i.e., where $g(m_{0},0)$ changes sign, in order to use a perturbative
framework in terms of the small parameter

$u_{0}\equiv g(m_{0},0).$ (33)

We assume $\nabla\phi=O(u_{0})$ since quadratic terms in $\nabla\phi$ have to
balance the contribution in $u_{0}\dot{m}\partial_{\dot{m}}\phi$. At order
$u_{0}$, Eq. (LABEL:eq:phi:quadrat) reduces to

$\dot{m}\partial_{m}\phi-V^{\prime}(m)\partial_{\dot{m}}\phi=0.$ (34)

The general solution of Eq. (34) reads [42]

$\phi(m,\dot{m})=f\big{(}H(m,\dot{m})\big{)}+f_{0}$ (35)

where the function $H(m,\dot{m})$ takes a form similar to a Hamiltonian,

$H(m,\dot{m})=\frac{\dot{m}^{2}}{2}+V(m)\,.$ (36)

The minimum value of $V(m)$ is set to $V=0$, and $f$ is at this stage an
arbitrary function, satisfying for convenience $f(0)=0$. The constant $f_{0}$
in Eq. (35) ensures that the minimal value of $\phi(m,\dot{m})$ is zero.
Contributions of order $u_{0}^{2}$ to Eq. (LABEL:eq:phi:quadrat) yield a
condition determining the derivative $f^{\prime}(H)$ (see [42] for a detailed
derivation)

$f^{\prime}(H)=-\frac{\int_{m_{1}}^{m_{2}}dm\,\sqrt{2[H-V(m)]}\,g\left(m,\sqrt{2[H-V(m)]}\right)}{\int_{m_{1}}^{m_{2}}dm\,\bigg{[}D_{11}\frac{V^{\prime}(m)^{2}}{\sqrt{2[H-V(m)]}}+2D_{12}V^{\prime}(m)+D_{22}\sqrt{2[H-V(m)]}\bigg{]}}\,,$
(37)

where $m_{1}$ and $m_{2}$ are such that $V(m_{1})=V(m_{2})=H$ and $V(m)\leq H$
for $m_{1}\leq m\leq m_{2}$.

The form Eq. (35) of the large deviation function $\phi(m,\dot{m})$ can be
interpreted as giving a statistical weight to deterministic trajectories
determined by the Hamiltonian dynamics

$\frac{dm}{dt}=\frac{\partial
H}{\partial\dot{m}},\qquad\frac{d\dot{m}}{dt}=-\frac{\partial H}{\partial m},$
(38)

valid at order $\varepsilon$, where the Hamiltonian $H(m,\dot{m})$ is defined
in Eq. (36). Denoting $m_{0}$ a minimum of $V(m)$ (we assume here for
simplicity that $V(m)$ has a single minimum or two symmetric minima) the case
$H=V(m_{0})$ corresponds to a fixed point ($m_{0},0$) of the deterministic
dynamics, whereas values $H>V(m_{0})$ correspond to closed orbits, and thus to
oscillations. The most probable value of $H$, and thus the macroscopically
observed behavior, is determined by the global minimum $H^{*}$ of $f(H)$. Note
that the method used here follows similar lines as the determination of
nonequilibrium potentials in dissipative dynamical systems [49, 50, 51].

When the three phases meet at the critical point $(T_{c},\mu_{c})$, the
ferromagnetic points, noted $\pm m_{0}$, have a small amplitude
($m_{0}\leq\varepsilon$), so that $g(m_{0},0)\sim\varepsilon$ with
$\varepsilon=(T_{c}-T)/T_{c}$. Therefore, close to the critical point, i.e.,
for small $\varepsilon$, the framework described above can be used to obtain
the large deviation function $\phi(m,\dot{m})$ and thus the probability
density $P_{N}(m,\dot{m})$ for large $N$.

##  IV Continuous transition from a paramagnetic phase to an oscillating phase

In this section, we briefly recall and extend results presented in [42] for
the continuous transition observed when $J_{1}$ and $J_{2}$ are positive, at
$T=T_{c}$, for $\mu>\mu_{c}$, from a high-$T$ paramagnetic phase to a low-$T$
oscillating phase [vertical green line in Fig. 1(a), (b)], corresponding to a
Hopf bifurcation at the deterministic limit. Above $T_{c}$ ($\varepsilon<0$),
the system is in a paramagnetic phase, whereas below $T_{c}$ ($\varepsilon>0$)
it is in an oscillating phase.

###  IV.1 Transition from a paramagnetic phase to an elliptic limit cycle

####  IV.1.1 Large deviation function

Figure 6: Large deviation function around the paramagnetic-oscillating
transition. (a) Examples of $f(H)$ for $\varepsilon=10^{-3}$ (blue curve) and
for $\varepsilon=-10^{-3}$ (orange curve). The inset represents the shape of
$V(m)$. (b) Colormap of $\phi(m,\dot{m})$ in the space $(m,\dot{m})$ for
$\varepsilon=10^{-3}$. Scalings with $\varepsilon$: $m\sim\varepsilon^{1/2}$,
$\dot{m}\sim\varepsilon^{1/2}$, $H\sim\varepsilon$ and
$f=\phi\sim\varepsilon^{2}$. Parameters: $J_{1}=0.6$, $J_{2}=0.4$,
$(\mu-\mu_{c})/\mu_{c}=1$.

As in the paramagnetic phase $m$ and $\dot{m}$ are small, we use a power-
series expansion of $V(m)$ and $g(m,\dot{m})$ in $m$ and $\dot{m}$. At the
lowest order required to describe the transition, one has

$\displaystyle V(m)=\frac{\mu-\mu_{l}(T)}{2T^{2}}m^{2},$ (39) $\displaystyle
g(m,\dot{m})=a_{0}\varepsilon-a_{1}m^{2}-a_{2}m\dot{m}-a_{3}\dot{m}^{2},$ (40)

where $\varepsilon=(T_{c}-T)/T_{c}$; $\mu_{l}(T)$ and $a_{0}$ were introduced
previously in Eq. (17) and their expressions are recalled in Appendix B along
with the expressions of $a_{1}$ and $a_{2}$, which are all positive
quantities. Compared to Eq. (16), we only keep the quadratic term in $V(m)$,
which is positive around $T_{c}$ for $\mu>\mu_{c}=\mu_{l}(T_{c})$. Higher
order terms are necessary only when the quadratic term is negative, in order
to describe ferromagnetic order, or when the quadratic term is small, which is
discussed in the next section. An illustration of the quadratic potential
$V(m)$ is plotted in the inset of Fig. 6(a).

In the deterministic limit, the oscillating phase appears for $\varepsilon>0$
when the paramagnetic point $(0,0)$ looses stability. The small perturbative
parameter $u_{0}$ introduced in Sec. III [see Eq. (33)] is proportional to
$\varepsilon$, since here $m_{0}=0$ and $u_{0}=g(0,0)=a_{0}\varepsilon$ from
Eq. (40). Hence the formalism introduced in the previous section to obtain the
large deviation function $\phi(m,\dot{m})=f(H)$ can be used to describe the
phase transition for small $\varepsilon$, when the perturbative approach is
valid.

We find from Eq. (37), after integration,

$f(H)=-\varepsilon aH+bH^{2},$ (41)

with

$a=\frac{T^{2}a_{0}}{T^{2}D_{22}+D_{11}[\mu-\mu_{l}(T)]}$ (42)

and

$b=\frac{a_{1}T^{4}+3T^{2}(\mu-\mu_{l})a_{3}}{4(\mu-\mu_{l})[T^{2}D_{22}+D_{11}(\mu-\mu_{l})]}.$
(43)

When $\varepsilon<0$, $f(H)$ is minimal for $H=0$, which corresponds to the
paramagnetic phase. When $\varepsilon>0$, $f(H)$ has a minimum in
$H^{*}=\varepsilon a/2b$, see Fig. 6(a) for examples of $f(H)$ around the
paramagnetic-oscillating transition. The equation $H(m,\dot{m})=H^{*}$
describes an ellipse in the phase space $(m,\dot{m})$ as depicted in Fig.
6(b). The period $\tau$ of a limit cycle described by
$V(m)+\dot{m}^{2}/2=H^{*}$ is obtained as

$\tau=2\int_{-m^{*}}^{m^{*}}\frac{dm}{\dot{m}}=\sqrt{2}\int_{-m^{*}}^{m^{*}}\frac{dm}{\sqrt{H^{*}-V(m)}}\,,$
(44)

where $m^{*}$ is such that $H^{*}=V(m^{*})$. Using expression (39) of $V(m)$,
we find

$\tau=\frac{2\sqrt{2}\pi T}{\sqrt{\mu-\mu_{l}(T)}}.$ (45)

####  IV.1.2 Order parameters

The paramagnetic to oscillating phase transition is characterized by two order
parameters, $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$, where
$\langle x\rangle=\int dmd\dot{m}P_{N}(m,\dot{m})\,x(m,\dot{m})$, and where
the observable $x$ stands for $m^{2}$ or $\dot{m}^{2}$ [or, more generally,
any even function $x(m,\dot{m})$]. Using Eqs. (28) and (35),
$P_{N}(m,\dot{m})$ can be approximated by its properly normalized large
deviation form,

$P_{N}(m,\dot{m})\approx\frac{\exp\big{[}-Nf\big{(}H(m,\dot{m})\big{)}\big{]}}{\int
dm^{\prime}d\dot{m}^{\prime}\,\exp\big{[}-Nf\big{(}H(m^{\prime},\dot{m}^{\prime})\big{)}\big{]}}.$
(46)

Then from Eq. (36), $\int d\dot{m}$ can be replaced by $\int
dH/\sqrt{2[H-V(m)]}$, so that $\langle x\rangle$ becomes, making the
integration intervals explicit:

$\langle
x\rangle=\frac{\int_{-1}^{1}dm\int_{V(m)}^{\infty}\frac{dH}{\sqrt{H-V(m)}}\,x\,e^{-Nf(H)}}{\int_{-1}^{1}dm\int_{V(m)}^{\infty}\frac{dH}{\sqrt{H-V(m)}}e^{-Nf(H)}}\,.$
(47)

From Eqs. (39) and (41) one can then compute the values of $\langle
m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ using Eq. (47). In the limit of
large system sizes, in the paramagnetic phase ($\varepsilon<0$),
$\langle\dot{m}^{2}\rangle=1/(|\varepsilon|aN)$ which vanishes in the limit
$N\to\infty$. By contrast, in the oscillating phase ($\varepsilon<0)$,
$\langle\dot{m}^{2}\rangle=H^{*}\sim\varepsilon$ is constant in the limit
$N\to\infty$. At the transition, $\varepsilon=0$,
$\langle\dot{m}^{2}\rangle=1/\sqrt{\pi bN}$. For smaller system sizes
$N\ll|\varepsilon|^{-2}$, $f(H)$ can be approximated as $f(H)=bH^{2}$ so that
one finds

$\langle\dot{m}^{2}\rangle\sim N^{-1/2},$ (48)

which is independent of $\varepsilon$.

###  IV.2 Non-elliptic limit cycle near the tricritical point

####  IV.2.1 Large deviation function

Figure 7: Large deviation function close to the tricritical point. (a) Example
of $f(H)$ from Eq. (52); the inset represents the shape of $V(m)$. (b)
Colormap of $\phi(m,\dot{m})$ in the space $(m,\dot{m})$. Note the different
scalings $m\sim\varepsilon^{1/2}$ and $\dot{m}\sim\varepsilon$. Parameters:
$J_{1}=0.6$, $J_{2}=0.4$, $\mu=\mu_{c}$ and $\varepsilon=10^{-3}$.

The three phases (paramagnetic, oscillating and ferromagnetic) meet at the
tricritical point $T=T_{c}$ and $\mu=\mu_{c}$. We now look at how the limit
cycle changes when approaching the tricritical point at $\mu=\mu_{c}$ for
$T<T_{c}$. We have

$\mu_{c}-\mu_{l}(T)=T_{c}^{2}\varepsilon^{2}\,,$ (49)

such that at $\mu=\mu_{c}$, the quadratic term $\propto(\mu_{c}-\mu_{l})m^{2}$
in the expression of $V(m)$ scales as $\varepsilon^{2}m^{2}$. Keeping only the
quadratic term in $V(m)$, we would find $m\sim\varepsilon^{1/2}$ by following
the same reasoning as in Sec. IV.1.1. Hence, the quadratic term becomes of
order $\varepsilon^{3}$ whereas the next-order term $\propto m^{4}$ scales as
$\varepsilon^{2}$, so that the assumption of neglecting the $m^{4}$-term in
the expansion of $V(m)$ is inconsistent. To obtain the right behavior at
$\mu=\mu_{c}$, it is thus necessary to include the contribution of $m^{4}$ in
$V(m)$. At $\mu=\mu_{c}$, we have

$V(m)=\frac{\varepsilon^{2}T_{c}^{2}}{2T^{2}}m^{2}+\frac{v_{1}(T,\mu)}{4}m^{4},$
(50)

with $v_{1}(T,\mu)$ given in Appendix B. In the following, we assume that the
scaling $m\sim\varepsilon^{1/2}$ remains valid close to $\mu_{c}$, and we
check below that the assumption is consistent. [Note that the reason why the
scaling $m\sim\varepsilon^{1/2}$ will eventually prove valid is different from
the one mentioned in the previous paragraph, which relies only on the
quadratic term in $V(m)$]. Under this assumption, the term in
$\varepsilon^{2}m^{2}\sim\varepsilon^{3}$ in Eq. (50) is negligible compared
to the term of order $m^{4}$. Hence, to leading order, we can use the simple
form

$V(m)=\frac{v_{1}(T,\mu)}{4}m^{4}.$ (51)

We also consider that $v_{1}(T_{c},\mu_{c})>0$ as discussed in Sec. II.2 (see
Fig. 2 for possible values of $J_{1}$, $J_{2}$) and that $v_{1}(T,\mu_{c})$
remains positive for small $\varepsilon$. An example of $V(m)$ is plotted in
the inset of Fig. 7(a). Using the same expression of $g(m,\dot{m})$ as before,
Eq. (40), we obtain from Eq. (37) after integration,

$f(H)=-\frac{\varepsilon a_{0}}{D_{22}}H+cH^{3/2}$ (52)

with

$c=\frac{8\,\Gamma\left(\frac{3}{4}\right)^{4}a_{1}}{5\pi^{2}D_{22}\sqrt{v_{1}}}\,$
(53)

where $\Gamma(x)=\int_{0}^{\infty}t^{x-1}e^{-t}dt$. One finds that $f(H)$ has
a minimum for

$H^{*}=\left(\frac{2\varepsilon a_{0}}{3cD_{22}}\right)^{2}$ (54)

[see Fig. 7(a)] corresponding to a limit cycle in the phase space
$(m,\dot{m})$. However, $H$ is no longer quadratic in $m$, because of the
quartic form (51) of $V(m)$. Hence the limit cycle is no longer elliptic, see
Fig. 7(b) for a colormap of $\phi(m,\dot{m})$ in the phase space
$(m,\dot{m})$. We recall that for $\mu\gg\mu_{c}$, one had
$m\sim\dot{m}\sim\varepsilon^{1/2}$, whereas now one finds distinct scalings
$m\sim\varepsilon^{1/2}$ and $\dot{m}\sim\varepsilon$. These scalings are
obtained by using $H^{*}\sim\varepsilon^{2}$ from Eq. (54), and the expression
(36) of $H(m,\dot{m})$ together with the quartic form (51) of $V(m)$. Finally,
evaluating the oscillation period using Eq. (44), we find

$\tau=\frac{4\sqrt{\pi}\Gamma\\!\left(\frac{5}{4}\right)}{\Gamma\left(\frac{3}{4}\right)(v_{1}H^{*})^{1/4}}\,.$
(55)

As $H^{*}\sim\varepsilon^{2}$, the period diverges as $\varepsilon^{-1/2}$
when $\varepsilon\to 0$.

Figure 8: (a) $P_{N}(m,\dot{m})$ obtained from stochastic numerical
simulations. (b) Theoretical $P_{N}(m,\dot{m})$ evaluated by including
leading-order corrections, as given in Eq. (60). Parameters:
$\varepsilon=10^{-2}$, $\mu=\mu_{l}(T)$, $N=10^{7}$.

###  IV.3 Comparison with stochastic simulations and need for higher order
corrections

The method used to obtain analytically the large deviation function, developed
in Sec. III, relies on two main assumptions: $N$ is large and $\varepsilon$ is
small. We now compare the analytical results with numerical simulations of the
stochastic spin model. We use the Gillespie algorithm [52] to simulate the
stochastic dynamics with the rates given by Eq. (26) for a time-interval
$\tau$. The initial condition for the simulations are $m=0$ and $h=0$. In this
algorithm, time-steps are of $O(N^{-1})$ such that the number of steps
required to have $t=1$ is of $O(N)$. To observe an non-elliptic limit cycle
close to $T_{c}$, one needs to have $N|f(H^{*})|\gg 1$ which corresponds to
$N\varepsilon^{3}\gg 1$ [see Eqs. (52) and Eq. (54)]. For example, a value
$\varepsilon=10^{-3}$ would require simulations with at least $N\sim 10^{11}$.
To obtain data with converged statistics depicting the transition, we make
simulations for larger $\varepsilon$ where the approximations made in Sec III
are no longer expected to be quantitatively valid. We now discuss the notable
differences observed in numerical simulations due to larger $\varepsilon$
values and smaller system sizes $N$.

In Fig. 8, we plot $P_{N}(m,\dot{m})$ obtained from stochastic simulations for
$(T_{c}-T)/T_{c}=\varepsilon=10^{-2}$ and for $\mu=\mu_{l}(T)$. We take
$\mu=\mu_{l}(T)$ (instead of $\mu=\mu_{c}$) so that the term in front of
$m^{2}$ in $V(m)$ is exactly zero. Significant discrepancies are observed
between the simulation results and the theoretical predictions of the
perturbative approach described in Sec. III. From leading order calculations
in $\varepsilon$, we obtain that $\phi(m,\dot{m})=f(H)$ with
$H=V(m)+\dot{m}^{2}/2$ which has, in particular, two consequences. First, $m$
and $\dot{m}$ are decoupled, and the symmetries $m\,\mapsto\,-m$ and
$\dot{m}\,\mapsto\,-\dot{m}$ hold independently. Second, the probability
$P_{N}(m,\dot{m})$ is uniform along the limit cycle, corresponding in a
dynamical view (and in the deterministic limit) to a constant ‘speed’ along
the limit cycle. Both of these characteristics are not observed in the
stochastic simulations, see Fig. 8(a).

These differences come from higher order corrections, in $\varepsilon$ and in
$N$ in the probability density $P_{N}(m,\dot{m})$. Similar corrections were
studied in [49] in the context of noisy dynamical systems. We now give an
example of the first corrections for $\mu=\mu_{l}(T)$. The detailed steps of
the derivation are given in [53]. To perform a systematic
$\varepsilon$-expansion, we introduce rescaled variables $\tilde{H}=H/H^{*}$
and $x=m/m_{0}$ with $H^{*}$ given in Eq. (54) and
$m_{0}=(4H^{*}/v_{1})^{1/4}$ consistently with Eq. (19). At lowest order in
$\varepsilon$, from Eq. (52), one finds
$\phi(m,\dot{m})\sim\varepsilon^{3}\tilde{\phi}_{1}(x,\tilde{H})$, where
$\tilde{\phi}_{1}$ is a rescaled function independent of $\varepsilon$. We
introduce $C(m,\dot{m})$ the correction of $O(N^{0})$ to $\ln
P_{N}(m,\dot{m})$,

$P_{N}(m,\dot{m})\propto\exp[-N\phi(m,\dot{m})+C(m,\dot{m})],$ (56)

and we expand $\phi(m,\dot{m})$ and $C(m,\dot{m})$ in power series of
$\varepsilon^{1/2}$ [53],

$\phi(m,\dot{m})=\varepsilon^{3}\sum_{i=1}^{\infty}\varepsilon^{(i-1)/2}\tilde{\phi}_{i}(x,\tilde{H})$
(57)

and

$C(m,\dot{m})=\sum_{i=0}^{\infty}\varepsilon^{i/2}\tilde{C}_{i}(x,\tilde{H}).$
(58)

Injecting these expressions into the master equation on $P_{N}(m,\dot{m})$,
one finds equations on $\tilde{\phi}_{i}$ and $\tilde{C}_{i}$ at each order
$i$ [53]. For the lowest order, we find
$\phi_{1}(m,\dot{m})=\varepsilon^{3}\tilde{\phi}_{1}(x,\tilde{H})=f(H)$ [Eq.
(52)] and
$C_{0}(m,\dot{m})=\varepsilon^{1/2}\tilde{C}_{0}(x,\tilde{H})=\varepsilon^{1/2}c_{0}$
a constant given by the normalization of $P_{N}(m,\dot{m})$. For
$\tilde{\phi}_{2}(x,\tilde{H})$ and $\tilde{C}_{1}(x,\tilde{H})$ one finds:

$\displaystyle\tilde{\phi}_{2}(x,\tilde{H})$
$\displaystyle=A\left(1-\sqrt{\tilde{H}}\right)\tilde{H}x\left[a_{0}\,{}_{2}F_{1}\left(-\frac{1}{2},\frac{1}{4},\frac{5}{4},x^{4}\right)-\frac{2}{3}\frac{D_{22}}{\alpha}x^{2}\,{}_{2}F_{1}\left(-\frac{1}{2},\frac{3}{4},\frac{7}{4},x^{4}\right)\right]\,,$
(59) $\displaystyle\tilde{C}_{1}(x,\tilde{H})$
$\displaystyle=B\tilde{H}^{1/4}\left[-6a_{0}x\sqrt{1-x^{4}}+30a_{0}\sqrt{v_{1}}x\,{}_{2}F_{1}\left(-\frac{1}{2},\frac{1}{4},\frac{5}{4},x^{4}\right)-8\frac{D_{22}}{\alpha}x^{3}{}_{2}F_{1}\left(\frac{1}{2},\frac{5}{4},\frac{7}{4},x^{4}\right)\right]\,,$

where $\alpha$, $A$ and $B$ depend on $v_{1}$, $a_{0}$ and $a_{1}$ and are
given in Appendix B, and ${}_{2}F_{1}(a,b,c,x)$ denotes the hypergeometric
function. The correction $\tilde{\phi}_{2}$ changes the orientation and the
shape of the limit cycle, whereas the correction $\tilde{C}_{1}$ breaks the
uniformity of the probability $P_{N}(m,\dot{m})$ along the limit cycle. We
plot in Fig. 8 the following expression of $P_{N}(m,\dot{m})$ that includes
leading corrections,

$P_{N}(m,\dot{m})=\exp[-N(\phi_{1}+\phi_{2})+C_{0}+C_{1}]$ (60)

with the definitions
$\phi_{i}(m,\dot{m})=\varepsilon^{3+(i-1)/2}\tilde{\phi_{i}}(x,\tilde{H})$ and
$C_{i}(m,\dot{m})=\varepsilon^{i/2}\tilde{C}_{i}(x,\tilde{H})$. The main
features of the probability density obtained from the simulations are captured
by these leading corrections.

##  V Type-I discontinuous transition between ferromagnetic and oscillating
phases

In this section, we investigate the properties, near the tricritical point
$(T_{c},\mu_{c})$, of the ferromagnetic to oscillating phase transition where
a limit cycle appears around the ferromagnetic points, called coexistence of
Type I. This case corresponds to $v_{1}(T_{c},\mu_{c})>0$, see Fig. 2 as well
as the phase diagram of Fig. 1(a) and the trajectories displayed in Figs. 1(c)
and 1(e).

###  V.1 Large deviation function and phase diagram

####  V.1.1 Validity of the perturbative approach

We start from the generic expansion of the potential $V(m)$ given in Eq. (16),
recalled here for clarity,

$V(m)=\frac{\mu-\mu_{l}(T)}{2T^{2}}m^{2}+\frac{v_{1}(T,\mu)}{4}m^{4}+V_{0},$
(61)

where $V_{0}$ is chosen such that $V(m)\geq 0$ and its minimal value is zero.
An example of $V(m)$ for $\mu<\mu_{l}(T)$ is given in Fig. 9. We recall that
$g(m,\dot{m})$ is given in Eq. (40). As discussed in Sec. II.3, ferromagnetic
points $m_{0}^{2}=(\mu_{l}(T)-\mu)/T^{2}v_{1}$ exist for $\mu<\mu_{l}(T)$, and
are locally stable for $\mu\leq\mu_{F}$ with
$\mu_{F}(T)=\mu_{l}(T)-\varepsilon a_{0}T^{2}v_{1}/a_{1}$. In this section, we
focus on the region where the ferromagnetic points loose stability
($\mu\approx\mu_{F}$), thus:

$\mu_{l}(T)-\mu\sim\varepsilon\,$ (62)

and one has $m_{0}^{2}\sim\varepsilon$ and thus
$u_{0}=g(m_{0},0)\sim\varepsilon$. For small $\varepsilon=(T_{c}-T)/T_{c}$,
the main assumption made in Sec. III, i.e., that $u_{0}$ is small, is
verified, and we can use the method developed in this section to obtain the
large deviation function and study the phase transition from a ferromagnetic
phase to an oscillating phase.

####  V.1.2 Typical $f(H)$ and phase diagram

Figure 9: Large deviation function for Type-I discontinuous transition between
ferromagnetic and oscillating phases. (a) $f(H)$ determined numerically from
Eq. (37) (full line); the dashed line corresponds to local approximations
described in Sec. V.1.3. Inset: shape of $V(m)$ with two minima. (b) Colormap
of $\phi(m,\dot{m})$ in the plane $(m,\dot{m})$. Note the different scalings
$m\sim\varepsilon^{1/2}$ and $\dot{m}\sim\varepsilon$. Parameters:
$J_{1}=0.6$, $J_{2}=0.4$, $(\mu-\mu_{c})/\mu_{c}=-3.18\times 10^{-5}$ and
$\varepsilon=10^{-3}$.

In general, except for particular cases as the one described in the previous
section, one cannot obtain explicit analytical expressions of $f^{\prime}(H)$
from Eq. (37), and one needs to perform a numerical integration of the
integrals in Eq. (37) to determine $f(H)$. An example of $f(H)$, for
$\mu\leq\mu_{l}(T)$, numerically obtained from Eq. (37), is plotted in Fig.
9(a). We observe that $f(H)$ has two local minima: one in $H=0$ corresponding
to the ferromagnetic points $m=m_{0}$ and $\dot{m}=0$ [since $V(m_{0})=0$],
and one for $H=H^{*}>0$ corresponding to a limit cycle in the phase space
$(m,\dot{m})$. We numerically obtain (not shown) that

$H^{*}\sim\varepsilon^{2}.$ (63)

An example of colormap of $\phi(m,\dot{m})=f\big{(}H(m,\dot{m})\big{)}$ in the
phase space $(m,\dot{m})$ is displayed in Fig. 9(b). Here, the most stable
phase is the oscillating phase as $f(0)>f(H^{*})$. Contrary to Sec. IV, no
analytical expression of $H^{*}$ is available in the present case.

The transition from the ferromagnetic phase to the oscillating phase takes
place when $f(0)=f(H^{*})$; we note $\mu_{t}(T)$ the value of $\mu$ at the
transition. The value of $H$ jumps from $H=0$ to the nonzero value $H^{*}$ at
the transition, meaning that the latter is discontinuous. We obtain
numerically that $(\mu_{c}-\mu_{t})/\mu_{c}\sim\varepsilon$ where
$\varepsilon=(T_{c}-T)/T_{c}$ [see Fig. 10 with
$(\mu_{c}-\mu_{F})/\mu_{c}\sim\varepsilon$].

From the numerical determination of $f(H)$, one obtains a phase diagram in the
space $(\varepsilon,\mu)$ with the determination of the different phases:
ferromagnetic phase (F), oscillating phase (O) or the phase where both
coexist, with one being more stable than the other. A close up on the phase
diagram near the tricritical point for $J_{1}=0.6$ and $J_{2}=0.4$ is plotted
in Fig. 10, where we represent $(\mu-\mu_{F})/\mu_{F}$, with $\mu_{F}$ given
in Eq. (20), in order to visualize the different phases. The limits of
existence of the ferromagnetic and oscillating states can also be obtained in
the deterministic limit, but for the determination of $\mu_{t}$ (which
characterizes the most stable phase) it is necessary to consider finite system
sizes using the large deviation approach.

Figure 10: Close up of the phase diagram of Fig. 1(a) near the tricritical
point, for $J_{1}=0.6$ and $J_{2}=0.4$, in the reduced parameters
$\varepsilon=(T_{c}-T)/T_{c}$ and $(\mu-\mu_{F}(T))/\mu_{F}(T)$, where
$\mu_{F}$ is defined in Eq. (20) and
$(\mu_{c}-\mu_{F}(T))/\mu_{c}\sim\varepsilon$. O corresponds to the
oscillating phase, F to the ferromagnetic phase. In the hatched area, both
phases are locally stable. The most stable phase, given by the global minima
of $f(H)$, is either the ferromagnetic phase [(O)+F] or the oscillating phase
[O+(F)], separated by the transition line $\mu_{t}$. $\mu_{l}(T)$ corresponds
to the limit of existence of the ferromagnetic points in the deterministic
limit, $\mu_{F}$ to their linear stability limit, and $\mu_{O}$ to the limit
of existence of the oscillating state at deterministic level.

####  V.1.3 Local analytical expressions of $f(H)$

In most cases, keeping only the first orders of the series expansions of
$V(m)$ and $g(m,\dot{m})$ is not enough to obtain an analytical expression of
$f^{\prime}(H)$. Still, local approximations can be obtained. Near a minimum
$m_{0}$ of $V(m)$, a quadratic expansion of $V$ gives

$f(H)=f_{F}(H)\equiv-\frac{g(m_{0},0)}{D_{22}+2v_{1}m_{0}^{2}D_{11}}H,$ (64)

where $f_{F}(H)$ stands for the local approximate expression of $f(H)$ in the
ferromagnetic state. We recover that the point $(m,\dot{m})=(m_{0},0)$ is
stable when $g(m_{0},0)<0$. This expression of $f(H)$ is valid for small $H$
only. For $H\gg H^{*}$, we recover that $f(H)\approx cH^{3/2}$, with $c$
defined in Eq. (53). For intermediate values of $H$ ($H\sim H^{*}$), we do not
have an analytical expression of $f(H)$. However the regime $H\gg H^{*}$ is
similar to the one obtained for $\mu=\mu_{c}$, which suggests that the form of
$f(H)$ obtained for $\mu\approx\mu_{c}$ in Eq. (52) remains approximately
valid up to a redefinition of coefficient values. One can perform a local fit
of the form $f(H)=f_{O}(H)+f(H^{*})$ with

$f_{O}(H)=\tilde{c}\left(H^{3/2}-\frac{3}{2}\sqrt{H^{*}}H+\frac{1}{2}H^{*3/2}\right),$
(65)

where the parameters $\tilde{c}$ and $H^{*}$ are fitted on the numerically
evaluated $f(H)$ to get a local approximation of $f(H)$ near $H^{*}$ (see Fig.
9 for an example of a fit of $f(H)$ close to its minimum). The functional form
(65) provides a reasonable description of the large $H$ behavior of $f(H)$,
and is more accurate than a simple parabolic fit around the minimum $H=H^{*}$.

###  V.2 Scalings of order parameters with system size at the transition

####  V.2.1 Large-$N$ scaling at $\mu=\mu_{t}$

Using the two local approximations of $f(H)$ given in Eqs. (64) and (65), we
study the behaviors of $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$
in the large-$N$ limit when approaching the critical point $\varepsilon=0$
where the three phases meet. In the ferromagnetic phase ($\mu<\mu_{t}$), using
Eq. (64) and $g(m_{0},0)=a_{0}\varepsilon-a_{1}m_{0}^{2}$, one finds in the
large-$N$ limit,

$\displaystyle\langle m^{2}\rangle$
$\displaystyle=m_{0}^{2}=\frac{\mu_{l}(T)-\mu}{T^{2}v_{1}},$ (66)
$\displaystyle\langle\dot{m}^{2}\rangle$
$\displaystyle=\frac{D_{22}+2v_{1}m_{0}^{2}D_{11}}{(a_{1}m_{0}^{2}-a_{0}\varepsilon)N}.$
(67)

We recover the results of the deterministic limit for $N\to\infty$, $\langle
m^{2}\rangle=m_{0}^{2}\sim\varepsilon$ and $\langle\dot{m}^{2}\rangle=0$. For
large but finite $N$, we obtain that
$\langle\dot{m}^{2}\rangle\sim\varepsilon^{-1}N^{-1}$. In the oscillating
phase ($\mu>\mu_{t}$), $f(H)$ is minimal in $H^{*}$, so that for large enough
$N$ we can replace $e^{-Nf(H)}$ by $\delta(H-H^{*})$ in Eq. (47), yielding

$\displaystyle\langle m^{2}\rangle$
$\displaystyle=\frac{\int_{-m^{*}}^{m^{*}}dm\frac{m^{2}}{\sqrt{H^{*}-V(m)}}}{\int_{-m^{*}}^{m^{*}}dm\frac{1}{\sqrt{H^{*}-V(m)}}},$
(68) $\displaystyle\langle\dot{m}^{2}\rangle$
$\displaystyle=\frac{\int_{-m^{*}}^{m^{*}}dm\sqrt{2(H^{*}-V(m))}}{\int_{-m^{*}}^{m^{*}}dm\frac{1}{\sqrt{2(H^{*}-V(m))}}},$
(69)

where $m^{*}$ is such that $H^{*}=V(m^{*})$. Both $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ reach constant values at large $N$. We obtain
numerically that $\langle m^{2}\rangle\sim\varepsilon$ and
$\langle\dot{m}^{2}\rangle\sim\varepsilon^{2}$, which are the same scalings as
the one observed for the nonelliptic limit cycle for $\mu=\mu_{c}$ (see Sec.
IV.2).

####  V.2.2 Moderate-$N$ scaling at $\mu=\mu_{t}$

Figure 11: Order parameters $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ as a function of system size $N$, for
$\varepsilon\in[10^{-3},10^{-4},7\times 10^{-6},6\times 10^{-7},4.6\times
10^{-8},3.6\times 10^{-9}$] at $\mu=\mu_{t}(\varepsilon)$. Order parameters
are evaluated numerically from Eqs. (37) and (47). (a) $\langle m^{2}\rangle$
vs $N$. (b) $\langle\dot{m}^{2}\rangle$ vs $N$. The red dashed line
corresponds to the scaling predictions for moderate values of $N$ from Eqs.
(70) and (71), with $\langle m^{2}\rangle\sim N^{-1/3}$ and
$\langle\dot{m}^{2}\rangle\sim N^{-2/3}$. (c) $\langle
m^{2}\rangle/\varepsilon$ and (d) $\langle\dot{m}^{2}\rangle\varepsilon^{2}$
vs the rescaled system size $N/\varepsilon^{-3}$, showing data collapse.
Parameters: $J_{1}=0.6$, $J_{2}=0.3$.

Unlike for large values of $N$, for intermediate values of $N$, $\langle
m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ are found not to depend much on
$\varepsilon$ and on $\mu$.

In Figs. 11(a) and 11(b), we plot $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ for different $\varepsilon$ at the transition,
$\mu=\mu_{t}(\varepsilon)$. We observe that $\langle m^{2}\rangle\sim
N^{-1/3}$ and $\langle\dot{m}^{2}\rangle\sim N^{-2/3}$ for moderate $N$
values. These scaling behaviors can be understood as follows. The expression
of the average $\langle x\rangle$ contains an integral over $H$ with the
factor $\exp[-Nf(H)]$, see Eq. (47). For moderate $N$, the integral is
dominated by the $H^{3/2}$-term in $f(H)$. This can be justified by performing
the change of variable $z=NH^{3/2}$ in the integral. One then finds that
higher order powers of $H$ in $f(H)$ are negligible when $N\gg 1$, while
linear contributions in $H$ are also negligible as long as
$\sqrt{H^{*}}N^{1/3}\ll 1$, i.e., $N\ll(H^{*})^{-3/2}$, where $H^{*}$ is small
for small $\varepsilon$ [see Eq. (63)]. The integral is then dominated by the
contribution of the region $z\sim 1$, i.e., which corresponds to $H\gg H^{*}$.
In addition, Eq. (47) also contains an integral over $m$. Due to the factor
$[H-V(m)]^{-1/2}$ in the integrals, values of $m$ which contribute the most
are where $V(m)\approx v_{1}m^{4}/4$. In a similar way as above, this can be
justified by performing the change of variable $z^{\prime}=m/H^{1/4}$ in the
integral with the expression of $V(m)$ given in Eq. (61). One then finds that
higher order powers of $m$ are negligible when $H\ll 1$ and the quadratic
order is negligible when $v_{0}(\mu-\mu_{l})\ll H^{1/4}$ which is verified for
$H\gg H^{*}$ as $\mu-\mu_{l}\sim\varepsilon$ [Eq. (62)] and
$H^{*}\sim\varepsilon^{2}$ [Eq. (63)]. Using these two approximations on
$f(H)$ and $V(m)$, we find:

$\langle
m^{2}\rangle=\frac{4\Gamma\left(\frac{5}{6}\right)\Gamma\left(\frac{3}{4}\right)^{4}}{\pi^{5/2}\sqrt{v_{1}}c^{1/3}}N^{-1/3},$
(70)

and

$\langle\dot{m}^{2}\rangle=\frac{4\Gamma\left(\frac{7}{6}\right)}{3\sqrt{\pi}c^{2/3}}N^{-2/3},$
(71)

where $c$ is given in Eq. (53). We plot these quantities in red in Figs. 11(a)
and 11(b) alongside the numerical values obtained from Eq. (47) using the
numerical evaluation of $f(H)$ from Eq. (37).

We note $N^{*}$ the crossover value of $N$ between the moderate- and large-$N$
regimes. The crossover takes place when the value of $\langle m^{2}\rangle\sim
N^{-1/3}$ in the moderate-$N$ approximation is comparable to the one in the
large-$N$ approximation, $\langle m^{2}\rangle\sim\varepsilon$. One thus finds
that $N^{*}$ behaves as

$N^{*}\sim\varepsilon^{-3}.$ (72)

[Note that, according to the integration argument above,
$N^{*}\sim(H^{*})^{-3/2}$, implying $H^{*}\sim\varepsilon^{2}$, consistently
with Eq. (63)]. A similar argument for $\dot{m}$ yields the same scaling for
$N^{*}$: in the moderate-$N$ regime ($N\ll N^{*}$),
$\langle\dot{m}^{2}\rangle\sim N^{-2/3}$, while in the large-$N$ approximation
($N\gg N^{*}$), one has $\langle\dot{m}^{2}\rangle\sim\varepsilon^{2}$ in the
oscillating phase and $\langle\dot{m}^{2}\rangle\sim\varepsilon^{-1}N^{-1}$ in
the ferromagnetic phase, which both give a crossover at
$N^{*}\sim\varepsilon^{-3}$.

Focusing on the oscillating phase, these scaling behaviors of $\langle
m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ can be encompassed into two
scaling functions

$\langle
m^{2}\rangle=\varepsilon\,\mathcal{F}_{m}(\varepsilon^{3}N),\quad\langle\dot{m}^{2}\rangle=\varepsilon^{2}\,\mathcal{F}_{\dot{m}}(\varepsilon^{3}N),$
(73)

with asymptotic behaviors $\mathcal{F}_{m}(x)\sim x^{-1/3}$ and
$\mathcal{F}_{\dot{m}}(x)\sim x^{-2/3}$ for $x\to 0$, while both functions go
to constant values for $x\to\infty$. Figs. 11(c) and 11(d) show the data
collapse obtained by plotting the rescaled variables $\langle
m^{2}\rangle/\varepsilon$ and $\langle\dot{m}^{2}\rangle/\varepsilon^{2}$
versus the rescaled system size $N/\varepsilon^{-3}$, for different values of
$\varepsilon$ at $\mu=\mu_{t}(\varepsilon)$.

###  V.3 Detailed study of the crossover regime

We reported above two distinct scaling regimes $N\ll N^{*}$ and $N\gg N^{*}$
of the observables $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ as a
function of system size $N$ for $\mu=\mu_{t}(\varepsilon)$, and we identified
the scaling with $\varepsilon$ of the crossover size $N^{*}$. We now
investigate in more details the behavior of these observables in the crossover
regime $N\sim N^{*}$, now focusing on the effect of the variations of $\mu$
close to the transition value $\mu_{t}$, for a fixed $\varepsilon$. We find in
particular that in the crossover regime, $\langle m^{2}\rangle$ has a non-
monotonic behavior as a function of $N$, whose details significantly depend on
$\mu$. The behavior of $\langle\dot{m}^{2}\rangle$, while monotonic as a
function of $N$, is found to strongly depend on $\mu$.

####  V.3.1 Influence of $\mu$ on the crossover regime

Figure 12: (a) $\langle m^{2}\rangle$ and (b) $\langle\dot{m}^{2}\rangle$ as a
function of system size $N$, for different values of $\mu$ close to $\mu_{c}$,
corresponding to $(\mu_{c}-\mu)/\mu_{c}\in[4.0,3.6,3.3,3.1,2.5]\times 10^{-5}$
from darker to lighter colors. The transition between ferromagnetic and
oscillating states takes place at $\mu_{t}$ given by
$(\mu_{c}-\mu_{t})/\mu_{c}\approx 3.2\times 10^{-5}$. Parameters: $J_{1}=0.6$,
$J_{2}=0.4$, and $\varepsilon=10^{-3}$. Order parameters are evaluated
numerically in the same way as in Fig. 11.

As mentioned above, the observables $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ are seen to have a very weak dependence on
$\varepsilon$ and $\mu$ in the moderate-$N$ regime. For large $N$, the value
of $\langle\dot{m}^{2}\rangle$ is found to be significantly different in the
oscillating phase $(\mu>\mu_{t})$ where
$\langle\dot{m}^{2}\rangle\sim\varepsilon^{2}$ and in the ferromagnetic phase
$(\mu<\mu_{t})$ where $\langle\dot{m}^{2}\rangle\sim\varepsilon^{-1}N^{-1}$.
In contrast, the value of $\langle m^{2}\rangle$ is similar in both phases,
with $\langle m^{2}\rangle\sim\varepsilon$. In Fig. 12, $\langle m^{2}\rangle$
and $\langle\dot{m}^{2}\rangle$ are plotted as a function of system size $N$
in the crossover regime $N\sim N^{*}$, for different values of $\mu$ across
the transition, keeping $\varepsilon$ fixed. After an initial decay for $N\ll
N^{*}$, we observe that $\langle m^{2}\rangle$ slightly increases before
reaching a constant value, as expected in the limit $N\to\infty$. For some
values of $\mu$, like for $(\mu-\mu_{c})/\mu_{c}=-3.3\times 10^{-5}$, a second
decay is observed before reaching the asymptotic constant value. The behavior
of $\langle\dot{m}^{2}\rangle$ is significantly different from that of
$\langle m^{2}\rangle$ as the large-$N$ limit yields two different behaviors
in the ferromagnetic or in the oscillating phase. We observe that for
$\mu<\mu_{t}$, $\langle\dot{m}^{2}\rangle$ first reaches a plateau for a
significant range of $N$, before steeply decreasing to eventually reach the
large-$N$ scaling $\langle\dot{m}^{2}\rangle\sim N^{-1}$.

####  V.3.2 Interpretation as a finite-size phase coexistence

The observed non-trivial behaviors can be given a simple interpretation in
terms of finite-size phase coexistence and metastability. For a finite-size
system, a metastable state has a finite probability to be visited, and this
probability decreases exponentially with system size. Based on this idea, we
introduce a simple decomposition of average values into contributions of each
phase, and show that such a decomposition is sufficient to account for most of
the observed behaviors.

Figure 13: (a) $\langle m^{2}\rangle$ and (b) $\langle\dot{m}^{2}\rangle$ vs
system size $N$. In black, data from Fig. 12 for
$(\mu_{c}-\mu)/\mu_{c}=3.3\times 10^{-5}$. The dashed lines correspond to
$\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ computed for $f=f_{H}$
[Eq. (64)] (orange) and for $f=f_{O}$ [Eq. (65)] (blue). The plain red lines
correspond to $\langle x(m,H)\rangle_{\mathrm{approx}}$ defined in Eq. (74).
Here, $f(H^{*})=1.1\times 10^{-12}$.

Starting from the expression of an average observable given in Eq. (47), we
split the semi-axis $H\geq 0$ into two regions, separated by the value $H_{0}$
corresponding to the local maximum of $f(H)$. For small $H$, i.e., near the
ferromagnetic points, $f(H)$ is linear [see Eq. (64)]. For $H$ around $H^{*}$,
we write $f(H)=f_{O}(H)+f(H^{*})$ with $f_{O}$ from Eq. (65) where
$\tilde{c}$, $H^{*}$ and $f(H^{*})$ are parameters fitted on the numerically
evaluated $f(H)$. We consider $N$ large enough so that the integration
interval can be extended to the entire real axis due to the rapidly decaying
factor $\exp[-Nf(H)]$. Then, for any quantity $x(m,H)$, we use the approximate
expression of the average value $\langle x(m,H)\rangle$,

$\langle x(m,H)\rangle_{\mathrm{approx}}=\frac{\langle x\rangle_{F}+\langle
x\rangle_{O}C_{1}\sqrt{N}e^{-Nf(H^{*})}}{1+C_{1}\sqrt{N}e^{-Nf(H^{*})}}$ (74)

with

$C_{1}=\frac{1}{\sqrt{N}}\frac{\int_{-1}^{1}dm\int_{V(m)}^{\infty}dH\frac{1}{\sqrt{H-V(m)}}e^{-Nf_{O}(H)}}{\int_{-1}^{1}dm\int_{V(m)}^{\infty}dH\frac{1}{\sqrt{H-V(m)}}e^{-Nf_{F}(H)}}.$
(75)

In Eq. (74), $\langle x\rangle_{F}$ (resp. $\langle x\rangle_{O}$) corresponds
to the ‘pure-state’ average computed in the ferromagnetic state with
$f(H)=f_{F}(H)$ defined in Eq. (64) [resp. $f(H)=f_{O}(H)$ in the oscillating
state, see Eq. (65)]. We obtain

$C_{1}=-\frac{4g(m_{0},0)m_{0}\sqrt{v_{1}}H^{*1/4}}{\sqrt{3\pi\tilde{c}}(D_{22}+2v_{1}m_{0}^{2}D_{11})}\int_{m*}^{1}\frac{dm}{\sqrt{H^{*}-V(m)}},$
(76)

with $m^{*}$ such that $V(m^{*})=H^{*}$. The oscillating phase has a
contribution weighted with the factor
$C_{1}\sqrt{N}e^{-Nf(H^{*})}/(1+C_{1}\sqrt{N}e^{-Nf(H^{*})})$. When
$f(H^{*})>0$ (i.e., the ferromagnetic phase is the most stable one), the
contribution of the oscillating phase disappears at large $N$ but is non-
negligible for $N\sim f(H^{*})$. In Fig. 13, the ‘pure-state’ averages
$\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ evaluated using either
$f=f_{H}$ for the ferromagnetic state, or $f=f_{O}$ for the oscillating state,
as well as the ‘mixed-state’ approximation Eq. (74) are compared to the values
of $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ obtained from the
numerically evaluated $f(H)$ (same data as on Figs. 11 and 12). $\langle
m^{2}\rangle$ increases due to the influence of the oscillating phase as
$\langle m^{2}\rangle$ is higher for $f=f_{O}$ than for $f=f_{H}$; then for
larger $N$, it decreases to its expected value in the ferromagnetic state.
When $f(H^{*})<0$ (i.e., the oscillating phase is the most stable one), we
observe a monotonous increase between the moderate-$N$ decay the asymptotic
large-$N$ value (see Fig. 12). The influence of the oscillating phase is even
more pronounced for $\langle\dot{m}^{2}\rangle$, as we observe that for
moderate $N$, $\langle\dot{m}^{2}\rangle$ is almost constant and equal to the
value expected in the oscillating state (obtained using $f=f_{O}$), before
eventually steeply decreasing when $Nf(H^{*})\sim 1$.

####  V.3.3 Comparison with stochastic simulations

Figure 14: Comparison with stochastic simulations of the spin model. (a)
$\langle m^{2}\rangle$ and (b) $\langle\dot{m}^{2}\rangle$ obtained from
stochastic numerical simulations, as a function of $N$, for
$(\mu_{c}-\mu)/\mu_{c}\in[4.44,4.04,3.64,3.24]\times 10^{-3}$ (from darker to
lighter colors). Parameters: $J_{1}=0.6$, $J_{2}=0.4$, and
$\varepsilon=10^{-1}$.

In Fig. 14, $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ obtained
from stochastic simulations of the spin model are plotted for different system
sizes $N$, for $\varepsilon=(T_{c}-T)/T_{c}=10^{-1}$. Qualitatively, the
moderate- and large-$N$ regimes are visible on the data. However, the decay of
$\langle m^{2}\rangle$ in the moderate-$N$ regime is significantly slower than
the theoretically predicted behavior $\langle m^{2}\rangle\sim N^{-1/3}$. This
is most likely due to the fact that $\varepsilon$ is not small enough to enter
the asymptotic low-$\varepsilon$ regime, as discussed below. The moderate-$N$
decay of $\langle\dot{m}^{2}\rangle$ seems better described by the theoretical
prediction $\langle\dot{m}^{2}\rangle\sim N^{-2/3}$, although significant
deviations are also visible. In the large-$N$ regime,
$\langle\dot{m}^{2}\rangle$ reaches a constant value $\sim\varepsilon^{2}$ in
the oscillating phase ($\mu>\mu_{t}$), or decreases as
$\varepsilon^{-1}N^{-1}$ in the ferromagnetic phase ($\mu<\mu_{t}$). The
transition between the moderate- and large-$N$ regimes takes place around
$N^{*}\approx 10^{3}\approx\varepsilon^{-3}$, as expected.

To understand the discrepancies found between stochastic simulations data and
theoretical predictions, we note that the main approximation made to obtain
the power laws $\langle m^{2}\rangle\sim N^{-1/3}$ and
$\langle\dot{m}^{2}\rangle\sim N^{-2/3}$ is the assumption
$\sqrt{H^{*}}N^{1/3}\ll 1$ (see Sec. V.2.2), which is not valid here for
moderate $N$ values. Moreover, we showed in Sec. IV.3 that there are
discrepancies between the theory and the simulations for low values of
$\varepsilon$ when $N$ is not large enough. We discuss this issue in more
details in the next subsection.

####  V.3.4 Discussion on the low-$\varepsilon$ and large-$N$ approximations

Figure 15: (a) $P_{N}(m,\dot{m})$ obtained from stochastic simulations of the
spin model for $N=5\times 10^{5}$. (b) Trajectories ($m(t),\dot{m}(t))$ in the
phase space $(m,\dot{m})$ obtained in the deterministic limit. The color
corresponds to $v(m,\dot{m})^{-1}$, where $v(m,\dot{m})$ is the local speed
along the limit cycle, and can be interpreted as the local density
$p(m,\dot{m})\propto v(m,\dot{m})^{-1}$. The two black dots are added for
visual purpose and correspond to the ferromagnetic points. Parameters:
$J_{1}=0.6$, $J_{2}=0.4$, $\varepsilon=10^{-1}$ and
$(\mu_{c}-\mu)/\mu_{c}=3.5\times 10^{-3}$. At the deterministic limit, both
the limit cycle and the ferromagnetic points are linearly stable.

In Fig. 14, we compared $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$
obtained from the theoretical results of Sec. III with stochastic simulations
of the spin model. We observed that the behavior is qualitatively the same,
but we did not obtain quantitative results. We now discuss the
low-$\varepsilon$ approximation and its consequences. We plot in Fig. 15(a) an
example of $P_{N}(m,\dot{m})$ obtained from numerical simulations for
$N=5\times 10^{5}$, $\varepsilon=10^{-1}$ and $(\mu_{c}-\mu)/\mu_{c}=3.5\times
10^{-3}$. In the deterministic limit, both the limit cycle and the
ferromagnetic points are linearly stable for these parameter values. We
observe a significant difference with the results obtained in this section,
similar to what was observed in Sec. IV.3: there is no individual symmetry
$m\,\mapsto\,-m$ or $\dot{m}\,\mapsto\,-\dot{m}$, and the probability density
is not constant along the limit cycle. In Sec. IV, we showed that higher order
corrections [both in $\phi(m,\dot{m})$ and $C(m,\dot{m})$] to the large
deviation form of $P_{N}(m,\dot{m})$ may account for discrepancies between
analytical predictions and numerical results of stochastic simulations. We
recall that corrections in $\varepsilon$ to the large deviation function
$\phi$ lead to changes in the shape of the limit cycle, and are responsible
for the breaking of the individual reversal symmetry in $m$ and $\dot{m}$. In
contrast, corrections in $N^{-1}$ and $\varepsilon$ to the large deviation
function, given by the function $C(m,\dot{m})$ [Eq. (56)], break the
uniformity of the probability density along the limit cycle. In Sec. IV, we
computed the first corrections analytically for $\mu=\mu_{l}(T)$. However,
these corrections are more complicated to compute for any $\mu$, and we thus
propose here a different way to determine corrections to the large deviation
function. In the deterministic limit $N\to\infty$, Eq. (56) can be rewritten
in the form

$P_{N}(m,\dot{m})\underset{N\to\infty}{\to}\exp[C(m,\dot{m})]\,\delta(\phi(m,\dot{m})),$
(77)

and the probability density on the limit cycle can be obtained from the local
speed

$v(m,\dot{m})=[\left(dm/dt\right)^{2}+\left(d\dot{m}/dt\right)^{2}]^{1/2},$
(78)

since $P(m,\dot{m})\propto v(m,\dot{m})^{-1}$. The deterministic limit
provides information on the location of the minima of the function
$\phi(m,\dot{m})$, which corresponds to the limit cycle, as well as the value
of $C(m,\dot{m})$ on the limit cycle. To obtain these corrections, we compute
the trajectory $(m(t),\dot{m}(t))$ in the deterministic limit [using Eq. (8)]
and we plot, in Fig. 15(b) the trajectories in the phase space $(m,\dot{m})$
where the color can be interpreted as the local density $P(m,\dot{m})\propto
v(m,\dot{m})^{-1}$ along the limit cycle. We observe the same shape of limit
cycle as in the stochastic simulations, and we recover a higher probability
density near the axis $\dot{m}=0$ (close to the ferromagnetic points), in
qualitative agreement with numerical results.

###  V.4 Entropy production

Figure 16: Entropy production $\Sigma$ vs. $N$, for different values of $\mu$
such that
$(\mu_{c}-\mu)/\mu_{c}\in[4.0,3.6,3.4,3.33,3.29,3.28,3.275,3.272,3.271]\times
10^{-5}$ from darker to lighter colors. The inset represents $N_{\Sigma}$, the
value of $N$ for which $\Sigma$ is maximal, as a function of $\mu-\mu_{t}$.
The dashed line indicates an exponent $-1$. The entropy production is
evaluated numerically from Eq. (80) where the averages are computed as in Fig.
11. Parameters: $J_{1}=0.6$, $J_{2}=0.4$, and $\varepsilon=10^{-3}$;
$(\mu_{c}-\mu_{t})/\mu_{c}=3.271\times 10^{-5}$.

Beyond the order parameter $\langle\dot{m}^{2}\rangle$, the transition to an
oscillating phase may also be characterized thermodynamically as a transition
from microscopic to macroscopic irreversibility [31, 37], by introducing the
entropy production density $\sigma=\Sigma/N$ in the limit $N\to\infty$, where
the steady-state entropy production $\Sigma$ identifies with the entropy flux
[54, 55],

$\Sigma=\frac{1}{2}\sum_{\mathcal{C},\mathcal{C}^{\prime}}\big{[}W(\mathcal{C}^{\prime}|\mathcal{C})P(\mathcal{C})-W(\mathcal{C}|\mathcal{C}^{\prime})P(\mathcal{C}^{\prime})\big{]}\,\ln\frac{W(\mathcal{C}^{\prime}|\mathcal{C})}{W(\mathcal{C}|\mathcal{C}^{\prime})}\,.$
(79)

We briefly investigate the influence of the bistability of the system on the
entropy production. In the large-$N$ and small-$\varepsilon$ limits, one has
(see Appendix C and [42])

$\frac{\Sigma}{N}=\left[1+(T-J_{1})^{2}\right]\left\langle\dot{m}^{2}\right\rangle+T^{2}\left\langle
V^{\prime}(m)^{2}\right\rangle.$ (80)

In the paramagnetic phase or ferromagnetic phase, one finds $\Sigma=O(N^{0})$,
and in an oscillating phase $\Sigma=O(N)$. Using Eq. (80), we compute the
entropy production numerically for different system sizes; the results are
plotted in Fig. 16. For large $N$, we recover that $\Sigma$ is independent of
$N$ in the ferromagnetic phase, whereas $\Sigma\sim N$ in the oscillating
phase. However, for moderate values of $N$, one has $\Sigma\sim N^{1/3}$, due
to the scaling $\langle\dot{m}^{2}\rangle\sim N^{-2/3}$ obtained in Sec. V.3.
Note that this scaling is different from the scaling of the usual transition
to an oscillating phase with an elliptic limit cycle, where one finds for
moderate $N$, $\Sigma\sim N^{1/2}$, as a consequence of the scaling
$\langle\dot{m}^{2}\rangle\sim N^{-1/2}$ [see Eq. (48)]. In the ferromagnetic
phase, due to the influence of the oscillating phase, the entropy production
increases before having a steep decrease to its constant value. Interestingly,
this ‘overshoot’ effect is still present for $\mu$ slightly below $\mu_{O}$
(see Fig. 10), that is when the limit cycle no longer exists at the
deterministic level. In this situation, the fluctuations described by the
large deviation function keep track of the nearby existence of the limit cycle
in parameter space, and are still able to generate a non-monotonous behavior.
We introduce $N_{\Sigma}$ the value of $N$ where the entropy production is
maximal. The evolution of $N_{\Sigma}$ with $\mu-\mu_{t}$, where $\mu_{t}$ is
the value of $\mu$ at the transition between the oscillating and the
ferromagnetic phases, is plotted in the inset of Fig. 16 over a range of small
values of $\mu-\mu_{t}$. Numerical data can be approximately described by a
power-law decay $N_{\Sigma}\sim(\mu-\mu_{t})^{-1}$, although no theoretical
prediction is available to support this scaling relation. Accordingly, for $N$
larger than $N_{\Sigma}$, the entropy production drops by an amount
$\Delta\Sigma\sim N_{\Sigma}\sim(\mu-\mu_{t})^{-1}$, before reaching its
asymptotic constant value.

##  VI Type-II discontinuous transition between ferromagnetic and oscillating
phases

In this section, we investigate the properties of the transition of Type II
between the ferromagnetic and oscillating phases, near the tricritical point
$(T_{c},\mu_{c})$. In this case, obtained for $J_{1}=J_{2}$, a small, almost
elliptic limit cycle around the center is observed. The Type II scenario is
illustrated in the phase diagram of Fig. 1(b) and on the trajectories of Figs.
1(d) and 1(f). All figures in this section are obtained with
$J_{1}=J_{2}=0.5$.

###  VI.1 Large deviation function and phase diagram

####  VI.1.1 Validity of the perturbative approach

For $J_{1}=J_{2}$, one has $\mu_{c}=1$ and $T_{c}=J_{1}$. The main difference
with the previous case is that $v_{1}(T,\mu)$, the factor in front of $m^{4}$
in $V(m)$ [Eq. (61)] vanishes at $(T_{c},\mu_{c})$. We have to leading order
in an expansion in $\varepsilon$ and $\mu-\mu_{c}$,

$v_{1}(T,\mu)=-\frac{J^{2}}{12}\varepsilon-\frac{1}{4J^{2}}(\mu-\mu_{c})$ (81)

with $\varepsilon=(T_{c}-T)/T_{c}$. Corrections to Eq. (81) include terms
proportional to $(\mu-\mu_{c})^{2}$, $\varepsilon(\mu-\mu_{c})$ and
$\varepsilon^{2}$. Numerically, we observe that the transition between the
oscillating and ferromagnetic phases takes place for
$(\mu-\mu_{c})\ll\varepsilon$ so that the term in $\mu-\mu_{c}$ can be
neglected and we write $v_{1}=-\alpha\varepsilon$ with $\alpha=J^{2}/12$. As
$v_{1}<0$, higher order terms in the expansion of $V(m)$ are necessary to
compensate for the term $-\varepsilon m^{4}$ and thus to describe the
ferromagnetic points. Numerically, we observe ferromagnetic points whose
amplitude goes to zero with $\varepsilon$. The coefficients of the terms
proportional to $m^{6}$ and $m^{8}$ in the expansion of $V(m)$ obtained from
Eq. (105) scale as $\varepsilon$ for small $\varepsilon$, which would give
ferromagnetic points independent of $\varepsilon$ in this limit, if only these
terms were retained. Expanding $V(m)$ further, we find that the coefficient of
the term proportional to $m^{10}$ is independent of $\varepsilon$. Assuming
that the ferromagnetic point $m_{0}$ results from the balance of the terms in
$m^{4}$ and in $m^{10}$, i.e., $\varepsilon m_{0}^{4}\sim m_{0}^{10}$, yields
$m_{0}\sim\varepsilon^{1/6}$. We now check a posteriori that the assumption to
neglect the terms in $m^{6}$ and $m^{8}$ was valid. For
$m\sim\varepsilon^{1/6}$, one has $\varepsilon m^{6}\sim\varepsilon^{2}$ and
$\varepsilon m^{8}\sim\varepsilon^{7/3}$ which are both much smaller than the
term $m^{10}\sim\varepsilon^{5/3}$ for $\varepsilon\to 0$, so that neglecting
the terms in $m^{6}$ and $m^{8}$ was justified for $m\sim m_{0}$. We thus
write the following minimal form for $V(m)$,

$V(m)=\frac{\mu-\mu_{l}(T)}{2T^{2}}m^{2}-\frac{\alpha\varepsilon}{4}m^{4}+\frac{v_{4}}{10}m^{10}+V_{0}$
(82)

where $v_{4}=J^{2}/81$. $V_{0}$ is such that $V(m)\geq 0$ and the minimal
value of $V(m)$ is zero. An example of the shape of $V(m)$ is given in Fig.
17(a). Until now, we have considered only $V(m)$ with one or two minima,
whereas now it can have three of them. As we now show, this has important
consequences which makes this case of interest, and quite different from the
previous ones. We found that the ferromagnetic points are
$m_{0}\sim\varepsilon^{1/6}$ when $\mu-\mu_{c}\ll\varepsilon$, so that
$u_{0}=g(m_{0},0)\sim\varepsilon^{1/3}$ is small. Thus, for small
$\varepsilon$ and close to $\mu_{c}$, one can use the perturbative method
described in Sec. III to obtain the large deviation function.

####  VI.1.2 Typical $f(H)$ and phase diagram

Figure 17: Type-II transition between ferromagnetic and oscillating states.
(a) Potential $V(m)$; the dashed lines correspond to $m=\pm m_{l}$. (b)
Representation of the different areas in the plane $(m,\dot{m})$. (c) $f(H)$
in the different areas defined in Eq. (83). (d) Colormap of $\phi(m,\dot{m})$
in the plane $(m,\dot{m})$. Note the different scalings
$m\sim\varepsilon^{1/2}$ and $\dot{m}\sim\varepsilon$. Parameters:
$J_{1}=J_{2}=0.5$, $(\mu-\mu_{l})/\mu_{l}=4.55\times 10^{-4}$ and
$\varepsilon=10^{-2}$.

The main difference with the previous case is as follows. As illustrated in
Fig. 17(a), the condition $V(m)=H$ may correspond to six values of $m$ instead
of only two or four previously, when considering values of $H$ close to the
local minima of $V(m)$ with $m\neq 0$ (ferromagnetic points). In the
definition of $f^{\prime}(H)$ given in Eq. (37), we integrate over $m_{1}(H)$
and $m_{2}(H)$ such that $V(m_{1})=V(m_{2})=H$ and $V(m)\leq H$ for
$m\in[m_{1},m_{2}]$. Thus, for a given value of $H$, $f$ can have different
values depending on the range of values of $m$ over which the integral is
computed. We note $m_{l}$ the positive value of $m$ where $V(m)$ has a local
maximum. In the phase space $(m,\dot{m})$ there are four different areas,
which are represented in Fig. 17(b). A first area around the center
corresponds to small $m$ and $\dot{m}$, where $|m|<|m_{l}|$ and
$V(m)<V(m_{l})$, which is denoted area 1 and is represented in green. Two
symmetric domains situated around the ferromagnetic points, where $|m|>m_{l}$
and $V(m)<V(m_{l})$, correspond to area 2 and are represented in red. A last
area for higher values of $H$, is denoted area 3 and is represented in blue.
We define three different functions $f$, one for each area:

$\phi(m,\dot{m})\\!=\\!\left\\{\\!\begin{array}[]{ll}f_{1}(H(m,\dot{m}))\\!&\\!\text{if
}|m|\\!<\\!m_{l}\text{ and }V(m)\\!<\\!V(m_{l})\\\
f_{2}(H(m,\dot{m}))\\!&\\!\text{if }|m|\\!>\\!m_{l}\text{ and
}V(m)\\!<\\!V(m_{l})\\\
f_{3}(H(m,\dot{m}))\\!&\\!\text{otherwise.}\end{array}\right.$ (83)

As $f$ is defined up to a constant in every area, we impose that $f_{1}(0)=0$
and we assume $f(H)$ to be continuous at the border between two different
areas. In Fig. 17(c), (d), examples of $f(H)$ and $\phi(m,\dot{m})$ are
plotted. A limit cycle around the center, and two ferromagnetic points are
locally stable. Once again, the most stable phase is given by the global
minima of $f$, here the oscillating phase.

From the numerical determination of $f(H)$ and its minima, one obtains the
phase diagram in the parameter space ($T,\mu$). We plot in Fig. 18 the phase
diagram for $J_{1}=J_{2}$ close to the tricritical point $(T_{c},\mu_{c})$. We
introduce the line $\mu_{F}$ indicating the existence of the ferromagnetic
points and the line $\mu_{O}$ indicating the existence of the oscillating
state. We also introduce $\mu_{t}(T)$ the value of $\mu$ at the transition,
such that $f_{2}(H(m_{0},0))=f_{1}(H^{*})$ (where $m_{0}$ corresponds to the
ferromagnetic point and $H^{*}$ is where $f_{1}(H)$ is minimal). Numerically,
we obtain (see Fig. 18) that

$(\mu_{t}-\mu_{c})/\mu_{c}\sim\varepsilon^{4/3}.$ (84)

Indeed, the transition almost takes place when the ferromagnetic points
disappear, meaning that $\mu_{t}\approx\mu_{F}$. The ferromagnetic points
disappear when the term in $m^{2}$ balances the $m^{4}$ term in $V(m)$ at
$m_{0}\sim\varepsilon^{1/6}$, so that
$(\mu_{F}-\mu_{l})m_{0}^{2}\sim\varepsilon m_{0}^{4}$. As $\mu_{l}\sim\mu_{c}$
[Eq. (49)], this gives $(\mu_{F}-\mu_{c})/\mu_{c}\sim\varepsilon^{4/3}$. We
observe that when the ferromagnetic phase and the oscillating phase coexist,
the ferromagnetic phase is almost always the most stable one.

####  VI.1.3 Approximate local analytical expressions of $f(H)$

To go beyond the numerical evaluation of $f(H)$, we now try to obtain an
approximate analytical expression of $f(H)$, which will be helpful in
particular to determine the scaling regimes of the order parameters $\langle
m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$. As in Sec. V, one cannot obtain
a full analytic expression of $f(H)$, and we thus focus on local
approximations. For an expression of $f(H)$ around the ferromagnetic points, a
quadratic expansion of $V$ around one of its local minima $m_{0}$ gives
$f_{2}(H)=f_{F}(H)+\bar{f}_{2}$ where

$f_{F}(H)=-\frac{g(m_{0},0)}{D_{22}+6v_{4}m_{0}^{8}D_{11}}\left[H-V(m_{0})\right],$
(85)

and $\bar{f}_{2}$ is a constant such that $f$ is continuous in $H=V(m_{l})$.
We note here that $H(m_{0},0)\neq 0$ for the ferromagnetic points unlike in
the previous section as now $V(m_{0})$ can be nonzero if the global minimum of
$V$ is for $m=0$.

For the area around the center, the leading term of $V(m)$ is the $m^{2}$
term, leading to the same $f(H)$ as for the Hopf bifurcation [see Eq. (41)],
$f_{1}(H)=f_{O}(H)+\bar{f}_{1}$ where

$f_{O}(H)=-\varepsilon aH+bH^{2}+\frac{(\varepsilon a)^{2}}{b}$ (86)

with $a$ and $b$ given in Eqs. (42) and (43), and $\bar{f}_{1}$ a constant
chosen such that $f(H)$ is continuous in $H=V(m_{l})$. The minimum of
$f_{O}(H)$ corresponds to an elliptic limit cycle around the center, as
depicted in Fig. 17(d), with

$H^{*}=\frac{a\varepsilon}{2b}\sim\varepsilon[\mu-\mu_{l}(T)]\,,$ (87)

as $b\sim[\mu-\mu_{l}(T)]^{-1}$ from Eq. (43). At the transition
($\mu=\mu_{t}$), we find [Eqs. (49) and (84)],

$\mu_{t}-\mu_{l}(T)\sim\varepsilon^{4/3}$ (88)

and thus $H^{*}\sim\varepsilon^{7/3}$.

Figure 18: Close up of the phase diagram of Fig. 1(b) for $J_{1}=J_{2}=0.5$
for $T<T_{c}$. $\mu_{F}$ corresponds to the limit of existence of the
ferromagnetic points (obtained from $V(m)$), and $\mu_{O}$ of the limit cycle
(obtained in the deterministic limit). The transition between the F and O
phases takes place for $\mu=\mu_{t}$. In the hatched area both the
ferromagnetic points and the limit cycle are local minima of $f(H)$. As
$\mu_{F}$ and $\mu_{t}$ are very close to each other, in most of the
coexistence region, the ferromagnetic phase is the most stable one. Scalings
of transition lines: $(\mu_{t}-\mu_{c})/\mu_{c}\sim\varepsilon^{4/3}$,
$(\mu_{F}-\mu_{c})/\mu_{c}\sim\varepsilon^{4/3}$ and
$(\mu_{O}-\mu_{c})/\mu_{c}\sim\varepsilon^{2}$, with
$\varepsilon=(T_{c}-T)/T_{c}$.

In the phase diagram of Fig. 18, we observe that the transition line $\mu_{t}$
is very close to the line $\mu_{F}$ indicating the limit of existence of the
ferromagnetic points. Hence in most of the coexistence region, the
ferromagnetic phase is the most stable phase. This can be explained with the
following argument. Around the ferromagnetic points,
$f^{\prime}(H)\sim\varepsilon^{1/3}$ whereas near the limit cycle,
$f^{\prime}(H)\sim\varepsilon$. The slope of $f$ near the ferromagnetic points
is much steeper that around the limit cycle (see Fig. 17 for an example of
$f$). When the area around the ferromagnetic points exists, it rapidly becomes
the global minimum of $f$ when varying $\mu$ at fixed $\varepsilon$.

###  VI.2 Multiple scalings of order parameters with $N$ at the transition

In Figs. 19(a) and 19(b) we plot $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ as a function of $N$ at the transition, for
$\mu=\mu_{t}(\varepsilon)$. We observe three different regimes depending on
the value of $N$ and $\varepsilon$, that are described below.

####  VI.2.1 Large-$N$ scaling at $\mu=\mu_{t}$

Using the local approximations, we obtain the behaviors of $\langle
m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ in the large-$N$ limit when
approaching the critical point $\varepsilon=0$ where the three phases meet. In
the ferromagnetic phase, using Eq. (64), we find that $\langle
m^{2}\rangle=m_{0}^{2}\sim\varepsilon^{1/3}$ and

$\langle\dot{m}^{2}\rangle=\frac{D_{22}+6v_{4}m_{0}^{8}D_{11}}{-g(m_{0},0)N},$
(89)

so that $\langle\dot{m}^{2}\rangle\sim\varepsilon^{-1/3}N^{-1}$ as
$|g(m_{0},0)|\sim m_{0}^{2}\sim\varepsilon^{-1/3}$. Due to the $m^{10}$ term
in $V(m)$, we obtain power laws in $\varepsilon$ with critical exponents quite
different from the corresponding values previously obtained. In the
oscillating phase, one finds $\langle m^{2}\rangle\sim\varepsilon$ similarly
to the elliptic limit cycle obtained in Sec. IV, and
$\langle\dot{m}^{2}\rangle\sim\varepsilon[\mu_{t}-\mu_{l}(T)]$. One has
$\mu_{t}-\mu_{l}\sim\varepsilon^{4/3}$ [Eq. (88)] such that one finds
$\langle\dot{m}^{2}\rangle\sim\varepsilon^{7/3}$ in the oscillating phase.

####  VI.2.2 Moderate-$N$ scaling at $\mu=\mu_{t}$

Figure 19: Order parameters $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ as a function of the system size $N$, for
$\varepsilon\in[10^{-4},10^{-6},10^{-8},10^{-10}]$ from lighter to darker
colors at $\mu=\mu_{t}(\varepsilon)$. Order parameters are evaluated
numerically from Eqs. (37) and (47). (a) $\langle m^{2}\rangle$ vs $N$. (b)
$\langle\dot{m}^{2}\rangle$ vs $N$. The red dashed line corresponds to the
theoretical prediction for low $N$ [Eqs. (90) and (91)], with a scaling of
$N^{-1/6}$ for $\langle m^{2}\rangle$ and $N^{-5/6}$ for
$\langle\dot{m}^{2}\rangle$. (c) $\langle m^{2}\rangle$ rescaled by
$\varepsilon^{1/3}$ and (d) $\langle\dot{m}^{2}\rangle$ rescaled by
$\varepsilon^{5/3}$ versus rescaled system size $N/\varepsilon^{-2}$,
highlighting the crossover between the moderate and intermediate-$N$ regimes.
(e) $\langle m^{2}\rangle$ rescaled by $\varepsilon$ and (f)
$\langle\dot{m}\rangle$ rescaled by $\varepsilon^{7/3}$ versus rescaled system
size $N/\varepsilon^{-10/3}$, highlighting the crossover to the large-$N$
regime. Parameters: $J_{1}=J_{2}=0.5$.

For moderate values of $N$, we observe that $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ decrease with $N$ and are independent of
$\varepsilon$. Similarly to Sec. V.2.2, for low values of $N$ we can keep only
the first $\varepsilon$-independent term in $V(m)$, namely here
$V(m)=v_{4}m^{10}/10$. Using Eq. (37), this approximation gives
$f(H)=AH^{5/6}$, leading to

$\displaystyle\langle m^{2}\rangle\sim N^{-1/6},$ (90)
$\displaystyle\langle\dot{m}^{2}\rangle\sim N^{-5/6},$ (91)

where the exact asymptotic relations including prefactors are given in
Appendix D, and are plotted in dashed red lines in Fig. 19(a,b).

####  VI.2.3 Intermediate-$N$ scaling at $\mu=\mu_{t}$

For intermediate values of $N$, we observe that for all $\varepsilon$,
$\langle m^{2}\rangle\sim N^{-1/2}$ and $\langle\dot{m}^{2}\rangle\sim
N^{-1/2}$. Indeed, for values of $N$ such that $f(H(m_{l}))\sim N^{-1}$ where
$m_{l}$ is the positive junction point between the different areas (see Fig.
17), because the ferromagnetic areas are small, the main contribution to the
integrals corresponds to $m$ in area 1 (in the center). In this area, one has
$f(H)\sim f_{O}(H)$ where $f_{O}(H)$ is given in Eq. (86). The leading
correction in $N$ of $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ is
given by

$\langle\dot{m}^{2}\rangle=\frac{\mu-\mu_{l}(T)}{T^{2}}\langle
m^{2}\rangle=H^{*}+\frac{e^{-NbH^{*2}}}{2\sqrt{\pi bN}},$ (92)

with $H^{*}=a\varepsilon/b$. We recall that
$b\sim(\mu-\mu_{l})^{-1}\sim\varepsilon^{-4/3}$ [Eq. (88)] at the transition,
thus we find $\langle m^{2}\rangle\sim\varepsilon^{-2/3}N^{-1/2}$ and
$\langle\dot{m}^{2}\rangle\sim\varepsilon^{2/3}N^{-1/2}$ for
$N\varepsilon^{2}\ll 1$. We note $N_{1}$ the crossover value of $N$ between
those the moderate- and intermediate-$N$ regimes, and $N_{2}$ the crossover
value between the intermediate- and large-$N$ regimes.

For $N\ll N_{1}$, $\langle m^{2}\rangle\sim N^{-1/6}$ while for $N_{1}\ll N\ll
N_{2}$, $\langle m^{2}\rangle\sim\varepsilon^{-2/3}N^{-1/2}$, so that
$N_{1}\sim\varepsilon^{-2}$. The same argument holds for
$\langle\dot{m}^{2}\rangle$: for $N\ll N_{1}$, $\langle\dot{m}^{2}\rangle\sim
N^{-5/6}$ and for $N_{1}\ll N\ll N_{2}$,
$\langle\dot{m}^{2}\rangle\sim\varepsilon^{2/3}N^{-1/2}$, also implying
$N_{1}\sim\varepsilon^{-2}$. These different scaling behaviors for $\langle
m^{2}\rangle$ around the first crossover regime $N\sim N_{1}$ can be
encompassed into a single scaling function

$\langle m^{2}\rangle=\varepsilon^{1/3}\,\mathcal{F}_{m,1}(\varepsilon^{2}N),$
(93)

where $\mathcal{F}_{m,1}(x)$ asymptotically behaves as
$\mathcal{F}_{m,1}(x)\sim x^{-1/6}$ for $x\to 0$ and $\mathcal{F}_{m,1}(x)\sim
x^{-1/2}$ for $x\to\infty$. In a similar way, $\langle\dot{m}^{2}\rangle$ can
be expressed in terms of a scaling function,

$\langle\dot{m}^{2}\rangle=\varepsilon^{5/3}\,\mathcal{F}_{\dot{m},1}(\varepsilon^{2}N),$
(94)

with asymptotic behaviors $\mathcal{F}_{\dot{m},1}(x)\sim x^{-5/6}$ for $x\to
0$ and $\mathcal{F}_{\dot{m},1}(x)\sim x^{-1/2}$ for $x\to\infty$. We plot in
Fig. 19(c,d) $\langle m^{2}\rangle/\varepsilon^{1/3}$ and
$\langle\dot{m}^{2}\rangle/\varepsilon^{5/3}$ as a function of
$N/\varepsilon^{-2}$, which is proportional to the rescaled system size
$N/N_{1}$. As expected, the different curves corresponding to different values
of $\varepsilon$ collapse for moderate up to intermediate values of $N$.

We now turn to the second crossover $N\sim N_{2}$ between the intermediate-
and large-$N$ regimes. For $N_{1}\ll N\ll N_{2}$, one has $\langle
m^{2}\rangle\sim\varepsilon^{-2/3}N^{-1/2}$ while for $N\gg N_{2}$, $\langle
m^{2}\rangle\sim\varepsilon$. Balancing the two contributions thus gives
$N_{2}\sim\varepsilon^{-10/3}$ (note that $N_{2}\gg N_{1}$). The same argument
holds for $\langle\dot{m}^{2}\rangle$: for $N_{1}\ll N\ll N_{2}$, one finds
$\langle\dot{m}^{2}\rangle\sim\varepsilon^{2/3}N^{-1/2}$ and for $N\gg N_{2}$,
$\langle\dot{m}^{2}\rangle\sim\varepsilon^{7/3}$, which also gives
$N_{2}\sim\varepsilon^{-10/3}$. These scaling behaviors of $\langle
m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ can be encompassed into two
scaling functions

$\langle
m^{2}\rangle=\varepsilon\,\mathcal{F}_{m,2}(\varepsilon^{10/3}N),\quad\langle\dot{m}^{2}\rangle=\varepsilon^{7/3}\,\mathcal{F}_{\dot{m},2}(\varepsilon^{10/3}N),$
(95)

with asymptotic behaviors
$\mathcal{F}_{m,2}(x)\sim\mathcal{F}_{\dot{m},2}(x)\sim x^{-1/2}$ for $x\to
0$, while both functions go to constant values for $x\to\infty$. In Fig.
19(e,f), we plot $\langle m^{2}\rangle/\varepsilon$ and
$\langle\dot{m}^{2}\rangle/\varepsilon^{7/3}$ as a function of the rescaled
system size $N/\varepsilon^{-10/3}\sim N/N_{2}$. As expected, for different
values of $\varepsilon$, the different curves collapse for intermediate up to
large values of $N$.

###  VI.3 Crossover between intermediate and large $N$ regimes

We reported above three distinct scaling regimes, separated by $N_{1}$ and
$N_{2}$, for the observables $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ as a function of the system size $N$ for
$\mu=\mu_{t}(\varepsilon)$, and we identified the scalings with $\varepsilon$
of the crossover sizes $N_{1}$ and $N_{2}$. We now investigate in more details
the behavior of the observables in the crossover regimes, focusing on the
effect of the variable $\mu$ close to the transition value $\mu_{t}$ for a
fixed $\varepsilon$. We numerically find that in the moderate- and
intermediate-N regimes, the two observables only weakly depend on $\mu$,
unlike for large-$N$ values. Thus, we now focus on the dependence of the
observables $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ on $\mu$
for a fixed $\varepsilon$, in the crossover regime between intermediate- and
large-$N$ values ($N\sim N_{2}$).

####  VI.3.1 Influence of $\mu$ on the second crossover regime

Figure 20: (a) $\langle m^{2}\rangle$ and (b) $\langle\dot{m}^{2}\rangle$ as a
function of the system size $N$, for
$(\mu-\mu_{c})/\mu_{c}\in[8.310,8.340,8.348,8.350]\times 10^{-7}$ from darker
to lighter colors. Parameters: $J_{1}=J_{2}=0.5$ and $\varepsilon=10^{-4}$.
The transition takes place at $(\mu_{t}-\mu_{c})/\mu_{c}\approx 8.349\times
10^{-7}$. Order parameters are evaluated numerically in the same way as in
Fig. 19. The black dashed line correspond to $\langle
x(m,H)\rangle_{\mathrm{approx}}$ defined in Eq. (96).

In the large-$N$ limit, the observables $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ are discontinuous at the transition, whereas for
moderate and intermediate-$N$ values, they do not depend much on the value of
$\mu$. In Fig. 20, $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ are
plotted as a function of system size $N$ for different values of $\mu$ along
the transition, at fixed $\varepsilon$. Like in Sec. V, we observe jumps in
$\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$ in the ferromagnetic
phase ($\mu<\mu_{t}$) at a finite system size, which takes place for higher
values of $N$ when approaching the transition. The two values of $\langle
m^{2}\rangle$ in the different phases are nonzero whereas
$\langle\dot{m}^{2}\rangle$ goes from a nonzero value (in the oscillating
phase) to a value decreasing as $N^{-1}$. The main difference with Sec. V is
that the jump in $\langle m^{2}\rangle$ is now much more pronounced because of
its different dependence on $\varepsilon$. Indeed, we had in Sec. V that
$\langle m^{2}\rangle\sim\varepsilon$ in both phases, whereas now $\langle
m^{2}\rangle\sim\varepsilon$ in the oscillating phase and $\langle
m^{2}\rangle\sim\varepsilon^{1/3}$ in the ferromagnetic phase, leading for
small $\varepsilon$ to a strong mismatch of $\langle m^{2}\rangle$ between the
two phases.

####  VI.3.2 Approximation in terms of phase coexistence

Here again, the observed behaviors can be given a simple interpretation in
terms of finite-size phase coexistence and metastability. Hence, as in Sec. V,
we introduce a simple decomposition of average values into contributions of
each phases. For any quantity $x(m,H)$, we introduce the approximate average
value obtained by taking into account the statistical weight of each phase,

$\langle x(m,H)\rangle_{\rm approx}=\frac{\langle x\rangle_{F}+\langle
x\rangle_{O}C_{2}\sqrt{N}e^{-N(\bar{f}_{2}-\bar{f}_{1})}}{1+C_{2}\sqrt{N}e^{-N(\bar{f}_{2}-\bar{f}_{1})}}$
(96)

where the ‘pure-state’ averages $\langle x\rangle_{F}$ and $\langle
x\rangle_{O}$ are respectively obtained from the ferromagnetic state large
deviation function $f_{F}(H)$ given in Eq. (85), and from the oscillating
state large deviation function $f_{O}(H)$ given in Eq. (86); $\bar{f}_{i}$ is
the minimum of $f_{i}(H)$ ($i=1$, $2$). The constant $C_{2}$, whose expression
is similar to that of the constant $C_{1}$ given in Eq. (75), but with
different expressions for $f_{F}$ and $f_{O}$ from the ones found in Sec. V,
now becomes

$C_{2}=\frac{-g(m_{0},0)}{(D_{22}+6v_{4}m_{0}^{8}D_{11})}\sqrt{\frac{6\pi
v_{4}}{bv_{0}}}m_{0}^{4}\,.$ (97)

In Sec. V, we could not have a local expression of $f_{O}(H)$ in the
oscillating phase, so we used fitting parameters for the coefficients
$\tilde{c}$ and $H^{*}$. Here, both $f_{F}$ and $f_{O}$ are known
analytically. The only quantity which is not known analytically and is a
fitted parameter, obtained from the numerical evaluation of $f(H)$, is
$\bar{f}_{2}-\bar{f}_{1}$.

This decomposition is plotted for different values of $\mu$ in Fig. 20 in
black dashed lines. The jump in $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ is well described by this simple decomposition.
When increasing $N$, the influence of the oscillating phase dominates until
$f(H^{*})\sim N^{-1}$.

####  VI.3.3 Comparison with stochastic simulations

Figure 21: Comparison with stochastic simulations of the spin model. (a)
$\langle m^{2}\rangle$ and (b) $\langle\dot{m}^{2}\rangle$ obtained from
numerical simulations as a function of the system size $N$, for
$(\mu-\mu_{c})/\mu_{c}\in[2,3,4]\times 10^{-3}$ (from darker to lighter
colors). Parameters: $J_{1}=J_{2}=0.5$ and $\varepsilon=5\times 10^{-2}$.

We plot in Fig. 21 $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$
computed from stochastic simulations for different system sizes and different
$\mu$, with $\varepsilon=(T_{c}-T)/T_{c}=5\times 10^{-2}$. We observe the
different expected behaviors described above. At large $N$, a jump of $\langle
m^{2}\rangle$ is observed between the ferromagnetic phase (high values of
$\langle m^{2}\rangle\sim\varepsilon^{1/3}$ for low values of $\mu-\mu_{c}$)
and the oscillating phase (low values of $\langle m^{2}\rangle\sim\varepsilon$
for higher values of $\mu-\mu_{c}$). A jump of $\langle\dot{m}^{2}\rangle$ is
also observed: $\langle\dot{m}^{2}\rangle$ is constant and of order
$\varepsilon^{7/3}\sim 10^{-3}$ for sufficiently high values of $\mu-\mu_{c}$,
while it decreases as $N^{-1}$ for lower values of $\mu-\mu_{c}$.

Similarly to Sec. V, we are able to describe the qualitative behavior of the
observables $\langle m^{2}\rangle$ and $\langle\dot{m}^{2}\rangle$. However,
for the restricted range of values of $\varepsilon$ and $N$ accessible in
stochastic simulations, we are not able to reach a quantitative agreement with
analytical predictions obtained in the small-$\varepsilon$, large-$N$ limit.
We now briefly discuss the effect of considering finite values of
$\varepsilon$ and $N$.

Figure 22: (a) $P_{N}(m,\dot{m})$ obtained from stochastic simulations for
$N=2\times 10^{4}$. (b) Trajectories ($m(t),\dot{m}(t))$ in the phase space
$(m,\dot{m})$ obtained in the deterministic limit. The color corresponds to
$v(m,\dot{m})^{-1}$ where $v(m,\dot{m})$ is the local velocity, as the local
density $p(m,\dot{m})\sim v(m,\dot{m})^{-1}$. The black red dots are added for
visual purpose and correspond to the ferromagnetic points. Parameters:
$J_{1}=J_{2}=0.5$, $\varepsilon=5\times 10^{-2}$ and
$(\mu-\mu_{c})/\mu_{c}=3.981\times 10^{-3}$.

In Fig. 22(a), we plot an example of $P_{N}(m,\dot{m})$ obtained from
stochastic simulations for $\varepsilon=10^{-1}$ and $N=5\times 10^{5}$, for a
value of $\mu$ where both the limit cycle and the ferromagnetic points are
linearly stable in the deterministic limit. We observe important discrepancies
with the theory described in Sec. VI.1, which assumed a small elliptic limit
cycle in the center with uniform probability along the cycle, and
ferromagnetic points outside the cycle. From a dynamical viewpoint, a non-
uniform probability along the cycle in an ensemble approach means that an
individual system goes along the cycle at a non-uniform speed. Here, as shown
in Fig. 22(a), the limit cycle is hardly visible, meaning the probability is
strongly non-uniform along the cycle. To understand this result, we plot in
Fig. 22(b) trajectories, obtained in the deterministic limit, in the phase
space $(m,\dot{m})$ where the color codes for $v(m,\dot{m})^{-1}$, with
$v(m,\dot{m})=[\left(dm/dt\right)^{2}+\left(d\dot{m}/dt\right)^{2}]^{1/2}$ the
local speed on the cycle. The quantity $v(m,\dot{m})^{-1}$ is proportional the
local probability density $p(m,\dot{m})$ along the limit cycle. We observe
that the limit cycle is not elliptic, and that the speed $v(m,\dot{m})$ is far
from uniform along the cycle. The probability density $p(m,\dot{m})$ is much
higher close to the ferromagnetic points, in qualitative agreement with
stochastic simulations. This discrepancy with theoretical predictions derived
in Sec. VI.1 comes from both the small-$\varepsilon$ and large-$N$
approximations made to obtain analytical results. The small-$\varepsilon$
approximation gives an elliptic limit cycle, and the large-$N$ one gives a
constant speed along the limit cycle [see Sec. IV.3 and Sec. V].

###  VI.4 Entropy production

Figure 23: Entropy production $\Sigma$ vs. $N$, for
$(\mu-\mu_{c})/\mu_{c}\in[6,7.2,8.0,8.2,8.3,8.330,8.343,8.348,8.352]\times
10^{-7}$ from darker to lighter colors. The entropy production is evaluated
numerically from Eq. (80). Inset: $N_{\Sigma}$, the value of $N$ where
$\Sigma$ is maximal, as a function of $\mu-\mu_{t}$, with
$(\mu_{t}-\mu_{c})/\mu_{c}=8.352\times 10^{-7}$. A slope $-1.3$ is indicated
by a dashed line, as a guide to the eye. Parameters: $J_{1}=J_{2}=0.5$ and
$\varepsilon=10^{-4}$.

Similarly to Sec. V, we discuss the behavior of the entropy production Eq.
(80) with system size at the transition, which is plotted in Fig. 23 for
different values of $\mu$ across the transition ($\mu\approx\mu_{t}$). For low
values of $N$, the entropy production increases as $N^{1/6}$ as
$\langle\dot{m}^{2}\rangle\sim N^{-5/6}$. For larger values of $N$ and for
$\mu>\mu_{t}$, $\Sigma\sim N$ as expected in an oscillating phase. For
$\mu<\mu_{t}$, as in Sec. V, we observe that $\Sigma$ increases like in the
oscillating phase before steeply decreasing to a constant value. This behavior
is similar to the one observed in Sec. V and is a consequence of the proximity
of the oscillating phase. We again denote as $N_{\Sigma}$ the value of $N$
when $\Sigma$ is maximal. In the inset of Fig. 23, we plot $N_{\Sigma}$ as a
function of the distance to the transition $\mu-\mu_{t}$, showing an
approximate power-law divergence of $N_{\Sigma}$ when $\mu-\mu_{t}\to 0$. As a
result, getting closer to the transition, the drop of $\Sigma$ takes place for
larger $N$ and thus becomes bigger, since it eventually decays to
approximately the same asymptotic large-$N$ value for all $\mu-\mu_{t}$.

###  VI.5 Comment on the continuous transition for $J_{2}=\pm 2+J_{1}$

In this section, we investigated the properties of the discontinuous
transition of Type II between the ferromagnetic phase and the oscillating
phase taking place for $J_{1}=J_{2}$ where $v_{1}(T_{c},\mu_{c})=0$. As seen
in Fig. 2, the condition $v_{1}(T_{c},\mu_{c})=0$ is also satisfied for
$J_{2}=\pm 2+J_{1}$, and it would thus be natural to also study the transition
in this case. However, in the deterministic limit we find for $J_{2}=\pm
2+J_{1}$ that the transition between the ferromagnetic phase and the
oscillating phase is of a different type: the transition is continuous, and
the ferromagnetic points turn into a limit cycle with a infinite period at the
transition. We now comment further on the difference between the two cases
$J_{1}=J_{2}$ and $J_{2}=\pm 2+J_{1}$ and why the method presented in this
paper does not allow for a characterization of continuous transitions between
ferromagnetic points and a limit cycle with infinite period.

In Sec. VI.1, we computed the large deviation locally using Eq. (37) and
assumed that the large deviation function is continuous in $m=m_{l}$ in order
to obtain the large deviation function in all three areas, and thus in the
whole plane ($m,\dot{m}$). However, an important issue is that at the point
where the different areas meet, $m=m_{l}$, we have $V^{\prime}(m_{l})=0$ and
$g(m_{l},0)\neq 0$. Thus, the assumption $V^{\prime}(m)\gg\dot{m}g(m,\dot{m})$
made to split the different orders of Eq. (LABEL:eq:phi:quadrat) is not valid
close to $m=m_{l}$. Therefore, we expect to get corrections close to $m_{l}$
which are not taken into account here. When the ferromagnetic points and the
limit cycle are well separated and far from $m_{l}$, which is true for
$J_{1}=J_{2}$ as we have $m\sim\varepsilon^{1/6}$ for the ferromagnetic phase
and $m\sim\varepsilon^{1/2}$ for the oscillating phase, the corrections close
to $m_{l}$ do not affect the qualitative behavior at the transition, and thus
the method presented in this section describes well the phase transition.
However, if one or both of them are close to $m_{l}$, corrections, that are
not taken into account in this paper and which would require a different
approach, are necessary.

For $J_{2}=\pm 2+J_{1}$, we find that the $m^{6}$-term in $V(m)$ is
independent of $\varepsilon$ (whereas for $J_{1}=J_{2}$ it scales as
$\varepsilon$) and thus is enough to describe the ferromagnetic points, which
turn out to scale as $m\sim\varepsilon^{1/2}$ (as seen by balancing the terms
in $\varepsilon m^{4}$ and $m^{6}$). If we blindly apply the method described
in this section, we find a limit cycle with $m\sim\varepsilon^{1/2}$ in
between the ferromagnetic points, which have the same scaling. The limit cycle
is very close to the point $m=m_{l}$ where $V^{\prime}(m_{l})=0$ and where the
theory breaks downs. Furthermore, from the definition of the period of the
limit cycle Eq. (44), we find that when $H^{*}\approx V(m_{l})$, the period is
very large, and diverges when $H^{*}=V(m_{l})$. Thus, by applying the method
without enough care, we would still recover some qualitative properties of the
transition: the ferromagnetic points and the limit cycle are very close to
each other and the period of the limit cycle is very large. Yet, this would
not be a correct description of the transition as we would find a
discontinuous transition instead of a continuous one as observed numerically.

##  VII Conclusion

We have considered a mean-field spin model with a dynamics breaking detailed
balance due to the non-reciprocal couplings between spins and auxiliary
dynamic fields. The presence of ferromagnetic interactions between spins on
one side, and between dynamic fields on the other side, allows for the
presence of both ferromagnetic and spontaneously oscillating phases. We have
characterized in details the transition between these two phases, showing that
it is discontinuous with the coexistence of both ferromagnetic and oscillating
states, one state being stable and the other one metastable. The relative
stability of both states is determined by a large deviation function,
generalizing the Landau free energy, that we evaluated explicitly in different
cases thanks to a perturbative framework. Two main scenarios are discussed,
whether the ferromagnetic points turn out to be inside or outside the limit
cycle. In addition, we found that the entropy production is peaked as a
function of system size, leading to a maximally dissipative system for an
optimal finite system size.

A natural generalization of this work may be to try to extend these results
beyond mean-field, by considering finite-dimensional systems, with the goal to
formulate a Ginzburg-Landau theory extending the present Landau framework
based on a large deviation principle. Such a theory might then be amenable to
a renormalization group treatment, extending the results of [3, 4] which
considered the synchronization of coupled oscillators. Here, we have started
from more basic ingredients, in the sense that the microscopic degrees of
freedom of the model (i.e., the spins and dynamic fields) do not oscillate in
the absence of interaction. Connecting these types of models to previous
results obtained on the synchronization transition is thus an interesting
challenge for future work.

###### Acknowledgements.

L. G. acknowledges funding from the French Ministry of Higher Education and
Research.

##  Appendix A Derivation of the deterministic evolution equations

In this appendix, we derive the deterministic evolution equations Eqs. (6) and
(7) from the microscopic spin and field dynamics. The dynamics of the system
is determined from the master equation [see Eq. (5)]. As the average is
defined as $\langle x\rangle=\sum_{\mathcal{C}}x(\mathcal{C})P(\mathcal{C})$,
we find after rearranging terms,

$\displaystyle d_{t}\langle
m\rangle=\sum_{\mathcal{C}}P(\mathcal{C})\sum_{\mathcal{C}^{\prime}\neq\mathcal{C}}\left[m(\mathcal{C}^{\prime})-m(\mathcal{C})\right]W(\mathcal{C}^{\prime}|\mathcal{C}),$
(98) $\displaystyle d_{t}\langle
h\rangle=\sum_{\mathcal{C}}P(\mathcal{C})\sum_{\mathcal{C}^{\prime}\neq\mathcal{C}}\left[h(\mathcal{C}^{\prime})-h(\mathcal{C})\right]W(\mathcal{C}^{\prime}|\mathcal{C}),$
(99)

with the shorthand notation $d_{t}\equiv d/dt$. From a configuration
$\mathcal{C}=\\{s_{1},...,s_{N},h_{1},...,h_{N}\\}$ with magnetization
$m=N^{-1}\sum_{i=1}^{N}s_{i}$ and average field $h=N^{-1}\sum_{i=1}^{N}h_{i}$,
there are $(1\pm m)N/2$ possibilities to flip a spin $s_{i}=\pm 1$ and $(1\pm
h)N/2$ possibilities to flip a field $h_{i}=\pm 1$. For a flip of a spin
$s_{i}=\pm 1$, $m(\mathcal{C}^{\prime})-m(\mathcal{C})=\mp 2/N$ and for a flip
of a field $h_{i}=\pm 1$, $h(\mathcal{C}^{\prime})-h(\mathcal{C})=\mp 2/N$.
Thus, using the definition of the transition rates given in the main text [Eq.
(2)], we find:

$\displaystyle d_{t}\langle m\rangle$
$\displaystyle=\big{\langle}-m+\tanh[\beta(J_{1}m+h)]\,\big{\rangle},$ (100)
$\displaystyle d_{t}\langle h\rangle$
$\displaystyle=\big{\langle}-h+\tanh[\beta(J_{2}h+(1-\mu)m)]\,\big{\rangle}.$
(101)

Assuming that the law of large numbers applies in the limit $N\to\infty$, $m$
and $h$ obey the following deterministic equations:

$\displaystyle d_{t}m$ $\displaystyle=-m+\tanh[\beta(J_{1}m+h)],$ (102)
$\displaystyle d_{t}h$ $\displaystyle=-h+\tanh[\beta(J_{2}h+(1-\mu)m)].$ (103)

These deterministic equations can be used to determine the macroscopic phase
when a single solution exists for given values of the control parameters
$\beta=T^{-1}$ and $\mu$. When two solutions exist, the most stable one has to
be determined from the large deviation function approach, as explained in the
main text.

##  Appendix B Values of the different functions and coefficients of the model

In this appendix, we give the values of the different functions and
coefficients introduced in the main text. The function $Y(m,\dot{m})$
introduced in Eq. (8) as $d_{t}\dot{m}=Y(m,\dot{m})$ has been split into a
$\dot{m}$-independent part, $V^{\prime}(m)=Y(m,0)$ and a $\dot{m}$-dependent
part $\dot{m}g(m,\dot{m})=Y(m,\dot{m})-Y(m,0)$. From Eqs. (6) and (7), we
have:

$\displaystyle Y(m,\dot{m})=$ $\displaystyle\beta J_{1}m+(-1+\beta
J_{1})\dot{m}-\tanh^{-1}(m+\dot{m})+\beta\tanh[J_{2}\tanh^{-1}(m+\dot{m})+\beta(1-\mu-
J_{1}J_{2})m]$ (104)
$\displaystyle+(m+\dot{m})^{2}\left[\tanh^{-1}(m+\dot{m})-\beta\tanh[J_{2}\tanh^{-1}(m+\dot{m})+\beta(1-\mu-
J_{1}J_{2})m]-\beta J_{1}(m+\dot{m})\right]\,,$ $V^{\prime}(m)=-\beta
J_{1}m+\beta
J_{1}m^{3}+(1-m^{2})\tanh^{-1}(m)-\beta(1-m^{2})\tanh[J_{2}\tanh^{-1}(m)+\beta(1-\mu-
J_{1}J_{2})m]\,.$ (105)

In the main text, we introduced the following expansions of $V(m)$ and
$g(m,\dot{m})$ [see Eqs. (16) and (17)],

$V(m)=\frac{v_{0}}{2}m^{2}+\frac{v_{1}}{4}m^{4},$ (106)

and

$g(m,\dot{m})=a_{0}\varepsilon-a_{1}m^{2}-a_{2}m\dot{m}-a_{3}\dot{m}^{2}.$
(107)

The coefficients appearing in these expansions are given by

$\displaystyle v_{0}=(\mu-1)/T^{2}+(1-J_{1}/T)(1-J_{2}/T),$ (108)
$\displaystyle\begin{aligned}
v_{1}=&-2/3+(2J_{2}+3J_{1})/3T-(\mu-1+J_{1}J_{2})/T^{2}\\\
&-(\mu-1-J_{2}T+J_{1}J_{2})^{3}/3T^{4},\end{aligned}$ (109)
$\displaystyle\varepsilon=(T_{c}-T)/T_{c},$ (110) $\displaystyle
a_{0}=2T_{c}/T,$ (111) $\displaystyle\begin{aligned}
a_{1}&=-2+(2J_{2}+J_{2}^{3}+3J_{1})/T\\\
&+J_{2}(-1+J_{1}J_{2}+\mu)^{2}/T^{3}\\\
&-2(1+J_{2}^{2})(-1+J_{1}J_{2}+\mu)/T^{2},\end{aligned}$ (112)
$\displaystyle\begin{aligned} a_{2}&=-2+2\beta J_{2}+\beta J_{2}^{3}+3\beta
J_{1}\\\ &+\beta^{2}(1+J_{2}^{2})(-1+J_{1}J_{2}+\mu),\end{aligned}$ (113)
$\displaystyle a_{3}=-2/3+(2J_{2}+J_{2}^{3}+3J_{1})/3T.$ (114)

In addition, we introduced in Eq. (31) the coefficients $D_{11}(m,\dot{m})$,
$D_{12}(m,\dot{m})=D_{21}(m,\dot{m})$ and $D_{22}(m,\dot{m})$, which read:

$\displaystyle D_{11}=1-m(m+\dot{m}),$ (115) $\displaystyle
D_{12}=[1-m(m+\dot{m})][-1+\beta J_{1}(1-(m+\dot{m})^{2})]$ (116)
$\displaystyle\begin{aligned} D_{22}&=[1-m(m+\dot{m})][-1+\beta
J_{1}(1-(m+\dot{m})^{2})]^{2}\\\
&+[1-h(h+\dot{h})]\beta^{2}[1-(m+\dot{m})^{2}]^{2}\end{aligned}$ (117)

with

$h=T\tanh^{-1}(m+\dot{m})-J_{1}m$ (118)

and

$\dot{h}(m,\dot{m})=T\frac{Y(m,\dot{m})+\dot{m}}{1-(m+\dot{m})^{2}}-J_{1}\dot{m}.$
(119)

In the main text, we use the coefficients evaluated at $m=0$ and $\dot{m}=0$,
which simplify to:

$\displaystyle D_{11}(0,0)=1,$ (120) $\displaystyle D_{12}(0,0)=-1+J_{1}/T,$
(121) $\displaystyle D_{22}(0,0)=1/T^{2}+(J_{1}/T-1)^{2}.$ (122)

##  Appendix C Entropy production

In steady-state, the entropy production can be identified with the entropy
flux, which is defined from the microscopic configurations
$\mathcal{C}=\\{s_{1},...,s_{N},h_{1},...,h_{N}\\}$ as

$\Sigma=\sum_{\mathcal{C},\mathcal{C}^{\prime}}W(\mathcal{C}^{\prime}|\mathcal{C})P(\mathcal{C})\,\ln\frac{W(\mathcal{C}^{\prime}|\mathcal{C})}{W(\mathcal{C}|\mathcal{C}^{\prime})}\,.$
(123)

Note that Eq. (123) is equivalent to the definition (80) given in the main
text, up to a symmetrization of expression (123). To compute the entropy
production, we aim at changing the sum over the configurations into integrals
over $m$ and $\dot{m}$. We thus replace $\sum_{\mathcal{C}}$ by $\int
dmd\dot{m}\sum_{\mathcal{C}\in S(m,\dot{m})}$, where $S(m,\dot{m})$ denotes
the ensemble of configurations $\mathcal{C}$ with $m(\mathcal{C})=m$ and
$\dot{m}(\mathcal{C})=\dot{m}$. We transform the integral over
$\mathcal{C}^{\prime}$ into a sum over all possible transitions. We recall
that spin reversals are labelled with $k=1$ and field reversals with $k=2$,
keeping track of the sign $\sigma=\pm$ of the variable prior to reversal. We
denote as $W_{k}^{\sigma}$ the coarse-grained transition rates given in Eq.
(26) and $n_{k}^{\sigma}$ the fraction of possible transitions,

$n_{1}^{\pm}=\frac{1}{2}(1\pm m),\quad n_{2}^{\pm}=\frac{1}{2}(1\pm h).$ (124)

We find

$\Sigma=N\sum_{k,\sigma}\left\langle
W_{k}^{\sigma}(\mathbf{x})\ln\left[\\!\frac{W_{k}^{\sigma}(\mathbf{x})n_{k}^{-\sigma}(\mathbf{x}+\frac{\sigma\mathbf{d}_{k}}{N})}{W_{k}^{-\sigma}(\mathbf{x}+\frac{\sigma\mathbf{d}_{k}}{N})n_{k}^{\sigma}(\mathbf{x})}\right]\\!\right\rangle$
(125)

with the shorthand notation $\mathbf{x}=(m,\dot{m})$. To further simplify
notations, we make the dependence on $\mathbf{x}$ implicit in what follows. At
leading order in $N$, we find

$\Sigma=N\sum_{k}\left\langle\left(W_{k}^{+}-W_{k}^{-}\right)\ln\left[\frac{W_{k}^{+}n_{k}^{-}}{W_{k}^{-}n_{k}^{+}}\right]\right\rangle.$
(126)

We define

$\Sigma_{k}=\left\langle\left(W_{k}^{+}-W_{k}^{-}\right)\ln\left[\frac{W_{k}^{+}n_{k}^{-}}{W_{k}^{-}n_{k}^{+}}\right]\right\rangle.$
(127)

For $k=1$, we have $W_{1}^{+}-W_{1}^{-}=\dot{m}/2$ and

$\ln\left[\frac{W_{k}^{+}n_{k}^{-}}{W_{k}^{-}n_{k}^{+}}\right]=2\tanh^{-1}(m+\dot{m}),$
(128)

thus we find that

$\Sigma_{1}=N\left\langle\dot{m}\tanh^{-1}(m+\dot{m})\right\rangle.$ (129)

In the main text, we showed that close to a transition, one generically has a
scaling behavior $\dot{m}\sim\varepsilon^{\alpha}$ (with $\alpha>0$ an
exponent depending on the specific transition considered) and $\langle
m\dot{m}\rangle=0$. Thus, keeping the lowest order in $\varepsilon$, we find

$\Sigma_{1}/N=\langle\dot{m}^{2}\rangle.$ (130)

For $k=2$, we recall that we have
$h(m,\dot{m})=-J_{1}m+T\tanh^{-1}[m+\dot{m}]$ and we introduce
$\dot{h}=-h+\tanh[\beta(J_{2}h+(1-\mu)m)]$, the equivalent of $\dot{m}$ for
the fields variables $h_{i}$, given in Eq. (119) as a function of $m$ and
$\dot{m}$. We find $\Sigma_{2}=N\langle\dot{h}\tanh^{-1}(h+\dot{h})\rangle$.
As spins and fields play symmetric roles, one also has $\langle
h\dot{h}\rangle=0$ close to a transition. Using that
$Y(m,\dot{m})=-V^{\prime}(m)$ and $\langle V^{\prime}(m)\dot{m}\rangle=0$ at
the lowest order in $\varepsilon$, we find

$\frac{\Sigma_{2}}{N}=\langle
T^{2}V^{\prime}(m,\dot{m})^{2}+(T-J_{1})^{2}\dot{m}^{2}\rangle.$ (131)

Gathering contributions, one finds

$\frac{\Sigma}{N}=\left[1+(T-J_{1})^{2}\right]\left\langle\dot{m}^{2}\right\rangle+T^{2}\left\langle
V^{\prime}(m)^{2}\right\rangle$ (132)

at leading order in $\varepsilon$ and $N$.

##  Appendix D Moderate-$N$ approximation of Sec. VI

In this appendix, we give the exact expressions of $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ for moderate values of $N$ for the Type-II
discontinuous transition of Sec. VI. For intermediate values of $N$, only high
values of $V(m)$ contribute, thus we consider that $V(m)=v_{4}m^{10}/10$.
Using Eq. (37), we find

$f(H)=AH^{6/5}\,,$ (133)

with

$A=\frac{9\Gamma\left(\frac{3}{10}\right)^{2}a_{1}}{(2^{6}5^{4}v_{4})^{1/5}\sqrt{\pi}\Gamma\left(\frac{1}{10}\right)}\,.$
(134)

Thus, using the definition of $\langle m^{2}\rangle$ and
$\langle\dot{m}^{2}\rangle$ from Eq. (47) we find:

$\displaystyle\langle m^{2}\rangle=\frac{c_{m}}{(a_{1}v_{4})^{1/6}}N^{-1/6},$
(135)
$\displaystyle\langle\dot{m}^{2}\rangle=c_{\dot{m}}\frac{v_{4}^{1/6}}{a_{1}^{5/6}}N^{-5/6},$
(136)

with

$\displaystyle
c_{m}=\frac{5^{1/3}\Gamma\left(\frac{3}{10}\right)\Gamma\left(\frac{2}{3}\right)}{3\times
2^{2/5}\pi}\left(\frac{\Gamma\left(\frac{13}{10}\right)}{\Gamma\left(\frac{11}{10}\right)}\right)^{5/6}\\!\left(\frac{\Gamma\left(\frac{9}{5}\right)}{\Gamma\left(\frac{8}{5}\right)}\right)^{1/6}\\!\approx
0.53,$ (137) $\displaystyle
c_{\dot{m}}=\frac{5^{2/3}\Gamma\left(\frac{3}{5}\right)\Gamma\left(\frac{4}{3}\right)\left(\Gamma\left(\frac{11}{10}\right)\Gamma\left(\frac{9}{5}\right)\right)^{5/6}}{2\sqrt{\pi}\Gamma\left(\frac{8}{5}\right)^{11/6}\Gamma\left(\frac{13}{10}\right)^{5/6}}\approx
1.33.$ (138)

## References

  * Acebrón  _et al._ [2005] J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005).
  * Gupta  _et al._ [2014] S. Gupta, A. Campa, and S. Ruffo, J. Stat. Mech.: Theor. Exp. , R08001 (2014).
  * Risler  _et al._ [2004] T. Risler, J. Prost, and F. Jülicher, Phys. Rev. Lett. 93, 175702 (2004).
  * Risler  _et al._ [2005] T. Risler, J. Prost, and F. Jülicher, Phys. Rev. E 72, 016130 (2005).
  * Nicolis [1986] G. Nicolis, Rep. Prog. Phys. 49, 873 (1986).
  * Kamino  _et al._ [2017] K. Kamino, Y. Kondo, A. Nakajima, M. Honda-Kitahara, K. Kaneko, and S. Sawai, Proc. Natl. Acad. Sci. USA 114, E4149 (2017).
  * Wang and Tang [2019] S.-W. Wang and L.-H. Tang, Nat. Commun. 10, 5613 (2019).
  * Saha  _et al._ [2020] S. Saha, J. Agudo-Canalejo, and R. Golestanian, Phys. Rev. X 10, 041009 (2020).
  * You  _et al._ [2020] Z. You, A. Baskaran, and M. C. Marchetti, Proc. Natl. Acad. Sci. USA 117, 19767 (2020).
  * Cao  _et al._ [2015] Y. Cao, H. Wang, Q. Ouyang, and Y. Tu, Nat. Phys. 11, 772 (2015).
  * Nguyen  _et al._ [2018] B. Nguyen, U. Seifert, and A. C. Barato, J. Chem. Phys. 149, 045101 (2018).
  * Aufinger  _et al._ [2022] L. Aufinger, J. Brenner, and F. C. Simmel, Nat. Commun. 13, 2852 (2022).
  * Devailly  _et al._ [2015] C. Devailly, C. Crauste-Thibierge, A. Petrosyan, and S. Ciliberto, Phys. Rev. E 92, 052312 (2015).
  * Andrae  _et al._ [2010] B. Andrae, J. Cremer, T. Reichenbach, and E. Frey, Phys. Rev. Lett. 104, 218102 (2010).
  * Duan  _et al._ [2019] D. Duan, B. Niu, and J. Wei, Chaos, Solitons and Fractals 123, 206 (2019).
  * Gualdi  _et al._ [2015] S. Gualdi, J.-P. Bouchaud, G. Cencetti, M. Tarzia, and F. Zamponi, Phys. Rev. Lett. 114, 088701 (2015).
  * Yi  _et al._ [2015] S. D. Yi, S. K. Baek, G. Chevereau, and E. Bertin, J. Stat. Mech.: Theor. Exp. , P11001 (2015).
  * Collet  _et al._ [2016] F. Collet, M. Formentin, and D. Tovazzi, Phys. Rev. E. 94, 042139 (2016).
  * De Martino and Barato [2019] D. De Martino and A. C. Barato, Phys. Rev. E 100, 062123 (2019).
  * Dai Pra  _et al._ [2020] P. Dai Pra, M. Formentin, and P. Guglielmo, J. Stat. Phys. 179, 690 (2020).
  * Crawford [1991] J. D. Crawford, Rev. Mod. Phys. 63, 991 (1991).
  * Fei  _et al._ [2018] C. Fei, Y. Cao, Q. Ouyang, and Y. Tu, Nat. Commun. 9, 1434 (2018).
  * Gaspard [2002] P. Gaspard, J. Chem. Phys. 117, 8905 (2002).
  * Barato and Seifert [2016] A. C. Barato and U. Seifert, Phys. Rev. X 6, 041053 (2016).
  * Barato and Seifert [2017] A. C. Barato and U. Seifert, Phys. Rev. E 95, 062409 (2017).
  * Oberreiter  _et al._ [2022] L. Oberreiter, U. Seifert, and A. C. Barato, Phys. Rev. E 106, 014106 (2022).
  * Remlein  _et al._ [2022] B. Remlein, V. Weissmann, and U. Seifert, Phys. Rev. E 105, 064101 (2022).
  * Sagués  _et al._ [2007] F. Sagués, J. M. Sancho, and J. García-Ojalvo, Rev. Mod. Phys. 79, 829 (2007).
  * Xu  _et al._ [2020] H.-Y. Xu, Y.-P. Luo, J.-W. Wu, and M.-C. Huang, Physica D 411, 132612 (2020).
  * Crochik and Tomé [2005] L. Crochik and T. Tomé, Phys. Rev. E 72, 057103 (2005).
  * Xiao  _et al._ [2008] T. J. Xiao, Z. Hou, and H. Xin, J. Chem. Phys. 129, 114506 (2008).
  * Xiao  _et al._ [2009] T. Xiao, Z. Hou, and H. Xin, J. Phys. Chem. B 113, 9316 (2009).
  * Barato and Hinrichsen [2012] A. C. Barato and H. Hinrichsen, J. Phys. A: Math. Theor. 45, 115005 (2012).
  * Tomé and de Oliveira [2012] T. Tomé and M. J. de Oliveira, Phys. Rev. Lett. 108, 020601 (2012).
  * Noa  _et al._ [2019] C. E. F. Noa, P. E. Harunari, M. J. de Oliveira, and C. E. Fiore, Phys. Rev. E 100, 012104 (2019).
  * Martynec  _et al._ [2020] T. Martynec, S. H. L. Klapp, and S. A. M. Loos, New J. Phys. 22, 093069 (2020).
  * Seara  _et al._ [2021] D. S. Seara, B. B. Machta, and M. P. Murrell, Nat. Commun. 12, 392 (2021).
  * Le Bellac [1992] M. Le Bellac,  _Quantum and Statistical Field Theory_ (Oxford Science Publications, Oxford, 1992).
  * Meibohm and Esposito [2022] J. Meibohm and M. Esposito, Phys. Rev. Lett. 128, 110603 (2022).
  * Holtzman and Raz [2022] R. Holtzman and O. Raz, Commun. Phys. 5, 280 (2022).
  * Aron and Chamon [2020] C. Aron and C. Chamon, SciPost Phys. 8, 074 (2020).
  * Guislain and Bertin [2023] L. Guislain and E. Bertin, Phys. Rev. Lett. 130, 207102 (2023).
  * Fruchart  _et al._ [2021] M. Fruchart, R. Hanai, P. Littlewood, and V. Vitelli, Nature 592, 363 (2021).
  * Martin  _et al._ [2023] D. Martin, D. Daniel Seara, Y. Avni, M. Fruchart, and V. Vitelli, (2023), arXiv:2307.08251 .
  * Collet [2014] F. Collet, J. Stat. Phys. 157, 1301 (2014).
  * Sinelschikov  _et al._ [2023] D. Sinelschikov, A. Poggialini, M. F. Abbate, and D. De Martino, Emergence of collective self-oscillations in minimal lattice models with feedback (2023), arXiv:2306.01823 .
  * Touchette [2009] H. Touchette, Physics Reports 478 (2009).
  * Knessl  _et al._ [1985] C. Knessl, B. J. Matkowsky, Z. Schuss, and C. Tier, SIAM J. Appl. Math. 46, 1006 (1985).
  * Graham and Tél [1987] R. Graham and T. Tél, Phys. Rev. A 35, 1328 (1987).
  * Graham and Tél [1984] R. Graham and T. Tél, Journal of Statistical Physics 35, 729 (1984).
  * Graham [1989] R. Graham, Journal of Statistical Physics 54, 1207 (1989).
  * Gillespie [2007] D. T. Gillespie, Annu. Rev. Phys. Chem. 58, 35 (2007).
  * [53] See Supplemental Material at XXX.
  * Schnackenberg [1976] J. Schnackenberg, Rev. Mod. Phys. 48, 571 (1976).
  * Gaspard [2004] P. Gaspard, J. Stat. Phys. 117, 599 (2004).

Supplementary Information: Discontinuous phase transition from ferromagnetic
to oscillating states in a nonequilibrium mean-field spin model

## Corrections to the large deviation function for $\mu=\mu_{l}$

Stochastic simulations are made for values of $\varepsilon$ and $N$ for which
the assumptions of small $\varepsilon$ and large $N$ are not fully satisfied.