CONTENTS

Synchronization phenomena in large populations of interacting elements are the subject of intense research efforts in physical, biological, chemical, and social systems. A successful approach to the problem of synchronization consists of modeling each member of the population as a phase oscillator. In this review, synchronization is analyzed in one of the most representative models of coupled phase oscillators, the Kuramoto model. A rigorous mathematical treatment, specific numerical methods, and many variations and extensions of the original model that have appeared in the last few years are presented. Relevant applications of the model in different contexts are also included.


I. INTRODUCTION
Time plays a key role for all living beings. Their activity is governed by cycles of different duration which determine their individual and social behavior. Some of these cycles are crucial for their survival. There are biological processes and specific actions which require precise timing. Some of these actions demand a level of expertise that can only be acquired after a long period of training, but others take place spontaneously. How do these actions occur? Possibly through synchronization of individual actions in a population. A few examples follow. Suppose we attend a concert. Each member of the orchestra plays a sequence of notes that, properly combined according to a musical composition, elicit a deep feeling in us as listeners. The effect can be astonishing or a fiasco ͑apart from other technical details͒ simply depending on the exact moment when the sound was emitted. In the meantime, our hearts are beating rhythmically because thousands of cells synchronize their activity. The emotional character of the music can accelerate or decelerate our heartbeats. We are not aware of the process, but the cells themselves manage to change coherently, almost in unison. How? We see the conductor rhythmically moving his arms. Musicians know exactly how to interpret these movements and respond with the appropriate action. Thousands of neurons in the visual cortex, sensitive to specific space orientations, synchronize their activity almost immediately while the baton describes a trajectory in space. This information is transmitted and processed through some remarkably fast mechanisms. What more? Just a few seconds after the last bar, the crowd occupying the auditorium starts to applaud. At the beginning the rhythm may be incoherent, but the wish to get an encore can transform incoherent applause into perfectly synchronized applause, despite the different strength in beating or the location of individuals inside the concert hall.
These examples illustrate synchronization, one of the most captivating cooperative phenomena in nature. Synchronization is observed in biological, chemical, physical, and social systems, and it has attracted the interest of scientists for centuries. A paradigmatic example is the synchronous flashing of fireflies observed in some South Asian forests. At night, myriad fireflies rest in the trees. Suddenly, several fireflies start emitting flashes of light. Initially they flash incoherently, but after a short period of time the whole swarm is flashing in unison, creating one of the most striking visual effects ever seen. The relevance of synchronization has been stressed frequently although it has not always been fully understood. In the case of the fireflies, synchronous flashing may facilitate the courtship between males and females. In other cases, the biological role of synchronization is still under discussion. Thus perfect synchronization could lead to disaster and extinction, and therefore different species in the same trophic chain may develop different circadian rhythms to enlarge their probability of survival. Details about these and many other systems, together with many references, can be found in the recent excellent book by Strogatz ͑2003͒.
Research on synchronization phenomena focuses inevitably on ascertaining the main mechanisms responsible for collective synchronous behavior among members of a given population. To attain a global coherent activity, interacting oscillatory elements are required. The rhythmical activity of each element may be due to internal processes or to external sources ͑external stimuli or forcing͒. Even if the internal processes responsible for rhythmicity have different physical or biochemical origins and can be very complex, one may hope to understand the essence of synchronization in terms of a few basic principles. What might these principles be?
There are different ways to tackle this problem. Suppose that the rhythmical activity of each element is described in terms of a physical variable that evolves regularly in time. When such a variable reaches a certain threshold, the element emits a pulse ͑action potential for neurons͒, which is transmitted to the neighborhood. Later on, a resetting mechanism initializes the state of this element. Then, a new cycle starts. Essentially the behavior of each element is similar to that of an oscillator. Assuming that the rhythm has a certain period, it is convenient to introduce the concept of phase, a periodic measure of the elapsed time. The effect of the emitted pulse is to alter the current state of the neighbors by modifying their periods, lengthening or shortening them. This disturbance depends on the current state of the oscillator receiving the external impulse, and it can also be studied in terms of a phase shift. The analysis of the collective behavior of the system can be carried out in this way under two conditions: ͑i͒ the phase shift induced by an impulse is independent of the number of impulses arriving within an interspike interval, and ͑ii͒ the arrival of one impulse affects the period of the current time interval, but memory thereof is rapidly lost and the behavior in future intervals is not affected.
There is another scenario in which synchronization effects have been studied extensively. Let us consider an ensemble of nonlinear oscillators moving in a globally attracting limit cycle of constant amplitude. These are phase-or limit-cycle oscillators. We now couple them weakly to ensure that no disturbance will take any of them away from the global limit cycle. Therefore only one degree of freedom is necessary to describe the dynamic evolution of the system. Even at this simple level of description it is not easy to propose specific models. The first scenario of pulse-coupled oscillators is perhaps more intuitive, more direct, and easier to model. However, the discrete and nonlinear nature of pulse coupling gives rise to important mathematical complications. While the treatment of just a few pulse-coupled elements can be done within the framework of dynamical systems, the description becomes much more complicated for a large number of such elements. Proposing a model within the second scenario of coupled limit-cycle oscillators leaves ample room for imagination. We are forced to consider models that are mathematically tractable, with continuous time and specific nonlinear interactions between oscillators. Our experience tells us that models with the latter property are exceptional. Nevertheless some authors have been looking for a "solvable" model of this type for years. Winfree, for example, persistently sought a model with nonlinear interactions ͑Winfree, 1967, 1980͒. He realized that synchronization can be understood as a threshold process. When the coupling among oscillators is strong enough, a macroscopic fraction of them synchronizes to a common frequency. The model he proposed was hard to solve in its full generality, although a solvable version has been recently found ͑Ariaratnam and . Hence research on synchronization proceeded along other directions.
The most successful attempt was due to Kuramoto ͑1975͒, who analyzed a model of phase oscillators running at arbitrary intrinsic frequencies and coupled through the sine of their phase differences. The Kuramoto model is simple enough to be mathematically tractable, yet sufficiently complex to be nontrivial. The model is rich enough to display a large variety of synchronization patterns and sufficiently flexible to be adapted to many different contexts. This "little wonder" is the object of this review. We have reviewed the progress made in the analysis of the model and its extensions during the last 28 years. We have also tried to cover the most significant areas where the model has been applied, although we realize that this is not an easy task because of its ubiquity.
The review is organized as follows. The Kuramoto model with mean-field coupling is presented in Sec. II. In the limit of infinitely many oscillators, we discuss the characterization of incoherent, phase-locked, and partially synchronized phases. The stability of the partially synchronized state, finite-size effects, and open problems are also considered. Section III concerns the noisy mean-field Kuramoto model, resulting from adding external white-noise sources to the original model. This section deals with the nonlinear Fokker-Planck equation describing the one-oscillator probability density in the limit of infinitely many oscillators ͑which is derived in Appendix A͒. We study synchronization by first analyzing the linear stability of the simple nonsynchronized state called incoherence, in which every phase of the oscillators is equally probable. Depending on the distribution of natural frequencies, different synchronization scenarios can occur in parameter regions where incoherence is unstable. We present a complete analysis of these scenarios for a bimodal frequency distribution using bifurcation theory.
Our original presentation of bifurcation calculations exploits the Chapman-Enskog method to construct the bifurcating solutions, which is an alternative to the method of multiple scales for degenerate bifurcations and is simpler than using center-manifold techniques. Section IV describes the known results for the Kuramoto model with couplings that are not of the meanfield type. They include short-range and hierarchical couplings, models with disorder, time-delayed couplings, and models containing external fields or multiplicative noise. Extensions of the original model are discussed in Sec. V. Section VI discusses numerical solutions of the noisy Kuramoto model, both for the system of stochastic differential equations and for the nonlinear Fokker-Planck equation describing the one-oscillator probability density in the limit of infinitely many oscillators. Applications of the Kuramoto model are considered in Sec. VII. They include neural networks, Josephson junctions and laser arrays, and chemical oscillators. These applications are often directly inspired by the original model, share its philosophy, and represent an additional step toward the development of new ideas. The last section contains our conclusions and discusses some open problems and hints for future work. Some technical details are collected in five appendixes.

II. THE KURAMOTO MODEL
The Kuramoto model consists of a population of N coupled phase oscillators i ͑t͒ having natural frequencies i distributed with a given probability density g͑͒, and whose dynamics are governed by Thus each oscillator tries to run independently at its own frequency, while the coupling tends to synchronize it to all the others. By making a suitable choice of rotating frame, i → i − ⍀t, in which ⍀ is the first moment of g͑͒, we can transform Eq. ͑1͒ to an equivalent system of phase oscillators whose natural frequencies have zero mean. When the coupling is sufficiently weak, the oscillators run incoherently, whereas beyond a certain threshold collective synchronization emerges spontaneously. Many different models for the coupling matrix K ij have been considered such as nearest-neighbor coupling, hierarchical coupling, random long-range coupling, or even state-dependent interactions. All of them will be discussed in this review.
In this section, we introduce the Kuramoto model with mean-field coupling among phase oscillators. For this model, synchronization is conveniently measured by an order parameter. In the limit of infinitely many oscillators, N = ϱ, the amplitude of the order parameter vanishes when the oscillators are out of synchrony, and it is positive in synchronized states. We first present Kuramoto's calculations for partial synchronization of oscillators and bifurcation from incoherence, a state in which the oscillator phase takes values on the interval ͓− , ͔ with equal probability. The stability of incoherence is then analyzed in the limit N = ϱ. For coupling constant K Ͻ K c , a critical value of the coupling, incoherence is neutrally stable because the spectrum of the operator governing its linear stability lies on the imaginary axis. This means that disturbances from incoherence decay similarly to the Landau damping in plasmas ͑Strogatz et al., 1992͒. When K Ͼ K c and unimodal natural frequency distributions are considered, one positive eigenvalue emerges from the spectrum. The partially synchronized state bifurcates from incoherence at K = K c , but a rigorous proof of its stability is still missing. Finally, finite-size effects ͑N Ͻϱ͒ on oscillator synchronization are discussed.

A. Stationary synchronization for mean-field coupling
The original analysis of synchronization was accomplished by Kuramoto in the case of mean-field coupling, that is, taking K ij = K / N Ͼ 0 in Eq. ͑1͒ ͑Kuramoto, 1975. The model of Eq. ͑1͒ was then written in a more convenient form, defining the ͑complex-valued͒ order parameter Here r͑t͒ with 0 ഛ r͑t͒ ഛ 1 measures the coherence of the oscillator population, and ͑t͒ is the average phase. With this definition, Eq. ͑1͒ becomes i = i + Kr sin͑ − i ͒, i = 1,2, ... ,N, ͑3͒ and it is clear that each oscillator is coupled to the common average phase ͑t͒ with coupling strength given by Kr. The order parameter ͑2͒ can be rewritten as

͑4͒
In the limit of infinitely many oscillators, they may be expected to be distributed with a probability density ͑ , , t͒, so that the arithmetic mean in Eq. ͑2͒ now becomes an average over phase and frequency, namely, This equation illustrates the use of the order parameter to measure oscillator synchronization. In fact, when K → 0, Eq. ͑3͒ yields i Ϸ i t + i ͑0͒, that is, the oscillators rotate at angular frequencies given by their own natural frequencies. Consequently, setting Ϸ t in Eq. ͑5͒, by the Riemann-Lebesgue lemma, we obtain that r → 0 as t → ϱ and the oscillators are not synchronized. On the other hand, in the limit of strong coupling, K → ϱ, the oscillators become synchronized to their average phase, i Ϸ , and Eq. ͑5͒ implies r → 1. For intermediate couplings, K c Ͻ K Ͻϱ, part of the oscillators are phase locked ͑ i =0͒, and part are rotating out of synchrony with the locked oscillators. This state of partial synchronization yields 0 Ͻ r Ͻ 1 and will be further explained below. Thus synchronization in the mean-field Kuramoto model ͑with N = ϱ͒ is revealed by a nonzero value of the order parameter. The concept of order parameter as a measure of synchronization is less useful for models with short-range coupling. In these systems, other concepts are more appropriate to describe oscillator synchronization since more complex situations can happen ͑Strogatz and Mirollo, 1988a. For instance, it could happen that a finite fraction of the oscillators have the same average frequency i , defined by while the other oscillators might be out of synchrony, or that the phases of a fraction of the oscillators could change at the same speed ͑and therefore partial synchronization would occur͒, while different oscillator groups had different speeds ͑and therefore their global order parameter could be zero and incoherence would result͒. See Sec. IV for details. A continuity equation for the oscillator density can be found by noting that each oscillator in Eq. ͑1͒ moves with an angular or drift velocity v i = i + Kr sin͑ − i ͒. Therefore the one-oscillator density obeys the continuity equation to be solved along with Eq. ͑5͒, with the normalization condition ͵ − ͑,,t͒d = 1, ͑8͒ and an appropriate initial condition. The system of Eqs. ͑5͒-͑8͒ has the trivial stationary solution =1/͑2͒, r =0, corresponding to an angular distribution of oscillators having equal probability in the interval ͓− , ͔. Then, the oscillators run incoherently, and hence the trivial solution is called the incoherent solution, or simply incoherence. Let us now try to find a simple solution corresponding to oscillator synchronization. In the strongcoupling limit, we have global synchronization ͑phase locking͒, so that all oscillators have the same phase, i = ͓= i t + i ͑0͔͒, which yields r = 1. For a finite coupling, a lesser degree of synchronization with a stationary amplitude, 0 Ͻ r Ͻ 1, may occur. How can this smaller value of r arise? A typical oscillator running with velocity v = − Kr sin͑ − ͒ will become stably locked at an angle such that Kr sin͑ − ͒ = and − /2ഛ ͑ − ͒ ഛ / 2. All such oscillators are locked in the natural laboratory frame of reference. Oscillators with frequencies ͉͉ Ͼ Kr cannot be locked. They run out of synchrony with the locked oscillators, and their stationary density obeys v = C ͑constant͒, according to Eq. ͑7͒. We have obtained a stationary state of partial synchronization, in which part of the oscillators are locked at a fixed phase while all others are rotating out of synchrony with them. The corresponding stationary density is therefore Here H͑x͒ =1 if x Ͼ 0 and H͑x͒ = 0 otherwise, that is, H͑x͒ is the Heaviside unit step function. Note that one can write equivalently = ͱ K 2 r 2 − 2 ␦͓ − Kr sin͑ − ͔͒H͑cos ͒ for −Kr Ͻ Ͻ Kr. The normalization condition ͑8͒ for each frequency yields C = ͱ 2 − ͑Kr͒ 2 / ͑2͒.
We can now evaluate the order parameter in the state of partial synchronization by using Eqs. ͑5͒ and ͑9͒, Let us assume that g͑͒ = g͑−͒. Then, the symmetry relation ͑ + ,−͒ = ͑ , ͒ implies that the second term in this equation is zero. The first term is simply cos g͑Kr sin ͒Kr cos d, that is, This equation always has the trivial solution r = 0 corresponding to incoherence, = ͑2͒ −1 . However, it also has a second branch of solutions corresponding to the partially synchronized phase ͑9͒, satisfying This branch bifurcates continuously from r = 0 at the value K = K c obtained by setting r = 0 in Eq. ͑12͒, which yields K c =2/͓g͑0͔͒. Such a formula and the argument leading to it were first found by Kuramoto ͑1975͒. Considering, as an example, the Lorentzian frequency distribution g͑͒ = ␥/ ␥ 2 + 2 , ͑13͒ allows an explicit evaluation of the integrals above to be accomplished, which was done by Kuramoto ͑1975͒. Using Eq. ͑13͒, he found the exact result r = ͱ1−͑K c / K͒ for all K Ͼ K c =2␥. For a general frequency distribution g͑͒, an expansion of the right-hand side of Eq. ͑11͒ in powers of Kr yields the scaling law as K → K c . Throughout this review, we use the following definitions of the symbol ϳ ͑asymptotic͒ which compares two functions or one function and an asymptotic series in the limit as ⑀ → 0 ͑Bender and Orszag, 1978͒: According to Eq. ͑14͒, the partially synchronized phase bifurcates supercritically for K Ͼ K c if gЉ͑0͒ Ͻ 0, and subcritically for K Ͻ K c if gЉ͑0͒ Ͼ 0; see Figs. 1͑a͒ and 1͑b͒. Notice that Kuramoto's calculation of the partially synchronized phase does not indicate whether this phase is stable, either globally or locally.

Synchronization in the limit N = ؕ
Kuramoto's original construction of incoherent and partially synchronized phases concerns purely stationary states. Moreover, he did not establish any of their stability properties. The linear stability theory of incoherence was published by Strogatz et al. ͑1992͒ and interesting work on the unsolved problems of nonlinear stability theory was carried out by Balmforth and Sassi ͑2000͒. To ascertain the stability properties of the incoherent and partially synchronized solutions, it is better to work with the probability density ͑ , , t͒. Let us explain first what is known in the limit of infinitely many oscillators described by Eqs. ͑5͒-͑7͒. The linearized stability problem for this case is obtained by inserting =1/͑2͒ + ͑ , t ; ͒ with ͑ , t ; ͒ = exp͑t͒͑ , ͒ in Eqs. ͑5͒-͑8͒, and then ignoring terms nonlinear in : If Re Ͻ0 for all possible , incoherence is linearly stable, while it is unstable if some admissible has a positive real part. The periodicity condition implies = ͚ n=−ϱ ϱ b n ͑͒e in , which, inserted in Eq. ͑17͒, yields

͑19͒
Here we have used b −n = b n ͑a bar over a symbol denotes taking its complex conjugate͒, and defined the scalar product shows that n ͑͒ =−in, n = ±1, ±2,..., belong to the continuous spectrum describing the linear stability problem, provided is in the support of g͑͒.
͑Balmforth and Sassi, 2000͒. The function contains a term proportional to e −it , which is nondecaying and nonseparable, and does not correspond to a normal mode. As time elapses, this term becomes progressively more crenellated, and through increasing cancellations, integral averages of decay. Besides this, Eq. ͑21͒ contain two exponentially decaying terms which contribute to the order parameter ͑22͒. If g͑͒ has bounded support, the order parameter may decay more slowly, • Regularizing the Kuramoto model by adding a diffusive term to Eq. ͑7͒ and constructing the stationary probability density in the limit of vanishing diffusivity by means of boundary layer methods. This regularization corresponds to adding white-noise forcing terms to the Kuramoto model ͑1͒. The corresponding equation for the probability density is Eq. ͑7͒ with a diffusive term in its right-hand side, which is called the nonlinear Fokker-Planck equation. The nonlinear Fokker-Planck equation will be studied in Sec. III.
• A mixture of multiple scales and boundary layer ideas in the same limit of vanishing diffusivity of the nonlinear Fokker-Planck equation, but for K near values corresponding to the bifurcation of synchronized states from incoherence.
In the small-diffusivity limit of the nonlinear Fokker-Planck equation, D → 0 + , the calculations by Balmforth and Sassi ͑2000͒ indicate that the probability density with unimodal frequency distribution and K Ͼ K c tends toward a stationary phase that is concentrated around the curve = Kr sin for 2 Ͻ K 2 r 2 . The peak of the probability density is attained at ͱKr/͑2D͒. It would be useful to have consistent perturbation results for the evolution of the probability density near bifurcation points, in the low-noise limit D → 0 + of the nonlinear Fokker-Planck equation, and also for the Kuramoto model with D = 0. Both cases are clearly different, as shown by the fact that the synchronized phase is a generalized function ͑a distribution͒ when D = 0, while it is a smooth function for D Ͼ 0. In particular, it is clear that Kuramoto's partial synchronization solution ͑9͒ ͑involving a delta function͒ cannot be obtained by smallamplitude bifurcation calculations about incoherence, the laborious attempt by Crawford and Davies ͑1999͒ notwithstanding.
Thus understanding synchronization in the hyperbolic limit of the mean-field Kuramoto model ͑with N → ϱ͒ requires the following. First, a fully consistent asymptotic description of the synchronized phase and the synchronization transition as D → 0+ should be found. As indicated by Balmforth and Sassi ͑2000͒, the necessary technical work involves boundary layers and matching. These calculations could be easier if one works directly with the equations for ln , as suggested in the early paper ͑Bonilla, 1987͒. Second, and most likely harder, the same problems should be tackled for D = 0, where the stable synchronized phase is expected to be Kuramoto's partially synchronized state ͑which is a distribution, not a smooth function͒. Third, as pointed out by Strogatz ͑2000͒, the problem of proving stability of the partially synchronized state as a solution of the Kuramoto model remains open.

