The mean field analysis for the Kuramoto model on graphs I.
The mean field equation and transition point formulas
Hayato Chiba∗ and Georgi S. Medvedev†
December 20, 2016

Abstract
In his classical work on synchronization, Kuramoto derived the formula for the critical value of thecoupling strength corresponding to the transition to synchrony in large ensembles of all-to-all coupledphase oscillators with randomly distributed intrinsic frequencies. We extend the Kuramoto’s result to alarge class of coupled systems on convergent families of deterministic and random graphs. Specifically,
we identify the critical values of the coupling strength (transition points), between which the incoherent
state is linearly stable and is unstable otherwise. We show that the transition points depend on the largest
positive or/and smallest negative eigenvalue(s) of the kernel operator defined by the graph limit. Thisreveals the precise mechanism, by which the network topology controls transition to synchrony in the
Kuramoto model on graphs. To illustrate the analysis with concrete examples, we derive the transition
point formula for the coupled systems on Erdős-Rényi, small-world, and k-nearest-neighbor families ofgraphs. As a result of independent interest, we provide a rigorous justification for the mean field limit
for the Kuramoto model on graphs. The latter is used in the derivation of the transition point formulas.

1 Introduction
Synchronization of coupled oscillators is a classical problem of nonlinear science with diverse applicationsin science and engineering [3, 38]. Physical and technological applications of synchronization include
power, sensor, and communication networks [9], mobile agents [35], electrical circuits [1], coupled lasers[18], and Josephson junctions [44], to name a few. In biological and social sciences, synchronization is
studied in the context of flocking, opinion dynamics, and voting [14, 30]. Synchronization plays a prominentrole in physiology and in neurophysiology, in particular. It is important in the information processing in thebrain [36] and in the mechanisms of several severe neurodegenerative diseases such as epilepsy [43] andParkinson s Disease [22]. This list can be continued.
Identifying common principles underlying synchronization in such diverse models is a challenging task.
In seventies Kuramoto found an elegant approach to this problem. Motivated by problems in statistical∗

Institute of Mathematics for Industry, Kyushu University / JST PRESTO, Fukuoka, 819-0395, Japan,
chiba@imi.kyushu-u.ac.jp†
Department of Mathematics, Drexel University, 3141 Chestnut Street, Philadelphia, PA 19104, medvedev@drexel.edu

1physics and biology, he reduced a system of weakly coupled limit cycle oscillators to the system of equations
for the phase variables only1 . The resultant equation is called the Kuramoto model (KM) [20]. Kuramoto’s
method applies directly to a broad class of models in natural science. Moreover, it provides a paradigm forstudying synchronization. The analysis of the KM revealed one of the most striking results of the theoryof synchronization. For a system of coupled oscillators with randomly distributed intrinsic frequencies,
Kuramoto identified the critical value of the coupling strength, at which the gradual buildup of coherence
begins. He introduced the order parameter, which describes the degree of coherence in a coupled system.
Using the order parameter, Kuramoto predicted the bifurcation marking the onset of synchronization.
Kuramoto’s analysis, while not mathematically rigorous, is based on the correct intuition for the transition to synchronization. His discovery initiated a line of fine research (see [40, 41, 39, 5] and references
therein). It was shown in [40, 41] that the onset of synchronization corresponds to the loss of stability of
the incoherent state, a steady state solution of the mean field equation. The latter is a nonlinear hyperbolic
partial differential equation for the probability density function describing the distribution of phases on theunit circle at a given time. The bifurcation analysis of the mean field equation is complicated by the presenceof the continuous spectrum of the linearized problem on the imaginary axis. To overcome this problem, in[5] Chiba developed an analytical method, which uses the theory of generalized functions and rigged Hilbert
spaces [12].
The Kuramoto’s result and the analysis in [40, 41, 5] deal with all-to-all coupled systems. Real world
applications feature complex and often random connectivity patterns [34]. The goal of our work is to describe the onset of synchronization in the KM on graphs. To this end, we follow the approach in [27, 28].
Specifically, we consider the KM on convergent families of simple and weighted graphs, for which we derive and rigorously justify the mean field limit. Our framework covers many random graphs widely used
in applications, including Erdős-Rényi and small-world graphs. With the mean field equation in hand, we
derive the transition point formulas for the critical values of the coupling strength, where the incoherent
state loses stability. Thus, we identify the region of linear stability of the incoherent state. In the follow-up
work, we will show that in this region the incoherent state is asymptotically stable (albeit in weak topology),
analyze the bifurcations at the transition points, and develop the center manifold reduction.
The original KM with all-to-all coupling and random intrinsic frequencies has the following form:K∑sin(θj − θi ).
θ̇i = ωi +
n
n

(1.1)

j=1

Here, θi : R → S := R/2πZ, i ∈ [n] := {1, 2, . . . , n} is the phase of the oscillator i, whose intrinsic
frequency ωi is drawn from the probability distribution with density g(ω), n is the number of oscillators,
and K is the strength of coupling. The sum on the right-hand side of (1.1) describes the interactions of theoscillators in the network. The goal is to describe the distribution of θi (t), i ∈ [n], for large times andn ≫ 1.
Since the intrinsic frequencies are random, for small values of the coupling strength K > 0, the dynamics of different oscillators in the network are practically uncorrelated. For increasing values of K > 0,
however, the dynamics of the oscillators becomes more and more synchronized. To describe the degree of1

For related reductions predating Kuramoto’s work, see [25, 26, 3] and the discussion before Theorem 9.2 in [15].

2synchronization, Kuramoto used the complex order parameter:
r(t)e

iψ(t)

:= n

−1

n
∑

eiθj (t) .

(1.2)

j=1

