Exact Results for the Kuramoto Model with a Bimodal Frequency Distribution
Abstract
We analyze a large system of globally coupled phase oscillators whose natural frequencies are bimodally distributed. The dynamics of this system has been the subject of longstanding interest. In 1984 Kuramoto proposed several conjectures about its behavior; ten years later, Crawford obtained the first analytical results by means of a local center manifold calculation. Nevertheless, many questions have remained open, especially about the possibility of global bifurcations. Here we derive the system’s stability diagram for the special case where the bimodal distribution consists of two equally weighted Lorentzians. Using an ansatz recently discovered by Ott and Antonsen, we show that in this case the infinitedimensional problem reduces exactly to a flow in four dimensions. Depending on the parameters and initial conditions, the longterm dynamics evolves to one of three states: incoherence, where all the oscillators are desynchronized; partial synchrony, where a macroscopic group of phaselocked oscillators coexists with a sea of desynchronized ones; and a standing wave state, where two counterrotating groups of phaselocked oscillators emerge. Analytical results are presented for the bifurcation boundaries between these states. Similar results are also obtained for the case in which the bimodal distribution is given by the sum of two Gaussians.
pacs:
05.45.XtI Introduction
i.1 Background
Large systems consisting of many coupled oscillatory units occur in a wide variety of situations [[1]]. Thus the study of the behaviors that such systems exhibit has been an active and continuing area of research. An important early contribution in this field was the introduction in 1975 by Kuramoto [[2], [3]] of a simple model which illustrates striking features of such systems. Kuramoto employed two key simplifications in arriving at his model: (i) the coupling between units was chosen to be homogeneous and alltoall (i.e., ‘global’), so that each oscillator would have an equal effect on all other oscillators; and (ii) the oscillator states were solely described by a phase angle , so that their uncoupled dynamics obeyed the simple equation , where is the intrinsic natural frequency of oscillator , is the number of oscillators, and . The natural frequencies are, in general, different for each oscillator and are assumed to be drawn from some prescribed distribution function .
Much of the research on the Kuramoto model has focused on the case where is unimodal (for reviews of this literature, see [[4], [5], [6]]). Specifically, is usually assumed to be symmetric about a maximum at frequency and to decrease monotonically and continuously to zero as increases. In that case, it was found that as the coupling strength between the oscillators increases from zero in the large limit, there is a continuous transition at a critical coupling strength . For below , the average macroscopic, timeasymptotic behavior of the system is such that the oscillators in the system behave incoherently with respect to each other, and an order parameter (defined in Sec. II) is correspondingly zero. As increases past , the oscillators begin to influence each other in such a way that there is collective global organization in the phases of the oscillators, and the timeasymptotic order parameter assumes a nonzero constant value that increases continuously for [[3], [4], [5], [6], [7]].
It is natural to ask how these results change if other forms of are considered. In this paper we will address this question for what is perhaps the simplest choice of a nonunimodal frequency distribution: we consider a distribution that has two peaks [[8], [9]] and is the sum of two identical unimodal distributions , such that . We find that this modification to the original problem introduces qualitatively new behaviors. As might be expected, this problem has been previously addressed [[3], [10]]. However, due to its difficulty, the problem was not fully solved, and, as we shall show, notable features of the behavior were missed.
i.2 Reduction method
The development that makes our analysis possible is the recent paper of Ott and Antonsen [[11]]. Using the method proposed in Ref. [[11]] we reduce the original problem formulation from an integropartialdifferential equation [[4], [5], [7]] for the oscillator distribution function (a function of and ) to a system of just a few ordinary differential equations (ODEs). Furthermore, we analyze the reduced ODE system to obtain its attractors and the bifurcations they experience with variation of system parameters.
The reduced ODE system, however, represents a special restricted class of all the possible solutions of the original full system [[11]]. Thus a concern is that the reduced system might miss some of the actual system behavior. In order to check this, we have done numerical solutions of the full system. The result is that, in all cases tested, the timeasymptotic attracting behavior of the full system and the observed attractor bifurcations are all contained in, and are quantitatively described by, our ODE formulation. Indeed a similar result applies for the application of the method of Ref. [[11]] to the original Kuramoto model with unimodally distributed frequencies [[2], [3]] and to the problem of the forced Kuramoto model with periodic drive [[12]].
On the other hand, the reduction method has not been mathematically proven to capture all the attractors, for any of the systems to which it has been applied [[11], [12]]. Throughout this paper we operate under the assumption (based on our numerical evidence) that the reduction method is reliable for the bimodal Kuramoto model. But we caution the reader that in general the situation is likely to be subtle and systemdependent; see Sec. VI.4.1 for further discussion of the scope and limits of the reduction method.
i.3 Outline of the paper
The organization of this paper is as follows. In Sec. II we formulate the problem and reduce it to the abovementioned ODE description for the case where is a sum of CauchyLorentz distributions.
Sec. III provides an analysis of the ODE system. The main results of Sec. III are a delineation of the different types of attractors that can exist, the regions of parameter space that they occupy (including the possibility of bistability and hysteresis), and the types of bifurcations that the attractors undergo.
In Sec. IV, we establish that the attractors of the ODEs obtained in Section III under certain symmetry assumptions are attractors of the full ODE system. In Section V, we confirm that these attractors and bifurcations are also present in the original system. In addition, we investigate the case where is a sum of Gaussians, rather than CauchyLorentz distributions. We find that the attractors and bifurcations in the Lorentzian case and in the Gaussian case are of the same types and that parameter space maps of the different behaviors are qualitatively similar for the two distributions.
Ii Governing Equations
ii.1 Problem definition
We study the Kuramoto problem of oscillators with natural frequencies ,
(1) 
where are the phases of each individual oscillator and is the coupling strength. We study this system in the limit for the case in which the distribution of natural frequencies is given by the sum of two Lorentzian distributions:
(2) 
Here is the width parameter (halfwidth at halfmaximum) of each Lorentzian and are their center frequencies, as displayed in Fig. 1. A more physically relevant interpretation of is as the detuning in the system (proportional to the separation between the two center frequencies).
Note that we have written the distribution so that it is symmetric about zero; this can be achieved without loss of generality by going into a suitable rotating frame.
ii.2 Derivation
In the limit where , Eq. (1) can be written in a continuous formulation [[3], [4], [5]] in terms of a probability density . Here is defined such that at time , the fraction of oscillators with phases between and and natural frequencies between and is given by . Thus
(3) 
and
(4) 
by definition of .
The evolution of is given by the continuity equation describing the conservation of oscillators:
(5) 
where is the angular velocity of the oscillators. From Eq. (1), we have
(6) 
Following Kuramoto, we define a complex order parameter
(7) 
whose magnitude characterizes the degree to which the oscillators are bunched in phase, and describes the average phase angle of the oscillators. Expressing the velocity (6) in terms of we obtain
(8)  
(9) 
where the * denotes complex conjugate.
Following Ott and Antonsen [[11]], we now restrict attention to a special class of density functions. By substituting a Fourier series of the form
(10) 
where ‘c.c’ stands for the complex conjugate of the preceeding term, and imposing the ansatz that
(11) 
we obtain
(12) 
where
(13) 
We now consider solutions of (12) and (13) for initial conditions that satisfy the following additional conditions: (i) ; (ii) is analytically continuable into the lower half plane ; and (iii) as . If these conditions are satisfied for , then, as shown in [[11]], they continue to be satisfied by as it evolves under Eqs. (12) and (13). Expanding in partial fractions as
we find it has four simple poles at . Evaluating (13) by deforming the integration path from the real axis to , the order parameter becomes
(14) 
where
(15) 
Substitution of this expression into (12) yields two coupled complex ODEs, describing the evolution of two ‘sub’order parameters,
(16)  
(17)  
where we use dots to represent the time derivative from now on. (This system agrees with the results of [[11]] for the case of two equal groups of oscillators with uniform coupling strength and average frequencies and .)
ii.3 Reductions of the system
The system derived so far is fourdimensional. If we introduce polar coordinates and define the phase difference , the dimensionality can be reduced to three:
(18)  
(19)  
(20) 
To facilitate our analysis we now look for solutions of Eqs. (1820) that satisfy the symmetry condition
(21) 
In Sec. IV we will verify that these symmetric solutions are stable to perturbations away from the symmetry manifold and that the attractors of Eqs. (16, 17) lie within this manifold.
Our analysis of the problem thus reduces to a study in the phase plane:
(22)  
(23) 
Iii Bifurcation Analysis
Figure 2 summarizes the results of our analysis of Eqs. (22, 23). We find that three types of attractors occur: the wellknown incoherent and partially synchronized states [[2], [3], [4], [5], [6]] corresponding to fixed points of (22, 23), as well as a standing wave state [[10]] corresponding to limitcycle solutions. In addition, we will show that the transitions between these states are mediated by transcritical, saddlenode, Hopf, and homoclinic bifurcations, as well as by three points of higher codimension.
iii.1 Scaling
To ease the notation we begin by scaling Eqs. (22, 23). If we define and nondimensionalize the parameters and time such that
(24)  
we obtain the dimensionless system
(25)  
(26) 
Here the overdot now means differentiation with respect to dimensionless time, and we have dropped all the tildes for convenience. For the rest of this section, all parameters will be assumed to be dimensionless (so there are implicitly tildes over them) unless stated otherwise.
iii.2 Bifurcations of the incoherent state
The incoherent state is defined by , or by in the phase plane formulation. The linearization of the incoherent state, however, is most easily performed in Cartesian coordinates using the formulation in Eqs. (16) and (17). We find the degenerate eigenvalues
(27)  
(28) 
This degeneracy is expected because the origin is always a fixed point and because of the rotational invariance of that state. It follows that the incoherent state is stable if and only if the real parts of the eigenvalues are less than or equal to zero.
The boundary of stable incoherence therefore occurs when the following conditions are met:
These equations define the semicircle and the halfline shown in Fig. 2, labeled TC (for transcritical) and HB (for Hopf bifurcation), respectively. (Independent confirmation of these results can be obtained from the continuous formulation of Eq. (1) directly, as shown in the Appendix.) More precisely, we find that crossing the semicircle corresponds to a degenerate transcritical bifurcation, while crossing the halfline corresponds to a degenerate supercritical Hopf bifurcation.
In the latter case, the associated limitcycle oscillation indicates that the angle increases without bound; this reflects an increasing difference between the phases of the two ‘sub’order parameters of Eqs. (16, 17). In terms of the original model, this means that the oscillator population splits into two counterrotating groups, each consisting of a macroscopic number of oscillators with natural frequencies close to one of the two peaks of . Within each group the oscillators are frequencylocked. Outside the groups the oscillators remain desynchronized, drifting relative to one another and to the locked groups. This is the state Crawford [[10]] called a standing wave. Intuitively speaking, it occurs when the two humps in the frequency distribution are sufficiently far apart relative to their widths. In Kuramoto’s vivid terminology [[3]], the population has spontaneously condensed into “a coupled pair of giant oscillators.”
iii.3 Fixed point solutions and saddlenode bifurcations
Along with the trivial incoherent state , the other fixed points of Eqs. (25, 26) satisfy , and . Using trigonometric identities, we obtain
(29) 
or equivalently,
(30) 
Thus, the fixed point surface is defined implicitly. It can be single or doublevalued as a function of for fixed . To see this, consider how behaves as . We find that
(31) 
from which we observe that the behavior changes qualitatively at , as shown in Fig. 3.
The surface defined by can be plotted parametrically using and , as is seen in Fig. 4. The fold in the surface corresponds to a saddlenode bifurcation. Plots of the phase portrait of reveal that the upper branch of the doublevalued surface in Fig. 3 corresponds to sinks, and the lower branch to saddle points; see Fig. 2 (c), (d), and (g).
In physical terms, the sink represents a stable partially synchronized state, which is familiar from the classic Kuramoto model with a unimodal distribution [[3], [4], [5], [6]]. The oscillators whose natural frequencies are closest to the center of the frequency distribution become rigidly locked, and maintain constant phase relationships among themselves—in this sense, they act collectively like a “single giant oscillator,” as Kuramoto [[3]] put it. Meanwhile the oscillators in the tails of the distribution drift relative to the locked group, which is why one describes the synchronization as being only partial.
The saddle points also represent partially synchronized states, though of course they are unstable. Nevertheless they play an important role in the dynamics because they can annihilate the stable partially synchronized states; this happens in a saddlenode bifurcation along the fold mentioned above. To calculate its location analytically, we use (30) and impose the condition for a turning point, , which yields
(32) 
Eliminating from this equation using (30), we obtain the equation for the saddlenode bifurcation curve
(33) 
This curve is labeled SN in Fig. 2. Its intersection with the semicircle TC occurs at , and is labeled B in the figure. Note also that point C in the figure is not a TakensBagdanov point, as the saddlenode and Hopf bifurcations occur at different locations in the state space; see Figs. 2 (a) and (g).
iii.4 Bistability, homoclinic bifurcations, and SNIPER
An examination of the dynamics corresponding to the approximately triangular parameter space region ABC in Fig. 2 shows bistability. More specifically, we find that the stable incoherent fixed point coexists with the stable partially synchronized state produced by the saddlenode bifurcation described above, as shown in the statespace plot in Fig. 2(c).
Further study of these statespace plots led us to the homoclinic bifurcation curve marked HC, which was obtained numerically. The coexistence of states continues into region ACD, where we found that the stable partially synchronized state now coexists with the stable limit cycle created at the Hopf curve. (See Fig. 2(g).) This limit cycle is then destroyed by crossing the homoclinic curve, which is bounded by point A on one side and by point D on the other.
At point D, the homoclinic curve merges with the saddlenode curve. This codimensiontwo bifurcation, occurring at approximately (1.3589, 0.7483), is known as a saddlenodeloop [[13]]. Below D, however, the saddlenode curve exhibits an interesting feature: the saddlenode bifurcation occurs on an invariant closed curve. This bifurcation scenario is known as a saddlenode infiniteperiod bifurcation, or in short, SNIPER. If we traverse the SNIPER curve from left to right, the sink and saddle (the stable and unstable partially synchronized states) coalesce, creating a loop with infinite period. Beyond that, a stable limit cycle then appears—see Figs. 2 (d), (e), and (f).
In conclusion, we have identified six distinct regions in parameter space and have identified the bifurcations that occur at the boundaries.
Iv Transverse Stability
Our analysis so far has been based on several simplifying assumptions. First, we restricted attention to a special family of oscillator distribution functions and a bimodal Lorentzian form for , which enabled us to reduce the original infinitedimensional system to a threedimensional system of ODEs, Eqs. (1820). Second, we considered only symmetric solutions of these ODEs, by assuming ; this further decreased the dimensionality from three to two.
The next two sections test the validity of these assumptions. We begin here by showing that the nonzero fixed point attractor (the stable partially synchronized state) and the limit cycle attractor (the standing wave state) for Eqns. (25, 26) are transversely stable to small symmetrybreaking perturbations, i.e., perturbations off the invariant manifold defined by . This does not rule out the possible existence of attractors off this manifold, but it does mean that the attractors in the twodimensional symmetric manifold are guaranteed to constitute attractors in the threedimensional ODE system (1820).
Let and consider the reduced governing equations (1820) without symmetry. Introducing the longitudinal and transversal variables
(34) 
and substituting these into (1820), we derive the equation for the transversal component
which describes the order parameter dynamics off the symmetric manifold.
To simplify the notation, let and and scale the system using Eqs. (III.1), as before. Linearization and evaluation at the asymptotic solution denoted by , which may be either a fixed point or a limit cycle, yields the variational equation
(35) 
where
(36) 
Observe that and do not appear in linear order on the right hand side of (35). This decoupling implies that is the eigenvalue associated with the transverse perturbation , in the case where is a fixed point. Similarly, if is a limit cycle, the Floquet exponent associated with is simply , where the brackets denote a time average over one period. Hence the fixed point will be transversely stable if . The analogous condition for the limit cycle is .
iv.1 Fixed point stability
To test the transverse stability of sinks for the twodimensional flow, we solve Eq. (25) for fixed points and obtain
(37) 
Subtracting this from (36), we find
(38) 
Hence is a sufficient condition for transverse stability. But at a nontrivial fixed point,
(39) 
so the transverse stability condition is equivalent to .
We claim that this inequality holds everywhere on the upper branch of the fixed point surface (30). Obviously the inequality is satisfied at all points where . For all other cases, consider the turning point from Fig. 3 defined by . Since the function of interest, , has a global minimum with , and is independent of (at fixed ), it is a lower bound for all on the upper sheet of the fixed point surface, provided that is monotonically decreasing on the interval of . In fact, it is easier to establish that with ; the latter expression is positive, and whenever . Thus transverse stability for the nodes on the fixed point surface follows.
iv.2 Limit cycle stability
To examine the transverse linear stability of the limit cycle, we calculate the transverse Floquet exponent by averaging the eigenvalue over the period of one oscillation:
(40) 
In order to render this expression definite, we rewrite Eq. (25) in terms of the limit cycle solution :
(41) 
Periodicity on the limit cycle guarantees , and so we have
(42) 
which we subtract from the averaged eigenvalue to yield
(43) 
Although we are not able to analytically demonstrate that in (43) is negative, we have calculated and numerically for the limit cycle attractors of Eqs. (1820). This was done for parameter values corresponding to a grid in dimensionless parameter space, by sampling 50 evenly spaced values and . The simulations were run with oscillators. In all the cases that we tested, we found that .
V Numerical Experiments
All of the results described above were obtained using the reduced ODE models derived in Sec. II B and C, and are therefore subject to the restrictions described therein. It is therefore reasonable to ask if these results agree with the dynamics of the original system given in Eq. (1). To check this, a series of direct simulations of Eq. (1) using oscillators and fourthorder RungeKutta numerical integration were performed.
First, we compared solutions of Eq. (1) with those of our reduced system Eqs. (22, 23) in the region where we predicted the coexistence of attractors. For example, we show in Fig. 5 a bifurcation diagram computed along the line that traverses the region ABCD in Fig. 2. (Note that here and for the rest of the paper, we revert to using the original, dimensional form of the variables.) The vertical lines in Fig. 5 indicate the locations of the bifurcations that were identified using the ODE models. For each point plotted, the simulation was run until the order parameter exhibited its timeasymptotic behavior; this was then averaged over the subsequent 5000 time steps. Error bars denote standard deviation. Note in particular the hysteresis, as well as the point with the large error bar, indicating the predicted limit cycle behavior.
Next, we examined the behavior of Eq. (1) at 121 parameter values corresponding to an regular grid superimposed on Fig. 2, ranging from 0.1 to 2.1 at intervals of 0.2 on each axis. (In all cases, was set to 1, and and were varied.) An additional series was run using a smaller grid (from 0.6 to 1.6 at intervals of 0.1 on each axis), to focus on the vicinity of region ABCD in Fig. 2. Initial conditions were chosen systematically in 13 different ways, as follows:

