The mean field analysis of the Kuramoto model on graphs II.
Asymptotic stability of the incoherent state, center manifold
reduction, and bifurcations
Hayato Chiba∗ and Georgi S. Medvedev†
September 25, 2017

Abstract
In our previous work [3], we initiated a mathematical investigation of the onset of synchronizationin the Kuramoto model (KM) of coupled phase oscillators on convergent graph sequences. There, we
derived and rigorously justified the mean field limit for the KM on graphs. Using linear stability analysis,
we identified the critical values of the coupling strength, at which the incoherent state looses stability,
thus, determining the onset of synchronization in this model.
In the present paper, we study the corresponding bifurcation. Specifically, we show that similar tothe original KM with all-to-all coupling, the onset of synchronization in the KM on graphs is realized
via a pitchfork bifurcation. The formula for the stable branch of the bifurcating equilibria involves theprincipal eigenvalue and the corresponding eigenfunctions of the kernel operator defined by the limit of
the graph sequence used in the model. This establishes an explicit link between the network structure
and the onset of synchronization in the KM on graphs. The results of this work are illustrated with thebifurcation analysis of the KM on Erdős-Rényi, small-world, as well as certain weighted graphs on acircle.

1

Introduction

In 1970s, a prominent Japanese physicist Yoshiki Kuramoto described a remarkable effect in collective dynamics of large systems of coupled oscillators [11]. He studied all-to-all coupled phase oscillators withrandomly distributed intrinsic frequencies, the model which now bears his name. When the strength of coupling is small, the phases are distributed approximately uniformly around a unit circle, forming an incoherent
state. Kuramoto identified a critical value of the coupling strength, at which the the incoherent state looses
stability giving rise to a stable partially synchronized state. To describe the bifurcation corresponding to theonset of synchronization, Kuramoto introduced the order parameter - a scalar function, which measures the∗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degree of coherence in the system. He further showed that the order parameter undergoes a pitchfork bifurcation. Thus, the qualitative changes in the statistical behavior of a large system of coupled phase oscillatorsnear the onset to synchronization can be described in the language of the bifurcation theory.
Kuramoto’s discovery created a new area of research in nonlinear science [20, 21, 19]. The rigorous
mathematical treatment of the pitchfork bifurcation in the KM was outlined in [4] and was presented withall details in [1]. The analysis in these papers is based on the generalized spectral theory for linear operators[2] and applies to the KM with intrinsic frequencies drawn from a distribution with analytic or rational
probability density function. Under more general assumptions on the density, the onset of synchronizationin the KM was analyzed in [6, 8, 5]. These papers use analytical methods for partial differential equations
and build upon a recent breakthrough in the analysis of Landau damping [18].
In our previous work [3], we initiated a mathematical study of the onset of synchronization in the KM
on graphs. Following [14, 15], we considered the KM on convergent families of deterministic and randomgraphs, including Erdős-Rényi and small-world graphs among many other graphs that come up in applications. For this model, we derived and rigorously justified the mean field limit. The latter is a partial
differential equation approximating the dynamics of the coupled oscillator model in the continuum limit as
the number of oscillators grows to infinity. In [3], we performed a linear stability analysis of the incoherent
state and identified the boundaries of stability. Importantly, we related the stability region of the incoherent
state to the structural properties of the network through the spectral properties of the kernel operator definedby the limit of the underlying graph sequence [12]. In the present paper, we continue the study initiated in[3]. Here, we analyze the bifurcations at the critical values of the coupling strength, where the incoherent
state looses stability. In the remainder of this section, we review the model and the main outcomes of [3]and explain the main results of this work.
We begin with a brief explanation of the graph model that was used in [3] and will be used in this paper.
In [3], we adapted the construction of W-random graphs [13] to define a convergent sequence of graphs. The
flexible framework of W-random graphs allows us to deal with a broad class of networks that are of interestin applications. Let W be a symmetric measurable function on the unit square I 2 := [0, 1]2 with values in[−1, 1] and let Xn = {ξn1 , ξn2 , . . . , ξnn } form a triangular array of points ξni , i ∈ [n] := {1, 2, . . . n}, n ∈N, subject to the following condition
Z
n
X
−1
lim n
(1.1)
f (ξni ) = f (x)dx ∀f ∈ C(I).
n→∞

I

i=1

Γn = hV (Γn ), E(Γn ), (Wnij )i is a weighted graph on n nodes labeled by the integers from [n], whose edge
set isE(Γn ) = {{i, j} : W (ξni , ξnj ) 6= 0, i, j ∈ [n]} .
Each edge {i, j} ∈ E(Γn ) is supplied with the weight Wnij := W (ξni , ξnj ). In the theory of graph limits,
W is called a graphon [12]. It defines the asymptotic properties of {Γn } for large n.
Consider the initial value problem (IVP) for the KM on Γnθ̇ni = ωi + Kn−1nX

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

(1.2)

j=1

θni (0) =

0
θni
.

(1.3)
2

The intrinsic frequencies ωi , i ∈ [n], are independent identically distributed random variables. The distribution of ω1 has density g(ω). For the spectral analysis in Section 3, we need to impose the following
assumptions on g: a) g : R →RR+ ∪ {0} is an even unimodal function, and b) g is real analytic function withfinite moments of all orders: R |x|m g(x)dx < ∞, m ∈ N. For instance, the density of the Gaussian distribution satisfies these conditions. The KM on weighted graphs {Γn } (1.2), (1.3) can be used to approximate
the KM on a variety of random graphs (cf. § 4.2 [3]).
Along with the discrete model (1.2) we consider the IVP for the following partial differential equation∂∂
ρ(t, θ, ω, x) +{ρ(t, θ, ω, x)V (t, θ, ω, x)} = 0,∂t∂θρ(0, θ, ω, x) = ρ0 (θ, ω, x) ∈ S × R × I,
where

(1.4)
(1.5)

Z Z Z
W (x, y) sin(φ − θ)ρ(t, φ, ω, y)g(ω)dφdωdy.

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

R

(1.6)

S

Here, ρ(t, θ, ω, x) is the conditional density of the random vector (θ, ω) given ω, and parametrized by(t, x) ∈ R+ × I, and S = R/2πZ is a circle. In particular,Zρ(t, θ, ω, x)dθ = 1 ∀(t, ω, x) ∈ R+ × R × I.
(1.7)
S

It is shown in [3, Theorem 2.2] that
µnt (A) = n−1nX

δ(θni (t),ωi ,ξni ) (A)

(1.8)

i=1

interpreted as a probability measure on Borel sets A ∈ B(G), G = S × R × I, converges in the bounded
Lipschitz distance [7] uniformly on bounded time intervals to the absolutely continuous measureZµt (A) =ρ(t, θ, ω, x)g(ω)dθdωdx, A ∈ B(G),(1.9)
A

provided µn0 and µ0 are sufficiently close in the same distance. The latter can be achieved with the appropriate initial condition (1.3) and sufficiently large n (see [3, Corollary 2.3]). Therefore, the IVP (1.4),(1.5)
approximates the IVP (1.2),(1.3) on finite time intervals for sufficiently large n.
An inspection of (1.4) shows that ρu = 1/(2π), the density of the uniform distribution on S, is a steadystate solution of (1.4). It corresponds to the incoherent (mixing) state of the KM. Numerics suggests thatthe incoherent state is stable for small K ≥ 0. The loss of stability of the incoherent state is interpreted as
the onset of synchronization in the KM. This is the main focus of [3] and of the present paper. In [3], we
identified the boundaries of the region of stability of the incoherent state in (1.4). Specifically, we showedthat there exist Kc− ≤ 0 ≤ Kc+ such that ρu is linearly stable for K ∈ [Kc− , Kc+ ], and is unstable otherwise.
The critical values Kc− and Kc+ depend on the network topology through the eigenvalues of the compact
symmetric operator W : L2 (I) → L2 (I)Z
W[f ](x) =W (x, y)f (y)dy.
(1.10)
R

