Finite-size-induced transitions to synchrony in oscillator ensembles with nonlinear global coupling

# Finite-size-induced transitions to synchrony in oscillator ensembles with nonlinear global coupling

Maxim Komarov Department of Physics and Astronomy, University of Potsdam, Karl-Liebknecht-Str 24, D-14476, Potsdam, Germany Department of Cell Biology and Neuroscience, University of California Riverside, 900 University Ave. Riverside, CA 92521, USA Department of Control Theory, Nizhni Novgorod State University, Gagarin Av. 23, 606950, Nizhni Novgorod, Russia    Arkady Pikovsky Department of Physics and Astronomy, University of Potsdam, Karl-Liebknecht-Str 24, D-14476, Potsdam, Germany Department of Control Theory, Nizhni Novgorod State University, Gagarin Av. 23, 606950, Nizhni Novgorod, Russia
July 3, 2019
###### Abstract

We report on finite-sized-induced transitions to synchrony in a population of phase oscillators coupled via a nonlinear mean field, which microscopically is equivalent to a hypernetwork organization of interactions. Using a self-consistent approach and direct numerical simulations, we argue that a transition to synchrony occurs only for finite-size ensembles, and disappears in the thermodynamic limit. For all considered setups, that include purely deterministic oscillators with or without heterogeneity in natural oscillatory frequencies, and an ensemble of noise-driven identical oscillators, we establish scaling relations describing the order parameter as a function of the coupling constant and the system size.

###### pacs:
05.45.Xt, 05.45.-a

Collective synchronization phenomena are abundant in complex nonlinear systems, and onset of synchrony can be typically treated as a nonequilibrium phase transition. The Kuramoto model Kuramoto (1975); *Acebron-etal-05 of globally coupled phase oscillators is the simplest paradigmatic system, where this transition can be explored nearly in full details Strogatz (2000); *Ott-Antonsen-08, also a relation to equilibrium transitions is well studied Gupta et al. (2014). This model is universally applicable for ensembles of weakly coupled oscillators, possessing harmonic phase sensitivity (like, e.g., Josephson junctions Wiesenfeld et al. (1996)). In many cases one, however, needs to go beyond such a simple setup, allowing for couplings that include higher harmonics, this is relevant for electrochemical oscillators Kiss et al. (2005) and -Josephson junctions Goldobin et al. (2011). Moreover, as suggested in Rosenblum and Pikovsky (2007) and experimentally realized in Temirbayev et al. (2012); *Temirbayev_et_al-13, coupling terms can be nonlinear functions of the order parameters (mean fields).

In this letter we describe nontrivial properties of the synchronization transition in a model with simple nonlinear coupling, where the coupling is at the second harmonics of the phase, but is proportional to the the square of the Kuramoto order parameter. We will show that such an interaction on a microscopic level represents a fully connected hypernetwork. By performing the analysis in the thermodynamic limit, we will demonstrate, that for deterministic ensembles the asynchronous state with a uniform distribution of the phases never loses stability; for noisy oscillators it is possible to show that such an asynchropnous state is the only stationary solution. Nevertheless, direct numerical simulations of finite pupulations yield partially synchronous regimes. These regimes can be called finite-size-induced one; the main goal of this paper is to clarify their nature. Our main tool is the analysis of scaling of partial synchrony with the system size. We will establish scaling properties of this finite-size-induced transition for different setups: identical purely deterministic oscillators, identical noisy oscillators, and deterministic non-identical ensemble (the latter setup is mostly close to the Kuramoto model). From these scaling properties it follows, that in the thermodinamic limit the characteristic value of the observed order parameter tends to zero, while the critical value of the coupling strength at which partial synchrony is observed, diverges. Thus, this transition to synchrony can be called finite-size-induced one.

Let us start with formulation of general phase equations for an ensemble of nonlinear oscillators coupled via mean fields. In the limit of weak coupling (or weak external forcing), the dynamics of each oscillator can be described by its phase via the phase response function :

 ˙ϕ=Ω+ω+S(ϕ)F,

