Kuramoto model and mean field spin XY dynamics

# Dynamical aspects of mean field plane rotators and the Kuramoto model

## Abstract.

The Kuramoto model has been introduced in order to describe synchronization phenomena observed in groups of cells, individuals, circuits, etc… We look at the Kuramoto model with white noise forces: in mathematical terms it is a set of oscillators, each driven by an independent Brownian motion with a constant drift, that is each oscillator has its own frequency, which, in general, changes from one oscillator to another (these frequencies are usually taken to be random and they may be viewed as a quenched disorder). The interactions between oscillators are of long range type (mean field). We review some results on the Kuramoto model from a statistical mechanics standpoint: we give in particular necessary and sufficient conditions for reversibility and we point out a formal analogy, in the limit, with local mean field models with conservative dynamics (an analogy that is exploited to identify in particular a Lyapunov functional in the reversible set-up). We then focus on the reversible Kuramoto model with sinusoidal interactions in the limit and analyze the stability of the non-trivial stationary profiles arising when the interaction parameter is larger than its critical value . We provide an analysis of the linear operator describing the time evolution in a neighborhood of the synchronized profile: we exhibit a Hilbert space in which this operator has a self-adjoint extension and we establish, as our main result, a spectral gap inequality for every .

2010 Mathematics Subject Classification: 82C20, 35P15, 60K35

Keywords: Synchronization, Kuramoto model, (ir)reversibility, mean field Spin XY model, spectral gap inequality

## 1. Introduction

The Kuramoto model (e.g. [1, 21]) is defined by the set of coupled stochastic differential equations:

 \rm dφξj(t)=ξj\rm dt−KNN∑i=1sin(φξj(t)−φξi(t))\rm d% t+σ\rm dwj(t), (1.1)

for , where

1. is a family of independent and identically distributed standard Brownian motions. We refer to this source of randomness as thermal noise.

2. is a family of independent identically distributed random variables. This is another source of noise, and we refer to it as disorder.

3. is a real parameter and .

We stress from now that we consider the stochastic evolution (1.1) once a realization of the disorder variables is chosen, so the disorder is of quenched type. Moreover the law of and of the initial condition (specified below) does not depend on the values of the disorder variables.

###### Remark 1.1.

It will be at times interesting to discuss the role of the sine drift in the model. We will therefore refer to a -model when is replaced by a smooth, i.e. , -periodic function .

The variables are actually angles, so we focus on , which is an element of . The existence and uniqueness of a unique (strong) solution to the system (1.1), when the initial condition and are independent and are square integrable random variables, is a standard result. In our case we may therefore choose arbitrarily distributed provided that it is concentrated on .

The main result of this work is on the model in which there is no disorder, that is the case in which the law of is degenerate, so that for every , with is a real constant. In this case, with the change of variable we have

 \rm dφj(t)=−KNN∑i=1sin(φj(t)−φi(t))\rm dt+σ\rm dwj(t). (1.2)

Disregarding the disorder is actually a major simplification first of all because, if , the system (1.2) is reversible with respect to the (Gibbs) probability measure

 μN,K(\rm dφ):=1ZN,Kexp(−2Kσ2HN(φ))λN(\rm dφ), (1.3)

where , is the uniform probability measure on (that is the -fold product of Lebesgue measures normalized by ),

 HN(φ):=−12NN∑i=1N∑j=1cos(φj−φi), (1.4)

and is the partition function. In fact, the generator of the dynamics (1.2) acts on twice differentiable functions as

 LK,NF(φ)=σ22N∑i=1∂2F(φ)∂φ2i−KN∑i=1∂HN(φ)∂φi∂F(φ)∂φi, (1.5)

and one directly verifies the symmetry for , which implies that is invariant for the dynamics [20, 13].

The measure is the Gibbs measure of a classical statistical mechanics model: the mean field spin XY model with single spin state space , i.e. mean field plane rotators [18].

###### Remark 1.2.

It is important to notice that the system (1.1) is not reversible unless for every . Even the case for every is not reversible, but, as we have argued, it maps to a reversible system. Notice in fact that, unless , the transformation maps to a system with time dependent interactions. This strongly hints to the absence of reversibility and it is indeed the case, but proving such a statement is more delicate: we address this point in Section 4 below. The aspect that we want to stress here is the disorder induced non-equilibrium character of the full Kuramoto model.

### 1.1. Empirical measure and the large N limit.

Since we focus on the case, there is no loss in generality in choosing and we will do so from now on. We introduce the empirical measure

 νN,t(\rm dθ):=1NN∑j=1δφj(t)(\rm dθ), (1.6)