Here, 0 ≤ r(t) ≤ 1 and ψ(t) stand for the modulus and the argument of the order parameter defined by
the right-hand side of (1.2). Note that if all phases are independent uniform random variables and n ≫ 1then with probability 1, r = o(1) by the Strong Law of Large Numbers. If, on the other hand, all phasevariables are equal then r = 1. Thus, one can interpret the value of r as the measure of coherence in thesystem dynamics. Numerical experiments with the KM (with normally distributed frequencies ωi ’s) reveal
the phase transition at a certain critical value of the coupling strength Kc > 0. Specifically, numerics suggestthat for t ≫ 1 (cf. [39]){O(n−1/2 ),0 < K < Kc ,r(t) =
r∞ (K) + O(n−1/2 ), K > Kc .
Assuming that g is a smooth even function that is decreasing on ω ∈ R+ , Kuramoto derived the formula forthe critical value
2
.
(1.3)
Kc =
πg(0)
Furthermore, he formally showed that in the partially synchronized regime (K > Kc ), the steady-state value
for the order parameter is given by√−16 √r∞ (K) =K − Kc + O(K − Kc ).
(1.4)
πKc4 g ′′ (0)Recently, Chiba and Nishikawa [7] and Chiba [5] confirmed Kuramoto’s heuristic analysis with the rigorous
derivation of (1.3) and analyzed the bifurcation at Kc .
In this paper, we initiate a mathematical investigation of the transition to coherence in the Kuramotomodel on graphs. To this end, we consider the following model:K∑θ̇i = ωi +Wnij sin(θj − θi ),
n
n

(1.5)

j=1

Wn = (Wnij ) is an n × n symmetric matrix of weights. Note that in the classical KM (1.1), every oscillator
is coupled to every other oscillator in the network, i.e., the graph describing the interactions between theoscillators is the complete graph on n nodes. In the modified model (1.5), we supply the edges of thecomplete graph with the weights (Wnij ). Using this framework, we can study the KM on a variety ofdeterministic and random (weighted) graphs. For instance, let Wnij , 1 ≤ i < j ≤ n be independent
Bernoulli random variables
P(Wnij = 1) = p,for some p ∈ (0, 1). Complete the definition of Wn by setting Wnji = Wnij and Wnii = 0, i ∈ [n]. With
this choice of Wn , (1.5) yields the KM on Erdős-Rényi random graph.
The family of Erdős-Rényi graphs parameterized by n is one example of a convergent family of randomgraphs [23]. The limiting behavior of such families is determined by a symmetric measurable function on3

the unit square W (x, y), called a graphon. In the case of the Erdős-Rényi graphs, the limiting graphon isthe constant function W ≡ p. In this paper, we study the KM on convergent families of deterministic andrandom graphs. In each case the asymptotic properties of graphs are known through the limiting graphon
W . The precise relation between the graphon W and the weight matrix Wn will be explained below.
In studies of coupled systems on graphs, one of the main questions is the relation between the structureof the graph and network dynamics. For the problem at hand, this translates into the question of how thestructure of the graph affects the transition to synchrony in the KM. For the KM on convergent families ofgraphs, in this paper, we derive the formulas for the critical values
Kc+ =2
πg(0)ζmax (W)and Kc− =2,πg(0)ζmin (W)(1.6)

where ζmax (W) (ζmin (W)) is the largest positive (smallest negative) eigenvalue of the self-adjoint kernel
operator W : L2 (I) → L2 (I), I := [0, 1], defined by∫
W[f ] = W (·, y)f (y)dy, f ∈ L2 (I).
(1.7)
I

If all eigenvalues of W are positive (negative) then Kc− := −∞ (Kc+ := ∞). The main result of this workshows that the incoherent state is linearly (neutrally) stable for K ∈ [Kc− , Kc+ ] and is unstable otherwise.
The transition point formulas in (1.6) reveal the effect of the network topology on the synchronization
properties of the KM through the extreme eigenvalues ζmax (W) and ζmin (W). For the classical KM (W ≡1) ζmax (W) = 1 and there are no negative eigenvalues. Thus, we recover (1.3) from (1.6).
We derive (1.6) from the linear stability analysis of the mean field limit of (1.5). The latter is a partial differential equation for the probability density function corresponding to the distribution of the phasevariables on the unit circle (see (2.10), (2.11)). For the classical KM, the mean field limit was derived by
Strogatz and Mirollo in [40]. We derive the mean field limit for the KM on weighted graphs (1.5) andshow that its solutions approximate probability distribution of the phase variables on finite time intervals
for n ≫ 1. Here, we rely on the theory of Neunzert developed for the Vlasov equation [31, 32] (see also
[4, 8, 13]), which was also used by Lancellotti in his treatment of the mean field limit for the classical KM[21].
With the mean field limit in hand, we proceed to study transition to coherence in (1.5). As for theclassical KM, the density of the uniform distribution is a steady state solution of the mean field limit. The
linear stability analysis in Section 3 shows that the density of the uniform distribution is neutrally stable
for Kc− ≤ K ≤ Kc+ and is unstable otherwise. Thus, the critical values Kc± given in (1.6) mark theloss of stability of the incoherent state. The bifurcations at Kc± and the formula for the order parametercorresponding to (1.4) will be analyzed elsewhere using the techniques from [5, 6].
Sections 4 and 5 deal with applications. In the former section we collect approximation results, whichfacilitate application of our results to a wider class of models. Further, in Section 5, we discuss the KM forseveral representative network topologies : Erdős-Rényi, small-world, k-nearest-neighbor graphs, and theweighted ring model. We conclude with a brief discussion of our results in Section 6.

4

2

The mean field limit

Throughout this paper, we will use a discretization of I = [0, 1] :Xn = {ξn1 , ξn2 , . . . , ξnn },ξni ∈ I, i ∈ [n],(2.1)∀f ∈ C(I).

(2.2)

which satisfies the following property
lim n

n→∞

−1

n
∑

∫
f (ξni ) =f (x)dx,I

i=1

Example 2.1. The following two examples of Xn will be used in constructions of various graphs throughout
this paper.
1. The family of sets (2.1) with ξni = i/n, i ∈ [n] satisfies (2.2).
2. Let ξ1 , ξ2 , . . . be independent identically distributed (IID) random variables (RVs) with ξ1 having theuniform distribution on I. Then with Xn = {ξ1 , ξ2 , . . . , ξn } (2.2) holds almost surely (a.s.), by theStrong Law of Large Numbers.
Let W be a symmetric Lipschitz continuous function on I 2 :√|W (x1 , y1 ) − W (x2 , y2 )| ≤ LW (x1 − x2 )2 + (y1 − y2 )2∀(x1,2 , y1,2 ) ∈ I 2 .

(2.3)

The weighted graph Γn = G(W, Xn ) on n nodes is defined as follows. The node and the edge sets of
Γn are V (Γn ) = [n] andE(Γn ) = {{i, j} : W (ξni , ξnj ) ̸= 0, i, j ∈ [n]} ,(2.4)

respectively. Each edge {i, j} ∈ E(Γn ) is supplied with the weight Wnij := W (ξni , ξnj ).
On Γn , we consider the KM of phase oscillatorsθ̇ni = ωi + Kn

−1

n
∑

Wnij sin(θnj − θni ),i ∈ [n].

(2.5)

j=1

The phase variable θni : R → S := R/2πZ corresponds to the oscillator at node i ∈ [n]. Throughout thispaper, we identify θ ∈ S with its value in the fundamental domain, i.e., θ ∈ [0, 2π). Further, we equip S
with the distancedS (θ, θ′ ) = min{|θ − θ′ |, 2π − |θ − θ′ |}.
(2.6)
The oscillators at the adjacent nodes interact through the coupling term on the right hand side of (2.5).
The intrinsic frequencies ω1 , ω2 , . . . are IID RVs. Assume that ω1 has absolutely continuous probability
distribution with a continuous density g(ω). The initial condition
θni (0) = θi0 ,5i ∈ [n],(2.7)

are sampled independently from the conditional probability distributions with densities ρ̂0θ|ω (θ, ωi , ξni ), i ∈[n]. Here, ρ̂0θ|ω (θ, ω, ξ) is a nonnegative continuous function on G := S × R × I that is uniformly continuousin ξ, i.e., ∀ϵ > 0 ∃δ > 0 such that(ξ1 , ξ2 ∈ I) & |ξ1 − ξ2 | < δ =⇒

ρ̂0θ|ω (θ, ω, ξ1 ) − ρ̂0θ|ω (θ, ω, ξ2 ) < ϵuniformly in (θ, ω) ∈ S × R. In addition, we assume∫ρ̂0θ|ω (ϕ, ω, ξ)dϕ = 1S∀(ω, ξ) ∈ R × I.

(2.8)

(2.9)

We want to show that the dynamics of (2.5) subject to the initial condition (2.7) can be described interms of the probability density function ρ̂(t, θ, ω, x) satisfying the following Vlasov equation∂∂ρ̂(t, θ, ω, x) +{ρ̂(t, θ, ω, x)V (t, θ, ω, x)} = 0,∂t∂θ(2.10)

∫ ∫ ∫

where
V (t, θ, ω, x) = ω + K

I

R

S

W (x, y) sin(ϕ − θ)ρ̂(t, ϕ, λ, y)dϕdλdy(2.11)and the initial conditionρ̂(0, θ, ω, x) = ρ̂0θ|ω (θ, ω, x)g(ω).
By (2.9), ρ̂(0, θ, ω, x) is a probability density on (G, B(G)):}∫ ∫ ∫∫ ∫ {∫
ρ̂(0, θ, ω, x)dxdωdθ =ρ̂0θ|ω (θ, ω, x)dθ g(ω)dxdω
S R I
S
∫R I
=
g(ω)dω = 1.

(2.12)(2.13)

R

Here, B(G) stands for the Borel σ-algebra of G.
Below, we show that the solutions of the IVPs for (2.5) and (2.10), generate two families of Borel
probability measures parametrized by t > 0. To this end, we introduce the following empirical measureµnt (A) = n−1n∑δPni (t) (A)A ∈ B(G),(2.14)

i=1

where Pni (t) = (θni (t), ωi , ξni ) ∈ G.
To compare measures generated by the discrete and continuous systems, following [32], we use thebounded Lipschitz distance:∫∫(2.15)
f dµ −
f dν , µ, ν ∈ M,d(µ, ν) = sup
f ∈L

G

G

where L is the set of functions
L = {f : G → [0, 1] : |f (P ) − f (Q)| ≤ dG (P, Q), P, Q ∈ G}
6

(2.16)

and M stands for the space of Borel probability measures on G. Here,√dG (P, P ′ ) = dS (θ, θ′ )2 + (ω − ω ′ )2 + (x − x′ )2 ,for P = (θ, ω, x) and P ′ = (θ′ , ω ′ , x′ ). The bounded Lipschitz distance metrizes the convergence of Borel
probability measures on G [10, Theorem 11.3.3].
We are now in a position to formulate the main result of this section.
Theorem 2.2. Suppose W is a Lipschitz continuous function on I 2 . Then for any T > 0, there exists a
unique weak solution2 of the IVP (2.10), (2.11), and (2.12), ρ̂(t, ·), t ∈ [0, T ], which provides the density
for Borel probability measure on G:∫µt (A) =ρ̂(t, P )dP, A ∈ B(G),(2.17)
A

parametrized by t ∈ [0, T ]. Furthermore,d(µnt , µt ) → 0(2.18)

uniformly for t ∈ [0, T ], provided d(µn0 , µ0 ) → 0 as n → ∞.
Proof. We rewrite (2.5) as follows
θ̇ni = λi + Kn

−1

n
∑

W (xni , xnj ) sin(θnj − θni ),(2.19)

j=1

λ̇ni = 0,ẋni = 0,i ∈ [n],subject to the initial condition(θni (0), λni (0), xni (0)) = (θi0 , ωi , ξni ),i ∈ [n].

(2.20)As before, we consider the empirical measure corresponding to the solutions of (2.19), (2.20)µnt (A) = n−1n∑δPni (t) (A),A ∈ B(G),(2.21)

i=1

where Pni (t) = (θni (t), λni (t), xni (t)) ∈ G.
We need to show that µnt and µt are close for large n. This follows from the Neunzert’s theory [32].
Specifically, below we showd(µnt , µt ) ≤ Cd(µn0 , µ0 ), t ∈ [0, T ],(2.22)
for some C > 0 independent from n.
Below, we prove (2.22). Theorem 2.2 will then follow.
2

See [31, Remark 1] for the definition of the weak solution.

7

Let C(0, T ; M) denote the space of weakly continuous M-valued functions on [0, T ]. Specifically,µ. ∈ C(0, T ; M) means that∫t 7→

f (P )dµt (P )(2.23)

Gis a continuous function of t ∈ [0, T ] for every bounded continuous function f ∈ Cb (G).
For a given ν. ∈ C(0, T ; M), consider the following equation of characteristics:
dP= Ṽ [ν. ](t, P ),dt

P (s) = P 0 ∈ G.

(2.24)

where P = (θ, ω, x) and∫ω + K G W (x, y) sin(v − θ)dνt (v, ω, y).
0Ṽ [ν. ](t, P ) = 
0(2.25)

Under our assumptions on W , (2.24) has a unique global solution, which depends continuously on initial−1
data. Thus, (2.25) generates the flow Tt,s : G → G (Ts,s = id, Ts,t = Tt,s):
P (t) = Tt,s [ν. ]P 0 .
Following [31], we consider the fixed point equation:νt = ν0 ◦ T0,t [ν. ],t ∈ [0, T ],(2.26)

which is interpreted as
νt (A) = ν0 (T0,t [ν. ](A))∀A ∈ B(G).

It is shown in [31] that under the conditions (I) and (II) given below, for any ν0 ∈ M there is a uniquesolution of the fixed point equation (2.26) ν. ∈ C(0, T ; M). Moreover, for any two initial conditions ν01,2 ∈M, we havesup d(νt1 , νt2 ) ≤ exp{CT }d(ν01 , ν02 )(2.27)
t∈[0,T ]for some C > 0. By construction of Tt,s and (2.21), the empirical measure µn. satisfies the fixed point
equation (2.26). By [31, Theorem 1], νt , the solution of the (2.26), is an absolutely continuous measure withdensity ρ̂(t, ·) for every t ∈ [0, T ], provided ν0 is absolutely continuous with density ρ̂(0, ·) (cf. (2.13)).
Furthermore, ρ̂(t, P ) is a weak solution of the IVP for (2.10), (2.11), and (2.12). Therefore, since boththe empirical measure µn. and its continuous counterpart µ. (cf. (2.21) and (2.17)) satisfy the fixed point
equation (2.26), we can use (2.27) to obtain (2.22). It remains to verify the following two conditions on
the vector field Ṽ [ν. ], which guarantee the solvability of (2.26) and continuous dependence on initial data
estimate (2.27) (cf. [32]):(I) Ṽ [µ. ](t, P ) is continuous in t and is globally Lipschitz continuous in P with Lipschitz constant3 L1 ,which depends on W .
3

A straightforward estimate shows that L1 =√ ()
2 LW + ∥W ∥L∞ (I 2 ) + 1, where LW is the Lipschitz constant in (2.3).

8

(II) The mapping Ṽ : µ. 7→ Ṽ [µ. ] is Lipschitz continuous in the following sense:Ṽ [µ. ](t, P ) − Ṽ [ν. ](t, P ) ≤ L2 d(µt , νt ),for some L2 > 0 and for all µ. , ν. ∈ C(R, M) and (t, P ) ∈ [0, T ] × G. 4
For the Lipschitz continuous function W , it is straightforward to verify conditions (I) and (II). In particular, (I) follows from the weak continuity of µt (cf. (2.23)) and Lipschitz continuity of W and sin x. The
second condition is verified following the treatment of the mechanical system presented in [32] (see also[21]). We include the details of the verification of (II) for completeness.
Let P = (θ, ω, x) ∈ G be arbitrary but fixed and definef (ϕ, λ, y; P ) =W (x, y) sin(ϕ − θ) + ∥W ∥L∞ (I 2 ),2(∥W ∥L∞ (I 2 ) + LW )(2.28)

where LW is the Lipschitz constant of W (x, y) (cf. (2.3)). Then f ∈ L (cf. (2.15)). Further,∫Ṽ [ν. ](t, P ) − Ṽ [µ. ](t, P ) = K
W (x, y) sin(ϕ − θ) (dνt (ϕ, λ, y) − dµt (ϕ, λ, y))G∫= 2K(∥W ∥L∞ (I 2 ) + LW )f (ϕ, λ, y) (dνt (ϕ, λ, y) − dµt (ϕ, λ, y))G≤ L2 d(νt , µt ),L2 := 2K(∥W ∥L∞ (I 2 ) + LW ),which verifies the condition (II).
Corollary 2.3. For the empirical measure µnt (2.14) and absolutely continuous measure µt (2.17) definedon the solutions of the IVPs (2.5), (2.7), (2.9) and (2.10), (2.11), (2.12) respectively, we have
lim sup d(µnt , µt ) = 0

n→∞ t∈[0,T ]a.s..

Proof. In view of Theorem 2.2, we need to show
lim d(µn0 , µ0 ) = 0n→∞By [10, Theorem 11.3.3], it is sufficient to show∫limf (dµn0 − dµ0 ) = 0
n→∞ G

a.s..

∀f ∈ BL(G) a.s.,(2.29)

where BL(G) stands for the space of bounded real-valued Lipschitz functions on G with the supremum
norm. Since BL(G) is a separable space, let {fm }∞m=1 denote a dense set in BL(G). Using (2.7) and(2.14), we have
∫
n
n
∑
∑
fm dµn0 = n−1
fm (θi0 , ωi , ξni ) =: n−1
Ym,ni .
(2.30)
G
4

i=1

i=1

With these assumptions the estimate (2.27) holds with C = L1 + L2 (see [32]).

9

RVs Ym,ni , i ∈ [n], are independent and uniformly bounded. Further,
(
)
n
n ∫ ∫
∑
∑
−1
fm (ϕ, λ, ξni )ρ̂0θ|ω (ϕ, λ, ξni )g(λ)dλdϕ
E nYm,ni = n−1
i=1
−1

=: n

i=1 S
n
∑

R

(2.31)

Fm (ξni ),i=1

Because fm is Lipschitz and ρ0θ|ω is uniformly continuous in ξ (cf. (2.8)), the function∫ ∫Fm (ξ) =fm (ϕ, λ, ξ)ρ̂0θ|ω (ϕ, λ, ξ)g(λ)dλdϕ
SR

is continuous on I.
By (2.2), we have(−1

lim E n

n→∞

n
∑

)
Ym,ni

∫

∫

=

Fm (ξ)dξ =I

i=1

fm dµ0 .

By the Strong Law of Large Numbers5 , from (2.30), (2.31), and (2.32) we have
∫
n
∑
n
−1
lim
fm (dµ0 − dµ0 ) = lim n(Ym,ni − E Ym,ni ) = 0
n→∞ G

n→∞

Therefore,

{
P

lim

a.s..

(2.33)

i=1

}

∫

n→∞ G

(2.32)Gfm (dµn0 − dµ0 ) = 0 ∀m ∈ N= 1.

Using density of {fm }∞m=1 in BL(G), we have
{
}
∫
n
P lim
f (dµ0 − dµ0 ) = 0 ∀f ∈ BL(G) = 1.
n→∞ G

3

Linear stability

In the previous section, we established that the Vlasov equation (2.10), approximates discrete system (2.5)
for n ≫ 1. Next, we will use (2.10) to characterize the transition to synchrony for increasing K. To thisend, in this section, we derive the linearized equation about the incoherent state, the steady state solutionof the mean field equation. In the next section, we will study how the spectrum of the linearized equation
changes with K.
In this section, we assume that probability density g is a continuous and even function monotonicallydecreasing on R+ .
It is easy to adjust the proof of Theorem 6.1 [2] so that it applies to the triangular array Yni , i ∈ [n], n ∈ N (see [29,Lemma 3.1]).
5

10

3.1

Linearization

First, setting ρ̂(t, θ, ω, x) = ρ(t, θ, ω, x)g(ω), from (2.10) we derive the equation for ρ:∂∂ρ+{Vρ ρ} = 0,
∂t
∂θ

(3.1)

∫ ∫ ∫

where
Vρ (t, θ, ω, x) = ω + K

I

S

R

W (x, y) sin(ϕ − θ)ρ(t, ϕ, λ, y)g(λ)dλdϕdy.

By integrating (2.10) over S, one can see that∫∂
ρ̂(t, θ, ω, x)dθ = 0,∂t S
and, thus,∫∫ρ̂(t, θ, ω, x)dθ =

S

Thus,

∫
ρ̂(0, θ, ω, x)dθ =SSρ̂0θ|ω (0, θ, ω, x)g(ω)dθ = g(ω).

(3.2)

∫
S

ρ(t, θ, ω, x)dθ = 1 ∀(t, ω, x) ∈ R+ × R × I.

(3.3)

∫

In addition,g(ω)dω = 1,Rbecause g is a probability density function. The density of the uniform distribution on S, ρu = (2π)−1 , isa steady-state solution of the mean field equation (3.1). It corresponds to the completely mixed state. We
are interested in stability of this solution. In the remainder of this section, we derive the linearized equation
around ρu .
Let
ρ = ρu + z(t, θ, ω, x).
By (3.3),(3.4)∫
z(t, θ, ω, x)dθ = 0
S

By plugging in (3.4) into (3.1), we obtain∂∂z(t, θ, ω, x) +∂t
∂θ∀(t, ω, x) ∈ R+ × R × I.

{
(
)}
1
Vρu +z (t, θ, ω, x)
+z
= 0.
2π

(3.5)

(3.6)

The expression in the curly brackets has the following form:
(
) (
∫ ∫ ∫
1
Vρu +z (t, θ, ω, x)+z = ω+K
W (x, y) sin(ϕ − θ)×
2π
I S R
(
)
(
)
) 1
ω
−1
× (2π) + z(t, ϕ, λ, y) g(λ)dλdϕdy
+z =
+ ωz
2π
2π
∫ ∫ ∫
K
W (x, y) sin(ϕ − θ)z(t, ϕ, λ, y)g(λ)dλdϕdy + O(z 2 ).
+2π I S R
11

Thus,
∂
∂
z+
∂t
∂θ

{
}
K
ωz +
G[z] + O(z 2 ) = 0,
2πwhere the linear operator G is defined by∫ ∫ ∫G[z] :=W (x, y) sin(ϕ − θ)z(t, ϕ, λ, y)g(λ)dλdϕdy
I

S

R

(3.7)

We arrive at the linearized equation:
∂
∂
Z+
∂t
∂θ

3.2

{

}
K
ωZ +
G[Z] = 0.
2π

(3.8)

Fourier transform

We expand Z into Fourier seriesZ(t, θ, ω, x) =∞∑Ẑk (t, ω, x)e−ikθ +(∞∑)
Ẑk (t, ω, x)e−ikθ ,

(3.9)

k=1

k=1

where Ẑk stands for the Fourier transform of Z
Ẑk =

1
2π

∫
S

Z(θ, ·)eikθ dθ.

In (3.9), we are using the fact that Z is real and Ẑ0 = 0 (cf. (3.5)).
The linear stability of ρu is, thus, determined by the time-asymptotic behavior of Ẑk , k ≥ 1. To derive
the differential equations for Ẑk , k ≥ 1, we apply the Fourier transform to (3.8):
∂
Ẑk +
∂t

(

∂
∂θ

{ \
})
K
ωZ +
G[Z]
= 0.
2π
k

(3.10)

Using the definition of the Fourier transform and integration by parts, we have
(

∂
∂θ

{ \
})
∫
K
1
∂
K \
ωZ +
G[Z]
=
(. . . ) eikθ dθ = −ikω Ẑk − ik (G[Z])k.
2π2π S ∂θ2πk\ , k ≥ 1. To this, end we rewrite
It remains to compute (G[Z])
k
∫ ∫ ∫
(
)
1
G[Z] =W (x, y) ei(ϕ−θ) − e−i(ϕ−θ) Z(t, ϕ, λ, y)g(λ)dλdϕdy
2i I S R∫ ∫()π
W (x, y) e−iθ Ẑ1 − eiθ Ẑ−1 g(λ)dλdy.
=
i I R
12

(3.11)(3.12)

Using (3.12), we compute
∫ ∫ ∫
(
)
1 π
\
(G[Z])k =W (x, y) ei(k−1)θ Ẑ1 − ei(k+1)θ Ẑ−1 g(λ)dλdydθ
2π i{ π S I R
i P[Ẑ1 ], k = 1,=0,
k > 1,(3.13)∫ ∫where
P[Z] :=W (x, y)Z(t, λ, y)g(λ)dλdy.
IR(3.14)

The combination of (3.10), (3.11), and (3.13) yields the system of equations for Ẑk , k ≥ 1:∂
KẐ1 = iω Ẑ1 + P[Ẑ1 ],∂t
2∂Ẑk = ikω Ẑk , k > 1.
∂t

3.3

(3.15)(3.16)

Spectral analysis

In this section, we study (3.15), which decides the linear stability of the mixed state. We rewrite (3.15) as∂
Ẑ1 (t, ω, x) = T [Ẑ1 ](t, ω, x).
∂twhere
T [Z] = iωZ +

K
P[Z].
2

(3.17)(3.18)

Equations (3.14) and (3.18) define linear operators P and T on the weighted Lebesgue space L2 (X, gdωdx)
with X := R × I.
Lemma 3.1. T : L2 (X, gdωdx) → L2 (X, gdωdx) is a closed operator. The residual spectrum of T isempty and the continuous spectrum σc (T ) = i supp(g).
Proof. Consider the multiplication operator Miω : L2 (X, gdωdx) → L2 (X, gdωdx) defined byMiω z = iωz,ω ∈ R.

(3.19)

It is well known that Miω is closed and the (continuous) spectrum of Miω lies on the imaginary axis
σc (Miω ) = i · supp(g). Since W (x, y) is square-integrable by the assumption, P is a Hilbert-Schmidt
operator on L2 (X, gdωdx) and, therefore, is compact [46]. Then, the statement of the lemma follows fromthe perturbation theory for linear operators [17].
Similarly, the spectrum of the operator Mijω lies on the imaginary axis; σ(Mijω ) = ij ·supp(g). Hence,
the trivial solution Ẑj ≡ 0 of (3.16) for j = 2, 3, . . . is neutrally stable.
13

We define a Fredholm integral operator W on L2 (I) by∫
W[V ](x) =W (x, y)V (y)dy.
R(3.20)

Since W is compact, the set of eigenvalues σp (W) is a bounded countable with the only accumulation pointat the origin. Since W is symmetric, all eigenvalues are real numbers.
Lemma 3.2. The eigenvalues of T are given by{}2
σp (T ) = λ ∈ C\σc (T ) : D(λ) =, ζ ∈ σp (W)\{0} ,ζK(3.21)∫

where
D(λ) :=

1g(ω)dω.
R λ − iω(3.22)

Proof. Suppose v ∈ L2 (X, gdωdx) is an eigenvector of T corresponding to the eigenvalue λ:T [v] = λv.
Thenv = 2−1 K(λ − iω)−1 P[v].

By multiplying both sides of (3.23) by g(ω) and integrating with respect to ω, we have
∫ ∫
∫
∫
K
1
g(ω)dω ·W (x, y)v(ω, y)g(ω)dωdy.
v(ω, x)g(ω)dω =2 R λ − iωI RR

(3.23)(3.24)

Rewrite (3.24) asV =K
D(λ)W[V ],2(3.25)

∫

where
V :=

R

v(ω, ·)g(ω)dω ∈ L2 (I).

Equations (3.25) and (3.20) reduce the eigenvalue problem for T to that for the Fredholm operator W.
Suppose V ∈ L2 (I) is an eigenfunction of W associated with the eigenvalue ζ ̸= 0. Then (3.25) implies
D(λ) = 2/(ζK).
If 0 ∈ σp (W) and V is a corresponding eigenvector, then Equation (3.23) yields∫ ∫K 1
W (x, y)v(ω, y)g(ω)dωdyv =2 λ − iω I R
∫
K 1
K 1
=
W (x, y)V (y)dy =(W[V ])(x) = 0.
2 λ − iω I
2 λ − iωThus, ζ = 0 is not an eigenvalue of T .

14

For ζ ̸= 0 denote
K(ζ) =

2 1.
πg(0) |ζ|(3.26)

Lemma 3.3. For each ζ ∈ σp (W) and K > K(ζ) there is a unique eigenvalue of T λ = λ(ζ, K) satisfying
D(λ) = 2/(ζK).
∩
For ζ ∈ σp (W) R+ , λ(ζ, K) is a positive increasing function of K satisfyinglim

λ(ζ, K) = 0 + 0,

K→K(ζ)+0For ζ ∈ σp (W)∩lim λ(ζ, K) = ∞.

K→∞

(3.27)

R− , λ(ζ, K) is a negative decreasing function of K satisfyinglim
K→K(ζ)+0λ(ζ, K) = 0 − 0,lim λ(ζ, K) = −∞.

K→∞(3.28)

Finally,σp (T ) = {λ(ζ, K) :ζ ∈ σp (W)\{0}, K > K(ζ)} ⊂ R\{0}.

Proof. Since ζ ∈ R, setting λ = x + iy for the equation D(λ) = 2/(ζK) yields
 ∫
x
2

,
g(ω)dω =
22
ζK∫R x + (ω − y)
ω−yg(ω)dω = 0.
22
R x + (ω − y)(3.29)

With ζ = 1 this system of equations was analyzed in [5] in the context of the classical Kuramoto model(1.1).
The second equation of (3.29) gives
∫
∫ ∞
ω−y
ω
0=
g(ω)dω =(g(y + ω) − g(y − ω))dω.
2222
R x + (ω − y)0 x +ω

(3.30)

Since g is even, y = 0 is a solution of (3.30). Furthermore, since g is unimodal, there are no other solutionsof (3.30). Thus, λ is real.
With y = 0 the first equation of (3.29) yields∫
2x
g(ω)dω =.
2
2
ζK
Rx + ω

(3.31)

Since g is nonnegative and K > 0, from (3.31) we have ζx > 0. Further, the left-hand side of (3.31) satisfies∫xlim
g(ω)dω = ±πg(0),(3.32)
x→±0 R x2 + ω 2∫xlim
g(ω)dω = 0.
(3.33)
2x→±∞ R x + ω 2
15

The first identity follows from the Poisson’s integral formula for the upper half-plane [37]. The combinationof (3.31) and (3.32) implies that x → 0 as K → K(ζ) + 0, while that of (3.31) and (3.33) yields x → ±∞as K → ∞.
For the uniqueness, it is sufficient to show that the function∫x
x 7→g(ω)dω
2 + ω2xR
is monotonically decreasing in x except at the jump point x = 0. If this were not true, there would exist two
eigenvalues for some interval K1 < K < K2 , and two eigenvalues would collide and disappear at K = K1
or K = K2 . This is impossible, because D(λ) is holomorphic in λ.
To formulate the main result of this section, we will need the following notation:{∩∩max{ζ : ζ ∈ σp (W) R+ }, σp (W) R+ ̸= ∅,ξmax (W) =0 + 0, otherwise.
{∩∩min{ζ : ζ ∈ σp (W) R− }, σp (W) R− ̸= ∅,ξmin (W) =0 − 0, otherwise.
Further, let
Kc+ =2
πg(0)ξmax (W)and Kc− =2.
πg(0)ξmin (W)(3.34)(3.35)

Theorem 3.4. The spectrum of T consists of the continuous spectrum on the imaginary axis and possibly
one or more negative eigenvalues, if K ∈ [Kc− , Kc+ ], and there is at least one positive eigenvalue of T ,otherwise.
Theorem 3.4 follows from Lemma 3.3. It shows that the incoherent state is linearly (neutrally) stable forK ∈ [Kc− , Kc+ ] and is unstable otherwise. The critical values Kc± mark the loss of stability the incoherent
state. The detailed analysis of the bifurcations at Kc± will be presented elsewhere.
Remark 3.5. For the classical Kuramoto model on the complete graph, W ≡ 1 and ζmax (W) = 1. Thus, werecover the well-known Kuramoto’s transition formula (1.3) [20]. In the general case, the transition points
depend on the graph structure through the extreme eigenvalues of the kernel operator W.
Remark 3.6. For nonnegative graphons W , ζmax (W) coincides with the spectral radius of the limiting kernel
operator W (cf. [42]):
ϱ(W) = max{|ζ| : ζ ∈ σp (W)},This can be seen from the variational characterization of the eigenvalues of a self-adjoint compact operator,
which also implies
ϱ(W) = max (W[f ], f ).
(3.36)∥f ∥L2 (I) =1

4

Approximation

Equation (2.5) may be viewed as a base model. To apply our results to a wider class of deterministic andrandom networks, we will need approximation results, which are collected in this section.
16

4.1

Deterministic networks

Consider the Kuramoto model on the weighted graph Γ̃n = ⟨[n], E(Γ̃n ), W̃ ⟩∑˙θ̃ni = ωi + Kn−1
W̃nij sin(θ̃nj − θ̃ni ),ni ∈ [n],

(4.1)

j=1

where W̃n = (W̃nij ) is a symmetric matrix.
Denote the corresponding empirical measure byµ̃nt (A) = n−1n∑δP̃ni (t) (A),A ∈ B(G),

(4.2)

i=1

where P̃ni (t) = (θ̃ni (t), ωi , ξni ), i ∈ [n].
First, we show that if Wn and W̃n are close, so are the solutions of the IVPs for (2.5) and (4.1) with thesame initial conditions. To measure the proximity of Wn = (Wnij ) ∈ Rn×n and W̃n = (W̃nij ) and thecorresponding solutions θn and θ̃n , we will use the following norms:
v
v
u
u
n
n
∑
∑
u
u
2
2 ,
−2
−1
t
t
∥Wn ∥2,n = n
Wnij , ∥θn ∥1,n = n
θni
(4.3)
i,j=1

i=1

where θn = (θn1 , θn2 , . . . , θnn ) and Wn = (Wnij ).
The following lemma will be used to extend our results for the KM (2.5) to other networks.
Lemma 4.1. Let θn (t) and θ̃n (t) denote solutions of the IVP for (2.5) and (4.1) respectively. Suppose thatthe initial data for these problems coincide
θn (0) = θ̃n (0).

(4.4)

Then for any T > 0 there exists C = C(T ) > 0 such thatmax θn (t) − θ̃n (t)t∈[0,T ]1,n≤ C Wn − W̃n

2,n

,

(4.5)

where the positive constant C is independent from n.
Corollary 4.2.
sup d(µnt , µ̃nt ) ≤ C Wn − W̃nt∈[0,T ]

2,n

.

(4.6)

Proof. (Lemma 4.1) Denote ϕni = θni − θ̃ni . By subtracting (4.1) from (2.5), multiplying the result by17n−1 ϕni , and summing over i ∈ [n], we obtain(2K)−1d∥ϕn ∥21,n
dt

=

+

n (
∑

n−2
n−2

i,j=1
n
∑

)
Wnij − W̃nij sin(θnj − θni )ϕni()
W̃nij sin(θnj − θni ) − sin(θ̃nj − θ̃ni ) ϕin

i,j=1=: I1 + I2 .

(4.7)

Using an obvious bound |sin(θnj − θni )| ≤ 1 and an elementary inequality |ab| ≤ 2−1 (a2 + b2 ), we
obtain()|I1 | ≤ 2−1 ∥Wn − W̃n ∥22,n + ∥ϕn ∥21,n .
(4.8)
Further, from Lipschitz continuity of sin and the definition of ϕni , we havesin(θni − θnj ) − sin(θ̃ni − θ̃nj ) ≤ |ϕni − ϕnj | ≤ |ϕni | + |ϕnj |.
Therefore,|I2 | ≤ 2∥W̃ ∥L∞ (I 2 ) ∥ϕn ∥21,n .
The combination of (4.7), (4.8), and (4.9) yields()
dK −1 ∥ϕn (t)∥21,n ≤ 4∥W̃ ∥L∞ (I 2 ) + 1 ∥ϕn (t)∥21,n + ∥Wn − W̃n ∥22,n .
dt
By the Gronwall’s inequality [11, Appendix B.2.j],
(
)
∫ t
2
KC1 t
2
2
∥ϕn (t)∥1,n ≤ e∥ϕn (0)∥1,n +K∥Wn − W̃n ∥2,n ds
0≤ eKC1 t Kt∥Wn − W̃n ∥22,n ,()
where we used ϕn (0) = 0 and C1 := 4∥W̃ ∥L∞ (I 2 ) + 1 . Thus,
sup ∥ϕn (t)∥21,n ≤ eKC1 T KT ∥Wn − W̃n ∥22,n .
t∈[0,T ]Proof. (Corollary 4.2) Let f ∈ L (cf. (2.16)) and consider∫n∑f (θni (t)) − f (θ̃ni (t))
f (dµnt − dµ̃nt ) = n−1
G

≤ n−1

i=1
n
∑

θni (t) − θ̃ni (t)i=1≤ ∥θn (t) − θ̃n (t)∥1,n .
By Lemma 4.1,∫
max d(µnt , µ̃nt ) = supt∈[0,T ]f ∈Lf (dµnt − dµ̃nt ) ≤ C∥Wn − W̃n ∥2,n .
G

18

(4.9)

(4.10)

4.2

Random networks

We now turn to the KM on random graphs. To this end, we use W-random graph Γ̄n = Gr (Xn , W ), which
we define next. As before, Xn is a set of points (2.1),(2.2). Γ̄n is a graph on n nodes, i.e., V (Γ̄n ) = [n]. Theedge set is defined as follows:
P ({i, j} ∈ E(Γn )) = W (ξni , ξnj ).
(4.11)The decision for each pair {i, j} is made independently from the decisions on other pairs.
The KM on the W-random graph Γ̄n = Gr (Xn , W ) has the following form:θ̄˙ni = ωi + Kn−1n∑enij sin(θ̄nj − θ̄ni ),i ∈ [n] := {1, 2, . . . , n}(4.12)

j=1

where enij , 1 ≤ i ≤ j ≤ n are independent Bernoulli RVs:
P(enij = 1) = W (ξni , ξnj ),and enij = enji .
Lemma 4.3. Let θn (t) and θ̄n (t) denote solutions of the IVP for (2.5) and (4.12) respectively. Suppose thatthe initial data for these problems coincide
θn (0) = θ̄n (0).

(4.13)

Then
lim supn→∞ t∈[0,T ]θn (t) − θ̄n (t) 1,n = 0

a.s..

(4.14)

Proof. The proof follows the lines of the proof of Lemma 4.1. As before, we set up the equation for
ϕni := θ̄ni − θni :
−1 d

(2K)

dt

∥ϕn ∥21,n

=

+

n

−2

n−2

n
∑
i,j=1
n
∑

(enij − Wnij ) sin(θnj − θni )ϕni()
enij sin(θ̄nj − θ̄ni ) − sin(θnj − θni ) ϕin

i,j=1=: I1 + I2 .

(4.15)As in (4.9), we have bound|I2 | ≤ n−2n∑(|ϕni | + |ϕnj |) |ϕni | ≤ 2∥ϕn ∥21,n ,i,j=1

where we used 0 ≤ enij ≤ 1.

19

(4.16)

Next, we turn to the first term on the right hand side of (4.15). For this, we will need the following
definitions:
Zni (t) = n

−1

n
∑

anij (t)ηnij ,j=1

anij (t) = sin (θnj (t) − θni (t)) ,= enij − Wnij .

ηnijand Zn = (Zn1 , Zn2 , . . . , Znn ). With these definitions in hand, we estimate I1 as follows:|I1 | = |n−1

n
∑

(
)
Zni ϕni | ≤ 2−1 ∥Zn ∥21,n + ∥ϕn ∥21,n ,(4.17)

i=1

The combination of (4.15)-(4.17) yields,d∥ϕn (t)∥21,n ≤ 5K∥ϕn (t)∥21,n + K∥Zn (t)∥21,n .
dt(4.18)

Using the Gronwall’s inequality and (4.13), we have∥ϕn (t)∥21,n ≤ K exp{5Kt}∫ t∥Zn (s)∥21,n ds

0andsup ∥ϕn (t)∥21,n ≤ K exp{5KT }∫ T(4.19)

0t∈[0,T ]Our next goal is to estimate∥Zn (t)∥21,n dt.

∫T0∥Zn (t)∥21,n dt. Below, we show that

∫ T

∥Zn (t)∥21,n dt → 0

as n → ∞ a.s..

(4.20)

0

To this end, we will use the following observations. Note that ηnik and ηnil are independent for k ̸= l andE ηnij = E enij − Wnij = 0,

(4.21)

where we used P(enij = 1) = Wnij .
Further,22
E ηnij= E(enij − Wnij )2 = E e2nij − Wnij
2
= Wnij − Wnij
≤ 1/4.

(4.22)

and4
E(ηnij) = E(enij − Wnij )4
(
)
2
3
4
= E e4nij − 4e3nij Wnij + 6e2nij Wnij
− 4enij Wnij
+ Wnij2
3
4= Wnij − 4Wnij+ 6Wnij− 3Wnij
4= Wnij (1 − Wnij )4 + Wnij
(1 − Wnij ) ≤ 2−4 .

20(4.23)

Next,∫ T
2

Zni (t) dt = n

n
∑

−2

0

cnikl ηnik ηnil ,(4.24)

k,l=1

where

∫ T
cnikl =

anik (t)anil (t)dt and |cnikl | ≤ T.

(4.25)

0

Further,∫ T∥Zn (t)∥21,n dt = n−30n∑cnikl ηnik ηnil

(4.26)

i,k,l=1

and, finally,
E

(∫ T

)2
∥Zn (t)∥21,n dt = n−60n∑cnikl cnjpq E (ηnik ηnil ηnjp ηnjq ) .

(4.27)

i,k,l,j,p,q=1

We have six summation indices i, k, l, j, p, q ranging from 1 to n. Since E ηnik = 0 for i, k ∈ [n], and RVs
ηnik and ηnjp are independent whenever {i, k} ̸= {j, p}, the nonzero terms on the right-hand side of (4.27)
fall into two groups:
4c2nikk ηnik2 η2
2 η 2 (k ̸= l).
cnikk cnjpp ηnik
c2nikl ηnik
njp (i ̸= j) or
nil

I :
II :

There are n2 terms of type I and 3n3 (n − 1) terms of type II. Thus,
E

(∫ T

)2
∥Zn (t)∥21,n dt()≤ T 2 n−6 n2 + 3n3 (n − 1) = O(n−2 ),(4.28)

0where we used (4.27), (4.22), (4.23), and the bound on |cnikl | in (4.25).
Next, for a given ϵ > 0 we denote
{∫ T
Aϵn =

∥Zn (t)∥21,n dt}≥ϵ .

(4.29)0and use Markov’s inequality and (4.28) to obtain∞∑n=1

P(Aϵn ) ≤ ϵ−2

∞
∑
n=1

E

(∫ T

)2
∥Zn (t)∥21,n dt< ∞.

(4.30)

0

By the Borel-Cantelli Lemma, (4.29) and (4.30) imply (4.20). The latter combined with (4.19) concludes
the proof of Lemma 4.3.

21

5

Examples5.1

Erdős-Rényi graphs

Let p ∈ (0, 1), Xn be defined in (2.1), (2.2), and Wp (x, y) ≡ p. Then Γn,p = Gr (Xn , Wp ) is a family of
Erdős-Rényi random graphs. To apply the transition point formula (3.35), we need to compute the largest
eigenvalue of the self-adjoint compact operator Wp : L2 (I) → L2 (I) defined by∫∫
Wp [f ](x) = W (x, y)f (y)dy = p f (y)dy, f ∈ L2 (I).
(5.1)
I

I

Lemma 5.1. The largest eigenvalue of Wp is ζmax (Wp ) = p.
Proof. Suppose λ ∈ R\{0} is an eigenvalue of Wp and v ∈ L2 (I) is the corresponding eigenvector. Then∫p v(y)dy = ζv(x).
(5.2)
ISince the right-hand side is not identically 0, v = constant ̸= 0. By integrating both sides of (5.2), we findthat ζ = p.
Thus, for the KM on Erdős-Rényi graphs, we have
Kc+ =

5.2

2
,
πg(0)pKc− = −∞.

Small-world network

Small-world (SW) graphs interpolate between regular nearest-neighbor graphs and completely random
Erdős-Rényi graphs. They found widespread applications, because they combine features of regular symmetric graphs and random graphs, just as seen in many real-world networks [45].
Following [28, 29], we construct SW graphs as W-random graphs [24]. To this end let Xn be a set ofpoints from (2.1) satisfying (2.2), and define Wp,r : I 2 → I by{
1 − p, dS (2πx, 2πy) ≤ 2πr,
Wp,r (x, y) =(5.3)
p,otherwise,
where p, r ∈ (0, 1/2) are two parameters.
Definition 5.2. [28] Γn = Gr (Xn , Wp,r ) is called W-small-world graph.
The justification of the mean field limit in Section 2 relies on the assumption that the graphon W is
Lipschitz continuous. To apply Theorem 2.2 to the model at hand, we approximate the piecewise constantϵ :graphon Wp,r by a Lipschitz continuous Wp,rϵ
∥Wp,r − Wp,r∥L2 (I 2 ) < ϵ.

22

(5.4)

ϵ ). The approximation
Further, we approximate the SW graph Γn = Gr (Xn , Wp,r ) by Γϵn = Gr (Xn , Wp,rresults in Section 4 guarantee that the empirical measures generated by the solutions of the IVPs for theKMs on Γn and Γϵn (with the same initial data) are O(ϵ) close with high probability for sufficiently large n.
Thus, below we derive the transition point formula for the KM on (Γϵn ) and take the limit as ϵ → 0 to obtainthe critical values for the KM on SW graphs (Γn ).

Lemma 4.3 and Theorem 2.2 justify (2.10), (2.11) as the continuum limit for the KM on the sequenceof SW graphs (Γϵn ). Thus, (3.35) yields the transition point formula for the KM on SW graphs. To use thisformula, we need to compute the extreme eigenvalues of Wϵp,r : L2 (I) → L2 (I) defined by∫
ϵWϵp,r [f ] = Wp,r(·, y)f (y)dy, f ∈ L2 (I).
I

Lemma 5.3. The largest eigenvalue of Wϵp,r isζmax (Wϵp,r ) = 2r + p − 4pr + oϵ (1).

(5.5)

Proof. We calculate the largest eigenvalue of Wp,r . Using the definition (5.3), one can write∫
Wp,r [f ](x) = Kp,r (x − y)f (y)dy,S

where Kp,r is a 1−periodic function on R, whose restriction to the interval [−1/2, 1/2) is defined as follows{
1 − p, |x| ≤ r,
Kp,r (x) =(5.6)
p,otherwise,
The eigenvalue problem for Wp,r can be rewritten as
Kp,r ∗ v = ζv.
Thus, the eigenvalues of Wp,r are given by the Fourier coefficients of Kp,r (x) as∫
ζk = (K̂p,r )k := Kp,r (x)e−2πikx dx,
S

for k ∈ Z. The corresponding eigenvectors vk = ei2πkx , k ∈ Z, form a complete orthonormal set in L2 (I).
A straightforward calculation yields{
2r + p − 4rp,k = 0,
ζk =(πk)−1 (1 − 2p) sin(2πkr), k =
̸ 0.

(5.7)

Further,ζmax (Wp,r ) = 2r + p − 4rp.

(5.8)

The estimate (5.5) follows from (5.8) and (5.4) via continuous dependence of the eigenvalues with respect
to the parameter ε.
23Thus, for the KM on SW graphs, we obtain
Kc+ =

2
1
.
πg(0) 2r + p − 4pr

The smallest negative eigenvalue is given by ζk(r) for some k(r) ̸= 0. While it is difficult to obtain an
explicit expression k(r), from (5.7) we can find a lower bound for it:ξmin (W) >−1 + 2p
.
π

(5.9)

(
]
In particular, the incoherent state is stable for K ∈ 2 (g(0)(−1 + 2p))−1 , 2 (g(0)π(2r + p − 4pr))−1 .

5.3

Coupled oscillators on a ring

A common in applications type of network connectivity can be described as follows. Consider n oscillatorsplaced uniformly around a circle. They are labelled by integers from 1 to n in the counterclockwise direction.
To each potential edge {i, j} ∈ [n]2 we assign a weight
Wnij = G(ξni − ξnj ),(5.10)

where G is a 1-periodic even measurable bounded function.
Example 5.4. Let Xn = {1/n, 2/n, . . . , 1} and the restriction of the 1-periodic even function Gr on [0, 1/2]
is defined by{
1, 0 ≤ x ≤ r,
Gr (x) =(5.11)
0, x > r,where r ∈ (0, 1/2) is a parameter. With this choice of G := Gr , we obtain a k-nearest-neighbor model, in
which each node is connected to k = ⌊rn⌋ from each side.
Example 5.5. Another representative example was used by Kuramoto and Battogtokh [19]. Here, let Xn ={1/n, 2/n, . . . , 1} and the restriction of the 1-periodic even function G to [0, 1/2] is defined by G(x) :=e−κx , where κ > 0 is a parameter. With this choice of G, we obtain the KM where the strength of interactions
decreases exponentially with the distance between oscillators.
As in our treatment of the KM on SW graphs in §5.2, the integral operator W : L2 (I) → L2 (I) can be
written as a convolution∫
W[f ](x) = G(x − y)f (y)dy, f ∈ L2 (I).
(5.12)
I

It is easy to verify that the eigenvalues of W are given by the Fourier coefficients of G(x).
For instance, for the k-nearest-neighbor network in Example 5.4, by setting p = 0 in (5.8), for the
network at hand we obtain ζmax (W) = r and, thus,
Kc+ =

1
.
πg(0)r24

Note that, in accord with our physical intuition, the transition point is inversely proportional to the coupling
range.
For the network in Example 5.5, eigenvalues are given by2κ(1 − e−κ/2 ), (k : even),
2
2 k2
κ
+
4π
ζk :=
2κ


(1 + e−κ/2 ), (k : odd), k = 0, 1, 2, . . . .
2
κ + 4π 2 k 2(5.13)The largest positive eigenvalue is ζmax (W) = ζ0 , and we obtain
Kc+ =

κ
1
,
πg(0) 1 − e−κ/2

Kc− = −∞.

This recovers the classical result (1.3) as κ → 0.
If the explicit expression for the largest positive eigenvalue of W is not available, for a network with
nonnegative graphon W , the transition point can be estimated using the variational characterization of thelargest eigenvalue (3.36). Specifically, from (3.36), we have:∫ζmax (W) ≥W dxdy = ∥W ∥L1 (I 2 ) ,I2and, thus,
Kc+ ≤

6

2
.
πg(0)∥W ∥L1 (I 2 )Discussion

In this work, we derived and rigorously justified the mean field equation for the KM on convergent familiesof graphs. Our theory covers a large class of coupled systems. In particular, it clarifies the mathematicalmeaning of the mean field equation used in the analysis of chimera states (see, e.g., [33]). Moreover, weshow how to write the mean field equation for the KM on many common in applications random graphs
including Erdős-Rényi and small-world graphs, for which it has not been known before.
We used the mean field equation to study synchronization in the KM on large deterministic and randomgraphs. We derived the transition point formulas for the critical values of the coupling strength, at which
the incoherent state looses stability. The transition point formulas show explicit dependence of the stability
boundaries of the incoherent state on the spectral properties of the limiting graphon. This reveals the precise
mechanism by which the network topology affects the stability of the incoherent state and the onset ofsynchronization. In the follow-up work, we will show that the linear stability analysis of this paper canbe extended to show nonlinear stability of the incoherent state albeit with respect to the weak topology.
There we will also present the bifurcation analysis for the critical values Kc± . The analysis of the KM
on small-world networks shows that, unlike in the original KM (1.1), on graphs the incoherent state may
remain stable even for negative values of K, i.e., for repulsive coupling. In fact, the bifurcations at theleft and right endpoints can be qualitatively different. In the small-world case, the center manifold at Kc−
25

is two-dimensional, whereas it is one-dimensional at Kc+ . These first findings indicate that the bifurcation
structure of the KM on graphs (2.5) is richer than that of its classical counterpart (1.1) and motivates further
investigations of this interesting problem. In the future, we also plan to extend our analysis to the KM on
certain sparse graphs, including sparse power law networks considered in [16].
Acknowledgements. This work was supported in part by the NSF DMS 1412066 (GM).

26

References[1] V.S. Afraimovich, N.N Verichev, and M.I Rabinovich, Stochastic synchronization of oscillations indissipative system, Izy. Vyssh. Uchebn. Zaved. Radiofiz. 29 (1986).
[2] Patrick Billingsley, Probability and measure, third ed., Wiley Series in Probability and Mathematical
Statistics, John Wiley & Sons, Inc., New York, 1995, A Wiley-Interscience Publication. MR 1324786[3] I. I. Blekhman, Sinkhronizatsiya v prirode i tekhnike, “Nauka”, Moscow, 1971. MR 628458[4] W. Braun and K. Hepp, The Vlasov dynamics and its fluctuations in the 1/N limit of interactingclassical particles, Comm. Math. Phys. 56 (1977), no. 2, 101–113. MR 0475547[5] Hayato Chiba, A proof of the Kuramoto conjecture for a bifurcation structure of the infinitedimensional Kuramoto model, Ergodic Theory Dynam. Systems 35 (2015), no. 3, 762–834. MR
3334903[6], A spectral theory of linear operators on rigged Hilbert spaces under analyticity conditions,Adv. Math. 273 (2015), 324–379. MR 3311766[7] Hayato Chiba and Isao Nishikawa, Center manifold reduction for large populations of globally coupledphase oscillators, Chaos 21 (2011), no. 4, 043103, 10. MR 2918832[8] R. L. Dobrušin, Vlasov equations, Funktsional. Anal. i Prilozhen. 13 (1979), no. 2, 48–58, 96. MR
541637[9] F. Dorfler and F. Bullo, Synchronization and transient stability in power networks and non-uniform
Kuramoto oscillators, SICON 50 (2012), no. 3, 1616–1642.
[10] R. M. Dudley, Real analysis and probability, Cambridge Studies in Advanced Mathematics, vol. 74,
Cambridge University Press, Cambridge, 2002, Revised reprint of the 1989 original. MR 1932358[11] Lawrence C. Evans, Partial differential equations, second ed., Graduate Studies in Mathematics,
vol. 19, American Mathematical Society, Providence, RI, 2010. MR 2597943[12] I. M. Gel′ fand and N. Ya. Vilenkin, Generalized functions. Vol. 4, AMS Chelsea Publishing, Providence, RI, 2016, Applications of harmonic analysis, Translated from the 1961 Russian original [MR0146653] by Amiel Feinstein, Reprint of the 1964 English translation [ MR0173945]. MR 3467631[13] François Golse, On the dynamics of large particle systems in the mean field limit, Macroscopic andlarge scale phenomena: coarse graining, mean field limits and ergodicity, Lect. Notes Appl. Math.
Mech., vol. 3, Springer, [Cham], 2016, pp. 1–144. MR 3468297[14] Julien M. Hendrickx and Alex Olshevsky, On Symmetric Continuum Opinion Dynamics, SIAM J.
Control Optim. 54 (2016), no. 5, 2893–2918. MR 3562360[15] Frank C. Hoppensteadt and Eugene M. Izhikevich, Weakly connected neural networks, Applied Mathematical Sciences, vol. 126, Springer-Verlag, New York, 1997. MR 1458890[16] D. Kaliuzhnyi-Verbovetskyi and G. S. Medvedev, The semilinear heat equation on sparse randomgraphs, ArXiv e-prints (2016).
27

[17] Tosio Kato, Perturbation theory for linear operators, Classics in Mathematics, Springer-Verlag, Berlin,
1995, Reprint of the 1980 edition. MR 1335452[18] S. Yu. Kourtchatov, V. V. Likhanskii, A. P. Napartovich, F. T. Arecchi, and A. Lapucci, Theory of phaselocking of globally coupled laser arrays, Phys. Rev. A 52 (1995), 4089–4094.
[19] Y. Kuramoto and D. Battogtokh, Coexistence of coherence and incoherence in nonlocally coupledphase oscillators, Nonlinear Phenomena in Complex Systems 5 (2002), 380–385.
[20] Yoshiki Kuramoto, Self-entrainment of a population of coupled non-linear oscillators, International
Symposium on Mathematical Problems in Theoretical Physics (Kyoto Univ., Kyoto, 1975), Springer,
Berlin, 1975, pp. 420–422. Lecture Notes in Phys., 39. MR 0676492[21] C. Lancellotti, On the Vlasov limit for systems of nonlinearly coupled oscillators without noise, Transport Theory Statist. Phys. 34 (2005), no. 7, 523–535. MR 2265477 (2007f:82098)[22] Ron Levy, William D. Hutchison, Andres M. Lozano, and Jonathan O. Dostrovsky, High-frequency
synchronization of neuronal activity in the subthalamic nucleus of parkinsonian patients with limb
tremor, Journal of Neuroscience 20 (2000), no. 20, 7766–7775.
[23] L. Lovász, Large networks and graph limits, AMS, Providence, RI, 2012.
[24] L. Lovász and B. Szegedy, Limits of dense graph sequences, J. Combin. Theory Ser. B 96 (2006),
no. 6, 933–957. MR 2274085 (2007m:05132)[25] I. G. Malkin, Metody Lyapunova i Puankare v teorii nelineı̆nyh kolebaniı̆, OGIZ, Moscow-Leningrad,],1949. MR 0035895[26], Some problems of the theory of nonlinear oscillations, Gosudarstv. Izdat. Tehn.-Teor. Lit.,
Moscow, 1956. MR 0081402[27] Georgi S. Medvedev, The nonlinear heat equation on dense graphs and graph limits, SIAM J. Math.
Anal. 46 (2014), no. 4, 2743–2766. MR 3238494
[28]
[29]

, The nonlinear heat equation on W-random graphs, Arch. Ration. Mech. Anal. 212 (2014),
no. 3, 781–803. MR 3187677, Small-world networks of Kuramoto oscillators, Phys. D 266 (2014), 13–22. MR 3129708[30] Sebastien Motsch and Eitan Tadmor, Heterophilious dynamics enhances consensus, SIAM Rev. 56(2014), no. 4, 577–621. MR 3274797[31] H. Neunzert, Mathematical investigations on particle - in - cell methods, vol. 9, 1978, pp. 229–254.
[32], An introduction to the nonlinear Boltzmann-Vlasov equation, Kinetic theories and the Boltzmann equation (Montecatini, 1981), Lecture Notes in Math., vol. 1048, Springer, Berlin, 1984, pp. 60–
110. MR 740721 (87i:82061)[33] O.E. Omelchenko, Coherence-incoherence patterns in a ring of non-locally coupled phase oscillators,
Nonlinearity 26 (2013), no. 9, 2469.

28

[34] Mason A. Porter and James P. Gleeson, Dynamical systems on networks, Frontiers in Applied Dynamical Systems: Reviews and Tutorials, vol. 4, Springer, Cham, 2016, A tutorial. MR 3468887[35] Wei Ren, Randal W. Beard, and Timothy W. McLain, Coordination variables and consensus building
in multiple vehicle systems, Cooperative control, Lecture Notes in Control and Inform. Sci., vol. 309,
Springer, Berlin, 2005, pp. 171–188. MR 2106896[36] W. Singer, Synchronization of cortical activity and its putative role in information processing andlearning, Annual Review of Physiology 55 (1993), no. 1, 349–374.
[37] Elias M. Stein and Guido Weiss, Introduction to Fourier analysis on Euclidean spaces, Princeton
University Press, Princeton, N.J., 1971, Princeton Mathematical Series, No. 32. MR 0304972[38] Steven Strogatz, Sync, Hyperion Books, New York, 2003, How order emerges from chaos in the universe, nature, and daily life. MR 2394754[39] Steven H. Strogatz, From Kuramoto to Crawford: exploring the onset of synchronization in populationsof coupled oscillators, Phys. D 143 (2000), no. 1-4, 1–20, Bifurcations, patterns and symmetry. MR
1783382[40] Steven H. Strogatz and Renato E. Mirollo, Stability of incoherence in a population of coupled oscillators, J. Statist. Phys. 63 (1991), no. 3-4, 613–635. MR 1115806[41] Steven H. Strogatz, Renato E. Mirollo, and Paul C. Matthews, Coupled nonlinear oscillators below
the synchronization threshold: relaxation by generalized Landau damping, Phys. Rev. Lett. 68 (1992),
no. 18, 2730–2733. MR 1160290[42] Balázs Szegedy, Limits of kernel operators and the spectral regularity lemma, European J. Combin. 32(2011), no. 7, 1156–1167. MR 2825541[43] Roger D. Traub, Miles A. Whittington, and Mark O. Cunningham, Epileptic fast oscillations andsynchrony in vitro, Epilepsia 51 (2010), 28–28.
[44] S. Watanabe and S.H. Strogatz, Constants of motion for superconducting Josephson arrays, Physica
D: Nonlinear Phenomena 74 (1994), no. 34, 197 – 253.
[45] D.J. Watts and S.H. Strogatz, Collective dynamics of small-world networks, Nature 393 (1998), 440–442.
[46] Nicholas Young, An introduction to Hilbert space, Cambridge Mathematical Textbooks, Cambridge
University Press, Cambridge, 1988. MR 949693

29