Bifurcations in the Kuramoto model on graphs
Hayato Chiba, Georgi S. Medvedev, and Matthew S. Mizuhara

Citation: Chaos 28, 073109 (2018); doi: 10.1063/1.5039609
View online: https://doi.org/10.1063/1.5039609
View Table of Contents: http://aip.scitation.org/toc/cha/28/7
Published by the American Institute of Physics

CHAOS 28, 073109 (2018)

Bifurcations in the Kuramoto model on graphs
Hayato Chiba,1,a) Georgi S. Medvedev,2,b) and Matthew S. Mizuhara3,c)1

Institute of Mathematics for Industry, Kyushu University/JST PRESTO, Fukuoka 819-0395, Japan
Department of Mathematics, Drexel University, 3141 Chestnut Street, Philadelphia, Pennsylvania 19104, USA
3
Department of Mathematics and Statistics, The College of New Jersey, 2000 Pennington Road, Ewing, New Jersey 08628-0718,
USA
2(Received 9 May 2018; accepted 28 June 2018; published online 27 July 2018)
In his classical work, Kuramoto analytically described the onset of synchronization in all-to-all coupled networks of phase oscillators with random intrinsic frequencies. Specifically, he identified acritical value of the coupling strength, at which the incoherent state loses stability and a gradualbuild-up of coherence begins. Recently, Kuramoto’s scenario was shown to hold for a large class ofcoupled systems on convergent families of deterministic and random graphs [Chiba and Medvedev,“The mean field analysis of the Kuramoto model on graphs. I. The mean field equation and the transition point formulas,” Discrete and Continuous Dynamical Systems—Series A (to be published); “Themean field analysis of the Kuramoto model on graphs. II. Asymptotic stability of the incoherent state,
center manifold reduction, and bifurcations,” Discrete and Continuous Dynamical Systems—Series
A (submitted).]. Guided by these results, in the present work, we study several model problems illustrating the link between network topology and synchronization in coupled dynamical systems. First,
we identify several families of graphs, for which the transition to synchronization in the Kuramotomodel starts at the same critical value of the coupling strength and proceeds in a similar manner.
These examples include Erdős-Rényi random graphs, Paley graphs, complete bipartite graphs, andcertain stochastic block graphs. These examples illustrate that some rather simple structural properties such as the volume of the graph may determine the onset of synchronization, while finer structuralfeatures may affect only higher order statistics of the transition to synchronization. Furthermore, we
study the transition to synchronization in the Kuramoto model on power law and small-world randomgraphs. The former family of graphs endows the Kuramoto model with very good synchronizability:
the synchronization threshold can be made arbitrarily low by varying the parameter of the power law
degree distribution. For the Kuramoto model on small-world graphs, in addition to the transition tosynchronization, we identify a new bifurcation leading to stable random twisted states. The examples analyzed in this work complement the results in Chiba and Medvedev, “The mean field analysisof the Kuramoto model on graphs. I. The mean field equation and the transition point formulas,”Discrete and Continuous Dynamical Systems—Series A (to be published); “The mean field analysisof the Kuramoto model on graphs. II. Asymptotic stability of the incoherent state, center manifold
reduction, and bifurcations,” Discrete and Continuous Dynamical Systems—Series A (submitted).
Published by AIP Publishing. https://doi.org/10.1063/1.5039609
Understanding principles of collective dynamics and synchronization in complex networks is a challenging problem of nonlinear science with applications ranging fromneuronal networks to power grids. The Kuramoto modelof coupled phase oscillators provides a successful framework for studying synchronization in diverse biological,
physical, and social networks. In particular, it reveals auniversal mechanism for the transition to synchronizationin ensembles of coupled oscillators with random intrinsic frequencies. A pitchfork bifurcation responsible forthe onset of synchronization in the Kuramoto model of
all-to-all coupled phase oscillators is a classical result.
This paper explains the onset of synchronization in the
Kuramoto model on a broad class of dense and sparse
random graphs. It provides explicit formulas characterizing the onset of synchronization for several familiesa)chiba@imi.kyushu-u.ac.jp
medvedev@drexel.edu
c)mizuharm@tcnj.edub)1054-1500/2018/28(7)/073109/10/$30.00

of random graphs, which are important in applications,
including small-world, power law, Erdős-Rényi, and Paleygraphs, among others. For the Kuramoto model on smallworld graphs, in addition to the transition to synchronization, a new bifurcation leading to stable random twisted
states is identified. The examples analyzed in this work
elucidate the relation between the structure and dynamics
in complex networks.

I. INTRODUCTION

The Kuramoto model (KM) of coupled phase oscillators provides an important paradigm for studying collective
dynamics and synchronization in ensembles of interacting
dynamical systems. In its original form, the KM describes
the dynamics of all-to-all coupled phase oscillators with

28, 073109-1

Published by AIP Publishing.

073109-2

Chiba, Medvedev, and Mizuhara

Chaos 28, 073109 (2018)

randomly distributed intrinsic frequenciesθ̇i = ωi + Kn−1nsin(θj − θi ),where

i ∈ [n] := {1, 2, . . . , n}.

√r∞ (x) = r1 x + O(x),

j=1

