1 Introduction

Emerging periodic behavior in complex systems with a large number of interacting units is a commonly observed phenomenon in neuroscience [13], ecology [22], socioeconomics [4, 23] and life sciences in general. From a mathematical standpoint, when modeling such a phenomenon it is natural to consider large families of microscopic identical units evolving through noisy interacting dynamics. Single units have no natural tendency to behave periodically and oscillations are rather an effect of self-organization as they emerge in the macroscopic limit when the number of particles tends to infinity. Within this modeling framework mean field models have received much attention due to their analytical tractability. For mean field models one can obtain in closed form the stochastic dynamics of a single unit in the limit of infinite many units; the time evolution of the associated distribution is a, possibly infinite-dimensional, dynamical system. Periodic trajectories of this dynamical system correspond to the emergence of self-organized periodic oscillations in the interacting system; this will be referred to with the term self-sustained periodic behavior. One of the goals of the mathematical theory in this field is to understand which types of microscopic interactions and mechanisms can lead to or enhance the above type self-organization. Among others, we cite noise [10, 19, 21], dissipation in the interaction potential [1, 6, 7, 9], delay in the transmission of information and/or frustration in the interaction network [8, 11, 20]. In particular, in [11] the authors consider non-Markovian dynamics, studying systems of interacting nonlinear Hawkes processes for modeling neurons.

Although not proved in general, a strong belief in the literature is that, at least for Markovian dynamics, self-sustained periodic behavior cannot emerge if one does not introduce some time-irreversible phenomenon in the dynamics, as it is the case in all the above cited works (see e.g. [3, 6, 15]). The model treated here, in which the limit dynamics is still reversible with respect to the stationary distribution around which cycles emerge (see Remark 1 below), suggests that this paradigm could be false for the non-Markovian case.

Specifically, we give numerical and mathematical evidence of the emergence of self-sustained periodic behavior of the empirical magnetization in a mean field spin system related to the Curie–Weiss model, which happens to belong to the following universality class: it features the presence of a unique stable zero mean phase for values of the parameters corresponding to high temperatures, the emergence of periodic oscillations in an intermediate range of the parameter values, and a subsequent ferromagnetic ordered phase for increasingly lower temperatures. Our recipe consists in replacing the Poisson distribution of the spin-flip times with another renewal process, thus making the individual spin dynamics non-Markovian. In details, we consider the distribution of the interarrival times to have tails proportional to \(e^{-t^{\gamma +1}}\), for \(\gamma =1,2\). Then, we introduce an interaction among the spins via a time rescaling depending on the magnetization of the system. The specific choice of interarrival time distribution makes the computations developed in Sect. 4 as easiest as possible (to our knowledge), allowing for an explicit characterization of the discrete spectrum of the linearization of the limiting Fokker–Planck equation. A question which can arise naturally is whether similar results can be found for other classes of waiting times. Although we do not have a general answer to this, we want to remark that simulations with different distributions highlighted the same characteristics (e.g. Gamma distribution, tails proportional to \(e^{-t^{\gamma +1}}\) with \(\gamma \in {\mathbb {R}}\) such that \(\gamma \ge 1\)). All the working examples we considered feature exponentially or super-exponentially decaying tails. On the other hand, we have examples of polynomial tails (e.g. inverse Gamma distribution) where no oscillatory behavior was experienced.

The paper is organized as follows: in Sect. 2 we describe the model and the results obtained. In particular, before introducing the model (Sect. 2.2), we start by recalling basic facts about the Curie–Weiss model and its phase transition (Sect. 2.1); we then proceed with the results on the propagation of chaos (Sect. 2.3), and on the linearization of the Fokker–Planck equation around a zero mean equilibrium, for two different choices of renewal dynamics (Sect. 2.4). Notably, we determine the discrete spectrum of the linearized operator in terms of the zeros of two holomorphic functions. Numerical investigations of the discrete spectrum, studied as a function of the interaction parameters, show the presence of a Hopf bifurcation: a pair of complex conjugate eigenvalues of the linearized Fokker–Planck equation around a stationary solution crosses the complex plane imaginary axis (see Sects. 2.4.1 and  2.4.2). Such a bifurcation of the mean field limit is reflected in the behavior of the finite particle system. Indeed, numerical simulations of the empirical magnetization in Sect. 3 confirm the transition from an incoherent state to self-organized rhythmic oscillations. Section 4 contains the proofs of the results of Sect. 2.

2 Model and Results

2.1 Motivation

As we mentioned above, the model we consider can be seen as a proper modification of the Curie–Weiss dynamics. When we refer to the latter, we mean a spin-flip type Markovian dynamics for a system of N interacting spins \(\sigma _i \in \left\{ -1,1\right\} \), \(i =1,\ldots ,N\). Such dynamics is reversible with respect to the equilibrium Gibbs probability measure on the space of configurations \(\left\{ -1,1\right\} ^N\),

$$\begin{aligned} P_{N,\beta }(\varvec{\sigma }) := \frac{1}{Z_N(\beta )} \exp \left[ -\beta H_N(\varvec{\sigma })\right] , \end{aligned}$$
(1)

with \(\varvec{\sigma } := (\sigma _1,\ldots ,\sigma _N) \in \left\{ -1,1\right\} ^N\), \(\beta > 0\) (ferromagnetic case), \(Z_N(\beta )\) is a normalizing constant, and \(H_N\) is the Hamiltonian

$$\begin{aligned} H_N(\varvec{\sigma }) := - \frac{1}{2N}\left( \sum _{i=1}^N \sigma _i\right) ^2. \end{aligned}$$
(2)

Denote also the empirical magnetization as \(m^N := \frac{1}{N}\sum _{i=1}^N \sigma _i\). Note that the distribution (1) gives higher probability to the configurations with minimal energy, which by (2) are the ones where the individual spins are aligned in the same state. The equilibrium model undergoes a phase transition tuned by the interaction parameter \(\beta >0\), which can be recognized by proving a Law of Large Numbers for the equilibrium empirical magnetization

