From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators

  • Center for Applied Mathematics and Department of Theoretical and Applied Mechanics, Kimball Hall, Cornell University, Ithaca, NY 14853, USA

Abstract

The Kuramoto model describes a large population of coupled limit-cycle oscillators whose natural frequencies are drawn from some prescribed distribution. If the coupling strength exceeds a certain threshold, the system exhibits a phase transition: some of the oscillators spontaneously synchronize, while others remain incoherent. The mathematical analysis of this bifurcation has proved both problematic and fascinating. We review 25 years of research on the Kuramoto model, highlighting the false turns as well as the successes, but mainly following the trail leading from Kuramoto’s work to Crawford’s recent contributions. It is a lovely winding road, with excursions through mathematical biology, statistical physics, kinetic theory, bifurcation theory, and plasma physics.

Keywords

  • Kuramoto model;
  • Coupled oscillators;
  • Kinetic theory;
  • Plasma physics

1. Introduction

In the 1990s, Crawford wrote a series of papers about the Kuramoto model of coupled oscillators [1], [2] and [3]. At first glance, the papers look technical, maybe even a bit intimidating.

For instance, take a look at “Amplitude expansions for instabilities in populations of globally coupled oscillators”, his first paper on the subject [1]. Here, Crawford racks up 200 numbered equations as he calmly plows through a center manifold calculation for a nonlinear partial integro-differential equation.

Technical, yes, but a technical tour de force. Beneath the surface, there is a lot at stake. In his modest, methodical way, Crawford illuminated some problems that had appeared murky for about two decades.

My goal here is to set Crawford’s work in context and to give a sense of what he accomplished. The larger setting is the story of the Kuramoto model [4], [5], [6], [7], [8] and [9]. It is an ongoing tale full of twists and turns, starting with Kuramoto’s ingenious analysis in 1975 (which raised more questions than it answered) and culminating with Crawford’s insights. Along the way, I will point out some problems that remain unsolved to this day, and tell a few stories about the various people who have worked on the Kuramoto model, including how Crawford himself got hooked on it.

2. Background

The Kuramoto model was originally motivated by the phenomenon of collective synchronization, in which an enormous system of oscillators spontaneously locks to a common frequency, despite the inevitable differences in the natural frequencies of the individual oscillators [10], [11], [12] and [13]. Biological examples include networks of pacemaker cells in the heart [14] and [15]; circadian pacemaker cells in the suprachiasmatic nucleus of the brain (where the individual cellular frequencies have recently been measured for the first time [16]); metabolic synchrony in yeast cell suspensions [17] and [18]; congregations of synchronously flashing fireflies [19] and [20]; and crickets that chirp in unison [21]. There are also many examples in physics and engineering, from arrays of lasers [22] and [23] and microwave oscillators [24] to superconducting Josephson junctions [25] and [26].

Collective synchronization was first studied mathematically by Wiener [27] and [28], who recognized its ubiquity in the natural world, and who speculated that it was involved in the generation of alpha rhythms in the brain. Unfortunately Wiener’s mathematical approach based on Fourier integrals [27] has turned out to be a dead end.

A more fruitful approach was pioneered by Winfree [10] in his first paper, just before he entered graduate school. He formulated the problem in terms of a huge population of interacting limit-cycle oscillators. As stated, the problem would be intractable, but Winfree intuitively recognized that simplifications would occur if the coupling were weak and the oscillators nearly identical. Then one can exploit a separation of timescales: on a fast timescale, the oscillators relax to their limit cycles, and so can be characterized solely by their phases; on a long timescale, these phases evolve because of the interplay of weak coupling and slight frequency differences among the oscillators. In a further simplification, Winfree supposed that each oscillator was coupled to the collective rhythm generated by the whole population, analogous to a mean-field approximation in physics. His model is

View the MathML source
where θi denotes the phase of oscillator i and ωi its natural frequency. Each oscillator j exerts a phase-dependent influence X(θj) on all the others; the corresponding response of oscillator i depends on its phase θi, through the sensitivity function Z(θi).

