Solvable model of the collective motion ofheterogeneous particles interacting on a sphere
Takuma Tanaka
Department of Computational Intelligence and Systems Science, Interdisciplinary
Graduate School of Science and Engineering, Tokyo Institute of Technology,
Yokohama, Japan
E-mail: tanaka.takuma@gmail.com
Abstract. I propose a model of mutually interacting particles on an M -dimensional
unit sphere. I derive the dynamics of the particles by extending the dynamics ofthe Kuramoto–Sakaguchi model. The dynamics include a natural-frequency matrix,
which determines the motion of a particle with no external force, and an external
force vector. The position (state variable) of a particle at a given time is obtainedby the projection transformation of the initial position of the particle. The same
projection transformation gives the position of the particles with the same naturalfrequency matrix. I show that the motion of the centre of mass of an infinite numberof heterogeneous particles whose natural-frequency matrices are obtained from a classof multivariate Lorentz distribution is given by an M -dimensional ordinary diﬀerential
equation in closed form. This result is an extension of the Ott–Antonsen theory.

PACS numbers: 05.45.Xt,05.65.+b

Submitted to: New J. Phys.

1. Introduction
Systems consisting of many elements, such as a society, a school of fish, a flock of birds
and the network of neurons in the brain, have attracted much research attention overyears [1, 2, 3]. Oscillation is a ubiquitous phenomenon found in systems such as thesethat consist of a large number of mutually interacting elements. Other examples of thesesystems include reaction–diﬀusion systems [4], electrochemical reactions [5], electronic
circuits [6, 7], spiking neurons [8], the human gait [9], flashing of fireflies [10], the female
menstrual cycle [11], a linear array of CO2 waveguide lasers [12] and Josephson-junction
arrays [13]. Dynamics of oscillatory element i can be described byẋi = χi (xi ),where xi is the state vector and χi (xi ) is the function defining the behaviour of element
i. The dynamics of the system consisting of these elements can be described by∑ẋi = χi (xi ) +ψij (xi , xj ),1≤j≤NCollective motion of particles on a sphere

2

where ψij (xi , xj ) is the function characterising the interaction from j to i and N is
the number of oscillators in the system. In many cases, these systems are so diﬃcult
to solve that one can only perform numerical simulations and observe the results. To
gain insight into the collective phenomena of these complicated systems, we need todevelop methods to analyse these systems. The phase description or phase reductionof limit-cycle oscillators introduced by Kuramoto [14] is the most widely used method
to analyse the synchronisation phenomena that lead to observed oscillations. The phase
description describes the state of a limit-cycle oscillator by a variable ϕ (called thephase) whose dynamics areϕ̇ = ω + Z(ϕ) · p(t),where ω is the natural frequency, p(t) is an external force and Z(ϕ) is the phasesensitivity function. If there is no external force, this system oscillates with constantfrequency ω. The external force p(t) advances or delays the phase. Z(t) determines thesensitivity of the oscillator to the external force. If two oscillators with phases ϕ1 andϕ2 are coupled to each other, the force exerted by oscillator 2 on oscillator 1 is givenby the phase ϕ2 , i.e. p(ϕ2 ). Assuming the coupling is weak, the long-time average ofthe mutual coupling Z(ϕ1 ) · p(ϕ2 ) can be regarded as a function of the phase diﬀerence
ϕ2 − ϕ1 . Thus, the dynamics of two oscillators are described byϕ̇1 = ω1 + Γ12 (ϕ1 − ϕ2 ),ϕ̇2 = ω2 + Γ21 (ϕ2 − ϕ1 ),where Γij is the coupling function characterising the interaction from j to i.
The Kuramoto–Sakaguchi model is the simplest model of weakly coupled oscillatorsand is described byK ∑ϕ̇i = ωi +sin(ϕj − ϕi − α),N 1≤j≤N
where K is the coupling strength (figure 1 (a)). In the Kuramoto–Sakaguchi model, thecoupling function is identical among all oscillator pairs. In such a system, the natural
frequency ωi diﬀers from oscillator to oscillator; the most intensively studied is thesystem in which the natural frequencies follow a Lorentz distribution. The system has acritical coupling strength Kc at which the synchronisation transition occurs. If K < Kc ,the system remains desynchronised and the order parameter r (i.e. the centre of mass)defined by
1 ∑r=exp(iϕi )N 1≤i≤Nis zero. For K > Kc , the system is synchronised and has non-zero r. This model has
been studied intensively and has contributed to the understanding of the synchronisation
phenomena in general [15].
Recently, two powerful methods to analyse the Kuramoto–Sakaguchi model haveattracted attention. First, Watanabe and Strogatz [16] showed that the trajectory of aCollective motion of particles on a sphere

(a)

3

(b)

Figure 1. Kuramoto–Sakaguchi model (a) and proposed model (b). The Kuramoto–
Sakaguchi model may be regarded as the dynamical model of N particles interacting
on a circle. The proposed model describes the dynamics of N particles interacting on∑
an M -dimensional sphere. Here we set g = K
j xj .
N

system of N -coupled oscillators is on three-dimensional invariant manifold if the system
is described byK ∑ϕ̇i = ω +sin(ϕj − ϕi + α) + J sin(ϕi + β),(1)
N 1≤j≤N
where ω is the natural frequency common to all oscillators; K, J, α, and β aretime-dependent parameters. The phase of a single oscillator in the system at time
t is determined by its initial phase and a nonlinear transformation involving threeparameters. An important point is that the transformation with the same parameters
can be used to obtain the phase of any oscillator in the system. It has been shownthat the Möbius transformation describes the i-th phase at time t [17, 18]. Thistransformation is a one-variable projection transformation in the complex planezi (t) = exp(iσ)

a − zi (0),1 − āzi (0)(2)

where zi (t) = exp[iϕi (t)] and both of σ and a are time-dependent and oscillatorindependent parameters. This type of Möbius transformation maps the unit circle on
the complex plane to itself. Because the Möbius transformation (2) has a real parameterσ and a complex parameter a, all the trajectories of the system are confined to manifolds
with three dimensions. In other words, the dynamics of the system can be described by
the dynamics of σ and a. This method simplifies the analysis of a system of oscillatorswith identical natural frequencies and enables a precise analysis of a wider range ofsynchronisation phenomena such as chimera states [19] and persistent fluctuations insynchronisation rates [20] than before.

Collective motion of particles on a sphere

4