here stands for the population mean natural frequency and is an individual deviation from the mean. The term is an overall force produced by mean field coupling, i.e. it is a general function of mean fields (Daido order parameters Daido (1993); *Daido-96) (here means averaging over the ensemble) which can be represented as an expansion (where are constants). Next, let us introduce a slowly varying phase and the corresponding slow order parameters . Using a Fourier representation of the phase response function we get the following general equation for the slow phase:

 ˙φ=ω+∑n,k,msne−inφhk,m¯Zmkei(km−n)Ωt. (1)

Performing an averaging of the equation (1) over the fast time scale is equivalent to keeping only terms on the r.h.s. that do not contain explicit time dependence , i.e. we have to set :

 ˙φ=ω+∑k,mskme−ikmφhk,m¯Zmk. (2)

This is the general form of the phase equation derived for a weakly coupled oscillatory ensemble with generically nonlinear (due to terms with ) mean-field interaction (cf. Rosenblum and Pikovsky (2007)). Recall that the simplest case where gives us the term , and (2) reduces to the standard Kuramoto-Sakaguchi model Sakaguchi and Kuramoto (1986). Considering further terms with , we get a generalized bi-harmonic coupling on the r.h.s of the system (2). The case where the bi-harmonic coupling depends linearly on the order parameters () was extensively studied in Refs. Komarov and Pikovsky (2013); *Komarov-Pikovsky-14; *Vlasov-Komarov-Pikovsky-15.

In this letter we focus on the effects produced by purely nonlinear second-harmonic coupling (see Skardal et al. (2011) for the analysis of linear second-harmonic coupling ) and show that it is responsible for finite-size induced transitions to synchrony with nontrivial scaling on the ensemble size. We will demonstrate that while synchrony disappears in the thermodynamic limit, it is observed for finite ensembles. Thus, throughout this paper we consider the following model, where also external noise is taken into account for completeness:

 ˙φk=ωk+εR2sin(2Θ−2φk)+√Dηk(t),   k=1,...,N. (3)

Here we denote , is the size of population, and the noise is Gaussian delta-correlated . Qualitatively, nontrivial features of this model can be understood as follows. Because the interaction is proportional to the second harmonics of the phase , it supports formation of two clusters, with phase difference . However, the coupling term is determined by the order parameter which vanishes if two symmetric clusters are formed and is non-zero only due to asymmetry of the clusters. This asymmetry, as we shall demonstrate below, is due to finite-size fluctuations at the initial stage when the clusters are formed from the disorder.

Before proceed with the main analysis, we discuss physical relevance of the model (3). First, the purely second-harmonic coupling () appears when the force acts on nearly harmonic oscillators parametrically (a typical example here are mechanical pendula suspended on a vertically oscillating common beam). Another situation where the second harmonic in dominates is that of period-doubled oscillations. Noteworthy, due to nonlinear coupling model (3) represents a hypernetwork Wang, Jian-Wei et al. (2010); *Bilal-Ramaswamy-14 of oscillators. Indeed, substituting the expression for the mean field in (3), one can see that the microscopic coupling terms can be written as . This means that effective interactions are not pairwise (as in the standard Kuramoto model and its numerous generalizations) but via triplets; this is exactly the definition of a hypernetwork coupling structure.

Furthermore, it is worth discussing the rôle of different order parameters in the problem. The order parameter governs the force acting on the oscillators and is therefore of major importance. Because this force contains the second harmonics only (), the appearing order is of “nematic” type and corresponds to large absolute values of the order parameter ; at these states the order parameter may be small. In the disordered, asynchronous states both order parameters are small.

Linear stability analysis of a disordered state () in model (3) is straightforward, beacuse coupling is nonlinear in and thus does not contribute. The solution is the same as for the Kuramoto model with vanishing coupling Strogatz and Mirollo (1991); *Crawford-Davies-99: the disordered state is either neutrally (without noise ), or asymptotically stable. Thus, all the transitions described below are due to nonlinear and finite-size effects.

We start with the simplest case where all oscillators have identical frequencies () and are not affected by noise (). In this case, for any there are two stable positions for the phases: and . Any distribution with oscillators in the first state is possible, with order parameter

 R=2n1/N−1. (4)

Only the symmetric distribution with is not a solution, because the mean field vanishes. Therefore, in the thermodynamic limit , the stationary two-cluster distributions can be written as

 ρ(φ)=Sδ(φ−Θ)+(1−S)δ(φ−Θ+π), (5)