$$\begin{aligned} \text {Law}(m^N) \xrightarrow {N \rightarrow +\infty } {\left\{ \begin{array}{ll} \delta _0, \ \ \ &{} \,\text { if } \beta \le 1,\\ \frac{1}{2}\delta _{+m_\beta } + \frac{1}{2}\delta _{-m_\beta }, \ \ \ &{} \, \text { if } \beta > 1, \end{array}\right. } \end{aligned}$$
(3)

where \(m_\beta > 0\) is the so-called spontaneous magnetization [2, 12]. When we turn to the dynamics, different choices can be made in order to satisfy the above-mentioned reversibility with respect to (1). The prototype is a continuous-time spin-flip dynamics defined in terms of the infinitesimal generator L, applied to a function \(f : \left\{ -1,1\right\} ^N \rightarrow {\mathbb {R}}\),

$$\begin{aligned} Lf(\varvec{\sigma }) = \sum _{i=1}^N e^{-\beta \sigma _i m^N}\left[ f(\varvec{\sigma }^i)-f(\varvec{\sigma })\right] , \end{aligned}$$
(4)

where \(\varvec{\sigma }^i \in \left\{ -1,1\right\} ^N\) is obtained from \(\varvec{\sigma }\) by flipping the i-th spin. Dynamics (4) induces a continuous-time Markovian evolution for the empirical magnetization process \(m^N(t)\), which is given in terms of a generator \({\mathcal {L}}^N\) applied to a function \(g :[-1,1] \rightarrow {\mathbb {R}}\):

$$\begin{aligned} {\mathcal {L}}^Ng(m) = N \frac{1+ m}{2} e^{-\beta m}\left[ g\left( m-\frac{2}{N}\right) -g(m)\right] + N \frac{1-m}{2}e^{\beta m}\left[ g\left( m+\frac{2}{N}\right) - g(m)\right] . \end{aligned}$$
(5)

It is easy to obtain the weak limit of the sequence of processes \(\big (m^N(t)\big )_{t \ge 0}\), by studying the uniform convergence of the generator (5) as \(N \rightarrow +\,\infty \) (see e.g. [14]). The limit process \((m(t))_{t \ge 0}\) is deterministic and solves the Curie–Weiss ODE

$$\begin{aligned} {\left\{ \begin{array}{ll} {\dot{m}}(t) = 2\sinh (\beta m(t)) - 2 m(t) \cosh (\beta m(t)),\\ m(0) = m_0 \in [-1,1]. \end{array}\right. } \end{aligned}$$
(6)

The presence of the phase transition highlighted in (3) can be recognized as well in the out-of-equilibrium dynamical model (6). Indeed, studying the long-term behavior of (6), one finds that:

  • for \(\beta \le 1\), (6) possesses a unique stationary solution, globally attractive, constantly equal to 0;

  • for \(\beta > 1\), 0 is still stationary but it is unstable; two other symmetric stationary locally attractive solutions, \(\pm m_\beta \), appear: the two non-zero solutions to \(m = \tanh (\beta m)\). The dynamics m(t) gets attracted for \(t \rightarrow +\infty \) to the polarized stationary state which has the same sign as the initial magnetization \(m_0\).

Another concept which we refer to in what follows is that of renewal process, a generalization of the Poisson process. We identify a renewal process with the sequence of its interarrival times (also commonly referred to as sojourn times or waiting times in the literature) \(\left\{ T_n\right\} _{n=1}^\infty \), i.e. the holding times between the occurrences of two consecutive events. The Poisson process is characterized by having independent and identically distributed interarrival times, where each \(T_i\) is exponentially distributed. In particular, the memoryless property \({\mathbb {P}}(T_i> s+t | T_i> t) = {\mathbb {P}}(T_i > s)\), holds for any \(s,t \ge 0\). The interarrival times of a renewal process are still independent and identically distributed, but their distribution is not required to be exponential. We recall that a continuous-time homogeneous Markov chain can be identified by a Poisson process, modeling the jump times, and a stochastic transition matrix, identifying the possible arrival states at each jump time. Due to the lack of the memoryless property, when one replaces the Poisson process in the definition of the spin-flip dynamics with a more general renewal process, the resulting evolution is thus non-Markovian. In the literature, the associated dynamics is referred to as semi-Markov process, first introduced by Lévy in [18].

2.2 The Dynamics

Glauber dynamics as the Curie–Weiss dynamics in (4) can be constructed in two stages. First one system of independent spin-flips: at the times of a Poisson process of intensity 1, the spin in a given site flips; different sites have independent Poisson processes. Then, the interaction can be introduced as a spin-dependent time scale in the waiting time for updates. In this section we illustrate this procedure, and generalize it by replacing the Poisson Process with a more general renewal process.

For the moment we focus on a single spin \(\sigma (t) \in \left\{ -1,1\right\} \). If driven by a Poisson process of intensity 1, its dynamics has infinitesimal generator

$$\begin{aligned} {\mathcal {L}}f(\sigma ) = f(-\sigma ) - f(\sigma ), \end{aligned}$$
(7)

with \(f : \left\{ -1,1\right\} \rightarrow {\mathbb {R}}\). If the Poisson process is replaced by a renewal process, the spin dynamics is not Markovian. In what follows, we refer to the resulting dynamics as a spin-valued renewal process, that is an example of q two-state semi-Markov process. We can associate a Markovian description to the latter: define y(t) as the time elapsed since the last spin-flip occured up to time t. Suppose that the waiting times \(\tau \) (interchangeably referred to as interarrival times) of the renewal satisfy

$$\begin{aligned} {\mathbb {P}}(\tau > t) = \varphi (t), \end{aligned}$$
(8)

for some smooth function \(\varphi : [0,+\infty ) \rightarrow {\mathbb {R}}\). Then, the pair \((\sigma (t),y(t))_{t \ge 0}\) is Markovian with generator

$$\begin{aligned} {\mathcal {L}}f(\sigma ,y) = \frac{\partial f}{\partial y}(\sigma ,y) + F(y)[f(-\sigma ,0) - f(\sigma ,y)], \end{aligned}$$
(9)

for \(f : \left\{ -1,1\right\} \times {\mathbb {R}}^+ \rightarrow {\mathbb {R}}\), with

$$\begin{aligned} F(y) := -\frac{\varphi '(y)}{\varphi (y)}. \end{aligned}$$
(10)

This is equivalent to say that the couple \((\sigma (t), y(t))_{t \ge 0}\) evolves according to

$$\begin{aligned} {\left\{ \begin{array}{ll} (\sigma (t), y(t)) \mapsto (-\sigma (t),0), \ \ \ \text { with rate } F(y(t)), \\ dy(t) = dt, \ \ \ \text { otherwise.} \end{array}\right. } \end{aligned}$$
(11)

Note that the dynamics in (11) is well defined whenever F is continuous and nonnegative. Moreover, expression (10) for the jump rate follows by observing that, for an interarrival time \(\tau \) of the jump process \(\sigma (t)\), we have

$$\begin{aligned} {\mathbb {P}}(\sigma (t+h) = -\sigma | \sigma (t) = \sigma ) = 1 - {\mathbb {P}}(\tau> t+h | \tau > t) = 1 - \frac{\varphi (t+h)}{\varphi (t)}, \end{aligned}$$

for any \(h > 0\). Observe that when the \(\tau \)’s are exponentially distributed \(F(y) \equiv 1\), so we get back to dynamics (7). Dynamics (9) can be perturbed by allowing the distribution of the waiting time for a spin-flip to depend on the current spin value \(\sigma \); the simplest way is to model this dependence as a time scaling:

$$\begin{aligned} {\mathbb {P}}(\tau > t | \sigma ) = \varphi (a(\sigma ) t). \end{aligned}$$
(12)

Under this distribution for the waiting times the generator of \((\sigma (t),y(t))_{t \ge 0}\) becomes:

$$\begin{aligned} {\mathcal {L}}f(\sigma ,y) = \frac{\partial f}{\partial y}(\sigma ,y) + a(\sigma ) F(a(\sigma )y)[f(-\sigma ,0) - f(\sigma ,y)]. \end{aligned}$$

On the basis of what seen above, it is rather simple to define a system of mean field interacting spins with non-exponential waiting times. For a collection of N pairs \((\sigma _i(t),y_i(t))_{i=1,\ldots ,N}\), we set \(m^N(t) := \frac{1}{N}\sum _{i=1}^N \sigma _i(t)\) to be the magnetization of the system at time t. The interacting dynamics is

$$\begin{aligned} {\left\{ \begin{array}{ll} (\sigma _i(t), y_i(t)) \mapsto (-\sigma _i(t),0), \ \ \ \text { with rate } \ \ F\left( y_i(t) e^{-\beta \sigma _i(t) m^N(t)}\right) e^{-\beta \sigma _i(t) m^N(t)}, \\ dy_i(t) = dt, \ \ \ \text { otherwise.} \end{array}\right. } \end{aligned}$$
(13)

where \(\beta >0\) is a parameter tuning the interaction between the particles.

Denoting \(\varvec{\sigma }:= (\sigma _1,\ldots , \sigma _N) \in \left\{ -1,1\right\} ^N\), \(\varvec{y} := (y_1,\ldots ,y_N) \in ({\mathbb {R}}^+)^N\), \(m^N := \frac{1}{N}\sum _{i=1}^N \sigma _i\), the associated infinitesimal generator is

$$\begin{aligned} {\mathcal {L}}^N f(\varvec{\sigma },\varvec{y}) = \sum _{i=1}^N \frac{\partial f}{\partial y_i}(\varvec{\sigma },\varvec{y}) + \sum _{i=1}^N F\left( y_i e^{-\beta \sigma _i m^N}\right) e^{-\beta \sigma _i m^N}\left[ f(\varvec{\sigma }^i,\varvec{y}^i) - f(\varvec{\sigma },\varvec{y})\right] ,\nonumber \\ \end{aligned}$$
(14)

where \(\varvec{\sigma }^i\) is obtained from \(\varvec{\sigma }\) by flipping the i-th spin, while \(\varvec{y}^i\) by setting to zero the i-th coordinate. The additional factor \(e^{-\beta \sigma _i(t) m^N(t)}\) in the jump rate in (13) follows from the observation we made in (12) and the definition of \(F(y) = -\frac{\varphi '(y)}{\varphi (y)}\). Note that, for \(F \equiv 1\), we retrieve the Curie–Weiss dynamics (4) for the spins.

2.3 Propagation of Chaos

The macroscopic limit and propagation of chaos for the above class of models should be standard, although some difficulties may arise for general choices of F not globally Lipschitz. For computational reasons which will be made clear below, we focus on the case \(F(y) = y^\gamma \), for \(\gamma \in {\mathbb {N}}\), which corresponds to considering, in the single spin model (9), the tails of the distribution of the interarrival times to be \(\varphi (t) \propto e^{-\frac{t^{\gamma +1}}{\gamma +1}}\).

When \(F(y) = y^\gamma \), (13) becomes

$$\begin{aligned} {\left\{ \begin{array}{ll} (\sigma _i(t), y_i(t)) \mapsto (-\sigma _i(t),0), \ \ \ \text { with rate } \ \ y_i^{\gamma }(t) e^{-(\gamma +1)\beta \sigma _i(t) m^N(t)}, \\ dy_i(t) = dt, \ \ \ \text { otherwise.} \end{array}\right. } \end{aligned}$$
(15)

As for the Curie–Weiss model, dynamics (15) is subject to a cooperative-type interaction: the spin-flip rate is larger for particles which are not aligned with the majority. Assuming propagation of chaos, at the macroscopic limit \(N \rightarrow +\infty \) the representative particle \((\sigma (t), y(t))\) has a mean-field dynamics

$$\begin{aligned} {\left\{ \begin{array}{ll} (\sigma (t), y(t)) \mapsto (-\sigma (t),0), \ \ \ \text { with rate } \ \ y^{\gamma }(t) e^{-(\gamma +1)\beta \sigma (t) m(t)}, \\ dy(t) = dt, \ \ \ \text { otherwise,} \end{array}\right. } \end{aligned}$$
(16)

with \(m(t) = {\mathbb {E}}[\sigma (t)]\). To this dynamics we can associate (see [17]) the non-linear infinitesimal generator \({\mathcal {L}}(m(t))\) with depends on the value of m(t) and acts on a function \(f(\sigma ,y)\) as follows:

$$\begin{aligned} {\mathcal {L}}(m(t)) f(\sigma ,y) = \frac{\partial f}{\partial y}(\sigma ,y) + y^\gamma e^{-(\gamma +1)\beta \sigma m(t)}\left[ f(-\sigma ,0) - f(\sigma ,y)\right] , \end{aligned}$$
(17)

where the non-linearity is due to the dependence of the generator on m(t), a function of the joint law at time t of the processes \((\sigma (t),y(t))\). In Sect. 4 we study rigorously the well-posedeness of the pre-limit and limit dynamics and the propagation of chaos. The main result is collected in the following

Theorem 1

(Propagation of chaos) Fix \(\gamma \in {\mathbb {N}}\), and let \(T>0\) be the final time in (15) and (16). Assume that \((\sigma _i(0), y_i(0))_{i=1,\ldots ,N}\) are \(\mu _0\)-chaotic for some probability distribution \(\mu _0\) on \(\left\{ -1,1\right\} \times {\mathbb {R}}^+\). Then, the sequence of empirical measures \((\mu _t^N)_{t \in [0,T]}\) converges in distribution on the path space (in the sense of weak convergence of probability measures) to the deterministic law \((\mu _t)_{t \in [0,T]}\) of the unique solution to Eq. (16) with initial distribution \(\mu _0\).

2.4 Local Analysis of the Fokker–Planck

In this section we illustrate the results on the local analysis of the Fokker–Planck equation for the mean field limit dynamics (16) with \(\gamma = 1\) and \(\gamma = 2\). We adopt the following approach: we find a stationary solution of interest, then we linearize formally the dynamics around that equilibrium. The specific choice of the values of \(\gamma \) allows for an explicit characterization of the discrete spectrum of the linearization of the Fokker–Planck equation. Indeed, we compute the discrete spectrum of the associated linearized operator, which we show to be given by the zeros of an explicit holomorphic function \(H_{\beta ,\gamma }(\lambda )\). We then study numerically the character of the eigenvalues when \(\beta \) varies: for both \(\gamma =1,2\), we find that for all \(\beta < \beta _c(\gamma )\) all eigenvalues have negative real part; at \(\beta _c(\gamma )\) two eigenvalues are conjugate and purely imaginary, suggesting the possible presence of a Hopf bifurcation in the limit dynamics. These critical values of \(\beta \) are then compared to the ones estimated from the simulations of the finite particle system in Sect. 3.

The Fokker–Planck equation associated to (16) is a PDE describing the time evolution of the density function \(f(t,\sigma ,y)\) of the limit process \((\sigma (t),y(t))\). It is given by

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial }{\partial t} f(t,\sigma ,y) + \frac{\partial }{\partial y} f(t,\sigma ,y) + y^\gamma e^{-(\gamma + 1)\beta \sigma m(t)} f(t,\sigma ,y) = 0,\\ f(t,\sigma ,0) = \int _{0}^{+\infty } y^\gamma e^{(\gamma + 1)\beta \sigma m(t)} f(t,-\sigma ,y) dy, \\ m(t) = \int _{0}^{\infty } [f(t,1,y) - f(t,-1,y)]dy,\\ 1 = \int _{0}^{\infty } [f(t,1,y) + f(t,-1,y)]dy,\\ f(0,\sigma ,y) = f_0(\sigma ,y), \text { for } \sigma \in \left\{ -1,1\right\} , \ y \in {\mathbb {R}}^+. \end{array}\right. } \end{aligned}$$
(18)

We postpone a derivation of the Fokker–Planck equation to Sect. 4.1. Here we just spend a few words on the boundary integral condition - second line of System (18) - which is specific to our model and may not be clear at first sight. It is a mass-balance condition: heuristically, at each time t, the mass of the spins at state \((\sigma ,0)\) (i.e. \(f(t,\sigma ,0)\)) equals the spins just jumped from \((-\sigma ,y)\) weighted with their probabilities (i.e. \(y^\gamma e^{(\gamma + 1)\beta \sigma m(t)} f(t,-\sigma ,y) dy\)) integrated over all the possible jump times. While a general study of (18) is beyond the scope of this work, here we just observe that (18) can be seen as a system of two quasilinear PDEs (one for \(\sigma = 1\) and another for \(\sigma = -1\)), where the non-linearity enters in an integral form through m(t) in the exponent of the rate function. Moreover, the boundary integral condition in the second line poses additional challenges. Nevertheless, it is easy to exhibit a particular stationary solution to (18):

Proposition 1

The function

$$\begin{aligned} f^*(\sigma ,y) = \frac{1}{2\varLambda }e^{-\frac{y^{\gamma +1}}{\gamma +1}}, \end{aligned}$$
(19)

with \(\varLambda := \int _{0}^{+\infty } e^{-\frac{y^{\gamma +1}}{\gamma +1}}dy\), is a stationary solution to Sys. (18) with \(m = 0\).

Remark 1

Let \(g^*(\sigma )\) be the marginal of \(f^*(\sigma ,y)\) with respect to the first coordinate. Then, \(g^*(\sigma )\) is a stationary reversible distribution for the limit renewal process \((\sigma (t))_{t \ge 0}\). Indeed, by choosing \(\sigma (0) \sim g^*\), \(g^*(1) = g^*(-1) = \frac{1}{2}\), we have that \(m(t) \equiv 0\) and \((\sigma (t))_{t \ge 0}\) is a renewal process with interarrival times \(\tau \) such that \({\mathbb {P}}(\tau > t) \propto e^{-\frac{t^{\gamma +1}}{\gamma +1}}\) independently of the value of \(\sigma \), so its law is invariant by time reversal.

The linearization of the operator associated to Sys.  (18) around the equilibrium (19), yields the following eigen-system

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial }{\partial y} g(\sigma ,y) + y^\gamma g(\sigma ,y) - \frac{\beta \sigma k(\gamma + 1)}{2\varLambda }y^\gamma e^{-\frac{y^{\gamma +1}}{\gamma +1}} = -\lambda g(\sigma ,y),\\ g(\sigma ,0) = \frac{\beta \sigma k (\gamma + 1)}{2\varLambda } + \int _0^\infty g(-\sigma ,y) y^\gamma dy,\\ \int _0^\infty [g(\sigma ,y) + g(-\sigma ,y)]dy = 0, \ \ \ (\sigma ,y) \in \left\{ -1,1\right\} \times {\mathbb {R}}^+, \end{array}\right. } \end{aligned}$$
(20)

where \(k = 2\int _0^\infty g(1,y) dy\), and \(\varLambda = \int _0^\infty e^{-\frac{y^{\gamma +1}}{\gamma + 1}} dy\). For a formal derivation see Sect. 4.3.1. We work out the computations of the discrete spectrum of the linearized operator for the two cases \(\gamma = 1\), \(\gamma = 2\).

2.4.1 Case \(\gamma = 1\)

In this case, \(\varLambda = \sqrt{\frac{\pi }{2}}\), and the eigen-system (20) becomes

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial }{\partial y} g(\sigma ,y) + y g(\sigma ,y) + \lambda g(\sigma ,y) = \beta \sigma k \left( \sqrt{\frac{\pi }{2}}\right) ^{-1}y e^{-\frac{y^2}{2}} \\ g(\sigma ,0) = \beta \sigma k \left( \sqrt{\frac{\pi }{2}}\right) ^{-1} + \int _0^{\infty } y g(-\sigma ,y)dy,\\ \int _0^{\infty } [g(\sigma ,y) + g(-\sigma ,y)]dy = 0, \end{array}\right. } \end{aligned}$$
(21)

where \(k = 2\int _0^\infty g(1,y)dy\). The eigenvalues of (21) are given by the zeros of an explicit holomorphic function.

Proposition 2

The solutions in \(\lambda \in {\mathbb {C}}\) to (21) are the zeros of the holomorphic function

$$\begin{aligned} H_{\beta ,1}(\lambda ) := H_1(\lambda )\left[ -4\beta - \lambda ^3 \sqrt{\frac{\pi }{2}}\right] + \sqrt{2\pi }\lambda ^2 - 4\beta \lambda + 2\beta \sqrt{2\pi }, \end{aligned}$$
(22)

with

$$\begin{aligned} H_1(\lambda ) := \int _0^\infty e^{-\frac{y^2}{2}}e^{-\lambda y}dy. \end{aligned}$$
(23)

Moreover, it holds

$$\begin{aligned} H_1(\lambda ) = \sqrt{\frac{\pi }{2}}\sum _{m=0}^\infty \frac{\lambda ^{2m}}{(2m)!!} - \lambda \sum _{m=0}^\infty \frac{(2\lambda )^{2m}m!}{(2m + 1)!}\frac{1}{2^m}. \end{aligned}$$
(24)

Equation \(H_{\beta ,1}(\lambda ) = 0\) can be numerically investigated. We used the power series expansion (24) and a numerical root finding built-in function of the software Mathematica, specifically FindRoot, starting the search from different initial points of the complex plane and from different values of \(\beta \). Here we report the results:

  1. (1.1)

    we find two conjugate purely imaginary solutions to \(H_{\beta ,1}(\lambda ) = 0\), for \(\lambda = \pm \,\lambda _c(1):= \pm i (1.171)\) and

    $$\begin{aligned} \beta = \beta _c(1) := 0.769; \end{aligned}$$
    (25)
  2. (1.2)

    iterating the search around \((\beta _c(1),\lambda _c(1))\), the resulting complex eigenvalue goes from having a negative real part for \(\beta < \beta _c(1)\) to a positive real part for \(\beta > \beta _c(1)\);

  3. (1.3)

    no other purely immaginary solution \(\lambda = \pm \,i x\) is found for \(0 \le x \le 500\) and \(0 \le \beta \le 20 \);

  4. (1.4)

    for \(\beta < \beta _c(1)\) all the eigenvalues \(\lambda = ix + y\) are such that \(y < 0\). This was verified for \(-100\le x \le 100\), \(-100\le y \le 100\).

2.4.2 Case \(\gamma = 2\)

In this case the eigen-system is given by

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial }{\partial y} g(\sigma ,y) + y^2 g(\sigma ,y) +\lambda g(\sigma ,y) = \frac{3}{2\varLambda }\beta \sigma k y^2 e^{-\frac{y^3}{3}} ,\\ g(\sigma ,0) = \frac{3}{2\varLambda } \beta \sigma k + \int _0^{\infty } y^2 g(-\sigma ,y)dy,\\ \int _0^{\infty } [g(\sigma ,y) + g(-\sigma ,y)]dy = 0, \end{array}\right. } \end{aligned}$$
(26)

where \(\varLambda = \int _0^\infty e^{-\frac{y^3}{3}}dy= \frac{\varGamma \left( \frac{1}{3}\right) }{3^{2/3}}\), \(k = 2 \int _0^{\infty } g(1,y) dy\), and \(\varGamma (\cdot )\) is the Gamma function. Analogously to the case \(\gamma =1\), the eigenvalues of (26) can be computed as the zeros of an explicit holomorphic function.

Proposition 3

The solutions in \(\lambda \in {\mathbb {C}}\) to (26) are the zeros of the holomorphic function

$$\begin{aligned} \begin{aligned} H_{\beta ,2}&(\lambda ):= H_2(\lambda )\Bigg [12\beta - \lambda ^4 \varLambda + 6\beta \lambda \varLambda - 6\beta \lambda 3^{1/3}\varGamma (4/3) + 3\beta \lambda ^2 3^{2/3}\varGamma (5/3)\\&\quad \quad \quad - 6\beta \lambda ^2\frac{\varGamma (2/3)}{3^{1/3}}\Bigg ]+\Bigg [2\varLambda \lambda ^3 - 12\beta \varLambda + 12\beta \frac{\varGamma (2/3)}{3^{1/3}}\lambda - 6\beta \lambda ^2 \Bigg ], \end{aligned} \end{aligned}$$
(27)

with

$$\begin{aligned} H_2(\lambda ) := \int _0^\infty e^{-\lambda y} e^{-\frac{y^3}{3}} dy. \end{aligned}$$
(28)

Moreover, it holds

$$\begin{aligned} H_2(\lambda ) = \sum _{n=0}^{\infty } (-1)^n \frac{\lambda ^n}{n!} 3^{\frac{1}{3}(n-2)}\varGamma \left( \frac{n+1}{3}\right) . \end{aligned}$$
(29)

We solved numerically \(H_{\beta ,2}(\lambda ) = 0\) for different values of the parameters. Apart from being sensibly slower, it seems the numerical root finding for \(\gamma = 2\) suffers from numerical instability issues. This is why we were able to check the results for much smaller intervals in this case. Our results are the following:

  1. (2.1)

    we find two conjugate purely imaginary solutions \(\lambda = \pm \lambda _c(2) := \pm i(1.978)\) for

    $$\begin{aligned} \beta = \beta _c(2):= 0.362; \end{aligned}$$
    (30)
  2. (2.2)

    analogous to (1.2);

  3. (2.3)

    analogous to (1.3), verified for \(0 \le x \le 10\) and \(0\le \beta \le 5\);

  4. (2.4)

    analogous to (1.4), verified for \(-25 \le x \le 25\), \(-25 \le y \le 25\).

Fig. 1
figure 1

Simulation of the finite particle system dynamics (15) for \(\gamma =1\) (left) and \(\gamma = 2\) (right). The plot shows the time evolution of the empirical magnetization of \(N = 1500\) spins. In the picture is depicted the evolution of magnetization, with initial data \(\sigma _i(0) = 1\) for every \(i=1,\ldots ,N\). Simulations suggest that the zero mean trajectory is attractive and the magnetization goes to zero regardless of the initial datum

Fig. 2
figure 2

Simulation of the finite particle system dynamics (15) for \(\gamma =1\) (left) and \(\gamma = 2\) (right), with number of spins \(N = 1500\). Increasing the value of \(\beta \) oscillations of the empirical magnetization appear. Initially not very regular, the rhythmic behavior becomes evident for larger \(\beta \) as shown in Fig. 3 and 4

3 Finite Particle System Simulations

Fig. 3
figure 3

Above a critical value of \(\beta \) the dynamics defined in (15) is periodic. In the figure the time evolution of the pair \((m^N(t),y^N(t))\) is plotted for \(N=1500\) particles. Here \(m^N(t)\) is the empirical magnetization of the spins, whereas \(y^N(t)\) indicates the empirical mean of the \(y_i\)’s

Fig. 4
figure 4

Simulation of the finite particle system dynamics (15) for \(\gamma =1\) (left) and \(\gamma = 2\) (right), with number of spins \(N = 1500\). For \(\beta \) above a critical value the empirical magnetization has a clear periodic behavior

Fig. 5
figure 5

Increasing \(\beta \) further above the first critical value, the system (15) crosses a second critical value above which oscillations disappear, and the system magnetizes, i.e. the magnetization stabilizes to a non-zero value. The simulation in the figure shows the evolution of the empirical magnetization for a system comprised of \(N = 1500\) spins

We ran several simulations of the particle system with \(N = 1500\) spins for \(\gamma = 1, 2\). Simulations are in accordance with the above numerical results on the eigenvalues. In particular, above the critical values of \(\beta \) (see (25) and (30)) periodic behavior of the finite interacting particle system appears. Description of the evidences is the following:

  • For \(\beta \) small the magnetization goes to zero regardless of the initial datum (Fig. 1).

  • There is a critical \(\beta \) (around 0.75 for \(\gamma = 1\), around 0.35 for \(\gamma = 2\)) above which the magnetization starts oscillating. Close to the critical point oscillations do not look very regular (corrupted by noise?), but they soon become very regular if \(\beta \) is not too close to the critical value (see Fig. 2). We also made joint plots of the magnetization with the empirical mean of the \(y_i\)’s (Fig. 3). Periodic oscillations seems to emerge.

  • As \(\beta \) increases, the amplitude of the oscillation of the magnetization increases, while the period looks nearly constant (see Fig. 4). As \(\beta \) crosses another critical value (around 1.3 for \(\gamma = 1\), around 1.65 for \(\gamma = 2\)) oscillations disappear, and the system magnetizes, i.e. the magnetization stabilizes to a non-zero value, actually close to \(\pm 1\) (Fig. 5).

  • The oscillations are lasting for a wider interval of \(\beta \)’s for \(\gamma =2\) (from \(\beta \approx 0.35\) until \(\beta \approx 1.65\)) than \(\gamma = 1\) (from \(\beta \approx 0.75\) until \(\beta \approx 1.3\)). The period is instead smaller for \(\gamma = 2\) than for \(\gamma = 1\).

  • For both \(\gamma = 1,2\), the appearance of oscillations in the particle system dynamics does not seem to depend on the initial data. This behavior of the finite volume system may reflect the presence of a global Hopf bifurcation in the mean field limit dynamics.

4 Proofs

4.1 Formal Derivation of the Fokker–Planck Equation

We give a sketch of the derivation of the Fokker–Planck equation (18). Let the process \((\sigma (t),y(t))_{t\ge 0}\) be distributed at time t according to some distribution \(dm(t,\sigma ,y) \in {\mathcal {P}}\left( \left\{ -1,1\right\} \times [0,+\infty )\right) \), which we assume to be regular in y, i.e. absolutely continuous with respect to the Lebesgue measure on \([0,+\infty )\), with a smooth density \(f(t,\sigma , y)dy\). Let \(h : \left\{ -1,1\right\} \times [0,+\infty ) \rightarrow {\mathbb {R}}\) be a smooth function, and compute

$$\begin{aligned} \frac{d}{dt}{\mathbb {E}}\left[ h(\sigma (t),y(t))\right] = {\mathbb {E}}\left[ {\mathcal {L}}(m(t))h(\sigma (t),y(t))\right] , \end{aligned}$$

with \(m(t) := \int _{[0,+\infty )}[f(t,1,y) - f(t,-1,y)]dy\). Under our regularity assumptions, the above expression is equivalent to

$$\begin{aligned}&\int _{[0,+\infty )} h(1,y)\frac{\partial }{\partial t}f(t,1,y)dy + \int _{[0,+\infty )} h(-1,y)\frac{\partial }{\partial t}f(t,-1,y)dy \\&\quad = \int _{[0,+\infty )} \frac{\partial }{\partial y} h(1,y) f(t,1,y)dy + \int _{[0,+\infty )} \frac{\partial }{\partial y} h(-1,y) f(t,-1,y)dy \\&\quad \quad + \int _{[0,+\infty )}y^{\gamma } e^{-(\gamma +1)\beta m(t)}\left[ h(-1,0) - h1,y)\right] f(t,1,y)dy \\&\quad \quad + \int _{[0,+\infty )}y^{\gamma } e^{(\gamma +1)\beta m(t)}\left[ h(1,0) - h(-1,y)\right] f(t,-1,y)dy. \end{aligned}$$

