# Synchronization transition in ensemble of coupled phase oscillators with coherence-induced phase correction

###### Abstract

We study synchronization phenomenon in a self-correcting population of noisy phase oscillators with randomly distributed natural frequencies. In our model each oscillator stochastically switches its phase to the ensemble-averaged value at a typical rate which is proportional to the degree of phase coherence . The system exhibits a continuous phase transition to collective synchronization similar to classical Kuramoto model. Based on the self-consistent arguments and on the linear stability analysis of an incoherent state we derive analytically the threshold value of coupling constant corresponding to the onset of a partially synchronized state. Just above the transition point the linear scaling law is found. We also show that nonlinear relation between rate of phase correction and order parameter leads to non-trivial transition between incoherence and synchrony. To illustrate our analytical results, numerical simulations have been performed for a large population of phase oscillators with proposed type of coupling. The results of this work could become useful in designing distributed networked systems capable of self-synchronization.

## I Introduction

The spontaneous synchronization of mutually coupled oscillators is observed in many complex biological, chemical, physical and sociological systems with different origins of rhythmical activity and different mechanisms of coupling. Depending on the context the term ”oscillator” may refer to neurons Hopfield_1995 (), cardiac pacemaker cells Peskin_1975 (), yeast cells Ghosh_1971 (), flashing fireflies Buck_1976 (), chirping crickets Sismondo_1990 (), applauding spectators Nikitin2001 (), nano-mechanical resonators Bagheri_2003 (), lasers Jiang_1993 (), superconducting Josephson junctions Wiesenfeld_1996 (), etc. The concept of spontaneous synchronization has proved to be useful in designing of artificial networked systems capable of self-organization in the absence of any centralized control mechanism Arenas_2008 (). A common example of these systems is the so-called wireless sensor networks that typically consist of a number of spatially distributed devices with computational, sensing, memory and wireless communication capabilities Hekmat_2006 (). Many industrial and civil applications of sensor networks require global clock synchronization of sensors Wu_2011 () or their distributed consensus on certain quantities of interest Olfati-Saber_2004 (). A good synchronization scheme should provide the robustness with respect to failure of an individual sensors and take into account the diversity of sensor’s properties and the presence of noise. The nature-inspired approach to this problem is completely distributed synchronization strategies based on the mutual coupling between sensors Simeone_2007 (); Carli_2008 (); Scutari_2008 (); Schenato_2011 (); Monti_2014 (). In self-organizing sensor network the nodes are coupled through the exchange of information: each node transmit its own state and receive the signals from the neighbors. The knowledge of this information allows to single node to update itself occasionally in accordance with some correction rule so as to facilitate the global synchrony. Another important technological application of synchronization phenomenon lies in the field of smart power grids Rohden_2012 (); Motter_2013 (). The modern electrical grids consist of thousands of generators and substations linked across large distances. To ensure stable operation and efficiency of the entire system, all the components that generate electricity must operate at the same frequency. In a decentralized smart power grid, the nodes could communicate among each other and correct their parameters in order to achieve and maintain the synchronized state.

Motivated by the discussed applications, in this work, we investigate an ensemble of phase oscillators for which the mutual coupling has a form of phase correction. Our model implies that each oscillator evolves autonomously except for discrete times when the instantaneous correction events occur. We focus our attention on the simple correction rule which prescribe to oscillator to reset its phase to the ensemble-averaged value at a typical rate which is determined by the product of the coupling constant and the degree of global phase coherence. This basic model is inspired by the famous Kuramoto model in which the effective coupling force is proportional to the amplitude of the mean field Kuramoto_1975 (). The central result of the work is a synchronization transition upon the change of coupling constant. We demonstrate that for sufficiently strong coupling a partially synchronized state continuously bifurcates from incoherence. It is also shown how the synchronization properties of oscillators change in the case of non-linear relation between the correction rate and the degree of synchrony.

## Ii Model formulation

Consider a diverse ensemble of noisy oscillators which are solely characterized by the phase variables , where . We will use the Kuramoto’s order parameter Kuramoto_book (); Kuramoto_1975 ()

(1) |