3

The eigenvalues of W are real with the only accumulation point at 0. Denote the largest positive andsmallest negative eigenvalues of W by µmax and µmin respectively. If all eigenvalues are nonnegative(nonpositive), we set µmin = −∞ (µmax = ∞). The main stability result of [3] yields explicit expressions
for the transition points
2
2
Kc− =
and Kc+ =.
(1.11)
πg(0)µminπg(0)µmax
Thus, the region of linear stability of ρu depends explicitly on the spectral properties of the limiting graphon
W . Recall that W represents the graph limit of {Γn }. Thus, (1.11) links network topology to synchronization in (1.2). For the classical KM (all-to-all coupling), W (x, y) = 1 and µmin = −∞, µmax = 1, which
recovers the known Kuramoto’s formula.
In the present paper, we study the onset of synchronization in (1.2) in more detail. After some preliminaries and preparatory work in Sections 2 and 3, we revisit linear stability of the incoherent solution. Thistime, we show that despite the lack of eigenvalues with negative real part and the presence of the continuous
spectrum on the imaginary axis, the incoherent state is an asymptotically stable solution of the linearized
problem (cf. Theorem 4.1). This is a manifestation of the Landau damping in the KM.
In Section 5, we study the bifurcation at Kc+ with a one-dimensional center manifold. To this end, werecall the order parameter
Z Z Zh(t, x) =W (x, y)eiθ ρ(t, θ, ω, y)g(ω)dθdωdy,(1.12)
I

R

S

which was introduced in [3] as a measure of coherence in the KM on graphs. This is a continuous analog of
the local order parametern

1X
Wnij eiθnj (t)n
j=1

for the discrete model (1.2). The order parameter generalizes the original order parameter used by Kuramoto
for the all-to-all coupled model. Note that (1.12) depends on x and contains information about the structureof the network through W . As will be clear below, the order parameter plays an important role in the analysisof the mean field equation. In particular, it can be used to locate nontrivial steady state solutions. To thisend note that the velocity field (1.6) can be conveniently rewritten in terms of the order parameterK  −iθV (t, θ, ω, x) = ω +e h(t, x) + eiθ h(t, x) .
(1.13)
2i
In particular, for a given steady state of the order parameter written in the polar formh∞ (x) = R(x)eiΦ(x) ,(1.14)

the velocity field takes the following form
V (t, θ, ω, x) = ω − KR(x) sin(θ − Φ(x)).
Setting ∂θ (V ρ) = 0, we find the corresponding steady state solution of (1.4)
 δ θ − Φ(x)p − arcsin(ω/KR(x)) , |ω| ≤ KR(x),ρ(θ, ω, x) =
ω 2 − K 2 R(x)21
,|ω| > KR(x),2π |ω − KR(x) sin(θ − Φ(x))|
4(1.15)

where δ stands for the Dirac delta function. The stationary solution (1.15) has the following interpretation:
the first line describes phase-locked oscillators, while the second line yields the distribution of the drifting
oscillators. Thus, solutions of this form may combine phase-locked oscillators and those moving irregularly.
Such solutions are called partially phase-locked or partially synchronized. The phase of a phase-locked
oscillator at x with a natural frequency ω, is given byωθ = Φ(x) + arcsin,(1.16)
KR(x)provided |ω| ≤ KR(x). In Sections 5 and 6, we will identify branches of stable equilibria bifurcating from
Kc± in terms of the corresponding values of the order parameter. Then Equation (1.16) will be used todescribe the corresponding stable phase-locked solutions.
In Section 5, assuming that µmax is a simple eigenvalue of W, we show that the coupled system (1.2) undergoes a supercritical pitchfork bifurcation at Kc+ . Specifically, we derive an ordinary differential equation
for the order parameter h and show that the trivial solution of this equation looses stability at Kc+ and gives
rise to a stable branch of (nontrivial) equilibria, corresponding to partially synchronized state (cf. (5.1)). In
Section 6, we consider the onset of synchronization in networks with certain symmetries (cf. (6.1)). This
leads to the bifurcation with a two-dimensional center manifold. The bifurcation analysis in Sections 5 and
6 is illustrated with the analysis of the KM on Erdős-Rényi, small-world graphs, and to a class of weighted
graphs on a circle.

2

Preliminaries

In the remainder of this paper, we will assume that K ≥ 0. The case of negative K is reduced to thatabove by switching to K := −K and W := −W . Furthermore, without loss of generality we assume that
µmax > 0.

2.1

Fourier transform

We rewrite (1.4) in terms of the complex Fourier coefficientsZzj = eijθ ρ(t, θ, ω, x)dθ,j ∈ Z.

(2.1)

S

Applying the Fourier transform to (1.4) and using integration by parts, we obtainZ∂zj∂= − eijθ{V (t, θ, ω, x)ρ(t, θ, ω, x)} dθ∂t∂θZS= ij eijθ V (t, θ, ω, x)ρ(t, θ, ω, x)dθ
S
Z


jK
= ijωzj +eijθ e−iθ h(t, x) − eiθ h(t, x) ρ(t, θ, ω, x)dθ,
2 S

5

(2.2)

whereZ Z Z

W (x, y)eiθ ρ(t, θ, ω, y)g(ω)dθdωdyh(t, x) =
ZI ZR

S

(2.3)

W (x, y)z1 (t, ω, y)g(ω)dωdy.

=
I

R

From (2.2), we further obtain∂zjjK
hzj−1 − hzj+1 ,= ijωzj +∂t2j ∈ Z.

(2.4)

By (1.7), z0 = 1. Further, z−j = z̄j , because ρ is real. Thus, in (2.4) we can restrict to j ∈ N.
Let
Z
Pf (ω, x) =W[f ](ω, x)g(ω)dω
ZR Z=
W (x, y)f (ω, y)g(ω)dωdy.
I

(2.5)

R

Combining these observations, we rewrite (2.4):
∂K
z1 = iωz1 +Pz1 − Pz1 z2 ,
∂t
2

∂
jK
zj = ijωzj +Pz1 zj−1 − Pz1 zj+1 ,
∂t
2

(2.6)
j = 2, 3, . . . .

(2.7)

Note that the trivial solution Z := (z1 , z2 , . . . ) ≡ 0 is a steady state solution of (2.6), (2.7). It corresponds to the uniform distribution ρu = 1/(2π), a constant steady state solution of (1.4). Linearizing aroundZ ≡ 0, we arrive at∂
z1 = Tz1 ,∂t∂zj = ijωzj ,
∂t

(2.8)
j = 2, 3, . . . ,(2.9)

where T is a linear operator on H := L2 (R × I, g(ω)dωdx)
Tf = iωf +

2.2

K
Pf.
2

(2.10)

ω∈R(2.11)The eigenvalue problem

The multiplication operator Miω : H → H defined byMiω f = iωf,6

is a closed operator. The continuous spectrum of Miω fills the imaginary axis
σc (Miω ) = isupp(g) = iR.

(2.12)

Since P is compact (as a Hilbert-Schmidt operator), T : H → H is closed and σc (T) = iR.
Next we turn to the eigenvalue problem
Tf = λf,(2.13)

