PHYSlCAELSEVIER

Physica D 100(1997)279-300

Self-synchronization of coupled oscillators with hysteretic responses
Hisa-Aki Tanaka a,b,., Allan J. Lichtenberg b, Shin'ichi Oishi aa Department of Information and Computer Sciences, Waseda Universi~, Tokyo 169, Japanb Electronics Research Laborato~ and Department of Electrical Engineering and Computer Sciences, University of California,
Berkeley, CA 94720, USA

Received 21 May 1996; revised 26 July 1996; accepted 26 July 1996
Communicated by Y. Kuramoto

Abstract
We analyze a large system of nonlinear phase oscillators with sinusoidal nonlinearity, uniformly distributed natural frequencies and global all-to-all coupling, which is an extension of Kuramoto's model to second-order systems. For small coupling,
the system evolves to an incoherent state with the phases of all the oscillators distributed uniformly. As the coupling isincreased, the system exhibits a discontinuous transition to the coherently synchronized state at a pinning threshold of thecoupling strength, or to a partially synchronized oscillation coherent state at a certain threshold below the pinning threshold.
If the coupling is decreased from a strong coupling with all the oscillators synchronized coherently, this coherence can persist
until the depinning threshold which is less than the pinning threshold, resulting in hysteretic synchrony depending on theinitial configuration of the oscillators. We obtain analytically both the pinning and depinning threshold and also expalin the
discontinuous transition at the thresholds for the underdamped case in the large system size limit. Numerical explorationshows the oscillatory partially coherent state bifurcates at the depinning threshold and also suggests that this state persists
independent of the system size. The system studied here provides a simple model for collective behaviour in damped driven
high-dimensional Hamiltonian systems which can explain the synchronous firing of certain fireflies or neural oscillators withfrequency adaptation and may also be applicable to interconnected power systems.
PACS: 02.50.-r; 05.40.+j; 05.90.+m; 84.30.Jc; 87.10.+e
Keywords: Phase model; Mutual entrainment; Hysteresis; Bifurcation; Adaption

I. Introduction
There has been a continuing effort to understand the collective synchronous behaviour in dynamical systems withmany degrees of freedom. Particularly, for the spontaneous emergence of synchronization, extensive studies havebeen devoted to large populations of limit-cycle oscillators ([ 1-3,5,6,15,16,20-24,26,29-31,38], see [27] for manyother interesting references). The studies have been motivated by intrinsic scientific interest [4,8,9,12,13,17,36,37]
along with the applications to engineering [7,11,13,32,33,35,39]. Recently, a breakthrough has been made in [35]* Corresponding author. E-mail: tanaka@oishi.info.waseda.ac.jp.
0167-2789/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved
Pll SO 167-2789(96)00 193-5

280

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

that relates the study of a practical current-biased series array of non-identical Josephson junctions and Kuramoto's
theory on spontaneous synchrony which was originally a product of an ingenious and bold insight in nonlinear
physics [15,16,24].
Certain interacting limit-cycle oscillators can be reduced to a phase model through averaging, allowing them tobe approximated by a set of first-order phase equations of the formdOi-~- O)i -}- 2..* nji(Oj -- Oi),
d--~-

(1)

J

in order to analyze their synchrony. On the other hand, systems with many degrees of freedom of the form
d20i-- Z Hji (Oj - Oi)
dt 2
)

(2)

can be Hamiltonian, having special dynamical characteristics [ 14,19,34]. A damped, driven generalization of (1)and (2)
d20i
dOi
m-~
+ dt : ~2i -1- Z(3)aji(Oj - Oi)
J

is a practical model for systems with many degrees of freedom in various fields of science. The question naturally
arises to the behavior of (3): do the oscillators synchronize as in (1) or do they exhibit dynamical patterns which
are a remnant of the dynamics in (2)?
Another motivation of this work is related to the progress in modeling the synchronous firing of the southeast Asia
fireflies - Pteroptyx malaccae. Unlike other species such as P.cribellata, P.malaccae has an ability to alter its firing
frequency as much as 15% to higher or lower frequency in response to the external signals, and this adaptation isconsidered to be responsible for the synchronous firing with very small phase lags [ 12]. Based on the physiologicalreasoning as well as the data from field work by Hanson [12], Ermentrout [8] developed a phase model explaining
the synchronous firing in swarms of P.malaccae. His idea is to introduce an adaptive effect of the firing frequency
to the phase model (1), as the simplest model for synchronous firing with small phase lags, in the following way:
Oi = O)i,N
O9i = E(O)i0 -- O)i) q- y ~ nji(Oj -- Oi).
j=l

(4)

Choosing small E and a particular function Hji (which can be nearly sinusoidal if o)i is in the range of adaptive
frequency), phase model (4) can simulate the synchrony observed in swarms of P.malaccae [8]. Interestingly, (4) isequivalent to (3) if (3) has an all-to-all, global coupling. More specifically, (4) can be transformed toN

mOi + Oi = ,.Qi @ ( K / N ) Zh(Oj - Oi),(5)

i = 1. . . . . N

j=l

by assuming that all Hji are identical and letting
told ~

E/new,~f2i =

EO)i,O,

Hji(Oj-- Oi) =(K/N)h(Oj

--

Oi) ,

m = E-2.

(6)Here, we study the particular model (5) with sinusoidal nonlinearity h(.) = sin(.). In addition to the interest in
firefly synchrony, this model may be significant for understanding the synchrony in power systems modeled by theswing equation [25], and also to extend the analysis of the Hamiltonian system studied by Konishi and Kaneko [ 14]to a damped, driven (continuous time) system.

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

281

1.O-

f

r

coherentlocking

0.5

0.0

incoherentstate

Kl°wer

Kupper

K

1.0-

f

r

0.5

1

1

0.0

Klc°wetK

KupperFig. 1. Hysteretic synchrony observed in (5) with system size N = 500. (a) As K exceeds the depinning threshold KD the system switches
discontinuously from the incoherent state to the coherent locked state. The data were obtained for m ~ 0.85 and uniformly distributed~Qi in [-5.0, 5.0]. (b) For larger m ( ~ 0.95) and the same natural frequency distribution if2M as in (a), the system exhibits an oscillationof the coherency r in the shaded region. The dotted line shows a time average of r.