with an arbitrary indicator constant , the order parameter is and is arbitrary.

A nontrivial question here is: which of possible synchronous states establishes if one starts from a fully disordered initial configuration with uniformly distributed initial phases . Numerically obtained distributions of the final states, for different sizes , are shown in Fig 1(a). Here the order parameter can attain only discrete values according to (4), and are probabilities of these states. Remarkably, these distributions collapse perfectly after rescaling , , as is shown in Fig. 1(b). This means that the stationary order parameter scales as , i.e. it disappears in the thermodynamic limit. To this scaling corresponds also the scaling of the characteristic transient time from initial disorder to a final synchronous configuration: as one can see from (3), this time is , what is confirmed by numerics (not shown).

Next, we consider the case when the oscillators have different natural frequencies and are not affected by noise . We assume a Gaussian distribution (without loss of generality the width of the distribution is set to 1 and the mean frequency to 0). First, following Refs. Komarov and Pikovsky (2013); *Komarov-Pikovsky-14; *Vlasov-Komarov-Pikovsky-15, we find stationary solutions in the thermodynamical limit by virtue of a self-consistent scheme. We will see, that although the analysis of the thermodynamic limit does not provide a transition to partial synchrony, it allows us to find states close to that observed in simulations of finite systems. For the sake of brevity of presentation we restrict to symmetric, non-rotating solutions only, then without loss of generality one can set . For such states the conditional distribution density is stationary, and the order parameter can be defined as follows:

 R=∬dωdφg(ω)ρ(φ|ω)cosφ. (6)

It is convenient to introduce an auxiliary parameter (the overall amplitude of the coupling function) and the rescaled frequency . It easy to see from (3) that all the oscillators can be divided into those locked by the mean field () and the unlocked (rotating) ones . The distribution of the latter ones is inversely proportional to the phase velocity (here is a normalization constant) and, because of the symmetry, it does not contribute to the order parameter in (6). The distribution of the locked oscillators is in fact the same as in (5), but frequency-dependent:

 ρL(φ|x) =S(x)δ(φ−Φ)+(1−S(x))δ(φ−Φ−π) (7)

where . Similar to the case of identical oscillators, the indicator function is arbitrary (due to assumed symmetry we restrict ourselves to the case , asymmetric functions lead generally to rotating solutions), it describes redistribution of oscillators between two stable branches and . Below we consider the simplest case of constant indicator function . In order to get the final closed self-consistent scheme, we need to substitute the distribution function (7) into the definition of the order parameter (6). This yields the order parameter as a function of the coupling constant in a parametric (in dependence on auxillary parameter ) analytic form:

 R =A(2σ−1)H(A),   ε=AR−2, (8) H(A) =∫π4−π42cosφcos2φg(Asin2φ)dφ.

Figure 2(a) illustrates stationary solutions at different indicator constants . The black curves denote the main solution at . At a critical coupling two branches, stable (black solid) and unstable (dashed line) arise (stability is determined by direct finite-size numerical simulations). Note that these lines are separated from the disordered solution , although as at the unstable branch. Solutions for can be easily found from the rescaling of the main dependence (at ) according to (8). So, the curves at have qualitatively similar structure, however, they are shifted to larger values of for smaller values of (see dotted red lines in Fig. 2(a)). In particular, the critical points scale as . This blue bold line in Fig. 2(a) together with the black solid line at define the region of possible synchronous solutions, characterized by different indicator constants .

It is worth mentioning that the incoherent solution exists at any value of coupling and it is always stable in the thermodynamical limit. However, in the finite-size simulations of the system (3), we found that the incoherent state has a finite lifetime: after a long transient a synchronous state from the described above solution family (i.e. between the blue bold and black solid curves in Fig. 2(a)) sets on, and this state remains further stable. Blue markers in the Fig. 2(a) denote averaged value of obtained from direct numerical simulations of (3) with . The averaging was performed over distinct simulation runs (until final time ) with disordered initial conditions. The final state to which the system jumps from the incoherence is not always the same and has a deviation range depicted by the gray area in Fig. 2(a).