to describe the collective rhythm produced by the whole system. The magnitude of this complex parameter measures the degree of macroscopic coherence, while is the average phase. In the absence of coupling the dynamics of the th oscillator is given by

(2) |

where the intrinsic frequency is randomly chosen from some probability density and is the Gaussian white noise,

(3) |

The the angular brackets denote averages over statistics of the noise and the non-negative parameter measures its intensity. Dropping the index, we introduce the one-oscillator probability density of a phase distribution as , where is a solution of the stochastic equation (2) for a fixed realization of noise. The time evolution of this function obeys the following Fokker-Planck equation Risken ()

(4) |

Note that is nonnegative, -periodic in , and satisfies the normalization condition

(5) |

In the thermodynamic limit , the order parameter (1) can be expressed as

(6) |

For an arbitrary initial distribution the solution of Eq. (4) is given by the convolution , where the Green function

(7) |

obeys the Fokker-Plank equation (4) together with the initial condition . It follows from (7) that any initial state relaxes exponentially fast to the uniform distribution , which corresponds to zero order parameter (6). Predictably, the diverse population of non-interacting noisy units behaves incoherently.

Next, let us assume that oscillators are coupled and update itself in order to improve the degree of global synchrony. Namely, the th oscillator evolves accordingly to the equation (2), but from time to time it instantaneously shifts the phase to the current ensemble-averaged value . We treat the updating time instants as completely random and statistically independent for different oscillators. The intensity of phase correction is characterized by a typical rate of updates per oscillator which, in general, may change over time. Under these assumptions, the phase dynamics can be viewed as a drifting diffusion process on a circle in the presence of stochastic resetting Evans_2011 () with and playing a role of resetting position and resetting rate, respectively. We, thus, write the following master equation in the limit

(8) |

where is determined by (6) and the one-oscillator probability density implies the averaging of phase dynamics over both the noise and stochastic phase correction. The structure of the equation (8) is quite transparent: the first and second terms in the right hand side correspond to a diffusion and uniform rotation at natural frequency (cf. with Eq. (4)), while the third and fourth terms are related to updates and represent a negative probability flux out of each point and a corresponding positive probability flux into . Taking into account the -periodicity of function in it is easy to show the conservation of total probability (5).

The equation (8) constitutes the basis of our statistical model for coupled oscillators. The focus of this study is the synchronization properties of the system in the case when the correction rate is a function of the mean-field amplitude . It is useful to start with the description of the steady-state dynamics arising from the Eq. (8) with constant parameter .

## Iii Steady-state analysis

Suppose that the rate of updates is time-independent. Since the phase correction promotes the consistent behaviour of oscillators, the incoherent state does not solve the Eq. (8) at . If frequency distribution has a single maximum at frequency and is symmetric, it is reasonable to seek in the form of some stationary profile uniformly rotating at frequency . Then passing into the rotating frame and neglecting the time fluctuations of the order parameter at one can set without loss of generality. In result, the problem reduces to the stationary equation

supplemented by the periodic condition , where . It is straightforward to find the following solution

(9) |

in which and

(10) | |||

(11) |

It is noteworthy, that for this stationary distribution is nonequilibrium since phase correction produces ongoing circulation of probability.

We can now calculate the order parameter (6) by using (9-11). Since by assumption, one readily obtains

(12) |

In a particular case of a Lorentzian distribution we find exactly

(13) |

whereas for a normal distribution with standard deviation one obtains

(14) |

where erfc denotes the complementary error function.

Importantly, the result (12) is extracted from the equation (8) which is justified by the assumption that updates are the random Poisson events. There is a more general way to derive the steady-state coherence which is applicable to the case of arbitrary stationary statistics of updating time instants as long as these instants are independent for different oscillators. As previously, we suppose that density profile of oscillators in the phase circle is stationary in the the rotating frame, , and set . Then the oscillator dynamics may be thought of as a renewal process: each update resets the phase of oscillators to , whereas in the inter-update time intervals the evolution of phase is given by superposition of diffusion and detuning drift. The Green function (7) rewriting in the rotating frame yields the probability density of phase distribution for a randomly chosen oscillator

(15) |

where denotes the time elapsed since the last update of this oscillator. To obtain the steady-state distribution one should average the Eq. (15) over different realization of the phase correction process, i.e.