(1)
Here, θi : R → T := R/2π Z is the phase of the oscillatori ∈ [n]. The intrinsic frequencies ωi , i ∈ [n], are independent identically distributed random variables. Throughout thispaper, we assume that the density of the absolutely continuous probability distribution of ω1 , g, is an even smoothfunction. The sum on the right-hand side of (1) models theinteractions between oscillator i and the rest of the network.
Parameter K controls the strength of the interactions. The
sign of K determines the type of interactions. The couplingis attractive if K > 0, and is repulsive otherwise. Sufficiently
strong attractive coupling favors synchrony.
For small positive values of K, the KM shows little coherence. The phases fill out the unit circle approximately uniformly [Fig. 1(a)]. This dynamical regime is called
an incoherent state. It persists for positive values of K
smaller than the critical value Kc = 2[π g(0)]−1 . For values of
K > Kc , the system undergoes a gradual build-up of coherence approaching complete synchronization as K → ∞[Fig. 1(b)]. Kuramoto identified the critical value Kc anddescribed the transition to synchrony, using the complex orderparameterh(t) = n−1neiθj (t) .

(2)

j=1

Below, we will often use the polar form of the order parameterh(t) = r(t)eiψ(t) .

(3)Specifically, he showed that for t  1 (cf. Ref. 20)

0 < K < Kc ,O(n−1/2 ),r(t) =
r∞ (K − Kc ) + O(n−1/2 ), K > Kc ,
(a)

(b)

(4)

x ≥ 0,

r1 =

Kc2

4

. (5)
−π g (0)Here, 0 ≤ r(t) ≤ 1 and ψ(t) stand for the modulus and theargument of the order parameter, respectively, defined by theright-hand side of (2). The value of r is interpreted as a measure of coherence in the system dynamics. Indeed, if all phases
are independent random variables distributed uniformly on T(complete incoherence), then r = o(1) with probability one,by the Strong Law of Large Numbers. If, on the other hand,
all phase variables are equal, then r = 1. Equations (4) and(5) suggest that the system undergoes a pitchfork bifurcation en route to synchronization. This bifurcation is clearlyseen in the numerical experiments [see Fig. 1(c)]. Furthermore, Kuramoto’s scenario of the transition to synchrony wasrecently confirmed by rigorous mathematical analysis of theKM.3,6,7
In this paper, we investigate bifurcations in the KM on
a variety of graphs ranging from symmetric Cayley graphs
to random small-world and power law graphs. To visualize
the connectivity of a large graph, we use the pixel pictureof a graph. This is a square black and white plot, where
each black pixel stands for entry 1 in the adjacency matrixof the graph (see Fig. 2). As we show below, the graphstructure plays a role in the transition to synchrony. Furthermore, some graphs may force bifurcations to spatial patterns
other than synchrony. To highlight the role of the network
topology, we present numerical experiments of the KM on
small-world graphs.21 These graphs are formed by replacingsome of the edges of a regular graph with random edges [seeFig. 2(a)]. These results demonstrate that some graphs may
exhibit bifurcations to spatial patterns other than synchrony.
Figure 3(a) shows that for the KM on small-world graphs
like in that on complete graphs, the order parameter undergoes a smooth transition. Near the critical value Kc+ ≈ 3.2,
the asymptotic in time value of the order parameter starts to(c)FIG. 1. The distribution of the phases of coupled oscillators is shown on the unit circle in the complex plane: θk → eiθk ∈ C, k ∈ [n]. The strength of coupling
is below the critical value Kc in (a) and is above Kc in (b). The black arrow depicts the order parameter, as a vector in the complex plane [cf. (2)]. The bigger size
(modulus) of the order parameter corresponds to the higher degree of coherence. In (a), the modulus of the order parameter is close to zero, and the distributionof the oscillators is close to the uniform distribution. In contrast, in (b), the distribution develops a region of higher density. The complex order parameter points
to the expected value (center of mass) of the distribution of the oscillators around T. (c) The modulus of the order parameter is plotted for different values ofK. For the KM on graphs, the order parameter must be suitably redefined [see (21)]. The maximum value of the order parameter for the KM on graphs depends
on the edge density of the graph [cf. (21)]. The gradual change of the modulus of the order parameter from values around zero to those close to 1 marks thetransition to synchronization. This and all other numerical bifurcation diagrams in this paper were computed using the KM with 4001 coupled oscillators. The
realization of ωi ’s was generated for each value of K separately.

073109-3(a)Chiba, Medvedev, and Mizuhara

Chaos 28, 073109 (2018)(b)

FIG. 2. Pixel pictures of small-world(a) and power law (b) graphs on 100
nodes.

