Oscillatory behavior in a model of non-Markovian mean field interacting spins

We analyze a non-Markovian mean field interacting spin system, related to the Curie--Weiss model. We relax the Markovianity assumption by replacing the memoryless distribution of the waiting times of a classical spin-flip dynamics with a distribution with memory. The resulting stochastic evolution for a single particle is a spin-valued renewal process, an example of two-state semi-Markov process. We associate to the individual dynamics an equivalent Markovian description, which is the subject of our analysis. We study a corresponding interacting particle system, where a mean field interaction is introduced as a time scaling, depending on the overall magnetization of the system, on the waiting times between two successive particle's jumps. Via linearization arguments on the Fokker-Planck mean field limit equation, we give evidence of emerging periodic behavior. Specifically, numerical analysis on the discrete spectrum of the linearized operator, characterized by the zeros of an explicit holomorphic function, suggests the presence of a Hopf bifurcation for a critical value of the temperature, which is in accordance with the one obtained by simulating the N-particle system.


Introduction
Emerging periodic behavior in complex systems with a large number of interacting units is a commonly observed phenomenon in neuroscience ( [14]), ecology ( [23]), socioeconomics ( [4,24]) 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 where each individual particle has 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.Throughout the paper we refer to the emergence of self-organized periodic oscillations 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 self-organization.Among others, we cite noise ( [10], [20], [22]), dissipation in the interaction potential ( [1], [6], [7], [9]), delay in the transmission of information and/or frustration in the interaction network ( [8], [12], [21]).In particular, in [12] 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 timeirreversible phenomenon in the dynamics, as it is the case in all the above cited works (see e.g.[2], [16]).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 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 neutral phase for values of the parameters corresponding to high temperatures, the emergence of periodic orbits 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 γ+1 , for γ = 1, 2. Then, we introduce an interaction among the spins via a time rescaling depending on the overall magnetization of the system.The specific choice of interarrival time distribution makes the computations developed in Section 4 the easiest as possible (to our knowledge), allowing for an explicit characterization of the discrete spectrum of the linearized operator.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 γ+1 with γ ∈ R such that γ ≥ 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 Section 2 we describe the model and the results obtained.In particular, before introducing the model (Subsection 2.2), we start by recalling basic facts about the Curie-Weiss model and its phase transitions (Subsection 2.1); we then proceed with the results on the propagation of chaos (Subsection 2.3), and on the linearized Fokker-Planck equation around a neutral equilibrium, for two different choices of renewal dynamics (Subsection 2.4).Notably, we determine the discrete spectrum of the linearized operator in terms of the zeros of two holomorphic functions.Section 3 contains the numerical results on the discrete spectrum, studied as a function of the interaction parameters (Subsection 3.1).These results are then compared in Subsection 3.2 with the ones obtained by simulating the finite particle system, finding a precise accordance between the two approaches.Section 4 contains the proofs of the results of Section 2.

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 σ i ∈ {−1, 1}, i = 1, . . ., N , which is reversible with respect to the equilibrium Gibbs probability measure on the space of configurations with σ := (σ 1 , . . ., σ N ) ∈ {−1, 1} N , β > 0 (ferromagnetic case), Z N (β) is a normalizing constant, and H is the Hamiltonian, which in the Curie-Weiss setting is given by Denote also the empirical magnetization as m N := 1 N N i=1 σ 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 β > 0, which can be recognized by proving a Law of Large Numbers for the equilibrium empirical magnetization where m β > 0 is the so-called spontaneous magnetization.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 : where N is obtained from σ 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 L applied to a function g : [−1, 1] → R: It is easy to obtain the weak limit of the sequence of processes m N (t) t≥0 , by studying the uniform convergence of the generator (5) as N → +∞ (see e.g.[15]).The limit process (m(t)) t≥0 is deterministic and solves the Curie-Weiss ODE The presence of the phase transition highlighted in (3) can be recognized as well in the outof-equilibrium dynamical model (6).Indeed, studying the long-term behavior of (6), one finds that: • for β ≤ 1, (6) possesses a unique stationary solution, globally attractive, constantly equal to 0; • for β > 1, 0 is still stationary but it is unstable; two other symmetric stationary locally attractive solutions, ±m β , appear: the two non-zero solutions to m = tanh(βm).The dynamics m(t) gets attracted for t → +∞ 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 a 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) {T n } ∞ n=1 , 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 P(T i > s + t|T i > t) = P(T i > s), holds for any s, t ≥ 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 continuoustime 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 Levy in [19].