(16) |

where is the probability density of random variable . After the substitution the expressions (15), (16) and into (6) we derive the following simple formula

(17) |

which allows us to readily determine the stationary coherence once is known. For stochastic phase correction with Poisson rate of updates we have . Then integration of (17) over yields exactly (12).

As an example, let us apply the Eq. (17) to the case of independent periodical correction processes. We assume that the sequence of updating events of th oscillator is given by , where is updating period and is a random variable uniformly distributed over the interval . Then for and otherwise. Performing integration in (17) with one obtains

(18) |

In agreement with intuition, Eqs. (13), (14) and (18) give (incoherence) and (perfect synchrony) at and respectively.

## Iv State-dependent correction rate

We now turn to the more non-trivial case when the correction rate is not fixed but depends itself on the state of population. Let us assume that updates occur the more frequently, the more higher the current degree of phase synchrony . Apparently, this coupling scenario sets up a positive feedback so that as the population becomes more coherent, grows and so the correction rate increases, which tends to synchronize the oscillators even stronger. The situation is quite similar to classical Kuramoto model where the mean-field coupling becomes stronger with the growth of global synchrony. It is well-known that Kuramoto-like coupled oscillators exhibit spontaneous synchronization provided the interaction constant is large enough. The major question we wish to address here is whether there is a similar phenomenon in the discussed case of coupling through the coherence-induced phase correction. The answer is positive, below we demonstrate that for sufficiently strong coupling the phase correction compensates the destructive effects of noises and frequencyâs diversity and a partially synchronized state emerges. We focus our attention on the linear relation between coherence and correction rate, i.e. , where is the interaction constant. The nonlinear generalizations of the model will be discussed briefly in Section V.

### iv.1 Self-consistency analysis and numerics

One may expect that starting from an arbitrary initial condition the population of oscillators will settle into a steady state with constant coherence and average phase running at frequency . Then the correction rate becomes time-independent and the results of the previous section can be straightforwardly applied. Replacing by in Eq. (17), we find the following self-consistency condition for stationary coherence

(19) |

A trivial zero solution , corresponding to a incoherent state, holds for any value of . A second brunch of solutions bifurcates continuously at critical coupling

(20) |

obtained by letting , and corresponds to a partially synchronized state with non-zero .

Remarkably, formula (20) has the same structure as that for the critical coupling in the noisy Kuramoto model Sakaguchi_1988 (). At the same time, the models differ by the critical behaviour of the order parameter. An expansion of the right-hand side of Eq. (19) in powers of yields the linear scaling law in the vicinity of transition. In contrast, in Kuramoto model the order parameter of the bifurcating branch obeys the square-root law near threshold. Note also that how grows with strictly depends on the statistics of updates. For the sake of illustration, suppose that is a Lorentzian distribution. Then, for stochastic phase correction we find from (13) exactly at . However, for periodic phase correction with coherence-dependent frequency of updates , one obtains from (18) the non-trivial solution at .

The emergence of spontaneous synchronization is confirmed by simulations of large but finite population of oscillators. Our numerical scheme use temporal discretization of Eq. (2) with time step so that is the approximation of , where . Then, the complex order parameter is . The random forces are modeled by the telegraph processes whose correlation times are equal to time step duration . Specifically, the values of inside the th time step are chosen to be independent random constants with identical normal distributions. To ensure a given value of the diffusion coefficient , one should choose . At each time step we generate the set of random numbers having distribution with . If , we switch the phase of th oscillator to the current ensemble-averaged value . Thus, the rules for the evolution of the system are as follows

In numerics the initial phases of oscillators were uniformly distributed over the interval and the discrete time step was used. The diffusivity was and the natural frequencies were chosen from normal distribution with expected value and standard deviation . Figure 1 represents numerical results for the time-averaged coherence as a function of coupling constant in comparison with theoretical curve for infinite- system. It can be seen that level of synchrony in population remains low until reaches a critical value , above which monotonically increases towards its asymptotic value of unity. The non-zero values of below in the simulation reflect fluctuations due to finite .

### iv.2 Global dynamics in homogeneous system