where T and P are operators on H (cf. (2.10) and (2.5)).
We will locate the eigenvalues of T through the eigenvalues of W (cf. (1.10)). Since W is a compact
symmetric operator on L2 (I), it has a countable set of real eigenvalues with the only accumulation point atzero. All nonzero eigenvalues have finite multiplicity.
Suppose λ is an eigenvalue of T and v ∈ H is the corresponding eigenfunction. Then a simple calculation yields (cf. [3])K
w = D(λ)Ww,(2.14)
2whereZg(ω)dω
D(λ) =,(2.15)
λ − iω
ZRw =v(ω, ·)g(ω)dω ∈ L2 (I).
(2.16)
R

Equation (2.14) yields the equation for eigenvalues of T
D(λ) =

2,Kµ(2.17)

where µ is a nonzero eigenvalue of W.
Using (2.17), we establish a one-to-one correspondence between the eigenvalues of W and those of T.
Specifically, for every positive eigenvalue of W, µ, there is a branch of eigenvalues of T,
λ = λ(µ, K),

2,πg(0)µ(2.18)

lim λ(µ, K) = ∞.

(2.19)K ≥ K(µ) :=such thatlim

λ(µ, K) = 0+,

K→K(µ)+0

K→∞

Recall that µmax stands for the largest positive eigenvalue of W. Then for K ∈ [0, K(µmax )) there are no
eigenvalues with positive real part. Furthermore, for small ε > 0 and K ∈ (K(µmax ), K(µmax ) + ε) there
is a unique positive eigenvalue of T, λ(K, µmax ), which vanishes as K → K(µmax ) + 0 (see [3] for moredetails).

7

3

The generalized spectral theory

The major obstacle in studying stability and bifurcations of the incoherent state is the continuous spectrumof the linearized problem on the imaginary axis (cf. (2.12)). To deal with this difficulty, we develop thegeneralized spectral theory following the treatment of the classical KM in [1]. Below, we outline the key
steps in the analysis of the generalized eigenvalue problem referring the interested reader to [1] for missing
proofs and further details.

3.1

The rigged Hilbert space

Here, we define the rigged Hilbert space which will be used in the spectral analysis below. Let Exp(β, n)be the set of holomorphic functions on the region Cn := {z ∈ C | Im(z) ≥ −1/n} such that the norm||φ||β,n :=supe−β|z| |φ(z)|(3.1)

Im(z)≥−1/n

