A linear reformulation of the Kuramoto model of self-synchronizing oscillators

0.11 MB·English

by  David C. Roberts

A linear reformulation of the Kuramoto model of self-synchronizing coupled
oscillators
David C. Roberts
Theoretical Division and Center for Nonlinear Studies,
Los Alamos National Laboratory, Los Alamos, NM
8 (Dated: February1,2008)
0
0 Abstract
2

n The present paper introduces a linear reformulation of the Kuramoto model describing a self-
a
J
synchronizing phase transition in a system of globally coupled oscillators that in general have different

3
2
characteristic frequencies. The reformulated model provides an alternative coherent framework through

]
S whichonecananalytically tacklesynchronization problemsthatarenotamenabletotheoriginalKuramoto
P
. analysis. It allows one to solve explicitly for the synchronization order parameter and the critical point
n
i
l of 1) the full phase-locking transition for a system with a finite number of oscillators (unlike the original
n
[
Kuramotomodel,whichissolvableimplicitlyonlyinthemean-fieldlimit)and2)anewclassofcontinuum

5
v systems. Italsomakesitpossibletoprobethesystem’sdynamicsasitmovestowardsasteadystate. While
6
6 discussion inthis paper isrestricted tosystems withglobal coupling, thenew formalism introduced bythe
1
1
linearreformulation alsolendsitselftosolving systemsthatexhibitlocalorasymmetriccoupling.
.
4
0
7
0
:
v
i
X
r
a
1
I. INTRODUCTION
Spontaneous synchronization of coupled oscillators with different natural frequencies is at the
core of many striking phenomena in dynamical systems in realms from biology and physics to
social dynamics [1, 2]. The complexity of these systems defied attempts to encapsulate them in
tractable mathematical formulations until, in 1975, Kuramoto produced a model that he was able
to solve exactly for systems containing an infinite number of globally weakly coupled nonlinear
oscillators [3, 4, 5, 6]. In such a system, the global coupling strength across the oscillators and
the width of the oscillators’ initial characteristic frequency distribution determine whether or not
thesystemwillself-synchronize,andwiththeseKuramotowas abletodescribeaphasetransition
intoaself-synchronizingsystem. TheKuramotomodelhassinceprovidedthebasisformanylater
effortstoexplorespontaneoussynchronization. Moreimportantlystill,itisavaluableexplanatory
model in its own right for understanding numerous situations involvingsynchronization. In addi-
tiontoexamplesinbiology(see[1,6]and reference therein)suchasneuralfiring patterns,clouds
offireflies flashing as one, and thecoordinatedaction ofcardiac pacemaker cells, manyexamples
of theKuramotomodel havebeen recently discoveredthroughoutthephysicalsciences including
in such diverse systems as Josephson junction arrays [7], collective atomic recoil lasing [8], fla-
vor evolution of oscillating neutrinos [9], and phase locking of oscillations in coupled chemical
reactions [10].
This paper reformulates the Kuramoto model in terms of linear dynamics, permitting its so-
lution through an eigenvalue/eigenvector approach. The analysis here is restricted to solving for
the critical point of the fully locking transition and the synchronization order parameter beyond
this transition;the properties ofpartially locked states are not considered. Withinthis regime, the
reformulation of the Kuramoto model has a number of advantages over the original model. In
particular, in the continuumlimit,this linear model makes it possibleto solveexactly a new class
ofsynchronizationmodelsthatisdistinctfromtheclassofsystemsforwhichtheKuramotomodel
in its original form has an analytic solution. Furthermore, whereas in the original version of the
Kuramotomodel the synchronizationorder parameter is only solvablein thecontinuumlimit and
then only implicitly,thealternativeform presented below allowstheorder parameter to be solved
explicitly and for any number of oscillators. In addition, the linearity of the reformulation makes
it possible to investigate the time evolution of a system’s self-synchronization and lends itself to
adaptationtosystemsthatexhibitlocal and/orasymmetriccouplingbetween oscillators.
2
The present paper will begin with an introduction to the original Kuramoto model and its im-
plicit analyticsolutionin the thermodynamiclimit. Then the linear reformulationand its solution
will be presented, and the mapping between the Kuramoto model in its currently accepted form
and linear version presented here will be shown. This will be followed by some examples to
demonstratepropertiesofthelinearreformulation.
II. KURAMOTOMODEL
TheKuramotomodel[3,4,5,6]describesacollectionofN oscillatorsthatareweaklycoupled:
N
θ˙ = ω + K sin(θ θ ) (1)
k k jk j k
−
j=k
X6
where K is the coupling constant, ω is the characteristic frequency of the kth oscillator, and θ
jk k k
isitsphase. AlthoughKuramoto’smethodonlyyieldsan analyticsolutionforauniformcoupling
scheme (further discussion below), note that the model’s coupling constant is general and does
not imply uniform coupling. Kuramoto was able to exactly solve eq. (1) for a system of globally
and uniformly coupled oscillators,i.e. K = K/N, with theconstraintsthat the systembe in the
jk
thermodynamic limit, i.e. N , and have reached a steady-state; in so doing, he showed the
→ ∞
existenceof anontrivialself-synchronizationphasetransitioninsuch a system. In thispaper, as a
measure of the synchronization of a system at a point in time, we will use the amplitude r of the
complexphaseorderparameter,
N
1
reiφ eiθj, (2)
≡ N
j=1
X
(where φ represents the collective phase of the synchronized state) which ranges between 0 (no
synchronization) and 1 (perfect synchronization); steady-state synchronization of the system will
be specified as such. Kuramoto’s result was an implicit equation for r in the synchronized state
givenby
π
2
r = Kr cos2θg(Krsinθ)dθ (3)
π
Z−2
whereg = g(ω )isthedistributionofω andtheframeofreferenceistherotatingframeinwhich
k k
the frequency of the synchronized solution is zero. We will refer hereafter to eq. (3), with the
abovementioned accompanying constraints, as the Kuramoto solution. The Kuramoto model eq.
(1), alongwithitsmanyvariations,has beenstudiedingreat detail[6], andisseenas thestandard
modelforsynchronizationofgloballyweaklycoupled nonlinearoscillators.
3
III. LINEARREFORMULATIONANDSOLUTION
WenowproposealinearreformulationoftheKuramotomodelofspontaneoussynchronization:
N
ψ˙ = (iω γ)ψ + Ω ψ , (4)
k k k jk j
−
j=k
X6
where Ω is the coupling constant of this linear model and γ is the decay constant to be tuned to
jk
bring theamplitudeof ψ — a complexvariable— to a steady state(details givenbelow). As we
k
willlatershow,theargumentofψ correspondstoθ ineq. (1). Forsimplicity,weassumeγ > 0,
k k
and ω isreal. In thisreformulationweconsiderglobalcouplingbecauseitisconvenientandwill
k
permiteasycomparisonofourresultswiththoseofpreviouslystudiedmean-fieldmodels,butthis
approach lendsitselfequally welltootherlinearcouplingschemes [11].
Thelinearmodel(eq. 4)has thesimplesolution
N
ψ~ = a ~v eλjt, (5)
j j
j=1
X
wherea areconstantsdeterminedbytheinitialconditions,and~v andλ aretheeigenvectorsand
j j j
eigenvaluesassociatedwiththematrixdefinedbytheRHSofeq. (4). Thesynchronizingbehavior
of the system in the long-time limit is dependent on the eigenvalue(s) with the greatest real part.
Thus in this analysis we will adopt the convention of ordering all eigenvalues by their real part,
from λ (least) to λ (greatest). Distinct eigenvalues with the same real part shall arbitrarily be
1 N
assigned consecutivesubscripts within the larger sequence. In our discussion we will assume λ
N
isnot degenerate.
Wetuneγ so that [λ ] = 0. If [λ ] = 0 thenthesolutionbecomes
N N 1
ℜ ℜ − 6
lim ψ~ = a ~v eiωrt (6)
N N
t
→∞
wherethecollectivefrequency ofthefully lockedstate,ω , isgivenby
r
ω = iλ = λ (7)
r N N
− | |
and is related to the collective phase of the locked state by φ = ω t. As a result, r will tend to
r
a steady-state value between 0 and 1, which (in a finite-N system) indicates full locking, where
the entire oscillator population is locked to one particular frequency. The ensuing analysis will
focus on thetransitionto thisfully lockedstate(thefull lockingtransition). Forfinite-N systems,
4
the transition may be from partial locking — where there are subpopulations locked to different
frequencies — tofulllocking,orfromincoherence tofulllocking.
One can calculate the properties of the steady-state synchronization phase, where only one
eigenvector remains in the long-time limit. To find an expression for r one must determine the
eigenvector associated with λ . Hereafter, for simplicity we assume Ω = Ω > 0 unless other-
N jk N
wise specified. The general form of thecorresponding eigenvector for 1 j N, where j is the
≤ ≤
indexforthecomponentsof~v , isgivenby
N
i(ω ω ) Ω/N γ
(v ) = N − r − − . (8)
N j
i(ω ω ) Ω/N γ
j r
− − −
The general explicit expressionfor r of a fully phase-locked system with an arbitrary distribution
ofω overfiniteN in thelong-timelimitisthen
k
N N N
1 1 ψ 1 (v )
r eiθj = j = N j
≡ N N ψ N (v )
(cid:12) (cid:12) (cid:12) j (cid:12) (cid:12) N j (cid:12)
(cid:12) Xj=1 (cid:12) (cid:12)Xj=1 | |(cid:12) (cid:12)Xj=1 | |(cid:12)
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
(cid:12)1 N i(ω(cid:12) ω (cid:12)) Ω/N(cid:12) γ (cid:12)(ω ω )2 +(cid:12) (Ω/N +γ)2
= (cid:12) (cid:12)N − r(cid:12) − (cid:12)− (cid:12) j − r (cid:12) , (9)
N (cid:12) i(ωj ωr) Ω/N γ s(ωN ωr)2 +(Ω/N +γ)2(cid:12)
(cid:12)Xj=1 − − − − (cid:12)
(cid:12) (cid:12)
which, it should be str(cid:12)essed, is independent of initial conditions. If we assume th(cid:12)e distributionof
(cid:12) (cid:12)
frequenciestobesymmetricaboutω ,i.e. g(ω ω ) = g(ω +ω ),thenwecansimplifyeq. (9)
r r k r k
−
toarriveat
1 N ω ω 2 −1/2
r = 1+ j − r . (10)
N Ω/N +γ
" #
j=1 (cid:18) (cid:19)
X
If r goes to a nonzero steady-state value, there is steady-state synchronization of the system. In a
finitesystem,thefulllockingtransitiontakesplacewherer goesbetweenhavingandnothavinga
steady-state value. This transition point occurs for a given frequency distribution g(ω ) when the
k
followingrelationshipwiththecriticalvalueofthecouplingconstantΩ issatisfied:
c
[λ (Ω )] = 0. (11)
N 1 c
ℜ −
To see the connection between the linear reformulation, eq. (4), and the Kuramoto model, eq.
(1), weperform thenonlineartransformationψ (t) = R (t)eiθk(t) on eq. (4)toarriveat
k k
N
Ω R
θ˙ = ω + j sin(θ θ ), (12)
k k j k
N R −
k
j=k
X6
N
Ω
R˙ (t) = γR + R cos(θ θ ). (13)
k k j j k
− N −
j=k
X6
5
By tuningγ so that (λ ) = 0 wecan force each R togo toasteady state,namely
N k
ℜ
(ω ω )2 +(Ω/N +γ)2
lim R (t) = a (v ) = a N − r . (14)
t→∞ k | N N k| | N|s(ωk −ωr)2 +(Ω/N +γ)2
If all R go to a steady state for large times, then eq. (12) becomes eq. (1). Therefore, this linear
k
versionmapsontoeq. (1)withtheeffectivecouplingconstant
Ω R Ω (ω ω )2 +(Ω/N +γ)2
K˜ = lim j = k − r . (15)
jk t→∞ N Rk Ns(ωj −ωr)2 +(Ω/N +γ)2
It is important to note that K˜ is independent of initialconditionsas a cancels out, and that the
jk N
mappingholdsonlyintheregimeinwhichthereissteady-statesynchronization(where [λ ] =
N 1
ℜ − 6
0). Withthis,eq. (12)can berewrittenas
N
θ˙ = ω + K˜ sin(θ θ ). (16)
k k jk j k
−
j=k
X6
In other words, by introducing the amplitude R properly constrained by the decay constant γ,
k
we would be able to perform a nonlinear transformation of the Kuramoto model eq. (1) (which
only has an implicit solution in the infinite-N limit) with an effective coupling constant of K˜
jk
into our linear version eq. (4) that can be solved exactly for any N. One can conceive of much
more general mappings between eq. (1) and eq. (4), such as by allowing each oscillator to have
its own independent γ γ , or by breaking the assumption of uniform coupling in the linear
k
→
reformulation, e.g. Ω Ω (which does not break the linearity) to accommodate a larger class
jk
→
ofcouplingsK .
jk
IV. EXAMPLES
To emphasize the mathematical properties of our linear model, we work through three groups
of examples in detail below. In each, we solve for r using the linear approach and compare the
result to theKuramoto solution,eq. (3), in thethermodynamiclimit(the regimeof validityof the
Kuramotosolution).
A. Identicaloscillators
In our first example we consider a system where all characteristic frequencies are equal, i.e.
ω = ω. Thisisasimplesystemwhichwillexhibitsteady-stateperfectsynchronizationregardless
k
6
of initial conditions (as long as Ω = 0). In this case, to solve the linear model (eq. (4)), we first
6
observethatthereareN 1degenerateeigenvalues,each equaltoiω γ Ω/N andoneunique
− − −
eigenvalue,
λ = iω γ +Ω(N 1)/N, (17)
N
− −
that has an associated eigenvector v = 1,1,1,...,1 . (Note that the normalization of the
N
{ }
eigenvector does not matter for the calculation of the order parameter because r is defined such
that every element has a modulus of 1). In order for each ψ to go to a steady state, we set
k
| |
γ = Ω(N 1)/N (consistent with eq. (13) for a perfectly synchronized state) so that as t
− → ∞
all of the degenerate eigenvectors decay away on a time scale N , which for large N reduces
(N 2)Ω
−
to 1/Ω. AssumingΩ = 0, thesolutionto thereformulationwithω = ω isthen
k
∼ 6
lim ψ~ = a ~v eiωt. (18)
N N
t
→∞
Thesynchronizationorderparametercannowbecomputedeasily,givingthesteady-statevalue
ofr = 1forlongtimes. Becauseeachelementinthesumisnormalizedbyitsmagnitude,thelong-
time steady-state behavior of r is independent of the initial conditions represented by a . Since
j
K˜ Ω/N in the long-time limit in this example, in the infinite-N limit the linear model’s
jk
→
solution,eq. (10), and theKuramoto solution, eq. (3), givethe identical result of r = 1 assuming
the frequency distribution ω = ω. It is also instructive to take small perturbations around the
k
characteristicfrequencies,i.e. ω = ω+ǫη where0 < ǫ 1andη aresymmetricallydistributed
k k k
≪
around 0. In thiscase, to leadingorder
∆2
r 1 ǫ2 , (19)
≈ − 2Ω2
where∆2 isthevariancegivenby ∆2 = 1 N η2.
N j=1 j
P
B. Bimodaldistributionofcharacteristic frequencies
Foroursecondexample,welookatasysteminwhichbothsteady-statepartialsynchronization
and a full-locking transition occur. Consider a bimodal distribution where ω = ω ∆ for
j 0 − 2
1 j N/2 and ω = ω + ∆ for N/2 < j N where ∆ 0 is the width of the distribution
≤ ≤ j 0 2 ≤ ≥
(we will assume N is even for simplicity). The solution to this problem maps onto a system of
twocoupled oscillatorsbecauseofthepecularityofthedistribution.
7
Solvingthelinearreformulation,weseethatthereare(N 2)/2degenerateeigenvaluesgiven
−
by i(ω ∆) Ω γ and (N 2)/2 degenerateeigenvaluesgivenby i(ω + ∆) Ω γ. The
0 − 2 − N − − 0 2 − N −
final twoeigenvaluesλ and λ are givenby
N N 1
−
Ω(N 2) √Ω2 ∆2
λ = iω γ + − − (20)
N,N 1 0
− − 2N ± 2
wherethe‘+’refers toλ , and the‘ ’refers toλ .
N N 1
− −
First let us take the situation where Ω < ∆. Here we shall see that r does not go to a steady-
state value and therefore that the system will not fully lock and no steady-state synchronization
occurs. Unlikethepreviousexample,therearenotonebuttwodistincteigenvalueswiththelargest
real part, namely ℜ[λN−1] = ℜ[λN] = −γ + Ω(N2N−2). So even after setting γ = Ω(N2N−2) we find
that
√∆2−Ω2 √∆2−Ω2
lim ψ~ = a ~v ei„ω0− 2 «t +a ~v ei„ω0+ 2 «t. (21)
(N 1) (N 1) N N
t − −
→∞
It is clear that as long as Ω < ∆, r will never go to a steady-state value and hence the whole
system will never fully lock or reach steady-state synchronization. However, independent of all
initial conditions and assuming Ω > 0, the population of oscillators will split into two perfectly
synchronizedsubsetswithdifferentcollectivefrequencies givenby [λ ] = ω √∆2 Ω2.
ℑ N,N−1 0 ± 2−
Now let us consider the situation where Ω > ∆. If we set γ = Ω(N−2) + √Ω2−∆2 (note
2N 2
γ Ω(N 1)/N for large Ω, consistent with eq. (13) for a perfectly synchronized state) then,
→ −
on atimescalegivenby 1/√Ω2 ∆2,all butoneeigenvaluedieoutovertime:
−
lim ψ~ = a ~v eiω0t (22)
N N
t
→∞
wherethecomponentsoftheeigenvectorare givenby
i∆ 2Ω/N √Ω2 ∆2
(v ) = − − − (23)
N j
i∆ 2Ω/N √Ω2 ∆2
− − − −
for1 j N/2 and (v ) = 1 forN/2 < j N. Thusr goesto asteady-statevalue:
N j
≤ ≤ ≤
1+ 1 ∆ 2
r = v − Ω (24)
u q 2
(cid:0) (cid:1)
u
t
for Ω > ∆. It is important to note that although this equation is true for finite N — unlike the
Kuramotosolutioneq. (3)— theexpressionforr isindependentofN.
Since when Ω < ∆ thesystem does not fully lock and no steady-statesynchronization occurs,
and when Ω > ∆thereisfull lockingandsteady-statesynchronization,itis clearthat when Ω =
c
8
∆ a “first-order” full-locking transition occurs, marked by steady-state partial synchronization
beyondr = 1/√2. Usingeq. (24)wecandeduceasquare-rootscalinglawoftheorderparameter
c
near the full-locking transition (which in this case is also the transition into a steady-state self-
synchronizingsystem),i.e.
1 Ω Ω
r r − c. (25)
c
− ≈ 2 Ω
r
As in the previous example, this bimodal example has the property of K˜ Ω/N in the
jk
→
long time limit, so the linear form maps onto the Kuramoto model (eq. (1)) with the simple
transformationΩ K. Fromtheabovediscussiontherefore,ifoneassumesabimodalfrequency
→
distribution, i.e. g(ω) = δ(ω−∆/2) + δ(ω+∆/2), then solving eq. (3), which is only valid in the
2 2
infinite-N limit,givesourresult,eq. (24).
C. Continuumlimit
In this final example, we work exclusively in the continuum limit and demonstrate that, even
within this regime where the Kuramoto solution is valid, the linear reformulation still has the ad-
vantageofasimplersolutionfortheorderparameterandopensupanewclassofexactlysolvable
synchronizationmodels. Wefirstdeterminethecriticalpointofthelockingtransitionandthesyn-
chronization order parameter given by the linear model for a system of oscillators with a general
symmetric characteristic frequency distribution and a particular effective coupling. Systems with
such coupling, which we will consider below, go from completely incoherent to fully locked at
one point in parameter space, i.e. there are no partially locked states; the partial-locking and full-
lockingtransitionsmerge [15]. Next,wesolvefor thecritical pointand thesynchronizationorder
parameter for a Lorentzian frequency distribution and a uniform frequency distribution. We then
lookat howtheresult differsfrom thetraditionaluniformcouplingsolution.
Webeginbysolvingforr,γ, and thelockingtransitionpoint(also thesynchronizationcritical
point), each in terms of Ω and a general symmetric frequency distribution. As N , eq. (10)
→ ∞
becomes
1/2
ω ω 2 −
r = ∞ dωg(ω) 1+ − r , (26)
γ
" #
Z (cid:18) (cid:19)
−∞
which holds for Ω > Ω . The anomalous scaling properties of this solution are discussed in
c
[17]. The linear reformulation eq. (4) in the continuum limit has the same form as the equation
describingtheevolutionofthefundamentalmodeinthestabilityanalysisin[14]oftheincoherent
9
stateintheKuramotomodel:
ψ˙(ω,t) = (iω γ)ψ(ω,t)+Ω ∞ g(ω )ψ(ω ,t)dω , (27)
′ ′ ′
−
Z
−∞
whichhas an effectivecouplinginthecontinuumKuramotomodelof
(ω ω )2 +γ2
K(ω,ω ) = Ω ′ − r . (28)
′
s(ω ωr)2 +γ2
−
Followingtheanalysisin[14],thespectrumofthelinearoperatoroftheRHSofeq. (27)comprises
an eigenvalue λ at the origin of the complex plane and a continuous line of eigenvalues along
N
γ oftherealaxis. IfΩ > Ω , γ > 0and λ liesapart fromtheothereigenvalues,so thereisfull
c N
−
locking and synchronization in the long-time limit. However if Ω Ω , then γ = 0, λ merges
c N
≤
withthecontinuumalongtheimaginaryaxis,andcontributionfromeveryfrequencyremainsinthe
long-timelimit,whichresultsinr = 0asalloscillatorsare“drifting”. Thisimpliesasecond-order
phasetransition. Inthiscase,thepartial-lockingtransition(oftenreferredtoasthesynchronization
phasetransition)andthefull-lockingtransitionoccuratthesamepoint. Furthermore,thecoupling
described by eq. (28) is such that, at the onset of the transition, there is a uniform distribution of
phases from 0 to 2π across theoscillatorpopulation(i.e. r 0 as Ω Ω ). We can determineγ
c
→ →
usingtheself-consistencyequation[14]
g(ω)
∞
1 = Ω dω. (29)
γ iω
Z−∞ −
At most one solution for γ exists and γ 0 [12]. Therefore the system will reach steady-state
≥
synchronization for γ > 0, and the critical point Ω occurs as γ 0+. This can easily be shown
c
→
tobe
1
Ω = . (30)
c
πg(0)
Now, considering a Lorentzian distribution of frequencies about ω , i.e. g(ω ω ) =
r r
−
∆ , from eq. (30) we obtain Ω = ∆. Integrating eq. (29) gives γ = Ω ∆. So
π[∆2+(ω ωr)2] c lor −
−
whereΩ > Ω ,
c
1/2
2 Ω Ω 2 −
r = cos 1 c 1 c (31)
lor −
π Ω Ω − Ω Ω
(cid:18) − c(cid:19)" (cid:18) − c(cid:19) #
inthelong-timelimit. Similarlyforauniformdistributionoffrequenciesaboutω ,i.e. g(ω ω ) =
r r
−
1 for ω ω < π∆/2 and 0 otherwise, the same steps give us Ω = ∆, γ = ∆π cot(π∆),
π∆ | − r| c unif 2 2Ω
and
πΩ πΩ
r = cot c sinh 1 tan c (32)
unif −
2 Ω 2 Ω
(cid:18) (cid:19) (cid:20) (cid:18) (cid:19)(cid:21)
10