A more detailed description of the final synchronous state is presented in the Fig. 2(b). Here the top panel shows dependence of the phases on the internal frequency . The area in the center clearly shows two stable branches of locked oscillators. Outside this area one can see the clouds of points, which depict unlocked oscillators. The unlocked phases rotate in relation to the mean field phase, and, therefore, constitute an asynchronous part of the ensemble. The bottom panel shows statistics of the function , which is calculated in the range where oscillators are typically locked to the mean field. The function has certain profile depicted in the bottom panel of the Fig. 2(b). As one can see, the distribution of locked oscillators over the branches remains close to a constant value in the center of the range, however, it drops significantly close to the boundaries of the coherent region.

The averaged lifetime of the incoherent state drastically increases with decrease of coupling constant , what is shown in the inset of Fig. 2(a). Thus, below it is impossible to collect any reasonable statistics with a finite simulation time . However, even for relatively small values of , the transition from the incoherent state is possible, what is shown by the blue markers to the left of they gray colored area.

It is instructive to characterize this finite-size induced transition to synchrony in dependence on the system size . The dependence of averaged lifetime on the rescaled coupling , plotted in the Fig. 2(c), demonstrates a nice collapse of data points. This scaling follows from the fact that the characteristic amplitude of the coupling term is and in the disordered state. Furthermore, in the inset of the Fig. 2(c) we show a rescaled order parameter, obtained at the end of a fixed integration time ; one can see that it scales as , what can be explained as follows. For sufficiently small values of the coupling the system always exhibits finite-size fluctuations and remains in the asynchronous state. With increase of , the transition to synchronous state becomes more probable what leads to an increase of the averaged final value of . The upper branch reflects the scaling of the lowest border of synchronous states mentioned above. Note that the critical value of coupling resulting from this scaling is , so the transition effectively disappears in the thermodynamic limit.

Finally, we describe the finite-size-induced transitions to synchrony in the ensemble of identical oscillators () with noise (3). Without lose of generality, we can take , so that the only parameters are and . In the thermodynamic limit, when , the system does not have any non-trivial coherent solution, because due to the symmetry of the coupling function, the stationary density that follows from the Fokker-Planck equation is also symmetric (in the small noise limit it tends to (5) with ), thus the only stationary solution is the incoherent one with . However, similar to the situations described above, for small system sizes a transition to synchronous two-cluster configurations is observed (cf. Ref. Pikovsky et al. (1994)). In contradistinction to the noise-free case, here also back transitions to disorder are possible due to noise, so that at a long run the process looks like an intermittent order-disorder dynamics.

Qualitatively, this dynamics can be understood as effect of noise on the multiplicity of synchronous states described for the noise-free case above (cf. discussion around Eq. (5)): due to small noise now transitions between these states occur. The transitions rates one can estimate using the Kramers’ formula, they are exponentially small in the potential barrier which is where . For small , the Kramers’ rate does not apply, here one can phenomenologically set the transition rate to a constant. As a result, one obtains for the order parameter a random walk model, which can be described by the corresponding master equation. Without going into details, which will be presented elsewhere, we present here the main results of this statistical model. The stationary distribution (Fig. 3(a)) shows a transition to synchrony at , so that for larger couplings we get , while below this threshold only finite-size fluctuations around disordered state with are observed. The characteristic time scale of the time evolution from asynchrony to synchrony is however extremely large, because the Kramers’ rate at large is exponentially small. Thus, direct simulations of finite ensembles, started from a disordered state and performed over a finite time interval , allow to reveal only order parameters for which , thus . At this stage the evolution becomes effectively “frozen”. We illustrate this in Fig. 3(b): only for one observes a saturation of the order parameter as predicted by the random walk model, while for larger , values of order parameter close to one are never achieved during available integration times. Of course, if one starts from a state with close to one, it remains practically frozen as well.

Summarizing, in this letter we considered a model of oscillators, globally coupled via a nonlinear function (in our case, square) of the Kuramoto mean field. Equivalently, on the microscopic level such a coupling can be described as a fully connected hypernetwork. While the disordered state remains stable in the linear approximation and in the thermodynamic limit, a transition to synchrony is observed due to finite-size effects: the characteristic critical coupling parameter value scales typically as , also the transient time from disorder to order diverges as . For the deterministic ensembles we demonstrated scaling properties of the transition in form of dependence of the order parameter on the coupling strength and the ensemble size. For the noisy case, the system demonstrates effective breaking of ergodicity, being trapped in frozen metastable states due to exponentially small hopping rates. While we focused on the purely quadratic in mean field coupling, the described framework allows one also to consider a general combination of linear and nonlinear couplings. The approach based on the master equation provides a framework for a description of finite-size transitions not only in the context of phase oscillator networks, but in other types of mean-field coupled systems demonstrating finite-size-induced transition Pikovsky et al. (1994).