Finite-size effects
Another way to regularize the hyperbolic equation ͑7͒ is to study a large population of finitely many phase oscillators, which are globally coupled. The analysis of this large system may shed some light on the stability properties of the partially synchronized state. The question can be posed as follows. What is the influence of finitesize effects on Kuramoto's partially synchronized state as N → ϱ?
One issue with kinetic equations describing populations of infinitely many elements is always that of finitesize effects. This issue was already raised by Zermelo's paradox, namely, that a system of finitely many particles governed by reversible classical Hamiltonian mechanics was bound to have recurrences according to Poincaŕe's recurrence theorem. Then, this system would come back arbitrarily close to its initial condition infinitely many times. Boltzmann's answer to this paradox was that the recurrence times would become infinite as the number of particles tend to infinite. Simple model calculations illustrate the following fact. A nonrecurrent kinetic description for a system of infinitely many particles approximates the behavior of a system with a large but finite number of particles during finite time intervals, after which recurrences set in ͑Keller and Bonilla, 1986͒. The same behavior, denoting the noncommutativity of the limits N → ϱ and t → ϱ, is also present in the Kuramoto model. For instance, Hemmen and Wreszinski ͑1993͒ used a Lyapunov-function argument to point out that a population of finitely many Kuramoto oscillators would reach a stationary state as t → ϱ. Our derivation of the nonlinear Fokker-Planck equation in Appendix A suggests that fluctuations scale as N −1/2 as N → ϱ, a scaling that, for the order parameter, is confirmed by numerical simulations ͑Kuramoto and Nishikawa, 1987;Daido, 1990͒. More precise theoretical results are given by DaiPra and den Hollander ͑1996͒, for rather general mean-field models that include Kuramoto's and also spin systems. DaiPra and den Hollander ͑1996͒ obtained a central limit theorem, which characterizes fluctuations about the one-oscillator probability density, for N = ϱ, as Gaussian fields having a certain covariance matrix and scaling as N −1/2 . Near bifurcation points, a different scaling is to be expected, similarly to Dawson's results for related meanfield models ͑Dawson, 1983͒. Daido ͑1987b, 1989͒ explored this issue by dividing the oscillator phase and the order parameter in Eq. ͑2͒ into two parts: their limits achieved when N → ϱ and their fluctuating parts ͑which were regarded as small͒. In the equations for the phase fluctuations, only terms linear in the fluctuation of the order parameter were retained. The result was then inserted into Eq. ͑2͒, and a self-consistent equation for the fluctuation of the order parameter was found. For unimodal frequency distributions, Daido found the scaling ͓͑K c − K͒N͔ −1/2 for the rms fluctuation of the order parameter as K → K c − ͑from below͒. For coupling strengths larger than K c , he found that the fluctuation of the order parameter was consistent with the scaling ͑K − K c ͒ −1/8 N −1/2 as K → K c + ͑Daido, 1989͒. Balmforth and Sassi ͑2000͒ carried out numerical simulations to investigate finite-size effects, and discussed how sampling the natural frequency distribution affects the one-oscillator probability density. In particular, they found that sampling may give rise to unexpected effects, such as timeperiodic synchronization, even for populations with unimodal frequency distribution, g͑͒. Note that such effects do not appear in the limit N = ϱ. Further work on this subject would be interesting, in particular, finding a formulation similar to Daido's fluctuation theory for the order parameter but for the one-oscillator probability density instead.

III. THE MEAN-FIELD MODEL INCLUDING WHITE-NOISE FORCES
In this section, we analyze the mean-field Kuramoto model with white-noise forcing terms. This generalization makes the model more physical, in that unavoidable random imperfections can be taken into account. At the same time, the ensuing model turns out to be mathematically more tractable. In fact, in the limit of infinitely many oscillators, the one-oscillator probability density obeys a parabolic equation ͑the nonlinear Fokker-Planck equation͒, instead of the hyperbolic equation ͑7͒, which is harder to analyze. The nonlinear Fokker-Planck equation is derived in Appendix A. Its simplest solution is also = ͑2͒ −1 , corresponding to incoherence. First, we shall study its linear stability properties. When K Ͻ K c the incoherence turns out to be linearly stable, unlike what happens in the Kuramoto model, for which incoherence is neutrally stable. When the coupling strength is greater than the critical coupling, incoherence becomes linearly unstable, since the real part of one of the eigenvalues in the discrete spectrum of the linearized problem becomes positive. Then, different bifurcation scenarios and phase diagrams will occur, depending on the distribution of natural frequencies g͑͒. Synchronized phases branching off from incoherence have been constructed by using different singular perturbation techniques. In particular, a powerful technique, the Chapman-Enskog method, has been used to study in detail these synchronization transitions for a bimodal frequency distribution which has a very rich phase diagram.

A. The nonlinear Fokker-Planck equation
The result of adding white-noise forcing terms to the mean-field Kuramoto model is the system of stochastic differential equations Here the i 's are independent white-noise stochastic processes with expected values Introducing the order parameter ͑4͒, the model equations ͑23͒ and ͑24͒ can be written as The Fokker-Planck equation for the one-oscillator probability density ͑ , , t͒ corresponding to this stochastic equation is provided the order parameter re i is a known function of time and we ignore the subscript i. In the limit N → ϱ and provided all oscillators are initially independent, we can derive Eq. ͑5͒: the periodicity condition, ͑ +2 , , t͒ = ͑ , , t͒, and an appropriate initial condition for ͑ , ,0͒. In most works ͑an exception is , the natural frequency distribution g͑͒ is a non-negative even function, to be considered later.

B. Linear stability analysis of incoherence
The trivial solution of the nonlinear Fokker-Planck equation, 0 =1/͑2͒, with order parameter r = 0, represents incoherent or nonsynchronized motion of all oscillators. A natural method for studying how synchronized phases with r Ͼ 0 may branch off from incoherence is to analyze its linear stability as a function of the parameters of the model, and then construct the possible solutions bifurcating from it. The first results were obtained by Strogatz and Mirollo ͑1991͒. They studied the linear stability problem setting =1/͑2͒ + ͑ , t ; ͒ with ͑ , t ; ͒ = exp͑t͒͑ , ͒ in Eqs. ͑26͒ and ͑8͒, and then neglecting terms nonlinear in , Incoherence is linearly stable as long as Re Ͻ0, and it becomes unstable if some admissible has positive real part. The periodicity condition implies = ͚ n=−ϱ ϱ b n ͑͒e in , which, inserted in Eq. ͑30͒, yields ͑ + in + n 2 D͒b n = K 2 ͑␦ n,1 + ␦ n,−1 ͒͗1,b n ͘.

͑32͒
Here we have used the relation b −n = b n , and the scalar product defined in Eq. ͑20͒. As in the Kuramoto model, the numbers n ͑͒ = − Dn 2 − in, n = ± 1, ± 2, ... , ͑33͒ with belonging to the support of g͑͒, form the continuous spectrum relevant to the linear stability problem. Note that the continuous spectrum lies to the left side of the imaginary axis when D Ͼ 0 ͓b 0 = 0 because of the normalization condition ͑31͔͒. Then the eigenvalues ͑33͒ have negative real parts and therefore correspond to stable modes. The case n = ± 1 is special for two reasons. First, the right-hand side of Eq. ͑32͒ does not vanish, and thus b 1 = ͑K /2͒͗1,b 1 ͘ / ͑ + i + D͒. Then, provided ͗1,b 1 ͘ 0, we find the following equation for ͑Strogatz and Mirollo, 1991͒: The solutions of this equation are the eigenvalues for the linear stability problem in Eq. ͑30͒. Clearly, they are independent of . Since the continuous spectrum lies on the left half plane, the discrete spectrum determines the linear stability of the incoherence. Second, the nonlinear Fokker-Planck equation and therefore the linear stability equation ͑30͒ are invariant under the reflection symmetry, → −, → −, assuming g͑͒ to be even. This implies that two independent eigenfunctions exist, corresponding to each simple solution of Eq. ͑34͒, Note that these two linearly independent eigenfunctions are related by the reflection symmetry. When is real, these eigenfunctions are complex conjugates of each other. When is a multiple solution of Eq. ͑34͒, the eigenvalue is no longer semisimple ͑Crawford, 1994͒.

C. The role of g͑͒: Phase diagram of the Kuramoto model
The mean-field Kuramoto model for infinitely many oscillators may have different stable solutions ͑also called phases͒ depending on the natural frequency distribution g͑͒, the values of the coupling strength K, and the diffusion constant D. Many phases appear as stable solutions bifurcating from known particular solutions, which lose their stability at a critical value of some parameter. The trivial solution of the nonlinear Fokker-Planck equation is incoherence, and therefore much effort has been devoted to studying its stability properties, as a function of K, D, and parameters characterizing g͑͒. As explained in Sec. II, we can always consider the first moment of g͑͒ to be equal to zero, shifting the oscillator phases if necessary. Most of the work reported in the literature refers to an even function, g͑͒, g͑−͒ = g͑͒. In addition, if g͑͒ has a single maximum at = 0, we call it a unimodal frequency distribution. For general even unimodal frequency distributions, Strogatz and Mirollo ͑1991͒ proved that the eigenvalue equation ͑34͒ has at most, one solution, which is necessarily real, and it satisfies + D Ͼ 0. Explicit calculations can be carried out for discrete ͓g͑͒ = ␦͔͑͒ and Lorentzian ͓g͑͒ = ͑␥ / ͒ / ͑ 2 + ␥ 2 ͔͒ frequency distributions. We find =−D − ␥ + K / 2, with ␥ = 0 for the discrete distribution. Clearly, incoherence is linearly stable for points ͑K , D͒ above the critical line D =−␥ + K / 2, and unstable for points below this line; see Fig. 1͑c͒. In terms of the coupling strength, incoherence is linearly stable provided that K Ͻ K c ϵ 2D +2␥, and unstable otherwise. This conclusion also holds for D = 0 in the general unimodal case, for which Strogatz and Mirollo ͑1991͒ recovered Kuramoto's result K c =2/͓g͑0͔͒ ͑K c =2␥ in the Lorentzian case͒. When D = 0, the stability analysis is complicated by the fact that the continuous spectrum lies on the imaginary axis.
For even or asymmetric multimodal frequency distributions, the eigenvalues may be complex ͑Bonilla et . The simple discrete bimodal distribution g͑͒ = ͓␦͑ − 0 ͒ + ␦͑ + 0 ͔͒ / 2 has been studied extensively ͑Bonilla et Bonilla, Pérez-Vicente, and Spigler, 1998͒. In this case, Eq. ͑34͒ has two solutions, ± =−D + ͓K ± ͱ K 2 −16 0 2 ͔ / 4. The stability boundaries for the incoherent solution can be calculated by equating to zero the greater of Re + and Re − . The resulting phase diagram on the plane ͑K , D͒ is depicted in Fig. 2 ͑Bonilla et al., 1992͒. When the coupling is small enough ͑K Ͻ 2D͒, incoherence is linearly stable for all 0 , whereas it is always unstable when the coupling is sufficiently strong, K Ͼ 4D. For intermediate couplings, 2D Ͻ K Ͻ 4D, incoherence may become unstable in two different ways. For 0 Ͻ D, ± are real and incoherence is linearly stable provided that K Ͻ K c =2D͓1+͑ 0 / D͒ 2 ͔, and unstable when K Ͼ K c . For 0 Ͼ D, ± are complex conjugate and have zero real parts at K c =4D.
What happens in regions of the phase diagram where incoherence is unstable? Typically there appear stable solutions with r Ͼ 0, which correspond to synchronized phases. As discussed below, their study has been based on bifurcation theory, for parameter values close to critical lines of the phase diagram where incoherence is neutrally stable. These analytical results are supplemented by numerical simulations of the nonlinear Fokker-Planck equation far from critical lines, or by numerical continuation of synchronized solutions bifurcating from incoherence ͑Acebrón, Perales, and Spigler, 2001͒. Besides this, Acebrón and Bonilla ͑1998͒ have developed a singular perturbation description of synchronization for multimodal g͑͒, arbitrarily far from critical lines. Their idea is to consider a g͑͒ with m maxima located at 0 ⍀ l , where 0 → ϱ, g͑͒d Ϸ͚ l=1 m ␣ l ␦͑⍀ − ⍀ l ͒d⍀, and then to use a method of multiple scales. The main result is that the solution of the nonlinear Fokker-Planck equation splits ͑at lowest order͒ into l phases, each obeying a nonlinear Fokker-Planck equation with a discrete unimodal distribution centered at ⍀ l and a coupling strength ␣ l K. Depending on the value of ␣ l K, the lth phase turns out to be either synchronized or incoherent. The overall order parameter is the weighted sum ͑with weight ␣ l ͒ of the order parameters of the corresponding phases. If first-order terms are included, the method of multiple scales describes incoherence, as well as oscillator synchronization for multimodal frequency distributions, rather well. These results hold for general frequency distributions ͑both with or without reflection symmetry͒ and for relatively low values of 0 ͑Acebrón and Bonilla, 1998͒. Related work on multimodal frequency distributions include that of Acebrón, Perales, and An interesting open problem is to generalize the method of Acebrón and Bonilla ͑1998͒ so that both the location and the width of the peaks in g͑͒ are taken into account.

D. Synchronized phases as bifurcations from incoherence, D 0
At the parameter values for which incoherence ceases to be linearly stable, synchronized phases ͓stable solutions of the nonlinear Fokker-Planck equation, ͑ , , t͒, with r Ͼ 0͔ may bifurcate from it. In this rather technical subsection, branches of these bifurcating solutions will be constructed in the vicinity of the bifurcation point by means of the Chapman-Enskog method; see Bonilla ͑2000͒. We shall study the Kuramoto model with the discrete bimodal natural frequency distribution, whose phase diagram is depicted in Fig. 2. The stability boundaries in this rich phase diagram separate regions where incoherence becomes unstable, undergoing a transition to either a stationary state, when 0 Ͻ D and K Ͼ K c =2D͓1+͑ 0 / D͒ 2 ͔, or to an oscillatory state, when 0 Ͼ D and K Ͼ K c =4D. The bifurcating solutions are as follows: ͑1͒ When 0 Ͻ D / ͱ 2, the synchronized phases bifurcating from incoherence are stationary and stable. The bifurcation is supercritical, hence the synchronized phases exist for K Ͼ K c .
An unstable branch of synchronized stationary solutions bifurcates for K Ͻ K c , reaches a limit point at a smaller coupling constant, and there coalesces with a branch of stable stationary solutions having larger r.
͑3͒ When 0 Ͼ D, the synchronized phases bifurcating from incoherence are oscillatory and the corresponding order parameter is time periodic. Two branches of solutions bifurcate supercritically at K c =4D, a branch of unstable rotating waves and a branch of stable standing waves.
͑4͒ At the special point 0 = D / ͱ 2 and K c =3D, the bifurcation to stationary solutions changes from supercritical to subcritical. Near this point, the bifurcation analysis can be extended to describe analytically how the subcritical branch of stationary solutions turns into a branch of stable solutions at a limit point.
͑5͒ At the special point 0 = D and K c =4D, which can be called the tricritical point, a line of Hopf bifurcations coalesces with a line of stationary bifurcations and a line of homoclinic orbits. The study of the corresponding O͑2͒-symmetric Takens-Bogdanov bifurcation shows how the oscillatory branches die at a homoclinic orbit of an unstable stationary solution.
The Chapman-Enskog method is flexible enough to analyze all these bifurcations and, at the same time, simpler than alternatives such as constructing the center manifold ͑Crawford, 1994͒. Other than at the two special bifurcation points, the simpler method of multiple scales explained in Appendix B yields the same results ͑Bonilla et al., 1992; Bonilla, Pérez-Vicente, and Spigler, 1998͒. The Chapman-Enskog method ͑Chapman and Cowling, 1970͒ was originally employed by Enskog ͑1917͒ in the study of the hydrodynamic limit of the Boltzmann equation. It becomes the averaging method for nonlinear oscillations ͑Boltzmann and Mitropolsky, 1961͒ and is equivalent to assuming a center manifold in bifurcation calculations ͑Crawford, 1994͒.
FIG. 2. Linear stability diagram for the incoherent solution 0 =1/͑2͒ and the discrete bimodal frequency distribution in the parameter space ͑K / D , 0 / D͒. 0 is linearly stable to the left of the lines K =4D, 0 Ͼ D ͑where Hopf bifurcations take place͒ and K / ͑2D͒ =1+ 0 2 / D 2 , 0 Ͻ D ͓where one solution of Eq. ͑34͒ becomes zero͔. To the right of these lines, the incoherent solution is unstable. At the tricritical point, K =4D, 0 = D, two solutions of Eq. ͑34͒ become simultaneously zero. The dashed line separates the region where eigenvalues are real ͑below the line͒ from that where they are complex conjugate ͑above the line͒. From Bonilla, Pérez-Vicente, and Spigler, 1998.

Bifurcation of a synchronized stationary phase
Let 0 Ͻ D and consider K close to its critical value K c =2D͓1+͑ 0 / D͒ 2 ͔. The largest eigenvalue satisfies ϳ͑K − K c ͒ / ͓2͑1− 2 / D 2 ͔͒ as K → K c . As indicated in Eq. ͑35͒, there are two eigenfunctions associated with this eigenvalue, e i / ͑D + i͒ and its complex conjugate. Except for terms decaying exponentially fast in time, the solution of the linearized stability problem at K = K c is therefore Ae i / ͑D + i͒ + c.c., where A is a constant and c.c. means complex conjugate of the preceding term. Let us suppose now ͑and justify later͒ that K = K c + 2 K 2 , where is a small positive parameter. The probability density corresponding to initial conditions close to incoherence will have the form ϳ͑2͒ −1 + A͑͒e i / ͑D + i͒ + c.c. The correction to incoherence will be close to the solution of the linearized stability problem, but now we can assume that the complex constant A varies slowly with time. How slowly? The linearized solution depends on time through the factor e t , and = O͑K − K c ͒ = O͑ 2 ͒. Thus we assume = 2 t. The probability density can then be written as ͑,,t;͒ ϳ 1 The corrections to 1 / ͑2͒ can be telescoped into an exponential with a small argument, which ensures that the probability density is always positive. Typically, by the exponential ansatz, the parameter region, where the asymptotic expansion is a good approximation to the probability density, is widened. The functions n and n are linked by the relations and so on. They depend on a fast scale t corresponding to stable exponentially decaying modes, and on a slow time scale through their dependence on A. All terms in Eq. ͑36͒ that decrease exponentially in time will be omitted. In Eq. ͑36͒, the slowly varying amplitude A obeys the equation The functions F ͑n͒ ͑A , Ā ͒ are determined from the conditions that n or n be bounded as t → ϱ ͑on the fast time scale͒, for fixed A, and periodic in . Moreover, they cannot contain terms proportional to the solution of the linearized homogeneous problem, e ±i / ͑D ± i͒, because all such terms can be absorbed in the amplitudes A or Ā ͑Bonilla, 2000͒. These two conditions imply that The normalization condition together with Eq. ͑36͒ yields ͵ − n ͑,t,;A,Ā ͒d = 0, n ജ 2. ͑40͒ To find n , we substitute Eqs. ͑36͒ and ͑38͒ into Eq. ͑26͒ and use Eq. ͑39͒ to simplify the result. This yields the following hierarchy of linear nonhomogeneous equations: and so on. Clearly, 1 = A͑t͒e i / ͑D + i͒ + c.c. obeys the linearized stability problem ͑30͒ with =0, L 1 = 0 up to terms of order . Thus it is not obvious that each linear nonhomogeneous equation of the hierarchy has a bounded periodic solution. What is the necessary solvability condition to ensure that the linear nonhomogeneous equation, has a solution of the required features? To answer this question, we assume that, in fact, Eq. ͑43͒ has a bounded periodic solution of the form n = Pe i +¯. Then P is given by from which we obtain the nonresonance condition first found by Bonilla et al. ͑1992͒. Note that we obtain Eq. ͑45͒ even when ͗1,P͘ 0. e i times the term proportional to ͗1,P͘ in Eq. ͑44͒ is a solution of L = 0 and should therefore be absorbed in the definitions of A and Ā . This implies Eq. ͑39͒. Inserting 1 into the right side of Eq. ͑41͒, we find whose solution is We see that 1 contains odd harmonics and 2 contains even harmonics ͑a possible -independent term is omitted in 2 because of the normalization condition͒. This is actually true in general: 2n contains harmonics e i2j , j =0, ±1, ... , ±n, and 2n+1 contains harmonics e i͑2j+1͒ , j =−͑n +1͒ , ... ,n. The nonlinearity of the Fokker-Planck equation is responsible for the appearance of resonant terms in the equations for 2n+1 , which should be eliminated through the terms containing F ͑2n͒ . Then, we can set F ͑2n+1͒ = 0 and we only need the scaling K − K c = O͑ 2 ͒. The ultimate reason for these cancellations is, of course, the O͑2͒ symmetry of our problem, i.e., reflection symmetry and invariance under constant rotations, → + ␣ ͑Crawford, 1994͒. Similarly, the nonresonance condition for Eq. ͑42͒ yields Keeping this term in Eq. ͑38͒, we obtain dA / d ϳ F ͑0͒ , which is a reduced equation with the stationary solution ͉A͉ = ͱ͑K 2 D /4͓͒4+͑ 0 / D͒ 2 ͔ / ͓1−2͑ 0 / D͒ 2 ͔. The corresponding order parameter is which was obtained by Bonilla et al. ͑1992͒ using a different procedure. The solution ͑49͒ exists for K Ͼ K c ͑supercritical bifurcation͒ provided that 0 Ͻ D / ͱ 2, whereas it exists for K Ͻ K c ͑subcritical bifurcation͒ when 0 Ͼ D / ͱ 2. The amplitude equation ͑38͒ implies that the supercritical bifurcating solution is stable and that the subcritical solution is unstable. We can describe the transition from supercritical to subcritical bifurcation at 0 = D / ͱ 2, K c =3D, by evaluating F ͑4͒ and adding it to the right-hand side of Eq. ͑38͒. The result is see Appendix C. The stationary solutions of this equation are the stationary synchronized phases. We see that stable phases bifurcate supercritically for K 2 Ͼ 0 if K 2 Ͼ 2 ͱ 2 2 , whereas a branch of unstable stationary solutions bifurcates subcritically for This branch of unstable solutions coalesces with a branch of stable stationary synchronized phases at the limit point K 2 Ϸ −19 2 ͑7K 2 −4 ͱ 2 2 ͒ 2 / 612.