The dynamics.
In order to introduce the model, we start by observing that the Curie-Weiss dynamics (4), as any spin-flip Glauber dynamics, can be obtained by adding interaction to a 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.Our aim here is to replace Poisson processes by more general renewal processes, otherwise keeping the structure of the interaction.
For the moment we focus on a single spin σ(t) ∈ {−1, 1}.If driven by a Poisson process of intensity 1, its dynamics has infinitesimal generator (7) Lf 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 two-states 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 τ (interchangeably referred to as interarrival times) of the renewal satisfy ( 8) for some smooth function ϕ : [0, +∞) → R.Then, the pair (σ(t), y(t)) t≥0 is Markovian with generator ( 9) for This is equivalent to say that the couple (σ(t), y(t)) t≥0 evolves according to Expression (10) for the jump rate follows by observing that, for an interarrival time τ of the jump process σ(t), we have for any h > 0. Observe that when the τ 's are exponentially distributed F (y) ≡ 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 σ; the simplest way is to model this dependence as a time scaling: Under this distribution for the waiting times the generator of (σ(t), y(t)) t≥0 becomes: 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 (σ i (t), y i (t)) i=1,...,N , we set m N (t) := 1 N N i=1 σ i (t) to be the magnetization of the system at time t, and a parameter β > 0 tuning the interaction between the particles.The interacting dynamics is where σ i is obtained from σ by flipping the i-th spin, while y i by setting to zero the i-th coordinate.
The additional factor e −βσ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) = − ϕ (y) ϕ(y) .Note that, for F ≡ 1, we retrieve the Curie-Weiss dynamics (4) for the spins.

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 γ , for γ ∈ N, which corresponds to considering, in the single spin model, the tails of the distribution of the interarrival times to be ϕ(t) ∝ e − t γ+1 γ+1 .When F (y) = y γ , ( 13) becomes (15) (σ i (t), y i (t)) → (−σ i (t), 0), with rate y γ i (t)e −(γ+1)βσi(t)m N (t) , dy i (t) = dt, otherwise.
As for the Curie-Weiss model, dynamics ( 15) is subject to a cooperative-type interaction: the spinflip rate is larger for particles which are not aligned with the majority.Assuming propagation of chaos, at the macroscopic limit N → +∞ the representative particle (σ(t), y(t)) has a mean-field dynamics To this dynamics we can associate (see [18]) the non-linear infinitesimal generator ( 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 (σ(t), y(t)).In Section 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 γ ∈ N, and let T > 0 be the final time in (15) and (16).Assume that (σ i (0), y i (0)) i=1,...,N are µ 0 -chaotic for some probability distribution µ 0 on {−1, 1} × R + .Then, the sequence of empirical measures (µ N t ) t∈[0,T ] converges in distribution (in the sense of weak convergence of probability measures) to the deterministic law (µ t ) t∈[0,T ] on the path space of the unique solution to Eq. ( 16) with initial distribution µ 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 γ = 1 and γ = 2. Our approach is the following: we find a neutral stationary solution of interest, we linearize formally the dynamics around that equilibrium and 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 β,γ (λ).In Subsection 3.1 we then study numerically the character of the eigenvalues when β varies: for both γ = 1, 2, we find that for all β < β c (γ) all eigenvalues have negative real part; at β c (γ) two eigenvalues are conjugate and purely imaginary, suggesting the possible presence of a Hopf bifurcation in the limit dynamics.These critical values of β are then compared to the ones obtained by simulating the finite particle system in Subsection 3.2.
The Fokker-Planck equation associated to ( 16) is a PDE describing the time evolution of the density function f (t, σ, y) of the limit process (σ(t), y(t)).It is given by ( 18) 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 σ = 1 and another for σ = −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): , is a stationary solution to Sys. (18) with m = 0.
The linearization of the operator associated to Sys. (18) around the neutral equilibrium (19), yields the following eigen-system (20) where k = 2 Moreover, it holds In this case the eigen-system is given by (25) e −λy e − y 3 3 dy.
Moreover, it holds