and observe that, by Itô’s rule, for every and

 ∫SF(θ)νN,t(\rm dθ)−∫SF(θ)νN,0(\rm dθ)=−K∫t0∫S2F′(θ)sin(θ−θ′)νN,s(\rm dθ)νN,s(\rm dθ′)\rm ds+12∫t0∫SF′′(θ)νN,s(\rm dθ)\rm ds+MN,F(t), (1.7)

where is a continuous martingale with quadratic variation at time , , equal to . Therefore, by Doob’s inequality, for every we have that is bounded by : this guarantees that the thermal noise disappears as , so that the limit of the empirical measure, if it exists, is not random. To make this precise we introduce the space , where are the probability measures on equipped with the topology of the weak convergence, and observe that if a subsequence of (of elements of ) converges to a limit , we have that for and every

 ∫SF(θ)νt(\rm dθ)−∫SF(θ)ν0(\rm dθ)=12∫t0∫SF′′(θ)νs(\rm dθ)\rm ds−K∫t0∫S2F′(θ)sin(θ−θ′)νs(\rm dθ)νs(\rm dθ′)\rm ds. (1.8)

This is a weak form of the equation

 (1.9)

when . So that if we assume that converges to a non random limit and if there is a unique solution to (1.8), then the evolution is non random and determined by (1.8).

More precisely, we have the following:

###### Proposition 1.3.

If there exists such that for every and every we have

 limN→∞P(∣∣∣∫SF(θ)νN,0(\rm dθ)−∫SF(θ)ν0(\rm dθ)∣∣∣>ε)=0, (1.10)

then for every we have that for every and

 limN→∞P(∣∣∣∫SF(θ)νN,t(\rm dθ)−∫SF(θ)νt(\rm dθ)∣∣∣>ε)=0, (1.11)

where is the unique solution of (1.8). Moreover, for every the measure is absolutely continuous with respect to the Lebesgue measure with (strictly) positive density and the function , from to , is smooth and solves (1.9).

Proposition 1.3 is a particular (and particularly easy) case of far more general results (see for example [8, 15]). The derivation goes along proving tightness of and then proving uniqueness for the limiting equation (1.8). In our case such an equation is particularly nice and the evolution is smoothing, so that even if the initial datum is not a function (i.e is not absolutely continuous with respect to the Lebesgue measure) or it is not smooth, and for every . These analytic aspects are taken up with more details in Section 3 (Proposition 3.1)

###### Remark 1.4.

It is however important to recall here that Proposition 1.3 can be generalized to cover the disordered case (1.1). We refer to [7, 10] for precise statements, but, roughly, if the law of the random variable is denoted by (let us for example assume that is bounded), the empirical average at time converges as to a measure with density , where is the unique solution to

 ∂tqt(θ;ξ)=12∂2qt(θ;ξ)∂θ2+∂∂θ[(∫R(∫SKsin(θ−θ′)qt(θ′;ξ′)\rm dθ′)μ(\rm dξ′)+ξ)qt(θ;ξ)],q0(θ;ξ)=\rm dν0(\rm dθ)\rm dθ, (1.12)

for every in the support of . We have of course assumed, for simplicity, that is absolutely continuous with respect to the Lebesgue measure.

###### Remark 1.5.

The Kuramoto limit evolution (1.9) comes up also as mesoscopic scaling limit for the density of Kač models with conservative dynamics [5, 9], when the interaction potential is chosen equal to . This is just a formal analogy, but it is going to be crucial for the sequel.

### 1.2. Stationary profiles