Integrating by parts in the first two integrals on the right hand side of the equality, and regrouping the terms, one finds

$$\begin{aligned}&\int _{[0,+\infty )}h(1,y)\left[ \frac{\partial }{\partial t} f(t,1,y) + \frac{\partial }{\partial y} f(t,1,y) + y^\gamma e^{-(\gamma + 1) \beta m(t)} f(t,1,y)\right] dy\\&\quad + \int _{[0,+\infty )}h(-1,y)\left[ \frac{\partial }{\partial t} f(t,-1,y) + \frac{\partial }{\partial y} f(t,-1,y) + y^\gamma e^{(\gamma + 1) \beta m(t)} f(t,-1,y)\right] dy\\&\quad + h(1,0)\left[ f(t,1,0) - \int _{[0,+\infty )} y^\gamma e^{(\gamma + 1)\beta m(t)} f(t,-1,y)dy\right] \\&\quad + h(-1,0)\left[ f(t,-1,0) - \int _{[0,+\infty )} y^\gamma e^{-(\gamma + 1)\beta m(t)} f(t,1,y)dy\right] = 0, \end{aligned}$$

which, by the arbitrariness of h implies the first two equations in (18); the definition of m(t) and the requirement for f to be a density function over \(\left\{ -1,1\right\} \times [0,+\infty )\) complete the equations in (18).