Bifurcation of synchronized oscillatory phases
The bifurcation in the case of complex eigenvalues can be easily described by the same method. The main difference is that the solution to the linearized problem is now where ⍀ 2 = 0 2 − D 2 , K c =4D, the eigenvalues with zero real part being ͑K c ͒ = ±i⍀. For the two slowly varying amplitudes, A + , A − , we assume that equations of the form hold. Following the previous method, the nonresonance conditions for Eq. ͑42͒, with 1 given by Eq. ͑51͒, yield F + ͑0͒ and F − ͑0͒ . The corresponding amplitude equations are ͑Bonilla, Pérez-Vicente, and Spigler, 1998͒ where the overdot denotes differentiation with respect to , and To analyze the amplitude equations ͑53͒, we define the new variables By using Eq. ͑53͒, we obtain for u and v the system Clearly, u = v or u =−v correspond to traveling-wave solutions with only one of the amplitudes A ± being nonzero. The case v ϵ 0 corresponds to standing-wave solutions, which are a combination of rotating and counterrotating traveling waves with the same amplitude. We can easily find the phase portrait of Eqs. ͑56͒ corresponding to ␣, ␤, and ␥ given by Eqs. ͑54͒ ͑see Fig. 3͒. Up to, possibly, a constant phase shift, the explicit solutions ͓or A + ͑͒ϵ0 and A − ͑͒ as A + ͑͒ above͔ are obtained in the case of the traveling-wave solutions, while are obtained in the case of the standing-wave solutions. Note that both standing wave and traveling wave bifur-cate supercritically with ʈr SW ʈ / r TW Ͼ 1. Re͑␤ + ␥͒ and Re ␥ are both positive when K 2 = 1, whereas the square roots in Eqs. ͑57͒ and ͑58͒ become pure imaginary when K 2 = −1. This indicates that the bifurcating branches cannot be subcritical. An analysis of the phase portrait corresponding to Eq. ͑56͒ shows that the standing-wave solutions are always globally stable, while the travelingwave solutions are unstable. This result was first pointed out by Crawford ͑1994͒.

Bifurcation at the tricritical point
At the tricritical point, K =4D, 0 = D, a branch of oscillatory bifurcating phases coalesces with a branch of stationary bifurcating phases and a branch of homoclinic orbits, in an O͑2͒-symmetric Takens-Bogdanov bifurcation point. Studying the bifurcations in the vicinity of such a point shows how the stable and unstable branches of oscillatory phases, standing wave and traveling wave, respectively, end as the coupling is changed. Analyzing transitions at the tricritical point is a little more complicated because it requires changing the assumptions on the amplitude equation ͑Bonilla, Pérez-Vicente, and Spigler, 1998; . First of all, at the tricritical point, ͗1,͑D + i͒ −2 ͘ =Re͑D + iD͒ −2 = 0. This innocent-looking fact implies that the term −F ͑0͒ ‫ץ‬ A 1 on the right-hand side of Eq. ͑42͒ dissapears in the nonresonance condition, and therefore using the same ansatz as in Eqs. ͑36͒ and ͑38͒ will not deliver any amplitude equation. Second, the oscillatory ansatz ͑51͒ breaks down FIG. 3. Phase planes ͑a͒ ͑u , v͒ and ͑b͒ ͉͑A + ͉ , ͉A − ͉͒ showing the critical points corresponding to traveling wave ͑TW͒ and standing wave ͑SW͒ solutions. From Bonilla, Pérez-Vicente, and Spigler, 1998. too, because ⍀ = 0 at the tricritical point, and the factor multiplying A − is simply the complex conjugate of the factor multiplying A + . Therefore only one independent complex amplitude exists, and we are brought back to Eq. ͑36͒. How should one proceed?
In order to succeed, one should recognize that a basic slow time scale different from does exist near the tricritical point. The eigenvalues with largest real part are along with their complex conjugates, provided that K −4D = K 2 2 , 0 − D = 2 2 with 2 Ͼ K 2 / 4. Therefore the time-dependent factors e t appearing in the solution of the linearized problem indicate that the perturbations about incoherence vary on a slow time scale T = t near the tricritical point. This leads to the Chapman-Enskog ansatz ͑,,t;͒ = 1 The equation for A is second order ͓rather than first order as in Eq. ͑38͔͒ because resonant terms appear at O͑ 3 ͒ for the first time. These are proportional to A TT = d 2 A / dT 2 . The quantities F ͑0͒ and F ͑1͒ are evaluated in Appendix D. The resulting amplitude equation is is in the scaled normal form studied by Dangelmayr and Knobloch ͓1987, cf. their equations ͑3.3͒, p. 2480͔. Following these authors, we substitute in Eq. ͑61͒, separate real and imaginary parts, and obtain the perturbed Hamiltonian system

͑63͒
Here L = R 2 T is the angular momentum, and is the potential. This system has the following special solutions: ͑i͒ The trivial solution, L =0, R = 0, which corresponds to the incoherent probability density, =1/2. This solution is stable for which exist provided that K 2 Ͼ 0 and 2 Ͼ 19K 2 / 56. These solutions bifurcate from the trivial solution at K 2 = 2 = 0. When 2 =19K 2 / 56, the branch of traveling waves merges with the steady-state solution branch. This solution is always unstable.
͑iv͒ The standing-wave solutions, L =0, R = R͑T͒ periodic. Such solutions have been found explicitly ͑see Sec. V.A of Dangelmayr and Knobloch, 1987͒. The standing waves branch off from the trivial solution at K 2 = 2 = 0, exist for 2 Ͼ 11K 2 /19Ͼ 0, and terminate by merging with a homoclinic orbit of the steady state ͑ii͒ on the line 2 =11K 2 /19 ͓see Eq. ͑5.8͒ of Dangelmayr and Knobloch, 1987͔. This solution is always stable.
All these results are depicted in Fig. 4, which corresponds to Fig. 4, IV, in the general classification of stability diagrams reported by Dangelmayr and Knobloch ͑1987, p. 267͒. For fixed 0 Ͼ D, the bifurcation diagram near the tricritical point is depicted in Fig. 5. Note that Eq. ͑59͒ yields, to leading order, ͑,,t;͒ ϳ 1 and hence re i ϳ Re −i / ͑2D͒. It follows that r ϳ R / ͑2D͒ and ϳ −, which shows that, essentially, the In closing, the Chapman-Enskog method can be used to calculate any bifurcations appearing for other frequency distributions and related nonlinear Fokker-Planck equations. As discussed in Sec. II.B, nontrivial extensions are needed in the case of the hyperbolic limit, D → 0 + .

IV. VARIATIONS OF THE KURAMOTO MODEL
We have seen in the preceding section how the longrange character of the coupling interaction in the Kuramoto model allows us to obtain many analytical results. Yet one might ask how far these results extend beyond the mean-field limit in finite dimensions. Also, one might ask how synchronization effects in the Kuramoto model are modified by keeping long-range interactions but including additional sources of quenched disorder, multiplicative noise, or time-delayed couplings. Unfortunately, many of the analytical techniques developed in the preceding section hardly cover such new topics. In particular, the treatment of short-range couplings ͑oscillators embedded in a lattice with nearest-neighbor inter-actions͒ presents formidable difficulties at both the analytical and the numerical level. This challenges our current understanding of the mechanisms lying behind the appearance of synchronization. The next sections are devoted to discussing a number of these cases. The present knowledge of such cases is still quite modest, and major work remains to be done.

A. Short-range models
A natural extension of the Kuramoto model discussed in Sec. III includes short-range interaction effects ͑Sakaguchi et al., 1987;Daido, 1988;Strogatz and Mirollo, 1988a have considered the case in which oscillators occupy the sites of a d-dimensional cubic lattice and interactions occur between nearest neighbors, where the pair ͑i , j͒ stands for nearest-neighbor oscillators and the i 's are independent random variables chosen according to the distribution g͑͒. Compared to the Kuramoto model ͑23͒, the coupling strength K does not need to be scaled by the total number of oscillators. However, convergence of the model ͑66͒ in the limit of large d requires that K scale as 1 / d. Although this model can be extended so as to include stochastic noise, i.e., finite temperature T, most of the work on this type of model has been done at T = 0. Solving the short-range version of the Kuramoto model is a hopeless task ͑except for special cases such as one-dimensional modelssee below-or Cayley-tree structures͒ due to the difficulty of incorporating the randomness in any sort of renormalization-group analysis. In short-range systems, one usually distinguishes among different synchronization regimes. Global synchronization, which implies that all the oscillators are in phase, is rarely seen except for K ϰ N → ϱ. Phase locking or partial synchronization is observed more frequently. This is the case in which a local ensemble of oscillators verify the condition i = const, for every i. A weaker situation concerns clustering or entrainment. Although this term sometimes has been used in the sense of phase locking, usually it refers to the case in which a finite fraction of the oscillators have the same average frequency i defined by There is no proof that such a limit exists. However, if it does not exist, no synchronization whatsoever is possible. The condition of clustering is less stringent than phase locking, and therefore its absence is expected to preclude the existence of phase locking. Note that the previous definitions do not exhaust all possible types of synchronized stationary solutions, such as, for instance, the existence of moving traveling-wave structures. The concept of synchronization ͑either global or partial͒ is different from the concept of phase coherence introduced in the context of the Kuramoto model in Sec. II; see Eq. ͑2͒. Phase coherence is a stronger condition than synchronization as it assumes that all phases i are clustered around a given unique value, as are their velocities i . The contrary is not necessarily true, as phases can change at the same speed ͑synchronization takes place͒ while having completely different values ͑incoherence͒. Coherence seems less general than synchronization as the former bears connection to the type of ferromagnetic ordering present in the Kuramoto model. Although this type of ordering is expected to prevail in finite dimensions, synchronization seems more appropriate for discussing oscillator models with structural disorder built in.
For short-range systems, one would like to understand several issues: • The existence of a lower critical dimension above which any kind of entrainment is possible. In particular, it is relevant to prove the existence of phase locking and clustering for large enough dimensionality and the differences between both types of synchronization.
• The topological properties of the entrained clusters and the possibility of defining a dynamical correlation length describing the typical length scale of these clusters.
• The existence of an upper critical dimension above which the synchronization transition is of the meanfield type.
• The resulting phase diagram in the presence of thermal noise.
Sakaguchi et al. ͑1987͒ have proposed some heuristic arguments showing that any type of entrainment ͑global or local͒ can occur only for d ജ 2. This conjecture is supported by the absence of entrainment in one dimension. Strogatz andMirollo ͑1988a, 1988b͒, however, have shown that no phase locking can occur at any finite dimension. As phase locking occurs in mean-field theory, this result suggests that the upper critical dimension in the model is infinite. Particularly interesting results were obtained in one dimension ͑chains of oscillators͒. In this case, and for the case of a normal distribution of natural frequencies i , it can be proven that the probability of phase locking vanishes as N → ϱ, while it is finite whenever K Ϸ ͱ N. The same result can be obtained for any distribution ͑not necessarily Gaussian͒ of independently distributed natural frequencies. The proof consists of showing that the probability of phase locking is related to the probability that the height of a certain Brownian bridge is not larger than some given value which depends on K and the mean of g͑͒. Recall that by Brownian bridge we mean a random walk described by n moves of length x i extracted from a given probability distribution, where the end-to-end distance l͑n͒ = ͚ i=1 n x i is constrained to have a fixed value for a given number n of steps. However, one of the most interesting results in these studies is the use of block renormalization-group techniques to show whether clustering can occur in finite dimensions. Nearly at the same time, Strogatz andMirollo ͑1988a, 1988b͒ and Daido ͑1988͒ presented similar arguments but leading to slightly different yet compatible conclusions. For Strogatz andMirollo ͑1988a, 1988b͒, the goal was to calculate the probability P͑N , K͒ that a cube S containing a finite fraction ␣N ͑␣ Ͻ 1͒ of the oscillators could be entrained to have a single common frequency. Following Strogatz andMirollo ͑1988a, 1988b͒, let us assume the macroscopic cluster S to be divided into cubic subclusters S k of side l, the total number of subclusters being N s = ␣N / l d which is of order N. For each subcluster k, the average frequency ⍀ k and phase ⌰ k are defined as Summing Eq. ͑66͒ over all oscillators contained in each subcluster S k we get where the sum in the right-hand side runs over all links ͑i , j͒ crossing the surface ‫ץ‬S k delimiting the region S k . Considering that there are 2dl d−1 terms in the surface ͓and hence in the sum in Eq. ͑69͔͒, this implies that If S is a region of clustered oscillators around the frequency , then, after time averaging, the limit ͑67͒ gives ͉ − ⍀ k ͉ ഛ 2dK / l for all 1 ഛ k ഛ N s . Since the ⍀ k are uncorrelated random variables, the probability that such a condition is simultaneously satisfied for all N s oscillators is p N s , p being the typical probability value such that ͉ − ⍀ k ͉ ഛ 2dK / l is satisfied for a given oscillator. This probability is therefore exponentially small when the number of subclusters N s is of order O͑N͒, Here c is a constant, and the factor N in front of the exponential is due to the number of all possible ways the cluster S can be embedded in the lattice. Therefore the probability vanishes in any dimension in the limit of large population. This proof assumes that entrained clusters have a compact structure such as a cubical shape. However, this need not necessarily be true. Had the clusters a noncompact shape ͑such as space-filling sponges or lattice-animal treelike structures͒, the proof would not hold anymore, since the number of subclusters N s does not have to scale necessarily as 1 / l d . Therefore such a result does not prevent the existence of a macroscopic entrainment in noncompact clusters.
This finding does not appear to contradict that reported by Daido ͑1988͒, who has shown the existence of a lower critical dimension d l , depending on the tails of a class of frequency distributions g͑͒. When then normalization requires ␣ Ͼ 0. Moreover, Daido considers that ␣ ഛ 2 for distributions with an infinite variance, the limiting case ␣ = 2 corresponding to the case of a distribution with finite variance ͑such as the Gaussian distribution͒. The argument put forth by Daido resorts to a similar block decimation procedure as that outlined in Eqs. ͑68͒ and ͑69͒. However he reports different conclusions from those by Strogatz andMirollo ͑1988a, 1988b͒. Having defined the subcluster or block frequencies ⍀ k and phases ⌰ k in Eq. ͑68͒, Daido shows that only for d Ͻ ␣ / ͑␣ −1͒, and in the limit when l → ϱ, does Eq. ͑70͒ yield the fixed-point dynamical equation for every k, which shows that no clustering can occur for 0 Ͻ ␣ ഛ 1 in any dimension. However, for 1 Ͻ ␣ ഛ 2, macroscopic entrainment should be observed for dimensions above d l = ␣ / ͑␣ −1͒. For the Gaussian case, ␣ = 2, entrainment occurs above d l = 2 as was also suggested by Sakaguchi et al. ͑1987͒. Numerical evidence in favor of a synchronization transition in dimension d =3 ͑Daido, 1988͒ is not very convincing. The correct value of the upper critical dimension remains to be solved. The Kuramoto model in an ultrametric tree has also been studied by Lumer and Huberman ͑1991, 1992͒. They considered a general version of the model in Eq. ͑66͒, where N oscillators sit in the leaves of a hierarchical tree of branching ratio b and L levels; see Fig. 7. The coupling among the oscillators K ij is not uniform but depends on their ultrametric distance l ij , i.e., the number of levels in the tree separating the leaves from its common ancestor, where d͑x͒ is a monotonically decreasing function of the distance. The existence of a proper thermodynamic limit requires, for d͑x͒ in the limit of large L , N, for all j. For a given value of K, as the ultrametric distance l ij increases, entrainment fades away. However, as K increases, more and more levels tend to synchronize. Therefore this model introduces in a simple way the clusterization of synchronization, thought to be relevant in the perception problem at the neural level ͑see Sec. VII.A͒. The simplest function d͑x͒ that incorporates such effects has an exponential decay d͑x͒ϳ1/a x , where the coupling strength decreases by a factor a at consecutive levels. It can be shown ͑Lumer and Huberman, 1992͒ that a cascade of synchronization events occurs In such a regime, the model displays a nice devil's staircase behavior when plotting the synchronization parameter as a function of K. This is characteristic of the emergent differentiation observed in the response of the system to external perturbations. Other topologies beyond the simple cubic lattice structure have also been considered. Niebur, Schuster, et al. ͑1991͒ analyze spatial correlation functions in a square lattice of oscillators with nearest-neighbor, Gaussian ͑the intensity of the interaction between two sites decays with their distance according to a Gaussian law͒, and sparse connections ͑each oscillator is coupled to a small and randomly selected subset of neighbors͒. Overall, they find that entrainment is greatly enhanced with sparse connections. Whether this result is linked to the supposed noncompact nature of the clusters is yet to be checked. An intermediate case between the long-range Kuramoto model and its short-range version ͑66͒ occurs when the coupling among oscillators decays as a power law, 1/r ␣ , r denoting their mutual distance. The intensity of the coupling is then properly normalized in such a way that the interaction term in Eq. ͑66͒ remains finite in the limit of large population. For the normalized case, in one dimension, it has been shown ͑Rogers and Wille, 1996͒ that a synchronization transition occurs when ␣ Ͻ ␣ c ͑K͒ with K͑␣ =0͒ =2/͓g͑ =0͔͒, corresponding to the Kuramoto model ͓see the paragraph just after Eq. ͑11͔͒. It is found that ␣ Ͻ 2 is required for a synchronizing transition to occur at finite K. Note that the same condition is found for one-dimensional Ising and XY models in order to have a finite-temperature transition. For the non-normalized case results are more interesting ͑Marodi et al., 2002͒, as they show a transition in the population size ͑rather than in the coupling constant K͒ for ␣ Ͻ d. In such a case, synchronization occurs provided that the population is allowed to grow above some critical value N c ͑K , ␣ Ͻ d͒. The relevance of such a result stems from the fact that sufficiently large threedimensional populations, interacting through a signal whose intensity decays as 1 / r 2 , such as sound or light, can synchronize whatever the value of the coupling K might be. FIG. 7. A hierarchical tree with N = 9 two-level ͑L =2͒ oscillators having branching ratio b = 3. The distance l ij between two oscillators is given by the number of levels between them and their closest common ancestor. In this example, l 12 = l 46 = l 78 =1 and l 15 = l 58 =2.
Before ending the present overview of short-range models, it is worth mentioning that the complexity of the synchronization phenomenon in two dimensions has been emphasized in another investigation conducted by Kuramoto and co-workers ͑Sakaguchi et al., 1988͒. Here they studied a model where the coupling function was slightly modified to account for an enhancement of the oscillator frequencies due to their interaction, with − /2ഛ ␣ ഛ / 2. This is a nonvariational model, since the interaction term does not correspond to the gradient of a two-body potential. The effect of the parameter ␣ is to decrease the entrainment among oscillators that have different natural frequencies, giving rise to a higher value of the critical coupling. Notably this model has been found to describe synchronization of Josephson-junction arrays ͓see Eq. ͑155͒ and the ensuing discussion͔. The parameter ␣ already has an interesting effect in the nondisordered model i =0 ͑for ␣ = 0 this is the XY model͒, where neighboring phases inside a vortex can differ greatly, thereby inducing an increase in the local frequency. The vortex acts as a pacemaker, enhancing entrainment in the surrounding medium.

