Nonequilibrium firstorder phase transition in coupled oscillator systems with inertia and noise
Abstract
We study the dynamics of a system of coupled oscillators of distributed natural frequencies, by including the features of both thermal noise, parametrized by a temperature, and inertial terms, parametrized by a moment of inertia. For a general unimodal frequency distribution, we report here the complete phase diagram of the model in the space of dimensionless moment of inertia, temperature, and width of the frequency distribution. We demonstrate that the system undergoes a nonequilibrium firstorder phase transition from a synchronized phase at low parameter values to an incoherent phase at high values. We provide strong numerical evidence for the existence of both the synchronized and the incoherent phase, treating the latter analytically to obtain the corresponding linear stability threshold that bounds the firstorder transition point from below. In the limit of zero noise and inertia, when the dynamics reduces to the one of the Kuramoto model, we recover the associated known continuous transition. At finite noise and inertia but in the absence of natural frequencies, the dynamics becomes that of a wellstudied model of longrange interactions, the Hamiltonian meanfield model. Close to the firstorder phase transition, we show that the escape time out of metastable states scales exponentially with the number of oscillators, which we explain to be stemming from the longrange nature of the interaction between the oscillators.
pacs:
05.70.Fh, 05.70.Ln, 05.45.XtI Introduction
Collective synchronization refers to the remarkable phenomenon of a large population of coupled oscillators spontaneously synchronizing to oscillate at a common frequency, despite each constituent having a different natural frequency. This manybody cooperative effect is observed in many physical and biological systems, pervading length and time scales of several orders of magnitude. Some examples are metabolic synchrony in yeast cell suspensions Bier:2000 (), synchronized firings of cardiac pacemaker cells Winfree:1980 (), flashing in unison by groups of fireflies Buck:1988 (), voltage oscillations at a common frequency in an array of currentbiased Josephson junctions Wiesenfeld:1998 (), phase synchronization in electrical power distribution networks Filatrella:2008 (); Rohden:2012 (); Dorfler:2013 (), rhythmic applause Neda:2000 (), animal flocking behavior Ha:2010 (); see Ref. Strogatz:2003 () for a recent survey.
A paradigmatic model to study synchronization is the Kuramoto model comprising phaseonly oscillators of distributed natural frequencies that are globally coupled through the sine of their phase differences Kuramoto:1984 (); Strogatz:2000 (). Specifically, the system involves interacting oscillators . The th oscillator has natural frequency , and is characterized by its phase which is a periodic variable of period . The ’s have a common probability distribution given by . The phase evolves in time according to the equation
(1) 
where is the coupling constant, while the factor makes the model well behaved in the continuum limit .
In this work, we study a generalization of the dynamics Eq. (1) that includes inertial terms parametrized by a moment of inertia and stochastic noise parametrized by a temperature Acebron:1998 (); Hong:1999 (); Acebron:2000 (). Noise accounts for the temporal fluctuations of the natural frequencies Sakaguchi:1988 (), while inertia elevates the firstorder Kuramoto dynamics to secondorder Ermentrout:1991 (). For a general unimodal distribution of the natural frequencies, we report here the complete phase diagram of the model in the space of dimensionless moment of inertia, temperature, and width of the frequency distribution, showing that the system in the steady state may exist in either of two possible phases, namely, a synchronized phase and an unsynchronized or incoherent phase. We show that a nonequilibrium firstorder transition occurs from the synchronized phase at low parameter values to the incoherent phase at high values. While strong numerical evidence is provided to support the existence of both the synchronized and the incoherent phase, only the latter could be treated analytically to obtain the corresponding linear stability threshold that bounds the firstorder transition point from below. In proper limits of the dynamics, we recover the known continuous phase transitions in the Kuramoto model and in its noisy extension Sakaguchi:1988 (), and an equilibrium continuous transition in a related model of longrange interactions, the Hamiltonian meanfield model Inagaki:1993 ().
The Kuramoto model has been almost exclusively studied within the field of synchronization and nonlinear dynamical systems. On the other hand, there has been much recent activity within the community of statistical physicists to study nonequilibrium stationary states (NESSs) and develop a general framework akin to the one due to Boltzmann and Gibbs for equilibrium that allows analysis of nonequilibrium states on a general footing Mallick:2009 (). Unfortunately, there are few examples of NESSs for which one knows the probability measure of configurations exactly, so that the bulk of studies have relied on numerical simulations and approximate analysis Privman:1997 ().
Our work interprets the dynamics of the Kuramoto model to be of true nonequilibrium character. Moreover, quenched disorder in the form of natural frequencies of the oscillators provides a very rich setting to study the interplay of the nonequilibrium character of the dynamics with the disorder. In this rich backdrop, we are able to characterize the nature of the NESS and ascertain under quite general conditions the whole spectrum of phase transitions.
The paper is organized as follows. In the following section, we describe the model of interest and briefly review previous studies of the model. In Sec. III, we present the complete phase diagram of the model, providing numerical simulation results in support. In Sec. IV, we present an analytical treatment of the properties of the incoherent phase, based on the Kramers equation for the singleoscillator distribution. This is followed in Sec. V by a comparison of our analytical predictions with numerical simulations. The paper ends with conclusions. Some of the technical details are relegated to the two appendices.
Ii The model
We now give a precise definition of the generalized dynamics that we study in this paper. In addition to phase , we associate with the th oscillator another dynamical variable, namely, the angular velocity . With a Gaussian noise force and the natural frequency , the dynamics is Acebron:1998 (); Acebron:2000 ()
(2) 
where is the oscillator moment of inertia, is the friction constant, while is the synchronization order parameter:
(3) 
Here, we have
(4) 
with temperature in units of the Boltzmann constant. We consider a unimodal (that is, symmetric about mean , and decreases to zero with increasing ), and denote its width by . In the absence of inertia, the dynamics with the redefinition reduces at to that of the Kuramoto model Kuramoto:1984 (); Strogatz:2000 () and at to that of its extension studied by Sakaguchi in Ref. Sakaguchi:1988 ().
The dynamics Eq. (2) also describes motion of particles with an interaction on a unit circle, with and being respectively the angular coordinate, velocity and external torque. In the absence of ’s, Eq. (2) for is the microcanonical dynamics of the Hamiltonian meanfield model Inagaki:1993 (), a prototype of longrange interacting systems Campa:2009 (). In this case, the equations of motion are the Hamilton equations associated with the Hamiltonian
(5) 
with the momentum of the th particle. The dynamics of this system is microcanonical, conserving energy and total momentum. With no ’s, but , the dynamics of the resulting Brownian meanfield (BMF) model is canonical, mimicking the interaction of the HMF system with a heat bath Chavanis:2013 ().
The dynamics Eq. (2) is invariant under , and the effect of may be made explicit by replacing in the second equation with . We thus consider from now on the dynamics Eq. (2) with the substitution . In the resulting model, we take to have zero mean and unit width, without loss of generality.
For , using dimensionless variables
(6)  
(7)  
(8)  
(9)  
(10)  
(11) 
the dynamics becomes
where
(13) 
For , using dimensionless time , the dynamics becomes the overdamped motion
(14) 
where
(15) 
From now on, we will consider in place of dynamics Eq. (2) the reduced dynamics Eq. (LABEL:eomscaled) [that reduces for to the overdamped dynamics Eq. (14)] involving three dimensionless parameters, ; we will drop overbars for simplicity of notation. With (i.e. Acebron:1998 (),Acebron:2000 ()), the resulting BMF dynamics has an equilibrium stationary state Chavanis:2013 (). For other , the dynamics Eq. (LABEL:eomscaled) violates detailed balance due to the external driving by the set of torques , yielding a NESS. We demonstrate this in Appendix A.
Several stationary state aspects of the dynamics Eq. (LABEL:eomscaled) in the continuum limit are known. For the Kuramoto dynamics (), the system exhibits a continuous transition from a low synchronized [] to a high incoherent () phase across the critical point Kuramoto:1984 ()
(16) 
extending to , the point becomes a secondorder critical line on the plane, given, on using the results of Sakaguchi in Ref. Sakaguchi:1988 (), by solving
(17) 
For the BMF dynamics (), the synchronization transition is again continuous, occurring at the critical temperature given by Chavanis:2013 ()
(18) 
Although there have been some numerical studies of the full dynamics for nonzero Tanaka:1997 (); Acebron:1998 (); Acebron:2000 (), the complete synchronization phase diagram for a general unimodal has not been addressed before, a question we take up and answer in this paper. In the next section, we describe the complete phase diagram that emerges out of our analysis.
Iii Phase diagram
The complete phase diagram is shown schematically in Fig. 1(a), where the thick red secondorder critical lines stand for the continuous transitions mentioned above. For nonzero , the synchronization transition becomes firstorder, occurring across the shaded blue transition surface; this surface is bounded by the secondorder critical lines on the and planes, and by a firstorder transition line on the plane.
The phase diagram in Fig. 1(a) is a generalization of the one for typical fluids where a firstorder transition line ends in a critical point, while we have here a firstorder transition surface ending in critical lines. All transitions for are in NESS, and we interpret them to be of dynamical origin, accounted for by stability considerations of stationary solutions of equations describing evolution of phasespace distribution. Showing that the phases extremize a freeenergylike quantity (e.g., a large deviation functional Touchette:2009 ()) in NESS is a daunting task in the absence of a general framework akin to that for equilibrium Derrida:2005 (). For , the different phases actually minimize the equilibrium free energy Campa:2009 ().
To confirm the firstorder transition, we performed body simulations involving integrations of Eq. (LABEL:eomscaled) for a representative , i.e., a Gaussian. Details of the simulation procedure are given in the Appendix B. For given and and an initial state with oscillators at and angular velocities sampled from a Gaussian distribution with zero mean and width , we let the system equilibrate at . We then tune adiabatically to high values and back in a cycle. Figure 2 shows the behavior of the synchronization order parameter for several ’s at a fixed less than the BMF transition point , illustrating sharp jumps and hysteresis behavior expected of a firstorder transition. With decrease of , the jump in becomes less sharp and the hysteresis loop area decreases, both consistent with the transition becoming secondorderlike as decreases, see Fig. 1. For , Fig. 2 shows and , the stability thresholds for the incoherent and the synchronized phase, respectively; the phase transition point lies in between the two thresholds (see Fig. 1(d)). Figure 2 shows that the thresholds decrease and approach zero with the increase of ; it also suggests, together with Fig. 3, that and coincide both on the second order critical lines and as at a fixed .
For given and and between and , versus time in the stationary state shows bistability, with the system switching back and forth between incoherent and synchronized states [Fig. 4(a)]. To have nottoolarge switching times, these simulations have been performed with a relatively small number of oscillators, , causing large fluctuations in the order parameter . Therefore, in Fig. 4(a) the synchronized and the unsynchronized state are characterized by values of fluctuating above and below , respectively; however, this does allow for a clear visualization of the switches. The distribution in Fig. 4(b) is bimodal with a peak around or as varies between and , consistent with the transition being firstorder. Indeed, a firstorder transition point is characterized by two equally likely values of the order parameter, while at a secondorder phase transition point, the order parameter has its value equal to zero Goldenfeld:1992 ().
Iv Analytical treatment
We now turn to an analytical treatment of the firstorder transition. In the continuum limit , the dynamics Eq. (LABEL:eomscaled) is described by the singleoscillator distribution which gives at time and for each the fraction of oscillators with phase and angular velocity . The distribution is periodic in , and obeys the normalization , while evolving following the Kramers equation Acebron:2000 ()
(19) 
where . We now give the derivation of Eq. (19), followed by a discussion of its stationary solution corresponding to the incoherent phase.
iv.1 The Kramers equation for the singleoscillator distribution: Incoherent stationary state
Here, we start with deriving the BogoliubovBornGreenKirkwoodYvon (BBGKY) hierarchy equations for the dynamics Eq. (LABEL:eomscaled) for any number of oscillators , and, from there, by considering the limit , obtain the Kramers Eq. (19). For simplicity, we first discuss the derivation of the BBGKY equations for the case of a bimodal , and then generalize it to a general .
Consider a given realization of , in which there are oscillators with frequencies , and oscillators with frequencies , where . We then define the oscillator distribution function as the probability density at time to observe the system around the values . In the following, we use the shorthand notations and . Note that satisfies the normalization . We assume that