4.2 Proof of Theorem 1

Here we prove rigorously a propagation of chaos property for the N-particle dynamics to its mean-field limit, for any \(\gamma \in {\mathbb {N}}\). Actually, we establish the proofs for \(\gamma = 1\), where the rates enjoy globally Lipschitz properties, and then we generalize them to any \(\gamma \in {\mathbb {N}}\) in Remark 2. The generalization to non-Lipschitz rates is possible because of the a-priori bound on the variables \(y_i\)’s which, by definition, are such that \(0 \le y_i \le T\), where \(T < \infty \) is the final time horizon of the dynamics. For the convenience of the reader, we write again the dynamics

$$\begin{aligned} {\left\{ \begin{array}{ll} (\sigma _i(t), y_i(t)) \mapsto (-\sigma _i(t),0), \ \ \ \text { with rate } \ \ y_i^{\gamma }(t) e^{-(\gamma + 1)\beta \sigma _i(t) m^N(t)}, \\ dy_i(t) = dt, \ \ \ \text { otherwise,} \end{array}\right. } \end{aligned}$$
(31)

and the mean-field version

$$\begin{aligned} {\left\{ \begin{array}{ll} (\sigma (t), y(t)) \mapsto (-\sigma (t),0), \ \ \ \text { with rate } \ \ y^{\gamma }(t) e^{-(\gamma +1) \beta \sigma (t) m(t)}, \\ dy(t) = dt, \ \ \ \text { otherwise,} \end{array}\right. } \end{aligned}$$
(32)

with \(m(t) = {\mathbb {E}}[\sigma (t)]\). We represent both the microscopic and the macroscopic model as solutions of certain stochastic differential equations driven by Poisson random measures, in order to apply the results in [16]. As anticipated, in the proof we restrict to a finite interval of time [0, T].

To begin with, let us fix a filtered probability space \(\big ((\Omega , {\mathcal {F}},{\mathbb {P}}),({\mathcal {F}}_t)_{t \in [0,T]}\big )\) satisfying the usual hypotheses, rich enough to carry an independent and identically distributed collection \({\mathcal {N}}\), \(({\mathcal {N}}_i)_{i \in {\mathbb {N}}}\) of stationary Poisson random measures on \([0,T] \times \varXi \), with intensity measure \(\nu \) on \(\varXi := [0,+\infty )\) equal to the restriction of the Lebesgue measure on \([0, +\infty )\). For any N, consider the system of Itô-Skorohod equations

$$\begin{aligned} {\left\{ \begin{array}{ll} \sigma _i(t) = \sigma _i(0) + \int _0^t \int _{\varXi }f_1(\sigma _i(s^-), \xi , m^N(s^-), y_i(s^-)){\mathcal {N}}_i(ds,d\xi ),\\ y_i(t) = y_i(0) + t + \int _0^t \int _{\varXi }f_2(\sigma _i(s^-), \xi , m^N(s^-), y_i(s^-)){\mathcal {N}}_i(ds,d\xi ), \end{array}\right. } \end{aligned}$$
(33)

and the corresponding limit non-linear reference particle’s dynamics

$$\begin{aligned} {\left\{ \begin{array}{ll} \sigma (t) = \sigma (0) + \int _0^t \int _{\varXi }f_1(\sigma (s^-), \xi , m(s^-), y(s^-)){\mathcal {N}}(ds,d\xi ),\\ y(t) = y(0) + t + \int _0^t \int _{\varXi }f_2(\sigma (s^-), \xi , m(s^-), y(s^-)){\mathcal {N}}(ds,d\xi ). \end{array}\right. } \end{aligned}$$
(34)

The functions \(f_1, f_2 : \left\{ -1,1\right\} \times {\mathbb {R}}^+ \times [-1,1] \times {\mathbb {R}}^+ \rightarrow {\mathbb {R}}\), modeling the jumps of the process, are given by

$$\begin{aligned} f_1(\sigma , \xi , m, y) := - 2\sigma \mathbb {1}_{]0,\lambda [}(\xi ), \quad f_2(\sigma , \xi ,m,y):= -y\mathbb {1}_{]0,\lambda [}(\xi ), \end{aligned}$$
(35)

with \(\lambda := \lambda (\sigma ,m,y)\) being the rate function \(\lambda (\sigma ,m,y) = y^{\gamma } e^{-(\gamma +1) \beta \sigma m}\).

Proposition 4

For \(\gamma = 1\), Eqs. (33) and (34) possess a unique strong solution for \(t \in [0,T]\).

Proof

With the choices in (35), the well-posedeness of Eqs. (33) and (34) follows by Theorems 1.2 and 2.1 in [16]. Indeed, even though the function \(f_2\) is not globally Lipschitz continuous in y, the \(L^1\) Lipschitz assumption of the theorem still holds, by noting that

$$\begin{aligned}&\int _{\varXi } \Big |f_2(\sigma ,\xi ,m,y) - f_2({\tilde{\sigma }},\xi ,{\tilde{m}},{\tilde{y}})\Big | d\xi = \int _{\varXi } \Big | y \mathbb {1}_{]0,\lambda (\sigma ,m,y)[}(\xi ) - {\tilde{y}}\mathbb {1}_{]0,\lambda ({\tilde{\sigma }},{\tilde{m}},{\tilde{y}})[}(\xi )\Big |d\xi \\&\quad \le |y| \big |\lambda (\sigma ,m,y) - \lambda ({\tilde{\sigma }},{\tilde{m}},{\tilde{y}})\big | + \big |\lambda ({\tilde{\sigma }},{\tilde{m}},{\tilde{y}})\big | \big |y - {\tilde{y}}\big |\\&\quad \le |y| \big [\big |\lambda (\sigma ,m,y) - \lambda ({\tilde{\sigma }},{\tilde{m}},{\tilde{y}})\big | + \big |\lambda ({\tilde{\sigma }},{\tilde{m}},{\tilde{y}})\big ||y - {\tilde{y}}|\big ]\\&\quad \le CT \big [|m- {\tilde{m}}| + |y - {\tilde{y}}| + |\sigma - {\tilde{\sigma }}|\big ], \end{aligned}$$

where in the last step we have used that, by construction, the processes \(y_i(t) \le T\) for every \(t \in [0,T]\), so that the rates are a priori bounded and the Lipschitz properties of \(y e^{-2\beta \sigma m}\) for \((y,\sigma ,m) \in {\mathbb {R}}^+ \times \left\{ -1,1\right\} \times [-1,1]\). \(\square \)

Now, define the empirical measures \(\mu ^N := \frac{1}{N}\sum _{i=1}^N \delta _{(\sigma _i,y_i)}\), and their evaluation along the paths of (33),

$$\begin{aligned} \mu ^N_t := \frac{1}{N}\sum _{i=1}^N \delta _{(\sigma _i(t), y_i(t))}. \end{aligned}$$
(36)

The measures \((\mu _t^N)_{t \in [0,T]}\) can be viewed as random variables with values in \({\mathcal {P}}(D)\), the space of probability measures on D, where \(D:= {\mathcal {D}}\big ([0,T];\left\{ -1,1\right\} \times {\mathbb {R}}^+\big )\) is the space of \(\left\{ -1,1\right\} \times {\mathbb {R}}^+\)-valued càdlàg functions equipped with the Skorohod topology.

Proof of Theorem 1

Consider the i.i.d. processes \(({\tilde{\sigma }}_i(t), {\tilde{y}}_i(t))_{i=1,\ldots ,N}\), coupled with the N-particle dynamics \((\sigma _i(t),y_i(t))_{i=1,\ldots ,N}\),

$$\begin{aligned} {\left\{ \begin{array}{ll} {\tilde{\sigma }}_i(t) = {\tilde{\sigma }}_i(0) + \int _0^t \int _{\varXi }f_1({\tilde{\sigma }}_i(s^-), \xi , m(s^-), {\tilde{y}}_i(s^-)){\mathcal {N}}_i(ds,d\xi ),\\ {\tilde{y}}_i(t) = {\tilde{y}}_i(0) + t + \int _0^t \int _{\varXi }f_2({\tilde{\sigma }}_i(s^-), \xi , m(s^-), {\tilde{y}}_i(s^-)){\mathcal {N}}_i(ds,d\xi ), \end{array}\right. } \end{aligned}$$
(37)

with \(m(t) = {\mathbb {E}}[{\tilde{\sigma }}_i(t)]\). Let \(({\tilde{\mu }}_t^N)_{t \in [0,T]}\) be the empirical measure associated to (37). Clearly, one has \(({\tilde{\mu }}_t^N)_{t \in [0,T]} \rightarrow (\mu _t)_{t \in [0,T]}\) in the weak convergence sense (by a functional LLN, see [16] for e.g.) where \((\mu _t)_{t \in [0,T]}\) is the deterministic law of the unique solution to Eq. (16) with initial distribution \(\mu _0\). We are thus left to show

$$\begin{aligned} \varvec{d_1}\Big (\text {Law}\big ((\mu _t^N)_{t \in [0,T]}\big ),\text {Law}\big (({\tilde{\mu }}_t^N)_{t \in [0,T]}\big )\Big ) \xrightarrow {N \rightarrow +\infty } 0, \end{aligned}$$

with \(\varvec{d_1}\) being the 1-Wasserstein distance (which metrizes the weak convergence of probability measures) on \({\mathcal {P}}({\mathcal {P}}(D))\). Since (for instance see [9])

$$\begin{aligned} \varvec{d_1}\Big (\text {Law}\big ((\mu _t^N)_{t \in [0,T]}\big ),\text {Law}\big (({\tilde{\mu }}_t^N)_{t \in [0,T]}\big )\Big ) \le \frac{1}{N}\sum _{i=1}^N {\mathbb {E}}\big [d_{Sko}\big ((\sigma _i,y_i),({\tilde{\sigma }}_i,{\tilde{y}}_i)\big )\big ], \end{aligned}$$

with \(d_{Sko}\) the Skorohod metric on D, it is enough to show that

$$\begin{aligned} \frac{1}{N}\sum _{i=1}^N {\mathbb {E}}\Bigg [\sup _{t \in [0,T]} \Big (|\sigma _i(t) - {\tilde{\sigma }}_i(t)| + | y_i(t) - {\tilde{y}}_i(t)|\Big )\Bigg ] \xrightarrow {N \rightarrow +\infty } 0. \end{aligned}$$
(38)

The proof of (38) starts with the following inequalities which make use of the bounds for \(f_2\) in Proposition 4. In the subsequent inequalities the value of the constant C may change from line to line. We have:

$$\begin{aligned}&{\mathbb {E}}\left[ \sup _{s \in [0,t]}|y_i(s) - {\tilde{y}}_i(s)|\right] \le {\mathbb {E}}\Big [|y_i(0) - {\tilde{y}}_i(0)|\Big ] \\&\quad + C \int _0^t {\mathbb {E}}\Bigg [|m^N(s) - m(s)| + |y_i(s) - {\tilde{y}}_i(s)| + |\sigma _i(s) - {\tilde{\sigma }}_i(s)|\Bigg ]ds \le C \int _0^t {\mathbb {E}}\\&\qquad \Bigg [\sup _{r \in [0,s]}|m^N(r) - m(r)| + \sup _{r \in [0,s]}|y_i(r) - {\tilde{y}}_i(r)| + \sup _{r \in [0,s]} |\sigma _i(r) - {\tilde{\sigma }}_i(r)|\Bigg ]ds + C(N), \end{aligned}$$

with \(C(N) \xrightarrow {N \rightarrow +\infty } 0\) because of the chaoticity assumption on the initial datum. We proceed similarly for the \(\sigma _i\)’s, using the the Lipschitz continuity of \(f_1\), we obtain

$$\begin{aligned}&{\mathbb {E}}\left[ \sup _{s \in [0,t]}|\sigma _i(s) - {\tilde{\sigma }}_i(s)|\right] \le {\mathbb {E}}\Big [|\sigma _i(0) - {\tilde{\sigma }}_i(0)|\Big ] \\&\quad + C \int _0^t {\mathbb {E}}\Bigg [|m^N(s) - m(s)| + |y_i(s) - {\tilde{y}}_i(s)| + |\sigma _i(s) - {\tilde{\sigma }}_i(s)|\Bigg ]ds \le C \int _0^t {\mathbb {E}}\\&\qquad \Bigg [\sup _{r \in [0,s]}|m^N(r) - m(r)| + \sup _{r \in [0,s]}|y_i(r) - {\tilde{y}}_i(r)| + \sup _{r \in [0,s]} |\sigma _i(r) - {\tilde{\sigma }}_i(r)|\Bigg ]ds+ C(N). \end{aligned}$$

Denoting \({\tilde{m}}^N (t):= \frac{1}{N}\sum _{i=1}^N {\tilde{\sigma }}_i(t)\), we find

$$\begin{aligned}&{\mathbb {E}}\Bigg [\sup _{s \in [0,t]} |m^N(s) - m(s)|\Bigg ] \le {\mathbb {E}}\Bigg [\sup _{s \in [0,t]} |m^N(s) - {\tilde{m}}^N(s)|\Bigg ] + {\mathbb {E}}\Bigg [\sup _{s \in [0,t]} |{\tilde{m}}^N(s) - m(s)|\Bigg ]\\&\quad \le \frac{1}{N}\sum _{j=1}^N {\mathbb {E}}\Bigg [\sup _{s \in [0,t]}|\sigma _j(s) - {\tilde{\sigma }}_j(s)|\Bigg ] + {\mathbb {E}}\Bigg [\sup _{s \in [0,t]} |{\tilde{m}}^N(s) - m(s)|\Bigg ]\\&\quad = {\mathbb {E}}\Bigg [\sup _{s \in [0,t]}|\sigma _i(s) - {\tilde{\sigma }}_i(s)|\Bigg ] + C(N), \end{aligned}$$

with \(C(N) \xrightarrow {N \rightarrow +\infty } 0\) because of the chaoticity of the i.i.d. processes \(({\tilde{\sigma }}_i(t), {\tilde{y}}_i(t))_{i=1,\ldots ,N}\), and where in the equality we have used the exchangeability of the processes \((\sigma _i, {\tilde{\sigma }}_i)_{i=1,\ldots ,N}\). Collecting the estimates, we have shown, for any \(t \in [0,T]\),

$$\begin{aligned}&\frac{1}{N}\sum _{i=1}^N\left\{ {\mathbb {E}}\Bigg [\sup _{s \in [0,t]} |\sigma _i(s) - {\tilde{\sigma }}_i(s)|\Bigg ]+ {\mathbb {E}}\Bigg [\sup _{s\in [0,t]}| y_i(s) -{\tilde{y}}_i(s)|\Bigg ]\right\} \\&\quad \le C(N) + C\int _0^t \frac{1}{N}\sum _{i=1}^N{\mathbb {E}}\Bigg [\sup _{r \in [0,s]} |\sigma _i(r) - {\tilde{\sigma }}_i(r)|+ \sup _{r \in [0,s]}|y_i(r) - {\tilde{y}}_i(r)| \Bigg ]ds, \end{aligned}$$

which by the Gronwall’s lemma applied to

$$\begin{aligned} \phi (t) := \frac{1}{N}\sum _{i=1}^N\left\{ {\mathbb {E}}\Bigg [\sup _{s \in [0,t]} |\sigma _i(s) - {\tilde{\sigma }}_i(s)|\Bigg ]+ {\mathbb {E}}\Bigg [\sup _{s\in [0,t]}| y_i(s) -{\tilde{y}}_i(s)|\Bigg ]\right\} , \end{aligned}$$

implies (38), because \(\phi (T)\) is an upper bound for the left hand side of (38). \(\square \)

Remark 2

Proposition 4 and the proof of Theorem 1 can be generalized to any \(\gamma \in {\mathbb {N}}\). Indeed, analogous Lipschitz \(L^1\) estimates on the rates of Proposition 4 (used also in Theorem 1) hold by estimating

$$\begin{aligned}&\big |\lambda (\sigma ,m,y) - \lambda ({\tilde{\sigma }},{\tilde{m}},{\tilde{y}})\big |= \big | y^\gamma e^{-(\gamma + 1)\beta m \sigma } - {\tilde{y}}^\gamma e^{-(\gamma + 1)\beta {\tilde{m}}{\tilde{\sigma }}}\big |\\&\quad \le \big |y^\gamma e^{-(\gamma + 1)\beta m \sigma } - {\tilde{y}}^\gamma e^{-(\gamma + 1)\beta m\sigma } \big | + \big | {\tilde{y}}^\gamma e^{-(\gamma + 1)\beta m \sigma } - {\tilde{y}}^\gamma e^{-(\gamma + 1)\beta {\tilde{m}}{\tilde{\sigma }}}\big | \\&\quad \le \big |e^{-(\gamma + 1)\beta m \sigma }\big | \big |y^\gamma - {\tilde{y}}^\gamma \big | + {\tilde{y}}^\gamma \big |e^{-(\gamma + 1)\beta m \sigma } - e^{-(\gamma + 1)\beta {\tilde{m}}{\tilde{\sigma }}}\big | \\&\quad \le C\big |y - {\tilde{y}}\big | \big |p(y,{\tilde{y}})\big | + {\tilde{y}}^\gamma \Big [C |m - {\tilde{m}}| + C|\sigma - {\tilde{\sigma }}|\Big ]\le C\Big [|y - {\tilde{y}}| + |m - {\tilde{m}}| + |\sigma - {\tilde{\sigma }}|\Big ], \end{aligned}$$

with \(p(y,{\tilde{y}})\) a polynomial of degree \(\gamma - 1\). In the last step we have used the a priori bounds on \(y \le T\) to get \(|p(y,{\tilde{y}})| \le C(T)\) and the Lipschitz properties of \(e^{-(\gamma +1)\beta m \sigma }\) for \((\sigma ,m) \in \left\{ -1,1\right\} \times [-1,1]\).

4.3 Proofs of the Local Analysis of Sect. 2.4

In this section we address the proofs of the results illustrated in Sect. 2.4.

Proof of Proposition 1

Setting \(m= 0\) in Sys. (18), the stationary version of the first equation becomes

$$\begin{aligned} \frac{\partial }{\partial y} f(\sigma ,y) + y^\gamma f(\sigma ,y) = 0, \end{aligned}$$
(39)

whose solution is of the form \(f^*(\sigma ,y) = c(\sigma ) f(\sigma ,0) e^{-\frac{y^{\gamma +1}}{\gamma +1}}\). Denoting \(\varLambda := \int _{0}^{+\infty } e^{-\frac{y^{\gamma +1}}{\gamma +1}} dy\), it is easy to see that the integral conditions imply \(c(\sigma ) = c(-\sigma ) = \frac{1}{\varLambda }\) and \(f(\sigma ,0) = f(-\sigma ,0) = \frac{1}{2}\). \(\square \)

4.3.1 Formal Derivation of Sys. (20)

We now compute formally the linearization of the operator associated to Sys. (18) around the solution (19) with \(m = 0\). Namely, if we write the first equation in (18) in operator form

$$\begin{aligned} \frac{\partial }{\partial t} f(t,\sigma ,y) - {\mathcal {L}}_{\gamma }^{nl}f(t,\sigma ,y) = 0, \end{aligned}$$

with \({\mathcal {L}}_{\gamma }^{nl}f(t,\sigma ,y) := -\frac{\partial }{\partial y} f(t,\sigma ,y) - y^\gamma e^{-(\gamma + 1)\beta \sigma m(t)} f(t,\sigma ,y)\), we want to find the linearized version of the operator \({\mathcal {L}}_{\gamma }^{nl}\).

For the purpose, we express a generic stationary solution to (18) as

$$\begin{aligned} f(\sigma ,y) = f^*(\sigma ,y) + \varepsilon g(\sigma ,y), \end{aligned}$$

where \(f^*\) is the stationary solution corresponding to \(m=0\) and g satisfies the condition

$$\begin{aligned} \int _{0}^{\infty } [g(1,y) + g(-1,y)]dy = 0, \end{aligned}$$
(40)

so that \(\int _{0}^{\infty } [f(1,y) + f(-1,y)]dy = 1\) holds. We also denote \(m_f := \int _{0}^{\infty } [f(1,y) - f(-1,y)]dy\), which by the above consideration satisfies

$$\begin{aligned} m_f = 2\varepsilon \int _0^{\infty } g(1,y) dy =: \varepsilon k. \end{aligned}$$
(41)

The stationary version of the first equation in (18) becomes

$$\begin{aligned} \frac{\partial }{\partial y} f^*(\sigma ,y) + \varepsilon \frac{\partial }{\partial y} g(\sigma ,y) + y^\gamma e^{-\beta \sigma \varepsilon k (\gamma + 1)}[f^*(\sigma ,y) + \varepsilon g(\sigma ,y)] = 0. \end{aligned}$$

By expanding at the first order in \(\varepsilon \) the term \(e^{-\beta \sigma \varepsilon k (\gamma + 1)} \approx 1 - (\gamma + 1)\beta \sigma \varepsilon k\), and by considering only the resulting linear terms in \(\varepsilon \), we get

$$\begin{aligned} \frac{\partial }{\partial y} f^*(\sigma ,y) + \varepsilon \frac{\partial }{\partial y} g(\sigma ,y) + y^\gamma f^*(\sigma ,y) + y^\gamma \varepsilon g(\sigma ,y) - y^\gamma (\gamma +1)\beta \sigma \varepsilon k f^*(\sigma ,y) = 0. \end{aligned}$$

Finally, using that \(f^*\) solves (39) and substituting its expression (19), we get

$$\begin{aligned} \frac{\partial }{\partial y} g(\sigma ,y) + y^\gamma g(\sigma ,y) - \frac{\beta \sigma k(\gamma +1)}{2\varLambda } y^\gamma e^{-\frac{y^{\gamma + 1}}{\gamma + 1}} = 0. \end{aligned}$$

We can define the linearized operator as

$$\begin{aligned} {\mathcal {L}}_{\gamma }^{\text {lin}}g(\sigma ,y) := -\frac{\partial }{\partial y} g(\sigma ,y) - y^\gamma g(\sigma ,y) + \frac{\beta \sigma k(\gamma +1)}{2\varLambda } y^\gamma e^{-\frac{y^{\gamma + 1}}{\gamma + 1}}. \end{aligned}$$
(42)

We proceed with the linearization of the integral condition in the second line of Sys. (18):

$$\begin{aligned}&f^*(\sigma ,0) + \varepsilon g(\sigma ,0) = \int _0^\infty [f^*(-\sigma ,y) + \varepsilon g(-\sigma ,y)] y^\gamma e^{\beta \sigma \varepsilon k (\gamma + 1)}dy\\&\quad \approx \int _0^\infty f^*(-\sigma ,y) y^\gamma ( 1+ \beta \sigma \varepsilon k (\gamma + 1))dy + \varepsilon \int _0^\infty g(-\sigma ,y) y^\gamma ( 1+ \beta \sigma \varepsilon k (\gamma + 1))dy\\&\quad \approx \int _0^\infty f^*(-\sigma ,y) y^\gamma dy+ \beta \sigma \varepsilon k (\gamma + 1)\int _0^\infty f^*(-\sigma ,y) y^\gamma dy+ \varepsilon \int _0^\infty g(-\sigma ,y) y^\gamma dy. \end{aligned}$$

Using again that \(f^*\) solves (39) and its expression in (19), we get

$$\begin{aligned} g(\sigma ,0) = \frac{\beta \sigma k (\gamma + 1)}{2\varLambda } + \int _0^\infty g(-\sigma ,y) y^\gamma dy. \end{aligned}$$
(43)

In order to gain indications on the stability properties of the stationary solution to (18) with \(m = 0\), we study the discrete spectrum of \({\mathcal {L}}_{\gamma }^{\text {lin}}\) defined in (42), i.e., we search for the eigenfunctions g and the eigenvalues \(\lambda \in {\mathbb {C}}\), satisfying the linearized integral conditions (40) and (43) found above, and such that

$$\begin{aligned} {\mathcal {L}}_{\gamma }^{\text {lin}} g(\sigma ,y) = \lambda g(\sigma ,y), \end{aligned}$$
(44)

which is equivalent to

$$\begin{aligned} \frac{\partial }{\partial y} g(\sigma ,y) + y^\gamma g(\sigma ,y) - \frac{\beta \sigma k(\gamma + 1)}{2\varLambda }y^\gamma e^{-\frac{y^{\gamma +1}}{\gamma +1}} = -\lambda g(\sigma ,y). \end{aligned}$$
(45)

The eigen-system around \(m = 0\) is thus given by (20), where, recall by (41), \(k = 2\int _0^\infty g(1,y) dy\), and \(\varLambda = \int _0^\infty e^{-\frac{y^{\gamma +1}}{\gamma + 1}} dy\).

Remark 3

The derivation of the linearized operator (44) was formal. One could think to define it more rigorously, by indicating an Hilbert space where \({\mathcal {L}}_{\gamma }^{\text {lin}}\) acts on. The natural choice appears to be (a subspace of) \(\Big (L^2_{\mu _{\gamma }}\big ({\mathbb {R}}^+\big )\Big )^2\) satisfying conditions (40) and (43), where the outer square comes from the explicitation of the spin variable \(\sigma = \pm 1\), and the measure \(\mu _\gamma \) is defined as

$$\begin{aligned} \mu _\gamma (dy) := f^*(\sigma ,y)dy = \frac{1}{2\varLambda }e^{-\frac{y^{\gamma + 1}}{\gamma +1}} dy. \end{aligned}$$
(46)

As in the computations we do not use the particular choice of domain of the operator or its properties, we do not investigate further on this.

Proof of Proposition 2

Recall that in this case we have \(\gamma =1\). In order to solve the first equation in (21), we set \(h(\sigma ,y) := g(\sigma ,y) e^{\frac{y^2}{2}}\). It holds

$$\begin{aligned} \frac{\partial }{\partial y}h(\sigma ,y) = -\lambda h(\sigma ,y) + \frac{y \beta \sigma k}{\sqrt{\frac{\pi }{2}}}, \end{aligned}$$

whose solution is

$$\begin{aligned} h(\sigma ,y) = e^{-\lambda y}\left[ h(\sigma ,0) + \frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}}\int _{0}^y u e^{\lambda u} du\right] . \end{aligned}$$

Noting that \(\int _0^y u e^{\lambda u} du = \frac{1}{\lambda ^2} - \frac{e^{\lambda y}}{\lambda ^2} + \frac{e^{\lambda y}}{\lambda }y\), we obtain

$$\begin{aligned} g(\sigma ,y) = e^{-\frac{y^2}{2}} e^{-\lambda y} \left[ g(\sigma ,0) + \frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}}\left( \frac{1}{\lambda ^2} - \frac{e^{\lambda y}}{\lambda ^2} + \frac{e^{\lambda y}}{\lambda }y\right) \right] . \end{aligned}$$
(47)

We now impose the integral conditions. First, we note that \(\int _0^\infty [g(\sigma ,y) + g(-\sigma ,y)]dy = 0\) is equivalent to \(g(\sigma ,y) + g(-\sigma ,y) = 0\) for every \(y \in {\mathbb {R}}^+\) because of expression (47). For the computation of k, recalling notation (23), we find

$$\begin{aligned} k = 2\int _0^\infty g(1,y) dy = 2g(1,0) H_1(\lambda ) + 2\frac{\beta k}{\sqrt{\frac{\pi }{2}}}\frac{1}{\lambda ^2} H_1(\lambda ) - 2\frac{\beta k}{\lambda ^2}\frac{1}{\sqrt{\frac{\pi }{2}}}\sqrt{\frac{\pi }{2}} + 2\frac{\beta k}{\sqrt{\frac{\pi }{2}}}\frac{1}{\lambda }, \end{aligned}$$

so that

$$\begin{aligned} k = \frac{2 g(1,0) H_1(\lambda )}{1 - 2\frac{\beta }{\lambda \sqrt{\frac{\pi }{2}}} - 2\frac{\beta H_1(\lambda )}{\lambda ^2 \sqrt{\frac{\pi }{2}}} + 2\frac{\beta }{\lambda ^2}}. \end{aligned}$$
(48)

The integral condition in the second line of (21) gives

$$\begin{aligned} g(\sigma ,0)&= \frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}} + \int _0^\infty y\left[ e^{-\frac{y^2}{2}}e^{-\lambda y}\left( g(-\sigma ,0) - \frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}} \left( \frac{1}{\lambda ^2} - \frac{e^{\lambda y}}{\lambda ^2} + \frac{e^{\lambda y}}{\lambda } y\right) \right) \right] dy\\&= \frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}} - g(\sigma ,0)(1-\lambda H_1(\lambda )) - \frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}} \frac{(1-\lambda H_1(\lambda ))}{\lambda ^2} \\&\quad + \frac{1}{\lambda ^2} \frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}} - \frac{1}{\lambda } \frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}}\int _0^\infty y^2 e^{-\frac{y^2}{2}}dy \\&= \frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}} - g(\sigma ,0)(1-\lambda H_1(\lambda )) - \frac{(1-\lambda H_1(\lambda ))}{\lambda ^2}\frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}} + \frac{1}{\lambda ^2}\frac{\beta \sigma k}{\sqrt{\frac{\pi }{2}}} - \frac{1}{\lambda }\beta \sigma k. \end{aligned}$$

