Exactly Solvable Phase Oscillator Models with Synchronization Dynamics

Populations of phase oscillators interacting globally through a general coupling function fsxd have been considered. We analyze the conditions required to ensure the existence of a Lyapunov functional giving close expressions for it in terms of a generating function. We have also proposed a family of exactly solvable models with singular couplings showing that it is possible to map the synchronization phenomenon into other physical problems. In particular, the stationary solutions of the least singular coupling considered, fsxd ­ sgnsxd, have been found analytically in terms of elliptic functions. This last case is one of the few nontrivial models for synchronization dynamics which can be analytically solved. [S0031-9007(98)07451-1]

The dynamics of systems consisting of large sets of mutual interacting units is a very interesting problem in statistical physics. In some cases, these units can be thought of as subsystems characterized by hidden internal degrees of freedom obeying their own dynamics. This is the case of large populations of individuals forming complex systems of interest in interdisciplinary fields such as biology, economy, neurophysiology, and ecology [1,2]. On the other hand, the subsystems can be seen as entities interacting with each other in the presence of quenched disordered external fields. The effect of the external forces is to pump energy into the system thereby leading to nonlinear oscillatory behavior. This last case describes charged wave instabilities in plasmas [3], dynamical response of Josephson junction arrays in the presence of an external ac field [4] or nonlinear oscillations, and coherent motion of magnetized domains in strongly coupled magnetic systems [5]. It was realized by Kuramoto [6] that this rich dynamical behavior can appear in mean-field models of phase oscillators with randomly quenched natural frequencies. The crucial feature of the Kuramoto model is that the random precession frequencies play the role of external forces which pump energy into the system. This leads to dynamical instabilities and synchronization in low dissipation regimes. Moreover, the mean-field character of the model captures the essential mechanism which leads to synchronized dynamics while retaining mathematical simplicity.
The purpose of this Letter is twofold: (i) to find Lyapunov functionals, and (ii) to show exact results for the synchronization dynamics of simple oscillator models. Exact results are very difficult to obtain in models with general smooth oscillator coupling functions. Here we follow a different strategy and analyze singular coupling models which can be mapped into well-known physical problems. Later we will go further and consider models with smoother couplings such as Daido's [7] which we analytically solve. This last case describes the dynamics of a charged plasma in the mean-field approximation. Thus this model is related to Lennard's for one-dimensional plasmas with global charge neutrality. The stationary solution thereof was found exactly in Ref. [8] using techniques different from ours. While our analysis is restricted to models without disorder, the techniques developed here could be extended to disordered cases [9]. We consider a system of nonlinearly globally coupled phase oscillators with random frequencies v i taken from a distribution g͑v͒ and subject to external independent white noise sources h i (of strength p 2T ): where i 1, . . . , N. Here f i ͑t͒, K 0 . 0, and f denote the ith oscillator phase, coupling strength, and a 2pperiodic real function, respectively. In Fourier space, the latter can be decomposed as f͑x͒ Pǹ 2`an e inx , with a n a ‫ء‬ 2n . Important work on the Kuramoto model with a generalized coupling function was carried out by Daido [7] who introduced the concept of order function and Crawford [10] who proposed scaling behavior of bifurcating branches from the incoherent solution. These works introduced rather general theories and methods which, however, had inherent limitations (zero temperature for the order function theory and vicinity to bifurcation points in Crawford's work). Thus it seems desirable to have exact results for solvable models where such theories could be checked and extended. 0031-9007͞98͞81(17)͞3643(4)$15.00 The moment approach and general equations.-A simple way to study the dynamics of the Kuramoto model with general coupling functions (1) is to use the recently introduced moment approach [11] and define H m k 1 N P N j1 ͗e ikf j ͘v m j , where the brackets ͗· · ·͘ denote average with respect to the external noise and the overbar denotes average with respect to the random oscillator frequency. The equations of motion for the moments read (2) The following differential equation for the generating function, g͑x, y, t͒ Pk 2`Pm0 e 2ikx y m 2pm! H m k ͑t͒, is easy to derive from (2): Here y͑x, t͒ is the drift velocity, defined by It is easy to check that y͑x, t͒ 2K 0 H͑x 2 V e t, t͒, where V e and H͑x, t͒ are the synchronization frequency (zero for stationary states) and Daido's order function [7], respectively. On the other hand, the generating function g͑x, y, t͒ is related to the one-oscillator probability density r͑x, v, t͒ by g͑x, y, t͒ R2`e yv r͑x, v, t͒g͑v͒ dv. We derive from (3) the usual nonlinear Fokker-Planck equation for r͑x, v, t͒ wherever g͑v͒ fi 0: Existence of a Lyapunov functional.-The solution of (5) is generally quite complicated: only in special cases is it possible to work out explicit results. Here we want to analyze under which conditions it is possible to find a Lyapunov functional for the dynamics. This is an important result since it explicitly yields a functional equation for the stationary states as well as stability results. In particular, we will find a Lyapunov functional when the coupling f͑x͒ is an odd function (i.e., if detailed balance is satisfied) which is symmetric about x p͞2, and there is no frequency disorder [i.e., g͑v͒ d͑v͒]. Then there exists a partition function which characterizes the thermodynamic properties of the stationary states. Previous work in this direction [12] was done for f͑x͒ sin x.
Let us define the potential function V ͑x, t͒ R x 2p y͑s, t͒ ds and restrict ourselves to finding stationary solutions; rotating wave solutions may be reduced to this case after moving to a rotating frame. It is possible to show that the stationary solutions should have the form where Z and J are independent of x. We now impose the condition that r be a 2p-periodic function of x, use the symmetry properties of the drift, and find the probability flux J as a function of Z: J͞T 2Z 21 sinh͑2pT 21 v͒͞ R p 2p e 2͓͑p1x͒v1V ͑x͔͒͞T dx. Z can be found from the normalization condition R p 2p rdx 1. An interesting case corresponds to the case without disorder, g͑v͒ d͑v͒, for which J 0, r͑x, 0, t͒ g͑x, y, t͒ ϵ g͑x, t͒. We can obtain the stationary solutions (subscript zero below) by solving the system of equations: It is possible to show that, for this type of potential solutions, there is a Lyapunov functional [13] defined by the relative entropy, g͑x, t͒ e ͓V ͑x,t͒2m͑t͔͒͞T , Direct calculation (to be presented elsewhere) shows that h͑t͒ is bounded from below, h 0 ͑t͒ # 0, and that g tends to a stationary solution proportional to g. We then conclude that, for odd coupling functions and in absence of disordered frequencies, the stationary state (7) corresponds to the thermodynamic equilibrium state of a system of N oscillators interacting through the Hamiltonian where´͑x͒ is a pair interaction energy 2p-periodic function defined by´͑x͒ R x 2p f͑s͒ds. Note that V ͑x͒ 2K 0 R p 2p g͑x͒´͑x͒dx 1 const. Thus a computation of the partition function for (12) yields the stationary states. It is important to stress that potential solutions of the type (7) are no longer stationary solutions of the dynamics in the presence of disordered frequencies v i [11] or when detailed balance is violated [i.e., f͑x͒ is not an odd function]. The physical meaning of these two conditions is quite clear: suppose a single oscillator performs a global rotation of angle 2p. The global energy of the system, given by (12), does not change because´͑x͒ is 2p periodic. Thus the coupling strength does not exert work into the system of oscillators when a rotation path ͑f i ! f i 1 2p͒ is performed. Indeed, the amount of work is W Singular coupling models.-For general coupling functions, with infinitely many nonvanishing Fourier modes a n , it is difficult to find explicit analytical expressions for the stationary states. But calculations turn out to be much simpler for singular coupling functions. Here we analyze three different cases: (a) f͑x͒ d 0 ͑x͒; (b) f͑x͒ d͑x͒, and (c) f͑x͒ sgn͑x͒. In case (a), we have y͑x͒ 2K ≠g͞≠x, This equation is related to porous media systems with the additional 2p-periodicity condition for g. It is easy to check that h͑t͒ ͑1͞2͒ R p 2p g 2 ͑x, t͒dx $ 0 is a Lyapunov functional because h 0 ͑t͒ 2 R p 2p ͑T 1 K 0 g͒ ͑≠g͞≠x͒ 2 dx # 0. Then the incoherent solution g͑x͒ 1͞2p is globally stable and model (a) does not synchronize at any T .
Case (b) is described by the Burgers equation, for a 2p-periodic g, and it can be solved by the Hopf-Cole transformation. h͑t͒ defined for case (a) is also a Lyapunov functional for case (b), but now h 0 ͑t͒ 2T R p 2p ͑≠g͞≠x͒ 2 dx # 0. At finite T , incoherence is again stable. But, at zero temperature, initially synchronized oscillators remain forever synchronized.
Analytical solution of the model with Daido coupling.-Model (c) was originally proposed by Daido in the general disordered case but we will solve it for g͑v͒ d͑v͒. Since f͑x͒ is odd and symmetric about x p͞2, it is in principle possible to find the stationary states by solving the functional Eqs. (7) and (8). Such a calculation is quite involved and here we follow a different and novel approach. The dynamical equations for model (c) are nonlocal because the drift velocity (4) is but they become local after inserting the definitions r͑x, t͒ 2K 0 ͓g͑x, t͒ 2 g͑x 1 p, t͔͒ , s͑x, t͒ 2K 0 ͓g͑x, t͒ 1 g͑x 1 p, t͔͒ . (16) The new system is The drift velocity y͑x, t͒ is given by ≠y͞≠x 2r. All functions are 2p periodic and it is easy to check that both r and y are antiperiodic if x ! x 1 p, while s is p periodic. The stationary solutions of (17) and (18) satisfy where L is a positive constant. Equation (20) may be integrated once yielding where w͑j͒ y͑x͒͞ p L , x T j͞ p L, and E is a new constant. Equation (21) can be easily solved by quadratures in terms of Jacobi elliptic functions [14]. The periodicity properties of the drift y yield finitely many solutions (characterized by their value of E and an odd integer n) for a fixed value of K 0 ͞T . The relation among E , n, and K 0 ͞T is Here K͑m͒ and E͑m͒ are the complete elliptic integrals of the first and second kind with parameter m ͑1 2 p 1 2 2E ͒͑͞1 1 p 1 2 2E ͒, respectively. For every admissible value of E and n, with K 0 ͞T fixed, the drift velocity and the probability density can be found from (16) Here u 2nK͑m͒ ͑x 2 x 0 ͒͞p, x 0 const, while sn, cn, and dn are Jacobi elliptic functions defined in [14]. The details of these calculations will be presented elsewhere.
As K 0 ͞T !`the number of possible stationary solutions becomes infinite. Each different stationary state belongs to a different solution branch which bifurcates from incoherence at a critical value of the coupling K 0 n 2 T p͞2, where n 2p 2 1 and p $ 1 is an integer (see Fig. 1). The first branch bifurcates from incoherence at K 0 ͞T p͞2 and remains stable for all larger K 0 ͞T. In the limit K 0 ͞T !`, y͑x͒͞K 0 sgn͑x 2 x 0 ͒, g͑x͒ d͑x 2 x 0 2 p͒, and full synchronization is achieved. A convenient synchronization order parameter r is defined through the global magnetization M re ia ͑1͞N͒ P N j1 e if j . The order parameter can be calculated from Monte Carlo simulations of the Hamiltonian (12) with´͑x͒ jxj and K 0 1. The simulations use the heat-bath algorithm with a random sequential updating of the phases. The transition temperature corresponding to the first branch is easily obtained through standard finite-size scaling methods. Consider the kurtosis (or Binder parameter) for the synchronization parameter r, g where ͗· · ·͘ is the standard configurational average [weighted with the usual Boltzmann-Gibbs factor, exp͑2bH ͒, and H is given by (12)]. The curves for g are shown in Fig. 2 for different sizes. Note that these curves (especially for N 50, 100, and 500; data for N 1000 is more noisy) intersect at a common point characterizing the bifurcation temperature.
In summary, we have considered models of oscillators coupled through a general function f͑x͒. Without precessing frequencies and for odd-coupling functions there exists a Lyapunov functional. Then the probability den- sity evolves toward stable stationary states which can be described by an equilibrium measure. We have also proposed a family of exactly solvable models with singular couplings. The oscillators may synchronize more easily for less singular couplings. In particular, case (a) (the porous medium oscillators), f͑x͒ d 0 ͑x͒, never synchronizes, case (b) (Burgers oscillators), f͑x͒ d͑x͒, retains synchronization only at zero temperature, and case (c) (charged plasma in the mean-field approximation), f͑x͒ sgn͑x͒, synchronizes at a finite temperature. This last case is one of the few nontrivial models for synchronization dynamics which can be analytically solved. Stationary synchronized phases bifurcate from incoherence at a finite temperature which corresponds to a thermodynamic second order phase transition. It would be very interesting to extend our considerations to the case of random frequencies for models with a synchronization transition at finite temperature, where we expect dynamical instabilities as well as nonlinear behavior to appear. This work has been supported by DGES (