grow monotonically and approaches a value slightly greaterthan 0.5 for increasing values of K. This change in the orderparameter corresponds to the transition to synchronization.
Interestingly, in contrast to the original KM, the model on
small-world graphs exhibits another transition at a negativevalue of K, Kc ≈ −27 [Fig. 3(c)]. This time it grows monotonically for decreasing values of K and approaches the value
approximately equal to 0.0475. The corresponding state of
the network is shown in Fig. 3(c). In the pattern shown inthis figure, the oscillators are distributed randomly about a2-twisted state. A q-twisted state (q ∈ Z) is a linear function
on the unit circle: Tq (x) = 2πnqx mod 2π , x ∈ [n], q-twisted
states have been studied before as stable steady states in
repulsively coupled KM on Cayley graphs.17,22 By changing
parameters controlling small-world connectivity, we find different q-twisted states bifurcating from the incoherent state fordecreasing negative K [see Figs. 4(a) and 4(b)]. These numerical experiments show that in the KM on small-world graphs,
the incoherent state is stable for K ∈ (Kc− , Kc+ ). For values
K > Kc+ , the system undergoes the transition to synchrony,while for K < Kc− , it leads to random patterns localizedaround q-twisted states.
The loss of stability of the incoherent state in the KM on
spatially structured networks was analyzed in Refs. 4 and 5.
The linear stability analysis of the incoherent state yields thecritical values Kc± ,4 and the analysis of the bifurcations at Kc±
explains the emerging spatial patterns.5 In the present work,
using the insights from Refs. 4 and 5, we conduct numerical experiments illustrating the role of the network topology
in shaping the transition from the incoherent state to coherent structures like synchronous and twisted states in the KM
on graphs. In particular, we present numerical simulations
showing for many distinct graphs the onset of synchronization
takes place at the same critical value and the graph structure has only higher order effects. We identify the dominant
structural properties of the graphs shaping the transition tosynchronization in the KM.
The effective analysis of the KM on large graphs requires
an analytically convenient description of graph sequences. To
this end, we employ the approach developed in Refs. 14 and
15 inspired by the theory of graph limits11 and, in particular,
by W-random graphs.12 We explain the graph models used inthis paper in Sec. II. The W-random graph framework affords
an analytically tractable mean field description of the KM on abroad class of dense and sparse graphs. The mean field partial

differential equation for the KM on graphs is explained inSec. III. There, we also explain the generalization of the orderparameter suitable for the KM on graphs. After that we goover the main results of Refs. 4 and 5. In Sec. IV, we present
our numerical results. We use carefully designed examplesof graphs to highlight the structural properties affecting theonset of synchrony. For instance, we show that the KMs on
Erdős-Renyi and Paley graphs have the same critical value
Kc marking the onset of synchronization. Furthermore, themean field limits of the KM on these graphs coincide withthat for the KM on weighted complete graph. We demonstrate
the similarity in the transition to synchronization in the KM on
bipartite graphs and on the family of stochastic block graphsinterpolating a disconnected graph and a weighted completegraph. This time the mean field equations for the two models are different, but the critical value Kc remains the same.
Furthermore, we investigate the bifurcations in the KM on
small-world and power law graphs illuminating the effectsof these random connectivity patterns on the dynamics of the
KM.
II. THE KM ON GRAPHS

We begin with the description of graph models thatwill be used in this paper. Let n = V (n ), E(n ), An be aweighted directed graph on n nodes. V (n ) = [n] stands forthe node set of n . An = (an,ij ) is an n × n weight matrix. Theedge setE(n ) = {i, j} ∈ [n]2 : an,ij = 0 .
If n is a simple graph, then An is the adjacency matrix.
The KM on n is defined as follows:−1u̇n,i = ωi + (nαn )nan,ij sin(un,j − un,i ),i ∈ [n],(6)

j=1

un,i (0) = u0n,i .

(7)

The sequence αn = 0 is only needed and {n } is a sequenceof sparse graphs. Without proper rescaling, the continuum
limit (as n → ∞) of the KM on sparse graphs is trivial. The
sequence {αn } will be explained below, when we introduce
sparse W -random graphs. Until then, one can assume αn ≡ 1.
We will be studying solutions of (6) in the limit as
n → ∞. Clearly, such limiting behavior is possible only if
the sequence of graphs {n } is convergent in the appropriate

073109-4

Chiba, Medvedev, and Mizuhara(a)Chaos 28, 073109 (2018)(b)

FIG. 3. Transition to synchronization in
the KM on small-world graphs (a). For
negative K, the KM undergoes anothertransition at Kc− < 0 (b). The emerging
pattern is shown in Figs. 4(a) and 4(b).

sense. We want to define {n } in such a way that it includes
graphs of interest in applications and, at the same time, is convenient for deriving the continuum limit. To achieve this goal,
we use the ideas from the theory of graph limits.11 Specifically, we choose a square integrable function W on a unit
square I 2 := [0, 1]2 . W is used to define the asymptotic behavior of {n } as n → ∞. It is called a graphon in the theory ofgraph limits.11 Let

distributed (IID) random variables such that
P an,ij = 1 = W In,i ×In,j and
P an,ij = 0 = 1 − P an,ij = 1 .

(9)(RS) Sparse random graphs n = Hr (n, W , {αn }). Let W ∈ W
be a nonnegative function and positive sequence 1 ≥αn  0 satisfy nαn → ∞ as n → ∞W := W ∈ L2 (I 2 ) : W (x, y) = W (y, x) a.e. .

P (i → j) = αn W̃n In,i ×In,j

We define W0 ⊂ W as a set of measurable symmetric functions on I 2 with values in I. We discretize the unit interval
I: xn,j = j/n, j ∈ {0} ∪ [n] and denote In,i := (xn,i−1 , xn,i ], i ∈[n].
The following graph models are inspired by W-randomgraphs.2,12 Let n = V (n ), E(n ), An = (an,ij ) .

We illustrate the random graph models in (RD) and (RS) withthe following examples.
Example II.1. (SW ) Small-worldgraph
Sn,p,r =Hr (n, Wp,r ), r ∈ (0, 0.5), p ∈ (0, 0.5], is defined via