In the second equality we have used that \(\int _0^\infty y e^{-\frac{y^2}{2}}e^{-\lambda y}dy = 1-\lambda H_1(\lambda )\) which can be obtained by an integration by parts using the property \(g(\sigma ,0)=-g(-\sigma ,0)\). Solving for g(1, 0) in the above

$$\begin{aligned} g(1,0)[2 - \lambda H_1(\lambda )] = \beta k \left[ \frac{1}{\sqrt{\frac{\pi }{2}}} - \frac{(1-\lambda H_1(\lambda ))}{\lambda ^2}\frac{1}{\sqrt{\frac{\pi }{2}}} + \frac{1}{\lambda ^2}\frac{1}{\sqrt{\frac{\pi }{2}}} - \frac{1}{\lambda }\right] . \end{aligned}$$

Substituting the value of k we found in (48), we get

$$\begin{aligned} g(1,0)[2&- \lambda H_1(\lambda )] = \frac{2 \beta g(1,0) H_1(\lambda )}{1 - 2\frac{\beta }{\lambda \sqrt{\frac{\pi }{2}}} - 2\frac{\beta H_1(\lambda )}{\lambda ^2 \sqrt{\frac{\pi }{2}}} + 2\frac{\beta }{\lambda ^2}}\left[ \frac{1}{\sqrt{\frac{\pi }{2}}} - \frac{(1-\lambda H_1(\lambda ))}{\lambda ^2}\frac{1}{\sqrt{\frac{\pi }{2}}} + \frac{1}{\lambda ^2}\frac{1}{\sqrt{\frac{\pi }{2}}} - \frac{1}{\lambda }\right] , \end{aligned}$$