In this paper, we consider the phase model (5) to understand the effect of inertia, namely frequency adaptation
to synchrony of oscillators. The main new finding is that the system can exhibit hysteretic synchrony dependingon the initial conditions. A schematic illustration of this phenomenon is shown in Figs. l(a) and (b). The verticalaxis represents a degree of the synchrony of the oscillators which will be explicitly defined in Section 2. Fig. ! (a)shows the simplest pattem obtained numerically for m = 0.75, .('2 M = 5.0 and N = 500 where the hysteresis ischaracterized by two different sets of (K, r) values at which a discontinuous jump across between the incoherent
state and coherent synchronized state. This kind of hysteretic synchrony has been observed in the charged density
wave (CDW) model studied in [29]. While in Fig. 1(b), an oscillation of the coherence r can be observed in theshaded region for m = 0.85 and the other conditions as in Fig. l(a), it suggests, that the synchrony in (5) has amore involved nature than in the CDW model [29]. The common feature observed in Figs. l(a) and (b) can be

H.-A. Tanaka et al. / Physica D 100 (1997) 279-300

282

summarized as follows: (i) there is an upper/lower limit of K where the coherent synchrony onsets/disappears; and(ii) the coherent synchronized state and the incoherent state can coexist for a certain coupling strength K between
the upper and lower limits of K.
This paper aims to clarify these characteristics in the limit of large N and to examine the oscillating sate of
the coherence. The paper is organized as follows. In Section 2, to make our starting point clear, we review the
self-consistent theory of Kuramoto's model in its simplest form. In Section 3, we develop a self-consistent theory
for the synchrony of (5) in the limit of N --+ ~ . The main idea is an extension of Kuramoto's theory which
involves the self-consistent mean field theory taking different initial configurations of oscillations into account. In
Section 3.1, we introduce an order parameter r to decouple (5) into a swarm of single pendulum oscillators andconsider the hysteretic response of each oscillator. In Section 3.2, we start from a system with a distribution of thenatural frequencies 12M, for which two different self-consistent equations are derived, corresponding to the twodifferent initial configurations of oscillators. We then focus on one simplifying case in Section 3.3, that permits ananalytical treatment of the self-consistent equation- that of uniformly distributed natural frequencies in a finite range[--S2M, S2M] and m • K >> 1. From the self-consistent equation, which is perturbatively approximated under these
conditions, analytic results are obtained in Section 3.4, i.e. the lower and upper critical points and the associateddiscontinuous jumping between incoherent and coherent synchronized states are obtained. Specifically, the lower
critical point is shown to be the same as that from the original Kuramoto model, while the upper critical point isobtained from a set of nonlinear finite-dimensional equations and shown to have a greater K than the lower criticalpoint.
In the numerical part (Section 4), we numerically simulate the model (5), and verify the existence of hysteretic
synchrony, where the lower and upper critical points are obtained and show a nice agreement with the theoreticalprediction in Section 3.4. In addition to the hysteretic synchrony phenomena, we examine the oscillatory state of thecoherence r for different system size N and coupling strength K to explore the size dependence of the oscillatory
region.

2. Kuramoto's theory in its simplest form
The "phase transition"-like synchronization onset is a large population of interacting oscillators and was firstdescribed by Winfree in [36]. An analytically tractable model, the Kuramoto model, was then proposed by Kuramoto
to elucidate a subtle connection between collective synchronization and phase transitions [ 15,16,24]. The Kuramotomodel, which is derived by averaging from a certain class of weakly, globally coupled limit-cycle oscillators, takes
the following form in the simplest case:N

Oi = ~)i + ( K / N )Z sin(0j - Oi) ,i = 1. . . . . N.

(7)

j=l

H e r e , 0 i and ,.(2i a r e respectively, the instantaneous phase and the natural frequency of the ith oscillator, and K/N > 0

is the coupling strength. In the large N limit, we assume that the intrinsic frequencies have a certain distributiongff2). Since (7) has the arbitrariness in shifting all Oiby S2t, we can assume the distribution g (I-2) becomes zero-mean
after a shift. Here, we consider a unimodal, symmetric distribution g(I-2) with g(12) = g ( - 1 2 ) which includes theGaussian and also the uniform distribution, although certain bimodal distributions can be considered as in [ 16]. The"all-to-all" global coupling in (7) might seem unrealistic at the first sight, but this global coupling arises naturally incertain circumstances; for example, in realistic, electrical systems such as a Josephson array [35] and the interactionof quasi-optical oscillators with a cavity [39]. The uniform global coupling in (7) provides a simplifying property

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

283

that enables one to utilize the "mean field" description of macroscopic behavior. This can be realized by use of"order parameters" r and 4~ defined byN

re i4' = ( l / N ) ~

e i0/

(8)

j=l

in which r = r (t) and ~ = 4' (t) can be interpreted as the degree of synchrony and the mean phase angle, respectively.
From a trignometric identity and (8), one can transform (7) equivalently to the following set of equations:

Oi = F 2 i + K r s i n ( q ~ - O i ) ,i=1. . . . . N.

(9)

In (9), each oscillator is equivalent to an overdamped pendulum with torque £2i and restoring force proportional toK r. Thus, the order parameters r and q~determine the time evolution of Oi, and Oi determines r and q~ self-consistentlyby (8). Using the above formalism, Kuramoto obtained the following key insight:[AI] There are self-consistent steady solutions of (8) and (9). In particular, there are solutions with r(t) and qS(t)being constant, in which all fluctuations vanish in the large N limit.
In the present situation, q~ can be set to 0 without loss of generality from the arbitrariness of Oi in (7). One "trivial"
solution is found to be r(t) -~ 0, where all oscillators rotate at their intrinsic frequencies. This completely incoherent
solution exists for all K, but is not necessarily stable. A secondary family of steady solutions bifurcates at K = Kc
with

Kc =

2

7rg(O).

( 1O)

This critical value can be obtained from the self-consistent formalism which is briefly explained as follows. For
a given g ( ~ ) and constrants K >_ Kc and r > 0, the oscillators are separated into two groups. One group is themutually synchronized oscillators with I~] < Kr (denoted by [S]), and the other group is the drifting oscillatorswith t~[ > Kr (denoted by [D]). In the synchronized group [S], each oscillator is locked to the mean phase withphase lags determined by (9). A drifting oscillator in [D] with a natural frequency £2 has a velocity 0 also determinedby (9), Thus, given the natural frequency distribution g(~2), the normalized density of the oscillators in [S] and [D]can be, respectively, obtained explicitly as ns (0) and nD (0, ~Q) through [15,16,24]. This static density of oscillatorson the unit circle of 0 ~ [0, 2Jr) determines the contribution to the mean phase from the oscillators in IS] and [D].
Namely, the contribution of the locked oscillators [S] is expressed as
7r/2
rlock

C~ns(O)e iO d O = Kr [*J

cos20g(Krsin0)d0,

(11)

*]