Using numerical simulations and analytical approximations, Winfree discovered that such oscillator populations could exhibit the temporal analog of a phase transition. When the spread of natural frequencies is large compared to the coupling, the system behaves incoherently, with each oscillator running at its natural frequency. As the spread is decreased, the incoherence persists until a certain threshold is crossed — then a small cluster of oscillators suddenly freezes into synchrony.

This cooperative phenomenon apparently made a deep impression on Kuramoto. As he wrote in a paper with his student Nishikawa ([8], p. 570):

“…Prigogine’s concept of time order [29], which refers to the spontaneous emergence of rhythms in nonequilibrium open systems, found its finest example in this transition phenomenon… It seems that much of fresh significance beyond physiological relevance could be derived from Winfree’s important finding (in 1967) after our experience of the great advances in nonlinear dynamics over the last two decades.”

Kuramoto himself began working on collective synchronization in 1975. His first paper on the topic [4] was a brief note announcing some exact results about what would come to be called the Kuramoto model. In later years, he would keep wrestling with that analysis, refining and clarifying the presentation each time, but also raising thorny new questions too [5], [6], [7], [8] and [9].

3. Kuramoto model

3.1. Governing equations

Kuramoto [5] put Winfree’s intuition about phase models on a firmer foundation. He used the perturbative method of averaging to show that for any system of weakly coupled, nearly identical limit-cycle oscillators, the long-term dynamics are given by phase equations of the following universal form:

View the MathML source
The interaction functions Γij can be calculated as integrals involving certain terms from the original limit-cycle model (see Section 5.2 of [5] for details).

Even though the reduction to a phase model represents a tremendous simplification, these equations are still far too difficult to analyze in general, since the interaction functions could have arbitrarily many Fourier harmonics and the connection topology is unspecified — the oscillators could be connected in a chain, a ring, a cubic lattice, a random graph, or any other topology.

Like Winfree, Kuramoto recognized that the mean-field case should be the most tractable. The Kuramoto model corresponds to the simplest possible case of equally weighted, all-to-all, purely sinusoidal coupling:

View the MathML source
where K ≥ 0 is the coupling strength and the factor 1/N ensures that the model is well behaved as N → ∞.

The frequencies ωi are distributed according to some probability density g(ω). For simplicity, Kuramoto assumed that g(ω) is unimodal and symmetric about its mean frequency Ω, i.e., g(Ω + ω) = g(Ω − ω) for all ω, like a Gaussian distribution. Actually, thanks to the rotational symmetry in the model, we can set the mean frequency to Ω = 0 by redefining θi → θi + Ωt for all i, which corresponds to going into a rotating frame at frequency Ω. This leaves the governing equations

equation3.1
View the MathML source
invariant, but effectively subtracts Ω from all the ωi and therefore shifts the mean of g(ω) to zero. So from now on,
g(ω)=g(−ω)
for all ω, and the ωi denote deviations from the mean frequency Ω. We also suppose that g(ω) is nowhere increasing on [0,∞), in the sense that g(ω) ≥ g(v) whenever ω ≤ v; this formalizes what we mean by “unimodal”.

3.2. Order parameter

To visualize the dynamics of the phases, it is convenient to imagine a swarm of points running around the unit circle in the complex plane. The complex order parameter [5]

equation3.2
View the MathML source
is a macroscopic quantity that can be interpreted as the collective rhythm produced by the whole population. It corresponds to the centroid of the phases. The radius r(t) measures the phase coherence, and ψ(t) is the average phase (Fig. 1).

Full-size image (1 K)
Fig. 1. 

Geometric interpretation of the order parameter (3.2). The phases θj are plotted on the unit circle. Their centroid is given by the complex number r eiψ, shown as an arrow.

For instance, if all the oscillators move in a single tight clump, we have r ≈ 1 and the population acts like a giant oscillator. On the other hand, if the oscillators are scattered around the circle, then r ≈ 0; the individual oscillations add incoherently and no macroscopic rhythm is produced.

Kuramoto noticed that the governing equation