is finite. With this norm, Exp(β, n) is a Banach space [1]. Let Exp be their inductive limit with respect ton = 1, 2, . . . and β = 0, 1, 2, . . .
[ [Exp(β, n).
(3.2)
Exp = lim lim Exp(β, n) =−→ −→β≥0 n≥1β≥0 n≥1

With the inductive limit topology, Exp is a complete Montel space 1 . The properties of Exp are describedin detail in [1].
Next, we define X := Exp ⊗L2 (I) as a topological tensor product of Exp and L2 (I). X is a complete
Montel space. It is a dense subspace of H. The topology of X is stronger than that of H. For every f ∈ X
we have f (ω, ·) ∈ L2 (I) for each ω ∈ R. In addition, f (ω, x) is holomorphic in ω on the upper half plane,
where it can grow at most exponentially.
Let X 0 be the dual space of X. It is the space of continuous antilinear functionals on X. Let h·, ·i denotethe pairing between X 0 and X, i.e., for l ∈ X 0 and f ∈ X, hl, f i := l(f ) stands for the corresponding
antilinear functional.
The dual space X 0 is equipped with the weak dual topology (cf. [1])2 . H is a dense subspace of X 0 , andthe topology of H is stronger than that of X 0 . Hence,
X ⊂ H ⊂ X0
form a rigged Hilbert space (a.k.a. Gelfand triple) [10].
Note that if l ∈ X 0 ∩ L2 (R × I, g(ω)dωdx) thenZ Z∗
hl, f i := (l, f )L2 (R×I) =l(ω, x)f (ω, x)g(ω)dωdx,
I R
1
2

A locally convex topological vector space is called Montel if every bounded set in it is relatively compact.
{ln } ⊂ X 0 converges to l ∈ X 0 if hln , f i ∈ C tends to hl, f i for every f ∈ X.

8

for f ∈ X, where f ∗ (ω, x) := f (ω, x). Thus, T and the rigged Hilbert space defined above satisfy all
assumptions of the generalized spectral theory in [1].

3.2

The generalized eigenvalue problem

In this subsection we calculate the resolvent of T and spectral projections. With the rigged Hilbert space
defined above, we will view the resolvent as an operator from X to X 0 .
Below, we will need to construct analytic continuation for certain functions involving integrals of Cauchy
type. For this, we are going to use an implication of the Sokhotski formulas, which we formulate as aseparate statement for convenience.
Lemma 3.1. (Sokhotski, cf. [9]) Let f be a complex valued function on R. Suppose f has at most a finite
number of integrable discontinuities. ThenZf (ω)dωF (z) =
(3.3)
R z − iωis an analytic function in the right and left open half-planes of C. Furthermore, for z = x+iy, the following
formulas determine the limits of F (z) as x → 0±:Zf (−iφ)dφlim F (z) = ±πf (y) + i PV
x→0±
(3.4)
iR φ − iy= ±πf (y) − iπH[f ](y),where PV stands for the principal value in the sense of Cauchy and H[f ] denotes the Hilbert transform off.
Corollary 3.2. Suppose f is holomorphic on the real axis and admits the analytic continuation to the upper
half-plane. Thenx > 0,
 F (z),limx→0+ F (z),x = 0,(3.5)
F̃ (z) =F (z) + 2πf (−iz), x < 0,is an entire function.

3.3

The generalized resolvent

Our next goal is to compute the resolvent of T
R(λ) = (λ − T)−1 .

(3.6)

To this end, we first compute R(λ) for Re(λ) > 0 and extend it analytically to the left half-plane as an
operator from X to X 0 .

9In the right half-plane Re(λ) > 0, R(λ) can be rewritten as follows
−1KR(λ) = A(λ) I − PA(λ)
,
2

(3.7)

A(λ) = (λ − iω)−1 ,

(3.8)

where
and I stands for the identity operator. Note that A(λ) ceases to exist as the multiplication operator on H as
Re(λ) → 0 (recall that the imaginary axis is the continuous spectrum of Miω ). However, it can be extended
to the left half-plane as as an operator A : X → X 0 defined as follows∗Re(λ) > 0,
 (A(λ)u, v )H ,∗limRe(λ)→0+ (A(λ)u, v )H ,Re(λ) = 0, (3.9)hA(λ)u, vi =R(A(λ)u, v ∗ )L2 (R×I,gdωdx) + 2πg(−iλ) I u(−iλ, x)v(−iλ, x)dx, Re(λ) < 0.
By Corollary 3.2, hA(λ)u, vi is an entire function in λ for all u, v ∈ X. This suggests an appropriate
generalization of R(λ), R(λ) : X → X 0 defined byKR(λ) = A(λ) I − P× A(λ)
2

−1
,

(3.10)

where P× : X 0 → X 0 is the dual operator of P.
For each u ∈ X, R(λ)u is an X 0 -valued meromorphic function. For Re(λ) > 0, it coincides with therestriction of R(λ) to X. Thus, R(λ) is a meromorphic continuation of R(λ) from the right half-plane tothe left-half plane as an X 0 -valued operator. Note that since T has the continuous spectrum on the imaginary
axis, R(λ) can not be continued to the left-half plane as an operator on H.
We define the generalized eigenvalues of T as the singularities of the generalized resolvent R(λ).
Definition 3.3. λ ∈ C is called a generalized eigenvalue of T if there is a nonzero v ∈ X 0 such thatK ×
I − P A(λ) v = 0(3.11)
2
In this case, v is called a generalized eigenfunction.
It turns out that the generalized eigenvalues and the corresponding eigenfunctions of T are, in fact, theeigenvalues and eigenfunctions of the dual of T, T× .
Theorem 3.4. (cf. [2]) Let λ ∈ C be a generalized eigenvalue of T and v ∈ X 0 is the corresponding
eigenfunction. Then T× v = λv.
Remark 3.5. Using (3.9) and (3.11), one can see that the generalized eigenvalues λ = λ(µ, K) of T are theroots of the following equation
2= D(λ),(3.12)
Kµ10where µ is a nonzero eigenvalue of W and
Re(λ) > 0, D(λ),
lim D(λ),Re(λ) = 0,D(λ) =
Re(λ)→0+
D(λ) + 2πg(−iλ), Re(λ) < 0.

(3.13)The right hand side of (3.13) is an entire function (cf. Corollary 3.2). For Re(λ) > 0, (3.12) is reduced tothe equation for the eigenvalues of T (cf. (2.17)). In this case, the corresponding generalized eigenfunctionv is included in L2 (R × I, gdωdx), i.e., λ is an eigenvalue of T. On the other hand, for Re(λ) ≤ 0, thegeneralized eigenfunction v is not in H but is an element of the dual space X 0 .

3.4

The generalized Riesz projection

Let µ be a positive eigenvalue of W and w ∈ L2 (I) be the corresponding eigenfunction. The largest positive
eigenvalue of W and the corresponding eigenfunction are denoted by µmax and wmax respectively. For
every K > Kc+ = 2/(πg(0)µmax ) there is a real positive eigenvalue of T, λ = λ(µ, K). The corresponding
eigenfunction is given by
K w(x)v(ω, x) =.
(3.14)
2 λ − iωAs K approaches the critical value Kc+ from above, the eigenvalue λ(µmax , K) converges to 0+ along
the real axis and at K = Kc+ it hits the continuous spectrum on the imaginary axis. The corresponding
eigenfunction approaches the critical vector
X 0 3 vc+ :=

Kc+
wmax
lim*
,
2 λ→0+ λ − iω

(3.15)

where lim* stands for the limit in X 0 with respect to the weak dual topology, i.e., the action of vc+ ∈ X 0 on
u ∈ X is given byZK+
wmax (x)g(ω)
hvc+ , ui = c limu(ω, x)dωdx
2 λ→0+ R×I
λ − iω
(3.16)Z (w∗max , u (ω, ·))L2 (I) g(ω)dω
Kc+
=
lim
.
2 λ→0+ R
λ − iω
Let λ ∈ C be a generalized eigenvalue of T. Then the generalized Riesz projection Πλ : X → X 0 is defined
by
Z
1
Πλ =
R(z)dz,(3.17)
2πi γ(λ)
where γ(λ) is a simple closed curve around λ oriented counterclockwise that does not encircle or intersect
the rest of the spectrum. Below, we shall refer to such curves as contours. The image of Πλ gives thegeneralized eigenspace of λ [1].

11

Theorem 3.6. Suppose the algebraic and geometric multiplicities of µmax coincide. Then the generalized
Riesz projection of λ = 0, the generalized eigenvalue of T for K = Kc+ , has the following form
Π0 = − lim* D0 (λ)−1 A(λ)Π̃µmax D(λ)λ→0+
(3.18)= g1 lim* A(λ)Π̃µmax D(λ) ,λ→0+where g1 = − limλ→0+ D0 (λ)−1 is a positive constant, A(λ) was defined in (3.8), and Π̃µ stands for the
Riesz projection onto the eigenspace of W corresponding to the eigenvalue µ. The operator D(λ) on H isdefined byZv(ω, ·)g(ω)dωD(λ)v =, v ∈ L2 (R × I, gdωdx).
(3.19)
λ − iω
R
The proof of Theorem 3.6 relies on three technical lemmas. Below we state and prove these lemmas first
and then prove the theorem.
Lemma 3.7. Let Re(λ) > 0 thenR(λ)v = A(λ)v +


−1
K
K
A(λ)W I − D(λ)W
D(λ)v,
22v ∈ H.

(3.20)

Proof. By definition of R(λ) (3.6), for any v ∈ H, we haveK
λ − iω − P R(λ)v = v,
2
and, thus,
R(λ)v = A(λ)v +K
A(λ)PR(λ)v.
2(3.21)

Using Fubini theorem, from (2.5) we haveZ(R(λ)v)(ω, ·)g(ω)dωPR(λ)v = W
R(3.22)=: WQv.
On the other hand, integrating both sides of (3.21) against g(ω)dω, we obtain
Z
Qv = (R(λ)v)(ω, ·)g(ω)dωRZKW (·, y) (R(λ)v) (ω, y)g(ω)dωdy= D(λ)v + D(λ)
2
R×IK
= D(λ)v + D(λ)WQv.
2
and
Q=


−1
K
I − D(λ)W
D(λ).
2
12

(3.23)(3.24)Plugging (3.24) into (3.22), we have−1KPR(λ)v = W I − D(λ)W
D(λ)v.
2(3.25)

The combination of (3.21) and (3.25) yields (3.20).
Lemma 3.8. Let λ = λ(µ, K) > 0 be an eigenvalue of T corresponding to the positive eigenvalue of W,µ, and K > Kc+ , and suppose that the geometric and algebraic multiplicities of µ coincide.
ThenΠλ = −D0 (λ)−1 A(λ)Π̃µ D(λ),(3.26)

where Π̃µ is the Riesz projection onto the eigenspace of W corresponding to µ.
Proof. As before, let γ(λ) denote a contour around λ. From (3.20), we haveZ
R(z)dz =
γ(λ)

K
2


−1
K
A(z)W I − D(z)W
D(z)dz.
2γ(λ)

Z(3.27)We change variable in the integral on the right–hand side to ζ = 2(KD(z))−1 . By deforming the contour
γ(λ) if necessary, we can always achieve D0 (z) 6= 0 for z ∈ γ(λ), so that this change of variable ζ = ζ(z)is well defined. Under this transformation, γ(λ) is mapped to γ̃(µ), a contour around µ. Thus, we haveZZdζR(z)dz = −A(z(ζ))W (ζ − W)−1 D(z(ζ)) 0.
(3.28)
ζD (z(ζ))
γ(λ)
γ̃(µ)
Since the algebraic and geometric multiplicities of µ are equal, the singularity of (ζ − W)−1 at ζ = µ is a
simple pole, and the other factor in the integrand of the above is regular at ζ = µ. Therefore, the right–handside of (3.28) simplifies to
!
Z
Z
1
−1
R(z)dz = −A(λ)W(ζ − W) D(z(ζ))dζ .
(3.29)µD0 (λ)γ(λ)γ̃(µ)
By multiplying both sides of (3.29) by (2πi)−1 , we have
Πλ = −

1
A(λ)WΠ̃µ D(λ).
µD0 (λ)Finally, since Π̃µ is the projection on the eigensubspace of W,
WΠ̃µ D(λ) = µΠ̃µ D(λ).
Thus,Πλ = −(D0 (λ))−1 A(λ)Π̃µ D(λ).

13(3.30)

Lemma 3.9.

g(ω)dω
i n π  (n)
(n)
=
g
(0)
−
iH[g
](0)
.
n+1
n!
R (z − iω)

Z
lim

z→0+

(3.31)

Proof. Using integration by parts n times, we obtaing(ω)dω
in=n+1n!
R (z − iω)Zg (n) (ω)dω.
R z − iωZ

The application of Lemma 3.1 to the integral on the right-hand side yields (3.31).
Below will need the following implications of Lemma 3.9.
Corollary 3.10.
lim D0 (z) = −πH[g 0 ](0) < 0,

z→0+

Z
lim

z→0+

g(ω)dω
3
R (z − iω)

=

−π 00
g (0).
2

Proof. Differentiating D(z) and using (3.31), for z off the imaginary axis we haveZ
Z 0g(ω)dω
g (ω)dω0
D (z) = −
= −i
.
2
(z
−
iω)
R
R z − iω(3.32)(3.33)(3.34)The integral on the right–hand side is of Cauchy type and Lemma 3.1 applies. By (3.4),
lim D0 (z) = −iπg 0 (0) − πH[g 0 ](0).

z→0+(3.35)

Since g is even, g 0 (0) = 0 and g 0 is odd. Because g is also nonnegative and unimodal g 0 (x) ≤ 0, x > 0.
Thus,
Z ∞ 0
−1
g (s)ds0
H[g ](0) =
PV
π
s
−∞
(3.36)
Z ∞ 0−2
g (s)ds=lim> 0.
π →0+ s
The combination of (3.4) and (3.36) yields (3.32).
Likewise, (3.33) follows from Lemma 3.9 for n = 2 and Lemma 3.1.
Proof of Theorem 3.6. Theorem 3.6 follows from (3.30) and (3.32).

14

Figure 1: Deformation of the integral path for the Laplace inversion formula.

4

Asymptotic stability of the incoherent state

We now return to the problem of stability of the incoherent state. Recall that in the Fourier space theincoherent state corresponds to the trivial solution Z = (z1 , z2 , · · · ) = 0. The linearization about Z = 0
shows that it is a neutrally stable equilibrium of (2.8), (2.9) for 0 ≤ K < Kc+ . There are no eigenvaluesof T for these values of K and the continuous spectrum fills out the imaginary axis. Nonetheless, we showthat the incoherent state is asymptotically stable with respect to the weak dual topology.
Theorem 4.1. For K ∈ [0, Kc+ ) the trivial solution of (2.8), (2.9) is an asymptotically stable equilibrium
for initial data from X ⊂ X 0 with respect to the weak dual topology on X 0 .
Remark 4.2. The stability with respect to the weak dual topology is weaker than that with respect to the
topology of the Hilbert space H. Still it is a natural topology for the problem at hand. In particular,
Theorem 4.1 implies that the order parameter evaluated on the trajectories of the linearized problem tends
to 0 as t → ∞.
Proof. Integrating (2.9) subject to zj (0, ·) ∈ X, we havezj (t, ·) = eijωt zj (0, ·), j ≥ 2.
By the Riemann-Lebesgue lemma,
Z Z
hzj (t, ·), ψ(·)i =eijωt zj (0, ω, x)ψ(ω, x)dωdx → 0, as t → ∞,I∀ψ ∈ X.

(4.1)

R

We now turn to (2.8). Consider the semigroup generated by the operator T on H given by the Laplace
inversion formula
Z a+ib1tT
e = limeλt (λ − T)−1 dλ, t > 0,
(4.2)
b→∞ 2πi a−ib
where a > 0 is arbitrary. Thus, the (continuous) spectrum of T lies to the left of the integration path alongx = a (see Fig. 1a).
15

For arbitrary φ, ψ ∈ H, we have
Z a+ib
1e φ, ψ H = limeλt (λ − T)−1 φ, ψ H dλ.
(4.3)
b→∞ 2πi a−ibFor φ, ψ ∈ X, (λ − T)−1 φ, ψ H is an analytic function in the right half–plane, which can be extended tothe entire complex plane as a meromorphic function hR(λ)φ, ψi. Thus,
Z a+ib1
eλt hR(λ)φ, ψidλ ∀φ, ψ ∈ X.
(4.4)
hetT φ, ψi = limb→∞ 2πi a−ibtTLet K ∈ [0, Kc+ ) be fixed. Next we claim that one can choose ε = ε(K) > 0 such that there are no
generalized eigenvalues of T on or inside the contour
a − iR → a + iR → −ε + iR → −ε − iR → a − iR

Cε,R :(Fig. 1b)for every R > 0. To construct Cε,R with the desired property, we first fix δ > 0. Then we recall thatgeneralized eigenvalues of T satisfy (3.12). From (3.13), under our assumptions on g, there exists R0 =
R0 (δ) such that there are no roots of (3.12) in the region+DR= {z ∈ C : |z| ≥ R0 & − δ < Re(z) ≤ a},0 ,δ
+for R  1. On the other hand, D(λ)
because (3.12) can be reduced to 2/(Kµ) = O(1/|λ|) in DR,δis holomorphic. Thus, the set of roots of (3.12) (i.e., the set of generalized eigenvalues) does not haveaccumulation points in−DR= {z ∈ C : |z| ≤ R0 & − δ < Re(z) ≤ a}.
0 ,δ
−+. This completes∪ DR
Thus, we can choose ε > 0 such that there are no generalized eigenvalues in DR
0 ,ε
0 ,ε
the construction of Cε,R .

By the Cauchy Integral theorem,Ieλt hR(λ)φ, ψidλ = 0 ∀φ, ψ ∈ X,

(4.5)

Cε,R

for any R > 0, andZ a+iRλt

Z −ε+iRe hR(λ)φ, ψidλ =a−iRZ −ε+iR−−ε−iR

Z a−iR −a+iR

eλt hR(λ)φ, ψidλ.

(4.6)

−ε−iR

Below, we show that the last two integrals on the right–hand side of (4.6) tend to 0 as R → ∞. SendingR → ∞ in (4.6) and using (4.4), we arrive at
Z −ε+iR
1tT
he φ, ψi = lim
eλt hR(λ)φ, ψidλR→∞ 2πi −ε−iR
Z R
(4.7)
e−εt
=
lim
ieiλt hR(iλ − ε)φ, ψidλ
2πi R→∞ −R= O(e−εt ),∀φ, ψ ∈ X16

as t → ∞ because

Z R

ieiλt hR(iλ − ε)φ, ψidλ = O(1).

−R

It remains to prove the following lemma.
Lemma 4.3. For K ∈ [0, Kc+ ),Z −ε+iR

Z a−iRλte hR(λ)φ, ψidλ = lim

lim

R→∞ −ε−iRR→∞ a+iR

eλt hR(λ)φ, ψidλ = 0 ∀φ, ψ ∈ X.

(4.8)

R −ε+iR
Proof. We show that the integral a+iR eλt hR(λ)φ, ψidλ tends to zero as R → ∞. The second integral
R a−iR λt−ε−iR e hR(λ)φ, ψidλ can be treated in the same way. Further, we decompose the integral into two
integrals as
Z −ε+iReλt hR(λ)φ, ψidλ =

Z iR

a+iR

eλt hR(λ)φ, ψidλ +a+iRZ −ε+iReλt hR(λ)φ, ψidλ.

(4.9)

iR

We show that the first integral on the right hand side tends to zero as R → ∞. For Re(λ) > 0, we have

−1
K
K
hR(λ)φ, ψi = hA(λ)φ, ψi + hA(λ)W I − D(λ)W
D(λ)φ, ψi,22see (3.20). For the first term, we haveZ iR
eλt hA(λ)φ, ψidλ
a+iR

Z 0

= eiRt

eλt

a

Z Z
I

1
φ(ω, x)ψ(ω, x)g(ω)dωdxdλ.
λ
+
i(R
− ω)
R

Since the integral above is finite, for any ε0 > 0, there exists L > 0 such that
Z 0

λt

Z Z

e
a

I

1
φ(ω, x)ψ(ω, x)g(ω)dωdxdλ < ε0
|ω|>L λ + i(R − ω)On the other hand, the integrandeλt

1φ(ω, x)ψ(ω, x)g(ω) → 0,
λ + i(R − ω)as R → ∞ uniformly in x ∈ I, ω ∈ (−L, L) and λ ∈ (0, a). This implies that the integral
Z iR

eλt hA(λ)φ, ψidλ → 0,

a+iR

17

as R → ∞.

Next, define a new function φeλ by−1K
φeλ = W I − D(λ)W
D(λ)φ.
2
For 0 < K < Kc+ , the operator T has no generalized eigenvalues in the right half plane and on the imaginary
−1axis. This and (2.14) shows that I − K2 D(λ)Wis bounded uniformly in λ such that Re(λ) ≥ 0. Hence,
R iR λte
φλ is in X for any λ. Then the integral a+iR e hA(λ)φeλ , ψidλ tends to zero by the same estimate as above.
R iR
This proves that a+iR eλt hR(λ)φ, ψidλ decays to zero as R → ∞. The second integral in (4.9) is estimatedin a similar manner. This completes the proof of Lemma 4.3.

5

Bifurcation with a one-dimensional null space

In the previous section, we proved asymptotic stability of the equilibrium at the origin of the linearized
system (2.8), (2.9) for K ∈ [0, Kc+ ). On the other hand, for K > Kc+ there is a positive eigenvalue inspectrum of the linearized problem (cf. [3]). This signals a bifurcation at Kc+ . This bifurcation is analyzed
in this present section. As in the classical KM, the loss of stability of the incoherent state at Kc+ and thedevelopment of partial synchronization for K > Kc+ is best seen in terms of the order parameter.
Throughout this section, we assume that the largest positive eigenvalue µmax of W with the eigenfunction wmax is simple. Furthermore, we assume that at Kc+ there is a (one-dimensional) smooth center
manifold of the equilibrium at the origin of (2.6), (2.7)3 . Under these assumptions, below we show thatthe order parameter undergoes a supercritical pitchfork bifurcation at Kc+ . The stable branch of equilibria
bifurcating from 0 is given byspg(0)2 π 3/2 3/2
1 p
h∞ (K) = p 00 µmax
K − Kc+ + o( K − Kc+ ), K > Kc+ ,
(5.1)
C(x)
−g (0)
where
C(x) :=

Π̃µmax (|wmax |2 wmax ).
|wmax |2 wmax

(5.2)

Formula (5.1) generalizes the classical Kuramoto’s formula describing the pitchfork bifurcation in the allto-all coupled model to the KM on graphs. The network structure enters into the description of the pitchfork
bifurcation through the largest eigenvalue µmax and the corresponding eigenspace.
3

The proof of existence of the center manifold is a technical problem and is beyond the scope of this paper (see [1] for the proofof existence of the center manifold in the original KM).

18

5.1

Preparation

Throughout this section, we assume that µmax is a simple eigenvalue of W. Let K = Kc+ +  with
0 <   1 and rewrite (2.6),(2.7) as follows∂
K
z1 = T0 z1 + Pz1 − Pz1 z2 ,
∂t
2
2

∂
jK
Pz1 zj−1 − Pz1 zj+1 ,zj = ijωzj +
∂t
2

(5.3)
j = 2, 3, . . . ,(5.4)

where T0 is T estimated at K = Kc and T = T0 + P/2.
For small  > 0, the equilibrium of (5.3), (5.4) at the origin has a 1D unstable manifold. We reducethe dynamics on the 1D unstable manifold, which we approximate by the center manifold of the origin for
K = Kc+ , i.e., for  = 0. For the latter, we assume zk = hk (z1 ), k = 2, 3, . . . , on the center manifold,where hk are smooth functions such that hk (0) = h0k (0) = 0.
Let Π0 be the projection to the eigenspace of λ = 0 spanned by vc+ (cf. Section 3.4). To track theevolution on the slow manifold we adopt the following Ansatz:
z1 = Π0 z1 + (I −Π0 )z1 = αc(t)vc+ + O(α2 ),2zk = hk (z1 ) = O(α ),k = 2, 3, . . . ,2

 = α ,

(5.5)
(5.6)
(5.7)

where α > 0 is a small parameter, c(t) is the coordinate along the center manifold, and vc+ is the generalized
eigenfunction of T0 corresponding to the zero eigenvalue (cf. (3.15)). The Ansatz (5.5)-(5.7) follows rightaway once existence of the center manifold is shown.
We will start by deriving several auxiliary facts that follow from the Ansatz (5.5)-(5.7). First, using(5.5)-(5.7) and Theorem 3.4, from (5.3), we have+22ż1 = T0 z1 + O(α2 ) = T×0 (αc(t)vc ) + O(α ) = O(α ).

(5.8)

Next, we estimate the order parameter.
Lemma 5.1.
h(t, x) = αc(t)wmax (x) + O(α2 ).

(5.9)


h = Pz1 = P αc(t)vc+ + O(α2 )Z Z
Kc+
W (x, y)wmax (y)g(ω)= αc(t)lim
dydω + O(α2 ).
2 λ→0+ R I
λ − iω

(5.10)

Proof.

19

Applying the Fubini theorem, (3.12) and (2.15), we have
Kc+
(Wwmax )D(0+) + O(α2 )2K + µmax= αc(t) c
D(0+)wmax + O(α2 )2= αc(t)wmax + O(α2 ).

h = αc(t)Lemma 5.2.
z2 =αc(t)Kc+
2

2

2
wmax
+ O(α3 ).
λ→0+ (λ − iω)2

lim*

(5.11)(5.12)

Proof. Using (5.5)-(5.7) and (5.8), we obtainż2 = h02 (z1 )ż1 = O(α3 ),(Pz1 )z3 = O(α3 ).

(5.13)

By plugging (5.13) into (5.4) for j = 2, we obtain0 = 2iωz2 + K(Pz1 )z1 + O(α3 ).

(5.14)

Next we plug in the expressions for z1 , Pz1 , and z2 (see (5.5), (5.9), (5.12)) into (5.14) to verify that theysatisfy this equation up to O(α3 ) terms. Specifically, we haveαc(t)Kc+
2

2

2
wmax
lim*
λ→0+ (λ − iω)2
+ K αc(t)wmax + O(α2 ) αc(t)vc+ + O(α2 ) + O(α3 )2iωz2 + K(Pz1 )z1 = 2iω

(Kc+ )2(λ − iω) − λ 2
lim*
wmax
λ→0+
2
(λ − iω)2
K+
wmax+ Kc+ α2 c(t)2 wmax c lim*
+ O(α3 )2 λ→0+ λ − iω
(K + )2
w2(K + )2w2= −α2 c(t)2 c
lim* max + α2 c(t)2 c
lim* max + O(α3 )2 λ→0+ λ − iω
2 λ→0+ λ − iω3
= O(α ).
= −α2 c(t)2

5.2

The slow manifold reduction

Projecting both sides of (5.3) onto the center subspace, we have
Π0 ż1 = Π0 T0 z1 +K
Π0 h −
Π0 (hz2 ).
2220(5.15)

Using (5.5), we have
Π0 ż1 = αċ(t)vc+ ,× +Π0 T0 z1 = T×0 Π0 z1 = αc(t)T0 vc = 0.

(5.16)

Further,
Π0 h = g1 lim* (λ − iω)−1 Π̃µmax D(λ)h.
λ→0+(5.17)

To evaluate (5.17), we take the following steps
Z

αc(t)wmax
g(ω)dω + O(α2 )λ→0+ Rλ − iω
= αc(t)wmax D(0+) + O(α2 )2αc(t)wmax+ O(α2 )=+Kc µmax

lim* D(λ)h = lim*

λ→0+

and
lim Π̃µmax D(λ)h =λ→0+2αc(t)wmax+ O(α2 ).
Kc+ µmax

Finally,2αc(t)wmax+ O(α2 )g1 lim*
+
λ→0+
λ
−
iω
Kc µmax


αc(t)
2 2 +
vc + O(α2 ).
=
g1
+
µmax
Kc

Π0 h =(5.18)

Similarly, to evaluate
Π0 (hz2 ) = g1 lim* (λ − iω)−1 Π̃µmax D(λ)(hz2 ),λ→0+(5.19)

we first compute
2Zg(ω)dω
Kc+
2
lim* D(λ)(hz2 ) = α |c(t)| c(t)|wmax | wmax lim+ O(α4 )λ→0+
λ→0+ R (λ − iω)32
 + 2Kc
3
2= −α |c(t)| c(t)g2|wmax |2 wmax + O(α4 ),
2
3

2



where
g2 =

πg 00 (0).
2

(5.20)(5.21)

By plugging (5.20) into (5.19), we obtainΠ0 (hz2 ) = −α3 |c(t)|2 c(t)g1 g2

Kc+ Π̃µmax (|wmax |2 wmax ) +vc + O(α4 ).
2
wmax
21

(5.22)

By plugging (5.16), (5.18), and (5.22) into (5.15), dividing both sides by α and vc+ and keeping terms
up to O(α2 ) we have!
2g1 c(Kc+ )4 µmax g2 2 2 Π̃µmax (|wmax |2 wmax )ċ =++ O(α3 ).
(5.23)
α |c|8
wmax(Kc+ )2 µmax
It is instructive to recast (5.23) in terms of the order parameter h (cf. (2.3)). By Lemma 5.1,h(t, x) = αc(t)wmax + O(α2 ).
Thus, by multiplying both sides of (5.23) by αwmax and neglecting higher order terms, we obtain!
2g1(Kc+ )4 µmax g2 Π̃µmax (|wmax |2 wmax ) 2
h +
ḣ =
|h|
8
|wmax |2 wmax
(Kc+ )2 µmax

2g1
g 00 (0)2=h + 3
C(x)|h| .
π g(0)4 µ3max(Kc+ )2 µmax

(5.24)

Equation (5.24) shows that the trivial solution (the incoherent state) looses stability at  = 0 and for small
 > 0 there is a nonzero stable equilibriumspp|wmax |2 wmax
−8
+
+
o(
·
K
−
K
K − Kc+ )|h∞ | =c
(Kc+ )4 µmax g2 Π̃µmax (|wmax |2 wmax )(5.25)spg(0)2 π 3/2 3/2
1 p= p 00 µmax
K − Kc+ + o( K − Kc+ ), K > Kc+ .
C(x)−g (0)

5.3

ExamplesIn [3], we derived the transition formulas for the onset of synchronization in the KM on several networks.
We now return to these examples and describe the transition to synchronization in more detail using theresults of this section.
We start with the KM on the Erdős-Rényi graphs. To this end, let W ≡ p ∈ (0, 1). In [3], we showedthat the largest positive eigenvalue of W in this case is µmax = p. The corresponding eigenfunction wmax
is constant. This yields the critical value Kc+ = 2(πg(0)p)−1 . By plugging in these values into (5.1), we
obtainppg(0)2 π 3/2
h∞ (K) = p 00 p3/2 K − Kc+ + o( K − Kc+ ).
(5.26)−g (0)
We next turn to the KM on small-world graphs. This family of graphs is defined via the following
graphon:
1 − p, min{|x − y|, 1 − |x − y|} ≤ r,
(5.27)
Wp,r (x, y) =p,otherwise,
where p, r ∈ (0, 1/2) are two parameters. The former stands for the probability of long range random
connections and the latter is the range of regular local connections (cf. [16]).
22

The largest eigenvalue of Wp,r is equal to 2r + 2p − 4pr and the corresponding eigenfunction wmax isconstant (cf. [3]). This implies that the critical value is
Kc+ =

2
.
πg(0)(2r + p − 4pr)Using (5.1), we further haveppg(0)2 π 3/2
h∞ (K) = p 00 (2r + 2p − 4pr)3/2 K − Kc+ + o( K − Kc+ ).
−g (0)

6

Bifurcation with a two-dimensional null space

6.1

The slow manifold reduction(5.28)

Many networks in applications can be described with the limiting graphon of the following form
W (x, y) = G(x − y)(6.1)

for some G ∈ L2 (S) such that G(x) = G(−x). The graphons of this form are used in the description of thesmall-world and many other networks (cf. §5.3 [3]).
A graphon satisfying (6.1) admits Fourier series expansionX
W (x, y) =ck e2πik(x−y) , c−k = ck ∈ R.

(6.2)

k∈Z

By Parseval’s identity,X

c2k = (2π)−1 kGk2L2 (S) < ∞.

(6.3)

k∈Z

It follows from (6.2) that the eigenvalues of the kernel operator W coincide with the Fourier coefficients
ck , k ∈ Z. The Fourier modes e2πikx , k ∈ Z, yield the corresponding eigenfunctions.
We continue to assume that the largest eigenvalue of W is positive, i.e., there is at least one positivecoefficient ck . In view of (6.3), there is a finite set
M = {m ∈ Z : cm = sup{ck : k ∈ Z}} .

(6.4)

If M = {0} the null space of W is one-dimensional. This case was analyzed in the previous section. Here,
we assume |M | = 2, i.e., there exists a unique m ∈ N such that
µmax = sup{ck : k ∈ Z} = cm = c−m .
The corresponding eigenspace is spanned by w+ = e2πimx and w− = e−2πimx .

23

From now on, the slow manifold reduction proceeds along the lines of the analysis in Section 5. The
generalized center subspace of T0 is spanned by
v± =

Kc+
w± (x)lim*.
2 λ→0+ λ − iω

(6.5)

On the center manifold, we adopt the following Anzats
z1 = α (c− (t)v− + c+ (t)v+ ) + O(α2 ),2zk = hk (z1 ) = O(α ),k = 2, 3, . . . ,2

 = α ,

(6.6)
(6.7)
(6.8)

where (c− , c+ ) is the coordinate along the center manifold. Following the lines of Lemma 5.1 and 5.2, we
obtainh(t, x) = α (c− (t)w− (x) + c+ (t)w+ (x)) + O(α2 )
(6.9)
and

2
1
αKc+
(c− (t)w− (x) + c+ (t)w+ (x)) · lim*
+ O(α3 ).
z2 =λ→0+ (λ − iω)22


(6.10)In analogy to (5.15) and (5.16), projection of (5.3) onto the center subspace yieldsα(ċ− v− + ċ+ v+ ) =KΠ0 h −
Π0 (hz2 ).
22(6.11)(c− v− + c+ v+ ) + O(α2 ).

(6.12)As in (5.18) and (5.22), we further obtain
αg1
Π0 h =
µmax



2
Kc+

2

and
3



Π0 (h̄z2 ) = −α g1 g2

Kc+
2

2

1
Π̃µmax ((c− w− +c+ w+ )2 (c− w− + c+ w+ ))+O(α4 ). (6.13)
λ→0+ λ − iω

· lim*Taking into account that w± = e±2πimx , we computeΠ̃µmax ((c− w− + c+ w+ )2 (c− w− + c+ w+ )) = c− |c− |2 + 2c− |c+ |2 w− + c+ |c+ |2 + 2c+ |c− |2 w+ .
(6.14)
Plugging (6.14) into (6.13), we obtainΠ0 (h̄z2 ) = −α3 g1 g2 Kc+
2 
c− |c− |2 + 2c− |c+ |2 v− + c+ |c+ |2 + 2c+ |c− |2 v+ + O(α4 ).

(6.15)

Combining (6.11), (6.12), and (6.15), and by comparing the coefficients of v± on both sides of theresultant equation, we arrive at the following system of equationsċ− = p1 c−  − α2 p2 |c− |2 + 2|c+ |2  + O(α3 ),(6.16)
ċ+ = p1 c+  − α2 p2 |c+ |2 + 2|c− |2 + O(α3 ),24where


2g1
p1 =
,
+ 2
(Kc ) µmax

p2 = −2g2

Kc+
2

4
µmax

(6.17)

are positive constants (see Theorem 3.6).
Using the polar form for c± = r± eiφ± , we rewrite (6.16) as follows

2 p r 2 + 2r 2
ṙ
=
p
r

−
α
+ O(α3 ),
−
1
−
2

−
+



2 + 2r 2
ṙ+ = p1 r+  − α2 p2 r++ O(α3 ),
−
φ̇
= 0,


 −
φ̇+ = 0.

(6.18)

Neglecting the higher order terms, we locate the fixed pointspppp(r− , r+ ) = (0, 0), (0, /p2 α2 ), ( /p2 α2 , 0), ( /3p2 α2 , /3p2 α2 ).
The linearization of (6.18) about these fixed points yields








2
1 0
−1 0
−2 0
−1 −2
p1
, p1
, p1
, p1
,
0 1
0 −2
0 −1
−2 −1
3
respectively. Thus, the second and the third fixed points are stable for 0 < K − Kc+  1. This proves thatthe order parameter (6.9) tends tosp
K − Kc+ iφh(t, x) =e w+ (x) + o( K − Kc+ )p2
orsh(t, x) =p
K − Kc+ iφe w− (x) + o( K − Kc+ )p2

as t → ∞, where φ is a constant which depends on an initial condition.

6.2

ExamplesBefore presenting the main example of this section, we briefly revisit small-world networks. The graphon
Wp,r (see (5.27)) satisfies (6.1) with the largest positive eigenvalue
Z Z
c0 =W (x, y)dxdy.
I

I

The corresponding eigenspace is one-dimensional. Thus, the bifurcation at K = Kc+ falls into the settingof the previous section. However, the eigenspace corresponding to the smallest negative eigenvalue is twodimensional for generic values of the parameters p and r. Consequently, the bifurcation at the other endpointof the stability region, K = Kc− < 0, is of the type considered in the present section. Interestingly, thebifurcations at Kc± are qualitatively different. We do not present a numerical example for the bifurcation
25

θi

a

6

6

6

5

5

5

4

θi

3

4

θi

3

4
3

2

2

2

1

1

1

20

40

60

80

100

b

20

40

60

80

100

c

20

40

60

80

100

Figure 2: Partially synchronized state in the KM on a) complete graph, b) small-world graph (p = r = 0.25),
and c) weighted graph constructed using kernel (6.19). The horizontal axis denotes the lables of oscillators,
and the vertical axis denotes the phases of them in the partially locked state.
at Kc− , because it is difficult to locate the smallest negative eigenvalue in this case. Instead, we construct
an example of the kernel operator with a two-dimensional eigenspace corresponding to the largest positive
eigenvalue.
Specifically, we consider
W (x, y) = G(x − y) = cos (2π(x − y)) .