which is equivalent to

$$\begin{aligned} 2 - \lambda H_1(\lambda ) = \frac{2\beta H_1(\lambda )\left[ \lambda ^2 + \lambda H_1(\lambda ) - \lambda \sqrt{\frac{\pi }{2}}\right] }{\lambda ^2\sqrt{\frac{\pi }{2}} - 2\beta \lambda - 2\beta H_1(\lambda ) + 2\beta \sqrt{\frac{\pi }{2}}}. \end{aligned}$$
(49)

As a polynomial in \(\lambda \), (49) can be written as

$$\begin{aligned} -\lambda ^3 H_1(\lambda )\sqrt{\frac{\pi }{2}} + \lambda ^2\sqrt{2\pi } -4\beta \lambda - 4\beta H_1(\lambda ) + 2\sqrt{2\pi }\beta = 0, \end{aligned}$$

or, grouping for \(H_1(\lambda )\),

$$\begin{aligned} H_1(\lambda )\left[ -4\beta - \lambda ^3 \sqrt{\frac{\pi }{2}}\right] + \sqrt{2\pi }\lambda ^2 - 4\beta \lambda + 2\beta \sqrt{2\pi } = 0, \end{aligned}$$

i.e. the zeros of \(H_{\beta ,1}(\lambda )\). In fact, as defined in (23), \(H_1(\lambda )\) is a holomorphic function on \({\mathbb {C}}\), whose expression in series is