-7rt2
while the contribution from the drifting oscillators [D] is expressed asra,.ift=ffei°nD(O,,f2)g(S'2)d,.(2dO=O.

(12)

IS21>Kr

This particular form of rlock appears again for (5) as shown in Section 3.2. The sum of rtock and rdrift again determines
the mean phase such that the following self-consistent equation must be satisfied [ 15,16,24]:r ----rlock -4- rdrift.

(13)

We see in (12) that rdrift vanishes due to the symmetry of g(I2), g(£2) = g(-X2), i.e. in the Kuramoto model withsymmetric g(£2), a swarm of evenly distributed oscillators has no effect on r. However, this cancellation does not284

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

hold for (5), even for a symmetric g(O), as we see in Section 3. Thus, assuming K r sin 0 is small in g(.) of (l 1),
we are led to the following equation from (13):r = rlock q- rdrift = rlock --~ ½YrKrg(O) - lzrK3r3g(2)(O) h- O(r4).

(14)

Dividing (14) by r and letting r ~ 0, we find Kc = 2/yrg(0) and the coherence factor r grows as (K - Kc) 1/2
just above the transition at K = Kc. The exponent ½ is due to the pure sine nonlinearity in (7) [6], and thus shows
the analogy to the second-order phase transition in magnetic spins. It should be mentioned that much progress has
been made in this formalism; see for example [24] for modulated sinusoidal coupling function, and [6] for certain
general periodic functions.

3. Extension of Kuramoto's theory to a second-order system

In this section, we develop a self-consistent theory for the coherence r in (5) (hereafter, we call (5) the extended
Kuramoto model). The basic idea is the same as Kuramoto's, however, the analysis for the extended K ramotomodel becomes more involved because (a) Eq. (5) is genralized from a second-order system (a pendulum equation)
whose static states can be multi-stable and, therefore, dependent on the initial conditions, and (b) there is no longer
an exact solution to (5) for the drifting oscillators. We start from general m, K and g(O) with g(O) = g ( - O ) , andderive self-consistent for different configurations of initial conditions taking the above property (a) into account,then focus on one simplifying case of m, K and g(O) which enables us to give a good approximation of the solution
and the density of drifting oscillators on the unit circle.

3.1. Governing equation and hysteretic curve o f ( O , (0'})

Utilizing the identity ( K / N ) Y~7=l sin(0j - Oi) = K r sin(4~ - Oi) as obtained from (8) and a trigonometricidentity, (5) can be transformed to
mOi -~-6i --:--Oi - K r s i n O i ,i = 1. . . . . N,

(15)

which is the governing equation for a single damped driven pendulum. A damped driven pendulum is known to havehysteresis with respect to the initial conditions (Oi, Oi), which can be characterized by two different driving frequencies: the pinning frequency Op and the depinning frequency OD. As shown in Fig. 2, a pendulum starts whirling
once the applied torque O goes beyond a certain threshold S2D. This S'2 D is characterized by the disappearance ofthe equilibrium point determined by (15): 0 = sin - l ( O / K r ) = 4-½~r. Then I2D ( > 0) is given by S2D = K r . On
the other hand, once a pendulum starts whirling, it continues to rotate even though the torque O is varied down
below OD. This is due to the inertia of the pendulum and can be understood by the appearance of the second kindof periodic solution in a two-dimensional phase portrait of (15) which is well understood (for example, see [ 18]).
If O is again brought down, the frequency of ,O ~ /0) which is the average of 6 over one period, tends to 0
continuously as Op is approached. Specifically, if the damping is sufficiently small (underdamped case), (6) tends
to 0 as [ln(O - Op)] -r, which means the falling of the frequency of the running periodic solution 0, becomes
extremely steep as O tends to Op [28]. Furthermore, in this underdamped case an analytical approximation of Op
becomes possible, using Melnikov's method [10] (see Appendix B for the derivation of I2p).

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

285/z~ S"/
pinning frequency-\~pdepinning thresholdfrequency : ~D

Fig. 2. (/2, (0)) curve for a single oscillator.

, 1,f,,:..,.......................
m ~ l ~ ~ n ~....
m : locked oscillators" case [ l]• case [11],, : drifting oscillators

Fig. 3. Locked oscillators [S] and drifting oscillators [D] distributed on (~2, (#)) curve.

3.2. Self-consistent equation f o r two different cases
We now look for a solution of (15) satisfying [A 1] by deriving a self-consistent equation determining the coherence
r. There are two cases to be considered: case [I] in which the coupling strength K is increased from the incoherentr = 0 state; and the case [II] in which K is decreased from the coherent state with r ~ 1. In both cases, we start by
assuming the state where there are swarms of locked and drifting oscillators and they have some coherence qock andrdrift, respectively. In case [I], suppose a certain r is given at any moment, then r determines the pinning threshold~2p and the initially drifting oscillators with I2 < ~(2pcan be entertained to the locked oscillators [S] after a certaintransient. Thus, it is reasonable to assume that the swarms of [S] and the drifting oscillators [D] can be separated
at I2 -- ~2p in the natural frequency distribution as shown in Fig. 3. On the other hand, in case [II] initially locked

H.-A. Tanakaet al./Physica D 100 (1997) 279-300

286

oscillators are desynchronized and fall into the drifting state [D] once I2 exceeds I2D(= Kr). Hence, IS] and [D]can be assumed to be separated at the depinning frequency S-2 = 12D.
The coherence rlock for case [II] takes exactly the same form as (11) in Kuramoto's theory:,~/2
II = frlock

ns(O)eiOdO= fg(n)eiOdn = Kr f{.(~ {< .&"2D

COS2 0 g ( K r sin 0) dO,(16)--7r/2where ns is the normalized density of the synchronized oscillators with phase 0. In (16), the relations ns(0) =g(12) d I 2 / d O and Kr cos 0 = d ~ / d O were used. For case [1], oscillators with I2 > S2p are still drifting on theunit circle and the coherence of the locked oscillators takes the following form:
sin-1 [~2p/(Kr)]rock = f ns(O)ei°dO = f/

g(I2)ei°da'2=

IS21<S2v-

ei°g(Kr sin O)K cos 0 d0

sin-1 [~?p/(Kr)l0e= Kr f cosZOg(Kr sinO) dO,
-op

(17)

where 0p = sin -1 ~2p/(Kr) and 0p < lrr follows from 12e < .('2D = Kr.
Now we consider the effect from [D] on r:

I'I'ffei°nD(O,n)g(S2)dOdn,(18)

rdrift -~I.Ql>.f2p,owhere nD(0, $'2) is the density of the desynchronized oscillators [D]with phase 0 and given by driving frequency~2. As nD(O, I2) is proportional to 101 - l and f riD(O, I2) dO = f :nD(O, 12)0 d t = 1, we obtainriD(O, 12) = 7"-1101-1 = (~/2zr)101 -~.

(19)

I,I1
with T and ~ being the period and the frequency of the running periodic solution of 0, respectively. Thus, rdrift are,respectively,

1,,, =
rdriftff eiO(t,n)(~/2~r)lOl-,Og(12)dtdI2 = 2---~
1I£2I>£2p,D 0/

f eiO(t.n)l~lg(S2)dtd~"(20)

IF2I>f2p.D 0

Here 0 (t, I2) denotes the running solution of (15) with driving frequency I2, and the following properties were used:
0 > 0 for a'2 > 0, 0 < 0 for K2 < 0, and ~ ( - 1 2 ) = - ~ ( S 2 ) . In (20), cos 0 (t) and sin O(t) are periodic in t because
0 (t + T) =-- 0 (t) + 2re, Yt. Then, the integration of cos 0 (t) and sin 0 (t) over one period T does not depend on the spanover which the integration is made. Hence, we may assume that 0(0, a'2) = 0 without loss of generality. Also from(15) and the above initial condition 0(0, I2) = 0, we have m(-O(t, J2)) + (-O(t, £2)) = -a'2 - Kr s i n ( - 0 ( t , a'2))
such that

O(t, --a'-2) = --O(t, 0).
From g(I2) = g ( - ~ 2 ) , (20) takes the following form:(21)14.-A. Tanaka et al./ Physica D 100 (1997) 279-300

287

i"

I.,

,ff

rdrift = 2-'-~[COS0(t, ~ ) + COS0(t, --I-2)]l~lg(a'2) dt d$2$'2 > ff2p, D 0[sin O(t, $-2) + sin O(t, - I 2 ) ] J~lg($-2) dt dI2

+ ~

(22)

a"2 > ~f2p, D 0

and can be simplified by (21) tordrif t =

,,,

' f f

-~

,..,

.Q > ff2p,D 0

Thus, the two self-consistent equations for r I'II,
rI.ll = rl,lll.IIlock q- rdrift'(24)

are obtained from (16), (17) and (23).

3.3. Perturbation approximation of self-consistent equations
In this section we consider a condition that enables us to treat (24) analytically. For this purpose, the integrationin (23) is performed perturbatively. One way to do this without making the inertia term m~J small is to startfrom the underdamped case of (15) with uniform distribution of natural frequency $2. We consider this in detailbelow.
First, in a new time-scale: told = mtnew, (15) takes the form

Oi +Oi =m,Qi -mKrsinOi,=m' (-~,f2i) - m ' ( -m~ K )rsinOi.

(25)

This implies that the original (15) with given m is equivalent to a new system with arbitrary m' and an associatedrescaled distribution of I2i, together with the rescaled K. For instance, suppose that we have a coherence r for(15), with m = 1.0, uniformly distributed $2 E [---QM, $2M] with ,f2M = 5.0 and K = 10.0, then the same rwill be obtained for m = 0.1, $2M = 50.0 and K = !00.0. Thus we do not have to consider m and I2M(/K)
separately and we are naturally led to one possible way to introduce a small parameter 8 in the underdamped ( 15):8 ~ ( K m ) - I << 1.
By letting told -+ K-ltnew and r = 87 with ? being an arbitrary constant, (15) can be transformed to+ 80 + 827 sin0 = 8£2/K.

(26)

In the underdamped Josephson junction array, the exactly equivalent perturbation equation is used for a single
junction equation [32]. The running periodic solution in (26) is proved to be unique and globally exponentiallystable [18]. We look for a series expression of this running periodic solution 0 (t) using the Poincar6-Lindstead
mehtod as follows. First, by settingO(~t) = Oo(f2t) + 801 (~t) + . . . ,.O = I2o + 8~1 + . . .

(27)
(28)

288

H.-A. Tanaka et aL/Physica D 100 (1997) 279-300

with period T = 2 z r / ~ and the initial condition 0(0) = 0, we can obtain the following by matching terms at eachpower of 3:

O(t) = .0t + A2r sin I2t + A3?(cos ~ t -- 1) + Aalr 2 sin 2 ~ t - A47 sin S2t+ 3 A S ? Z ( c o s 2.0t -- 1) -- AS(~ + ?2)(cos ~ t -- 1) + . - .

(29)with A _= (mS2) - I and 7 -- 3 - 1 r = Kmr, and= S2/K- I;2 a33 + I ; 3 A 5 3 . . . .

(30)

In the above, ~(2n (n > 1) first becomes non-zero at the term 34 term due to the cancellation of the secular term.
Using .f2p obtained by the Melnikov method (see Appendix B), the following estimates are obtained, which will
be used later in this section:
A = (mI2) - I _< (m~2p) -1 ---- ¼7rV/~/r(31)

A27 = Ar(K/f2) < (¼7r) 2 ~ 0.62.

(32)

and

As we shall see in Sections 3.4 and 4, the coherence r takes a value, _> ¼rr as shown in Section 3.4, above critical
K. This implies that we do not have to consider small r. Thus, we only consider 8 = (Kin) -1 that satisfies 3 << rand A << 1 from (31). From the series solution of O(t) in (29), it can be verified that the coefficients of sin/cos
terms are multiplied by O(327) as the harmonics increased by one. This multiple A27 is estimated as A2r < (l:rr)2
even at ~ = I2p as in (32) and this can be much smaller for S2 within S2p < S2 < S2M. Also in (30), the second
term, --(~/2)A33, can be estimated from (31) and (32) as l r 2 A 3 3 = I ,~,,
( I rr'~3rl/2,~l/2+,v<< 1, and the multiplication
by A2~ between the higher-order sin/cos terms is verified as in (29). Thus, in the underdamped (26) with ,4 << 1,the lowest-order (harmonic) solution can be considered to be a reasonable approximation to the running periodic
solution 0 (t).
Under these conditions, we are now able to derive an analytic approximation or" r I'IIdriftas follows. By discarding
higher-order terms with A n coefficients (n > 4) in (29), a straightforward calculation (see Appendix A) leads tothe integrationf cos 0 dt ~ - 2 ( ~ r / ~ J1 (a) cos(C + A3~),(33)

0in which a = A2~/1 + A 2. ? and C = sin-I (1/~/-i- + A2). From (31) the oscillation amplitude a can be estimatedas a = A 2 ~• ~ < (¼Jr)2~/1 -4- (Tr/4)2(6/r). It should be noted that a becomes much smaller than (¼7r) 2
for I2 > f2p. Hence, Jl(a) in (33) can be approximated as Jl(a) ~ I a = A 2 ~
A < (rn~2p)-J and the above approximation, we have to lowest order
Jl (a) cos(C -f- A3r) ~'~ l r A 3• ~/2. Moreover, from(34)

(see Appendix A.) We use (34) to approximate (23) and finally obtain the following:
I2M
1,1I
1 f _2(zr/~)(~/2)A3l.(21g(I2)dJ 2 _rdrift = ~-~P,D

Kr [ _ 1A2]$2M

292M

2 ~ JY2p,D"(35)

H. -A. Tanaka et al./Physica D 1O0 (1997) 279-300

289

On the other hand, the uniform distribution of the natural frequency g (Y2) on [-S2M, I2M] is independent of S-2 andtherefore gives

g($2) = 1/212M for S'2 E [--I2M, 52M],(36)

which reduces (16) and (17), respectively, to the following forms (37) and (38):0v
r~ock -----Kr f

cos 2 0 g ( K r sin0)dO = K----f-r(0p + ½ sin20p),
2£2M

(37)

-0p
zr/2
n =Kr
rl°ck

f cos2Og(KrsinO)dO= rr Kr

2 212-----M'

(38)

-~/2

where 0p ~ sin -I [I2p/(Kr)] and g(Kr sin 0) = 1/2~M is used.
. upper
For this particular uniform g(I2) case, g($2) has a bounded support [ - I 2 m , g2M] and in the cases of K > ~c(in increasing K) and K > /(lower(in decreasing K) the contribution from the locked oscillators r lockm takes the--cform of Kr fo~,~ c°s 2 0 g ( K r sin 0) dO with sin -1 [~2M/(Kr)I.
We note that in Kuramoto's model with uniformly distributed natural frequencies, the rlock term in (11) has thesame form as (38). Thus, for Kurarnoto's model with uniformly distributed natural frequencies .(2 ~ [-S2M, S2M],
the self-consistent equation (14) for the synchronization onset can be written in the following form:

r--r, ock--rdrift=r(l

n"
= 2~M
0K " ) 2(39)

Summarizing, we have the following sets of self-consistent equations, respectively, for the increasing K case [ 1]and decreasing K case [II]:[I] increasing K case: if I2p < £2M (; K < Kcupper), r = r~ock + rlrift , then from (35) and (37)r{

~ K [ - IA2]I2M } = 0 ;+ ½ sin 20p) + 2S M
2-K

1 - 2S2 ,(40)

r2.. cos2Og(KrsinO)dO = (Kr/2I-2M)(OM+~sin2OM)andfrom,,., ---- Kr - - ~i f ~ p >_ I-2M (; K > ,,.
--cuppe,,, , r = rlock
(37)

r

1-

2~2M

(Op + ½ sin 201,-----O,(41)

where O~ = sin -I [ ~ / ( K r ) ] .
/(lower]., r = rH
[I1] decreasing K case: if $2D < $2M (; K < --clock + rlIdrift' anta~ Irom (35) and (38)
{

r

rr

1

K

--+
2 2£2M

K

12~12M }[ - 7 `4 Jot, =0;(42)( l o w e r ~,, r = flock,M leading to the same form as (41),
if I2D >_ I2M (; K > l--cr

{ 1 -- 2----~-~M
K (OM+½sin2OM)}=O.

(43)

290

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

3.4. Discontinuous phase transition at the critical point (Kc, rc)
In (40) and (42), r = 0 is just a formal solution since we assumed r >> ~ > 0 in Section 3.3. However, directlyfrom (8) and (15) the incoherent stater(t) ~ 0

Vt,

(44)

VK

is verified to be a possible (static) state be setting 0 i as follows. For any ~ , 6 [--~2M, 12M], assume n (> 2)
oscillators have this natural frequency and their initial state Oi (t = 0) satisfy ~ Oi (0) = 0, where the summation istaken over the n oscillators with S-2 = ~ . . This particular configuration of Oi leads to r(t) =- 0 and ~b = constantin t for any N. Of course, this situation is too restricted to model nearly incoherent states with certain (small)fluctuations around the above-mentioned completely incoherent state observed in numerical simulations and onlyshows the existence of the static r --- 0 solution in the limit of large N. (For a more general configuration, see thediscussion in [25].)On the other hand, for certain large K the coherence r can be determined by (41) or (43). This (K, r) curve canexist for K such that S2p > ~ M in case [I] and for 12D > S2M in case [II], respectively, as shown in Fig. l(a). The
end points of this curve correspond to I2p = 12M for case [I] and 12D = S2M for case [II]. Since r is shown to be
an increasing function of K (see Appendix C) and I2p and ~2D are, respectively, given by I2p = ¼zrq~-g-P/m and. . . upper~2D = K r , ~2p and I2D are increasing functions of K. Hence, we can determine these end points ttCc
, r upper)
for [I] and (g~ °wer, r~°wer) for [II] as follows.
For [II]+A/ ~ ~ / t,- upper~upper/ -~"~M = ~('2p = V-r/+~jW,, c"c/,,~ (end point of(K, r) curve),.. upper1 -- t C c [@ + ½ sin2@],(from (41))
212M

(45)
(46)

and

Kcu p prcuppere r sin Op = X?M(from the definition of Op).

(47).~ l o w e r l o w e r/kcrc(end point of(K, r) curve),
R"lower
1 --=c(from (43))
2 2~M

(48)

For [II]~(2M = "~D =(49)such that
Klower----_412M/Zr,
C

lower

rc

lg.

(50)

It should be noted that the above solutions (Kupper, rcUpper')and (K cl°wer, rl°wer)c lie on the same (K, r) curve defined by
the self-consistent equation r = Hock. From (46) and (49), Kg pper > K~°wet follows because 0r,+ ½ sin 2@ < ½zr for
0 < 0p < ½zr. Since the (K, r) curve is given by a continuous monotonic increasing function r ( K ) for K > --c](l°wer
as shown in Appendix C, rcupper > r clower = ¼zr follows.
We now consider the solution of the self-consistent equations if K is below -K- Cl°wer for [II] and/c'upper
for [I]. We
sLC
first assume non-zero r > 0 exists and show that this cannot be possible; the only possibility is r = 0. Given anyr > 0, the terms (K/212M)[. • -] from drifting oscillators in (40) and (42) have some non-zero value and this implies
that r is no longer on the coherent (K, r) curve for K < ](lower(/gupper). From A << 1 and (32), higher-order--C14.-A. Tanaka et al./Physica D 100 (1997) 279-300

291

terms in the (K/2~a4)[. -.] term in (40) and (42) are negligibly small and the lowest-order term is dominant. This
lowest-order term ( K / 2 ~ M ) [ - - ½ A 2] is clearly always positive, For case [II], the 1 - ½7r(K/2~M) term in (42) isalso positive because K < .i
cKl°wer = 4[2M/Tr. Also for [I], the 1 - (K/2~2M)(Oe + ½ sin 20p) term in (40) is shown. upper
to be positive for K < t~cfrom a straightforward calculation (see Appendix C). Thus, both for [I] and [II], the{. • .} terms in (40) and (42) are always positive, proving non-existence of any non-zero (static) coherence r > 0 for. upper
KlowerK < t%
a n d K <__c•
The above reasoning also applies to Kuramoto's model with uniformly distributed natural frequencies, whose
self-consistent equation for the synchronization onset is given by (39).

4. Numerical simulations
As we see in (25), in Section 3.3, we do not have to consider m and a2M (/K) separately. Thus, for numericalsimulation purpose we fix the distribution of Y2 as g2M = 5.0 and gradually change m in a region of large mK.

4.1. Numerical scheme
The numerical simulations of (5) were performed using a fourth-order Runge-Kutta integration scheme with atimestep 0.2. Several of the results were checked with a timestep 0.05. The number of oscillators usually used was500. In the numerical simulations, we choose N evenly spaced natural frequencies Y2i on the interval [ - - ~ M , ~M]
such that ~2i = --~2M + 2(i -- 1 ) ~ M / ( N -- 1). Two different types of initial conditions are considered: (a) Oi = 0
and O) i = 0, and (b) Oi uniformly distributed in [0, 27t) and wi also uniformly distributed with wi = ~i.
To obtain the (K, r) characteristics, we employed two different schemes (i) and (ii): (i) starting from (b), K isvaried up until all oscillators are coherently locked and r becomes r ~_ l, and then K is varied down without feeding
new initial condition for Oi and wi from outside; and (ii) in increasing/decreasing K, the initial conditions are given
as in (b)/(a), respectively for each update of K. The results are compared for several choices of m and a"2M, and
it is verified that the (K, r) characteristics obtained from (i) and (ii) do not show significant differences except forsmall fluctuations near the critical K of the coherence onset. Most of the numerical results presented in this section
were obtained by scheme (i).

4.2. Critical points and discontinuous jumps
Fig. 4 shows a typical example of the hysteretic synchrony observed numerically. These data were obtained using
the scheme (i) for N = 500 and m = 0.85. As we showed theoretically in Section 3, the discontinuous jumps
were observed at two different critical points of (K, r). From the enlargement of Fig. 4, these critical points are:(K, r) ~ (6.37, 0.785) and (K, r) "" (13.77, 0.966). The theoretical values from (45)-(47) and (50) lead to thecorresponding points: (K, r) ~ (6.366, 0.785) and (K, r) _~ (13.44, 0.975).
A non-stational oscillation of r was observed when we chose m larger than m ~ 0.90. For such m-values andin a range of K the oscillators have a partially coherent, oscillating state from the initial (nearly) incoherent state.
As shown in Fig. 5, this oscillating state onsets at a certain K and exists below the critical K above which all theoscillators coherently synchronize. Careful examination of the data from the enlargement of Figs. 5(a)-(c) leads tothe numerically obtained upper critical points:
for m = 0.90,(K, r) ~ (14.61,0.979),

form =0.95,( K , r ) ~ (15.32, 0.981),

for m = 1.0,

(K, r) --~ (15.92, 0.982).

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

292

1
0 .9
0 .8
0 .7

-

0 .6

-

0 •5

-

0 .4

-

0 •3

-

0 .2

-

0 .!

-

0
8

6

I0

12

14

16

K
Fig. 4. Typical example of the hysteretic synchrony; m = 0.85, N = 500 and .-QM = 5.0.

l

T

I

T

ll

12

13
K

0 9
0 8
0 7
0 6
0 5
0 4
0 .3

02?
0
I0
(a)

14

15

16

Fig. 5. Incoherent state, oscillating r state, and coherent locked state, observed with increasing coupling strength K: (a) m = 0.90; (b)m = 0.95; (c) m = 1.00.

293

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

1

I

"T

T

"f

~

ii

12

13
K

14

15

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.I
0

10
(b)

16

1
0.9

-

0.8

-

0.7

-

0.6

-

0.5

-

0.4

-

0.3

-

0.2

-

0.i
0
12

13

!

I

I

I

I

14

15

16

17

18

K

ic)
Fig. 5. Continued.

19

H.-A. Tanaka et al. IPhysica D 100 (1997) 279-300

294

1

I
,'

0.8

!

I
,

¢

I

,

" ,'I ',I,'/'"k. /If,,'
,

I

I

0.6

0.4

0.2

0

'

60

I

I

I

I

62

64

66

68

70

t
Fig. 6. Wave forms of oscillating r for rn = 0.95; K = 12.0 (solid curve), K = 13.0 (dashed curve), and K = 14.0 (dotted curve).

Solving (45)-(47) leads to the theoretical values of the upper critical points:
for m = 0.90,(K, r) ~ (14.19, 0.978),for m = 0.95,(K, r) "~ (14.95, 0.980),

for m = 1.0,

(K, r) --~ (15.70, 0.982).

These results, along with the above data for m = 0.85, show reasonable agreement between the simulations and, thetheoretical predictions, and the agreement improves as m becomes larger. Thus, we have verified the existence of•
upper
upper,
the upper critical pomt (Kc, rc) as obtained from the stationarity assumption [A 1] for the coherence r. The
oscillating r state, however, clearly breaks the assumption [A1] and therefore cannot be explained from the theory
in Section 3.

4.3. Oscillation of r
We now focus our attention on the oscillating state• Wave forms of r are shown in Fig. 6 for m = 0.95,
N = 500, for three values of K, which seem to be periodic. As we see in Figs. 5(a)-(c), the range of K showing
an oscillation of r increases as m is increased. More specifically, the oscillating region of K is defined by theregion of K values for which there is a spontaneous transition from an incoherent state to an oscillating state (upper
critical K upper) - (oscillation onset Kosc). In Fig. 7, we compare the oscillation onsets for different system size
N = 300, 500, and 700 using the scheme (i). The result shows the onset Kosc is increased as N becomes larger
upper
and Kosc seems to tend to Kc. However, this does not necessarily imply an oscillating state cannot exist, as N
gets large, for specific initial conditions• To see if such an oscillation of r can be a possible state for large N, we
numerically explore (5) for various N up to 2500, using the following initial conditions: (c) Oi and wi are given asin (b) for I~oi(0)] = ]I2i1 > 0.1, otherwise Oi(O) = 0 and O ) i ( 0 ) --- ~f2i. As shown in Fig. 8, the amplitude of theH.-A. Tanaka et al./Physica D 100 (1997) 279-300

1

I

0 .9

-

0 .8

-

0 .7

-

0 .6

-

0 .5

-

295

!

0 .4
0 .3
0 .2
0 .1

-

l0

ii

12

'!' ' ~ ' ~ ' ' - t ""l" " . . . . . . . .

13
K

14

'

15

16

Fig. 7. Oscillation onset for different s y s t e m sizes N ; N = 300 (solid curve), N = 500 (dashed curve), N ---- 700 (dotted curve).

0.5

I

I

I

I

0.45
0.4
0)

0.35
0.3

-H

0.25
0.2
0.15

",+
% ~-+- - - - + - - - -+ - - , ,:~" - ~ / ~v"""--w.........

,
"-

-

0.I
0.05

7-.iK

~,-D"-D---~~

0 ~ . u ~500.........

G ..........

I

I000

..........

I

I

1500

2000

.

2500

N
Fig. 8. Oscillation a m p l i t u d e for different s y s t e m sizes N with m = 0.95, and K = 12.0 (<>), K = 13.0 ( ÷ ) and K = 14.0 ([3).

296

H.-A. Tanaka et al./Physica D tO0 (1997) 279-300

oscillation of r does not decrease as N is increased to 2500. This result, along with the result in Fig. 7, seems to
support the following conjecture: an oscillating state can be a possible state in the large N limit for a range of mand K,for certain initial conditions. The transition from the incoherent state to this oscillating state can be triggeredby fluctuations in the (nearly) incoherent state, which increase with decreasing N.

5. Conclusion and discussion
In this paper, we studied the collective dynamics of coupled, many oscillator systems with hysteretic characteristics. To facilitate the analysis, we made two assumptions familiar from statistical mechanics: that the coupling
between oscillators is global and that the system is infinitely large. These assumptions made it possible to determine
the onset of collective synchrony in a self-consistent manner, which can be considered as a higher-order extensionof Kuramoto's theory. We have shown that, both for the extended Kuramoto's equation (5), as well as Kuramoto's
equation with uniformly distributed natural frequencies, infinitely many oscillators exhibit a discontinuous firstorder phase-transition-like change between the incoherent state and the coherent synchronized state. We have also
shown that, depending on initial conditions, this transition takes place at different critical values/(lower and/(upper•
The K~°wer value is the same as Kuramoto's, while the Kcupper value islarger than /(lower--c• Our numerical simulations
show that an oscillatory r state can exist and bifurcate at K - - - gcl°wer from the coherent state for non-small m.
Discontinuous jumps have also been found and examined analytically in a CDW model with a global coupling
structure [29] and also numerically observed even in randomly coupled neural elements [26]. A common featureof these systems is the existence of a uniformly distributed random variable in a finite interval, e.g. the random
pinning in [29] and the uniformly random natural frequencies in [26]. The question naturally arises if the uniform
distribution of the random variable is the cause of the discontinuous jumps. Comparison of a uniform distribution
with a non-uniform distribution of the natural frequencies in (5) is in progress.
Daido examined the first-order system (1) with matched, general periodic function nij (-) ~ h ( ' ) and a Lorentzian
distribution g(I2). He showed that hysteresis and bistability occurs in the "inverted bifurcation" for certain nonsinusoidal functions [6]. He called this phenomenon coupling induced hysteresis and bistability. Even for theKuramoto model with symmetric and double-peaked distribution g(Y2), such hysteresis appears [2]. Compared to
Daido's work [6], this might be calledfrequency distribution induced hysteresis. While, the hysteresis in the extended
Kuramoto model reflects the hysteretic response in each oscillator and therefore might be called adaptation induced
or inertia inducedhysteresis, although we only considered a non-generic uniformly distributed frequency distribution
in this article.
There are also several mathematical and experimetal open problem as follows.
(1) Although we have verified the oscillating state of r up to N = 2500, it is not yet clear how this state can becharacterized in the large N limit. Also, if such a state is possible in the large N limit, how can the stability ofthis state be related to m?
(2) We have shown that for a range of K the coherent locked state and the incoherent state can coexist. However,
we have not proved these are stable. Also, we have not proved whether the incoherent state is destabilized above
K > g upper or not.
(3) In the synchronous firing of a swarm of the fireflies, the strength of mutual entrainment K can be interpreted
as a ratio of the firing brightness to the background brightness. Thus, K should be varied up/down as it gets
dark/light in the laboratory experiment. It would be very interesting to observe if a swarm of P malaccae
exhibits a hysteresis in its synchrony as it gets drak/light.
(4) Recent studies on Josephson junction arrays revealed that their synchrony can be hysteretic depending on
initial condition when the second-order time derivative (/3tp) is not negligibly small in the governing equation.

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

297

Wiesenfeld pointed out that such a hysteretic synchrony may be related to the hysteresis found in the present
study. Can we explain the hysteresis in Josephson junction arrays in such a self-consistent theory? Work in thisdirection is in progress.

Acknowledgements
One of the authors (HT) deeply thanks Prof. M.A. Lieberman and Dr. M. de Sousa Vieira for their continuous
support and helpful discussions. He would like to thank Dr. S. Watanabe (Niels Bohr Institute) for many stimulating
discussions over the year and Prof. K. Wiesenfeld for his comments on the applications to the Josephson junction
arrays. Thanks are also due to Prof. S.H. Strogatz for fruitful discussions on the second-order phase model andKuramoto's theory at the Jackson Hole meeting. HT also thanks Prof. T. Endo (Meiji University) and Prof. R.
Hirota (Waseda University) for their helpful comments and encouragements. This work was initiated at the 1995
IIAS summer seminar organized by Prof. A. Hasegawa and Prof. N. Itoh. HT thanks Professors N. Itoh, K. Gohara,
M. Kikuchi, H. Nishimori, T. Ohira, T. Yanagita and H. Yuasa, and other participants for their helpful discussions.
The research of HT, was supported in part by JSPS Research Fellowships for Young Scientists.

Appendix A. Approximation of the cos 0 termThis appendix gives a calculation for the cos 0 term needed in Section 3.4. We employed the approximation of0: 0(r) : 00 + ~202 + ~303, discarding ~n(n > 4) terms. This yields a good approximation of cos0(r) because~ 4 ( K / f 2)4~2 = (A4~)2 becomes small compared to the coefficients in the lower-order terms. Then, cos 0 (~) canbe expressed by the Bessel functions as follows:cos0(r) ~ c o s [ r + A2~ sin r + A3~cosr -- A3~] = c o s i e r + a sin(~t + C + A3~)]-= Jo(a) cos ~ t - Jj (a) cos(C + A37) + J1 (a) cos(2~t + C + A3~),+ (higher harmonics)(A. l)

leading tof cos 0 (t) dt =Jl(a)cos(C + A3~),

(A.2)

0

where A = (mY2)-I , and a and C are determined froma sin(r + C) = `62f sin r + A37 cos r,
a=Ae~l+Aef,
sinC=l/~/l+,62,

cosC=,6/~/l+A

e.

(A.3)

If ,6 << 1, (A.2) can be further approximated as follows. In Section 3.3, we obtained J1 (a) -----½a = ~Aex/1 + ,62.
?/2. Using this approximation and (A.3), the Jl (a) cos(C + A37) term in (33) can be given asJl (a)cos(C + A37) ~ ½f,62 /~-_}_ A2[COS(`637)(`6/~/1 q._ A2) _ sin(`637)(1/~/1 + A2)]
= ½~[`63 Cos(A3~) _ `62 sin(,62f)]

(A.4)

298

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

As ,4 < ( m I 2 p ) - l = (4/Tr)vC~-/r) << 1 (see Appendix B) and ,43? : ,42?. ,4 << 1, COS(,437) ~ 1 - ½(,43?)2 +O([A37] 4) and sin(A3?) "- "43? + 0(['43?]3) follow. Neglecting 0([.43?] 4) terms in the cos term and 0(['43?] 3)
terms in the sin term, we finally obtain the following from (A.4):
JI ( a ) c o s ( C + ,43?) _ ½?[A 3 -- ?Zl 5 -- 172"49].

(A.5)

Appendix B. Analytic expression of S2p
This appendix gives an analytic expression of ~2p needed in Section 3.4. A new time-scale tnew = (~/-(-KT/m)told
transforms (15) to the following form:+ (1/mx/-m--K~r).0 -4- sin0 = 12/(Kr).

(B.1)

Using 1/v/-m-K-7 = 6(~/-~/r) and ~2/(Kr) = ~(12m/r), (B.1) becomes~ / = - sin0 + 6(mF2 / r - (v/-~r/r)O) = - sin0 + 3(1 - otO).

(B.2)

Eq. (B.2) has a suitable form for Melnikov's analysis, if8 << 1, I =_ J2m/r = O(1) and ct - v~r/r = O(1). The
homoclinic bifurcation curve in the (I, or) parameter space is tangent to the line I = 4ot/rr and this line is close to
the homoclinic bifurcation curve even if 3a is not so small (see [10,28]). Thus, we obtain the approximation of S2p
from I = 4ot/zr, i.e. m l 2 p / r = (4/zr)v/-~/r. The result is given as
4 K~

(B.3)

Appendix C. K, r curve defined by r ~D flockThis appendix gives an analytic expression for the (K, r) curve for K > /(lower--c, which is defined by (41) or (43).
We prove here K = K (r) is a continuous monotonic increasing function of r. Eliminating K from (43) and setting
sin Or, = ~2M/(Kr), we obtainr(0p) =

0p + (1/2) sin 20p(C.1)

2 sin Or,

Then, r'(0p) becomescos 0p (1 sin 20p - Or,),r~(0P) - 2 sin 2 0p

(C.2)

which is always negative for 0 < 0p < ½zr. Hence, r(0p) is a monotone decreasing function of 0p. Also, (43)shows that K(0p) is a monotone decreasing function of 0p. Thus, r ( K ) is a monotone increasing function of K. The
continuity of r ( K ) is clear from (43) and sin 0r, = J2M/(Kr).
We show below that the 1 - (K/212M)(Op + ½ sin 20p) term in (40) is positive. By letting K(0r, + ½ sin 20p) =f(K),df(K)dK = (0p + ½ sin 20r,) + K ( I + cos 20p) d0Pd____K_(C.3)

follows. From sin Or, = $2p/(Kr) = (4/Jr)(mrK) -1/2, the following relations are obtained:! sin 20p = (4/zr)2(mrK)-U2[(¼zr) 2 - ( m r K ) - l ] -U2,2

cos 20p = 1 - 2(4/yr )2 (mr K) -1,

H.-A. Tanaka et al./Physica D 100 (1997) 279-300

299

and
d0p = l ( m r ) _ l / 2 K _ 2 / 3 [ ( 1 ~ ) 2
dK_ (mrK)_,]-l/2.

U s i n g t h e s e i n ( C . 3 ) , d f ( K ) / d K c a n b e s i m p l i f i e d t o d f ( K ) / d K = 0p > 0. As 1--(K/212M)(Op+ 1 sin 20p) = 0
follows.
> 0, 1 - (K/212M)(Op + 1 sin 20p) > 0 for K < /(lower--cfor K = --cgl°werand d f ( K ) / d K

References
Ill Y. Aizawa, Synergetic approach to the phenomena of mode-locking in nonlinear systems, Prog. Theoret. Phys. 56 (1976) 703-716.
[21 L.L. Bonilla, J.C. Neu and R. Spigler, J. Stat. Phys. 67 (1992) 313.
[31 J.D. Crawford, Scaling and singularities in the entrainment of globally coupled oscillators, Phys. Rev. Lett. 74 (1995) 4341-4344.
[41 F. Crick and C. Koch, Toward a neurobiological theory of consciousness, Seminars Neurosci. 2 (1990) 263-275.
[5] H. Daido, Lower critical dimension for populations of oscillators with randomly distributed frequencies: A renormalization group
analysis, Phys. Rev. Lett. 61 (1988) 231-234.
[6] H. Daido, Onset of cooperative entrainment in limit-cycle oscillators with uniform all-to-all interactions: Bifurcation of the orderfunction, Physica D 91 (1996) 24-66.
[7] M. de Sousa Vieira, A.J. Lichtenberg and M.A. Lieberman, Self-synchronization of many coupled oscillators, Int. J. Bifurcation
Chaos 4 (1994) 1563-1577.
[8] G.B. Ermentrout, An adaptive model for synchrony in the firefly Pteroptyx malaccae, J. Math. Biol. 29 (1991) 571-585.
[9] C.M. Gray, P. K~nig, A.K. Engel and W. Singer, Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization
which reflect global stimulus properties, Nature 338 (1989) 334-337.
[101 J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (Springer, New
York. 1983).
[ 11 ] P. Hadley and M.R. Beasley, Dynamical states and stability of linear arrays of Josephson junctions, J. Appl. Phys. 50 (1987) 621-623.
[ 121 EE. Hanson, Pacemaker control of rhythmic flashing of fireflies, in: Cellular Pacemakers, ed. D.O. Carpenter, Vol. 2 (Wiley, New
York, 1982) pp. 81-100.
[ 13] J.J. Hopfield, Pattern recognition computation using action potential timing for stimulus representation, Nature 376 (1995) 33-36.
[141 T. Konishi and K. Kaneko, Clustered motion in symplectic coupled map systems, J. Phys. A. 25 (1992) 6283-6296.
[15] Y. Kuramoto, Self-entrainment of a population of coupled nonlinear oscillators, in: Int. Syrup. on Mathematical Problems in
Theoretical Physics, ed. H. Araki, Lecture Notes in Physics, Vol. 39 (Springer, New York, 1975) pp. 420-422.
[16] Y. Kuramoto, Chemical Oscillation, Waves and Turbulence (Springer, Berlin, 1984).
[ 171 G. Laurent and H. Davidowitz, Encoding of olfactory information with oscillating neural assemblies, Science 265 (1994) 1872-1875.
[18] M. Levi, EC. Hoppensteadt and W.L. Miranker, Dynamics of the Josephson junction, Quart. Appl. Math. (1978) 167-198.
[19] A.J. Lichtenberg and M.A. Lieberman, Regular and Chaotic Dynamics, 2nd Ed. (Springer, New York, 1992).
1201 P.C. Mattews, R.E. Mirollo and S.H. Strogatz, Dynamics of a large system of coupled nonlinear oscillators, Physica D 52 (1991)
293-331.
[211 R.E. Mirollo and S.H. Strogatz, Synchronization of pulse-coupled biological oscillators, SIAM J. Appl. Math. 50 (1990) 1645-1662.
[221 N. Nakagawa and Y. Kuramoto, From collective oscillations to collective chaos in a globally coupled oscillator systems, Physica D
75 (I 994) 74-80.
[231 E. Niebur, H.G. Schuster, D.M. Kammem and C. Koch, Oscillator-phase coupling from different two-dimensional network
connectivities, Phys. Rev. A 44 (1991) 6895-6904.
[241 H. Sakaguchi and Y. Kuramoto, A soluble active rotator model showing phase transitions via mutual entrainment, Progr. Theoret.
Phys. 76 (1986) 576-581.
1251 F.M.A. Salam, J.E. Marsden and P.P. Varaiya, Arnold diffusion in the swing equations of a power system, 1EEE Trans. CAS 32(1983) 784-796.
[26] K. Satoh, Computer experiment on the cooperative behaviour of a network of interacting nonlinear oscillators, J. Phys. Soc. Japan58 (1989) 2010-2021.
[27] S.H. Strogatz, Norbert Wiener's brain waves, in: Lecture Notes in Biomathematics, Vol. 100 (Springer, Berlin, 1994).
[28] S.H. Strogatz, Nonlinear Dynamics and Chaos (Addison-Wesley, Reading, MA, 1994) p. 273.
[29] S.H. Strogatz, C.M. Marcus, R.E. Mirollo and R.M. Westervelt, Collective dynamics of coupled oscillators with random pinning,
Physica D 36 (1989) 23-50.
130[ S.H. Strogatz and R.E. Mirollo, Phase-locking and critical phenomena in lattices of coupled nonlinear oscillators with random
intrinsic frequencies, Physica D 31 (1988) 143-168.
[311 S.H. Strogatz and R.E. Mirollo, Stability of incoherence in a population of coupled oscillators, J. Stat. Phys. 63 (1991) 613-635.

300

H.-A. Tanaka et al./Physica D 100 (1997) 279-300[32] S. Watanabe and S.H. Strogatz, Constants of motion for superconducting Josephson array, Physica D 73 (1994) 197-253.
[33] S. Watanabe and J.W. Swift, Stability of periodic solutions in series arrays of Josephson junctions with internal capacitance, J.
Nonlin. Sci., to appear.
[34] C.E. Wayne, Bounds on the trajectories of a system of weakly coupled rotators, Commun. Math. Phys. 104 (1986) 21-36.
[35] K. Wiesenfeid, P. Colet and S.H. Strogatz, Synchronization transitions in a disordered Josephson series array, Phys. Rev. Lett. 76(1996) 40n~ 497.
[36] A.T. Winfree, Biological rhythms and the behavior of populations of coupled oscillators, J. Theor. Biol. 16 (1967) 15-42.
[37] A.T. Winfree, The Geometry of Biological Time (Springer, New York, 1980).
[38] Y, Yamaguchi and H. Shimizu, Theory of self-synchronization in the presence of native frequency distribution and external noises,
Physica D 11 (1984) 212-226.
[39] R,A. York and R.C. Compton, Quasi-optical power combining using mutually synchronized oscillator arrays, 1EEE Trans. Microwave
Theor. and Tech. 39 (1991) 1000-1009.