Second, consider Kuramoto–Sakaguchi oscillators whose natural frequencies comefrom a Lorentz distribution1γp(ω) =,2π ω + γ2where γ is the scale parameter. Ott and Antonsen showed that the order parameter rof an infinite number of such oscillators obeys the Stuart–Landau equation
1
Kṙ = (K cos α − 2γ − iK sin α)r − exp(iα)|r|2 r2
2
under the condition that the initial phases are on a specified manifold [21]. Ott–
Antonsen theory has been used to describe systems with bimodal natural-frequency
distributions [22] and systems with two subpopulations [23]. In these studies, thesystems contain an infinite number of heterogeneous oscillators, but are describedby diﬀerential equations with a few dynamical variables, allowing the systems to beanalysed. Ott–Antonsen theory made the analysis of the synchronisation phenomena ofcomplicated systems much simpler and much more thorough than previous methods.
The Kuramoto model with α = 0 and ω = 0 can be regarded as the XY spinmodel if noise is injected into the oscillators. The Hamiltonian of the XY spin model∑is H = − ⟨i,j⟩ Jij cos(ϕj − ϕi ), where the summation runs over all connected i and j.
Spins in the XY spin model are two-dimensional unit vectors. The XY spin model isgeneralised to the n-vector model, in which spins are unit vectors of arbitrary dimension.
The XY spin model and n-vector model have been used to analyse frustrated magnets.
The n-vector models are highly simplified and abstract models for describing magnetic
spins; they do not necessarily correspond to physical substances. However, solving then-vector models has advanced our understanding of phase transition. Moreover, then-vector model with n = 0 was suggested to correspond to self-avoiding walks [24, 25].
Similarly, although percolation models in high-dimensional lattices and random graphs
do not necessarily correspond to physical materials, the results of infinite-dimensional
percolation were useful to interpret the results of percolation on complex networks suchas the attack tolerance of networks [26]. The generalisation of such models allows us topredict what happens in novel problems and to find relationships between two seemingly
unrelated physical systems.
The n-vector model, which is a high-dimensional extension of the XY spinmodel, has contributed to statistical physics by improving our understanding of phasetransitions. A high-dimensional extension of the Kuramoto model may similarly advancethe understanding of the collective phenomena of coupled systems and can be applied toother fields. If the extended model is solvable, it would have many possible applications.
In particular, if the methods similar to those of Watanabe and Strogatz and Ott
and Antonsen are applicable to the model, it would be quite useful for understanding
collective phenomena.
The study of collective phenomena goes back a long time. In the 17th century,
Huygens observed the antiphase locking of two pendulum clocks [27]; in the 18th
century, Kaempfer reported the synchronised flashing of fireflies at the banks of theCollective motion of particles on a sphere

5Chao Phraya [10]. Although these descriptions date back to the early modern period,
the study of synchronisation phenomena advanced only in the latter half of the 20th
century with the development of the phase-description method of limit-cycle oscillators[14, 28]. Similarly, a solvable generalisation of the Kuramoto–Sakaguchi model would
be useful for understanding the collective phenomena in a broader setting.
Motivated by this historical perspective, I propose an extension of the Kuramoto–
Sakaguchi model to describe the collective motion of particles interacting on a unit
sphere (figure 1 (b)). First, by extending the framework of the Watanabe–Strogatz
transformation, we derive the dynamics of the individual elements, which I call‘particles’. The present states of these particles, which I sometimes call ‘positions’,
are provided by the projection transformation (a high-dimensional extension of Möbiustransformation) of their initial values. This dynamics of a particle are given byẋ = − (g † x)x + Ωx + g= − xg † x + Ωx + g,
where † denotes the Hermitian conjugate, x is an M -dimensional real or complex
vector representing the state of the particle, Ω is an M × M matrix corresponding
to the natural frequency of phase oscillators and g determines the force exerted on
the particle. Assuming that x is a two-dimensional real vector and assuming thatx is a one-dimensional complex vector lead to the same dynamics as the Kuramoto–
Sakaguchi model. The position of the particle at time t is given by the same projection
transformation, irrespective of the initial position x0 . Second, we derive the centreof mass r of the particles that initially are uniformly distributed on the sphere. The
centre of mass corresponds to the order parameter of the Kuramoto–Sakaguchi model.
In particular, I show that the centre of mass is given by one of the parameter vectors of
the projection transformation if the state vector x is a complex vector. Third, we extend
the Ott–Antonsen theory to high-dimensional systems. By assuming that the naturalfrequency matrices of particles are obtained from a multivariate Lorentz distribution, I
show that the centre of mass of particles with complex variables is described by lowdimensional ordinary diﬀerential equations. In particular, I show that, if the particles
in a system are attracted to Kr, where K is an M × M matrix, the dynamics of r aredescribed byṙ = −rr † K † r + Ω̂r + Kr,
where Ω̂ is a matrix determined by the probability distribution of the natural frequenciesof the particles. We derive the limit cycle from these dynamics and show that a transition
emerges from a desynchronised state to a synchronised state that is similar to that foundin the Kuramoto model. I also show that this theoretically derived limit cycle agrees
well with the results of numerical simulation. Finally, I discuss the possible applicationsof the present model and problems to be addressed in the future.

Collective motion of particles on a sphere

6