View the MathML source
can be rewritten neatly in terms of the order parameter, as follows. Multiply both sides of the order parameter equation by e−iθi to obtain
View the MathML source
Equating imaginary parts yields
View the MathML source
Thus (3.1) becomes
equation3.3
View the MathML source
In this form, the mean-field character of the model becomes obvious. Each oscillator appears to be uncoupled from all the others, although of course they are interacting, but only through the mean-field quantities r and ψ. Specifically, the phase θi is pulled toward the mean phase ψ, rather than toward the phase of any individual oscillator. Moreover, the effective strength of the coupling is proportional to the coherence r. This proportionality sets up a positive feedback loop between coupling and coherence: as the population becomes more coherent, r grows and so the effective coupling Kr increases, which tends to recruit even more oscillators into the synchronized pack. If the coherence is further increased by the new recruits, the process will continue; otherwise, it becomes self-limiting. Winfree [10] was the first to discover this mechanism underlying spontaneous synchronization, but it stands out especially clearly in the Kuramoto model.

3.3. Simulations

If we integrate the model numerically, how does r(t) evolve? For concreteness, suppose we fix g(ω) to be a Gaussian or some other density with infinite tails, and vary the coupling K. Simulations show that for all K less than a certain threshold Kc, the oscillators act as if they were uncoupled: the phases become uniformly distributed around the circle, starting from any initial condition. Then r(t) decays to a tiny jitter of size O(N−1/2), as expected for any random scatter of N points on a circle (Fig. 2).

Full-size image (1 K)
Fig. 2. 

Schematic illustration of the typical evolution of r(t) seen in numerical simulations of the Kuramoto model (3.1).

But when K exceeds Kc, this incoherent state becomes unstable and r(t) grows exponentially, reflecting the nucleation of a small cluster of oscillators that are mutually synchronized, thereby generating a collective oscillation. Eventually r(t) saturates at some level r < 1, though still with O(N−1/2) fluctuations.

At the level of the individual oscillators, one finds that the population splits into two groups: the oscillators near the center of the frequency distribution lock together at the mean frequency Ω and co-rotate with the average phase ψ(t), while those in the tails run near their natural frequencies and drift relative to the synchronized cluster. This mixed state is often called partially synchronized. With further increases in K, more and more oscillators are recruited into the synchronized cluster, and r grows as shown in Fig. 3.

Full-size image (1 K)
Fig. 3. 

Dependence of the steady-state coherence r on the coupling strength K.

The numerics further suggest that r depends only on K, and not on the initial condition. In other words, it seems there is a globally attracting state for each value of K.

3.4. Puzzles

These numerical results cry out for explanation. A good theory should provide formulas for the critical coupling Kc and for the coherence r(K) on the bifurcating branch. The theory should also explain the apparent stability of the zero branch below threshold and the bifurcating branch above threshold. Ideally, one would like to formulate and prove global stability results, since the numerical simulations give no hint of any other attractors beyond those seen here. Even more ambitiously, can one formulate and prove some convergence results as N → ∞?

As we will see below, the first few of these problems have been solved, while the rest remain open. Specifically, Kuramoto derived exact results for Kc and r(K), Mirollo and I solved the linear stability problem for the zero branch, and Crawford extended those results to the weakly nonlinear case. But we still do not know how to show that the bifurcating branch is linearly stable along its entire length (if it truly is), and nobody has even touched the problems of global stability and convergence.

4. Kuramoto’s analysis

In his earliest work, Kuramoto analyzed his model without the benefit of simulations — he guessed the correct long-term behavior of the solutions in the limit N → ∞, using symmetry considerations and marvelous intuition. Specifically, he sought steady solutions, where r(t) is constant and ψ(t) rotates uniformly at frequency Ω. By going into the rotating frame with frequency Ω and choosing the origin of this frame correctly, one can set ψ  0 without loss of generality.

Then the governing equation (3.3) becomes

equation4.1
View the MathML source
Since r is assumed constant in (4.1), all the oscillators are effectively independent — that is the beauty of steady solutions. The strategy now is to solve for the resulting motions of all the oscillators (which will depend on r as a parameter). These motions in turn imply values for r and ψ which must be consistent with the values originally assumed. This self-consistency condition is the key to the analysis.

The solutions of (4.1) exhibit two types of long-term behavior, depending on the size of |ωi| relative to Kr. The oscillators with |ωi| ≤ Kr approach a stable fixed point defined implicitly by

equation4.2
View the MathML source
where View the MathML source. These oscillators will be called “locked” because they are phase-locked at frequency Ω in the original frame. In contrast, the oscillators with |ωi| > Kr are “drifting” — they run around the circle in a nonuniform manner, accelerating near some phases and hesitating at others, with the inherently fastest oscillators continually lapping the locked oscillators, and the slowest ones being lapped by them. The locked oscillators correspond to the center of g(ω) and the drifting oscillators correspond to the tails, as expected.

At this stage, Kuramoto has neatly explained why the population splits into two groups. But before we get too complacent, notice that the existence of the drifting oscillators would seem to contradict the original assumption that r and ψ are constant. How can the centroid of the population remain constant with all those drifting oscillators buzzing around the circle?

Kuramoto deftly avoided this problem by demanding that the drifting oscillators form a stationary distribution on the circle. Then the centroid stays fixed even though individual oscillators continue to move. Let ρ(θω) dθ denote the fraction of oscillators with natural frequency ω that lie between θ and θ + dθ. Stationarity requires that ρ(θω) be inversely proportional to the speed at θ; oscillators pile up at slow places and thin out at fast places on the circle. Hence

equation4.3
View the MathML source
The normalization constant C is determined by View the MathML source for each ω, which yields
View the MathML source
Next, we invoke the self-consistency condition: the constant value of the order parameter must be consistent with that implied by (3.2). Using angular brackets to denote population averages, we have
View the MathML source
Since ψ = 0 by assumption, 〈eiθ〉 = r eiψ = r. Thus,
View the MathML source
We evaluate the locked contribution first. In the locked state, sin θ* = ω/Kr for all |ω| ≤ Kr. As N → ∞, the distribution of locked phases is symmetric about θ = 0 because g(ω) = g(−ω); there are just as many oscillators at θ* as at −θ*. Hence 〈sin θlock = 0 and
View the MathML source
where θ(ω) is defined implicitly by (4.2). Changing variables from ω to θ yields
View the MathML source
Now, consider the drifting oscillators. They contribute
View the MathML source
It turns out that this integral vanishes. This follows from g(ω) = g(−ω) and the symmetry ρ(θ + π, −ω) = ρ(θω) implied by (4.3).

Therefore, the self-consistency condition reduces to

equation4.4
View the MathML source
Eq. (4.4) always has a trivial zero solutionr = 0, for any value of K. This corresponds to a completely incoherent state with ρ(θω) = 1/2π for all θ, ω. A second branch of solutions, corresponding to partially synchronized states, satisfies
equation4.5
View the MathML source
This branch bifurcates continuously from r = 0 at a value K = Kc obtained by letting r → 0+ in (4.5). Thus,
View the MathML source
which is Kuramoto’s exact formula for the critical coupling at the onset of collective synchronization. By expanding the integrand in (4.5) in powers of r, we find that the bifurcation is supercritical if g″(0) < 0 (the generic case for smooth, unimodal, even densities g(ω)) and it is subcritical if g″(0) > 0. Near onset, the amplitude of the bifurcating branch obeys the square-root scaling law:
equation4.6
View the MathML source
where
View the MathML source
is the normalized distance above threshold. For the special case of a Lorentzian or Cauchy density
equation4.7
View the MathML source
Kuramoto [4] and [5] integrated (4.5) exactly to obtain
View the MathML source
for all K ≥ Kc. This formula was later shown to match the results of numerical simulations [6] and [7].

5. Two unsolved problems

5.1. Finite-N fluctuations

In the last of her three Bowen lectures at Berkeley in 1986, Kopell pointed out that Kuramoto’s argument contained a few intuitive leaps that were far from obvious — in fact, they began to seem paradoxical the more one thought about them — and she wondered whether one could prove some theorems that would put the analysis on firmer footing. In particular, she wanted to redo the analysis rigorously for large but finite N, and then prove a convergence result as N → ∞.

But it would not be easy. Whereas Kuramoto’s approach had relied on the assumption that r was strictly constant, Kopell emphasized that nothing like that could be strictly true for any finite N. Think about the simple case K = 0. Then View the MathML source and every trajectory is dense on the N-torus, at least for the generic case where the frequencies are rationally independent. But then r(t) eventually passes through every possible value between 0 and 1, completely unlike the constant value r  0 implied by Kuramoto’s argument! Admittedly, r(t) would spend nearly all its time very close to zero, at r = O(N−1/2) ⪡ 1, and only blip up extremely rarely — in that sense r  0 is practically correct. But how can this rough idea be made precise? When K ≠ 0, the situation would become still more difficult, because now there would be three subpopulations of oscillators — locked and drifting ones as in Kuramoto’s analysis, but also some fuzzy oscillators between them, determined by the ever-fluctuating boundary ωi ≈ Kr(t).

Kopell’s suggestion was to try to prove something like this: For large N, for most initial conditions, and for most realizations of the ωi, the coherence r(t) approaches the Kuramoto value r(K) and stays within O(N−1/2) of it for a large fraction of the time. Around the same time, Daido [30], [31], [32] and [33], and Kuramoto and Nishikawa [8] and [9] began exploring the finite-N fluctuations using computer simulations and physical arguments. It appears that the fluctuations are indeed O(N−1/2) except very close to Kc, where they may be amplified [30], [31], [32] and [33].

Still, the issue of fluctuations remains wide open mathematically. As of March 2000, there are no rigorous convergence results about the finite-N behavior of the Kuramoto model.

5.2. Stability

The other major issue left unresolved by Kuramoto’s analysis concerns the stability of the steady solutions. It was in this arena that Crawford ultimately contributed so much, and so we will focus on it for the rest of this paper. Kuramoto was well aware of the stability problem; he writes [5] (p. 74):

“One may expect that negative μ (i.e., weaker coupling) makes the zero solution stable, and positive μ (i.e., stronger coupling) unstable. Surprisingly enough, this seemingly obvious fact seems difficult to prove. The difficulty here comes from the fact that an infinitely large number of phase configurations {θi,i=1,…,N} belong to an identical “macroscopic” state specified by a given value of r.”

He also remarks that it “appears to be difficult to prove” that the branch of partially synchronized states is stable when the bifurcation is supercritical, and unstable when it is subcritical.

6. Stability theories of Kuramoto and Nishikawa

Kuramoto and Nishikawa [8] and [9] were the first to tackle the stability problem. They proposed two different theories, both based on plausible physical reasoning, but neither of which ultimately turned out to be correct. Nevertheless, it is interesting to look back at their pioneering ideas, partly because they came tantalizingly close to the truth, and partly to remind us how subtle the stability problem appeared at the time.

6.1. First theory

In their first approach, Kuramoto and Nishikawa [8] tried to derive an evolution equation for r(t) in closed form, a dynamical extension of the earlier self-consistency equation (4.5). The hope was that this might be possible close to the bifurcation, where r(t) would be expected to evolve extremely slowly compared to the relaxation time of the individual oscillators. Then each oscillator would follow the order parameter almost adiabatically, allowing these rapid variables to be eliminated and causing a great reduction in the dynamics.

To push this strategy through, Kuramoto and Nishikawa [8] made several approximations whose validity was uncertain. As in the steady-state theory, they separated the population into locked and drifting groups; such a sharp division should be possible if r(t) varies slowly enough. The characteristic timescale of the locked oscillators was argued to be of order (Kr)−1, which is very slow since r(t) ⪡ 1 near the bifurcation. The theory also suggested that the drifting oscillators make a negligible contribution to the dynamics of r(t).

In the end, they were led to the following unconventional equation (see Eq. (3.36) in [8]):

equation6.1
View the MathML source
where ξs is an O(1) constant that arises in their theory, μ = (K − Kc)/Kc as before, and View the MathML source. Note the peculiar extra factor of r on the right-hand side as compared to the usual normal form near a pitchfork bifurcation. Eq. (6.1) predicts that the zero solution is stable below threshold (μ < 0), but with anomalously slow algebraic decay
View the MathML source
as t → ∞. Above threshold, the zero solution is unstable, though weakly so: r(t) initially grows only linearly in t, then eventually relaxes exponentially fast to View the MathML source.

6.2. Second theory

Kuramoto and Nishikawa soon realized that something was wrong. Two years later, they revisited the problem [9] and stated with admirable candor, “In the past, we seem to have held an erroneous view about the onset of collective oscillation…”. They now believed that the drifting oscillators are not negligible throughout the whole evolution of r(t) — rather, these oscillators play a decisive dynamical role in the earliest stages, thanks to their rapid response to fluctuations in r(t), though in the long run they still do not affect the steady value of r.

Kuramoto and Nishikawa [9] also proposed a new strategy for deriving an evolution equation for r(t). In the governing equation

View the MathML source
they pretend that r(t) is an external force, say h(t), and then derive the responses of the individual oscillators to h(t), restricting attention to the linear regime where h(t) ⪡ 1. These individual responses (which depend on the whole history of h(t)) can then be combined to yield the response of r(t). On general grounds, and without giving a derivation, Kuramoto and Nishikawa [9] guessed that r(t) should be a linear functional of h(t) of the form
View the MathML source
where M is a memory function to be determined. But since h is really r in disguise, the equation must be
equation6.2
View the MathML source
To calculate the kernel M, they consider the response to a step function
View the MathML source
and find that, for example, M(t) = et when the distribution is the Lorentzian g(ω) = [π(ω2 + 1)]−1. (The calculation of M is straightforward. The oscillators are initially distributed according to the stationary density ρ(θω) found in Section 4, where h0 plays the role of r in the earlier formulas. The density ρ is smooth in θ for the drifting oscillators and a delta function in θ for the locked oscillators. Then, since h(t) = 0 for t > 0, all the oscillators and their corresponding densities rotate rigidly and independently at their natural frequencies. The corresponding evolution of r(t) can be found by integrating eiθ with respect to these rotating densities, weighted by g(ω), and then M(t) can be extracted from the result.)

Within this revised framework, Kuramoto and Nishikawa [9] now found that r(t) grows exponentially above threshold, and decays exponentially below threshold. In other words, the zero solution was now predicted to change stability in the most standard way — it goes from linearly stable to linearly unstable as K increases through Kc.

But, should one really believe this prediction? Remember, the integral equation (6.2) was not derived in any systematic way from the governing equation (3.1). On the other hand, the intuitive argument for (6.2) looked plausible, and maybe even convincing.

7. Continuum limit of the Kuramoto model

It was against this confusing backdrop that Mirollo and I began thinking about the stability problem. At the time, it was unclear how to formulate the problem mathematically. We did not even know how to write down an infinite-N version of the Kuramoto model, let alone analyze the stability of its steady solutions.

We eventually realized that the continuum limit should be phrased in terms of densities, just as in traffic flow, kinetic theory, or fluid mechanics [34]. For each natural frequency ω, imagine a continuum of oscillators distributed on the circle. Let ρ(θtω) dθ denote the fraction of these oscillators that lie between θ and θ + dθ at time t. Then ρ is nonnegative, 2π-periodic in θ, and satisfies the normalization

equation7.1
View the MathML source
for all t and ω. The evolution of ρ is governed by the continuity equation
equation7.2
View the MathML source
which expresses conservation of oscillators of frequency ω. Here the velocity v(θtω) is interpreted in an Eulerian sense as the instantaneous velocity of an oscillator at position θ, given that it has natural frequency ω. From (3.3), that velocity is
equation7.3
View the MathML source
where r(t) and ψ(t) are now given by
equation7.4
View the MathML source
which follows from the law of large numbers applied to (3.2). Equivalently, these equations can be combined to yield a single equation for ρ in closed form:
equation7.5
View the MathML source
The expression in parentheses is v(θtω), written as the infinite-N version of (3.1).

Eq. (7.5) is the continuum limit of the Kuramoto model [34]. It is a nonlinear partial integro-differential equation for ρ. The virtue of (7.5) is that all questions about existence, stability, and bifurcation of various kinds of solutions can now be addressed systematically.

For instance, the stationary states of (7.5) are precisely the steady solutions that Kuramoto [4] and [5] wrote down intuitively. To see this, note that ∂ρ/∂t = 0 implies ρv = C(ω), where C(ω) is constant with respect to θ. If C(ω) ≠ 0, we recover the stationary density (4.3) for the drifting oscillators; if C(ω) = 0, we find that ρ is a delta function in θ, based at the locked phase found earlier.

The simplest state is the uniform incoherent state

View the MathML source
or what we earlier called the zero solution. As we will see in Section 8, its linear stability properties turn out to be stranger than anyone had expected.

Eqs. (7.2)–(7.5) had been studied previously by Sakaguchi [35], who extended the Kuramoto model to allow rapid stochastic fluctuations in the natural frequencies. The governing equations are

equation7.6
View the MathML source
where the variables ξi(t) are independent white noise processes that satisfy
View the MathML source
Here D ≥ 0 is the noise strength and the angular brackets denote an average over realizations of the noise. Sakaguchi argued intuitively that since (7.6) is a system of Langevin equations with mean-field coupling, as N → ∞ the density ρ(θtω) should satisfy the Fokker–Planck equation
equation7.7
View the MathML source
where v(θtω), r(t), and ψ(t) are given by (7.3) and (7.4). Thus Sakaguchi’s Fokker–Planck equation reduces to the continuum limit of the Kuramoto model when D = 0.

However, Sakaguchi [35] did not present a stability analysis of his model. Instead he solved for the stationary densities, and then extended Kuramoto’s self-consistency argument to determine where a branch of partially synchronized states bifurcates from the incoherent state. In this way he showed that the critical coupling is

equation7.8
View the MathML source
which reduces to Kuramoto’s formula Kc = 2/πg(0) as D → 0+.

8. Stability of the incoherent state

The linear stability problem for the incoherent state of Sakaguchi’s model was solved in [34]. Here is an outline of the approach and the results (for consistency with the rest of this paper, we will restrict attention to the Kuramoto model, where D = 0). Let

equation8.1
View the MathML source
where ε  1 and we write the perturbation η as a Fourier series in θ:
equation8.2
View the MathML source
Here c.c. denotes complex conjugate, and η contains the second and higher harmonics of η. (Note that η automatically has zero mean, because of (7.1).) We write the perturbation in this way because it turns out that the linearized amplitude equation for the first harmonic, c(tω), is the only one with nontrivial dynamics; that’s essentially because of the pure sinusoidal coupling in the Kuramoto model. Substituting for ρ into (7.5) yields
equation8.3
View the MathML source
The right-hand side of (8.3) defines a linear operator A, given by
equation8.4
View the MathML source

The spectrum of A has both continuous and discrete parts, as shown in [34]. Its continuous spectrum is pure imaginary, {iω: ω∈support(g)}, corresponding to a continuous family of neutral modes. These modes can be understood intuitively by imagining an initial perturbation η(θωt = 0) supported on a sliver of exactly one frequency, say ω = ω0. In other words, we disturb the slice of the oscillator population with intrinsic frequency ω0 and leave the rest alone in their perfectly incoherent state. The corresponding amplitude c(0, ω) would then take the form c(0, ω) = 0 for all ω ≠ ω0 (oscillators at those frequencies are not disturbed). As for ω = ω0, we can choose c(0, ω0) = 1 without loss of generality, since (8.4) is linear. The key point is that the integral in (8.4) vanishes for this strange sliver perturbation, and so (8.4) reduces to Ac = iω0c. Hence, c(0, ω) is (morally speaking) an eigenfunction with pure imaginary eigenvalue iω0, and that explains the form of the continuous spectrum. Of course, this argument is not strictly correct, because this sliver perturbation is equivalent in L2 to the zero perturbation, and so is not a valid eigenmode. But the intuition is right, and it agrees with the rigorous calculations given in [34].