(6.19)The Fourier expansion
1
1
W (x, y) = e−2πi(x−y) + e2πi(x−y)
22implies that the largest eigenvalue is 1/2, whose eignespace is spanned by
w± (x) = e±2πix .
The analysis of this section then yields the asymptotic state of the order parameter for 0 < K − Kc+  1 :s
pK − Kc+ iφ ±2πix
8g 00 (0)h∞ (x) =·e e+ o( K − Kc+ ), p2 = − 3,p2
π g(0)4
where φ ∈ R is a constant depending on the initial data. Using (1.16), we now can describe the phase-locked
solutions:
 r

ω
p2
,
θ = ±2πx + φ + arcsin
K K − Kc+q+cprovided |ω| < K K−K
p2 . In particular, in the KM with intrinsic frequencies having approximately thesame value the phase-locked solutions are approximately given by θ = ±2πx + constant. Such solutions
are called twisted states [17, 22].
Figure 6.2 presents numerical results for the KM (1.2) on complete graph (a), small-world graph (b),and the weighted graph constructed using kernel (6.19) (c). In these simulations, we used 100 coupled26

oscillators. The coupling strength was taken K = 5, which guaranteed stability of the partially locked
states. The horizontal axis denotes the lables of oscillators i = 1, · · · , 100, and the vertical axis denotes thephases θi of them in the partially locked state. In plots a and b the phases are distributed around θ = 0,
while in the plot c, the phases are close to the twisted state
θi = −2πi/N + π,

