Nonequilibrium first-order phase transition in coupled oscillator systems with inertia and noise
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 first-order 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 first-order 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 well-studied model of long-range interactions, the Hamiltonian mean-field model. Close to the first-order 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 long-range nature of the interaction between the oscillators.
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 many-body 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 , synchronized firings of cardiac pacemaker cells , flashing in unison by groups of fireflies , voltage oscillations at a common frequency in an array of current-biased Josephson junctions , phase synchronization in electrical power distribution networks , rhythmic applause , animal flocking behavior ; see Ref.  for a recent survey.
A paradigmatic model to study synchronization is the Kuramoto model comprising phase-only oscillators of distributed natural frequencies that are globally coupled through the sine of their phase differences . 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
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. (Equation 1) that includes inertial terms parametrized by a moment of inertia and stochastic noise parametrized by a temperature . Noise accounts for the temporal fluctuations of the natural frequencies , while inertia elevates the first-order Kuramoto dynamics to second-order . 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 first-order 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 first-order 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 , and an equilibrium continuous transition in a related model of long-range interactions, the Hamiltonian mean-field model .
The Kuramoto model has been almost exclusively studied within the field of synchronization and non-linear 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 . 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 .
Our work interprets the dynamics of the Kuramoto model to be of true non-equilibrium 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 Section 3, we present the complete phase diagram of the model, providing numerical simulation results in support. In Section 4, we present an analytical treatment of the properties of the incoherent phase, based on the Kramers equation for the single-oscillator distribution. This is followed in Section 5 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.
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 
where is the oscillator moment of inertia, is the friction constant, while is the synchronization order parameter:
Here, we have
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  and at to that of its extension studied by Sakaguchi in Ref. .
The dynamics Eq. (Equation 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. (Equation 2) for is the microcanonical dynamics of the Hamiltonian mean-field model , a prototype of long-range interacting systems . In this case, the equations of motion are the Hamilton equations associated with the Hamiltonian
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 mean-field (BMF) model is canonical, mimicking the interaction of the HMF system with a heat bath .
The dynamics Eq. (Equation 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. (Equation 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
the dynamics becomes
For , using dimensionless time , the dynamics becomes the overdamped motion
From now on, we will consider in place of dynamics Eq. (Equation 2) the reduced dynamics Eq. (Equation 4) [that reduces for to the overdamped dynamics Eq. (Equation 5)] involving three dimensionless parameters, ; we will drop overbars for simplicity of notation. With (i.e. ,), the resulting BMF dynamics has an equilibrium stationary state . For other , the dynamics Eq. (Equation 4) 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. (Equation 4) 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 
extending to , the point becomes a second-order critical line on the -plane, given, on using the results of Sakaguchi in Ref. , by solving
For the BMF dynamics (), the synchronization transition is again continuous, occurring at the critical temperature given by 
Although there have been some numerical studies of the full dynamics for non-zero , 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.
The complete phase diagram is shown schematically in Figure 1(a), where the thick red second-order critical lines stand for the continuous transitions mentioned above. For non-zero , the synchronization transition becomes first-order, occurring across the shaded blue transition surface; this surface is bounded by the second-order critical lines on the and planes, and by a first-order transition line on the -plane.
The phase diagram in Figure 1(a) is a generalization of the one for typical fluids where a first-order transition line ends in a critical point, while we have here a first-order 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 phase-space distribution. Showing that the phases extremize a free-energy-like quantity (e.g., a large deviation functional ) in NESS is a daunting task in the absence of a general framework akin to that for equilibrium . For , the different phases actually minimize the equilibrium free energy .
To confirm the first-order transition, we performed -body simulations involving integrations of Eq. (Equation 4) for a representative , i.e., a Gaussian. Details of the simulation procedure are given in the Appendix 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 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 first-order transition. With decrease of , the jump in becomes less sharp and the hysteresis loop area decreases, both consistent with the transition becoming second-order-like as decreases, see Fig. Figure 1. For , Figure 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 Figure 1(d)). Figure 2 shows that the thresholds decrease and approach zero with the increase of ; it also suggests, together with Figure 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 [Figure 4(a)]. To have not-too-large switching times, these simulations have been performed with a relatively small number of oscillators, , causing large fluctuations in the order parameter . Therefore, in Figure 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 Figure 4(b) is bimodal with a peak around or as varies between and , consistent with the transition being first-order. Indeed, a first-order transition point is characterized by two equally likely values of the order parameter, while at a second-order phase transition point, the order parameter has its value equal to zero .
We now turn to an analytical treatment of the first-order transition. In the continuum limit , the dynamics Eq. (Equation 4) is described by the single-oscillator 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 
where . We now give the derivation of Eq. (Equation 9), followed by a discussion of its stationary solution corresponding to the incoherent phase.
4.1The Kramers equation for the single-oscillator distribution: Incoherent stationary state
Here, we start with deriving the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy equations for the dynamics Eq. (Equation 4) for any number of oscillators , and, from there, by considering the limit , obtain the Kramers Eq. (Equation 9). 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 Fokker-Planck equation which may be straightforwardly derived from the equations of motion Eq. (Equation 4):
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 , and define the reduced distribution function , with and , as
Note that the following normalizations hold for the single-oscillator distribution functions: , and .
and similar equations for for oscillators of frequencies . The first equations of the hierarchy are
In the limit of large , we can write
In order to generalize Eqs. (Equation 12) and (Equation 13) to the case of a continuous , we denote for this case the single-oscillator distribution function as . The first equation of the hierarchy is then
In the continuum limit , we may neglect two-oscillator correlations and approximate as
The stationary solutions of Eq. (Equation 9) are obtained by setting the left hand side to zero. For , the stationary solution is
4.2Linear stability analysis of the incoherent state
with . The solution of the linearized equation yields that satisfies 
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. (Equation 17) for a general unimodal . The analysis for Lorentzian in Ref.  left untouched the crucial issue of the synchronization transition.
We rewrite Eq. (Equation 17) as
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. (Equation 18) by to obtain
Let us first look for pure imaginary solutions of this equation. Separating into real and imaginary parts, we have
In the second equation above, we make the change of variables , and exploit the parity in of the sum, to obtain
It can be shown that the sum on the right-hand 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 Figure 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 Figure 5(b). The point in Figure 5(b) is obtained for , in Figure 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
From the well-known theorem of complex analysis on the number of roots of a function in a given domain of the complex plane , 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. (Equation 22) [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
For and , one may see that the term within the integral in the last equation tends to . We thus obtain by examining Eq. (Equation 22) 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. (Equation 23), we have
The integral in can be easily performed. Making the change of variable , we arrive at the following equation:
The equation defines implicitly the function . We can show that this is a single-valued 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. (Equation 25), we readily obtain
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. (Equation 25) by to obtain
Let us consider the integral on the right-hand side
Since is negative for , we conclude that the derivative of this expression is negative, while its second derivative is positive. Then the right-hand side of Eq. (Equation 28) 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
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. (Equation 17) gives to be satisfying
In the space, Eq. (Equation 31) defines the stability surface . There will similarly be the stability surface (see Figure 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 first-order transition surface i.e., . We now show by taking limits that the surface meets the critical lines on the