The oscillator phases were uniformly distributed around the circle, so that the overall order parameter had magnitude .

The oscillators were all placed in phase at the same randomly chosen angle in , so that .

The remaining 11 initial conditions were chosen by regarding the system as composed of two subpopulations, one for each Lorentzian in the bimodal distribution of frequencies, as in [[8]]. In one of the subpopulations, the initial phases of the oscillators were chosen to be randomly spaced within the angular sector , where was chosen randomly in and was chosen at random such that the suborder parameter magnitude = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, or 0.9 (all approximately). The result was that had one of these magnitudes and its phase was random in . The same procedure was followed for the other subpopulation, subject to the constraint that . Our idea here was to deliberately break the symmetry of the system initially, to test whether it would be attracted back to the symmetric subspace defined by Eq. (21).
In all the cases we examined, no discrepancies were found between the simulations and the predicted behavior. Although these tests were not exhaustive, and certainly do not constitute a mathematical proof, they are consistent with the conjecture that no additional attractors beyond those described in Section III exist.
We then investigated the generality of our results by replacing the bimodal Lorentzian natural frequency distribution, Eq. (2), with the sum of two Gaussians:
(44) 
and computing the corresponding bifurcation diagram analogous to Fig. 2. The results are shown in Fig. 6. The transcritical (TC) and degenerate Hopf bifurcation (HB) curves were obtained using the continuous formulation of Eq. (1); see the Appendix for details. In addition, saddlenode, homoclinic, and SNIPER bifurcations were numerically observed at several parameter values, and based on these data, we estimated the location of the corresponding curves (dashed lines). All the features of Fig. 2 are preserved, but the curves are somewhat distorted.
Vi Discussion
We conclude by relating our work to three previous studies, and then offer suggestions for further research, both theoretical and experimental.
vi.1 Kuramoto’s conjectures
In his book on coupled oscillators, Kuramoto [[3]] speculated about how the transition from incoherence to mutual synchronization might be modified if the oscillators’ natural frequencies were bimodally distributed across the population. On pp.75–76 of Ref. [[3]], he wrote “So far, the nucleation has been supposed to be initiated at the center of symmetry of . This does not seem to be true, however, when is concave there.” His reasoning was that for a bimodal system, synchrony would be more likely to start at the peaks of . If that were true, it would mean that a system with two equal peaks would go directly from incoherence to having two synchronized clusters of oscillators, or what we have called the standing wave state, as the coupling is increased. The critical coupling at which this transition would occur, he argued, should be , analogous to his earlier result for the unimodal case. According to this scenario, the synchronized clusters would be tiny at onset, comprised only of oscillators with natural frequencies near the peaks of . Because of their small size, Kuramoto claimed these clusters “will behave almost independently of each other.” With further increases in , however, the clusters “will come to behave like a coupled pair of giant oscillators, and for even stronger coupling they will eventually be entrained to each other to form a single giant oscillator.” (This is what we have called the partially synchronized state.)
Let us now reexamine Kuramoto’s conjectures in light of our analytical and numerical results, as summarized in Fig. 7(a). For a fair comparison, we must assume that is concave at its center frequency ; for the bimodal Lorentzian (Eq. (2), this is equivalent to . (Otherwise is unimodal and incoherence bifurcates to partial synchronization as is increased, consistent with Kuramoto’s classic result as well as the lowest portion of Fig. 7(a).)
So restricting attention from now on to the upper part of Fig.7(a) where , what actually happens as increases? Was Kuramoto right that the bifurcation sequence is always incoherence standing wave partial synchronization?
No. For between and (meaning the distribution is just barely bimodal), incoherence bifurcates directly to partial synchronization—the “single giant oscillator” state—without ever passing through an intermediate standing wave state. In effect, the system still behaves as if it were unimodal. But there is one new wrinkle: we now see hysteresis in the transition between incoherence and partial synchronization, as reflected by the lower bistable region in Fig. 7(a).
Is there any part of Fig. 7(a) where Kuramoto’s scenario really does occur? Yes—but it requires that the peaks of be sufficiently well separated. Specifically, suppose , the value at the codimension2 saddlenodeloop point where the homoclinic and SNIPER curves meet (i.e., point D in Fig. 2). In this regime everything behaves as Kuramoto predicted.
An additional subtlety occurs in the intermediate regime where the peaks of are neither too far apart nor too close together. Suppose that . Here the system shows a different form of hysteresis. The bifurcations occur in the sequence that Kuramoto guessed as increases, but not on the return path. Instead, the system skips the standing wave state and dissolves directly from partial synchronization to incoherence as is decreased.
Finally we note that Kuramoto’s conjectured formula is incorrect, although it becomes asymptotically valid in the limit of widely separated peaks. Specifically, his prediction is equivalent to , which approaches the correct result as .
vi.2 Crawford’s center manifold analysis
Crawford [[10]] obtained the first mathematical results for the system studied in this paper. Using center manifold theory, he calculated the weakly nonlinear behavior of the infinitedimensional system in the neighborhood of the incoherent state. From this he derived the stability boundary of incoherence. His analysis also included the effects of white noise in the governing equations.
Figure 7(b), reproduced from Fig. 4 in Ref. [[10]], summarizes Crawford’s findings. Here is the noise strength (note: our analysis is limited to ), is the width of the Lorentzians (equivalent to in our notation), and are the center frequencies of the Lorentzians (as here). The dashed line in Fig. 7(b) shows Crawford’s schematic depiction of the unknown stability boundary between the standing waves and the partially synchronized state. He suggested a strategy for calculating this boundary, and highlighted it as an open problem, writing in the figure caption, “…the precise nature and location of this boundary have not been determined.” Our results, summarized in Figs. 2 and Fig. 7(a), now fill in the parts that were missing from Crawford’s analysis.
vi.3 Stochastic model of Bonilla et al.
In a series of papers (see [[6]] for a review), Bonilla and his colleagues have explored what happens if one replaces the Lorentzians in the frequency distribution with functions, and adds white noise to the governing equations. The resulting system can be viewed as a stochastic counterpart of the model studied here; in effect, the noise blurs the functions into bellshaped distributions analogous to Lorentzians or Guassians. And indeed, the system shows much of the same phenomenology as seen here: incoherence, partially synchronized states, standing waves, and bistability [[6]].
However, a complete bifurcation diagram analogous to Fig. 2 has not yet been worked out for this model. The difficulty is that no counterpart of the ansatz (11) has been found; the stochastic problem is governed by a secondorder FokkerPlanck equation, not a firstorder continuity equation, and the OttAntonsen ansatz (11) no longer works in this case. Perhaps there is some way to generalize the ansatz appropriately so as to reduce the stochastic model to a lowdimensional system, but for now this remains an open problem.
vi.4 Directions for future research
There are several other questions suggested by the work described here.
Validity of reduction method
The most important open problem is to clarify the scope and limits of the OttAntonsen method used in Sec. II.2. Under what conditions is it valid to assume that the infinitedimensional Kuramoto model can be replaced by the lowdimensional dynamical system implied by the OttAntonsen ansatz? Or to ask it another way, when do all the attractors of the infinitedimensional system lie in the lowdimensional invariant manifold corresponding to this ansatz?
This question has now become particularly pressing, because two counterexamples have recently come to light in which the OttAntonsen method [[11]] gives an incomplete account of the full system’s dynamics. When the method was applied to the problem of chimera states for two interacting populations of identical phase oscillators, it predicted only stationary and periodic chimeras [[14]], whereas subsequent numerical experiments revealed that quasiperiodic chimeras can also exist and be stable [[15]]. Likewise, chaotic states are known to emerge from a wide class of initial conditions for series arrays of identical overdamped Josephson junctions coupled through a resistive load [[16], [17]]. Yet the OttAntonsen ansatz cannot account for these chaotic states, because the reduced ODE system turns out to be only two dimensional [[18], [19]].
What makes this all the more puzzling is that the method works so well in other cases. It seems to give a full inventory of the attractors for the bimodal Kuramoto model studied here, as well as for the unimodal Kuramoto model in its original form [[2], [3], [11]] or with external periodic forcing [[11], [12]].
So we are left in the unsatisfying position of not knowing when the method works, or why. In some cases it (apparently) captures all the attractors, while in other cases it does not. How does one make sense of all this?
A possible clue is that in all the cases where the method has so far been successful, the individual oscillators were chosen to have randomly distributed frequencies; whereas in the cases where it failed, the oscillators were identical. Perhaps the mixing induced by frequency dispersion is somehow relevant here?
A resolution of these issues may come from a new analytical approach. Pikovsky and Rosenblum [[15]] and Mirollo, Marvel and Strogatz [[18]] have independently shown how to place the OttAntonsen ansatz [[11]] in a more general mathematical framework by relating it to the group of Mobius transformations [[18], [20]] or, equivalently, to a trigonometric transformation [[15]] originally introduced in the study of Josephson arrays [[17]]. This approach includes the OttAntonsen ansatz as a special case, but is more powerful in the sense that it provably captures all the dynamics of the full system, and it works for any , not just in the infinite limit. The drawback is that the analysis becomes more complicated. It remains to be seen what conclusions can be drawn—and, perhaps, what longstanding problems can be solved—when this new approach is unleashed on the Kuramoto model and its many relatives.
Even in those instances where the OttAntonsen ansatz doesn’t account for all the attractors of the full system, it can still provide useful information, for instance by giving at least some of the attractors and by easing the calculation of them. Moreover, the transient evolution from initial conditions off the OttAntonsen invariant manifold can yield interesting phenomena not captured by the ansatz, as discussed in Appendix C of [[21]].
Asymmetric bimodal distributions
Now returning to the specific problem of the bimodal Kuramoto model: What happens if the humps in the bimodal distribution have unequal weights? The analysis could proceed as in this paper, up to the point where we assumed symmetry between the two subpopulations. One would expect new phenomena such as traveling waves to arise because of the broken symmetry.
Finitesize effects
Comparison with experiment
Vii Acknowledgments
This research was supported in part by NSF grant DMS0412757 and ONR award N00014070734.
Appendix A Alternative calculation of the boundary of stability for the incoherent state
The system in Eq. (1), together with the bimodal natural frequency distribution given in Eq. (2), can be expressed using the formulation in [[8]] as two interacting populations of oscillators. In this case, each population has a separate Lorenzian frequency distribution of width and center frequency at or , and the twobytwo matrix describing the relative coupling weights (i.e, Eq. (1) in [[8]]) has in each entry. By postulating that a small perturbation to the incoherent state grows exponentially as , and setting for the marginally stable state, Eq. (9) of Ref. [[8]] gives the following expression for the critical coupling value :
(45) 
The boundary of stability of the incoherent state is obtained by requiring that this expression be strictly real. One solution is obtained for , resulting in , which is equivalent to
(46) 
This is the equation for the semicircle in Figure 2, corresponding to a transcritical bifurcation of the incoherent state. Another solution, obtained by assuming that in Eq. (45) and requiring , is . This is the equation for the halfline in Figure 2 corresponding to the degenerate Hopf bifurcation of the incoherent state.
If the bimodal natural frequency distribution is given by a sum of Gaussians of standard deviation and centers at , then the twopopulation approach outlined above leads to the following equation:
(47) 
where
(48) 
is known as the Faddeeva function and can be computed numerically [[26]]. Once again requiring that be real, two branches corresponding to being equal and not equal to zero can be obtained. These are the boundaries of stability of the incoherent state shown in Fig. 6.
References
 Some examples are the following: the coupling of oscillatory neurons in the suprachiasmatic nucleus of the brain governing circadian rhythms [S. Yamaguchi, et al., Science 302, 1408 (2002)]; the interaction of cells containing oscillatory chemically reacting constituents [I. Z. Kiss, Y. Zhai and J. L. Hudson, Science 296, 1676 (2002)]; Josephson junction circuits [K. Wiesenfeld and J. W. Swift, Phys. Rev. E 51, 1020 (1995)]; and neutrino oscillations [J. Pantaleone, Phys. Rev. D 58, 073002 (1998)].
 Y. Kuramoto, in International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics, Vol. 39, edited by H. Araki (SpringerVerlag, Berlin, 1975).
 Y. Kuramoto, Chemical Oscillations, Waves and Turbulence (Springer, 1984), pp. 75–76.
 S. H. Strogatz, Physica D 143, 1 (2000).
 E. Ott, Chaos in Dynamical Systems, second edition, chapter 6 (Cambridge University Press, 2002).
 J. A. Acebron, L. L. Bonilla, C. J. Perez Vicente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005).
 S. H. Strogatz, R. E. Mirollo and P. C. Matthews, Phys. Rev. Lett. 68, 2730 (1992).
 E. Barreto, B. Hunt, E. Ott and P. So, Phys. Rev. E 77, 036107 (2008).
 E. Montbrio, J. Kürths and B. Blasius, Phys. Rev. E 70, 056125 (2004).
 J. D. Crawford, J. Stat. Phys. 74, 1047 (1994).
 E. Ott and T. M. Antonsen, Chaos 18, 037113 (2008)
 T. M. Antonsen, R. Faghih, M. Girvan, E. Ott and J. Platig, Chaos 18, 0037112 (2008); L. M. Childs and S. H. Strogatz, arXiv:0807.4717 (2008).

J. Guckenheimer, Physica D 20, 1 (1986); J. Guckenheimer & P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations, (Springer, New York, 1983); J. Tyson,
http://mpf.biol.vt.edu/research/generic_model/ main/pp/intro.php
 D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. Wiley, Phys. Rev. Lett. 101, 084103 (2008).
 A. Pikovsky and M. Rosenblum, arXiv:0809.3700 (2008).
 D. Golomb, D. Hansel, B. Shraiman, and H. Sompolinsky, Phys. Rev. A 45, 3516 (1992)
 S. Watanabe and S. H. Strogatz, Physica D 74, 197 (1994).
 R. E. Mirollo, S. A. Marvel, and S. H. Strogatz (in preparation).
 S. A. Marvel and S. H. Strogatz (in preparation).
 C. J. Goebel, Physica D 80, 18 (1995).
 E. Ott, J. Platig, T. M. Antonsen, and M. Girvan, Chaos 18, 037115 (2008).
 E. J. Hildebrand, M. A. Buice, and C. C. Chow, Phys. Rev. Lett. 98, 054101 (2007).
 M. A. Buice, and C. C. Chow, Phys. Rev. E 76, 031118 (2007).
 I. Z. Kiss, W. Wang, and J. L. Hudson, Chaos 12, 252 (2002).
 A. S. Mikhailov, D. H. Zanette, Y. M. Zhai, I. Z. Kiss, and J. L. Hudson, Proc. Natl. Acad. Sci USA 101, 10890 (2004).
 J.A.C. Weideman, SIAM Journal on Numerical Analysis, 31, 14971518 (1994).