i = 1, 2, . . . , N.

Acknowledgements. This work was supported in part by the NSF DMS grants 1412066 and 1715161 (toGM).

References[1] 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[2], A spectral theory of linear operators on rigged Hilbert spaces under analyticity conditions,Adv. Math. 273 (2015), 324–379. MR 3311766[3] Hayato Chiba and Georgi S. Medvedev, The mean field analysis of the kuramoto model on graphs I.
The mean field equation and the transition point formulas, submitted.
[4] Hayato Chiba and Isao Nishikawa, Center manifold reduction for large populations of globally coupledphase oscillators, Chaos 21 (2011), no. 4, 043103, 10. MR 2918832[5] H. Dietert, Stability of partially locked states in the Kuramoto model through Landau damping with
Sobolev regularity, ArXiv e-prints (2017).
[6] Helge Dietert, Stability and bifurcation for the Kuramoto model, J. Math. Pures Appl. (9) 105 (2016),
no. 4, 451–489. MR 3471147[7] 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[8] Bastien Fernandez, David Gérard-Varet, and Giambattista Giacomin, Landau damping in the Kuramoto model, Ann. Henri Poincaré 17 (2016), no. 7, 1793–1823. MR 3510470[9] F. D. Gakhov, Boundary value problems, Dover Publications, Inc., New York, 1990, Translated fromthe Russian, Reprint of the 1966 translation. MR 1106850[10] I. M. Gel0 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[11] 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
27

[12] L. Lovász, Large networks and graph limits, AMS, Providence, RI, 2012.
[13] 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)[14] 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
[15]
[16]

, 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[17] G.S. Medvedev and X. Tang, Stability of twisted states in the Kuramoto model on Cayley and randomgraphs, Journal of Nonlinear Science (2015).
[18] Clément Mouhot and Cédric Villani, On Landau damping, Acta Math. 207 (2011), no. 1, 29–201. MR
2863910[19] 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[20] 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[21] 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[22] D.A. Wiley, S.H. Strogatz, and M. Girvan, The size of the sync basin, Chaos 16 (2006), no. 1, 015103,
8. MR 2220552 (2007e:37016)

28