Mean Field and Mean Ensemble Frequencies of a System of Coupled Oscillators
Spase Petkoski, Dmytro Iatsenko, Lasko Basnarkov, Aneta Stefanovska
aa r X i v : . [ n li n . AO ] F e b Mean Field and Mean Ensemble Frequencies of a System of Coupled Oscillators
Spase Petkoski , Dmytro Iatsenko , Lasko Basnarkov , , and Aneta Stefanovska ∗ Department of Physics, Lancaster University, Lancaster LA1 4YB, United Kingdom Ss. Cyril and Methodius University, Faculty of Computer Science and Engineering, P.O. Box 393, Skopje, Macedonia and Macedonian Academy of Sciences and Arts, Skopje, Macedonia (Dated: September 18, 2018)We investigate interacting phase oscillators whose mean field is at a different frequency fromthe mean or mode of their natural frequencies. The associated asymmetries lead to a macroscopictravelling wave. We show that the mean ensemble frequency of such systems differs from theirentrainment frequency. In some scenarios these frequencies take values that, counter-intuitively, liebeyond the limits of the natural frequencies. The results indicate that a clear distinction shouldbe drawn between the two variables describing the macroscopic dynamics of cooperative systems.This has important implications for real systems where a non-trivial distribution of parameters iscommon.
PACS numbers: 05.45.Xt, 89.75.-k
I. INTRODUCTION
Systems consisting of large numbers of interactingunits are common in science and nature, and have beenthe essential modelling tools in physics, biology, chem-istry and social science [1]. Here, the state of the wholesystem is characterized with macroscopic variables suchas the temperature, magnetization and so on. The valuesof these depend on both the microscopic laws governingthe dynamics of the units as well as their interaction.Generally, only macroscopic variables are accessible inexperiments. Thus their precise definition and interpre-tation is needed.In the case of populations of weakly interacting oscil-lators, application of the phase approximation leads tothe Kuramoto model (KM) for globally coupled phaseoscillators [2]. Although the model itself represents anidealized scenario, its analytical tractability makes it theprevailing approach in tackling a wide variety of impor-tant problems – from Josephson-junction arrays [3] tobrain dynamics under anaesthesia [4] and pedestrian in-duced oscillations in the Millenium bridge problem [5].This has led to many extensions of the basic model to al-low more realistic descriptions of actual systems, e.g. KMunder influence of external fields [6], or with time-varyingparameters [7] (for a review of generalizations and theproblems they address see [8] and references therein).A fundamental feature of this model is that for a largeenough coupling, synchronized behaviour emerges. De-pending on their inherent frequencies, some of the oscil-lators become locked, while the others continue to rotateasynchronously but with adjusted frequencies. The de-gree of the synchronization is usually characterized bysome order parameter. For example, in the paradigmaticexample of flashing fireflies [1], this parameter would de-scribe the fraction of fireflies that flash in synchrony. ∗ [email protected] Since the building units of the KM are defined by theirnatural frequencies, the macroscopic dynamics of the os-cillating system must be also characterized by some av-erage frequency. However, two quantities can be used:the effective frequency to which synchronized oscillatorsare locked and the average frequency of all the oscilla-tors, locked and unlocked, that belong to the observedsystem. The former represents the natural macroscopicfrequency, whilst the latter is the microscopically aver-aged mean frequency. We shall call these respectively the mean field frequency and the mean ensemble frequency .Returning to the example of fireflies, the frequency ofthose flashing in synchrony – the mean field frequency –can generally deviate from the observed mean frequencyof the whole population – the mean ensemble frequency.Similarly, some of the neurons in the brain are expectedto be mutually locked to a certain frequency, whereas anelectroencephalographic recording contains the mean ofall neurons in some area, not just the synchronized.In general, there is no reason for these frequency def-initions to coincide. Still, not enough attention hasbeen paid in formulating them for different parameters.Namely, due to the equality of the frequencies in the caseswith symmetry that were mostly studied, they were usedinterchangeably and without verification, even when theydo differ (e.g. see [9]). However, here we consider sce-narios which most closely resemble the actual physicalor natural phenomena; in particular, models with asym-metrically distributed frequencies, phase shifted couplingfunction, or asymmetric couplings of opposite sign. Forthem we show that these frequencies always differ andhave non-trivial values. Hence, one should be extremelycautious when the measured frequency of a populationis interpreted and then compared with the theoreticalmodel.We begin with a formulation of the model and its groupdynamics parameters. Section III describes the station-ary solutions of the KM, whilst all possible scenarios withnon-trivial mean field and mean ensemble frequencies aredescribed in Section IV. The summary of the work andits implications are discussed in Section V.
II. FORMULATION
The KM consists of phase oscillators running at arbi-trary intrinsic frequencies and coupled through the sineof their phase differences. In the case of heterogenouscoupling strengths, the dynamics of the phase ˜ θ i of the i th oscillator has the form˙˜ θ i = ˜ ω i + K i N N X j =1 sin(˜ θ j − ˜ θ i ) , i = 1 , . . . , N. (1)Here, K i is the coupling strength of each oscillator andit is drawn from a probability distribution Γ( K ). Simi-larly, the natural frequencies ˜ ω i are randomly distributedaccording to some ˜ g (˜ ω ). The tildes for the frequencies,phases and their distributions later are intentionally usedfor reasons to be explained below. Without loss of gen-erality, we assume ˜ g (˜ ω ) to have a mean < ˜ ω i > centeredat 0. Hereafter this frame of reference will be called nat-ural. Kuramoto introduced a complex order parameterfor this model, defined as a centroid of the complex rep-resentation of the oscillators˜ z ( t ) ≡ r ( t )e i ˜ ψ ( t ) = 1 N N X j =1 e i ˜ θ j . (2)It characterizes the macroscopic behaviour of the oscilla-tors, with r and ˜ ψ being the amplitude and the phase ofthe mean field respectively. The former shows the levelof synchronization, while the latter gives the position ofthe peak in the distribution of phases. Generally, thelong-term values for both can depend on time. ApplyingEq. (2) to the governing equation (1), it is rewritten as˙˜ θ i = ˜ ω i − K i r sin(˜ θ i − ˜ ψ ) . (3)For infinitely large populations N → ∞ , a probabilitydistribution function (PDF) ˜ f (˜ θ, ˜ ω, K, t ) is defined, suchthat R + π − π ˜ f (˜ θ, ˜ ω, K, t ) d ˜ θ = g (˜ ω )Γ( K ) . Thus, the complexmean field Eq. (2) becomes˜ z = Z π − π Z + ∞−∞ Z + ∞−∞ e i ˜ θ ˜ f (˜ θ, ˜ ω, K, t ) d ˜ θ d ˜ ω dK. (4)For convenience, infinite limits in all further definite in-tegrals will be omitted.As a consequence of the conservation of the number ofoscillators the evolution of the density function is gov-erned by a continuity equation ∂ ˜ f∂t = − ∂∂θ { [˜ ω + K i (˜ z e − i ˜ θ − ˜ z ∗ e i ˜ θ )] ˜ f } . (5)Here the right hand side of Eq. (3) is used and the sinefunction is expressed with complex exponents. Moreover, since ˜ f is 2 π periodic in ˜ θ it allows a Fourier expansionand can be written as˜ f = ˜ g (˜ ω )Γ( K )2 π { ∞ X k =1 [ ˜ f k (˜ ω, K, t ) e ik ˜ θ + c.c. ] } , (6)where c.c. are the complex conjugates and ˜ f − k = ˜ f ∗ k .In the limit t → ∞ , the ensemble described by (5)might settle into a stationary state for some rotatingframe. We define: i) systems that have stationary so-lutions in some frame of reference i.e. the complex meanfield and the distribution of phases rotate uniformly andthey have a constant mean field after the initial transi-tions; and ii) systems that experience complex non-stablebehavior, i.e. a time-varying amplitude of the order pa-rameter r ( t ). These definitions might differ from usualdescriptions found elsewhere which regard as stationaryonly those solutions that are fixed in the natural refer-ence frame. In this work, our attention is focused on theensembles with stationary solutions as described by i).For the case of identically coupled oscillators with uni-modal and symmetric distributions of their natural fre-quencies, above the critical coupling a phase locking ofthe oscillators takes place around the peak of ˜ g (˜ ω ) wherethe density of the oscillators is highest [2]. As a conse-quence of the symmetry, the group dynamics is station-ary and both mean frequencies are equal to the mode of˜ g (˜ ω ), which in this case is also its mean value.Nevertheless, introducing multimodal or asymmetric˜ g (˜ ω ), or distributed K leads to much richer dynam-ics. Thus, for multimodal ˜ g (˜ ω ) bistabilities and stand-ing waves emerge [10–12]. A standing wave is a macro-scopic solution where neither ˜ f (˜ θ, ˜ ω, K, t ), nor ˜ z ( t ) arestationary in any rotating frame. This further impliesnon-stationarity of r ( t ). Standing waves are also ob-served in systems with symmetrical bimodal distributionof natural frequencies, with the exact result for bifurca-tions between different states given in [13].Travelling waves (TW), as another peculiar group be-haviour, have also been observed for different parameterranges within the same systems [10, 11]. In our analysis,a TW state is considered to be any solution character-ized by long-term stationarity of the mean field ampli-tude, whereas the frequency of the locking Ω differs fromthe mean of the natural frequencies. In other words, thelocking of synchronized oscillators is in a frame differ-ent from the natural. This also represents a stationarysolution according to its definition above.A recent study [9] shows the occurrence of TW in mod-els with positive and negative coupling strengths, andidentifies so called conformists and contrarians. Simi-larly, a synchronization around a frequency that is dif-ferent from the mean or the peak of the distributionwas reported for ensembles that have an asymmetric uni-modal distribution of natural frequencies [14]. It is alsoworth mentioning here the Kuramoto-Sakaguchi model[15], where the phase shift is introduced into the couplingfunction, to allow synchronization at a frequency differ-ent from the mean of the natural frequencies. Hence, italways leads to TW states. Additionally, the whole classof models of non-isochronous oscillators with constantshear can be reduced to this model [2, 16].Stationary solutions for fully symmetric populationshave macroscopic frequencies that are equal to < ˜ ω i > .In asymmetric scenarios on the other hand, the syn-chronized cluster experience non-trivial phase velocity[9, 14, 15]. Thereafter, the focus of this work is inter-acting phase oscillators whose coherent behavior is char-acterized by a TW state, as defined earlier. Having cer-tain asymmetries, either in the frequencies, the couplingparameters, or in the coupling function itself, is a nec-essary condition for occurrence of this state. As a con-sequence, the influence of the unsynchronized oscillatorsto the entrainment frequency does not vanish [9, 14, 15].Additionally, we show that these oscillators also cause themean ensemble frequency to have a non-trivial value. Itgenerally differs from the entrainment frequency, but alsofrom the mean and the mode of the natural frequencies.Before proceeding with the analysis of all cases withTW state, let us define the macroscopic frequencies thatwe have already discussed. The mean field frequency rep-resents the velocity of the mean phase ˜ ψ and is obtainedfrom the time derivative of the complex mean field (2)˙˜ ze − i ˜ ψ = ˙ r + i ˙˜ ψr = 1 N N X j =1 i ˙˜ θ j e i (˜ θ j − ˜ ψ ) . (7)Taking into account that r and ψ are both real, the sameis true for their time-derivatives, so from Eqs. (3) and (7)one finally obtains following evolutions of the amplitudeand the phase of the complex mean field˙˜ ψ ≡ Ω = 1 rN N X j =1 [˜ ω j − K j r sin(˜ θ j − ˜ ψ )] cos(˜ θ j − ˜ ψ ) . (8)The expression for ˙˜ ψ represents the velocity of the meanphase, i.e. the frequency of the synchronized oscillators.For the other frequency parameter – the mean frequencyof the ensemble – used for characterizing these systems,its definition leads to˜ f ens ≡ N N X j =1 ˙˜ θ j = 1 N N X j =1 [˜ ω j − K j r sin(˜ θ j − ˜ ψ )] . (9)In the infinite limit, (6) is introduced into (4). By takingthe time derivative and applying the substitution (5), theevolution of the complex order parameter is obtained˙˜ z = Z Z [ i ˜ ω ˜ f ∗ + K z − ˜ z ∗ ˜ f ∗ )]˜ g (˜ ω )Γ( K ) d ˜ ωdK. (10)Similarly, the mean frequency of the ensemble becomes˜ f ens = Z π − π Z Z ˙˜ θ ˜ f (˜ θ, ˜ ω, K, t ) d ˜ θd ˜ ωdK. (11) III. STATIONARY SOLUTIONS FOR THEPHASE DISTRIBUTION
For a large class of problems, the long term coher-ent dynamics of an ensemble of phase oscillators is time-independent, as a consequence of the stationary distribu-tion of the phases. This implies existence of the station-ary solution of the continuity equation (5), which doesnot need to be in the natural reference frame. A recentwork [17] discusses generalized empirical stability condi-tions for these systems.From now on we consider that the ensemble has non-zero mean field, i.e. it is out of the incoherent state(a fully incoherent solution is also stationary), and thatstationary solutions for the phase distribution exist insome planes of reference. This is true for the simplestpossible scenario – unimodal symmetric frequency distri-bution and constant coupling strengths. The TW state isanother possible scenario with this property, despite thefact that in this situation the velocity of the mean phase(or the mean field frequency as defined here) differs fromthe mean of the natural frequencies.In our analysis the mean field frequency is allowed tobe nonzero despite assuming that < ˜ ω > = 0. Still, wework in the frame where ˜ f (˜ θ, ˜ ω, K, t ) is stationary – thereference frame rotating with the frequency Ω. Here thephases of the oscillators are θ = ˜ θ − Ω t and the phasecorresponding to the complex order parameter is ψ = ˜ ψ − Ω t . The distribution of the natural frequencies becomes g ( ω ) = ˜ g ( ω +Ω) with mean − Ω. In the same frame ψ = 0can be assumed after an appropriate phase shift.For any ensemble, if the stationary solution of (5) ex-ists, then it exhibits two types of long-term behavior,depending on the size of | ˜ ω − Ω | = | ω | relative to | Kr | and to the sign of K . The oscillators with | ω | < | Kr | approach a stable fixed point defined implicitly by θ = (cid:26) arcsin ω | K | r , if K > π + arcsin ω | K | r , if K < locked because they maintainconstant phase difference, while rotating at frequency Ωin the original frame. It is also assumed that synchro-nized oscillators with positive couplings have phases inthe interval ( − π/ , π/ π/ , π/ | ω i | > | K i r | rotate in a non-uniform manner. As ex-pected, the locked oscillators correspond to the center of g ( ω ) and the drifting oscillators correspond to the tails.For the synchronized oscillators the stationary distri-bution of the phases becomes f s ( θ, ω, K, t ) = (cid:26) g ( ω )Γ( K ) δ [ θ − arcsin( ω | K | r )] , if K > g ( ω )Γ( K ) δ [ θ − arcsin( ω | K | r ) − π ] , if K < . (13)Oscillators with frequencies in the interval | ω | > r | K | areout of synchrony with the mean phase. Their stationarydistribution is obtained from the continuity equation (5)and the normalization condition of f ( θ, ω, K, t ), such that f as ( θ, ω, K, t ) = g ( ω )Γ( K ) p ω − ( Kr ) π | ω − Kr sin( θ ) | . (14)Hence, the real and imaginary parts of the complex meanfield definition, Eq. (4), become r = Z π − π Z Z cos θf ( θ, ω, K, t ) dθ dω dK, (15)0 = Z π − π Z Z sin θf ( θ, ω, K, t ) dθ dω dK. (16)The latter is identified as the phase balance equation.The distribution of the phases for synchronized oscil-lators, Eq. (13), for each ω implies π difference betweenthe phases of the synchronized clusters with couplings ofopposite sign. Nevertheless, this holds for the distribu-tion f s ( ω, K, θ, t ) only if g ( ω ) is symmetric. This is notthe case for the TW state, where even if the oscillatorswere symmetrically distributed in the natural frame, thesymmetry would be broken when moving to the framerotating with Ω. Thereafter, the centroids of the phaseswill be shifted from the π mutual distance, a phenomenonthat was reported in [9].In the rotating reference frame the locked oscillatorsare frozen, i.e. by definition ˙ θ = ˙ ψ = 0. Thus only thedrifting ones need to be considered. This can be appliedin the expression for the mean frequency of the ensemble,Eq. (11), which becomes (see Appendix A) f ens = Z Z ∞| Kr | [ g ( ω ) − g ( − ω )]Γ( K ) p ω − ( Kr ) dωdK. (17)The definition (17) is not restricted to any distributionsof the natural frequencies and couplings.Going back to the rotating frame where the original˜ g (˜ ω ) has zero mean, frequency parameters of the TWstate are characterized by Ω and ˜ f ens , where˜ f ens = f ens + Ω . As a consequence of (16) Ω can equal 0 only if g ( ω )is symmetric in the boundaries of the integral in (17)[14]. This means that for any zero centered distributionof the natural frequencies ˜ g (˜ ω ) which is even in theintervals ˜ ω > | rK | (e.g. Lorentzian, Gaussian, etc) thefunction inside the integral is also even if Ω = 0. Thusthe integrals cancel each other and ˜ f ens = f ens = Ω = 0. IV. THE TRAVELLING WAVE STATES
In the following we discuss the possible scenarios thatlead to the TW state for which we also obtain mean fre-quency parameters described earlier.All further analyses are carried out in the rotatingframe of the entrainment frequency Ω, such that ˙ ψ ≡ g ( ω + Ω), where ˜ g (˜ ω ) has zero mean. However, for betterclarity of the figures, frequencies in the examples withasymmetric ˜ g (˜ ω ) are depicted as seen from the originaldistribution described in the captions. In other words,compared to the results from the analysis, the plots areshifted by the means of the given distributions ˜ g (˜ ω ). A. Traveling waves in the Kuramoto model withcontrarians and conformists
First we focus on the TW solution described in [9],which actually inspired this work. The model shows re-semblance to sociophysical models of opinion formation[18] and is also a continuation of the KM with distributedpositive couplings [19].The distribution of the couplings isΓ( K ) = (1 − p ) δ ( K − K ) + pδ ( K − K ) , where K < K >
0, and p denotes the probabilitythat a randomly chosen oscillator is a conformist, while q = 1 − p is the probability that a random oscillatoris a contrarian. The natural frequencies follow a zerocentered Lorentzian distribution with half-width γ ˜ g (˜ ω ) = γπ (˜ ω + γ ) . As stated in [9], if the absolute coupling strength is higherfor conformists than for the contrarians, then for someregion in the parameters space γ − p , the synchronizedoscillators will experience a TW. This means that bothpeaks in the phase distribution uniformly rotate in samedirection. The waves appear in symmetric pairs, withfrequencies ± Ω, since they result from the asymmetryin the coupling strengths. On the contrary, if the TW isdue to the asymmetry in the natural frequencies or in thecoupling function, the waves are not paired, as discussedlater.Following the definitions of Γ( K ) and ˜ g (˜ ω ), and usingthe substitution ω − ( K / r ) = u , , the integrals inEq.(17) are analytically solved, yielding f ens = −√ γ n − p qp ( γ + Ω + K r ) − K r Ω − Ω + K r + γ ++ p qp ( γ + Ω + K r ) − K r Ω − Ω + K r + γ o . (18) | Ω | , | ˜ f e n s | p r FIG. 1. (color online) The amplitude r and frequency Ω ofthe mean field, and the mean ensemble frequency ˜ f ens versusthe ratio p . Theoretical results are given with a solid line(black) for r , a dotted line (red) for Ω, and a dashed line(blue) for ˜ f ens . Results from numerical simulations are shownwith squares (black) for r , triangles (red) for Ω, and diamonds(blue) for ˜ f ens . Parameters: γ = 0 . K = − K = 2. This expression can be straightforwardly generalized formultimodal- δ distributed coupling strengths.It is obvious that Ω = 0 will imply ˜ f ens = 0. Henceonly in the presence of a TW may the mean frequency ofthe ensemble differ from the mean phase velocity for thismodel. Similarly, for the TW state f ens is non-zero andhas opposite sign from Ω. Additionally it can be shownthat the expression in the curly brackets is smaller than1 /γ , so that in the natural reference frame, | ˜ f ens | < | Ω | always holds. This is also evident from the numericalresults plotted in Fig. 1.Let us now derive the expression for obtaining the fre-quency Ω of the TW, as seen in the natural referenceframe. Ott and Antonsen, in their seminal work [21],showed that macroscopic evolution of large systems ofcoupled oscillators can be described by an explicit defi-nite set of nonlinear differential equations. They intro-duced an ansatz for the complex Fourier coefficients inEq. (6) ˜ f n (˜ ω, K, t ) = [˜ α (˜ ω, K, t )] n , (19)which exactly solves the governing equation (5), as longas ˜ α (˜ ω, K, t ) evolves following the nonlinear equation ∂ ˜ α∂t + i ˜ ω ˜ α + K z ˜ α − ˜ z ∗ ) = 0 . (20)When this ansatz is implemented in Eq.(4), the orderparameter reduces to˜ z = Z Z ˜ α ∗ ( ω, K, t )˜ g (˜ ω )Γ( K ) d ˜ ωdK. (21) Note that applying this ansatz in (10), leads to a sim-plified expression for the evolution of the complex orderparameter˙˜ z = Z Z [ i ˜ ω ˜ α ∗ + K z − ˜ z ∗ ˜ α ∗ )]˜ g (˜ ω )Γ( K ) d ˜ ωdK, (22)from where the frequency of entrainment can be ob-tained. Similarly, substituting Eqs. (3) and (6) and theansatz (19) into Eq. (11) transform the mean ensemblefrequency to˜ f ens = Z Z [˜ ω − K i (˜ zα − ˜ z ∗ ˜ α ∗ )]˜ g (˜ ω )Γ( K ) d ˜ ωdK. (23)We start by considering the low-dimensional evolution(20) in the reference frame of the TW. Similar analysis,but in the natural frame was performed in [9]. Firstthe bimodal- δ and Lorentzian distributions for Γ( K )and g ( ω ) respectively, are substituted into the integralEq. (21), such that it yields z ∗ = (1 − p ) α (Ω − iγ, K ) + pα (Ω − iγ, K ) . (24)This is then used in Eq. (20), which rewritten for both, α (Ω − iγ, K ) and α (Ω − iγ, K ), results in ∂α , ∂t +( i Ω + γ ) α , + K , zα , − z ∗ ) = 0 , (25)where we omit the dependencies of α ( ω, K, t ). As pre-viously stated, we are interested in a stationary solu-tion of the TW state in the t → ∞ limit. This im-plies time-independent distribution of the phases in thislimit, or from the ansatz (19), time-independent α , with ∂α , /∂t = 0 . Further, similar to [9], complex order pa-rameters for each of the subpopulations, and the differ-ence between their phases are defined r e − iψ = α , r e − iψ = α , δ = ψ − ψ = const . (26)In this way it is ensured that both synchronized pop-ulations, in-phase and antiphase, rotate with the samevelocity Ω in the natural frame, and preserve constantphase difference. In the rotating frame of the TW ˙ ψ = 0,so ψ ≡ { r , r , ψ , ψ } space.˙ r = − γr − K r − pr cos δ + qr )] = 0 , (27)˙ r = − γr − K r − pr + qr cos δ )] = 0 , (28)˙ ψ = Ω − K r pr sin δ ( r + 1) = 0 , (29)˙ ψ = Ω + K r qr sin δ ( r + 1) = 0 . (30)The low-dimensional parameters including Ω can be nowobtained self-consistently. The steady states (29-30) re-sult from the constant angle difference between the peaksin the phase distribution, i.e.˙ δ = − sin δ [ K r pr ( r + 1) + K r qr ( r + 1)] = 0 . Hence, when Im [ z ] = 0 is applied to Eq. (24), one canuse the expression qr sin ψ = − pr sin ψ (31)and the definition (26) of δ to obtain the values of ψ and ψ , such that the system will be fully described.The equations (27-30) and all further numerical inte-grations are performed using a Runge-Kutta 4th orderalgorithm. Once we have the low-dimensional parame-ters, Eq. (18) is applied to find the mean frequency of theensemble. Finally, results are compared with the valuesfor the order parameter r , the mean phase velocity of theensemble Ω and the mean frequency f ens , obtained fromthe numerical simulations of the ensemble Eq. (1), usingEq. (2) and Eq. (9). As in all later simulated scenarios,the number of oscillators was set to N = 100000, thetime step of the integration was 0 .
01. The simulationswere running for 10 time steps, with the initial 90% ofeach run discarded as possibly transient, while the restwere time averaged. The proportion of the conformists p is changed from 0 to 1 at 100 equally spaced points,and the obtained results given in Fig. 1 fully confirm thetheoretical analysis. B. Asymmetric unimodal frequency distribution
Another case of the KM characterized with a station-ary solution, where the mean frequency of locked oscilla-tors differs from the mean of the oscillators’ natural fre-quencies, are ensembles with asymmetric unimodal dis-tribution of the natural frequencies, and equal couplings.The asymmetric scenario is also more natural than thesymmetric, because any imperfection in the system, how-ever small, can destroy the ideal symmetry. Nevertheless,mostly due to the analytical difficulties, this case has ob-tained little attention. It was first examined by Sakaguchiand Kuramoto [15] and the self-consistent condition forthe mean field frequency was obtained in [14].It is known from [14], that the phase balance equation(16) for g ( ω ) asymmetric in the interval | ω | > rK , impliesthat the oscillators always lock to a frequency that differsfrom the mean of the natural frequencies. For the meanensemble frequency Eq. (17) was numerically integratedfor ensembles with triangular and lognormal frequencydistributions and constant coupling strength. We foundthat it always gives f ens = − Ω in the TW rotating frame.In other words, the mean ensemble frequency equals themean of the natural frequencies, although it differs fromthe mean field frequency. The results seemed unexpectedat first sight. However after careful analysis we realized Ω , ˜ f e n s c (a) r Ω , ˜ f e n s K (b) r Ω , ˜ f e n s σ (c) r FIG. 2. (color online) The amplitude r and frequency Ω ofthe mean field, and the mean ensemble frequency ˜ f ens forunimodal asymmetric natural frequencies. Results from thenumerical simulations are shown with squares (black) for r ,triangles (red) for Ω and diamonds (blue) for ˜ f ens . The lines(blue) are from theoretical results for ˜ f ens and they also matchthe means of the frequencies’ distributions – crosses (blue).(a-b) Triangular ˜ g (˜ ω ) in boundaries a = 1 and b = 10. (a)Coupling K = 4 . c ∈ [0 , K ∈ [4 ,
10] and c = 9; (c) log-normal distribution of natural frequencies, with µ = 0, σ ∈ [0 . ,
2] and coupling K = 4 .
2. The dashed lineshows the modes of the distributions. that the derivations for the equation of balance Eq. (16),as given in [14], in the frame rotating with Ω indeed leadsto the Eq. (17). Namely, applying the derivations in theAppendix A, Eq. (16) becomes0 = Ω +
Z Z ∞| Kr | g ( ω ) p ω − ( Kr ) dω − Z Z −| Kr |−∞ g ( ω ) p ω − ( Kr ) dω = Ω + f ens . (32)Thus, one can conclude that for the KM with asymmetricfrequency distribution, the mean ensemble frequency isalways equal to the mean of the natural frequencies.Next, a triangular distribution with limits a and b , andpeak at c was explored. In the scenario shown in Fig. 2(a) the peak is distributed in the interval [ a, b ], while inFig. 2 (b) the coupling strength is increasing from 4 to10 for fixed c . Similarly, for log-normal frequencies˜ g (˜ ω ) = 1˜ ωσ √ π e − (ln ˜ ω − µ ) / (2 σ ) , ˜ ω > ,µ is fixed to 0, while σ is logarithmically distributed, andthe results are given in Fig. 2 (b). The mean ensem-ble frequency values also match theoretical values for themeans of the given distributions, ( a + b + c ) / e µ + σ / respectively.For the triangular distribution shown in Fig. 2 (a) themode is c , while for (b) is fixed. So for the log-normal˜ g (˜ ω ) given in Fig. 2 (b) the mode is 1 for µ = 0 and σ → e µ − σ as shown. Hence, the presented results show that themean field frequency is always between the mode andthe mean of the distribution of natural frequencies, andreaches the second only when all oscillators become syn-chronized. That is to say, by increasing the couplingstrength for a given frequency distribution, the propor-tion of synchronized oscillators is also increased. Accord-ingly the value of Ω moves closer to the mean of the dis-tribution, until it eventually reaches it for r →
1. Forunbounded g ( ω ) the last can only occur when K → ∞ ,while in the case of bounded natural frequencies, for somevalue of rK all oscillators will be entrained. This canalso be deduced from Eq. (32). Namely, higher synchro-nization implies a smaller region for the integral on theright hand side. At the same time we assume g ( ω ) to bedecreasing left and right from the mode, meaning thathigher r leads to smaller value of f ens which eventuallybecomes 0, implying Ω = ˜ f ens . Furthermore, the valueunder the square root in the same integral also decreaseswith increasing r . These are confirmed in Fig. 2 (b-c),where we see that for smaller couplings and hence forsmaller mean field amplitudes, Ω is closer to the peak ofthe distribution and approaches ˜ f ens for larger coupling. C. Asymmetric multimodal frequency distribution
The KM with multimodal asymmetric distribution ofnatural frequencies is another candidate for the meanfield behavior described by a TW. Nevertheless, due todifficulties that arise in the mathematical analysis of thismodel, it was never fully solved, nor has a thorough dy-namical analysis of possible macroscopic solutions beenperformed. Still, following the analysis of the symmetricscenario [13] and qualitative descriptions of the dynamicsin the asymmetric case given in [2] and [13], some conclu-sions can be drawn. Namely, Kuramoto, in his seminalwork [2], discusses how transition from incoherence tomutual synchronization might be modified when the os-cillators’ natural frequencies are bimodally distributed.For sufficiently large coupling strength, he assumed thatthe clusters of synchronized oscillators “will eventually beentrained to each other to form a single giant oscillator”. For smaller coupling strengths compared to the distancebetween peaks, he envisaged that the synchronized nu-clei would be at the peaks of g ( ω ). Although some ofthe transitions between different states described in [2]for the symmetric case were shown to be wrong [13], thedescription of the mean field dynamics of a partially syn-chronized state is indeed correct.Hence, for a symmetrical g ( ω ) the ensemble should beeither partially synchronized for large enough couplingcompared to the peaks’ distance, or synchronized clustersshould exist near the peaks [2, 13, 22]. For asymmetri-cally and bimodally distributed frequencies the partialsynchronization will be characterized by TW, whereasstanding waves and other complex collective rhythms ap-pear for smaller couplings.Thereafter, for the partially synchronized state in theasymmetric bimodal case, one might expect a TW stateto occur. As a consequence of having a single synchro-nized cluster, the dynamic of the system is similar to themodels with unimodal asymmetric distribution describedin Section IV B. As a result the same conclusion drawnabout the asymmetric unimodal scenario, also holds forthe TW solution of the asymmetric bimodal case. Thismeans that the mean frequency of the ensemble is equalto the mean of ˜ g (˜ ω ), and as such differs from the fre-quency of the TW experienced by the mean field. Hence,in the TW frame, f ens = − Ω. This can be seen in Fig. 3,where results from numerical simulations confirmed thetheoretically expected values for ˜ f ens and Ω.As an example we analyze the dynamics of a popula-tion with bimodal Lorentzian distribution of frequencies˜ g (˜ ω ) = (1 − p ) γ π [(˜ ω − µ ) + γ ] + pγ π [(˜ ω − µ ) + γ ] . The subdistributions are peaked at µ and µ with thehalf-widths γ and γ respectively, where p is the pro-portion of the oscillators belonging to the second sub-distribution. The full low-dimensional dynamics for thissystem can be obtained using the ansatz Eq. (19) in asimilar manner as in Section IV A. The integral (21) forthe given Γ( K ) and g ( ω ) yields z ∗ = (1 − p ) α ( µ − γ , K, t ) + pα ( µ − γ , K, t ) . Hence, substituting α , in Eq. (20), the low-dimensionaldynamics are described by two ODEs˙ α , = − ( i µ , + γ , ) α , − K z ∗ + zα , ) . (33)In all the plots in Fig. 3 the values of the couplingstrength are large enough to induce TW instead of stand-ing wave states, i.e. one instead of two synchronizedclusters. For the Fig. 3 (a) both peaks are of equaldistance from the zero. This leads to zero mean for anywidth of the subdistributions and following previous dis-cussion and Eq. (32) it implies f ens = 0. For similar devi-ations of subdistributions the entrainment frequency willbe near zero, but closer to the narrower peak, since more −0.20.20.61 Ω , ˜ f e n s (a) r Ω , ˜ f e n s (b) r Ω , ˜ f e n s (c) r Ω , ˜ f e n s γ (d) r FIG. 3. (color online) The amplitude r and frequency Ω of themean field, and the mean ensemble frequency ˜ f ens versus thefrequency width of the first subpopulation γ . (a-d) Resultsfrom the simulations, Eq. (1), are shown with squares (black)for r , triangles (red) for Ω and diamonds (blue) for ˜ f ens . (a-c)The low-dimensional dynamics, Eq. (33), are shown with solidlines (black) for r , dotted lines (red) for Ω and dashed lines(blue) for ˜ f ens . (a-c) Bimodal Lorentzian ˜ g (˜ ω ). Parameters: γ ∈ [0 , γ = 0 . µ = − µ = 1, K = 5, and (a) p = 0 .
5, (b) p = 0 .
25 and (c) p = 0 .
75. (d) Bimodal Gaussian˜ g (˜ ω ) with γ ∈ [0 , γ = 0 . µ = − . µ = 1, coupling K = 4 .
25 and p = 0 . oscillators from that subdistribution will be entrained.However, for a large deviation of the first subdistribu-tion compared to the second one, a small number of theoscillators belonging to it will be entrained. These oscil-lators will be almost equally distributed on both sides ofthe second peak. Thus, the frequency of the entrained oscillators will asymptotically reach the peak of the sec-ond subdistribution in the limit case γ → ∞ . If in thesame limit case p → p = 0 .
25 and p = 0 .
75 as inFig. 3 (b-c), the mean of the frequencies is on 1 / / g (˜ ω ). D. Contrarians and conformists with asymmetricunimodal frequency distribution
The analysis naturally continues with the KM in thepresence of both previously analyzed conditions that arerequired for TW – contrarians and conformists, withasymmetric distribution of natural frequencies. Thephase locking in this scenario happens in frames otherthan the natural for all coherent solutions. But as it willbe shown, due to the distributed couplings, the mean en-semble frequency also has non-trivial values – it is notalways equal to < ˜ ω > . As another consequence of theasymmetry, the sign of the TW can no longer be expectedto be random, as is the case for the contrarians and con-formists over symmetrically distributed frequencies [9].However, the analysis of this case is largely more com-plicated compared to previous ones and we have takensome examples only to show the main points.The phase balance equation (16) for distributed K inthis case does not lead to a simple expression as wasEq. (32) for equal couplings. Still, following the sameprocedure, Eqs. (13, 14) are substituted into Eq. (16) andapplying the derivations from Appendix A, the balanceof phases becomes0 = Z Γ( K ) rK dK { Z ωg ( ω ) dω − Z ∞| Kr | ω [ g ( ω ) − g ( − ω )] p ω − ( Kr ) dω } = Z Γ( K ) rK dK [ − Ω + I ( K )] . (34)This expression self-consistently gives the frequency ofthe entrainment, whilst it no longer implies equality of˜ f ens and < ˜ ω > , as for constant coupling parameters.Results from numerical simulations shown in Fig. 4 con-firm the non-trivial nature of both global frequency pa-rameters. The most interesting phenomena is that Ωand ˜ f ens for some parameters can have values that notonly differ from either the mode and the mean of thenatural frequencies, but that are also outside of the in-terval between them as was case for equal K . Moreover,they can even be outside of the boundaries of their dis-tribution. For example for the triangular ˜ g (˜ ω ) ∈ [0 , ∼ − ∼
2, and similar happens with˜ f ens . These values are clearly out of the region of supportof ˜ ω , which is shaded in the same plot. Note that theseon first sight counterintuitive results are a consequenceof the very high coupling strengths with opposite sign,and they are still within the general boundaries of ˜ f ens and Ω. Namely, from Eqs. (11) and (3), it is clear thatˆ ω − max | K | < ˜ f ens < ˆ ω + max | K | , (35)while from the fact that the entrainment frequency can-not be out of the limits [min( ˙˜ θ ) , max( ˙˜ θ )] it followsmin(˜ ω ) − max | K | < Ω < max(˜ ω ) + max | K | . (36)Of course, these are broad limits and only for bounded˜ g (˜ ω ) do the boundaries for Ω not reach ±∞ .Hence, for bounded distributions, as the triangular,and for large enough couplings of opposite signs, the onlyway the balance equation (34) can be satisfied, is with theemergence of a TW with a large enough value, so that theintegral I ( K ) in the same equation will be non-vanishing.Contrary, if I ( K ) = 0 for the shown example, i.e.˜ ω − Ω < | rK | , ∀ ˜ ω for which ˜ g (˜ ω ) > , (37)then the phase balance becomes0 = (1 − p ) /K + p/K . (38)If (38) holds, then Ω can have any value that still obeyscondition (37), so each simulation would give a differentTW within these limits. Another consequence from con-dition (37) when applied to Eq. (17) is that f ens = 0.If (38) is not satisfied, then the oscillators rearrange,rendering the condition (37) also invalid. The result isthe appearance of non-zero I ( K ) in Eq. (34) to imposethe balance of the phases. There may be more than onevalue for Ω which impose this balance and if that was thecase one might expect multiple solutions or even hystere-sis behavior. In this work however, we only numericallyconfirm that for bounded natural frequencies, as in theexample in Fig. 4 (a), Ω, when observed in the naturalframe, can have two stable values with opposite signs. Inother words, numerical realizations could lead to any ofthe two stable values. This is illustrated in Fig. 4 (a)where one such realization is depicted, whilst in differentrealizations the values on both sides of the mean and themode of ˜ g (˜ ω ) appeared for p around (0 . , . p > .
65, then a state similar to the π state in the symmetricmodel in Section IV A is observed. Here, the equalityof Ω and ˜ f ens follows from the bounded distribution of Ω , ˜ f e n s p (a) r Ω , ˜ f e n s p (b) r −6 0 6 120.511.52 Ω , ˜ f e n s K (c) r FIG. 4. (color online) The amplitude r and frequency Ω of themean field, and the mean ensemble frequency ˜ f ens for asym-metric ensembles with distributed coupling strengths. Nu-merically obtained results, Eq. (1), are shown with squares(black) for r , triangles (red) for Ω and diamonds (blue) for˜ f ens . (a-c) The dashed lines (blue) are theoretically predictedresults for ˜ f ens , Eq. (17). Horizontal dashed and dotted linesare the modes and the means of ˜ g (˜ ω ) respectively, and its do-main is shaded for (a). (a) Triangular ˜ g (˜ ω ) within a = 0 and b = 1 and peak at c = 0 .
8; bimodal- δ Γ( K ) with K = − K = 8 and p ∈ [0 . , g (˜ ω ) with µ = 0and bimodal- δ Γ( K ). (b) σ = 1 . K = − K = 20; and(c) σ = 1, p = 0 . K ∈ [ − , K = 6. the oscillators’ natural frequencies, which all become en-trained. Similarly the integral I ( K ) vanishes, and therequirement for the phase balance holds only if Ω = 0.In this way only the first integral of the equation (34) sur-vives and it gives < ˜ ω > , which we set to be 0 (althoughin Fig. 4 the plots are for non-zero mean frequency dis-tributions, such that the plot (a) has mean 0 . I ( K ) =0, meaning that a TW always emerges. For the simula-tions we have performed, as seen in Figs. 4 (b-c), the sys-tem was always setting the same value for Ω that solvesEq. (34) for the given parameter range. Thus, for the log-normally distributed ω , Ω is always positive, although the0existence of another stable TW that also fulfills the phasebalance is not excluded. The entrainment of all oscilla-tors is never achieved in this scenario, and the mean fre-quencies will never reach < ˜ ω > . Still they will approachthis value when the number of entrained oscillators is in-creased, either by increasing the number of conformists,Fig. 4 (b), or by changing the coupling strength of one ofthe groups, as shown in Fig. 4 (c). The same plots alsoshow that even when the couplings are all positive, thegroup frequencies behave qualitatively differently fromthe case with constant K . This characteristic becomesless obvious once the modes are closer to each other. Tocheck the generality of the analysis, we also explored thecase when K is bimodal Gaussian distributed. The re-sults did not show any qualitative difference from thecase with bimodal- δ .It is expected that for multimodal couplings, cluster-ing emerges and plays a crucial role in defining the valuesof ˜ f ens and Ω. A similar phenomenon was also reportedin [19] for distributed positive and unimodal K , but asa finite-size effect only. That is to say, the distributionof the locked oscillators Eq. (13) for multimodal- δ Γ( K )will have peaks on different angles ψ i = arcsin ωK i r foreach mode K i . Further evidence for this explanation canbe seen in numerical simulations performed over differentΓ( K ) and g ( ω ). Fig. 5 presents obtained phase distribu-tions. As can be seen three qualitatively different typeswere identified.Firstly for very large coupling strengths compared tothe width of g ( ω ), and for bimodal Γ( K ) with modeshaving different signs, there are two separated peakscorresponding to the contrarians and conformists. Thedistance between them is smaller than π and the ob-served TW is similar as in Section IV A, although theasymmetry in g ( ω ) causes asymmetry in the distribution f ( θ, ω, K, t ). This is shown in Fig. 5 (a), and also in ex-amples (a-b) of Fig. 4, in the range for p around 0 . − . π statefrom [9], but the form of g ( ω ) will be mapped onto thepeaks which are separated by π . The asymmetry addi-tionally will influence the values of the ˜ f ens and Ω. How-ever, the latter will always be in the interval betweenthe mode and the mean of ˜ g (˜ ω ). This means that qual-itatively this corresponds to the case with constant K described in Section IV B, but with non-zero values for˜ f ens because of the phase balance equation (34).Lastly, if the couplings distribution have positive peaksonly, depending on their distance apart, the phases willhave either one or two close peaks, while still keeping theform of g ( ω ). Hence, the case with one peak will be thesame as having unimodally distributed K and numeri-cal results suggest that ˜ f ens ≈ h ˜ ω i as for constant K ,although this is not immediately obvious from Eqs. (17)and (34). Finally, if f ( θ, ω, K, t ) is no longer unimodal, f ( θ , ω , K , t ) − π π (a) θ f ( θ , ω , K , t ) − π π (b) θ f ( θ , ω , K , t ) − π π (c) θ f ( θ , ω , K , t ) − π π (d) θ FIG. 5. (color online) PDF of phases for different states ofthe ensemble with log-normally distributed natural frequen-cies and bimodal- δ distributed coupling strengths. Oscillatorswith coupling K are shown with a dashed line (red), thosewith K are following the dotted line (blue,) while the jointdistribution is shown with a solid line (black). Parameters: µ = 0, σ = 1 and p = 0 .
5. (a) K = − K = 20. Two peaksat a distance smaller than π and a TW occurs. (b) K = − K = 8. Two peaks at a distance π . (c) K = 2 K = 3. Thepeaks merge into a single peak in the joint distribution. (d) K = 2 K = 6. Two closely separated peaks. then ˜ f ens has values that differ from < ˜ ω > .The number and the nature of the emerging clustersis not so straightforward for non- δ Γ( K ) distributions.Still, from the PDF of the phases it is clear that for neg-ative modes in Γ( K ), the peaks in the PDF correspondingto those oscillators do not have the form of g ( ω ), but itis more symmetric. For positive modes, the peaks keepmapping the distribution of the natural frequencies. Thisis also an interesting peculiarity that requires further at-tention.The bottom line is that different and on first sightcounterintuitive values for the entrainment and mean en-semble frequency are obtained even for all positive butmultimodally distributed couplings, as long as the natu-ral frequencies are nonsymmetric. Nevertheless, a thor-ough analysis is still needed for this general case of theKM, with the emphasis put on the entrainment at differ-ent clusters and their influence to the group behavior. E. Phase shifted coupling function
The final model leading to TW of the macroscopic pa-rameters being analyzed here is actually the earliest gen-eralization of the KM [15]. This model was introduced toallow for entrainment in frames other than the natural.The so called Kuramoto-Sakaguchi model has the form˙˜ θ = ˜ ω − K r sin(˜ θ i − ˜ ψ − β ) , where β is a phase shift of the coupling function. Low-dimensional dynamics of this model with Lorentzian ˜ g (˜ ω )1
0 −4−2024 Ω , ˜ f e n s β − π / π / r FIG. 6. (color online) The mean field’s amplitude r , the phasevelocity Ω, and the mean frequency of the ensemble ˜ f ens ver-sus the phase shift β . Results from the theoretical predictionsfor the low-dimensional dynamics, Eqs. (39 - 40), are givenwith a solid line (black) for r , a dotted line (red) for Ω anda dashed line (blue) for ˜ f ens . Results from numerical simula-tions, Eq. (1), of the ensemble are shown with squares (black)for r , triangles (red) for Ω and diamonds (blue) for ˜ f ens . Pa-rameters: γ = 0 . K = 6 and α ∈ [ − π/ , π/ can be also obtained using the ansatz [21], with deriva-tions being similar as for the problem in Section IV A. Ityields ∂ ˜ z∂t + i (ˆ ω − iγ )˜ z + K z ˜ α e iβ − ˜ z ∗ e − iβ ) = 0 , which was also obtained in [20]. From here, one imme-diately obtains the long-term evolution of ψ , i.e. thefrequency of the TW d ˜ ψdt ≡ Ω = ˆ ω + K sin β − γ tan β. (39)For the mean ensemble frequency in this model, Eq. (17)can be analytically solved and with further substitutingstationary values for r and Ω, it transforms to˜ f ens = ˆ ω − sign (Ω) γ tan β. This directly leads to˜ f ens = ˆ ω + K sin β − γ tan β, (40)which can also be straightforwardly obtained fromEq. (23) after applying the residue theorem for the in-tegral over ω . The stationary values of the macroscopicparameters of this system are illustrated in Fig. 6. There,it is also obvious that | ˜ f ens | < | Ω | and that the frequen-cies have odd symmetry, features which are a direct con-sequence of Eqs. (39) and (40).It is worth noting that by taking the derivative over β from expressions (39, 40), extreme values for both meanfrequencies follow as β max/min = ± arccos p γ/K, for Ω and β max/min = ± arccos p γ/K, for ˜ f ens . Of course, since the population should be coher-ent, parameters have to be chosen such that K > K c =2 γ/ cos β holds. Thus, the values of the phase shift β thatproduce maximum deviation in both mean ensemble andmean field frequencies are directly obtained. V. SUMMARY
The TW state in the KM, occurs whenever the symme-try in either the natural frequencies or the coupling func-tion itself is broken. For asymmetric coupling strengthson the other hand, the wave occurs only for certain cou-pling parameters. The results of this study indicate thatin populations with an asymmetric, bell-shaped distribu-tion of the non-synchronized oscillators and equal cou-plings, the mean ensemble frequency is always the meanof the natural frequencies. In contrast, the TW of thesynchronized oscillators – the mean field frequency – hasa value between the mean and the mode of those. Whenthe asymmetry originates from the coupling strengths orthe coupling function, both mean frequencies have non-trivial values. In particular they differ from the mean orthe modes of the frequency distribution.The case of the asymmetric unimodal frequency dis-tribution has quite a straightforward explanation. Herethe frequency of synchronization is a result of the inter-play between the cluster of the locked oscillators and thewhole ensemble. Synchronized oscillators tend for fre-quency locking at the mode, where the density of similaroscillators is the highest, whilst the influence of all oscil-lators would be balanced if it is at the mean of all natu-ral frequencies. The compromise is achieved through theself-consistent system for the phase balance and Eq. (32)follows directly from this interplay. This balance wouldbe destroyed if the nucleus is either on the peak or onthe mean.By increasing the number of locked oscillators, themean field frequency shifts towards the less skewed sideof the distribution, i.e. towards the mean, because thenumber and the influence of the drifting oscillators de-creases. Finally, once all the oscillators get locked, thefrequency of the entrainment will be exactly the mean ofthe distribution. In this way, the influence of all oscilla-tors is equal. The same explanation can be applied forbimodal or multimodal cases, since again there is onlyone cluster of synchronized oscillators in the TW state.On the other hand, in the case of contrarians and con-formists, the physical explanation of the mean ensemblefrequency, the frequency of entrainment and their dif-ference, is not yet very clear. Previous work went asfar as showing that in some parameter spaces for oppo-site sign couplings, both peaks in the distribution of thephases lose their stability. Instead, they start to chaseeach other and are no longer a distance of π apart, henceproducing the TW. Here we have emphasized that thisbehavior simply acts to preserve the phase balance, whichdirectly gives the value of the entrainment frequency,2while its influence on arranging the non-synchronized os-cillators also makes it responsible for the mean ensemblefrequency. Despite this, we believe that even by show-ing, analytically and numerically, that the velocity of themean phase and the frequency of the ensemble differ inthis case, an important characteristic that seems to beunique for TW states only is immediately revealed. Thisexplanation was missing when the model was first in-troduced [9], and clear distinction between macroscopicfrequencies was not made. Moreover, “the mean phasevelocity” obtained through the low-dimensional dynam-ics was defined as the mean ensemble frequency.In the more general and realistic situation, with asym-metric frequencies and distributed couplings, the macro-scopic frequencies can have values in very unexpectedboundaries – even outside the limits of the natural fre-quencies. The reason for this complex behavior is againa result of the system’s self-arranging, such that the in-fluence of the non-synchronized oscillators to the meanfield amplitude vanishes and the phase balance is main-tained. In order for this to be achieved for multimodallydistributed couplings, there are different clusters of syn-chronized oscillators for each mode. This case also showspossibilities for further research either in eventual hys-teresis dynamics and seeking the number of stable so-lutions, or in the clustering phenomena for distributedcoupling strengths.When the asymmetry is induced through the couplingfunction itself the macroscopic frequencies always differ,in a similar way to the case when it stems from the cou-pling parameters. In that sense the mean ensemble fre-quency is of the same sign, but with lower absolute valuethan the mean field frequency. Also, their dependance onthe phase shift is nonlinear, and the ways they respondto the shift do not coincide, with both of them havingthe highest values for different phase shifts.Taken together, these results suggest that wheneverthe population is experiencing a TW, the locking of theoscillators is at a frequency different to the mean of theinstantaneous frequencies. Since in inverse problems orin experiments it is often only the macroscopic param-eters that can be obtained, a clear interpretation of theobserved mean frequency is always needed. Moreover, the asymmetric scenarios tend to be far more abundantin the real physical systems. That is to say, asymme-try of the natural frequencies can immediately make anyscenario more realistic, whilst few of the examples withopposite coupling strengths can be traced in inhibitoryand excitatory neurons [23] or in social dynamics [18].As for the phase-shifted coupling function, although itwas introduced to describe the formation of nonlinearwaves in non-oscillatory media [15], it can be used tomodel various phenomena, such as Josephson junctions[3], mammalian intestines and heart cells [24]. Hence, inthe models based on measurements, like those describ-ing the brain dynamics [4], one should precisely define towhich of the macroscopic frequencies, of the mean fieldor mean ensemble, the measured frequency corresponds.Finally, this work reveals the strong need for future re-search that should explain the physical link between theobserved mean frequency of any system with cooperativedynamics and the two macroscopic frequency parametersdescribed here for non-trivial cases. ACKNOWLEDGMENTS
The study was supported by the EPSRC, SP is sup-ported by the Lancaster University PhD grant and LBwas supported by Faculty of Computer Science and En-gineering at the SS.Cyril and Methodius University andERC Grant
Appendix A: General formula for the meanensemble frequency
We work in the reference frame rotating with Ω = ˙˜ ψ .In this frame, the locked oscillators are frozen and onehas to consider the drifting ones only, Eq. (14). This canbe applied in the expression for the mean frequency ofthe ensemble, Eq. (11), which becomes f ens = Z π − π Z Z | ω | > | Kr | ( ω − Kr sin θ ) g ( ω )Γ( K ) p ( ω ) − ( Kr ) π | ω − Kr sin( θ ) | dωdKdθ = I − I . (A1)The first integral is I = Z Z | ω | > | Kr | ω Γ( K ) g ( ω ) dKdω, (A2)because of the probability normalization Z π p ( ω ) − ( Kr ) π | ω − Kr sin θ | dθ = 1 . (A3) The second integral requires more calculation. Integrat-ing over the phases first, for positive frequencies ω > | Kr | it yields Z π p ω − ( Kr ) Kr sin θ π | ω − Kr sin θ | dθ = ω − p ω − ( Kr ) , (A4)3while for negative frequencies ω < −| Kr | Z π p ω − ( Kr ) Kr sin θ π | ω − Kr sin θ | dθ = ω + p ω − ( Kr ) . (A5)It is interesting to note that both integrals are even in K .Following this, the second integral, I , simply becomes I = Z ∞−∞ Z −| Kr |−∞ g ( ω )Γ( K )[ ω + p ω − ( Kr ) ] dωdK + Z ∞−∞ Z ∞| Kr | g ( ω )Γ( K )[ ω − p ω − ( Kr ) ] dωdK. (A6) After partial cancelation of I with I because of (A1)one obtains a simple formula for calculation of the meanensemble frequency in the reference frame of the orderparameter f ens = Z ∞| Kr | Z ∞−∞ g ( ω )Γ( K ) p ω − ( Kr ) dωdK − Z −| Kr |−∞ Z ∞−∞ g ( ω )Γ( K ) p ω − ( Kr ) dωdK = Z ∞| Kr | Z ∞−∞ [ g ( ω ) − g ( − ω )]Γ( K ) p ω − ( Kr ) dωdK. (A7)Note that g ( ω ) is centered in − Ω and hence is asymmetric, even for symmetric ˜ g (˜ ω ) which is often the case. Returningto the original frequencies ˜ ω , the last expression becomes f ens = Z ∞| Kr | Z ∞−∞ [˜ g ( ω + Ω) − ˜ g ( − ω + Ω)]Γ( K ) p ω − ( Kr ) dωdK. (A8) [1] S. Strogatz, Sync: The Emerging Science of SpontaneousOrder (Hyperion, New York, 2003).[2] Y. Kuramoto, in
International Symposium on Mathemat-ical Problems in Theoretical Physics , edited by H. Araki,Lecture Notes in Physics Vol. 39 (Springer, New York,1975), Y. Kuramoto,