Nonequilibrium first-order phase transition in coupled oscillator systems with inertia and noise

Nonequilibrium first-order phase transition in coupled oscillator systems with inertia and noise

Shamik Gupta, Alessandro Campa and Stefano Ruffo Laboratoire de Physique Théorique et Modèles Statistiques (UMR CNRS 8626), Université Paris-Sud, Orsay, France
Health and Technology Department, Istituto Superiore di Sanità, and INFN Sezione Roma1, Gruppo Collegato Sanità, Roma, Italy
Department of Physics and Astronomy and CSDC, University of Florence, CNISM and INFN, via G. Sansone, 1 50019 Sesto Fiorentino, Italy

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.

05.70.Fh, 05.70.Ln, 05.45.Xt

I 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 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 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 current-biased 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 phase-only 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


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 first-order Kuramoto dynamics to second-order 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 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 Sakaguchi:1988 (), and an equilibrium continuous transition in a related model of long-range interactions, the Hamiltonian mean-field model Inagaki:1993 ().

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 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 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 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 single-oscillator 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 ()


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 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 mean-field model Inagaki:1993 (), a prototype of long-range interacting systems Campa:2009 (). 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 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


the dynamics becomes



For , using dimensionless time , the dynamics becomes the overdamped motion




From now on, we will consider in place of dynamics Eq. (2) the reduced dynamics Eq. (LABEL:eom-scaled) [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:eom-scaled) 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:eom-scaled) 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 ()


extending to , the point becomes a second-order critical line on the -plane, given, on using the results of Sakaguchi in Ref. Sakaguchi:1988 (), by solving


For the BMF dynamics (), the synchronization transition is again continuous, occurring at the critical temperature given by Chavanis:2013 ()


Although there have been some numerical studies of the full dynamics for non-zero 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.

Figure 1: (Color online) (a) Schematic phase diagram of model Eq. (LABEL:eom-scaled) in terms of dimensionless moment of inertia , temperature , and width of the frequency distribution: the shaded blue surface is a first-order transition surface, the thick red lines are second-order critical lines. The system is synchronized inside the region bounded by the surface, and is incoherent outside. The limits (the Kuramoto model, the Sakaguchi model and the BMF model) in which known transitions are obtained are labeled. (b) The known transition line for the Sakaguchi model, given by Eq. (17), showing also the Kuramoto model transition point, Eq. (16), for a Gaussian with zero mean and unit width width (). (c) The known transition line for the BMF model, given by Eq. (18). The shaded blue surface in (a) is bounded from above and below by the dynamical stability thresholds and of the synchronized and the incoherent phase respectively. These thresholds may be estimated in -body simulations from hysteresis plots (see Fig. 2 for an example); Panel (d) shows the surfaces and obtained from -body simulations with for a Gaussian with zero mean and unit width, with cuts of the three-dimensional plot at shown in panel (e) and at shown in panel (f).
Figure 2: (Color online) vs. adiabatically tuned for different values at , showing also the stability thresholds, and , for . For a given , the branch of the plot to the right (left) marked with an arrowhead pointing down (up) corresponds to increasing (decreasing); for , the two branches almost overlap. The data are obtained in -body simulations with for a Gaussian with zero mean and unit width.
Figure 3: (Color online) vs. adiabatically tuned for different temperatures at a fixed moment of inertia . For a given , the branch of the plot to the right (left) marked with an arrowhead pointing down (up) corresponds to increasing (decreasing); for , the two branches almost overlap. The data are obtained in -body simulations with for a Gaussian with zero mean and unit width. Similar disappearance of the hysteresis loop with increase of was reported in Ref. Hong:1999 ().

Iii Phase diagram

The complete phase diagram is shown schematically in Fig. 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 Fig. 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 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 ().

Figure 4: (Color online) For , and a Gaussian with zero mean and unit width, (a) shows at , the numerically estimated first-order phase transition point, vs. time in the stationary state, while (b) shows the distribution at several ’s around . The data are obtained in -body simulations with .

To confirm the first-order transition, we performed -body simulations involving integrations of Eq. (LABEL:eom-scaled) 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 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. 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 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 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 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 Goldenfeld:1992 ().

Iv Analytical treatment

We now turn to an analytical treatment of the first-order transition. In the continuum limit , the dynamics Eq. (LABEL:eom-scaled) 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 Acebron:2000 ()


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 single-oscillator distribution: Incoherent stationary state

Here, we start with deriving the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy equations for the dynamics Eq. (LABEL:eom-scaled) 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

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

  2. , 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. (LABEL:eom-scaled):


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


Note that the following normalizations hold for the single-oscillator distribution functions: , and .

Using Eq. (20) in Eq. (21) and simplifying, we get the BBGKY hierarchy equations for oscillators with frequencies as


and similar equations for for oscillators of frequencies . The first equations of the hierarchy are




In the limit of large , we can write


and express Eqs. (23) and (24) in terms of .

In order to generalize Eqs. (23) and (24) 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


so that Eq. (LABEL:eq:final-eqn) 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


that corresponds to canonical equilibrium, with determined self-consistently Chavanis:2013 (). For , the incoherent stationary state is Acebron:2000 ()


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


with . The solution of the linearized equation yields that satisfies Acebron:2000 ()


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


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.

Figure 5: The loop in the complex -plane, (b), corresponding to the loop in the complex -plane, (a), as determined by the function in Eq. (IV.2).

Considering and strictly positive, we multiply for convenience the numerator and denominator of Eq. (IV.2) 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 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


From the well-known 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


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


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. (IV.2), 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. (IV.2) 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. (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


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


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 first-order 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