B. Models with disorder
The standard Kuramoto model already has disorder built in. However, one can include additional disorder in the coupling among the oscillators. Daido has considered the general mean-field model ͑Daido, 1992b͒ where the couplings K ij are Gaussian distributed, the natural frequencies being distributed according to a distribution g͑͒ and being Gaussian noise. The A ij are real-valued numbers lying in the range ͓− , ͔ and stand for potential vector differences between sites which amount to random phase shifts. The interest of this model lies in the fact that frustration, a new ingredient beyond disorder, is introduced. Frustration implies that the reference oscillator configuration, which makes the coupling term in Eq. ͑77͒ vanish, is incoherent, i.e., i − j 0. This is due to the competing nature of the different terms involved in that sum which contribute with either a positive or a negative sign. Frustration makes the reference configuration extremely hard to find using standard optimization algorithms. Compared to the models without disorder, few studies have been devoted to the disordered case so far, yet they are very interesting as they seem to display synchronization to a glassy phase rather than to a ferromagneticlike one. Moreover, as structural disorder tends to be widespread in many physical systems, disordered models seems to be no less relevant to realistic oscillator systems than their ordered counterparts.

Disorder in the coupling: the oscillator glass model
The model ͑77͒ with disorder in the coupling parameter K ij and A ij = 0 has been the subject of recent work ͑Daido, 1987a͑Daido, , 1992b͑Daido, , 2000Stiller andRadons, 1998, 2000͒. The model equation is where the K ij 's are given by Eq. ͑78͒. When i = 0 this corresponds to the XY spin-glass model ͑Sherrington and Kirkpatrick, 1975;Kirkpatrick and Sherrington, 1978͒, introduced to mimic the behavior of frustrated magnets. It is well known that the XY spin-glass model has a transition at a critical noise intensity D c = 1 from a paramagnetic phase ͑D Ͼ D c ͒ to a spin-glass phase ͑D Ͻ D c ͒. When natural frequencies exist, a synchronization transition is expected to occur from an incoherent to a synchronized glassy phase. The most complete study ͑yet with contradictory results; see below͒ was done without noise and for a Gaussian frequency distribution g͑͒ in two different works. One of these was done by Daido ͑1992b͒, the other by Stiller and Radons ͑1998͒. Daido ͑1992b͒ conducted a detailed numerical analysis of the distribution p͑h͒ of the local field h j acting on each oscillator defined by Daido showed how, as the value of K is increased, p͑h͒ changes from a pure Gaussian distribution, whose maximum, h * , is located at 0, to a non-Gaussian function where h * becomes positive. This establishes the value of the critical coupling as K c Ϸ 8. Glassy systems are known for their characteristic slow relaxational dynamics, which lead to phenomena such as aging in correlations and response functions and violations of the fluctuation-dissipation theorem ͑Crisanti and Ritort, 2003͒. Synchronization models are expected to show similar nonequilibrium relaxation phenomena, although aging is not expected to occur in the stationary state. In particular, Daido also considered the order parameter showing that Z͑t͒ decays exponentially with time in the incoherent-phase case ͑K Ͻ K c ͒, and as a power law in time in the synchronized-phase case. These results were criticized in a subsequent work by Stiller and Radons ͑1998͒. They considered the analytical solution of the dynamics of the same model, using the path-integral for-malism of Martin, Siggia, and Rose to reduce the N-oscillator problem to a single-oscillator problem with a correlated noise.
where C͑t , s͒ is the two-times complex-valued correlation function The order of the limits in Eq. ͑83͒ is important: the t 0 → ϱ limit assures that the stationary regime is attained first, the N → ϱ limit ensures ergodicity breaking, and the t → ϱ limit measures the equilibrium order within one ergodic component. Measurements of q reveal a neat transition across K c =24 ͑see Fig. 3

in Stiller and
Radons, 1998͒, nearly three times larger than the value reported by Daido ͑see Fig. 2 in Daido, 1992b͒. The origin of such a discrepancy is thus far unknown. Further work is needed to resolve this question. Before finishing this section, let us mention that the study of the model in Eq. ͑79͒ with site-disordered couplings has been considered by Bonilla et al. ͑1993͒ for the Kuramoto model with the van Hemmen-type interactions. Such interactions are characterized by the coupling parameters K ij = K 0 / N + ͑K 1 / N͒͑ i j + i j ͒, where i , i are random independent identically distributed quenched variables which take values ±1 with probability 1 / 2. Such a model is exactly solvable, so the resulting phase diagram can be obtained analytically. These au-thors found several phases, depending on the ratio between K 0 and K 1 , i.e., incoherence, synchronization, spin-glass phase, and mixed phase ͑where oscillators are partially coherently synchronized and partially in phase opposition͒.

The oscillator gauge glass model
The oscillator gauge glass model has seldom been studied and is included here for completeness. It corresponds to the model in Eq. ͑77͒, where K ij = K / N and frustration arises solely from the random phase shifts A ij ͓− , ͔. Note that this term by itself is enough to induce frustration. For instance, when A ij =0 or with probability 1 / 2, such a model coincides with the previous oscillator glass model where K ij = ± 1 with identical probability. In general, for nondisordered models ͑in the absence of natural frequencies͒, and in the absence of an external magnetic field, the vector potential components around a plaquette add to zero. In the presence of a field, this term is an integer number times the value of the quantized flux. Therefore the present oscillator gauge glass model in two dimensions can be seen as a simplified version of overdamped Josephson-junction arrays. Let us just mention a study by Park et al. ͑1998͒ of the phase diagram of the oscillator gauge glass model, in which the stationary states were taken as the equilibrium states of the corresponding Boltzmann system. The key technique was to map such a model into an equilibrium gauge glass model. However, as already explained in Sec. VI.C, this approach does not guarantee a full characterization of the stationary solutions.

C. Time-delayed couplings
One of the most natural and well-motivated extensions of the Kuramoto model concerns the analysis of time-delayed coupling among oscillators. In biological networks, electric signals propagate along neural axons at a finite speed. Thus delays due to transmission are natural elements in any theoretical approach to information processing. Time delays can substantially change the dynamical properties of coupled systems. In general, the dynamic behavior becomes much richer and, on occasion, even surprising. One might think that time delays tend to break coherence or to make it difficult in populations of interacting units, but this is not always the case. An important example has been the focus of intense research in the past few years: synchronization in chaotic systems. It has been proved that a delayed coupling is the key element for anticipating and controlling the time evolution of chaotic coupled oscillators by synchronizing their dynamics ͑Pikovsky et Voss, 2001͒. In excitable systems the situation is quite similar. It has been reported that delayed couplings may favor the existence of rapid phase-locked behavior in networks of integrate-and-fire oscillators ͑Gerstner, 1996͒.
How important are time delays for a population of coupled phase oscillators? This depends on the ratio of the time delay to the natural period of a typical oscilla-tor. In general, the delay should be kept whenever the transmission time lag is much longer than the oscillation period of a given unit ͑Izhikevich, 1998͒. On the other hand, the delay can be neglected whenever it is comparable to the period of the oscillators.
Let us start by discussing the dynamic properties of short-range coupled Kuramoto models with timedelayed couplings ͑Schuster and Wagner, 1989;Niebur, Schuster, and Kammen 1991;Nakamura et al., 1994͒: Here the sum is restricted to nearest neighbors, and is a constant time delay. In this model, each oscillator interacts with its neighbors in terms of the phase that they had at the time they sent a synchronizing signal. In a simple system of two oscillators, Schuster and Wagner ͑1989͒ found the typical fingerprint of delays in several relevant differences between this model and the standard Kuramoto model. When = 0, one synchronization state is stable. However, for nonzero delays and given values of and K, there are multiple stable solutions of Eq. ͑85͒ with more than one synchronization frequency and different basins of attraction. In a remarkable application, this model has been used recently to explain the experimentally observed oscillatory behavior of a unicellular organism ͑Takamatsu et al., 2000͒. Keeping these results in mind, it is not difficult to imagine what will happen for large ensembles of oscillators. In a 1D system, Nakamura et al. ͑1994͒ have shown that complex structures emerge spontaneously for N → ϱ. They are characterized by many stable coexisting clusters, each formed by a large number of entrained units, oscillating at different frequencies. Such structures have also been observed for a distribution of time delays ͑Zanette, 2000͒. Another peculiarity of the model has been observed for large coupling intensities and large delays. In this regime, the system exhibits frequency suppression resulting from the existence of a large number of metastable states. In two dimensions, Niebur, Schuster, and Kamman ͑1991͒ have shown that the system evolves towards a state with the lowest possible frequency, selected from all the possible solutions of Eq. ͑85͒, for given values of K and . More recently, Jeong et al. ͑2002͒ have shown that distance-dependent time delays induce various spatial structures such as spirals, traveling rolls, or more complex patterns. They also analyzed their stability and the relevance of initial conditions to select these structures.
The mean-field model has been studied theoretically by analyzing the corresponding nonlinear Fokker-Planck equation ͑Luzyanina, 1995; Yeung and Strogatz, 1999;Choi et al., 2000͒. As in the case of short-range couplings, the delay gives rise to multistability even for identical oscillators and without external white noise. The phase diagram ͑K , ͒ contains regions where synchronization is a stable solution of the dynamic equations, other regions where coherence is strictly forbidden, and still others where coherent and incoherent states coexist. The existence of all these regimes has been corroborated by simulations ͑Kim et al., 1997a;Yeung and Strogatz, 1999͒. The same qualitative behavior has been found for nontrivial distributions of frequencies. As in the standard Kuramoto model, in noisy systems there is a critical value of the coupling K c , above which the incoherent solution is unstable. The value of K c depends on the natural frequency distribution, the noise strength, and the delay. Yeung andStrogatz ͑1999͒ andChoi et al. ͑2000͒ carried out a detailed analysis of the bifurcation at K c for the nonlinear Fokker-Planck equation corresponding to Eq. ͑85͒. Kori and Kuramoto ͑2001͒ studied the same problem for more general phase oscillator models.

D. External fields
A natural extension of the original Kuramoto model is to add external fields, which gives rise to a much richer dynamical behavior. External fields can model the external current applied to a neuron so as to describe the collective properties of excitable systems with planar symmetry. For other physical devices, such as Josephson junctions, a periodic external force can model an oscillating current across the junctions. The Langevin equation governing the dynamics of the extended model is Shinomoto and Kuramoto ͑1986͒ studied the case h i = h, for every i. They analyzed the nonlinear Fokker-Planck equation associated with Eq. ͑86͒ and found two different regions of the phase diagram: a region of timeperiodic physical observables, and a region of stable stationary synchronized states. Sakaguchi ͑1988͒ replaced the last term in Eq. ͑86͒ by h sin͑ i − f t͒, where f is the frequency of the external force. Note that, by defining i = i − f t, we obtain again Eq. ͑86͒ for i , but with natural frequencies i − f . Therefore introducing a time-periodic external force amounts to modifying the statistical properties of the natural frequency distribution g͑͒. In general, there will be a competition between forced entrainment and mutual entrainment. When h is large, the oscillators tend to be entrained with the external force. On the other hand, when h is small there will be a macroscopic fraction of the population mutually entrained displaying a synchronous collective motion. Such a motion occurs with a frequency equal to the mean natural frequency of the population.
Arenas and Pérez-Vicente ͑1994a͒ studied the phase diagram of a Kuramoto model with a distribution of random fields, f͑h͒. They solved the nonlinear Fokker-Planck equation through a generating functional of the order parameters, and found analytical expressions thereof, which fully agreed with the numerical simulations. When f͑h͒ is centered at h = 0, they found a phase transition between an incoherent state with r = 0 and a synchronized state, which is similar to the transition in the static model as well as in the model of identical os-cillators ͑Arenas and Pérez-Vicente, 1993͒. The only difference here is that the critical value K c is larger. This is reasonable because larger values of the coupling strength are needed to counteract the effect of rotation due to the frequency distribution. When the field distribution is not centered at the origin, an effective force arises. This makes r 0 for any value of the ratio K / D ͑r ϰ ͗h͘ for small K / D and for ͗h͘ Ӷ 1͒, thereby precluding phase transitions, as in the static case.
Densities having time-dependent solutions with nonzero order parameters ͑Shinomoto and Arenas and Pérez-Vicente, 1994a͒ exist, provided that ͐hf͑h͒dh 0. Acebrón and Bonilla ͑1998͒ studied these solutions using the two-timescale asymptotic method mentioned in Sec. III.C. The probability density splits into independent components corresponding to different peaks in the multimodal frequency distribution. Each density component evolves toward a stationary distribution in a comoving frame, rotating at a frequency corresponding to the appropriate peak in g͑͒. Therefore the overall synchronous behavior can be determined by studying the synchronization of each density component ͑Acebrón and  Inspired by biological applications, Frank et al. ͑2000͒ analyzed the behavior of phase oscillators in the presence of forces derived from a potential with various Fourier modes. They studied in detail the transition from incoherence to phase locking. Finally, Coolen and Pérez-Vicente ͑2003͒ studied the case of identical oscillators with disordered couplings and subject to random pinning fields. This system is extremely frustrated and several spin-glass phases can be found in it. The equilibrium properties of such a model depend on the symmetries of the pinning-field distribution and on the level of frustration due to the random interactions among oscillators.

E. Multiplicative noise
To end this section, let us discuss the effect of multiplicative noise on the collective properties of phase oscillators. As far as is known, there have been two different approaches to this problem. Park and Kim ͑1996͒ studied the following rather complex version of the Kuramoto model: Here i is a zero-mean delta-correlated Gaussian noise with unit variance, measures the intensity of the noise, and is an integer. In this model, the phase oscillators are subject to an external pinning force and therefore they represent excitable units. Thus Eq. ͑87͒ describes the effect of multiplicative noise on a population of excitable units. Through analytical and numerical studies of the nonlinear Fokker-Planck equation, Park and Kim ͑1996͒ found the phase diagram of the model ͑h , ͒ for different values of . For identical oscillators, there are new phases in which two or more stable clusters of syn-chronized oscillators can coexist. This phenomenology is strictly induced by the multiplicative noise, without requiring time delays or high Fourier modes in the coupling. Kim et al. ͑1996;Kim, Park, Doering, and Ryu, 1997͒ supplemented Eq. ͑87͒ with additive noise. Rather surprisingly, the additive noise tended to suppress the effects triggered by the multiplicative noise, such as the bifurcation from one-cluster phase to the two-cluster state. The new phase diagram exhibited a very rich behavior, with interesting nonequilibrium phenomena such as reentrant transitions between different phases. The interesting physics induced by the combination of both types of noise caused Kim, Park, and Rhu ͑1997b͒ to propose a new mechanism for noise-induced current in systems under symmetric periodic potentials.
Reimann et al. ͑1999͒ tackled the problem from a different standpoint. They considered the standard meanfield equation ͑23͒ with a nonequilibrium Gaussian noise, characterized by where D͑͒ = D 0 + D 1 cos͑͒ and D 0 ജ D 1 ജ 0. The authors considered only one Fourier mode to make their analysis simpler, and studied the particular case D 0 = D 1 = Q / 2. Under such conditions, for identical oscillators and for arbitrarily weak couplings, time-dependent oscillatory synchronization appears for a certain value of the ratio K / Q, via spontaneous symmetry breaking. By means of numerical simulations, Reimann et al. also studied the effect of multiplicative noise in systems with short-range couplings. In such a case, even more complex behavior, including hysteretic phenomena and negative mobility, was found.
Recently, and only for identical oscillators, Kostur et al. ͑2002͒ studied Eq. ͑86͒ with h i = −1 for every i, with both additive and symmetric dichotomic noise. For the latter noise, they found a complex phase diagram with five different regions: incoherence, bistability, phase locking, hysteretic phenomena, and an oscillatory regime.

V. BEYOND THE KURAMOTO MODEL
Many generalizations of the Kuramoto model have been proposed to analyze synchronization phenomena in more complex situations. First of all, periodic coupling functions, which contain more harmonics than the simple sine function considered by Kuramoto, have been proposed. Moreover, there are oscillator models described by two angles ͑tops models͒, or by phase and amplitude ͑amplitude oscillators͒. Finally, the dynamical behavior described by the Kuramoto model changes if phase oscillators possess inertia, which makes synchronization harder but the emergence of spontaneous phase oscillations easier.

A. More general periodic coupling functions
An immediate generalization of the mean-field Kuramoto model is given by where h͑͒ is a general 2 periodic coupling function, and the natural frequencies i are distributed with probability density g͑͒, as usual. By shifting all phases, → + ⍀t, and selecting ⍀ appropriately, we can set ͐ −ϱ +ϱ g͑͒d = ͐ − h͑͒d = 0, without loss of generality. What can we say about oscillator synchronization in this more general context?
A particularly successful theoretical approach is due to Daido ͑1992a, 1993aDaido ͑1992a, , 1993bDaido ͑1992a, , 1994Daido ͑1992a, , 1995Daido ͑1992a, , 1996a, who generalized Kuramoto's idea on the order parameter and the partially synchronized state. Let us assume that the oscillators are phase locked at a common frequency e ͑t͒ = . Suppose that there exist the following order parameters: If the coupling function is expanded in Fourier series as then a simple calculation shows that we can write Eq. ͑89͒ in the following form: provided we define the order function H as Note the similarity between Eq. ͑92͒ and Kuramoto's expression, Eq. ͑3͒. Daido ͑1992a͒ derived a selfconsistent functional equation for H͑͒, assuming that H͑͒ possesses only one maximum and one minimum inside the interval − ഛ ഛ . This equation always has the trivial solution H͑͒ = 0, which corresponds to incoherence. Nontrivial solutions describe synchronized states of the oscillator population. Furthermore, one expects that, for K sufficiently large, the onset of mutual entrainment is a bifurcation of a nontrivial order function from incoherence.
Keeping this in mind, Daido considered several special cases ͑Daido, 1992a͑Daido, , 1993a. When the Fourier series of the coupling function contains only odd harmonics, there is spontaneous synchrony above the critical value of the coupling constant, This shows that K c depends on the first harmonic but not on the higher modes. Notice that the Kuramoto model belongs to this family of models ͓in fact h͑͒ = sin ͔ and the value given by Eq. ͑94͒ is consistent with the results given in Sec. II. The bifurcation to the synchronized state is supercritical provided that gЉ͑ e ͒ Ͻ 0 and h 1 s Ͼ 0. These results were confirmed by numerical simulations.
Daido ͑1994͒ worked out a bifurcation theory for the order function whose L 2 norm, ʈHʈϵ͓͐ − H 2 ͑͒d͔ 1/2 , can be represented in a bifurcation diagram as a function of the coupling strength K. Near the critical coupling, ʈHʈ ϰ ͑K − K c ͒ ␤ , which defines the critical exponent ␤. It turns out that the Kuramoto model describes the scaling behavior of a reduced family of phase oscillators for which ␤ =1/2, whereas ␤ = 1 for the vast majority of coupling functions ͑Daido, 1994, 1996b͒. This fact suggests the existence of different classes of universality. Crawford ͑1995͒ confirmed most of these findings by means of standard bifurcation theory, while Balmforth and Sassi ͑2000͒ gave a simple mode-coupling explanation for the different scalings in an example, where the only nonzero harmonics are h 1 s = 1 and h 2 s = . The previous theory can be generalized to order functions with several peaks in the interval ͑− , ͒ ͑Daido, 1995, 1996a͒. In this case, the oscillators may choose among different coexisting phase-locking states, and their resulting dynamical behavior is more complex than the standard case, where only one type of entrainment is possible ͑Daido, 1996a͒. When the order function has more than one local maximum and a local minimum, there is an overlap region with at least two branches where HЈ͑͒ Ͼ 0. Oscillators whose frequency lies in the overlap region may synchronize to the phase of one stable branch, depending on their initial condition. The number of such states is exponentially large, and the entropy per oscillator is an appropriate order parameter which characterizes the corresponding macroscopic state ͑Daido, 1996a͒.
Work carried out by other authors confirms the predictions made by means of the order-function formalism. For instance, Hansel et al. ͑1993b͒ considered coupling functions with two Fourier modes and a free parameter that controls the attractive/repulsive character of the interaction among oscillators. Besides the usual states typical of the Kuramoto model, the incoherent state and the synchronized state, they found more complex dynamical occurrences, according to suitable values of the control parameter. For instance, there was switching between two-cluster states connected by heteroclinic orbits. Each cluster contained a group of phase-locked oscillators running at a common frequency. Other studies discussing a large variety of clustering behavior for coupling functions having a large number of Fourier modes were conducted by Golomb et al. ͑1992͒, Okuda ͑1993͒, and Tass ͑1997͒.
Bonilla, Pérez-Vicente, Ritort, and Soler ͑1998͒ studied singular coupling functions such as h͑͒ = ␦Ј͑͒ , ␦͑͒ , sin͑͒, etc. which possess infinitely many nonvanishing Fourier modes. They showed that the dynamics of these models can be exactly solved using the moment approach discussed in Sec. VI.C. Considering the generating function for the moments, where the drift velocity v͑x , t͒ is defined by The moment-generating function and the oneoscillator probability density are related by and synchronization is possible only at zero temperature ͑D =0͒. Finally when h͑͒ = sin͑͒, the equation for can be recast as a pair of coupled nonlinear partial differential equations. Their stationary solutions can be evaluated explicitly in terms of elliptic functions, and therefore the associated bifurcation diagram can be constructed analytically for all K's. In this case, multiple solutions bifurcate from incoherence for different values of the coupling strength. It should be mentioned that these authors established a link between their approach and Daido's order-function method. One-dimensional chains of phase oscillators with nearest-neighbor interactions ͑and also beyond nearest neighbors͒ and arbitrary coupling function have also been studied recently by Ren and Ermentrout ͑2000͒. Given the complexity of such a problem, they studied general properties of the model such as the conditions required to ensure existence of phase-locking solutions. Numerical examples were provided to confirm their theoretical predictions.

B. Tops models
The Kuramoto model deals with interacting units which behave as oscillators that are described by only one variable, i.e., their phase. However, it is not difficult to imagine other cases in which the mutually interacting variables are not oscillators, but rather classical spins described by azimuthal and polar angles. Ritort ͑1998͒ introduced a tops model, and solved the associated phase diagram in the mean-field case.
Although the name "top" is a misnomer ͑tops are described by three Euler angles, rather than only two͒, it is perhaps a more appropriate word than the more correct term "spin." This choice may avoid some confusion in view of the overwhelming variety of spin models existing in the literature. The tops model might have experimental relevance whenever precession motion is induced by an external perturbation acting upon the orientational degrees of freedom of a given system. It can describe the synchronized response of living beings, such as bacteria, endowed with orientational magnetic properties, or complex resonance effects in random magnets or NMR.
The model consists of a population of N Heisenberg spins or tops ͕ ជ i ,1ഛ i ഛ N͖ ͑of unit length͒ which precess around a given orientation n i at a given angular velocity i . Larmor precession of the ith top is therefore described by a vector ជ i = i n i . Moreover, the tops in the population mutually interact, trying to align in the same direction. The tops-model equations of motion are where K is the coupling strength and E͑ ជ i , ជ j ͒ is the twobody energy term. Rotational invariance symmetry requires E to be a function of the scalar product ជ i ជ j , the simplest case being the bilinear form E͑ ជ i , ជ j ͒ = ជ i ជ j . Natural frequencies ជ i are chosen from a distribution P͑ ជ i ͒. The term ជ i ͑t͒ is a Gaussian noise of zero mean and variance equal to 6D, the factor 6 arising from the three degrees of freedom. As it stands, Eq. ͑101͒ is not yet well defined because the unit length of the vector ជ i is not constant, as can be seen by multiplying both sides of the equation by ជ i . In order to avoid such a problem, it is convenient to project the tops onto the surface of a unit radius sphere. This requires us to introduce the azimuthal angles i ͓0,͔ and the polar angles i ͓0,2͒, ជ i = ͓cos͑ i ͒sin͑ i ͒ , sin͑ i ͒sin͑ i ͒ , cos͑ i ͔͒. The main technical difficulty in this approach arises from the noise term, which should be described by Brownian motion constrained on a spherical surface ͑see, for instance, Coffey et al., 1996͒. The use of spherical coordinates also requires that the natural frequency vectors ជ i be specified in terms of their modulus i , and their azimuthal and polar angles ͑ i , i ͒. Note at this point that three parameters enter into the description of the disordered units, rather than only one ͑i.e., the natural frequency i ͒ in the Kuramoto model. The ensuing equations of motion are where the functions F , F , G , G can be easily obtained by transforming the first two terms on the righthand side of Eq. ͑101͒ to spherical coordinates, and where m ជ = ͑1/N͚͒ i=1 N ជ i is the average magnetization. The latter plays the role of the order parameter in the tops model, as does the parameter r in the Kuramoto model ͓see Eq. ͑2͒ in Sec. II͔. The noise terms i , i are Gaussian correlated with variance 2D and average D cot i and 0, respectively.
The solution of the tops model poses additional mathematical difficulties compared to the Kuramoto model, but most of the calculations can be done for the simplest disordered cases. Ritort ͑1998͒ accomplished the calculations were in the presence of only orientational disorder, where i = for every i. The natural frequency distribution is here given by a function p͑ , ͒. Perhaps the easiest approach to solving such a model is that of using the moment representation ͑see Sec. VI.C͒, by introducing the set of moments This formulation allows for a simple analysis of both the stationary solutions and their stability in the ͑K = K / D , = / D͒ plane, as well as providing an efficient method for the numerical integration of the model in the limit N → ϱ. Ritort ͑1998͒ investigated several axially symmetric disorder distributions where p͑ , ͒ϵp͑͒: ͑a͒ antiferromagnetically oriented precessing frequencies with p͑͒ = ͑1/4͒␦͑ −0͒ + ͑1/4͒␦͑ − ͒, ͑b͒ precessing orientations lying in the XY plane, p͑͒ = ͑1/2͒␦͑ − /2͒, and ͑c͒ isotropic disorder p͑͒ =1/4. While case ͑a͒ corresponds to a purely relaxational model ͑the mean-field antiferromagnet͒, displaying a single synchronization transition at K = 3, models ͑b͒ and ͑c͒ show a more complex dynamical behavior. In fact, in case ͑c͒ it was found that the incoherent solution is stable for K Ͻ 3, unstable for K Ͼ 9, and stable in the region 3 Ͻ K Ͻ 9 whenever 2 Ͼ ͑12K −36͒ / ͑9−K ͒, yielding a rich pattern of dynamical behavior.
More work still needs to be done on the tops model. Especially interesting would be an experimental verification of the orientational entrainment in magnetic systems or in magnetized living cells. Here, nonlinear effects are introduced by the inertial effects induced by the Larmor precession of magnetic moments in a magnetic field.

C. Synchronization of amplitude oscillators
As long as the attraction in each oscillator to its limit cycle dominates over the coupling among the members of the population, a phase model suffices to account for a number of phenomena. This is the case of the classical Kuramoto model. On the other hand, when the coupling is strong enough, the amplitude of each oscillator may be affected, and hence a more comprehensive model is required, in which the dynamics of the amplitudes as well as those of the phases is included.
Similarly to the case of the Kuramoto model, one can consider several interactions among the oscillators, most importantly, nearest-neighbor ͑also called diffusive͒ coupling ͑Bar-Eli, 1985͒, random coupling between each oscillator and an arbitrary number of neighbors ͑Satoh, 1989͒, and "all-to-all" ͑global͒ coupling ͑Ermentrout, 1990; Matthews and Strogatz, 1990;Matthews et al., 1991͒. Incidentally, it should be observed that in several papers the term "diffusive coupling" is used to refer to an "all-to-all" coupling instead of to nearest-neighbor interaction. It is the latter, in fact, that corresponds to the Laplace operator in a discretized form.
In this section, only amplitude models subject to a global interaction are discussed. These models, due to their simplicity, have been more extensively investigated so far. For instance, a case widely studied is given by the system of linearly coupled oscillators, each near a Hopf bifurcation,

͑105͒
Here z j is the position of the jth oscillator in the complex plane, j its natural frequency, picked up from a given frequency distribution, g͑͒, and K is the coupling strength. One should notice that, when the coupling is small, Eq. ͑105͒ also describes the behavior of a general dynamical system close to a Hopf bifurcation. In fact, strictly speaking, Eq. ͑105͒ with K = 0 represents a dynamical system in the Hopf normal form, whenever the frequency does not depend on the amplitude ͑Crawford, 1991͒.
The oscillators in Eq. ͑105͒ are characterized by two degrees of freedom ͑e.g., amplitude and phase͒. Thus the behavior is richer than that of phase oscillators governed by the Kuramoto model. In particular, amplitude death and chaos may appear in some range of parameters characterizing the oscillator populations, such as coupling strength and natural frequency spread ͑Matthews and Strogatz, 1990͒. When amplitudes besides phases are allowed to vary in time, it may happen that all the oscillator amplitudes die out, that is ͉z j ͉ = 0 for every j. Such a behavior is referred to as "amplitude death" or "oscillator death," and this phenomenon occurs when ͑a͒ the coupling strength is of the same order as the attraction to the limit-cycle, and ͑b͒ the frequency spread is sufficiently large ͑Ermentrout, 1990͒. In biological systems, this may be responsible for rhythmical loss.
Apparently, Yamaguchi and Shimizu ͑1984͒ were the first to derive the model equation ͑105͒. They considered a system of weakly globally coupled Van der Pol oscillators that included white-noise sources. Bonilla et al. ͑1987͒ also analyzed Eq. ͑105͒ ͑with noise͒ for a population of identical oscillators. They studied the nonlinear Fokker-Planck equation for the one-oscillator probability density, and found a transition from incoherence to a time-periodic state via a supercritical Hopf bifurcation Bonilla et al. ͑1988͒.
A comprehensive analysis of Eq. ͑105͒ is presented by Matthews and Strogatz ͑1990͒, Mirollo and Strogatz ͑1990͒, and Matthews et al. ͑1991͒. The phenomenon of amplitude death and its stability regions have been analyzed by Mirollo and Strogatz ͑1990͒, in terms of the coupling strength and the spread of natural frequencies.
In the same paper it was also shown that an infinite system gives a good description of large finite systems. Matthews and Strogatz ͑1990͒ and Matthews et al. ͑1991͒ present a detailed study of all possible bifurcations occurring in Eq. ͑105͒, discussing locking, amplitude death, incoherence, and unsteady behavior ͑Hopf oscillations, large oscillations, quasiperiodicity, and chaos͒. Matthews et al. ͑1991͒ survey some of the previous work, namely, the contributions of Shiino and Francowicz ͑1989͒ and Ermentrout ͑1990͒. Shiino and Frankowicz, using a selfconsistent equation, established the existence of partially locked solutions as well as of amplitude death states, but gave no analytical results on stability. On the other hand, Ermentrout was probably the first to point out the amplitude death phenomenon, and to study its stability analytically for certain frequency distributions and values of the coupling.
Recently, the dynamical behavior of this model has been revisited by Monte and D'ovidio ͑2002͒, focusing on the time evolution of the order parameter. This was done by a suitable expansion, valid for any population size, but only for strong couplings and narrow frequency distributions. Even though the truncation of the hierarchy yielding the order parameter was arbitrary, several known features of such systems were recovered qualitatively. Due to the limitation on the frequency spread, however, no amplitude death could be detected.

D. Kuramoto model with inertia
Considering that the Kuramoto model is merely a phase oscillator model and ignores the dynamical behavior of the corresponding frequencies, one might attempt to generalize the original model and cast such features in it. The simplest way to accomplish such a task is to in-clude inertial effects in the model. This leads to the system of second-order stochastic differential equations: Here ⍀ i represents the natural frequency of the ith oscillator, while i = i is the instantaneous frequency. These ideas have been developed by Acebrón and Spigler ͑1998, 2000͒, on the basis of the bivariate nonlinear Fokker-Planck equation, corresponding to Eq. ͑107͒, in the limit of infinitely many oscillators. Such an equation is again a nonlinear integro-differential Fokker-Planck-type equation, which somehow generalizes the Fokker-plank description of the Kuramoto model. In an earlier paper, Ermentrout ͑1991͒ revisited the special problem of self-synchronization in populations of certain types of fireflies, as well as certain alterations in circadian cycles in mammals. The point was that the Kuramoto model yielded too fast an approach to the synchronized state in comparison to the experimental observations and, moreover, an infinite value is needed for the coupling parameter in order to obtain a full synchronization. Consequently an adaptive frequency model was proposed for N nonlinearly coupled secondorder differential equations for the phases. In such a formulation, the striking difference from the Kuramoto model is that the natural frequency of each oscillator may vary in time. Moreover, a nonlinear, in general nonsinusoidal coupling function was considered in a noiseless framework. On the other hand, Tanaka et al. ͑1997a, 1997b͒ considered the same problem described by Ermentrout, but under the mean-field coupling assumption and with a sinusoidal nonlinearity, and still without any noise term. They observed hysteretic phenomena ͑bistability between incoherence and partially synchronized state͒ due to a first-order phase transition, also referred to as a subcritical bifurcation. Such phenomena, however, may also occur within the Kuramoto model when it is governed by bimodal frequency distributions ͑Bonilla et al., 1992͒.
Recently, the interplay between inertia and time delay and their combined effect on synchronization within the Kuramoto model has been investigated both analytically and numerically by Hong et al. ͑2002͒. The main feature here was the emergence of spontaneous phase oscillations ͑without any external driving͒. Such oscillations were also found to suppress synchronization, whereas their frequency was shown to decrease when inertia and delay decreased. Hong et al. also obtained the phase diagram, which shows the stationary and oscillatory regimes as a function of delay, inertia, and coupling strength. It appears clearly that, for any finite value of inertia and time delay, the system undergoes a transition from a stationary to an oscillatory state when the coupling is sufficiently large. Further investigation, however, seems desirable because in the strong-coupling limit the phase model itself breaks down.
Another line of research that should be pursued concerns the emergence of stochastic resonance phenomena ͑see, for instance, Bulsara and Gammaitoni, 1996;Gammaitoni et al., 1998͒ in oscillator systems characterized by many degrees of freedom. An investigation in this direction has been started by Hong and Choi ͑2000͒, who are studying synchronization as well as noiseinduced resonance phenomena and synchronization in systems of globally coupled oscillators, each having finite inertia. Hong, Choi, Yi, et al. ͑1999͒ considered synchronization phenomena described by Eq. ͑107͒, in the absence of noise, but under an external periodic driving force, that is, in systems governed by where ⍀ is the frequency of external driving, and A i is the amplitude of the ith force term. They found that in a system with inertia, unlocked oscillators besides those locked to the external driving contribute to the collective synchronization, whereas in the absence of inertia, only those oscillators locked to the external driving contribute ͑Choi et al., 1994͒.
Going back to the model equations in Eq. ͑107͒, it was established by Acebrón and Spigler ͑1998͒ and Acebrón et al. ͑2000͒ that such a model is also capable of synchronizing simultaneously both phase and frequency, as can be observed within the Kuramoto model ͑Acebrón and Spigler, 2000͒. Proceeding formally, as in Sec. III, to derive a nonlinear Fokker-Planck equation for the oneoscillator probability density, in the limit of infinitely many oscillators, they obtained the new equation ͑Acebrón and Spigler, 1998͒ Here the order parameter is given by where g͑⍀͒ is the natural frequency distribution. As in the nonlinear Fokker-Planck equation for the Kuramoto model, initial value and 2-periodic boundary conditions with respect to should be prescribed. Moreover, a suitable decay of ͑ , , ⍀ , t͒ as → ± ϱ is required. The initial profile ͑ , , ⍀ ,0͒ should be normalized according to ͐ −ϱ +ϱ ͐ 0 2 ͑ , , ⍀ ,0͒dd = 1. Setting m = 0 in Eq. ͑107͒, one recovers the Kuramoto model. On the other hand, the effects of a small inertia on synchronization can be ascertained by analyzing Eq. ͑109͒ in the hyperbolic limit as m → 0. In this limit, the first two terms on the right-hand side of Eq. ͑109͒ are of the same order and dominate the others. Scaling the equation accordingly, the leading-order approximation to the probability density is a shifted Maxwellian function of the fre-quency times a slowly varying density function of time and of the other variables. When one uses the Chapman-Enskog procedure, the Kuramoto nonlinear Fokker-Planck equation is recovered as the zeroth-order approximation for the slowly varying density ͑Bonilla, 2000͒. More interestingly, the first-order correction to this equation contains m and therefore consistently includes the effects of a small inertia. Similar results were obtained earlier by Hong, Choi, Yoon, et al. ͑1999͒, who used an arbitrary closure assumption, thereby omitting relevant terms that the systematic Chapman-Enskog procedure provides ͑Bonilla, 2000͒.
The incoherent solution to Eq. ͑109͒ is a -independent stationary solution. According to the definition of the order parameter in Eq. ͑109͒, in this case the result r = 0 is obtained. Such a solution is given by Following a procedure similar to that adopted to study linear stability in the Kuramoto model ͑see Sec. III͒, one can perturb the solution with a small term, which can be assumed depending on time as e t , around the incoherent solution 0 . Technical difficulties, however, more serious than in the case of the Kuramoto model, beset the general procedure. First of all, the equation satisfied by the perturbative term is now itself a nonlinear partial integro-differential equation. After rather lengthy and cumbersome calculations, the equation for small values of noise. Note that in the limit of vanishing mass, the result K c =2͑D + ͒, occurring in the Kuramoto model, is recovered. It should also be noted that when the population is characterized by several frequencies, the inertia does play a role, making it harder to synchronize the oscillator populations. In fact, the critical coupling in Eq. ͑112͒ increases with m. Finally, for the bimodal frequency distribution g͑⍀͒ = ͓␦͑⍀ − ⍀ 0 ͒ + ␦͑⍀ + ⍀ 0 ͔͒ / 2, see the stability diagram in Fig. 8.
To sum up, increasing any of the quantities m, D, ͑inertia, noise, frequency spread͒ makes it more difficult to synchronize the oscillator population.
As where the n 's are given in terms of parabolic cylinder functions ͑or, equivalently, Hermite polynomials͒, and the c n 's obey a certain system of coupled partial differential equations ͑Acebrón et al., 2000͒. Retaining a suit-able number of such coefficients, one can analyze the type of bifurcations branching off the incoherence. It can be found that both supercritical and subcritical bifurcations may occur, depending on the value of the inertia ͑see Tanaka et al., 1997aTanaka et al., , 1997b, concerning the subcritical behavior͒. For instance, in the case of a Lorentzian frequency distribution, this behavior differs from that observed in the Kuramoto model, where only the supercritical bifurcation appears. Indeed, Fig. 9 shows that apart from the Kuramoto model, a transition to a subcritical bifurcation takes place for any given spread when the inertia is sufficiently large ͑the smaller , the larger m͒. Therefore inertial effects favor the subcritical character of bifurcations in the transition from incoherence to a partially synchronized state. Concerning the numerical simulations, a number of tools have been used. The problem is clearly harder than treating the Kuramoto model, due to the double dimensionality of the system. Extensive Brownian simulations have been conducted on the system in Eq. ͑107͒, and a Fourier-Hermite spectral method has been implemented to solve the new nonlinear Fokker-Planck equation ͑Acebrón et al., 2000͒. The latter approach was suggested by the fact that the distribution is required to be 2 periodic in and decay to zero as → ± ϱ, with a natural weight e −m 2 /4D .

Numerical treatment of stochastic differential equations
In this section we review only some general ideas about the numerical treatment of stochastic differential FIG. 8. Discrete bimodal frequency distribution: Stability diagram of incoherence in the parameter space ͑⍀ 0 , K͒ for different mass values and D = 1. From Acebrón et al., 2000. FIG. 9. Stationary solutions of a Lorentzian frequency distribution with spread , separating the regions in parameter space ͑m , ͒ in which the transition to the synchronized state is either subcritcal or supercritical: dotted line, the analytical solution obtained using c n , n = 0,1,2; solid line, the solution computed numerically. From  equations. The reader is referred to the review papers of Platen ͑1999͒ and Higham ͑2001͒, and the book of Kloeden and Platen ͑1999͒ for more details. An autonomous stochastic differential equation has the form dX͑t͒ = f"X͑t͒…dt + g"X͑t͒…dW͑t͒, ͑116͒ which should be considered as an abbreviation of the integral equation for 0 Ͻ t ഛ T. Here f is a given n-dimensional vector function, g is a given n ϫ n matrix valued function, and X 0 represents the initial condition. W͑t͒ denotes an m-dimensional vector whose entries are independent standard scalar Brownian motions ͑Wiener processes͒, and the second integral on the right-hand side of Eq. ͑117͒ is intended as a stochastic integral in the sense of Itô. Recall that the formal derivative, dW / dt, of the Brownian motion is the so-called white noise. The solution X͑t͒ to Eqs. ͑116͒ and ͑117͒ is an n-dimensional stochastic process, that is, an n-dimensional random variable vector for each t. When g = 0 the problem becomes deterministic and then Eq. ͑116͒ is an ordinary differential equation. The simplest numerical scheme to compute the solution to Eq. ͑116͒ is the natural generalization of the Euler method, which here takes on the form for j =1, ... ,M. Here ⌬W j = ͑W j − W j−1 ͒, and X j represents an approximation of X͑t j ͒, t j = j⌬t, and ⌬t = T / M. As is well known, the Brownian increments W j − W j−1 are random variables distributed according to a Gaussian distribution with zero mean and variance ⌬t. Therefore Eq. ͑116͒ can be solved by generating suitable sequences of random numbers. As for the error one should attach to approximate solutions to Eq. ͑116͒, there are two different types, depending on whether one is interested in obtaining the paths or the moments. In the literature, these are referred to as "strong" and "weak" approximations, respectively. In the strong schemes, one should estimate the error which provides a measure of the closeness of the paths at the end of the interval. In the weak schemes, one is interested only in computing moments or other functionals of the process X͑t͒. For instance, in the case of the first moment, the error is given by Note that estimating the latter is less demanding. It can be shown that the Euler scheme is of order 1 / 2, that is s = O͑⌬t 1/2 ͒ for a strong method, but of order 1 ͓ w = O͑⌬t͔͒ for a weak method. Higher-order methods, however, do exist; for instance, the stochastic generalization of the Heun method, Others are based on the stochastic Taylor formula, e.g., the Taylor formula of order 3 / 2, which, for the scalar case, reads ͑Kloeden and Platen, 1999͒ where ⌬Z j are random variables, distributed according to a Gaussian distribution with zero mean, variance ͑⌬t͒ 3 / 3, and correlation E͑⌬W j ⌬Z j ͒ = ͑⌬t͒ 2 / 2. Using the Heun method to compute paths ͑strong scheme͒, the error can be shown to be s = O͑⌬t͒, while for the corresponding weak scheme it is O͑⌬t 2 ͒. It can be proved that computing paths by such a formula produces the inherent error s = O͑⌬t 3/2 ͒, while computing moments gives an error of order O͑⌬t 3 ͒.

The Kuramoto model
The Kuramoto model for N globally coupled nonlinear oscillators is given by the system made up of stochastic differential equations ͑23͒ and ͑24͒. System ͑23͒ can be solved numerically by the methods described in the previous subsection ͑Sartoretto et . In Fig. 10, plots of the amplitude of the order parameter as a function of time, obtained using Euler, Heun, and Taylor formulas of order 3 / 2, are compared. The reference solution was provided by solving the nonlinear Fokker-Planck equation ͑26͒ with a spectral ͑thus extremely accurate͒ method ͑Acebrón, Lavrentiev, and Spigler, 2001͒. Experiments were conducted for increasing values of N, which showed that results stabilize by N = 500.
The mean-field coupling assumption implies that in the thermodynamic limit ͑N → ϱ͒, averaging each single path i ͑t͒ over all the noise realizations yields the same results as averaging over the entire oscillator population. In practice, the population size N is assumed to be sufficiently large. Indeed, numerical simulations have shown that populations of a few hundred oscillators, in many instances, behave qualitatively almost the same as a population consisting of infinitely many members; see Fig. 11. In Fig. 12, the discrepancy between the orderparameter amplitude r N , computed with N oscillators, and that obtained when N → ϱ, in the long-time regime, was plotted versus N. The slope of log 10 ⌬r N , where ⌬r N = ͉r N − r͉, versus log 10 N, shows that ⌬r N ϳ N −1/2 as N grows.

B. Simulating infinitely many oscillators
In this section, the focus is on the numerical solution of the nonlinear Fokker-Planck equation ͑26͒, which FIG. 11. Comparison of the results of numerical simulations obtained solving the nonlinear Fokker-Planck equation with a spectral method with 40 harmonics, and solving the system of stochastic differential equations for two different population sizes, N = 500 and N = 5000. Simulations were performed using a Taylor 3 / 2 scheme with a time step ⌬t = 0.1. Parameters are the same as in Fig. 10. represents the thermodynamic limit N → ϱ, either by finite differences or by a spectral method. Let us rewrite this equation as

Finite differences
Explicit finite differences are easy to implement and yield where i n is an approximation to ͑i⌬ , n⌬t , ͒, and initial and boundary data are prescribed. In Eq. ͑126͒, forward time differences and space-centered finite differences have been used. The parameter in Eq. ͑126͒ should be picked up from the support of the natural frequency distribution g͑͒. In general, the integral terms I͓ i n ͔ and J͓ i n ͔ involve a quadrature over the frequencies, which in fact requires a suitable numerical treatment whenever g͑͒ is a continuous function. For instance, when g͑͒ is a Lorentzian frequency distribution, the Gauss-Laguerre quadrature has been successfully used ͑Acebrón et al., 2000͒.
If we use an explicit scheme, the size of the time step ⌬t has to be kept sufficiently small for stability reasons: ⌬tD / ͑⌬͒ 2 Ͻ 1 / 2. Using an implicit finite-difference scheme, such as Crank-Nicholson's ͑Acebrón and Bonilla, 1998; Bonilla, Pérez-Vicente, and Spigler, 1998͒, we can, in principle, increase the time step ͑for a given ⌬͒. The trouble is that our problem is nonlinear, and therefore we need to implement an additional iterative procedure at each time step to find i n+1 . In practice, this reduces the time step to rather small values, as illustrated by numerical experiments. In Fig. 13 is compared with the results of the numerical solution of the nonlinear Fokker-Planck equation evaluated by either implicit finite differences or by a spectral method. Note that inserting Eq. ͑127͒ into Eq. ͑128͒ yields a nonlinear equation for r 0 and 0 , whose solution can be computed with very high accuracy by the Brent method ͑Press et al., 1992͒.

Spectral method
The periodicity of the distribution in the nonlinear Fokker-Planck equation suggests that a suitable numerical method for solving such a partial differential equation could be based on expanding in a Fourier series with respect to ͑Shinomoto and Sartoretto et al., 1998͒. Such a spectral method is known to be very efficient, in that convergence is expected to be achieved exponentially fast, based on the number of harmonics ͑Fornberg, 1996͒. Expanding in the nonlinear Fokker-Planck equation leads to a system of infinitely many ordinary differential equations: which should be truncated to a suitable number of harmonics n , for each frequency in the support of g͑͒.
Such a system can be solved numerically by one of numerous extremely sophisticated packages that are available. Acebrón, Perales, and Spigler ͑2001͒ used a variable-step Runge-Kutta-Felhberg routine ͑Press et al.,

1992͒.
When, as is often the case, one is only interested in evaluating the order parameter, a variant of the spectral method can be used ͑Acebrón and Bonilla, 1998͒. The main difference consists of centering the phase in the basis functions around the mean ͑unknown͒ phase . The probability density is first expanded in a Fourier series around the mean phase, Then, differentiating R n and ⌿ n with respect to t, using the nonlinear Fokker-Planck equation, integrating by parts, and exploiting the 2 periodicity of , one obtains the hierarchy Here, In many instances, for example, when g͑͒ is an even function, R 1 turns out to be also even as a function of , and thus the last equation in Eq. ͑133͒ yields d / dt =0. In order to compare this approach with that of Eq. ͑130͒, let us write Eq. ͑130͒ in terms of real quantities. Setting n = c n + is n , we obtain ṡ n = − n 2 Ds n + nc n + nK

͑136͒
Note that R n = 2͓c n cos͑n͒ − s n sin͑n͔͒, ͑137͒ ⌿ n = 2͓s n cos͑n͒ + c n sin͑n͔͒. ͑138͒ Comparing the two hierarchies, it is now clear that Eq. ͑133͒ is simpler, and, in fact, the corresponding numerical treatment was observed to be faster. The spectral method of Eq. ͑130͒ has been analyzed by Acebrón, Lavrentiev, and Spigler ͑2001͒, who obtained explicit bounds for the space derivatives. Such bounds play a role in estimating the error term in the Fourier-series expansion. The number N of harmonics required to achieve a given numerical error has been investigated as a function of the nonlinearity parameter K, the noise strength D, and ͐ −ϱ +ϱ g͑͒d. Indeed, the L 2 norm ͑with respect to ͒ of the error can be estimated as where C p is an estimate for the pth derivative of with respect to . In practice, C p depends on the initial data and on the parameter ͑K / D͒͐ −ϱ +ϱ g͑͒d. Since C p grows rapidly with such a parameter, larger nonlinearities as well as lower noise levels require a higher number of harmonics. In Fig. 14, the harmonics amplitude ͉ n ͉ is plotted as a function of n for several values of K.
The performance of the spectral method is illustrated in Fig. 13, where a comparison with the result of finitedifference simulations is shown. In the latter, a mesh-size ⌬ = 0.04 was used, while the former algorithm was run with N =2,4,8 harmonics. In Fig. 15, the global L 2 error as a function of N, with N =1/⌬, is plotted in logarithmic scales, which shows that the spectral method outperforms the finite-difference method. The CPU time needed to implement N = 12 harmonics in the spectral method was approximately 25 times smaller than when using finite differences with ⌬ = 0.04, ⌬t =10 −4 .
FIG. 14. The coefficients ͉ n ͑ϱ , ͉͒ corresponding to the stationary solution are shown for various coupling strength parameters K. The frequency distribution here is g͑͒ = ␦͑͒, and the diffusion constant D = 1. From Acebrón, Lavrentiev, and Spigler, 2001.

Tracking bifurcating solutions
A powerful numerical code exists that is capable of tracking bifurcating solutions in general systems of ordinary differential equations. Such a code, called AUTO ͑Doedel, 1971͒, is based on continuing a given solution up to and beyond several bifurcating branches. In this way, one is able to assess the stability properties of each branch.
The code AUTO was applied to the system of nonlinear deterministic ordinary differential equations obtained through the spectral method described in the previous section ͑Acebrón, Perales, and Spigler, 2001͒. This made possible the analysis of the bifurcations occurring in infinite populations of globally coupled oscillators described by the noisy Kuramoto model in Eq. ͑23͒. In particular, the stability issue could be examined for the synchronized stationary states.
In Fig. 16, a bifurcation diagram for a unimodal frequency distribution is shown. The reader is referred to Acebrón, Perales, and Spigler ͑2001͒ for other more elaborate pictures corresponding to multimodal frequency distributions. Global stability of ͑partially͒ synchronized stationary solutions was conjectured and investigated for a rather long time by Strogatz ͑2000͒, and numerical evidence for such a stability was provided.

C. The moments approach
Pérez-Vicente and Ritort ͑1997͒ proposed an alternative derivation of the nonlinear Fokker-Planck equation for the mean-field Kuramoto model to that reported in Appendix A. The advantage of this approach is threefold: ͑1͒ it provides an efficient way to solve the meanfield equations numerically in the limit of large N's, free of finite-size effects; ͑2͒ it provides a simple proof that the stationary solutions of the dynamics are not Gibbsian, and therefore they cannot be derived within a Hamiltonian formalism ͑Park and Choi, 1995͒; and ͑3͒ it can be used as a basis for including fluctuations beyond the mean field, in the framework of certain approximate closure schemes. The idea of such an approach relies on the rotational symmetry of the Kuramoto model, i → i +2. The most general set of observables, invariant under these local transformations, are where k , m are integers with m ജ 0. Note that H k m ͑t͒ = ͓H −k m ͑t͔͒ * . Under the dynamics described by Eq. ͑23͒, these observables do not fluctuate in the limit of large N's, i.e., they are both reproducible ͑independent of the realization of noise͒ and self-averaging ͑independent of the realization of the random set of i ͒. 1 After averaging over the noise, the set of observables in Eq. ͑140͒ closes the dynamics in the limit N → ϱ, The order parameter in Eq. ͑4͒ was introduced in Eq. ͑141͒ through the relation H 1 0 = r exp͑i͒. Equation ͑141͒ leads immediately to the nonlinear Fokker-Planck equation ͑26͒, in terms of the one-oscillator probability density 1 Of course the specific set of natural frequencies for the realization must be considered as typical ͑atypical realizations such as the nondisordered choice i = are excluded͒.
Solving the hierarchy in Eq. ͑141͒ requires specification of the initial conditions H k m ͑t =0͒, where the values of H 0 m = m = ͐ m g͑͒d are determined solely by the natural frequency distribution g͑͒. The equations can be solved using standard numerical integration schemes. For many purposes a Heun scheme suffices. Obviously, the number of moments to be included must be finite, and in general a few tens of them has been shown to suffice. The method is particularly suited to models in which m takes a finite number of values. For instance, in the bimodal case ͑Sec. III.C͒, m = 2, and the set of moments H k m splits into two subsets H k + and H k − corresponding to m even or odd, respectively. In this case, the number of equations that have to be integrated is considerably reduced. Moreover, to limit "boundary effects," periodic boundary conditions within the set of moments should be implemented.  Figure  19 depicts the parameter r͑t͒ corresponding to the oscillatory synchronized phase, after the transients have died out. In this figure, the results obtained by the method of moments are compared with those given by Brownian simulations.
The moment formalism allows us to prove that the stationary distribution is not Gibbsian. It was shown by Hemmen and Wreszinski ͑1993͒ that the Hamiltonian is a Lyapunov functional for finite N, whose minima describe completely phase-locked states at D = 0. Subsequently, resorting to standard arguments of equilibrium statistical mechanics, such a Lyapunov functional was employed to characterize stationary states of the Kuramoto model in the thermodynamic limit, at either zero or finite "temperature" D ͑see Park and Choi, 1995͒. Pérez-Vicente and Ritort ͑1997͒ objected that the latter result applied to Gibbsian states, but not to general stationary states. In Eq. ͑144͒, it was necessary to assume − Ͻ i Ͻ in order that H be a univalued function. As the Langevin equation leads to the nonlinear Fokker-Planck equation ͑see Appendix A͒, this may suggest that, indeed, H in Eq. ͑144͒ possesses all the properties of an energy functional. In fact, it is well known that the stationary distribution of the quantities in Eq. ͑145͒ will be Gibbsian only when the probability currents across the boundaries i =− , vanish. Only in this case, the stationary one-oscillator probability density is described by the equilibrium distribution obtained from Eq. ͑144͒. In the general case ͑see Appendix E͒, it is shown that the moments of the Gibbsian equilibrium distribution function are not stationary ͑Pérez-Vicente and Ritort, 1997͒. Moreover, the Hamiltonian in Eq. ͑144͒ is not a Lyapunov functional because it is not bounded from below, although it decays in time. Only at D = 0 and for finite N do the local minima of H correspond to globally synchronized solutions.
The moment approach has also been applied to other models, such as random tops ͑Sec. V.B͒ and synchronization models without disorder ͑Sec. V.A͒. Such an approach provides an alternative and equivalent way to describe certain dynamical features of synchronization models ͑Aonishi and Okada, 2002͒.

VII. APPLICATIONS
The outstanding adaptability and applicability of the Kuramoto model makes it suitable to be studied in many different contexts ranging from physics to chemistry. Here, we present some of the most relevant examples discussed in the literature over the past few years, but its potential use is certainly still growing, and new applications will appear in the future.

A. Neural networks
Perception is an old, fascinating, and unsolved neurophysiological problem that has attracted the attention of neuroscientists for decades. The basic and fundamental task of sensory processing is to integrate stimuli across multiple and separate receptive fields. Such a binding process is necessary to create a complete representation of a given object. Perhaps the most illustrative example is the visual cortex. Neurons that detect features are distributed over different areas of the visual cortex. These neurons process information from a restricted region of the visual field, and integrate it through a complex dynamical process that allow us to detect objects, separate them from the background, identify their characteristics, etc. All these tasks combine to give rise to cognition.
What kind of mechanisms might be responsible for the integration of distributed neuronal activity? There is a certain controversy around this point. It is difficult to find a good explanation using exclusively models that encode information only from the levels of activity of individual neurons ͑Softky and Koch, 1993͒. There are theories that suggest that the exact timing in a sequence of firing events is crucial for certain perception tasks ͑Abeles, 1991͒. On the other hand, it has been argued ͑Tovee and Rolls, 1992͒ that long-time oscillations are irrelevant for object recognition. Notice that oscillations and synchrony need to be distinguished. Cells can synchronize their responses without experiencing oscillatory discharges and, conversely, responses can be oscillatory without being synchronized. The important point is the correlation between the firing pattern of simultaneously recorded neurons. In this context, the idea that global properties of stimuli are identified through correlations in the temporal firing of different neurons has gained support from experiments in the primary visual cortex of the cat, showing oscillatory responses coherent over long distances and sensitive to global properties of stimuli ͑Eckhorn et Gray et al., 1989͒. Oscillatory response patterns reflect organized temporal structured activity that is often associated with synchronous firing.
This experimental evidence has motivated intense theoretical research, which looks for models capable of displaying stimulus-dependent synchronization in neuronal assemblies. It would be a formidable task to enumerate and discuss all of them because they mainly concern populations of integrate-and-fire neurons, which are beyond the scope of this review. Here we shall focus exclusively on studies in which the processing units are modeled as phase oscillators. At this point, let us mention that there are two different theoretical approaches. The first line of reasoning is very much biologically oriented, trying as the primary goal to reproduce, at least qualitatively, experimental results. The second one is more formal and is connected with associative memory models of attractor neural networks, which have been extensively studied in the past few years.

Biologically oriented models
Phase oscillators can be used as elementary units in models of neural information processing. This fact can be accepted as a natural consequence of the previous discussion, but it is desirable to look for sound arguments supporting this choice. In order to describe the emergence of oscillations in a single column of the visual cortex, Schuster and Wagner ͑1990a, 1990b͒ proposed a model of excitatory and inhibitory neurons. The main ingredients of the model are ͑a͒ Inhomogeneous spatial distribution of connections.
Neurons are densely connected at a local scale but sparsely connected on a larger scale. This combination of short-and long-range couplings gives rise to the existence of local clusters of activity.
͑b͒ Measurement of the activity of the neurons in terms of their mean firing rate. For sufficiently weak couplings, it is possible to describe the dynamics of the whole population in terms of a single mean-field equation. Except for a more complex form of the effective coupling among units, this equation is analogous to that governing the evolution of the phases for a limit-cycle oscillator.
In the case of the model proposed by Schuster and Wagner ͑1990a, 1990b͒, the coupling strength depends on the activity of the two coupled clusters, and this remarkable feature is the key to reproducing stimulus-dependent coherent oscillations. Sompolinsky and co-workers refined and extended the previous idea in a series of papers ͑Sompolinsky et al Grannan et al., 1993͒. They proposed a similar model with more elaborate interneuron coupling ͑synapses͒, in which many calculations can be performed analytically. They considered a Kuramoto model with effective coupling among oscillators, given by K ij ͑r , rЈ͒ = V͑r͒W͑r , rЈ͒V͑rЈ͒, where V denotes the average level of activity of the presynaptic and post-synaptic cell and W͑r , rЈ͒ denotes the specific architecture of the connections. To be more precise, W͑r , rЈ͒ has two terms, one describing strong interactions among neurons in the same cluster ͑with large overlapping receptive fields͒ and another describing the weak coupling among neurons in different clusters ͑without common receptive fields͒. In addition, neurons respond to a preferred orientation as they do in certain regions of the visual cortex. With all these ingredients, the model displays an extremely rich range of behaviors. The results of Sompolinsky et al. agree remarkably well with experiments, even for small dynamic details, to the point that they have made this research a reference for similar studies. It provides a mechanism for linking and segmenting stimuli that span multiple receptive fields through coherent activity of neurons.
The idea of reducing the complexity of neuron dynamics to simplified phase-oscillator equations has been very fruitful in different contexts. For instance, in order to model the interaction of the septohippocampal region and cortical columns, Kazanovich and Borisyuk ͑1994͒ analyzed a system of peripheral oscillators coupled to a central oscillator. Their goal was to understand the problem of focusing attention. Depending on the parameters of the model, they found different patterns of entrainment between groups of oscillators. Hansel et al. ͑1993a͒ studied phase locking in populations of Hodgkin-Huxley neurons interacting through weak excitatory couplings.
They used the phase-reduction technique to show that, in order to understand synchronization phenomena, one must analyze the effective interaction among oscillators. They found that, under certain conditions, a weak excitatory coupling leads to an effective inhibition among neurons due to a decrease in their firing rates. A little earlier, Abbott ͑1990͒ had carried out a similar program using piecewise linear Fitzhugh-Nagumo dynamics ͑Fitzhugh, 1961͒, instead of the Hodgkin-Huxley equations. As in previous models, a convenient reduction of the original dynamical evolution led to a simplified model in which neurons could be treated as phase oscillators.
More recently, Seliger et al. ͑2002͒ have discussed mechanisms of learning and plasticity in networks of phase oscillators. They studied the long-time properties of the system by assuming a Hebbian-like ͑Hebb, 1949͒ principle. Neurons with coherent mutual activity strengthen their synaptic connections ͑long-term poten-tiation͒, while in the opposite situation they weaken their connections ͑long-term depression͒. The slow dynamics associated with synaptic evolution gives rise to multistability, i.e., coexistence of multiple clusters of different sizes and frequencies. The work of Seliger et al. ͑2002͒ is essentially numerical. More elaborate theoretical work can be carried out provided neurons and couplings evolve on widely separated times scales, fast for neurons and slow for couplings. An example is the formalism developed by Jongen et al. ͑2001͒ for an XY spin-glass model. Further work in this direction is desirable.

Associative memory models
The field of neural networks experienced a remarkable advance during the last 15 years of the twentieth century. One of the main contributions was the seminal paper published by Hopfield ͑1982͒, which was the precursor of the current studies in computational neuroscience. He discovered that a spin system endowed with suitable couplings could exhibit an appealing collective behavior which mimics some basic functions of the brain. To be more precise, Hopfield considered a system of formal neurons, modeled as two-state units, representing the active and passive states of real neurons. The interactions among units ͑synaptic efficacies͒ followed Hebbian learning ͑Hebb, 1949͒. The resulting model exhibits an interesting phase diagram with paramagnetic, spin-glass, and ferromagnetic phases, the latter having effective associative-memory properties. The dynamics of the Hopfield model is a heat-bath Monte Carlo process, governed by the Hamiltonian where S represents the two states of the neuron ͑±1͒ and the couplings J ij represent the synaptic strength between pairs of neurons. The latter is given by a Hebbian prescription, where each set i represents a pattern to be learned.
Hopfield showed that the configurations i are attractors of the dynamics: when the initial state is in the basin of attraction of a given pattern ͑partial information about a given memory͒, the system evolves toward a final state having a large overlap with this learned configuration. This evolution mimics a typical associative memory process.
Typical questions concern the number of patterns that can be stored in the system ͑in other words, the amount of information that can be processed by the model͒, the size of the basin of attraction of the memories, the robustness of the patterns in response to noise, as well as other minor aspects. Usually, all these issues are tackled analytically by means of a standard method in the theory of spin glasses, namely, the replica approach ͑Mezard et al., 1987͒. It is not the goal of this review to discuss such a technique, but simply to give the main results. The number of patterns p that can be stored scales with the size of the system ͑number of neurons N͒ as p Ϸ 0.138N. Therefore the storage capacity defined as the ratio ␣ = p / N is 0.138. A detailed analysis of the whole phase diagram ␣ − T can be found in many textbooks ͑Amit, 1989;Hertz et al., 1991;Peretto, 1992͒. The conventional models of attractor neural networks ͑ANNs͒ characterize the activity of the neurons through binary values, corresponding to the active and nonactive state of each neuron. However, in order to reproduce synchronization between members of a population, it is convenient to introduce new variables which provide information about the degree of coherence in the time response of active neurons. This can be done by associating a phase with each element of the system, thereby modeling neurons as phase oscillators. A natural question is whether large populations of coupled oscillators can store information after a proper choice of the matrix J ij , as in conventional ANNs. Cook ͑1989͒ considered a static approach ͑no frequencies͒ in which each unit of the system is modeled as a q-state clock, and the Hebbian learning rule ͓Eq. ͑147͔͒ is used as a coupling. Since the J ij 's are symmetric, it is possible to define a Lyapunov functional whose minima coincide with the stationary states of Eq. ͑146͒. Cook solved the model by deriving mean-field equations in the replica-symmetric approximation. In the limit q → ϱ and zero temperature, she found that the stationary storage capacity of the network is ␣ c = 0.038, much smaller than the storage capacity of the Hopfield model ͑q =2͒, ␣ c = 0.138, or of the model with q = 3 where ␣ c = 0.22. Instead of fixing the couplings J ij , Gerl et al. ͑1992͒ undertook to estimate the volume occupied in the space of couplings J ij according to the couplings obeying the stability condition i ͚ j J ij j Ͼ ജ 0. By using a standard formalism due to Gardner ͑1988͒, they found that, in the optimal case and for a fixed stability , the storage capacity decreases as q increases, and that the information content per synapse grows when scales as q −1 . Although this final result seems promising, it has serious limitations given by the size of the network ͑since N Ͼ q͒, and also because the time required to reach the stationary state is proportional to q, as corroborated by Kohring ͑1993͒. Other models with similar features display the same type of behavior ͑Noest, 1988͒.
In the same context, the first model of phase oscillators having an intrinsic frequency distribution was studied by Arenas and Pérez-Vicente ͑1994b͒. These authors considered the standard Kuramoto model dynamics discussed in previous sections with coupling strengths given by where K is the intensity of the coupling. This form preserves the basic idea of Hebb's rule adapted to the symmetry of the problem. Now the state of the system, described by an N-dimensional vector whose ith component is the phase of the ith oscillator, is changing continuously in time. However, this is not a problem. If there is phase locking, the differences between the phases of different oscillators remain constant in time.
Then it is possible to store information as a difference of phases between pairs of oscillators, which justifies the choice of the learning rule given in Eq. ͑148͒. Thus, if the initial state is phase locked with one of the embedded patterns, then the final state will also have a macroscopic correlation with the same pattern, at least for small p.
In the limit of ␣ → 0, Arenas and Pérez-Vicente ͑1994b͒ showed that the degree of coherence between the stationary state and the best retrieved pattern is Here m is analogous to the order parameter r defined in previous sections, I n is a modified Bessel function of the first kind and of order n, ␤ = J /2D, and ͗ ͘ means taking the average over the frequency distribution. In contrast with conventional ANN models, which, in the limit ␣ → 0, have a positive overlap below the critical temperature, in this model phase locking can be destroyed whenever the distribution of frequencies is sufficiently broad. From a linear analysis of Eq. ͑149͒, it is straightforward to show that there is no synchronization when Similarly, Hong et al. ͑2001͒ found that, for zero temperature and for a Gaussian frequency distribution with spread , a retrieval state can only exist above a critical coupling value K c / . The quality of the retrieval depends on the number of stored patterns. Their numerical simulations show a rather complex time evolution towards the stationary state. An analysis of the short-time dynamics of this network was performed by Pérez-Vicente et al. ͑1996͒.
The situation is more complex for networks that store an infinite number of patterns ͑finite ␣͒. There have been some attempts to solve the retrieval problem through a standard replica-symmetry formalism ͑Park and Choi, 1995͒. However, as has already been discussed in Sec. VI.C, there is no Lyapunov functional in such a limit. Therefore other theoretical formalisms are required to elucidate the long-time collective properties of associative memory models of phase oscillators.
The long-time properties of networks of phase oscillators without a Lyapunov functional have been successfully determined by the so-called self-consistent signalto-noise analysis ͑see Fukai, 1992, 1993͒. To apply this formalism it is necessary to have information about the fixed-point equations describing the equilibrium states of the system. To be more precise, the method relies on the existence of Thouless-Anderson-Palmer-like equations ͑Thouless et al., 1997͒, which, derived in the context of spin glasses, have been very fruitful in more general scenarios. These equations relate the equilibrium time average of spins ͑phase oscillators͒ to the effective local field acting on them. In general, the effective and the time-averaged local fields are different, and their difference coincides with the Onsager reaction term, which can be computed by analyzing the free energy of the network. Once the Thouless-Anderson-Palmer equations are available, the local field splits into the "signal" and "noise" parts. Then self-consistent signal-to-noise analysis yields expressions for the order parameters of the problem, and the ͑extensive͒ number of stored patterns is determined as a result of this. This method has been applied by several authors. For the standard Hebbian coupling given in Eq. ͑147͒, Aonishi ͑1998͒, Uchiyama and Fujisaka ͑1999͒, Yamana et al. ͑1999͒, Yoshioka and Shiino ͑2000͒, and Aonishi et al. ͑2002͒ found peculiar associative memory properties for some special frequency distributions. For instance, in the absence of noise and for a discrete three-mode frequency distribution, Yoshioka and Shiino ͑2000͒ observed the existence of two different retrieval regions separated by a window where retrieval is impossible. In the -␣ phase diagram, for sufficiently low temperature, a nonmonotonic functional relationship is found. This remarkable behavior, which is a direct consequence of having nonidentical oscillators, is not observed in standard ANN models. In the appropriate limit, the results given by Cook ͑1989͒ and Arenas and Pérez-Vicente ͑1994b͒ are recovered.
There are complementary aspects of phase-oscillator networks that are worth studying. For instance, how many patterns can be stored in networks with diluted synapses ͑Aoyagi and Kitano, 1997;Kitano and Aoyagi, 1998͒ or with sparsely coded patterns ͑Aoyagi, 1995; Aoyagi and Nomura, 1999͒? So far, these analyses have been carried out for systems with identical oscillators, so that a static approach can be used. The main result is that a network of phase oscillators is more robust against dilution than the Hopfield model. On the other hand, for low levels of activity ͑sparse by coded pat-terns͒, the storage capacity diverges as 1 / a log a for a → 0, a being the level of activity, as in conventional ANN models. It is also an open problem to analyze the effect of intrinsic frequency distributions on the retrieval properties of these networks.

B. Josephson junctions and laser arrays
In addition to the extensive development of synchronization models, in particular the Kuramoto model, because of their intrinsic interest, in recent years, several applications of superconducting Josephson-junction arrays and of laser arrays to physics and technology have been explored in detail. The main purpose of these applications is to synchronize a large number of single elements in order to increase the effective output power. As is well known, the Kuramoto model provides perhaps the simplest way to describe such collective behavior, and it has been shown that the dynamics of Josephson-junction arrays ͑Wiesenfeld, 1996͒, and laser arrays ͑Kozireff et al., 2000Vladimirov et al., 2003͒ can be conveniently mapped onto it. Even though a few other physical applications have been found, to isotropic gas of oscillating neutrinos ͑Pantaleone, 1998͒, beam steering in phased arrays ͑Heath et al., nonlinear antenna technology ͑Meadows et al., 2002͒ in this subsection only the first two subjects mentioned above will be discussed at some length.

Josephson-junction arrays
Josephson junctions are superconducting devices capable of generating high-frequency voltage oscillations, typically in the range 10 10 -10 12 Hz ͑Josephson, 1964;Duzer and Turner, 1999͒. They are natural voltage-tofrequency transducers. Applications to both analog and digital electronics have been made, to realize amplifiers and mixers for submillimeter waves, to detect infrared signals from distant galaxies, and to use as very sensitive magnetometers in SQUIDs ͑superconducting quantum interference devices͒.
A large number N of interconnected Josephson junctions may be combined in such a way as to achieve a large output power. This occurs because the power is proportional to V 2 , which turns out to be proportional to N 2 , provided that all members oscillate in synchrony. Moreover, the frequency bandwidth decreases as N −2 ͑Duzer and Turner, 1999͒. Networks of Josephsonjunction arrays coupled in parallel lead to nearestneighbor ͑that is diffusive͒ coupling, and more precisely to sine-Gordon discrete equations with soliton solutions. On the other hand, Josephson-junction arrays connected in series through a load exhibit "all-to-all" ͑that is glo-bal͒ coupling ͑Wiesenfeld et al., 1996͒. Moreover, the lat-ter configuration can be transformed into a Kuramoto model. The model equations are for the junctions, and for the load circuit ͑see Fig. 20͒. Here is the difference between the phases of the wave functions associated with the two superconductors, Q is the charge, I j is the critical current of the jth junction, C j and r j are its capacitance and resistance, respectively, and I B is the "bias current," while R, L, and C denote resistance, inductance, and capacitance of the load circuit. Note that the load provides a global coupling, and when Q = 0 all the junctions work independently. In such a case, assuming for simplicity C j = 0 and I B Ͼ I j , the jth element will undergo oscillations at its natural frequency, By using an averaging procedure Swift et al. ͑1992͒, Wiesenfeld ͑1996͒, and Wiesenfeld et al. ͑1996, 1998͒ have shown that Eqs. ͑151͒ and ͑152͒ with C j = 0 can be mapped onto the Kuramoto model provided coupling and disorder are weak. The precise assumptions are listed in Appendix F. An important step in the derivation is using the "natural" angles j 's, defined by instead of the phases j . These angles are called natural because they describe uniform rotations in the absence of the coupling, Q = 0, while the j 's do not. Using the transformation to natural angles and the averaging procedure, Eq. ͑152͒ yields, to leading order, Here K and ␣ depend on the Josephson parameters r j , I j , I B , R, L, and C. Strictly speaking, the model equation ͑155͒ represents a generalization of the classical Kuramoto model because of the presence of the parameter ␣ 0 ͑see Sec. IV.B͒. This case was studied by Sakaguchi et al. ͑1987͒, who observed that ␣ Ͼ 0 gives rise to a higher value of the critical coupling ͑for a given frequency spread͒, and a lower number of oscillators are involved in the synchronization.
Heath and Wiesenfeld ͑2000͒ and Sakaguchi and Watanabe ͑2000͒ pointed out that the model equation of the Kuramoto type in Eq. ͑155͒ does not explain the operation of Josephson-junction arrays described by Eqs. ͑151͒ and ͑152͒ in certain regimes with C j 0. In fact, both physical and numerical experiments using the latter equations showed the existence of hysteretic phenomena ͑Sakaguchi and Watanabe, 2000͒ which were not found in the model equations ͑155͒. This discrepancy was explained by Heath and Wiesenfeld, who recognized that a more appropriate averaging procedure needed to be used. Doing that, a model formally similar to that of Eq. ͑155͒ was found, the essential difference being that now K and ␣ depend on the dynamical state of the system.
When the capacitances are assumed to be nonzero, say C j = C 0 0 ͑see Sakaguchi and Watanabe, 2000͒, hysteresis and multistability are found within the framework provided by Eqs. ͑151͒ and ͑152͒. Proceeding as above and considering a sufficiently small value of C 0 , again a Kuramoto-type phase model like that in Eq. ͑155͒ is found, where now K and ␣ depend also on C 0 . A result is that the nonzero capacitance facilitates the mutual synchronization. At this point, we should stress that the essential parameter distinguishing the two regimes of negligible and non-negligible capacitance is given by the McCumber parameter, ␤ =2eI c r 2 C / ប. Depending on the properties of the Josephson junction ͑that is, r, I c , and C͒, ␤ can vary over a very broad range, say 10 −6 -10 7 ͑see Duzer and Turner, 1999͒.
Filatrella et al. ͑2000͒ considered a model for a large number of Josephson junctions coupled to a cavity. These authors were able to reproduce the synchronization behavior reported in the experiments conducted by Barbara et al. ͑1999͒. Even though these experiments concerned two-dimensional arrays of Josephson junc- tions, they found qualitatively no differences with respect to the case of one-dimensional arrays ͑Filatrella et al., 2001͒. In these studies, the capacitances of the Josephson junctions are nonzero ͑underdamped oscilla-tors͒, and the bias current is taken to be such that each junction is in the hysteretic regime. Depending on the initial conditions, the junctions may work in each of two possible states, characterized by zero or nonzero voltage. When the junctions work in the latter case, the phases of the oscillators describing them vary with time, and the oscillators are called "active." By a suitable averaging method, Filatrella et al. ͑2001͒ derived the modified Kuramoto model, Here N a is the number of the "active oscillators," and K i takes on two possible values, a larger value if the ith oscillator is frequency-locked, and a smaller value if it is not. These authors also predicted that some overall hysteretic behavior should be observed under certain circumstances, a feature that could not be observed in the experiments conducted by Barbara et al. ͑1999͒. Daniels et al. ͑2003͒ showed that the resistively shunted junction ͑RSJ͒ equations describing a ladder array of Josephson junctions, which are overdamped ͑zero capacitance͒ and different ͑with disordered critical cur-rents͒, can be taken into a Kuramoto-type model. Such a model exhibits the usual sinusoidal coupling, but the coupling mechanism is of the nearest-neighbor type. This mapping was realized by a suitable averaging method upon which the fast dynamics of the RSJ equations can be integrated out, the slow scale being that over which the neighboring junctions synchronize. The effect of thermal noise on the junctions has also been considered, finding a good quantitative agreement between the RSJ model and the locally coupled Kuramoto model, when a noisy term was added. However, the noise term now appearing in the Kuramoto model was not obtained from the RSJ noisy model by the aforementioned transformation procedure, which raises some doubt about the general applicability of this result. Locally coupled Kuramoto models like this are analyzed in more detail in Sec. IV.A.
In closing, another application should be mentioned, in the general category of Josephson-junction arrays. A network of superconducting wires provides an additional example of an exact mean-field system ͑Park and Choi, 1997͒. Such a network consists of two sets of parallel superconducting wires, mutually coupled by Josephson junctions at each crossing point. It turns out that each wire interacts with only half of the others, namely, with those perpendicular to it, which results in a semiglobal coupling. It was found that the relevant model equations consist of two sets of coupled phase oscillator equations and, under special conditions, these equations reduce to the classical Kuramoto equation.

Laser arrays
The idea of synchronizing laser arrays of various kinds and analyzing their collective behavior is appealing both technologically and theoretically. On the one hand, entrainment of many lasers results in a large power output from a high number of low-power lasers. When the lasers are phase locked, a coherent high-power beam can be concentrated in a single-lobe far-field pattern ͑Vladimirov et al., 2003͒. On the other hand, this setting provides an additional example of a physical realization of the Kuramoto phase model. Actually, there are a few other generalizations of the Kuramoto model, that have been obtained by investigating systems of laser arrays, as will be mentioned below.
It has been observed in the literature that solid-state laser arrays and semiconductor laser arrays behave differently, due to the striking differences in their typical parameters. In fact, for solid-state lasers the linewidth enhancement factor is about zero, and the upper-level fluorescence lifetime, measured in units of photon lifetime, is about 10 6 , while for semiconductor lasers they are respectively about 5 and 10 3 . It follows that they exhibit quite different dynamical behavior.
Both local and global coupling among lasers have been considered over the years ͑Li and Erneux, 1993; Silber et al., 1993͒. As might be expected, globally coupled lasers act more efficiently when stationary synchronized ͑i.e., in-phase͒ states are wanted. Global coupling is usually obtained by means of an optical feedback given by an external mirror. A widely used model, capable of describing the dynamical behavior of coupled lasers, is given by the socalled Lang-Kobayashi equations ͑Lang and Kobayashi, 1980;Wang and Winful, 1988;Winful and Wang, 1988͒, which were obtained using Lamb's semiclassical laser theory. They are where E j denotes the jth laser dimensionless electricfield envelope, Z j the excess free-carrier density ͑also called gain of the jth laser͒, = N −1 ͚ k k is the average frequency in the array, ␦ j = j − is the frequency mismatch, ␣ is the linewidth enhancement factor, and is the feedback rate. Other parameters are t D , the external cavity roundtrip time ͑thus t D is the mean optical phase shift between emitter and feedback fields͒, and p Ϸ 10 −12 s and c Ϸ 10 −9 s, which are, respectively, the photon lifetime and the free-carrier lifetime. The parameter P j represents the excess pump above threshold ͑Vladimirov et al., 2003͒. For instance, assuming that the Z j 's are given, the first equation in Eq. ͑158͒ is reminiscent of the amplitude Kuramoto model ͑see Sec. V.C͒, but with time delay.
When the coupling is realized through an external mirror located at a Talbot distance of the order of 1 mm, the time delay t D can be neglected. When, instead, the array and the feedback mirror are much larger, the time delay is important and should be taken into account. This was done by Kozireff et al. ͑2000, 2001͒. Here, the synchronization of a semiconductor laser array with a wide linewidth ␣ was studied, writing out E j = ͉E j ͉e i⌽ j , and obtaining an asymptotic approximation for the ⌽ j 's from the third-order phase equation: Here the scaled time s = ⍀ R t ͑s D = ⍀ R t D ͒, ⍀ R = ͱ 2P / c p P = N −1 ͚ j P j , has been introduced. The parameters appearing in Eq. ͑159͒ are = ͑2P +1͒ ͱ p /2P c , ⍀ j = ͑P j / P −1͒ / , ⌬ j = ␦ j c / ͑1+2P͒, and K = ␣ c / ͑1+2P͒. and are the amplitude and phase of the complexvalued order parameter, More details can be found in Vladimirov et al. ͑2003͒. Note that Eq. ͑159͒ represents a generalization of the Kuramoto model equation, reducing to it, but with delay, when the second-and third-order derivatives are neglected ͑see Sec. IV.C͒. Neglecting only the third-order derivative, the Kuramoto model with inertia is recovered ͑see Sec. V.D͒. Kozireff et al. ͑2000, 2001͒ and Vladimirov et al. ͑2003͒ showed that the time delay induces phase synchronization and may be used to control it in all dynamical regimes.
Oliva and Strogatz ͑2001, where the interested reader can find additional references concerning this subject͒ investigated large arrays of globally coupled solid-state lasers with randomly distributed natural frequencies.
Based on previous work on lasers ͑Jiang and McCall, 1993;Braiman et al., 1995;Kourtchatov et al., 1995͒, as well as on general amplitude oscillator models, Oliva and Strogatz considered a simplified form of the Lang-Kobayashi equations, in which the gain dynamics were adiabatically eliminated. The resulting model equation was Note that this is indeed an amplitude model, similar to the Kuramoto amplitude model studied in Sec. V.C. Analytical results, such as stability boundaries of a number of dynamical regimes, have been obtained, showing the existence of such diverse states as incoherence, phase locking, and amplitude death ͑when the system stops lasing͒. C. Charge-density waves Strogatz et al. ͑1988, 1989͒ and Marcus et al. ͑1989͒ have proposed and analyzed a model related to the Kuramoto model for charge-density-wave ͑CDW͒ transport in quasi-one-dimensional metals and semiconductors. Charge-density-wave conduction and dynamics have been reviewed by Grüner andZettl ͑1985͒ andGrüner ͑1988͒. An important consideration is that charge-density waves are pinned by impurities for zero applied electric field and they move for fields above a critical level. This is called the depinning transition of the charge-density wave. Experiments show hysteresis in the transition between pinned and sliding charge-density waves ͑Grüner and Zettl, 1985͒. The model by Strogatz et al. ͑1988͒ is as follows: Here the i are phase oscillators, ␣ i are random pinning angles, E is the applied electric field, and K is the oscillator coupling constant. If the sine function in the coupling term is linearized, we obtain the mean-field Fukuyama-Lee-Fisher model of CDW dynamics ͑Grüner, 1988͒, which has a continuous ͑nonhysteretic͒ transition from pinned to sliding charge-density waves. The analysis of Eq. ͑162͒ is similar to that of the Kuramoto model. For small values of E / h and K / h, there are stationary states in which the oscillator phases are locked to the values ␣ i . Let us keep K / h fixed. As E / h increases, the stationary state loses stability and a sliding state with a nonvanishing electric current ͑proportional to ͚ i i / N͒ appears. This occurs at a certain depinning threshold for the electric field, which has been evaluated analytically by Strogatz et al. ͑1989͒. The bifurcations of the system ͑162͒, in terms of the parameters E / h and K / h, have been studied. The pinning transition turns out to be a subcritical ͑discontinuous͒ bifurcation, and therefore hysteretic phenomena and switching between pinned and sliding change-density waves occur. Both discontinuous and hysteretic behavior have been observed experimentally in certain CDW materials ͑Grüner and Zettl, 1985͒.
We should note that system ͑162͒ can be recast into a Kuramoto model affected by external field and disorder at the same time. In fact, setting ⌰ i = i − ␣ i in Eq. ͑162͒, we obtain where A ij = ␣ j − ␣ i . Clearly, h sin ⌰ i plays the role of an external field ͑see Sec. IV.D͒, and the A ij 's represent disorder ͑see Sec. IV.B͒. All oscillators have the same natural frequency E.
The existence of oscillations in chemical reactions is well known. The Belousov-Zhabotinsky reaction is paradigmatic, and models such as the Brusselator and the Oregonator have been invented to understand its properties. These models are described in terms of a few coupled nonlinear reaction-diffusion equations, which may have time-periodic solutions and give rise to different spatiotemporal patterns for appropriate parameter values.
The relation between chemical oscillators and phase oscillator models has been a matter of discussion for many years. In 1975, Marek and Stuchl ͑1975͒ coupled two Belousov-Zhabotinsky reaction systems with different parameters and hence observed different periodic oscillations. Each reaction occurred in a separate stirred tank reactor, and both reactors could exchange matter through the common perforated wall. They observed phase locking when the oscillation frequencies in the two reactors were close to each other. For large frequency differences, the solution of the coupled system exhibited long intervals of slow variation in the phase difference separated by rapid fluctuations over very short intervals. These observations were qualitatively explained by Neu ͑1979͒. He considered two identical planar limit-cycle oscillators that were linearly and weakly coupled. In addition, one oscillator had a small imperfection of the same order as the coupling. A singular perturbation analysis showed that the phase difference between the oscillators evolved according to an equation reminiscent of the Kuramoto model. Its analysis revealed phase locking and rhythm splitting ͑Neu, 1979͒. Neu ͑1980͒ extended this idea to populations of weakly coupled identical chemical oscillators. Adding small imperfections to the oscillators, the resulting model equations became Here d ij = q j − q i is the spatial displacement between the oscillators and i is a random imperfection parameter. When ⑀ = 0, there are N identical coupled oscillators having a stable T-periodic limit cycle given by Neu's analysis yields the following equation for the evolution of the phases in the slow time scale = ⑀t: where P is a T-periodic function determined from the basic limit-cycle solution ͑165͒ that satisfies PЈ͑0͒ =1, and ␤ is a parameter ͑Neu, 1979͒. When X = cos͑t + ͒, Y = sin͑t + ͒, and P͑͒ = sin , Eq. ͑166͒ is essentially the Kuramoto model. Neu ͑1980͒ analyzed synchronization in the case of identical oscillators ͑ i =0͒ for both meanfield and diffusive coupling in the limit of infinitely many oscillators. His method involves finding an evolution equation for the time integral of the order parameter. This equation can be solved in special cases and provides information on how the oscillators synchronize as time elapses ͑Neu, 1980͒. Thus the Kuramoto model describes weakly coupled chemical oscillators in a natural way, as already discussed by Kuramoto ͑1984͒ and Bar-Eli ͑1985͒ and worked out by other authors later. Following previous experiments on two-coupled stirred tank reactors, Yoshimoto et al. ͑1993͒ have tried to test frustration due to disorder in oscillation frequencies, in a system of three coupled chemical reactors. Previous studies in systems of two coupled stirred tank reactors have shown that, depending on the coupling flow rate, different synchronization modes can emerge spontaneously: an inphase mode, an antiphase mode, and a phase-death mode. Yoshimoto et al. ͑1993͒ interpreted their experiment involving three reactors by using the numerical solutions of Kuramoto-type equations for phase oscillators with asymmetric couplings. Their numerical solutions exhibited different combinations of the previous three modes, as well as new complex multistable modes whose features depended on the level of asymmetry in the interaction among oscillators. More recently, Kiss et al. ͑2002͒ have confirmed experimentally the existence of all these patterns and a number of other predictions of the Kuramoto model by using an array of 64 nickel electrodes in sulfuric acid.

VIII. CONCLUSIONS AND FUTURE WORK
In this paper, we have extensively reviewed the main features of the Kuramoto model, which has been very successful in understanding and explaining synchronization phenomena in large populations of phase oscillators. The simplicity of the Kuramoto model allows for a rigorous mathematical analysis, at least for the case of mean-field coupling. Still, the long-time behavior of the Kuramoto model is nontrivial, displaying a large variety of synchronization patterns. Furthermore, this model can be adapted so as to explain synchronization behavior in many different contexts.
Throughout this review we have mentioned a number of open lines of investigation deserving of special attention. Let us summarize some of them. For the mean-field Kuramoto model, the recent work by Balmforth and Sassi ͑2000͒ raises interesting questions to be tackled in the future. At zero noise strength, D = 0, a stability analysis of the partially synchronized phase and a rigorous description of the synchronization transition are needed. The necessary work in this direction is expected to be technically difficult.
Much more work is needed to understand synchronization in the Kuramoto model with nearest-neighbor coupling. Recent work by Zheng et al. ͑1998͒ on the 1D case has shown that phase slips and bursting phenomena occur for couplings just below threshold. This fact is not surprising. It is well known that spatially discrete equations with overdamped or underdamped dynamics are models for dislocations and other defects that can either move or be pinned. Among these models, we can cite the 1D Frenkel-Kontorova model ͑Braun and Kivshar, 1998;Bonilla, 2001, 2003b͒ or the chain oscillator models for 2D edge dislocations ͑Carpio and Bonilla, 2003a͒. The latter is exactly the 2D Kuramoto model with asymmetric nearest-neighbor coupling and zero frequency on a finite lattice. A constant external field acts on the boundary and is responsible for depinning the dislocations if it surpasses a critical value. For these models, there are analytical theories of the depinning transition, and the effects of weak disorder on the transition are also understood. Perhaps this methodology could be useful for understanding either synchronization or its failure in the nearest-neighbor Kuramoto model. It would also be interesting to analyze the entrainment properties of large populations of phase oscillators connected in a scale-free network or in a smallworld network. Perhaps new clustering properties or multistability phenomena will come out in a natural way.
Extensions of the Kuramoto model to new scenarios should be considered. More work is needed to understand the stability properties of synchronized phases in models with general periodic couplings or models with inertia and time delay.
We hope the extensive discussion of the Kuramoto model in this review helps in finding new applications of this model. For the applications discussed here, Josephson-junction arrays with nonzero capacitances need to be better understood, given their frequent occurrence in real systems. More careful singular perturbation methods should yield more general, yet tractable, models of the Kuramoto type. On the other hand, quantum noise in the form of spontaneous emission and shot noise are important for certain laser systems ͑Wieczorek and Lenstra, 2004͒. Examining the role of noise in such systems suggests a new research line. Concerning biological applications, it would be interesting to investigate in depth adaptive mechanisms that go beyond the standard learning rules discussed in this review.
Finally, models of phase oscillators different from those discussed in this review should be explored. In this context, recent works on the Winfree model ͑Ariarat- For the sake of completeness, a derivation of the nonlinear Fokker-Planck equation ͑26͒ satisfied by the oneoscillator probability density is presented. Such a derivation follows the ideas of Bonilla ͑1987͒, adapted to the case of the noisy Kuramoto model. It is somewhat technical but it has an important advantage over other derivations. In fact, it shows that, in the limit as N → ϱ, the p-oscillator probability density factorizes into the product of p one-oscillator densities for all t Ͼ 0, provided all oscillators were statistically independent at time t =0. This result of "propagation of molecular chaos" is usually assumed in other simpler derivations which close a hierarchy of equations for p-oscillator densities ͑Crawford and Davies, 1999͒.
To derive the nonlinear Fokker-Planck equation, let us first write down the path-integral representation of the N-oscillator probability density N ͑t , គ , គ ͒ corresponding to the system of stochastic equations in Eq.
͑23͒. N is equal to a product of ␦͓ i ͑t͒ − ⌰ i ͑t ; គ ͔͒, averaged over the joint Gaussian distribution for the white noise i ͑t͒, and over the initial distribution of the oscillators. ⌰ i ͑t ; គ ͒ are the solutions of Eqs. ͑23͒ for a given realization of the noises. We have

͑A1͒
This expression is then transformed using the fact that the delta functions are the Fourier transforms of unity, ͟ t ␦͑f j ͒ = ͐exp͓͐ 0 t i⌿ j f j dt͔D⌿ j ͑t͒, and the Gaussian average ͗exp͓͐ 0 t i⌿ j j dt͔͘ = exp͓−2D͐ 0 t ⌿ j 2 dt͔. It follows that after transforming the functional determinant as indicated by Bausch et al. ͑1976͒, Phythian ͑1977͒, and Dominicis and Peliti ͑1978͒; and averaging over the initial conditions. In Eq. ͑A2͒, គ = 1 , ... , N are the oscillator angles and ⌿ គ = ⌿ 1 , ... ,⌿ N are their conjugate variables in phase space. N is 2 periodic in each of its phase arguments j . Initial data have been assumed to be independent and identically distributed, N ͑0, គ , គ ͒ = ͟ j=1 N ͑ j , j ͒ ͑the molecular chaos assumption͒. In addition, the normalization constant in the path integral will be omitted. Since we are interested in one-oscillator averages for systems of infinitely many oscillators, we should study the one-oscillator probability density ͑ , , t͒ such that ͑ 1 , 1 ,t͒ = lim To analyze this function, note that the exponential in Eq. ͑A2͒ contains double sums such as and others of a similar form. Here ⌿ j = i⌿ j . Note that the squares of the sums in the previous formulas can be eliminated by using Gaussian path integrals,

͑A12͒
Here the exponential terms coincide with that in the integrand in Eq. ͑A8͒. After inserting Eq. ͑A11͒ in these exponentials, and substituting the result into Eq. ͑A12͒, the terms containing k ͑k =1,2,3͒ become −K͗⌿ sin ͘cos − K͗cos ͘⌿ sin . The denominator in Eq. ͑A12͒ can be set equal to 1, upon appropriately definining the path integral so that ͗1͘ϵ1. We then obtain ͗⌿ sin ͘ = ␦͗1͘ / ␦͗cos ͘ = 0. The other functions k can be determined similarly, and we obtain ͑,,t͒ = ͵͵ ͵

͑A15͒
This is the path-integral representation of the solution of the nonlinear Fokker-Planck equation satisfying ͑ , ,0͒ = ͑ , ͒. Thus the one-oscillator density satisfies the nonlinear Fokker-Planck equation in the limit as N → ϱ ͑Bonilla, 1987͒. The same method can be used to show propagation of molecular chaos: the p-oscillator probability density is given by p ͑ 1 , 1 , ... , p , p ,t͒ = lim provided that the oscillators are independent and identically distributed initially and that ͑N − p͒ → ϱ ͑Bonilla, 1987͒.

͑B6͒
respectively. A͑͒ and B͑͒ are slowly varying amplitudes to be determined later. Inserting Eqs. ͑B5͒ and ͑B6͒ into ͑B4͒ and using the nonresonance condition ͑45͒, one obtains the relation dA / d = F ͑0͒ , where F ͑0͒ is given by Eq. ͑48͒. Thus, to leading order, the method of multiple scales and the Chapman-Enskog method yield the same amplitude equation. However, for more complicated bifurcations, such as the degenerate transition described by Eq. ͑50͒, the method of multiple scales still yields dA / d = F ͑0͒ and a linear inhomogeneous equation for the amplitude B͑͒. The reason for these unphysical results is that the method of multiple scales is severely limited by the fact that all terms in the reduced equations provided by it turn out to be of the same order. Equation ͑50͒ can still be derived from these two equations by an ad hoc ansatz, as was done by Bonilla, Pérez-Vicente, and Spigler ͑1998͒ in the case of the tricritical point. Note that Eq. ͑39͒ implies that B͑͒ =0 when the Chapman-Enskog method is used. In this appendix, we evaluate the term F ͑2͒ ͑A , Ā ͒, needed to describe the transition from supercritical to subcritical bifurcations at the parameters 0 = D / ͱ 2, K c =3D. The solution of Eq. ͑42͒ is The additional equations in the hierarchy that are needed are L 4 = − K c ‫ץ‬ ͕ 3 Im e −i ͗e −iЈ , 1 ͖͘ − K 2 ‫ץ‬ 1 Im e −i ͗e −iЈ , 1 ͘ − F ͑0͒ ‫ץ‬ A 2 + c.c., ͑C2͒ L 5 = − K c ‫ץ‬ ͕ 4 Im e −i ͗e −iЈ , 1 ͖͘ − F ͑2͒ ‫ץ‬ A 1 + c.c.

APPENDIX D: CALCULATION OF THE BIFURCATION AT THE TRICRITICAL POINT
Inserting Eq. ͑59͒ into Eq. ͑26͒ leads to the modified hierarchy L 2 = − K c ‫ץ‬ ͕ 1 Im e −i ͗e −iЈ , 1 ͖͘ − A T e i D + i + c.c.,

͑D4͒
Equation ͑60͒ has not yet been used. Using it leads to a correction in Eqs. ͑D2͒ and ͑D3͒. The second term in Eq. ͑D1͒ has the same form as 1 , but it is nonresonant because ͗1,͑D + i͒ −2 ͘ = 0 at the tricritical point. The solution of Eq. ͑D1͒ is We now insert this solution into Eq. ͑D2͒, and use the ansatz ͑60͒ because ‫ץ‬ T 2 contains a factor A TT in a truly resonant term. Recall that it is at this point that the routine Chapman-Enskog ansatz in Eq. ͑38͒ fails to deliver a resonant term. Note that the ansatz in Eq. ͑60͒ adds the term F ͑1͒ e i / ͑D + i͒ 2 + c.c. to the right-hand side of Eq. ͑D3͒. The nonresonance condition for Eq. ͑D2͒ yields F ͑0͒ = D 2 ͑K 2 − 4 2 ͒A + 2 5 ͉A͉ 2 A. ͑D6͒ The solution of Eq. ͑D2͒ is Applying the nonresonance condition to the right-hand side of Eq. ͑D3͒, we obtain Inserting Eqs. ͑D6͒ and ͑D8͒ into Eq. ͑60͒ yields the sought amplitude equation ͑61͒.

APPENDIX E: STATIONARY SOLUTIONS OF THE KURAMOTO MODEL ARE NOT EQUILIBRIUM STATES
In this appendix, we show how the stationary solutions of the Kuramoto model are not equilibrium states for the model defined by the Hamiltonian ͑144͒. Following Pérez-Vicente and Ritort ͑1997͒, the equilibrium value of the moments in Eq. ͑140͒, E k m , corresponding to Eq. ͑144͒, can be easily evaluated: Stationarity of the equilibrium solution requires that v m vanish. In the simple case of a symmetric frequency distribution g͑͒ = g͑−͒, it can be proven that v 2m = 0 for every m. However, odd moments do not vanish. The temperature dependence of v 2m−1 can be analytically evaluated in the high-temperature limit ͓v 2m−1 = 2m + O͑␤ 3 ͔͒ as well as in the low-temperature limit ͓v 2m−1 = 2m + O͑T͔͒. It is worth noting that, even in the zerotemperature limit ͑where a stability analysis reveals that the synchronized solution is neutrally stable; see Sec. II.B͒, the absolute minima of H are not stationary configurations of the dynamics. This shows that the assumption that the stationary solutions at T = 0 are local minima of H in Eq. ͑144͒ does not hold. For a symmetric frequency distribution, where v 0 = 0, this result does not guarantee that the Boltzmann distribution P eq ͕͑ i ͖͒ ϰ exp͓−␤H͕͑ i ͖͔͒ ͑which depends on the whole set of moments H k m ͒ is also stationary. In this case, in fact, both H k 0 = h k and the synchronization parameter r are stationary.

APPENDIX F: DERIVATION OF THE KURAMOTO MODEL FOR AN ARRAY OF JOSEPHSON JUNCTIONS
By using an averaging procedure, Swift et al. ͑1992͒, Wiesenfeld ͑1996͒, and Wiesenfeld et al. ͑1996, 1998͒ have shown that Eqs. ͑151͒ and ͑152͒ with C j =0 ͑the resistively-shunted-junction or RSJ case͒ can be mapped onto the Kuramoto model, provided coupling and disorder are weak. The main steps of the procedure are as follows. First, we change from the angular variables j to the natural angles j , which rotate uniformly in the absence of coupling: Direct integration of this equation yields I B − I j sin j = ͑I B 2 − I j 2 ͒ / ͑I B − I j cos j ͒, which is equivalent to Eq. ͑154͒. In terms of the new angular variables, Eq. ͑152͒ becomes Second, we want to estimate the orders of magnitude of the terms in our equations. Equation ͑152͒ yields the characteristic time scale over which the angular variables change, t = ប / ͑2erĪ͒, where r and Ī are the mean values of the resistances and critical currents in the junctions. Similarly, t Q = ͑R + Nr͒C is the characteristic time scale in the resistance-inductance-capacitance ͑RLC͒ load circuit. By assuming that t / t Q Ӷ 1, we can ignore the first two terms on the left side of Eq. ͑152͒. If we ignore Q into Eq. ͑152͒, and insert the resulting j into ͑152͒, we obtain the following order of magnitude for the charge: Q = O͑NCĪr͒. Then Q = O͑NCĪr / t ͒ = O͑eNCĪ 2 r 2 / ប͒ in Eq. ͑152͒. Third, we assume that the disorder in r j and I j is weak. More precisely, we assume that the fluctuations in the critical current of the junc-tions are of the same order as Q = O͑eNCĪ 2 r 2 / ប͒, and that both are much smaller than Ī. Thus we assume Last, and according to our assumptions, we proceed to ignore terms Q and Q in Eqs. ͑151͒ and ͑152͒. This yields Q as a function of the angular variables. We then insert this result into Eq. ͑F2͒, which is averaged over the fast time scale t . We obtain the modified Kuramoto model ͑155͒.