1 − p min{|x − y|, 1 − |x − y|} ≤ r,
Wp,r (x, y) =potherwise(11)[see Fig. 2(a)].
(ER) Erdős-Rényi graph Gn,p = Hr (n, Wp ), Wp ≡ p ∈ (0, 1).
(SER) Sparse Erdős-Rényi graph G̃n,p = Hr (n, Wp , {n−β }),β ∈ (0, 1).
(PL) Power law graph  = Hr (n, W γ , {n−β })W̃n (x, y) := min{αn−1 , W (x, y)}.

(DD) Weighted deterministic graph n = H(n, W ) is defined
as follows:
an,ij = W In,i ×In,j := n2In,i ×In,j

W (x, y)dxdy,(8)where W ∈ W .
(RD) Dense random graph n = Hr (n, W ). Let W ∈ W0
and an,ij , 1 ≤ i ≤ j ≤ n, be independent identically
(a)

(b)

and

W (x, y) = (xy)−γ ,0<γ <β<1

(10)

(12)

[see Fig. 2(b)].
(c)

FIG. 4. The oscillators in (a) are clustered around a 2-twisted state (p = 0.2, r = 0.3, K = −36). By changing one of the parameters controlling small-world
connectivity (namely, by reducing the range of local interactions), one can obtain q-twisted states for different values of q ∈ Z, such as the 4-twisted state shownin plot (b) (p = 0.2, r = 0.2, K = −40).

073109-5

Chiba, Medvedev, and Mizuhara

Chaos 28, 073109 (2018)

The graphon W carries all the information needed for thederivation of the continuum limit for the KM on deterministic
and random graphs (DD), (RD), and (RS). However, for theKM on random graphs (RD) and (RS) taking the continuum
limit involves an additional step. Namely, we first approximate the model on random graph by that on the averaged
deterministic weighted graph n = H(n, W )v̇n,i = ωn,i + Kn−1nW̄n,ij sin(vn,j − vn,i ),i ∈ [n],

(13)

j=1

vn,i (0) = u0n,i ,(14)where W̄n,ij := W̃n In,i ×In,j .
The approximate model (13) is formally derived fromthe original KM (6) on a random graph n = Hr (n, W , {αn })by averaging the right hand side of (6) over all possible
realizations of n . The justification of averaging is given in
Ref. 8.
III. THE MEAN FIELD LIMIT

Solutions of the KM with distributed frequencies such as
incoherent state and the solutions bifurcating from the incoherent state are best described in statistical terms. To this end,
suppose ρ(t, u, ω, x) is the conditional density of the random
vector (u, ω) given ω, and parametrized by (t, x) ∈ R+ × I.
Here, the spatial domain I = [0, 1] represents the continuumof the oscillators in the limit n → ∞ (see Ref. 4 for precise meaning of the continuum limit). Then, ρ satisfies thefollowing initial value problem:∂∂{ρ(t, u, ω, x)V (t, u, ω, x)} = 0,
ρ(t, u, ω, x) +∂t∂u
ρ(0, u, ω, x) = ρ 0 (u),

(15)
(16)

where initial density ρ 0 (u) independent of ω and x for simplicity (see Ref. 8 for the treatment of a more general case).
The velocity field V is defined via  V (t, u, ω, x) = ω + K
W (x, y) sin(φ − u)IRTρ(t, φ, ω, y)g(ω)dφdωdy.

(17)

The density ρ(t, ·) approximates the distribution of the oscillators around T at time t ∈ [0, T]. Specifically, in the large nlimit the empirical measures defined on the Borel σ -algebra
B (G), G := T × R:
μnt (A) = n−1nδ[uni (t),ωi ,xni ] (A),A ∈ B (G),R

Recall that W ∈ L (I ) is a symmetric function. Thus,
W : L2 (I 2 ) → L2 (I 2 ) is a self-adjoint compact operator. The
eigenvalues of W, on one hand, carry the information about
the structure of the graphs in the sequence {n }; on the other
hand, they appear in the stability analysis of the incoherent
state. Thus, through the eigenvalues of W, we can trace therelation between the structure of the network and the onset ofsynchronization in the KM on graphs.
Since W is self-adjoint and compact, the eigenvalues of
W are real and bounded with the only possible accumulation
point at 0. In all examples considered in this paper, the largest
eigenvalue μmax > 0 and the smallest μmin ≤ 0. The linear
stability analysis in Ref. 4 shows that the incoherent state isstable for K ∈ [Kc− , Kc+ ], where
2

Kc− =

2

2
π g(0)μmin

and Kc+ =2.
π g(0)μmax

Except for small-world graphs, for all other graphs consideredbelow, the smallest eigenvalue μmin = 0. Thus, the incoherent
state in the KM on these graphs is stable for K ≤ Kc +. In thissection, we comment on the bifurcation at Kc+ and postpone
the discussion of the bifurcation at Kc− until Sec. V, where wedeal with the KM on small-world graphs.
The analysis of the bifurcations in the KM on graphs
relies on the appropriate generalization of Kuramoto’s orderparameter (cf. Ref. 4)
hn (t) = [hn1 (t), hn2 (t), . . . , hnn (t)],1hni (t) =Wn,ij eiθnj (t) ,n j=1ni ∈ [n],

(21)

(18)

and its continuous counterpart  
W (x, y)eiθ ρ(t, θ, ω, y)g(ω)dθdωdy. (22)h(t, x) =(19)