Numerical results
3.1.Numerical evidence on the eigenvalues.We studied numerically the two eigenvalues equations (29) H β,1 (λ) = 0, and We used 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 β.Here we report the results: • Case γ = 1: ( (2.5) apart from being sensibly slower, the numerical root finding for γ = 2 suffers from numerical instability issues.This is why we were able to check the results for much smaller intervals in this case.
3.2.Finite particle system simulations.We made several simulations of the particle system (N large but finite, N = 1500) for γ = 1, 2, which seem in accordance with the above numerical results on the eigenvalues (compare with (31) and ( 32)).This is a description of the evidences: • For β small the system is stable, in particular the magnetization goes to zero regardless of the initial datum (Figure 1).• There is a critical β (around 0.75 for γ = 1, 0.35 for γ = 2) above which the magnetization starts oscillating.Close to the critical points oscillations (Figure 2) do not look very regular (corrupted by noise?), but they soon become very regular if β is not too close to the critical value.We also made joint plots of the magnetization with the empirical mean of the y i 's (Figure 3).A limit cycle seems to emerge. .Simulation of the finite particle system's dynamics for γ = 1 (left) and γ = 2 (right), with number of spins N = 1500.We plot the empirical magnetization of the spins against the empirical mean of the y i 's.
• As β increases, the amplitude of the oscillation of the magnetization increases (Figure 4), while the period looks nearly constant.As β crosses another critical value (around 1.3 for γ = 1, 1.65 for γ = 2) oscillations disappear, and the sistem magnetizes, i.e. the magnetization stabilizes to a non-zero value, actually close to ±1 (Figure 5).• The oscillations are lasting for a wider interval of β's for γ = 2 (from β ≈ 0.35 until β ≈ 1.65) than γ = 1 (from β ≈ 0.75 until β ≈ 1.3).The period is instead smaller for γ = 2 than for γ = 1.• For both γ = 1, 2, the appearance of the oscillations does not seem to depend on the initial data for the dynamics, suggesting the possible presence of a global Hopf bifurcation.Here we prove rigorously a propagation of chaos property for the N -particle dynamics to its mean-field limit, for any γ ∈ N. Actually, we establish the proofs for γ = 1, where the rates enjoy globally Lipschitz properties, and then we generalize them to any γ ∈ 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 ≤ y i ≤ T , where T < ∞ is the final time horizon of the dynamics.For the convenience of the reader, we write again the dynamics (33) (σ i (t), y i (t)) → (−σ i (t), 0), with rate y γ i (t)e −(γ+1)βσi(t)m N (t) , dy i (t) = dt, otherwise, and the mean-field version (34) (σ(t), y(t)) → (−σ(t), 0), with rate y γ (t)e −(γ+1)βσ(t)m(t) , dy(t) = dt, otherwise, with m(t) = E[σ(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 [17].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 (Ω, F, P), (F t ) t∈[0,T ] satisfying the usual hypotheses, rich enough to carry an inependent and identically distributed family (N i ) i∈N of stationary Poisson random measures N i on [0, T ] × Ξ, with intensity measure ν on Ξ := [0, +∞) equal to the restriction of the Lebesgue measure onto [0, +∞).For any N , consider the system of Itô-Skorohod equations and the corresponding limit non-linear reference particle's dynamics ( 36) Proof.With the choices in (37), the well-posedeness of Eqs. ( 35) and (36) follows by Theorems 1.2 and 2.1 in [17].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 where in the last step we have used that, by construction, the processes y i (t) ≤ T for every t ∈ [0, T ], so that the rates are a priori bounded and the Lipschitz properties of ye −2βσm for (y, σ, m

Now, define the empirical measures µ
, and their evaluation along the paths of ( 35), (38) The measures (µ N t ) t∈[0,T ] can be viewed as random variables with values in P(D), the space of probability measures on D, where D := D [0, T ]; {−1, 1}×R + is the space of {−1, 1}×R + -valued càdlàg functions equipped with the Skorohod topology.

Proofs of the local analysis of Subsection 2.4.
In this section we address the proofs of the results illustrated in Subsection 2.4.We start with a remark on the derivation of the Fokker-Planck mean-field limit equation: Remark 3.While the other equations in (18) are derived in a standard way from the expression of the generator (17), the boundary integral condition might need to be motivated.In words, it is a mass-balance between the spins that have just jumped (thus having y = 0).We reason heuristically by discretizing the state space [0, +∞) in small intervals of amplitude ε.The discretized version of y(t) takes values in {nε : n ∈ N}.The associated generator is, for n ∈ N and σ ∈ {−1, 1}, Denoting f (t, σ, 0) the density of the discretized process in (σ, 0) at time t, it follows from the expression of L ε , that is the discretized version of the integral condition in (17).
Proof of Proposition 1. Setting m = 0 in Sys.(18), the stationary version of the first equation becomes e − y γ+1 γ+1 , it is easy to see that the integral conditions imply c(σ) = c(−σ) = 1 Λ and f (σ, 0) = f (−σ, 0) = 1 2 .4.2.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 ∂ ∂t f (t, σ, y) − L nl γ f (t, σ, y) = 0, with L nl γ f (t, σ, y) := − ∂ ∂y f (t, σ, y)−y γ e −(γ+1)βσm(t) f (t, σ, y), we want to find the linearized version of the operator L nl γ .For the purpose, we express a generic stationary solution to (18) as Finally, using that f * solves (41) and substituting its expression (19) We proceed with the linearization of the integral condition in the second line of Sys.(18): Using again that f * solves (41) and its expression in (19), we get In order to gain indications on the stability properties of the stationary solution to ( 18) with m = 0, we study the discrete spectrum of L lin γ defined in (44), i.e., we search for the eigenfunctions g and the eigenvalues λ ∈ C, satisfying the linearized integral conditions (42) and ( 45 We now impose the integral conditions.First, we note that ∞ 0 [g(σ, y) + g(−σ, y)]dy = 0 is equivalent to g(σ, y) + g(−σ, y) = 0 for every y ∈ R + because of expression (49).For the computation of k, recalling notation (23), we find The integral condition in the second line of (21) gives In the second equality we have used that ∞ 0 ye − y 2 2 e −λy = 1 − λH 1 (λ) which can be obtained by an integration by parts.Solving for g(1, 0) in the above Substituting the value of k we found in (50), we get As a polynomial in λ, (51) can be written as or, grouping for H 1 (λ), i.e. the zeros of H β,1 (λ), provided we prove expression (24) for H 1 (λ).In fact, as defined in ( 23), H 1 (λ) is a holomorphic function on C, whose expression in series is We use these equalities, and reorder the terms of the absolutely convergent series of H 1 (λ) to finally get Proof of Proposition 3. We proceed as in the previous case, by setting h(σ, y) := g(σ, y)e
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.