By the regularizing character of the evolution (Proposition 3.1), the stationary solutions to (1.8) coincide with the stationary solutions to (1.9) (we are of course interested only in non-negative solutions of total mass equal to one: Proposition 3.1 guarantees also the positivity of stationary solutions). Let us notice moreover that if is a stationary solution, then is a stationary solution too, for arbitrary choice of . This is due to the invariance of (1.2) under rotations (that can of course be read also out of (1.4)). Note that is a solution to (1.2), regardless of the value of but there may be more solutions: in fact every stationary solution can be written as for some and

 ^q(θ):=exp(2Krcos(θ)∫Sexp(2Krcos(θ′)\rm dθ′, (1.13)

with a non-negative solution to

 r:=Ψ(2Kr),    with   Ψ(x):=∫Scos(θ)exp(xcos(θ))\rm dθ∫Sexp(xcos(θ))\rm dθ. (1.14)

In general, there is more than one solution to (1.14): in fact, there can be at most two, more precisely there is only the trivial solution for and there is also a second solution if . This is because and because is strictly concave [16]. In terms of stationary solutions, this means that for only the flat (incoherent) profile is stationary, while for also is a family of stationary solutions (they are the solutions that exhibit the coherence or synchronization of the system).

The result we just stated, that is (1.13)-(1.14), is a classical one in the sense that it is of course closely linked to the solution of the mean field planar rotator model [18] (result completed by the concavity result proven in [16]). It is however worthwhile recalling the proof: every stationary solution of (1.9) satisfies

 12^q′(θ)+K(∫Ssin(θ−θ′)^q(θ′)\rm dθ′)^q(θ)=C, (1.15)

for some constant . Since we know that , then (1.15) yields

 12(log^q(θ))′−K(∫Scos(θ−θ′)^q(θ′)\rm dθ′)′=C^q(θ), (1.16)

which implies . At this point, by playing on the rotation symmetry, we may assume that , so that any stationary non-negative solution with prescribed first Fourier cosine coefficient equal to satisfies

 12^q(θ)′−Krcos(θ)^q(θ)=0. (1.17)

A solution to (1.17) is proportional to . By normalizing ( is a probability density) and recalling the constraint on the first Fourier cosine coefficient we get to (1.13)–(1.14).

###### Remark 1.6.

Remarkably, a generalization of (1.13)–(1.14) holds also in the disordered case [19]. The key to such a derivation, like in the step above, is in the identification of the order parameter , that captures the degree of coherence (or synchronization) of the oscillators. In statistical mechanics terms this is nothing but the fact that the Hamiltonian may be rewritten as , with . However, if one considers an -model (cf. Remark 1.1), the Hamiltonian cannot be expressed any longer as a function of the total magnetization . For the identification of the order parameter in this more general context we refer to [6].

### 1.3. The gradient flow viewpoint

For our purposes the following fact is of crucial importance: (1.9) can be rewritten in the gradient form

 ∂tqt(θ)=∇[qt(θ)∇(δF(qt)δqt(θ))], (1.18)

where we use for for visual impact, is the standard Fréchet derivative of the functional and

 F(q):=12∫Sq(θ)logq(θ)\rm dθ−K2∫S2cos(θ−θ′)q(θ)q(θ′)\rm dθ\rm dθ′. (1.19)

Note that is Fréchet differentiable at for example if is continuous and and Proposition 3.1 guarantees that the evolution may be cast in the form (1.18) for . A direct consequence of (1.18) is that

 ∂F(qt)∂t=−∫Sqt(θ)(∇δF(qt)δqt(θ))2\rm dθ≤0. (1.20)

A first consequence of this observation is:

###### Proposition 1.7.

If there exists such that , then there exists a constant and a value satisfying (1.14) such that for every .

Proof. By Proposition 3.1 (guaranteeing smoothness and positivity of the solution) and by (1.20) we see that is constant for , so that for every and . But this implies that

 12∇(logqt(θ))+K(∫Ssin(θ−θ′)qt(θ′)\rm dθ′)=0, (1.21)

which is precisely (1.16) with . Therefore for every there exists a constant such that , with as in (1.13)-(1.14). Since is a stationary solution, the claim follows. ∎

Two observations are in order:

1. Proposition 1.7 generalizes to the non-disordered -model, when the latter is reversible (see Section 4), in the sense that the hypotheses imply and this condition identifies all the stationary solutions.

2. Proposition 1.7 shows that there is no non-trivial stationary solution to (1.1) when , a non-zero constant. This is simply because, otherwise, we would have a solution to (1.2) of the form , with non-constant, which violates Proposition 1.7. This is of interest also because it is not clear that Proposition 1.7 generalizes to disordered models. Clarifying the link between non-reversibility and coexistence of stationary and rotating solutions appears also to be an intriguing question.

### 1.4. On synchronization stability

The main result that we present addresses the important issue of the stability of the non-trivial stationary profiles , more precisely of the stability of the invariant manifold . In the literature we find a full analysis of incoherence stability [22] (also in presence of disorder) as well as an analysis of synchronized profiles as bifurcation from the incoherent profile (we refer to [1] and the several references therein). Our aim is to have a detailed non-perturbative analysis of the linearized evolution operator in the non disordered case, for every .

To address such an issue we observe that the linearized evolution around obeys the equation with a linear operator with domain defined as

 L^qu(θ)=12Δu(θ)+K∇[^q(θ)∫Ssin(θ−θ′)u(θ′)\rm dθ′+u(θ)∫Ssin(θ−θ′)^q(θ′)\rm dθ′]. (1.22)

It is easy to verify that , and this corresponds to the rotation invariance of the problem. However, what we are going to prove is that the remaining part of the spectrum is also real and it lies on the negative semi-axis. In order to make precise statements about we introduce the Hilbert space of distributions such that , with . Of course the derivative is taken in the sense of distributions and is determined, given , only up to a constant: we remove this uncertainty by stipulating that . The norm of is defined by

 ∥u∥2−1,1/^q:=∫SU(θ)2^q(θ)\rm dθ, (1.23)

and the scalar product of and is going to be denoted by : it is of course equal to , with definition of in analogy with . We will come back in the next section with more on , but what one can verify directly is that and are subsets of and that is symmetric as an operator on , that is

 ⟨⟨u,L^qv⟩⟩=⟨⟨v,L^qu⟩⟩, (1.24)

for every (for an explicit expression see (2.14)). We will actually prove (Proposition 2.6) that is essentially self-adjoint. Moreover:

###### Theorem 1.8.

The spectrum of (the self-adjoint extension of) is pure point and it lies in . The value is in the spectrum, with one-dimensional eigenspace (spanned, as we have seen, by ) and the distance between zero and the rest of the spectrum is of at least

 (1.25)

We stress that Theorem 1.8 holds as soon as there is a non-trivial solution to (1.14), that is for every . We have:

 λ(K)\lx@stackrelK↘1∼K−12    and    λ(K)\lx@stackrelK→∞∼exp(−8K+2)4K. (1.26)

Numerically increases till , where it reaches the value , and then it decreases.

The paper is organized as follows: in Section 2 we study and prove a spectral gap inequality, the essential self-adjontness of the operator and the fact that the spectrum is pure point. The nonlinear evolution properties mentioned in this introduction are treated in Section 3 and the (ir)reversibility issues are considered in Section 4.

## 2. Synchronization stability

In this section we prove the main result (Theorem 1.8). We assume and, for simplicity, we drop the hat from , so that a stationary solution is denoted by .

### 2.1. Some properties of the stationary profile

We first rewrite (1.13)-(1.14) by using the Bessel function notation:

 q(θ):=12πI0(2Kr)exp(2Krcos(θ)), (2.1)

and is the unique positive solution of

 r:=Ψ(2Kr),   with  Ψ(x):=I1(x)I0(x). (2.2)

We have used the standard notation for the modified Bessel functions of order and , explicitly

 Iν(x):=12π∫2π0(cos(θ))νexp(xcos(θ))\rm dθ,        for ν=0,1. (2.3)

As already mentioned before, uniqueness of is a non-trivial fact that follows from [16, Lemma 4], that establishes in particular the concavity of on the positive semi-axis. One can therefore define, via (2.2)), the function (one sets by continuity). We have that, for , or (equivalently)

 √1−1K

These bounds are easily checked for close to and large, and and the numerical plots of the three functions appearing in (2.4) ca be found in Figure 1. We could not find quick proofs of (2.4): a proof of the upper bound is a byproduct of one of the arguments that we develop below (see the proof of Lemma 2.2), while we prove here the lower bound by using Bessel functions properties.

Proof of (2.4), lower bound. By the change of variables we see that what we have to prove is equivalent to showing that

 Ψ2(y)+2yΨ(y)−1\lx@stackrely>0>0, (2.5)

Apply now the identity [23]

 I1(y)I0(y)=y2⎛⎜ ⎜⎝11+y2I2(y)I1(y)⎞⎟ ⎟⎠, (2.6)

so that

 1−2yΨ(y)=11+2yI1(y)I2(y)(>0), (2.7)

and therefore (2.5) is equivalent to

 Ψ2(y)(1+2yI1(y)I2(y))=I1(y)2I0(y)I2(y)>1, (2.8)

where the intermediate step follows from the identity [23]. But this is equivalent to for , a fact that is proven in [11, (3.8)]. ∎

In the sequel we will also use the notations

 J(θ):=−Ksin(θ)    and    ˜J(θ):=Kcos(θ), (2.9)

so that

 q′2q=J∗q, (2.10)

and with this change of notation (1.22) reads

 Lqu=12u′′−(q(J∗u)+u(J∗q))′. (2.11)

### 2.2. Rigged Hilbert spaces and H−1,1/q

We now introduce a rigged Hilbert spaces structure [4, pp. 81-82]. The pivot (Hilbert) space is (of course the scalar product is and the norm is denoted by ). The second Hilbert space we consider is , closure of the set of periodic functions such that with respect to the norm

 ∥u∥1,q:=√∫S(u′)2q, (2.12)

so that and the canonical injection of into is continuous (by the Poincaré inequality). Note that is dense in . We consider then the dual space of and the duality functional in defined by for every given (). We can define by setting . One can then show that is dense in and injects into in a continuous way [4, p. 82]. This injection allows considering as a subset of , by identifying and . Moreover if we have that can be made more explicit: given we call the primitive of such that and we observe that

 ∥u∥V′=supv∈V(u,v)∥v∥V=supv∈H1,q∫Uv′√∫q(v′)2=√∫U2q, (2.13)

where the last step follows on one hand by the Cauchy-Schwarz inequality (establishing the upper bound) and by choosing in the supremum (establishing the lower bound).

As already mentioned in the introduction, the scalar product in is denoted by . We observe also that these steps allow the precise identification of the functions in : if and only if (in the sense of distributions), with .

At this point it is crucial to observe that (recall that is the subset of periodic functions such that ) is dense in and that for , we have (use (2.10))

 ⟨⟨v,Lqu⟩⟩=⟨⟨Lqv,u⟩⟩=−12∫Suvq+(v,˜J∗u). (2.14)

In words: is a symmetric operator on .

### 2.3. Estimates on the Dirichlet form

It is going to be useful to introduce also the Hilbert space , which coincides with as a set of functions, but we equip it with the scalar product

 ⟨u,u⟩:=(u,u/q). (2.15)

We therefore introduce, for , the Dirichlet form . By (2.14) we have

 D(u)=12⟨u,u⟩−(u,˜J∗u). (2.16)

Our aim is to bound from below and we start with two technical lemmas. The first one yields the spectral decomposition of , viewed as a quadratic form on .

###### Lemma 2.1.

We have the orthogonal decomposition

 L21/q=V0⊕⊥V1/2⊕⊥VK−1/2, (2.17)

where

 V0:={θ↦a0+∑j≥2(ajcos(jθ)+bjsin(jθ)):∑ja2j+b2j<∞}, (2.18)

and both and are one dimensional subspaces generated respectively by and by . Moreover, when we have

 ˜J∗u=λqu. (2.19)

Proof. The -orthogonality statements and follow directly from the orthogonality in of the family . Instead because .

The validity of (2.19) follows by direct computation: for

 ˜J∗u(θ)=Kcos(θ)∫2π0cos(θ′)u(θ′)\rm dθ′+Ksin(θ)∫2π0sin(θ′))u(θ′)\rm dθ′, (2.20)