$$\begin{aligned} H_1(\lambda )&= \int _0^\infty e^{-\frac{y^2}{2}}e^{-\lambda y} dy = \sum _{n=0}^{\infty } (-1)^n \frac{\lambda ^n}{n!}\int _0^\infty y^n e^{-\frac{y^2}{2}}dy. \end{aligned}$$

The latter integral is known

$$\begin{aligned} \int _0^\infty y^n e^{-\frac{y^2}{2}}dy = 2^{\frac{1}{2}(n-1)}\varGamma \left( \frac{n+1}{2}\right) , \end{aligned}$$
(50)

where \(n\in {\mathbb {N}}\) and \(\varGamma (\cdot )\) is the Gamma function. When \(n = 2m + 1\) with \(m\in {\mathbb {N}}\) , for the properties of the Gamma function on \({\mathbb {N}}\), (50) reduces to

$$\begin{aligned} \int _0^\infty y^n e^{-\frac{y^2}{2}}dy = 2^{\frac{1}{2}(n-1)}\varGamma \left( \frac{n+1}{2}\right) = 2^m m!. \end{aligned}$$

For \(n = 2m\) with \(m\in {\mathbb {N}}\) instead we have, by the property \(\varGamma \left( l +\frac{1}{2}\right) = \frac{(2l -1)!!}{2^l}\sqrt{\pi }\) for any \(l \in {\mathbb {N}}\),

$$\begin{aligned} \int _0^\infty y^n e^{-\frac{y^2}{2}}dy = 2^{\frac{1}{2}(n-1)}\varGamma \left( \frac{n+1}{2}\right) = \sqrt{\frac{\pi }{2}}(2m-1)!!. \end{aligned}$$

We use these equalities, and reorder the terms of the absolutely convergent series of \(H_1(\lambda )\) to finally get

$$\begin{aligned} H_1(\lambda ) = \sqrt{\frac{\pi }{2}}\sum _{m=0}^{\infty }\frac{\lambda ^{2m}}{(2m)!!} - \lambda \sum _{m=0}^{\infty } \frac{(2\lambda )^{2m} m!}{(2m+1)!}\frac{1}{2^m}. \end{aligned}$$