The self-consistent arguments presented above do not provide any information about the stability of incoherent and partially synchronized states. The direct analysis of the stability properties is possible when the diversity of natural frequencies is negligible, i.e. . Then, using Eq. (8) one can derive the following closed-form equations associated with angular and radial motions of the order parameter

(21) | |||

(22) |

The equation (21) suggests that the average phase just uniformly rotates at frequency . Substituting we find easily the solution of amplitude equation (22)

(23) |

where is an initial condition. For the coherence always decays to zero as , but for incoherent state loses stability and one obtains another attracting point corresponding to partially synchronized state. Since by assumption, the critical value is consistent with (20).

If the natural frequencies of oscillators are non-identical, there is no closed description of global dynamics in terms of the order parameter. In the next subsection we propose perturbative approach which allows us to demonstrate that incoherent solution becomes linearly unstable as the strength of the coupling goes beyond the threshold (20).

### iv.3 Linear stability analysis of incoherence

In what follows we look at the oscillator dynamics in the frame moving with frequency . Let us consider a small perturbation of the uniform incoherent state

(24) |

where . At lowest (linear) order in the evolution of perturbation is governed by equation

(25) |

To proceed we write a -periodic function as a superposition of Fourier harmonics

(26) |

Note that normalization condition (5) automatically provides . When (26) is substituted into (25), we obtain the following equation for the amplitude

(27) |

The evolution equation for is just the conjugate of (27).

Next, substituting (26) into (1) yields

(28) |

Thus, only the first harmonic of the Fourier decomposition, which is the so-called fundamental mode, contributes to the order parameter.

From (27) and (28) one easily sees, that equation for the amplitude (and ) is uncoupled

(29) |

Let , then we obtain , where . The spectrum of the operator was constructed in Strogatz_1991 () in the context of noisy Kuramoto model Sakaguchi_1988 (). They show that the continuous part of spectrum lies in left half-plane for any , thus corresponding modes never cause instability. In contrast, the location of discrete spectrum depends strongly on interaction constant . When exceeds a threshold (20), the discrete real eigenvalue appears so that incoherent state becomes unstable.

Now we turn to the analysis of high order amplitudes. The equation (27) may be solved easily in terms of , and the initial condition . The result is

The order parameter is completely determined by the fundamental mode through Eq. (28), thus we put , where the eigenvalue belongs to the spectrum described above. After integration one obtains

(30) |

where and are the real and imaginary parts of , respectively. Notice that the denominator in (30) never vanishes since and for any Strogatz_1991 (). The second term in the rhs of (30) rapidly dies out, while the behaviour of the first one depends on the coupling strength . At the spectrum lies in the left half-plane so that always decays to zero. For the aforementioned positive discrete eigenvalue provides the exponential growth of .

Let us summarize the results of this subsection. The incoherent state goes linearly unstable to perturbation involving first harmonic for , where the critical coupling coincides with the prediction (20) of self-consistency analysis. Above the threshold the instability of the first harmonic induces also the growth of the high order amplitudes.

## V Nonlinear generalization

So far, we have considered the linear relation between correction rate and amplitude of the order parameter . One may expect more complex and interesting collective behaviour in system with non-linear dependence . Let us discuss as an illustration the quadratic model focusing on the case of identical oscillators, i.e. . Then, the amplitude equation (22) reads:

(31) |

At the only time-independent solution is and decays to zero as for any initial condition . However, for there are three fixed points: and . Thus, two partially synchronized branches bifurcate discontinuously from at , as it is shown in Fig. 2. The incoherence () and the upper synchronous state () are stable, and the lower synchronous one () is unstable. This means that the initial configurations with relax to incoherence as , while those with converge to the synchronized state at long-times. The described bistability is clearly observed in numerical simulations of finite- system, see Fig. 3. Note also that in more general case with the critical coupling is and two partially synchronized branches bifurcate from . Here again the system is bistable with respect to initial conditions.

For the second particular example the amplitude equation (22) yields

(32) |

The incoherent solution is unstable, while the only non-trivial solution is stable and exists for any . Thus infinitely small coupling leads to spontaneous synchronization of system. This feature is common for all models with .