is symmetric with respect to permutations of dynamical variables within the group of oscillators with the same frequency, and

, together with the derivatives , vanish on the boundaries of the phase space.
The evolution of follows the FokkerPlanck equation which may be straightforwardly derived from the equations of motion Eq. (LABEL:eomscaled):
(20) 
where we have defined the column vector whose first entries equal and the following entries equal , and where the superscript denotes matrix transpose operation: .
To proceed, we follow standard procedure Huang:1987 (), and define the reduced distribution function , with and , as
(21) 
Note that the following normalizations hold for the singleoscillator distribution functions: , and .
Using Eq. (20) in Eq. (21) and simplifying, we get the BBGKY hierarchy equations for oscillators with frequencies as
(22) 
and similar equations for for oscillators of frequencies . The first equations of the hierarchy are
(23) 
and
(24) 
In the limit of large , we can write
(25) 
In order to generalize Eqs. (23) and (24) to the case of a continuous , we denote for this case the singleoscillator distribution function as . The first equation of the hierarchy is then
In the continuum limit , we may neglect twooscillator correlations and approximate as
(27) 
so that Eq. (LABEL:eq:finaleqn) reduces to the Kramers Eq. (19).
The stationary solutions of Eq. (19) are obtained by setting the left hand side to zero. For , the stationary solution is
(28) 
that corresponds to canonical equilibrium, with determined selfconsistently Chavanis:2013 (). For , the incoherent stationary state is Acebron:2000 ()
(29) 
The existence of the synchronized stationary state is borne out by our simulation results discussed above (see Figs. 2,3, and 4), although its analytical form is not known.
iv.2 Linear stability analysis of the incoherent state
Let us now discuss the linear stability analysis of the incoherent state Eq. (29), pursued in Ref. Acebron:2000 () by linearizing Eq. (19) about the state by expanding as
(30) 
with . The solution of the linearized equation yields that satisfies Acebron:2000 ()
(31) 
The above equation contains valuable information about the range of values of the parameters for which the incoherent state is stable, and consequently, about the transition from the incoherent to synchronized phase. This warrants a detailed analysis of Eq. (31) for a general unimodal . The analysis for Lorentzian in Ref. Acebron:2000 () left untouched the crucial issue of the synchronization transition.
We rewrite Eq. (31) as
(32) 
where is unimodal. The incoherent state is unstable if there is a with a positive real part that satisfies the above eigenvalue equation. We will now prove that, depending on the values of the parameters appearing in the above equation, there can be at most one such that can be only real. In addition, for the case of a Gaussian explicitly used in simulations reported in this paper, we obtain the general shape of the surface in the space that defines the instability region of the incoherent state.
Considering and strictly positive, we multiply for convenience the numerator and denominator of Eq. (IV.2) by to obtain
(33) 
Let us first look for pure imaginary solutions of this equation. Separating into real and imaginary parts, we have
(34)  
(35) 
In the second equation above, we make the change of variables , and exploit the parity in of the sum, to obtain
(36) 
It can be shown that the sum on the righthand side is positive definite for any finite . Furthermore, for our class of distribution functions, one may see that the term in square brackets is positive (respectively, negative) definite for (respectively, for ). As a consequence, the last equation is never satisfied for and finite, and therefore, the eigenvalue equation does not admit pure imaginary solutions [the proof holds also for the particular case , as may be checked]. We also conclude that there can be at most one solution with positive real part. In fact, if in the complex plane, we perform the loop depicted in Fig. 5(a) (where it is meant that the points and represent , respectively, and the radius of the arc extends to ), then, in the complex plane, we obtain, due to the sign properties of just described, the loop qualitatively represented in Fig. 5(b). The point in Fig. 5(b) is obtained for , in Fig. 5(a), for values at points and and in the whole of the arc extending to infinity. The position of the point in the complex plane is determined by the value of , which is given by
(37) 
From the wellknown theorem of complex analysis on the number of roots of a function in a given domain of the complex plane Smirnov:1964 (), we therefore obtain that for , there is one and only one solution of the eigenvalue equation with positive real part; on the other hand, for , there is no such solution. When the single solution with positive real part exists, it is necessarily real, since a complex solution would imply the existence of its complex conjugate. The value of is readily seen to be equal to for . For positive , the value will depend on the particular form of the distribution function . However, it is possible to prove that the value is always smaller than ; this is consistent with the physically reasonable fact that if the incoherent state is stable for , which happens for , it is all the more stable for .
The surface delimiting the region of instability in the phase space is implicitly defined by Eq. (IV.2) [i.e. ], which, in principle, can be solved to obtain the threshold value of (denoted by ) as a function of : . On physical grounds, we expect that the latter is a single valued function, and that for any given value of , it is a decreasing function of for , reaching for . We are able to prove analytically these facts for the class of unimodal distribution functions considered in this work that includes the Gaussian case. However, we can prove in general for any that tends to for . This is done using the integral representation
(38) 
For and , one may see that the term within the integral in the last equation tends to . We thus obtain by examining Eq. (IV.2) that . Combined with the fact that , this shows that .
Let us now turn to the Gaussian case, . Denoting with a subscript in this case, and using Eq. (IV.2), we have
(39) 
The integral in can be easily performed. Making the change of variable , we arrive at the following equation:
(40) 
The equation defines implicitly the function . We can show that this is a singlevalued function with the properties and . We show this by explicitly computing the partial derivatives of with respect to and , and by evaluating the behavior with respect to changes in by adopting a suitable strategy. We begin by computing the derivative with respect to . From Eq. (IV.2), we readily obtain
(41) 
which is clearly negative. Second, the derivative with respect to gives
This derivative is negative, since is positive for . From the implicit function theorems, we then derive that . The study of the behavior with respect to a change in is a bit more complicated. Since we are considering , we multiply Eq. (IV.2) by to obtain
(43) 
Let us consider the integral on the righthand side
(44) 
Since is negative for , we conclude that the derivative of this expression is negative, while its second derivative is positive. Then the righthand side of Eq. (IV.2) can be zero, for , for at most one value of . Furthermore, since for fixed and the value of decreases if increases, the value for which decreases for increasing at fixed . This concludes the proof. Furthermore, for what we have seen before, and for .
From the above analysis, it should be clear that the proof is not restricted to the Gaussian case, but would work exactly in the same way for any such that
(45) 
is positive for any . However, on physical grounds, we are led to assume that the same conclusions hold for any unimodal .
On the basis of our analysis, it follows that at the point of neutral stability, one has , which when substituted in Eq. (31) gives to be satisfying
(46) 
In the space, Eq. (46) defines the stability surface . There will similarly be the stability surface (see Fig. 1(d) which shows the two surfaces as obtained in body simulations for for a Gaussian ). The two surfaces coincide on the critical lines on the and planes where the transition becomes continuous; outside these planes, the surfaces enclose the firstorder transition surface i.e., . We now show by taking limits that the surface meets the critical lines on the and planes, and also obtain its intersection with the plane. On considering at a fixed , only the term in the sum in Eq. (46) contributes, giving