\(\square \)

Proof of Proposition 3

Now we set \(\gamma =2\) and proceed as for \(\gamma =1\), by setting \(h(\sigma ,y) := g(\sigma ,y) e^{\frac{y^3}{3}}\), so that

$$\begin{aligned} \frac{\partial }{\partial y} h(\sigma ,y) = -\lambda h(\sigma ,y) + \frac{3}{2\varLambda }\beta \sigma k y^2. \end{aligned}$$

Thus,

$$\begin{aligned} h(\sigma ,y) = e^{-\lambda y}\left[ h(\sigma ,0) + \frac{3\beta \sigma k}{2\varLambda }\int _0^y u^2 e^{\lambda u} du\right] . \end{aligned}$$

Since \(\int _0^y u^2 e^{\lambda u} du = \frac{1}{\lambda ^3}[-2 + e^{\lambda y}(2 + \lambda y(-2 + \lambda y))]\), we can write

$$\begin{aligned} \begin{aligned} g(\sigma ,y) = e^{-\frac{y^3}{3}}e^{-\lambda y} \left[ g(\sigma ,0) + \frac{3\beta \sigma k}{2\varLambda }\frac{1}{\lambda ^3}\left( - 2 + 2e^{\lambda y} -2\lambda y e^{\lambda y} + \lambda ^2 y^2 e^{\lambda y}\right) \right] . \end{aligned} \end{aligned}$$
(51)

Recalling notation (28), we compute

$$\begin{aligned} k&= 2\int _0^\infty g(1,y) dy = 2 H_2(\lambda ) g(1,0) - 2H_2(\lambda ) \frac{3 \beta k}{\varLambda \lambda ^3} + 2\frac{3 \beta k}{\varLambda \lambda ^3}\int _0^\infty e^{-\frac{y^3}{3}}dy \\&\quad - 2\frac{3\beta k}{\varLambda \lambda ^2}\int _0^\infty y e^{-\frac{y^3}{3}}dy + \frac{3 \beta k}{\varLambda \lambda }\int _0^\infty y^2 e^{-\frac{y^3}{3}}dy\\&= 2 H_2(\lambda ) g(1,0) - 2H_2(\lambda ) \frac{3 \beta k}{\varLambda \lambda ^3} + 2 \frac{3 \beta k}{\lambda ^3} - 2 \frac{3\beta k}{\varLambda \lambda ^2} \frac{\varGamma (2/3)}{3^{1/3}} + \frac{3\beta k}{\varLambda \lambda }, \end{aligned}$$

which gives

$$\begin{aligned} k = \frac{2 g(1,0) H_2(\lambda )}{1 + 2H_2(\lambda )\frac{3\beta }{\varLambda \lambda ^3} - 2\frac{3\beta }{\lambda ^3} + 2\frac{3\beta }{\varLambda \lambda ^2}\frac{\varGamma (2/3)}{3^{1/3}} - \frac{3\beta }{\varLambda \lambda }}. \end{aligned}$$
(52)

As before, the condition \(\int _0^\infty [g(\sigma ,y) + g(-\sigma ,y)]dy = 0\) in (26) is equivalent to \(g(\sigma , y) + g(-\sigma , y) = 0\) for every \(y \in {\mathbb {R}}^+\) because of (51). Using this observation for \(y = 0\) in the other integral condition - see second line of system (26) - we compute

$$\begin{aligned} g(\sigma ,0)&= \frac{3 \beta \sigma k}{2\varLambda } + \int _0^\infty y^2 \Bigg [e^{-\frac{y^3}{3}}e^{-\lambda y}\Bigg (-g(\sigma ,0)\\&\quad - \frac{3\beta \sigma k}{2\varLambda } \frac{1}{\lambda ^3}(-2 + 2 e^{\lambda y} - 2\lambda y e^{\lambda y} + \lambda ^2 y^2 e^{\lambda y})\Bigg )\Bigg ]dy. \end{aligned}$$

Observing that, by integration by parts, \(\int _0^\infty y^2 e^{-\frac{y^3}{3}} e^{-\lambda y} dy = 1 -\lambda H_2(\lambda )\), we find

$$\begin{aligned} g(\sigma ,0)&= \frac{3 \beta \sigma k}{2\varLambda } - (1 -\lambda H_2(\lambda ))g(\sigma ,0) + \frac{3\beta \sigma k}{\varLambda \lambda ^3}(1-\lambda H_2(\lambda )) - \frac{3\beta \sigma k}{\varLambda \lambda ^3} \\&\quad + \frac{3\beta \sigma k }{\varLambda \lambda ^2}3^{1/3}\varGamma (4/3) - \frac{3\beta \sigma k}{2\varLambda \lambda }3^{2/3}\varGamma (5/3). \end{aligned}$$

Computing in \(\sigma = 1\) and grouping for g(1, 0),

$$\begin{aligned} g(1,0)[2 - \lambda H_2(\lambda )]&= k\left[ \frac{3 \beta }{2\varLambda } + \frac{3\beta }{\varLambda \lambda ^3}(1-\lambda H_2(\lambda )) - \frac{3\beta }{\varLambda \lambda ^3} + \frac{3\beta }{\varLambda \lambda ^2}3^{1/3}\varGamma (4/3) - \frac{3\beta }{2\varLambda \lambda }3^{2/3}\varGamma (5/3)\right] \\&= k \left[ \frac{3 \beta }{2\varLambda } - \frac{3\beta }{\varLambda \lambda ^2}H_2(\lambda ) + \frac{3\beta }{\varLambda \lambda ^2}3^{1/3}\varGamma (4/3) - \frac{3\beta }{2\varLambda \lambda }3^{2/3}\varGamma (5/3)\right] . \end{aligned}$$

Plugging expression (52) for k,

$$\begin{aligned} 2 - \lambda H_2(\lambda )&= \frac{2 H_2(\lambda )}{1 + 2H_2(\lambda )\frac{3\beta }{\varLambda \lambda ^3} - 2\frac{3\beta }{\lambda ^3} + 2\frac{3\beta }{\varLambda \lambda ^2}\frac{\varGamma (2/3)}{3^{1/3}} - \frac{3\beta }{\varLambda \lambda }}\left[ \frac{3}{2\varLambda }\beta - \frac{3\beta }{\varLambda \lambda ^2}H_2(\lambda ) \right. \\&\left. \quad + \frac{3\beta }{\varLambda \lambda ^2}3^{1/3}\varGamma (4/3) - \frac{3\beta }{2\varLambda \lambda }3^{2/3}\varGamma (5/3)\right] . \end{aligned}$$

This gives

$$\begin{aligned} 2 - \lambda H_2(\lambda ) = \frac{2\lambda H_2(\lambda )\left[ \frac{3}{2}\beta \lambda ^2 - 3\beta H_2(\lambda ) + 3\beta 3^{1/3}\varGamma (4/3) - \frac{3}{2}\beta \lambda 3^{2/3}\varGamma (5/3)\right] }{\varLambda \lambda ^3 + 6\beta H_2(\lambda ) - 6\beta \varLambda + 6\beta \lambda \frac{\varGamma (2/3)}{3^{1/3}} - 3\beta \lambda ^2}, \end{aligned}$$

which is equivalent to

$$\begin{aligned}&2 \varLambda \lambda ^3 + 12 \beta H_2(\lambda ) - 12\beta \varLambda + 12\beta \lambda \frac{\varGamma (2/3)}{3^{1/3}} \\&\qquad - 6\beta \lambda ^2 - \lambda ^4 \varLambda H_2(\lambda ) + 6\beta \lambda H_2(\lambda )\varLambda - 6\beta \lambda ^2\frac{\varGamma (2/3)}{3^{1/3}} H_2(\lambda ) \\&\qquad = 6\beta \lambda H_2(\lambda )3^{1/3}\varGamma (4/3) - 3\beta \lambda ^2 H_2(\lambda ) 3^{2/3}\varGamma (5/3). \end{aligned}$$

As a polynomial in \(\lambda \), this is

$$\begin{aligned}&-\varLambda H_2(\lambda )\lambda ^4+ 2\varLambda \lambda ^3 \\&\quad + \lambda ^2\left[ -6\beta -6\beta \frac{\varGamma (2/3)}{3^{1/3}}H_2(\lambda ) + 3\beta H_2(\lambda ) 3^{2/3}\varGamma (5/3)\right] \\&\quad +\lambda \left[ 6\beta \varLambda H_2(\lambda ) -6\beta H_2(\lambda )3^{1/3}\varGamma (4/3) + 12\beta \frac{\varGamma (2/3)}{3^{1/3}}\right] + 12\beta H_2(\lambda ) - 12\beta \varLambda = 0. \end{aligned}$$

Equivalently, in terms of \(H_2(\lambda )\) we have

$$\begin{aligned}&H_2(\lambda )[12\beta - \lambda ^4 \varLambda + 6\beta \lambda \varLambda - 6\beta \lambda 3^{1/3}\varGamma (4/3)\\&\quad + 3\beta \lambda ^2 3^{2/3}\varGamma (5/3)- 6\beta \lambda ^2\frac{\varGamma (2/3)}{3^{1/3}}]\\&\quad + \left[ 2\varLambda \lambda ^3 - 12\beta \varLambda + 12\beta \frac{\varGamma (2/3)}{3^{1/3}}\lambda - 6\beta \lambda ^2 \right] = 0, \end{aligned}$$

i.e. the zeros of \(H_{\beta ,2}(\lambda )\) in (27). As defined in (28), \(H_2(\lambda )\) is a holomorphic function on \({\mathbb {C}}\), which can be expressed in series as

$$\begin{aligned} H_2(\lambda )&= \int _0^\infty e^{-\lambda y}e^{-\frac{y^3}{3}}dy = \sum _{n=0}^{\infty } (-1)^n \frac{\lambda ^n}{n!}\int _0^\infty y^n e^{-\frac{y^3}{3}}dy\\&= \sum _{n=0}^{\infty } (-1)^n \frac{\lambda ^n}{n!} 3^{\frac{1}{3}(n-2)}\varGamma \left( \frac{n+1}{3}\right) , \end{aligned}$$

which is expression (29), where we have used the formula for \(\int _0^\infty y^n e^{-\frac{y^3}{3}}dy = 3^{\frac{1}{3}(n-2)}\varGamma \left( \frac{n+1}{3}\right) \). \(\square \)