which is equal to zero because does not contain the first harmonics. The other two cases follow by using the same trigonometric identity and the following two (clearly equivalent) identities:

 ∫2π0q(θ)sin2(θ)\rm dθ=12K,      ∫2π0q(θ)cos2(θ)\rm dθ=(1−12K). (2.21)

The second lemma is more technical and its interest will become clear in the proof of Proposition 2.3.

###### Lemma 2.2.

We have

 (2.22)

where

 ^u(θ)=−1−2π(1−12K)(r2−(1−12K))q(θ)+2πr(r2−(1−12K))q(θ)cos(θ), (2.23)

Proof. We have to minimize a quadratic functional under the linear constraints and . This corresponds to the three constraints:

 ∫2π0u(θ)\rm dθ=0,     ∫2π0cos(θ)u(θ)\rm dθ=0   and   ∫2π0sin(θ)u(θ)\rm dθ=0. (2.24)

The extrema (minima, by convexity) of such a problem can be found by the Lagrange multipliers method and they are of the form

 ^u(θ)=−1+λq(θ)+μq(θ)cos(θ)+ηq(θ)sin(θ), (2.25)

with , and three real numbers. The constraints (2.24), via (2.2) and (2.21), yield the linear system , and , which has a solution if and only if . Since the minimum exists for every and since for we see that for every (this is the upper bound in (2.4)). The proof is completed by making and explicit and using that the expression in (2.22) is equal to . ∎

The following is one of our main statements:

###### Proposition 2.3.

There exists such that, if is such that , then

 D(u)≥cK⟨u−u1/2,u−u1/2⟩, (2.26)

where . In particular, .

Proof. Lemma 2.1 shows that, if we write (according to the decomposition (2.17): of course ) we have

 (2.27)

We write and and, with this notations, one directly sees, using the definitions (2.1) and (2.2), that the constraint is equivalent to

 r˜c=−2πa0. (2.28)

Let us observe that we can assume : if