The value of the continuous order parameter (22) h(t, x)
carries the information about the (local) degree of coherence at point x and time t. It is adapted to a given network
connectivity through the kernel W .
The main challenge in the stability analysis of the incoherent state lies in the fact that for K ∈ [Kc− , Kc+ ], the linearized operator has continuous spectrum on the imaginary

i=1

converge weakly to the absolutely continuous measureμt (A) = ρ(t, u, ω, x)g(ω)dudωdx,

state in the mean field description corresponds to the stationary density ρu = g(ω)/2π . The stability analysis of ρu ,a steady state solution of (15), yields the region of stabilityof the incoherent state and the critical values of K, at which
it loses stability. Furthermore, the analysis of the bifurcations
at the critical values of K explains the phase transitions in theKM (cf. Ref. 5). Importantly, we can trace the role of the network topology in the loss of stability of the incoherent stateand emerging spatial patterns. In the remainder of this section,
we informally review some of the results of Refs. 4 and 5,
which will be used below.
The key ingredient in the analysis of the KM on graphs,
which was not used in the analysis of the original KM, is thefollowing kernel operator W : L2 (I) → L2 (I):W[f ](x) =W (x, y)f (y)dy.
(20)

I

A

provided the initial data (7) converge weakly to μ0 [cf. Ref. 4(Theorem 2.2), see also Ref. 8].
The mean field limit provides a powerful tool for studying
the KM. It gives a simple analytically tractable description ofcomplex dynamics in this model. In particular, the incoherentRT073109-6

Chiba, Medvedev, and Mizuhara

Chaos 28, 073109 (2018)

(a)

(b)

FIG. 5. Pixel pictures of Erdős-Rényi (a) and Paley (b) graphs.

axis and no eigenvalues (cf. Ref. 4). To overcome this difficulty, in Ref. 5, the generalized spectral theory was used to
identify generalized eigenvalues responsible for the instability
at Kc+ and to construct the finite-dimensional center manifold. The explanation of these results is beyond the scopeof this paper. An interested reader is referred to Ref. 5 fordetails. Here, we only present the main outcome of the bifurcation analysis. The center manifold reduction for the orderparameter near Kc+ yields a stable branch of solutions

g(0)2 π 3/2 3/2
h∞ (x, K) =  
μmax C(x) K − Kc+ + o( K − Kc+ ),−g (0)
0 < K − Kc+  1,

(23)

where C(x) is defined through the eigenfunction of W corresponding to μmax (see Ref. 5 for details). Formula (23)generalizes the classical Kuramoto’s formula describing thepitchfork bifurcation in the all-to-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.
IV. RESULTS

In this section, we present numerical results elucidatingsome of the implications of the bifurcation analysis outlined
in Sec. III.
(a)

(b)

A. Graphs approximating the complete graph

Our first set of examples deals with the KM on ErdősRényi and Paley graphs (see Fig. 5). The former is a randomgraph, whose edges are selected at random from the set of allpairs of vertices with fixed probability p (see Example II.1).
Before giving the definition of the Paley graphs, we first recall
the definition of a Cayley graph on an additive cyclic group
Zn = Z/nZ.
Definition IV.1. Let S be a symmetric subset of Zn (i.e.,
S = −S). Then, n = V (n ), E(n ) is a Cayley graph if
V (n ) = Zn and {a, b} ∈ E(n ), if b − a ∈ S. Cayley graphn is denoted Cay(Zn , S).
Definition IV.2. Let n = 1 (mod 4) be a prime anddenoteQn = {x2 (mod n) : x ∈ Zn /{0}}.
Qn is viewed as a set (not multiset), i.e., each element has multiplicity 1. Then, Qn is a symmetric subset of Zn× and |Qn | =2−1 (n − 1) [cf. Ref. 10 (Lemma 7.22)]. Pn = Cay(Zn , Qn ) is
called a Paley graph.10
As all Cayley graphs, Paley graphs are highly symmetric[Fig. 5(b)]. However, Paley graphs also belong to pseudorandom graphs, which share many asymptotic properties with
Erdős-Rényi graphs with p = 1/2.1 In particular, both ErdősRényi and Paley graphs have constant graph limit W1/2 ≡ 1/2
as n → ∞. The mean field limit for the KM for Erdős-Rényi
and Paley graphs has the same kernel constant W = W1/2 ≡1/2 [cf. (15)], exactly as in the case of the KM on weighted(c)FIG. 6. A large time asymptotic value of the order parameter (21) (in a scaled l2 -norm) is plotted for the KM on complete (a), Erdős-Rényi (b), Paley (c). The
transition from 0 to 1 takes place in the same region of K in all three plots. The critical value of the coupling strength Kc+ ≈ 3.2 is the same for all three models.

073109-7

Chiba, Medvedev, and Mizuhara

Chaos 28, 073109 (2018)

(a)

(b)

FIG. 7. Pixel pictures of the bipartite graph (a) and stochastic block graph with two weakly connected components (b).

complete graphs. Since μmax (W1/2 ) = 1/2, the transition tosynchronization takes place at
Kc+ =4≈ 3.183
π g(0)

C. Power law graphs

for all three networks (see Fig. 6).
B. Bipartite and disconnected graphs

Our second set of examples features another pair of distinct networks, which are not pseudorandom, and yet have thesame synchronization threshold.
Here, we consider the KM on bipartite graphs Bn,n =H(2n, Wbp ) [see Fig. 7(a)]
1, [0, 0.5) × [0.5, 1] ∪ [0.5, 1] × [0, 0.5),
Wbp (x, y) =0 otherwise.
Next, we consider the family of graphs interpolating between
the disconnected graph with two equal components and aweighted complete graphs, n,α = Hr (n, Wd,α ):1 − α,
Wd,α (x, y) =α[0, 0.5) × [0, 0.5) ∪ [0.5, 1] × [0.5, 1],
otherwise,α ∈ [0, 0.5].
A simple calculation yields largest positive eigenvalues
for each of these families of graphs
μmax (Wbp ) = μmax (Wd,α ) = 0.5 ∀ α ∈ [0, 0.5].
(a)Thus, as in the first set of examples, the KM on Bn,n and n,αundergoes transition to synchronization at Kc+ = 4[π g(0)]−1(see Fig. 8).

Next, we turn to the KM on power law graphs. To generate graphs with power law degree distribution, we use
the method of sparse W-random graphs.2 Specifically, n =Hr (n, W γ , n−β ) with graphon W γ defined in (12). For graphs
constructed using this method, it is known that the expected
degree of node i ∈ [n] is O(n1+γ −β i−γ ) and the edge density isO(n−β ) (cf. Ref. 9). Thus, n = Hr (n, W γ , n−β ) is a family ofsparse graphs with power law degree distribution. The mean
field limit for the KM on n is given by (15) with W := W γ .
The analysis of the integral operator W [cf. (20)] with kernel W := W γ shows that it has a single nonzero eigenvalue(cf. Ref. 13)
μmax = (1 − 2γ )−1 ,γ ∈ (0, 0.5).

Thus, the synchronization threshold for the KM on power law
graphs is
Kc+ =2(1 − 2γ ).
π g(0)

Note that as α → 1/2, Kc+ tends to 0. Thus, the KM on powerlaw graphs features remarkable synchronizability despite
sparse connectivity. Numerics in Fig. 9 illustrate the onset ofsynchronization in power law networks.

(b)FIG. 8. The onset of synchronizationin the KM on bipartite graph (a) andon the stochastic block graph with two
weakly connected components (b) with
α = 0.05.

073109-8

Chiba, Medvedev, and Mizuhara

Chaos 28, 073109 (2018)

eigenspace corresponding to μ = μ0 is spanned by w0 = 1.
For μ = μ0 , the eigenspace corresponding to μ is spanned by
wk = ei2π kx

andw−k = e−i2π kx ,k ∈ Iμ ,(27)

where Iμ = {k ∈ N : μk = μ}. The largest positive eigenvalue of Wp,r is
μ0 = Kp,r (x)dx = 2(r + p − 2pr).
I

Therefore, the onset of synchronization in the KM on smallworld graphs takes place at1
π g(0)(r + p − 2pr)

Kc+ =FIG. 9. The onset of synchronization in the KM on power law graphs.
The parameter values γ = 0.4 and β = 0.6 were used in these numericalexperiments.

V. SMALL-WORLD GRAPHS

Small-world graphs are obtained by random rewiring of
regular k-nearest-neighbor networks.21 Consider a graph on
n nodes arranged around a circle, with each node connected
to k = rn neighbors from each side for some r ∈ (0, 0.5).
With probability p ∈ [0, 0.5], each of these edges connecting k neighbors from each side is then replaced by a random
long-range (i.e., going outside the k-neighborhood) edge. The
pixel picture of a representative small-world graph is shownin Fig. 2(a). A family of small-world graphs parametrized
by p interpolates between regular k-nearest-neighbor graph(p = 0) and a fully random Erdős-Rényi graph (p = 0.5). Forintermediate values of p, small-world graphs combine features of regular and random connectivity, just as seen in many
real-world networks.21
It is convenient to interpret a small-world graph as a
W -random graph Sn,p,r = Hr (n, Wp,r ) [cf. (11)]. The mean
field limit of the KM on Sn,p,r is then given by (15) with
W := Wp,r . The corresponding kernel operator is given by
Wp,r (f ) =Kp,r (· − y)f (y)dy, T1 := R/Z,
(24)
T1

where Kp,r (x) = (1 − p)1|x|≤r (x) + p1|x|>r (x).
Using (24), we recast the eigenvalue problem for Wp,r as
follows:
Kp,r ∗ v = μv.

(25)

By taking the Fourier transform of both sides of (25), we findthe eigenvalues of Wp,rμk = Kp,r (x)e−2π ikx dxT
=2(r + p − 2rp),k = 0,(π k)−1 (1 − 2p) sin(2π kr), k ∈ Z/{0}.

(26)The corresponding eigenvectors are wk = ei2π kx , k ∈ Z. Notethat μk = μ−k , since Kp,r is even [cf. (26)]. Thus, the[see Fig. 3(b)].
Importantly, Wp,r also has negative eigenvalues. Since
Wp,r is a compact operator, 0 is the only accumulation pointof the spectrum of Wp,r . Thus, there is the smallest (negative)
eigenvalue μmin = min{μk : k ∈ N}. Let q ∈ N be such that
μq = μmin . Assuming that the multiplicity of μmin is 2, thecenter manifold reduction performed for the order parameterin Ref. 5 yields the following stable branch bifurcating fromthe incoherent state (h ≡ 0) at K = Kc− :h∞ (x, K) =√
κ ±i2π q(x+φ)
e+ o( κ),β0 < κ := Kc− − K  1,

β := −

8g (0),π 3 g(0)4(28)

where φ ∈ [0, 1) depends on the initial data.
Equation (28) implies that at K = Kc− , the KM on smallworld graphs undergoes a pitchfork bifurcation. Unlike thebifurcation at K = Kc+ considered earlier, in the present case,the center manifold is two-dimensional. It is spanned by
w±q = e±i2π qx . In the remainder of this section, we analyze
stable spatial patterns bifurcating from the incoherent state at
K = Kc− . To this end, we rewrite the velocity field (17) using
the order parameter
V (t, u, ω, x) = ω +K −iu
e h(t, x) + eiu h(t, x) .
2i

(29)

The velocity field in the stationary regime is then given by
V∞ (u, ω, x) = ω +K −iue h∞ (x, K) + eiu h∞ (x, K) . (30)2iUsing the polar form of the order parameterh∞ (x, Kc − κ) = R∞ (x, K)ei(x) ,R∞ (x, K) =√κ+ o( κ),β(31)

from (30), we obtainV∞ (u, ω, x) = ω − (Kc− − κ)√+ o( κ).

κ
sin[u − 2π q(x + φ)]β(32)

Since ρ∞ is a steady state solution of (15), it satisfies∂u {V∞ ρ∞ } = 0.

(33)

073109-9

Chiba, Medvedev, and Mizuhara(a)Chaos 28, 073109 (2018)

(b)

(c)

FIG. 10. The formation of 2-twisted states for decreasing values of coupling strength K. From (a) to (c), we have, respectively, K = −32, K = −36, andK = −50. In all simulations p = 0.2 and r = 0.3.