2. Results and discussion
2.1. Dynamics of the particles on the unit sphere induced by the projection
transformation
The Möbius transformation (2) has been shown to underlie the dynamics of (1) [18].
Marvel and collaborators showed that the phase of oscillator i at time t is given by theMöbius transformation (2) of its phase at time 0. Notably, given the initial conditions ofoscillators, the same parameters σ and a (i.e. the same Möbius transformation) can beused to calculate the phases at t of all oscillators in the system. To extend the Kuramoto–
Sakaguchi model to an M -dimensional system, we use the Möbius transformation.
Because this transformation is a one-dimensional projection transformation, we derive
the dynamics of variables whose present values are given by an M -dimensional projection
transformation of the initial values. We consider a system in which the M -dimensional
state vector x of a particle is given by the projection transformation
Ax0 + bx= †(3)
c x0 + dof the initial state x0 , where the M × M matrix A, M -dimensional vectors b and c andscalar value d are time dependent. The system can be real valued or complex valued(i.e. x can be a real or a complex vector). To derive the dynamics of particles whosetime evolution is described by (3), we diﬀerentiate x with respect to time to obtain
˙(Ȧx0 + ḃ)(c† x0 + d) − (Ax0 + b)(ċ† x0 + d)
ẋ =
.
(4)
(c† x0 + d)2
Substituting−1x0 = (A − xc† ) (dx − b)1=
(A−1 − c† A−1 xA−1 + A−1 xc† A−1 )(dx − b)
1 − c† A−1 x
into Qx0 + q, we have
Qx0 + q = ι−1 [Q(A−1 − c† A−1 xA−1 + A−1 xc† A−1 )(dx − b) + ιq]= ι−1 [QA−1 (dx − b) + Q(c† A−1 xA−1 − A−1 xc† A−1 )b + ιq],
where ι = 1 − c† A−1 x. Replacing Q and q with c† and d, respectively, yieldsc† x0 + d = κι−1 ,where κ = d − c† A−1 b. The first term in the numerator of (4) isκι−2 [ȦA−1 (dx − b) + Ȧ(c† A−1 xA−1 − A−1 xc† A−1 )b + ιḃ],and the second term in the numerator of (4) is˙κι−2 x[ċ† A−1 (dx − b) + ċ† (c† A−1 xA−1 − A−1 xc† A−1 )b + ιd].
Rearranging terms gives˙ † )A−1 xẋ = [x(−κċ† − ċ† A−1 bc† + dc
−1˙
+ (κȦ + ȦA−1 bc† − ḃc† + ċ† A−1 bA − dA)Ax

− ȦA−1 b + ḃ]/κ.

(5)

Collective motion of particles on a sphere

7

The dynamicsẋ = xh† x + Ωx + g(6)can be obtained from (5) if
˙ † )A−1 ,h† = κ−1 (−κċ† − ċ† A−1 bc† + dc
−1˙
Ω = κ−1 (κȦ + ȦA−1 bc† − ḃc† + ċ† A−1 bA − dA)A,g = κ−1 (−ȦA−1 b + ḃ).
The dynamics of the parameters of the projection transformation to realise (6) are givenbyd˙ = 0,(7a)† −1ċ† = κh† A(−κI − A−1 bc )−1

= − h† A(A + κ−1 bc† ) A)(
κ−1 A−1 bc† A−1†−1
A= −h A A −1 + κ−1 c† A−1 b= − h† A(A−1 − d−1 A−1 bc† A−1 )A= − h† (A − d−1 bc† ),(7b)Ȧ = κ−1 (κΩA − ȦA−1 bc† + ḃc† − ċ† A−1 bA)= κ−1 [κΩA + (−ȦA−1 b + ḃ)c† − ċ† A−1 bA]= κ−1 (κΩA + κgc† + h† bA − d−1 h† bc† A−1 bA)= ΩA + gc† + d−1 h† bA,(7c)ḃ = κg + ȦA−1 b= κg + Ωb + gc† A−1 b + d−1 h† bb= dg + Ωb + d−1 h† bb.

(7d)

Thus, there are two ways to obtain the state x at time t. First, to directly calculate (6)from the initial condition x0 . Second, to calculate the parameters of the transformation,
A, b, c and d, using (7a)–(7d) and the initial condition A = I, b = 0, c = 0 and
d = 1 and then map x0 to x using (3). Note that, with these initial conditions forthe parameters, A, b, c and d, the relationship x = x0 is satisfied at t = 0. The
first and the second method calculate the time evolution of the state vector and thetransformation, respectively. Thus, one can regard these methods as corresponding to
the Schödinger and Heisenberg representation of quantum mechanics.
The dynamics of phase oscillators can be regarded as those restricted to a unit
circle in a plane. Similarly, we restrict the dynamics of a particle to the sphere |x| = 1.
The conditiond|x|2= ẋ† x + x† ẋ
dt= (h† + g † )x + x† (Ω + Ω† )x + x† (h + g)
= 0,

Collective motion of particles on a sphere8

is satisfied by the constraints
g = − h,(8a)†Ω= −Ω .

(8b)

Thus, the dynamics are described byẋ = −xg † x + Ωx + g,

(9)where g is an arbitrary vector and Ω is an antisymmetric matrix in the real-valued
system and an anti-Hermitian matrix in the complex-valued system. This is a subclassof the diﬀerential equations called the Riccati matrix diﬀerential equation [29, 17].
For a real two-dimensional vector
[
]
x1
,
x=
x2
(9) reduces toẋ1 = − (g1 x1 + g2 x2 )x1 − ωx2 + g1 ,ẋ2 = − (g1 x1 + g2 x2 )x2 + ωx1 + g2 ,where we set

[
g=

and

]
g1
g2

(
Ω=

0 −ω
ω 0

)
.

Replacing x1 with cos ϕ and x2 with sin ϕ gives−ϕ̇ sin ϕ = − (g1 cos ϕ + g2 sin ϕ) cos ϕ − ω sin ϕ + g1
= − sin ϕ(ω + g2 cos ϕ − g1 sin ϕ)= − sin ϕ[ω + K sin(ψ − ϕ)],where we assume g1 = K cos ψ and g2 = K sin ϕ, which is the dynamics of a Kuramoto
oscillator under an external force. Similarly, for a complex variable x = [exp(iϕ)], (9)reduces toϕ̇ = iḡ exp(iϕ) + ω − ig exp(−iϕ)= ω + 2 Im[g exp(−iϕ)]= ω + K sin(ψ − ϕ),where Im denotes the imaginary part and we set Ω = (iω) and g = [g] = [K exp(iψ)/2],
which is the dynamics of a Kuramoto oscillator as well. Thus, these two systems have thesame dynamics. However, in general, a system of 2M -dimensional real state variables
and a system of M -dimensional complex state variables do not have the same dynamics
because the antisymmetric matrix Ω has 2M 2 − M real degrees of freedom in the former

Collective motion of particles on a sphere9case, whereas the anti-Hermitian matrix Ω has M 2 real degrees of freedom in the lattercase. These two examples suggest that g and Ω correspond to the external force andthe natural frequency, respectively, in the Kuramoto–Sakaguchi model. The naturalfrequency matrix Ω determines the direction and the speed of the rotation of the particle
on the unit sphere.
Next, we consider whether the interaction terms −xg † x + g of (9) can be regarded
as the tangential component of a central force on the surface of the unit sphere. Let us
assume that the central force is described byg−x
f (|g − x|),|g − x|whose origin is g. Assuming that x† x = 1, the tangential component is calculated byusing the projection operator I − xx† as(I − xx† )f (|g − x|)g−x
g − xx† g= f (|g − x|),|g − x||g − x|which gives −xg † x + g if x and g are real vectors and f (|g − x|) = |g − x|. Thus, theinteraction of particles in this system is the same as that for the system of objects thatare attracted to g by ideal springs and are restricted to lie on the M -dimensional unit
sphere.
Potential energy is obtained by integrating the force −xg † x + g. Assuming thatparticle i is attracted to particle j (i.e. g = xj ) the force exerted on particle i by particle
j is given by −xi x†j xi + xj . The potential energy is given by
∫ xi
E= −
(−xx†j x + xj ) · dxx∫ xj i= −xj · dx
xj= xj · xj − xj · xi
= − xj · xi + const.,
where we notice that the position vector x on the sphere is orthogonal to the tangent
vector dx. E is the potential energy of the n-vector model.
2.2. Dynamics of the parameters of the projection transformation
Constrained by (8a) and (8b), the dynamics of the parameters areȦ = ΩA + gc† − g † bA,

(10a)

†

ḃ = − bg b + Ωb + g,
†

†

(10b)

†

ċ = g (A − bc ),(10c)where we set d = 1. In this system, the parameters A, b and c of the projection
transformation are dependent variables. The projection transformation to and from theunit sphere satisfies†(Ax0 + b)† (Ax0 + b) = (c† x0 + 1) (c† x0 + 1)

Collective motion of particles on a sphere

10

under the constraint x0 † x0 = 1, i.e.
x0 † (A† A − cc† + |b|2 I − I)x0 + x0 † (A† b − c) + (b† A − c† )x0 = 0,
therefore,c†= b† A,
†

†

(11a)

A A = cc + (1 − |b| )I.
2

(11b)

From (11a), c is determined by A and b. Substituting (11a) into (11b), we obtainA† A = A† bb† A + (1 − |b|2 )I.
Rearranging the terms and multiplying by A†
right, we have−1from the left and by A−1 from the−1

I − bb† = (1 − |b|2 )(AA† )therefore,−1AA† = (1 − |b|2 )(I − bb† )
(
)
bb†
2
= (1 − |b| ) I +1 − b† b= (1 − |b|2 )I + bb† .
The matrix A satisfying this equation is defined by the polar decomposition
A = H 1/2 U ,where U is an arbitrary orthogonal matrix if x is a real vector and an arbitrary unitarymatrix if x is a complex vector andH 1/2 = V Σ1/2 V † ,√with Σ1/2 = ( Σij ) and V ΣV † being the singular value decomposition of
H = (1 − |b|2 )I + bb† .
V appears two times in the singular value decomposition of H because H is a Hermitianmatrix. Taken together, these relationships imply that the projection transformation
to and from the unit sphere is determined by a vector b and an orthogonal or unitarymatrix U . Because the degrees of freedom of an M -dimensional orthogonal matrix andan M -dimensional unitary matrix are M (M − 1)/2 and M 2 , respectively, the degree of
freedom of the transformation is M (M + 1)/2 in the real-valued system and M 2 + 2M
in the complex-valued system. Thus, the real-valued systems with M = 2 have threedegrees of freedom, and the complex-valued systems with M = 1 also have three degrees
of freedom. It is identical to the dimension of the invariant manifolds of Watanabe andStrogatz [16]. The vector b is an eigenvector of H with eigenvalue 1 because
Hb = (1 − |b|2 )b + bb† b= b.

Collective motion of particles on a sphere

11

Vectors orthogonal to b (i.e. vectors n satisfying b† n = 0) are eigenvectors of H witheigenvalues 1 − |b|2 because
Hn = (1 − |b|2 )n + bb† n= (1 − |b|2 )n.
From these observations, we obtain
H 1/2 b = b
and
H 1/2 n =(12)√
1 − |b|2 n.

(13)

H 1/2 approaches I in the limit |b| → 0, which indicates that the initial conditions
A = I and b = 0 are consistent with each other if we set U = I.
2.3. Centre of mass of particles
We set the initial conditions of the parameters of the transformation to A = I andb = 0 to satisfy x = x0 . The projection transformation of x0 to x is given as
H 1/2 U x0 + b
x = † 1/2
b H U x0 + 1
H 1/2 U x0 + b=.
b† U x0 + 1
We now derive the centre of mass
1 ∑
m=
xi
N 1≤i≤N(14)of the particle. We assume that, initially, the points are uniformly distributed on theunit sphere. This assumption simplifies the calculation of the centre of mass and enables
the low-dimensional description of the entire system. Assuming that the points are
uniformly distributed on the unit sphere in the initial conditions and that there are an
infinite number of particles, we have
∫
1
m=
x dx0
SM |x0 |=1
for real-valued systems and
∫
1
x dx0
m=
S2M |x0 |=1
M/2

2π
for complex-valued systems, where SM = Γ(M/2)
is the surface area of the M -dimensional
unit sphere. Equation (14) suggests that the density of points and the centre ofmass are completely determined by b because U does not aﬀect the distribution after
transformation if the particles are uniformly distributed on the unit sphere. Thus,
without loss of generality, we can replace U of (14) with I for calculating the centre ofmass m. Decomposing x0 intox0 = ηb1 + n,Collective motion of particles on a sphere

12

where b1 = |b|−1 b, η = b†1 x0 , and n is a vector orthogonal to b (i.e. n = x0 − ηb1 ) wehave
H 1/2 (ηb1 + n) + bx=b† (ηb1 + n) + 1√(η|b|−1 + 1)b + 1 − |b|2 n=,
η|b| + 1
where we use (12) and (13) and set U = I.
The centre of mass of the particles in real space is given by√∫(η|b|−1 + 1)b + 1 − |b|2 n1
m=dx0
SM |x0 |=1
η|b| + 1√∫ 1∫(η|b|−1 + 1)b + 1 − |b|2 n
1
1
dη
=
dn √
√
SM −1 |n|= 1−η2
η|b| + 1
1 − η2∫
SM −1 b 1 η|b|−1 + 1 (1 − η 2 )(M −2)/2
√
=
dη
SM
1 − η2−1 η|b| + 1∫
SM −1 b π |b|−1 cos θ + 1 M −2
=
sin
θ dθ
SM
|b| cos θ + 1
0()√ Γ[(M + 1)/2]
Γ(M/2)2bπ=√
2 F1 1/2, 1; (M + 2)/2; |b|
Γ[(M + 2)/2]πΓ[(M − 1)/2]
(
)
M −1
2
=
b,
(15)
2 F1 1/2, 1; (M + 2)/2; |b|M
where 2 F1 (a, b; c; x) is the ordinary hypergeometric function.
The centre of mass of the particles in complex space for M = 1 is given by√∫(η|b|−1 + 1)b + 1 − |b|2 n
1
m=dx0
S2 |x0 |=1
η|b| + 1∫ 2π −1|b| exp(iθ) + 1
b
dθ
=
2π 0
|b| exp(iθ) + 1= b,
because n = 0 in this case. For M ≥ 2, the centre of mass is given by√∫∫(η|b|−1 + 1)b + 1 − |b|2 n
1
1
m=
dn √
dη
√
S2M |η|≤1 |n|= 1−|η|2
η|b| + 1
1 − |η|2∫
S2M −2 b
η|b|−1 + 1 (1 − |η|2 )(2M −3)/2
√
dη
=
S2M
1 − |η|2|η|≤1 η|b| + 1∫ ∫S2M −2 b 1 2π |b|−1 r exp(iθ) + 1=(1 − r2 )M −2 r dθ dr
S2M|b|rexp(iθ)+100∫
S2M −2 b 1=
2π(1 − r2 )M −2 r dr
S2M
0
2πS2M −2 b=2(M − 1)S2M= b.
Thus, b is the centre of mass of the system with complex variables.

Collective motion of particles on a sphere13

2.4. Results of simulation for the system with identical particles
Figure 2 shows the results of simulation for particles described by a real threedimensional variable. In this system, all particles have the natural-frequency matrix
0 −1 1


Ω= 1
(16)
0 −1  ,−1 1
0and experience the external force


cos t


g =  sin t  .
0
No mutual coupling among particles is introduced in figure 2(a). The solid line infigure 2(a) is obtained by the direct simulation of (9) and shows the trajectory of one ofthe particles. The dashed line shows the trajectory obtained by simulating the dynamics
of parameters A and b by (10a) and (10b) and transforming the initial conditions of
the particle by the projection transformation (3) with these parameters. These twosimulation methods yield the same result.
Adding a mutual interaction and setting the external force to


cos t


g = Km +  sin t  ,(17)0where m is the centre of mass of N = 1000 particles and−1 0 11K =  0 −1 0  ,2
1
0 1
I perform the direct simulation with N = 1000 particles, which initially are randomly
distributed on the unit sphere. The same natural-frequency matrix as (16) is used forall particles. The solid line in figure 2(b) shows the trajectory of the centre of massof the particles in the system. Approximating m of (17) by (15) (Figure 2(b), dashedline), I simulate the dynamics of parameters A and b by (10a) and (10b) and obtainthe trajectory with the projection transformation (3). Figure 2(b) shows that the direct
simulation of the original system agrees well with this approximation. This agreement
occurs because the initial distribution of particles may be regarded as uniform when thenumber of particles is suﬃciently large.
2.5. Extension of the results of Ott and Antonsen
Next, we consider a heterogeneous system in which the natural-frequency matrix Ωvaries from particle to particle. Ott and Antonsen showed that the order parameter(i.e. the centre of mass) of an infinite number of Kuramoto–Sakaguchi oscillators obeys

Collective motion of particles on a sphere(a)Direct simulation
1
Watanabe–Strogatz
0.8
0.6
0.4
0.2
x3 0
–0.2
–0.4
–0.6
–0.8
–0.6
x2 0
0.6
1.2 –0.4 –0.2
0.2
0
x1

14

(b)

Direct simulation
1
Watanabe–Strogatz
0.8
0.6
0.4
0.2
x3 0
–0.2
–0.4
–0.6
–0.8
–0.4
x2 0
0.4
0.8 –0.6–0.4–0.2 0 0.2 0.4 0.6 0.8
0.4
x1

Figure 2. Dynamics of the system of identical particles. (a) Direct simulation of (9)(solid line) and the simulation by using (10a) and (10b) of a particle under an external
force (dashed line). (b) Direct simulation (solid line) and approximation of the centreof mass (dashed line) for the system with N = 1000 particles. Simulation time is
T = 300 for both (a) and (b).

the Stuart–Landau equation if their natural frequencies are obtained from a Lorentz
distribution [21]. To extend their results, we introduce two additional assumptions.
First, we assume a complex-valued system and that the state x of the particle is acomplex vector. Because the centre of mass m is approximated by b in a complexvalued system, this assumption facilitates the derivation of the dynamics of the centreof mass. Second, we assume that the natural-frequency matrix Ω is drawn from amultivariate Lorentz distribution. We assume that the entries of the natural-frequency
matrix Ω are given by the linear superposition of random variables drawn from Lorentz
distributions:∑Ω = Ω0 +ζz Ωz ,
1≤z≤Z

where ζz is obtained from the Lorentz distribution with location parameter µz and scaleparameter γz > 0, whose probability density is given by1γzpz (ζz ) =.
π (ζz − µz )2 + γz2
Because Ω must be an anti-Hermitian matrix, the eigenvalues of Ωz (0 ≤ z ≤ Z) arepurely imaginary. We assume that
∫ RZ
g(t) = K(t)m(Ω, t)p(ζ) dζ + f (t)= K(t)r(t) + f (t),where K(t) is the coupling strength, f (t) is the external force, and p(ζ) =∏1≤z≤Z pz (ζz ). Here, we show that, if Ω0 is an anti-Hermitian matrix and the eigenvaluesof Ωz are positive and imaginary for z > 0, the dynamics of this heterogeneous system
are described byṙ(t) = −r(t)[K(t)r(t) + f (t)]† r(t) + Ω̂r(t) + K(t)r(t) + f (t),(18)

Collective motion of particles on a sphere

15

where
r(t) = m(Ω̂, t)andΩ̂ = Ω0 +∑(19)(µz + iγz )Ωz .

1≤z≤Z

First, we derive the analytic continuation of m(Ω, t), assuming that ζz -s arecomplex numbers. Because m can be replaced by b in complex-valued systems, we
obtain from (10b) the dynamics of m(Ω, t), which isṁ(Ω, t) = −m(Ω, t)g(t)† m(Ω, t) + Ωm(Ω, t) + g(t).
Assuming that m(Ω, t) satisfies the Cauchy–Riemann equations∂∂m(Ω, t) = −im(Ω, t),∂uz∂vz
where ζz = uz + ivz and uz and vz are real numbers, then ṁ(Ω, t) satisfies the Cauchy–
Riemann equations∂∂ṁ(Ω, t) = −iṁ(Ω, t)∂uz
∂vz
because)∂ d∂ (†m(Ω, t) =−m(Ω, t)g(t) m(Ω, t) + Ωm(Ω, t) + g(t)∂uz dt∂uz∂∂m(Ω, t)g(t)† m(Ω, t) − m(Ω, t)g(t)†m(Ω, t)
= −
∂uz
∂uz
∂
∂
+
Ωm(Ω, t) + Ωm(Ω, t)∂uz∂uz∂∂m(Ω, t)g(t)† m(Ω, t) + im(Ω, t)g(t)†m(Ω, t)
=i
∂vz
∂vz
∂
∂
−i
Ωm(Ω, t) − iΩm(Ω, t)
∂vz
∂vz
∂ d
= −i
m(Ω, t),∂vz dtwhere we used
(
)
∑
∂
∂
Ωij =
Ω0ij +(uz + ivz )Ωzij
∂uz
∂uz
1≤z≤Z
= Ωzij
∂
= −i
∂vz

(
Ω0ij +

∂
Ωij .
= −i
∂vz

∑
1≤z≤Z

)
(uz + ivz )Ωzij

Collective motion of particles on a sphere

16

These relations indicate that m(Ω, t) is an analytic function of ζz if m(Ω, t′ ) is ananalytic function of ζz , where t′ < t [30]. Here, we have
∫ RZ−1 ∫ R
r(t) =m(Ω, t)pk (ζk ) dζk p̂k (ζ̂k ) dζ̂k
∫ RZ−1 ∫ R=

m(Ω, t)1γkdζk p̂k (ζ̂k ) dζ̂k ,π (ζk − µk )2 + γk2

where
ζ̂k

= [ζ1 , . . . , ζk−1 , ζk+1 , . . . , ζz ],∏p̂k (ζ̂k ) =pz (ζz ).
1≤z≤Z, z̸=kBecause m(Ω, t) is an analytic function, we have
∫ R1γkm(Ω, t)dζk
π (ζk − µk )2 + γk2∫ S1
m(Ω, t)m(Ω, t)= lim−dζk
S→∞ 2πi −S ζk − µk − iγkζk − µk + iγk∫ Cm(Ω, t)1m(Ω, t)+ lim−dζk
S→∞ 2πiζk − µk − iγk ζk − µk + iγk= m(Ω̂k , t),(20)

where C is a semicircle in the upper half of the complex plane with radius S and centred
at the origin, and∑ζz Ωz ,
Ω̂k = Ω0 + Ωk (µk + iγk ) +1≤z≤Z, z̸=k

if the second integral on the right-hand side converges to zero. Because Ωz is an antiHermitian matrix, the time evolution of |m(Ω, t)|2 on |m(Ω, t)| = 1 satisfies()
1 d|m(Ω, t)|2 = Re m(Ω, t)† ṁ(Ω, t)
2 dt()= Re m(Ω, t)† Ωm(Ω, t)
(
(
)
)
∑
= Re m(Ω, t)† Ω0 +ζz Ωz m(Ω, t)
1≤z≤Z

(
†

= Re m(Ω, t)(∑)(uz + ivz )Ωz m(Ω, t)

1≤z≤Z

)
∑
= Re m(Ω, t)†ivz Ωz m(Ω, t) ,1≤z≤Z

where Re denotes the real part. If all eigenvalues of the Hermitian matrix iΩz are
negative, thend|m(Ω, t)|2 ≤ 0,dtCollective motion of particles on a sphere17

on the sphere |m(Ω, t)| = 1 when ζk is in the upper-half complex plane (i.e. vz ≥ 0).
Thus, m(Ω, t) remains finite if Ωz are anti-Hermitian matrices with positive imaginary
eigenvalues and Im ζk ≥ 0. In the limit of large S, the dynamics of m(Ω, t), where|ζk | = S, can be approximated byṁ(Ω, t) = ζk Ωk m(Ω, t),where we approximate Ω by ζk Ωk and ignore the terms without Ω. If we use the same
assumption that all eigenvalues of Ωk are positive and imaginary, then the real part ofall eigenvalues of ζk Ωk are negative under the condition Im ζk > 0. In this case, m(Ω, t)approaches to zero as S increases, and so the second integral of (20) converges to zeroin the limit of large S. Thus, we have
∫ RZ−1r(t) =m(Ω̂k , t)p̂k (ζ̂k ) dζ̂k .
Because m(Ω̂k , t) is an analytic function of ζz (z ̸= k), the integration can be donerecursively to obtain (19), so the dynamics of the order parameter r(t) are given by(18). Note that |m(Ω , t)| remains finite and converges to zero in the limit of large
S = |ζz | (z ̸= k) even if Ω is replaced by Ω̂k .
2.6. Limit-cycle oscillation of the centre of mass
Assuming that there is no external force (i.e. f (t) = 0) and that the mutual interactions
among oscillators are constant (i.e. K(t) = K), the dynamics are described byṙ(t) = −r(t)r(t)† K † r(t) + Ω̂r(t) + Kr(t),(21)

which is a high-dimensional extension of the Stuart–Landau equation. The stability ofthe desynchronised state, r(t) = 0, is determined by the eigenvalues of Ω̂ + K. BecauseΩ0 is an anti-Hermitian matrix, Ωz (z > 0) are anti-Hermitian matrices with positive
imaginary eigenvalues and γz > 0, so the real part of∑x† Ω̂x = x† Ω0 x +(µz + iγz )x† Ωz x
1≤z≤Z

is non-positive. Therefore, all real parts of eigenvectors of Ω̂ are non-positive, so thesystem without any mutual interaction (K = 0) remains desynchronised. If Ω̂ + K has
an eigenvector ei with eigenvalue λi , where Re λi > 0, the limit cycle is given byri = Ri exp(iξi t)ei ,where
iξi = −Ri2 e†i K † ei + λi
and |ei | = 1. Because(Ω̂ + K)ei = λi ei ,we gete†i (Ω̂† + K † )ei = λ̄i .

(22)Collective motion of particles on a sphere18

The real part of the right-hand side is positive because Re λi > 0. Re(e†i Ω̂† ei ) ≤ 0
because the real part of x† Ω̂x is non-positive for any x. Thus, we have Re(e†i K † ei ) > 0,and so (22) is satisfied with a positive Ri2 and a λi with a positive real component.
Assuming that λ1 is the eigenvalue with the largest real component, we prove that
only the limit cycle r1 is stable. Let us examine the time evolution of the perturbed
solutionwhere d =∑ri = (Ri + ρ) exp(iξi t + iθ)(ei + d),
1≤j≤M, j̸=i dj ej , |di | ≪ 1, |ρ| ≪ 1 and |θ| ≪ 1. Because we havedri= ρ̇ exp(iξi t + iθ)(ei + d) + (Ri + ρ) exp(iξi t + iθ)d˙dt+ (Ri + ρ)i(ξi + θ̇) exp(iξi t + iθ)(ei + d)
≈ [ρ̇ + iθ̇Ri + iρξi ] exp(iξi t)ei + Ri exp(iξi t)d˙
+ iRi ξi exp(iξi t)d + iRi ξi exp(iξi t + iθ)ei ,− ri ri† K † ri + Ω̂ri + Kri= − (Ri + ρ)3 exp(iξi t + iθ)(ei + d)(ei + d)† K † (ei + d)
+ (Ri + ρ) exp(iξi t + iθ)(Ω̂ + K)(ei + d)
≈ [−3ρRi2 ei e†i K † ei − Ri3 de†i K † ei − Ri3 ei d† K † ei
− Ri3 ei e†i K † d + ρ(Ω̂ + K)ei + Ri (Ω̂ + K)d] exp(iξi t)+ [−Ri3 ei e†i K † ei + Ri (Ω̂ + K)ei ] exp(iξi t + iθ)= [−3ρRi2 ei e†i K † ei − Ri3 de†i K † ei − Ri3 ei d† K † ei
− Ri3 ei e†i K † d + ρ(Ω̂ + K)ei + Ri (Ω̂ + K)d] exp(iξi t)+ iRi ξi exp(iξi t + iθ)ei ,∑replacing d with 1≤j≤M, j̸=i dj ej , we obtain∑[ρ̇ + iRi θ̇ + iξi ρ]ei + Ri
1≤j≤M, j̸=i

= 3(iξi − λi )ρei + Ri (iξi − λi )− Ri3 ei e†i K †∑1≤j≤M, j̸=i∑d˙j ej + iRi ξi∑d j ej1≤j≤M, j̸=i∑

dj ej − Ri3 ei

1≤j≤M, j̸=i

dj ej + λi ρei + Ri∑

d¯j e†j K † ei

1≤j≤M, j̸=i

λj dj ej ,1≤j≤M, j̸=iwhere we used Ri2 e†i K † ei = λi − iξi and (Ω̂ + K)ei = λei . Thus, the dynamics of dj is
given byd˙j = (−λi + λj )dj
if i ̸= j. d1 is unstable if i > 1 because Re λ1 > Re λi . Hence, the limit cycle ri is
unstable if i > 1. Conversely, dj for j > 1 is stable if i = 1. This fact allows us to set
d = 0 in examining the stability of the limit cycle r1 . The resulting dynamics of ρ andθρ̇ + iR1 θ̇ = 2(iξ1 − λ1 )ρ

Collective motion of particles on a sphere19

reduce toρ̇= − 2 Re λ1 ρ,

R1 θ̇ = 2(ξ1 − Im λ1 )ρ.
Hence, ρ has a stable fixed-point at ρ = 0, and θ is neutrally stable.
2.7. Results of simulation of the heterogeneous particles with complex variables
Figure 3(a) and figure 3(b) compare the results of the direct simulation with theresults(of simulation( M = 2, )Z = 1:
) of the (reduced)dynamics (18) for a system with
−1 −i
i 0
−i 1
, Ω1 =, γ1 = 1, µ1 = 0 and K = k
. The
Ω0 =
i
1
0 i
−1 −i
solid lines are the trajectories of the order parameter r of N = 10 000 particles obtained
from the direct simulation ofẋi = −xi r † K † xi + Ωi xi + Kr,
where the order parameter is calculated by
1 ∑
r=
xi .
N 1≤i≤N
The numerical solution of xi is obtained by using the reference frame rotating with
Ωi . The dashed lines show the trajectories of r obtained from the reduced dynamics(21). The results from the reduced dynamics agree quite well with those of the direct
simulation. By varying the value of k, I obtain the fixed point and the limit-cycle
oscillation of r (figure 3(a) and (b)). This simple, low-dimensional behaviour of thelimit cycle is in sharp contrast with the complicated motions of individual particles.
Figure 3(c) shows the trajectories of two particles in the system of figure 3(b). One ofthe particles is entrained into the collective synchronisation whereas the other particle
has a complicated trajectory. Figure 3(d) shows how the radius of the limit cycle depends
on k. The radius obtained by solving (22) (dashed line) agrees well with that obtainedby the direct simulation (crosses) and by the reduced dynamics (circles).
3. Conclusions
In this study, we have extended the Kuramoto–Sakaguchi model to model particles
interacting on a high-dimensional unit sphere. The dynamics are described by a type ofmatrix Riccati diﬀerential equation, which is characterised by an external force vector
and a natural-frequency matrix. The position of a particle at a given time is obtained by
the projection transformation of the initial position. This result is an extension of the
Watanabe–Strogatz theory. By assuming that the particles are uniformly distributedon a unit sphere in the initial conditions, the centre of mass of the particles with thesame natural-frequency matrix is determined by the vector of the parameters of the
projection transformation. In a system of particles with complex variables and naturalfrequency matrices with a multivariate Lorentz distribution, the motion of the centreCollective motion of particles on a sphere(b)(a) Re(x2)(c)

Direct simulation
Ott–Antonsen

Direct simulation
Ott–Antonsen–0.3 –0.2 –0.1 0 0.1
Re(x1)

Re(x2)
1
0.5
0
–0.5
-0.5

0
Im(x1)

0.5

Re(x2)
0.6
0.4
0.2
0
–0.2
–0.4
–0.2 –0.1

0.2 0.3(d)0.9Re(x1)
–0.5
0
0.5
1 1

|r|

0.3
0.2
0.1
0
–0.1
–0.2
–0.3
–0.2
Im(x1) 0
0.2

20

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

–0.2
–0.1
0 Re(x )1
0.1
0.2
0 0.1 0.2 0.3 0.3
Im(x1)

Direct simulation
Ott–Antonsen
Theory
0

0.5

1

1.5
k

2

2.5

3

Figure 3. Dynamics of the system of N = 10 000 heterogeneous particles. (a,b) Direct
simulations of the centre of mass (solid line) and low-dimensional dynamics (dashedline) for (a) k = 0.5 and (b) k = 1. (c) Trajectory of two particles for k = 1. (d)Direct simulation of |r| at t = 100 (crosses), simulation of (18) at t = 100 (circles) andtheoretically derived radius of the limit cycle, R1 (dashed line).

of mass of all the particles can be described by a high-dimensional extension of theStuart–Landau equation. This result is an extension of the Ott–Antonsen theory. A
periodic solution of the extended Stuart–Landau equation agrees with the motion of
the centre of mass of the system. Thus, we have shown that the collective motion ofa system composed of elements with large degrees of freedom can be reduced to the
dynamics of a low-dimensional system.
Although the present model is an extension to high degrees of freedom of theKuramoto–Sakaguchi model, it diﬀers from the previously proposed extensions [31, 32].
For example, Ritort introduced a variable into the model to keep the oscillators on theunit sphere [31], whereas the particles of the present model (described by (9)) remainon the unit sphere without additional terms. In addition, unlike the dynamic variablesof the model by Gu and coworkers, which are represented by matrices [32], the dynamic
variables in the present model are represented by vectors. Finally, unlike the present
model with heterogeneous particles, neither of these two models has been reported tobe reducible to low-dimensional systems.
The present model provides a method to study new types of collective phenomena
by reducing the system behaviour of elements with large degrees of freedom to that ofa low-dimensional system. The present model and its low-dimensional description canbe applied to problems already studied from the viewpoint of Kuramoto oscillators [15],such as time-delay systems [33, 34], systems of multiple-peak natural frequencies [22],Collective motion of particles on a sphere21

non-local coupling [35], dynamics on complex networks [36], associative memory [37]and common-input synchronisation [38]. In particular, because the system of oscillatorswith heterogeneous interaction delays [33] and the system of oscillators whose natural
frequencies obey a mixture of Lorentz distributions [22] have been studied by using the
Ott–Antonsen theory, I expect that the high-dimensional extension of these systems canbe solved by the present method.
I have not found a physical system whose dynamics are described by the present
model. However, the phenomena observed in the present model might be usefulin interpreting experimental observations of the systems that can be regarded as apopulation of particles on a sphere. In particular, the normalised velocity of birds
in a flock approximated by the Heisenberg model [3] could be analysed by using thereal-valued systems with M = 3.
The results of the model suggest that the statetransition from the desynchronised state to synchronised state can occur in systems
with heterogeneous particles (figure 3(d)). The present results also suggest that the
trajectory of the centre of mass can be a very simple limit cycle (figure 3(b)) even if
individual particles in the system exhibit complicated trajectories (figure 3(c)). The
present model could be used to approximate and to analyse systems exhibiting theseproperties.
Acknowledgments
This work was supported by MEXT/JSPS KAKENHI Grant Numbers 24651184 and25115710.
References[1] Hoppensteadt F C and Izhikevich E M 1997 Weakly connected neural networks (Springer-Verlag)[2] Hayakawa Y 2010 Europhys. Lett. 89 48004[3] Bialek W, Cavagna A, Giardina I, Mora T, Silvestri E, Viale M and Walczak A M 2012 Proc. Natl.
Acad. Sci. 109 4786–4791[4] Zhabotinsky A M and Zaikin A N 1973 J. Theor. Biol. 40 45–61[5] Kiss I Z, Rusin C G, Kori H and Hudson J L 2007 Science 316 1886–1889[6] Appleton E V 1922 Automatic synchronization of triode oscillators Proc. Cambridge Phil. Soc.
vol 21 pp 231–248[7] Nana B and Woafo P 2006 Physical Review E 74 46213[8] Ermentrout B 1994 Neural Comput. 6 679–695[9] Strogatz S H et al. 2005 Nature 438 43–44[10] Buck J 1988 Q. Rev. Biol. 63 265–289[11] McClintock M K 1971 Nature 229 244–245[12] Glova A F, Kurchatov S Y, Likhanskii V V, Lysikov A Y and Napartovich A P 1996 Quantum
Electronics 26 500–502[13] Tsang K Y, Strogatz S H and Wiesenfeld K 1991 Phys. Rev. Lett. 66 1094–1097[14] Kuramoto Y 1984 Chemical Oscillations, Waves, and Turbulence (Springer-Verlag)[15] Acebrón J A et al. 2005 Rev. Mod. Phys. 77 137–185[16] Watanabe S and Strogatz S H 1994 Physica D 74 197–253[17] Goebel C J 1995 Physica D 80 18–20

Collective motion of particles on a sphere22[18] Marvel S A, Mirollo R E and Strogatz S H 2009 Chaos 19 043104[19] Laing C R 2009 Physica D 238 1569–1588[20] Atsumi Y and Nakao H 2012 Phys. Rev. E 85(5) 056207[21] Ott E and Antonsen T M 2008 Chaos 18 037113[22] Martens E A, Barreto E, Strogatz S H, Ott E, So P and Antonsen T M 2009 Phys. Rev. E 79
026204[23] Pikovsky A and Rosenblum M 2011 Physica D 240 872–881[24] de Gennes P G 1972 Phys. Lett. A 38 339–340[25] Bowers R G and McKerrell A 1973 J. Phys. C 6 2721
[26] Albert R, Jeong H and Barabási A L 2000 Nature 406 378–382[27] Pikovsky A, Rosenblum M and Kurths J 2001 Synchronization: a universal concept in nonlinear
sciences (Cambridge University Press)[28] Winfree A T 2001 The geometry of biological time (Springer-Verlag)[29] Reid W T and Binkley C 1972 Riccati Diﬀerential Equations (Academic Press New York)[30] Verhulst F 1996 Nonlinear diﬀerential equations and dynamical systems (Springer-Verlag)[31] Ritort F 1998 Phys. Rev. Lett. 80 6–9[32] Gu Z M, Zhao M, Zhou T, Zhu C P and Wang B H 2007 Phys. Lett. A 362 115–119[33] Lee W S, Ott E and Antonsen T M 2009 Phys. Rev. Lett. 103 044101[34] Abrams D M and Strogatz S H 2004 Phys. Rev. Lett. 93 174102[35] Kuramoto Y and Battogtokh D 2002 Nonlinear Phenom. Complex Syst. 5 380[36] Boccaletti S S, Latora V, Moreno Y, Chavez M and Hwang D U 2006 Phys. Rep. 424 175–308[37] Aoyagi T 1995 Phys. Rev. Lett. 74 4075–4078[38] Teramae J N and Tanaka D 2004 Phys. Rev. Lett. 93 204103