As it was stated in the previous section, if oscillators have different natural frequencies, there is no closed evolution equation for the order parameter. Nevertheless, the stationary values of can be found by substitution into (20) and solving the resulting self-consistency equation. It is easy to see from (13) that when the frequency diversity is described by a Lorentzian distribution one should just replace by in the above formulas for identical oscillators. Notice, however, that these arguments do not indicate which of the found branches are stable.

## Vi Conclusion

In this work, we proposed and investigated a new model of coupled non-identical noisy phase oscillators. We assumed that the mutual coupling manifests itself as a tendency of each oscillator to adjust its phase to the ensemble-averaged value through instantaneous phase shifts which occur at random updating time instants. Then in the thermodynamic limit the time evolution of the one-oscillator probability density function is governed by the equation (4). The key ingredient of our model is the coherence-dependent rate of updates: the more pronounced the coherence of the current state of population, the more frequently the phase shifts occur. We mainly concentrated on the simplest linear case when the typical rate of updates is given by the coupling constant multiplied by the coherence. The distribution of natural frequencies of oscillators was assumed to be unimodal and symmetric. Our study revealed that the system exhibits the transition to spontaneous synchronization in a fashion similar to the mean-field Kuramoto model. Specifically, using the self-consistency arguments for an infinitely large population, we evaluated analytically the critical coupling strength (20) that marks the continuous onset of a partially synchronized solution. The same threshold follows from the linear stability analysis of incoherent state. In the vicinity of the critical value, the linear scaling of the order parameter is found. For the sake of illustration, we performed the numerical simulations of large number of oscillators undergoing coherence-induced phase correction. The analytical predictions for the degree of phase coherence in dependence on the coupling strength turned out to be in good agreement with numerical results, see Fig. (1).

For homogeneous population of oscillators, we found a closed description of global dynamics, which consists of two uncoupled evolution equations (21) and (22) for the amplitude and phase of the complex order parameter. Based on this result, we were able to demonstrate easily the global stability of a partially synchronized state in the system of identical oscillators.

Finally we considered the direct generalization of the model to the case of a nonlinear coupling. When the relation between the correction rate and the coherence is given by a power law with exponent larger than unity, the synchronization transition is discontinuous (Fig. 2) and above the threshold the oscillators exhibit bistability of global activity between incoherence and a partial synchronization (Fig. 3). If the correction rate grows with the coherence slower than linearly, the system shows the non-zero level of synchrony for arbitrary small coupling constant (Fig. 2).

The model presented here can be further extended in many ways. Obviously, the knowledge of the current states of nearby nodes cannot be assumed perfect in a networked applications because of communication and processing time delays and information losses. From the perspective of the proposed theoretical framework this means that the delta-function in the rhs of our basic equation (8) should be replaced by some spread function which depends on the previous states of the system. Another natural extension of the model includes the short-range interaction effects for regular and complex networks. In future studies, it is important also to consider the synchronization properties of system in which the updating time instants of different oscillators are correlated.

Here we have restricted ourselves to the case of ”attractive” phase correction which attempts to mutually synchronize the oscillators. There are many practical situations in which the formation of a coherent state in ensemble of oscillatory units is unacceptable and should be avoided Louzada_2012 (); Degesys_2008 (). One may expect that ”repulsive” phase correction ( i.e. the phase shifts to rather than to ) is a good strategy to break the undesired synchronization of Kuramoto-like coupled oscillators .

## References

- (1) J.J. Hopfield, Pattern recognition computation using action potential timing for stimulus representation Nature 376, 33 (1995).
- (2) C.S. Peskin, Mathematical Aspects of Heart Physiology ( Courant Institute of Mathematical Sciences, NYU, New York, 1975).
- (3) Ghosh, A. K., Chance, B., & Pye, E. K. (1971). Metabolic coupling and synchronization of NADH oscillations in yeast cell populations. Archives of biochemistry and biophysics, 145(1), 319-331.
- (4) J. Buck and E. Buck, Synchronous fireflies Sci. Am. 234, 74 (1976).
- (5) E. Sismondo, Synchronous, alternating, and phase-locked stridulation by a tropical katydid Science 249, 55 (1990).
- (6) A. Nikitin, Z. NÃ©da, and T. Vicsek, Collective Dynamics of Two-Mode Stochastic Oscillators Phys.Rev.Lett. 87(2), 024101 (2001).
- (7) M. Bagheri, M. Poot, L. Fan, F. Marquardt, and H. X. Tang, Photonic cavity synchronization of nanomechanical oscillators Phys. Rev. Lett. 111(21), 213902 (2013).
- (8) Z. Jiang and M. McCall, Numerical simulation of a large number of coupled lasers. JOSA B 10(1), 155 (1993).
- (9) K. Wiesenfeld, P. Colet, and S. H. Strogatz, Synchronization transitions in a disordered Josephson series array Phys. Rev. Lett. 76, 404 (1996).
- (10) Arenas, A., DÃaz-Guilera, A., Kurths, J., Moreno, Y., & Zhou, C. (2008). Synchronization in complex networks. Physics Reports, 469(3), 93-153.
- (11) R. Hekmat, Ad-hoc Networks: Fundamental Properties and Network Topologies, Springer, Berlin, Germany, 2006
- (12) Wu, Yik-Chung, Q. Chaudhari and E. Serpedin, âClock synchronization of wireless sensor networksâ, IEEE Signal Processing Magazine, vol. 28, no. 1, (2011), pp. 124-138.
- (13) R. Olfati-Saber, R.M. Murray, Consensus problems in networks of agents with switching topology and time-delays, IEEE Trans. Automat. Control 49 (2004) 1520â1533.
- (14) Simeone, O., Spagnolini, U. (2007). Distributed time synchronization in wireless sensor networks with coupled discrete-time oscillators. EURASIP Journal on Wireless Communications and Networking, 2007.
- (15) Carli, R., Chiuso, A., Schenato, L., & Zampieri, S. (2008, July). A PI consensus controller for networked clocks synchronization. In IFAC world congress on automatic control (IFAC 08) (Vol. 3).
- (16) Scutari, G., Barbarossa, S., & Pescosolido, L. (2008). Distributed decision through self-synchronizing sensor networks in the presence of propagation delays and asymmetric channels. Signal Processing, IEEE Transactions on, 56(4), 1667-1684
- (17) Schenato, L., Fiorentin, F. (2011). Average TimeSynch: A consensus-based protocol for clock synchronization in wireless sensor networks. Automatica, 47(9), 1878-1886.
- (18) Monti, M., Sanguinetti, L., Tschudin, C. F., & Luise, M. (2014). A Chemistry-Inspired Framework for Achieving Consensus in Wireless Sensor Networks. Sensors Journal, IEEE, 14(2), 371-382.
- (19) Rohden, M., Sorge, A., Timme, M. & Witthaut, D. Self-organized synchronization in decentralized power grids. Phys. Rev. Lett. 109, 064101 (2012).
- (20) Motter, A. E., Myers, S. A., Anghel, M., & Nishikawa, T. (2013). Spontaneous synchrony in power-grid networks. Nature Physics, 9(3), 191-197.
- (21) Y. Kuramoto, in International Symposium on Mathematical Problems in Theoretical Physics, edited by H. Araki (Springer Lecture Notes Phys., v. 39, New York, 1975) p. 420.
- (22) Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence ( Springer, Berlin, 1984).
- (23) H. Risken, Fokker-Planck Equation (Springer Berlin Heidelberg, Berlin, 1984).
- (24) M. R. Evans and S. N. Majumdar, Diffusion with stochastic resetting Phys. Rev. Lett., 106(16), 160601 (2011).
- (25) H. Sakaguchi, Cooperative phenomena in coupled oscillator systems under external fields Progr. Theor. Phys. 79, 39 (1988).
- (26) S. H. Strogatz and R. E. Mirollo, Stability of incoherence in a population of coupled oscillators Journal of Statistical Physics 63(3-4), 613 (1991).
- (27) V. H. P. Louzada, N. A. M. AraÃºjo, Jr, J. S. Andrade, and H. J. Herrmann, How to suppress undesired synchronization Scientific reports 2, (2012).
- (28) Degesys, J., Rose, I., Patel, A., & Nagpal, R. (2008). Self-organizing Desynchronization and TDMA on Wireless Sensor Networks. In Bio-Inspired Computing and Communication (pp. 192-203). Springer Berlin Heidelberg.