From (32) and (33), we have
ρ∞ (u, ω, x)⎧δ [u ∓ 2π q(x + φ)
⎪
⎪
⎪
 √ 
⎨
√
−
arcsinω K √βκ + o( κ) ,
=
⎪
√
⎪
⎪
⎩1
ω2 −K 2 R∞ (x,K)2,
2π |ω−KR∞ (x,K) sin[u−(x)]||ω| ≤ KR∞ (x, K),|ω| > KR∞ (x, K).
(34)

Here, K = Kc− − κ and δ stands for the Dirac delta function.
Equation (34) describes partially phase locked solutions. They
combine phase locked oscillators (type I)√
u = ±2π q(x + φ) + Y (ω, κ) + O κ , √ βY (ω, κ) = arcsin ω √
,
(35)
K κ
and those whose distribution is given by the second line in(34) (type II). The oscillators of the first type form q-twisted
states subject to noise Y . Note that Y is a function of ωand, therefore, is random. Just near the bifurcation, where0 < κ  1, the probability of |ω| > KR∞ (x, K) is high and,thus, most oscillators are of the second type [Figs. 10(a)].
However, as we move away from the bifurcations, we see thenumber of oscillators of the first type increases and the noisy
twisted states become progressively more distinct [Figs. 10(b)and 10(c)]. Thus, the bifurcation analysis of the KM on smallworld graphs identifies a family of stable twisted states (cf.
Refs. 17 and 22). By changing parameters of small-world connectivity, one can control the winding number of the emerging
twisted states [see Fig. 4(c)].

and Paley graphs show more variability. The almost identical bifurcation plots for these models are due to the fact that
all three models result in the same mean field equation. Thismeans that the empirical measures generated by the trajectories of these models are asymptotically the same limit as
n → ∞. The differences seen in plots (a)–(c) of Fig. 6 are dueto the higher order moments, which are not captured by themean field limit. Other families of graphs considered in thispaper include the bipartite graph and a family of stochastic
block graphs interpolating between a graph with two disconnected components and Erdős-Rényi graph. The KM on allthese graphs (including the disconnected graph) feature thesame transition to synchronization. Finally, we studied thebifurcations in the KM on small-world and power law graphs
due to their importance in applications. A remarkable featureof the KM on small-world graphs is the presence of the bifurcation leading to stable noisy twisted states. Twisted states are
known as attractors in repulsively coupled KM with identical
intrinsic frequencies.16–18,22 For the KM with random frequencies from the Cauchy distribution, the transition from theVI. DISCUSSION

In this paper, we selected several representative familiesof graphs to illustrate the link between network topology andsynchronization and pattern formation in the KM of coupledphase oscillators. In particular, we showed that the transition to synchronization in the KM on pseudorandom graphs(e.g., Erdős-Rényi, Paley, and complete graphs) starts at thesame critical value of the coupling strength and proceeds inpractically the same way. The bifurcation plots shown in Fig.
6 are very similar, although plots for the KM on Erdős-Rényi

FIG. 11. Comparison of GPU and CPU computation time on Intel(R) Corei7-7700 CPU with NVIDIA Quadro K1200 GPU showing dramatic speed-up
due to parallelization of computations.

073109-10

Chiba, Medvedev, and Mizuhara

incoherent state to random twisted states was shown in Ref. 19
using the Ott-Antonsen ansatz. In this paper, we identified
the bifurcation in the repulsively coupled KM on smallworld graphs with intrinsic frequencies from arbitrary absolutely continuous distribution with even sufficiently smoothdensity.
To perform numerical experiments presented in thispaper, we had to overcome several challenges. Verificationof the bifurcation scenarios in the KM required a large number of simulations of large systems of ordinary differential
equations with random coefficients. Thus, the speed of computations was critical in this project. All computations werecompleted in MATLAB utilizing GPU computations for dramatic speed-up in the computational time (compared to CPU
computations), see Fig. 11. Time steps were performed using
Heun’s method with t = 0.01, which was sufficient for
stable simulation results.
Each simulation was initialized with a random statevector (un,1 , . . . , un,n )T with each component independently
chosen from a uniform distribution on [0, 2π ) representing the phase of each oscillator. We additionally initialize
a random vector of internal frequencies chosen from a normal distribution. Unless otherwise noted, in all simulations
we take n = 4001, which is a prime 1 modulo 4, so thatthe theory developed for Paley graphs applies. We observedthat taking Tfinal = 20 was sufficient for systems to exhibit
synchronization or q-twisted states.
Computationally solving (26) for the minimal negative eigenvalue is accomplished by observing the trivand iteratively computing μmin = μk0
ial bound μk > −1kπ
for k0 ∈ {1, . . . , M } for increasing values of M until
μmin ≤ M−1π .
ACKNOWLEDGMENTS

The work of the second author was supported in part by
the National Science Foundation Division of Mathematical
Sciences (NSF DMS) No. 1715161. The authors would liketo thank Nicholas Battista for helpful discussions related to
implementation of the numerical simulations.
1

N. Alon and J. H. Spencer, The Probabilistic Method, 4th ed., Wiley
Series in Discrete Mathematics and Optimization (John Wiley & Sons, Inc.,
Hoboken, NJ, 2016).

Chaos 28, 073109 (2018)2

C. Borgs, J. T. Chayes, H. Cohn, and Y. Zhao, “An Lp theory of sparsegraph convergence. I: Limits, sparse random graph models, and power law
distributions,” arXiv e-print (2014).
3
H. Chiba, “A proof of the Kuramoto conjecture for a bifurcation structureof the infinite-dimensional Kuramoto model,” Ergodic Theory Dyn. Syst.
35(3), 762–834 (2015).
4
H. Chiba and G. S. Medvedev, “The mean field analysis of the Kuramotomodel on graphs I. The mean field equation and the transition point formulas” Discrete and Continuous Dynamical Systems—Series A (to be
published).
5
H. Chiba and G. S. Medvedev, “The mean field analysis of the Kuramotomodel on graphs II. Asymptotic stability of the incoherent state, center
manifold reduction, and bifurcations” Discrete and Continuous Dynamical
Systems—Series A (submitted).
6
H. Dietert, “Stability and bifurcation for the Kuramoto model,” J. Math.
Pures Appl. 105(4), 451–489 (2016).
7
B. Fernandez, D. Gérard-Varet, and G. Giacomin, “Landau damping in the
Kuramoto model,” Ann. Henri Poincaré 17(7), 1793–1823 (2016).
8
D. Kaliuzhnyi-Verbovetskyi and G. S. Medvedev, “The mean field equation
for the Kuramoto model on graph sequences with non-Lipschitz limit,”
SIAM J. Math. Anal. 50(3), 2441–2465 (2018).
9
D. Kaliuzhnyi-Verbovetskyi and G. S. Medvedev, “The semilinear heat
equation on sparse random graphs,” SIAM J. Math. Anal. 49(2),1333–1355 (2017).
10
M. Krebs and A. Shaheen, Expander Families and Cayley Graphs (Oxford
University Press, Oxford, 2011) (A beginner’s guide).
11
L. Lovász, Large Networks and Graph Limits (AMS, Providence, RI,2012).
12
L. Lovász and B. Szegedy, “Limits of dense graph sequences,” J. Combin.
Theory Ser. B 96(6), 933–957 (2006).
13
G. S. Medvedev and X. Tang, “The Kuramoto model on power law graphs,”
arXiv e-prints (2017).
14
G. S. Medvedev, “The nonlinear heat equation on dense graphs and graphlimits,” SIAM J. Math. Anal. 46(4), 2743–2766 (2014).
15
G. S. Medvedev, “The nonlinear heat equation on W-random graphs,” Arch.
Ration. Mech. Anal. 212(3), 781–803 (2014).
16
G. S. Medvedev, “Small-world networks of Kuramoto oscillators,” Physica
D 266, 13–22 (2014).
17
G. S. Medvedev and X. Tang, “Stability of twisted states in the Kuramotomodel on Cayley and random graphs,” J. Nonlinear Sci. 25, 1169–1208(2015).
18
G. S. Medvedev and J. Wright, “Stability of twisted states in the continuum Kuramoto model,” SIAM J. Appl. Dyn. Syst. 16(1), 188–203(2017).
19
O. E. Omel’chenko, M. Wolfrum, and C. R. Laing, “Partially coherent
twisted states in arrays of coupled phase oscillators,” Chaos 24(2), 023102(2014).
20
S. H. Strogatz, “From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators,” Physica D 143(1-4),
1–20 (2000) (Bifurcations, patterns and symmetry).
21
D. J. Watts and S. H. Strogatz, “Collective dynamics of small-world
networks,” Nature 393, 440–442 (1998).
22
D. A. Wiley, S. H. Strogatz, and M. Girvan, “The size of the sync basin,”
Chaos 16(1), 015103 (2006).