###### Acknowledgements.
M.K. thanks Alexander von Humboldt foundation, NIH (R01 DC012943), ONR (N000141310672) and the Russian Science Foundation (Project 14-12-00811) for support. A. P. was supported by COSMOS ITN (EU Horizon 2020 research and innovation programme under Maria-Sklodowska-Curie grant agreement No 642563) and by the grant (agreement 02.Ð.49.21.0003 of August 27, 2013 between the Russian Ministry of Education and Science and Lobachevsky State University of Nizhni Novgorod).

## References

• Kuramoto (1975) 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.
• Acebrón et al. (2005) J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort,  and R. Spigler, Rev. Mod. Phys. 77, 137 (2005).
• Strogatz (2000) S. H. Strogatz, Physica D 143, 1 (2000).
• Ott and Antonsen (2008) E. Ott and T. M. Antonsen, CHAOS 18, 037113 (2008).
• Gupta et al. (2014) S. Gupta, A. Campa,  and S. Ruffo, J. Stat. Mech. - Theor. Exp. 8, R08001 (2014).
• Wiesenfeld et al. (1996) K. Wiesenfeld, P. Colet,  and S. H. Strogatz, Phys. Rev. Lett. 76, 404 (1996).
• Kiss et al. (2005) I. Z. Kiss, Y. Zhai,  and J. L. Hudson, Phys. Rev. Lett. 94, 248301 (2005).
• Goldobin et al. (2011) E. Goldobin, D. Koelle, R. Kleiner,  and R. G. Mints, Phys. Rev. Lett. 107, 227001 (2011).
• Rosenblum and Pikovsky (2007) M. Rosenblum and A. Pikovsky, Phys. Rev. Lett. 98, 064101 (2007).
• Temirbayev et al. (2012) A. A. Temirbayev, Z. Z. Zhanabaev, S. B. Tarasov, V. I. Ponomarenko,  and M. Rosenblum, Phys. Rev. E 85, 015204(R) (2012).
• Temirbayev et al. (2013) A. A. Temirbayev, Y. D. Nalibayev, Z. Z. Zhanabaev, V. I. Ponomarenko,  and M. Rosenblum, Phys. Rev. E 87, 062917 (2013).
• Daido (1993) H. Daido, Physica D 69, 394 (1993).
• Daido (1996) H. Daido, Physica D 91, 24 (1996).
• Sakaguchi and Kuramoto (1986) H. Sakaguchi and Y. Kuramoto, Prog. Theor. Phys. 76, 576 (1986).
• Komarov and Pikovsky (2013) M. Komarov and A. Pikovsky, Phys. Rev. Lett. 111, 204101 (2013).
• Komarov and Pikovsky (2014) M. Komarov and A. Pikovsky, Physica D 289, 18 (2014).
• Vlasov et al. (2015) V. Vlasov, M. Komarov,  and A. Pikovsky, J. Phys. A: Mathematical and Theoretical 48, 105101 (2015).
• Skardal et al. (2011) P. S. Skardal, E. Ott,  and J. G. Restrepo, Phys. Rev. E 84, 036208 (2011).
• Wang, Jian-Wei et al. (2010) Wang, Jian-Wei, Rong, Li-Li, Deng, Qiu-Hong,  and Zhang, Ji-Yong, Eur. Phys. J. B 77, 493 (2010).
• Bilal and Ramaswamy (2014) S. Bilal and R. Ramaswamy, Phys. Rev. E 89, 0622923 (2014).
• Strogatz and Mirollo (1991) S. H. Strogatz and R. E. Mirollo, J. Stat. Phys. 63, 613 (1991).
• Crawford and Davies (1999) J. D. Crawford and K. T. R. Davies, Physica D 125, 1 (1999).
• Pikovsky et al. (1994) A. S. Pikovsky, K. Rateitschak,  and J. Kurths, Z. Physik B 95, 541 (1994).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters