Instantaneous frequencies in the Kuramoto model
IInstantaneous frequencies in the Kuramoto model
Julio D. da Fonseca, Edson D. Leonel, and Hugues Chaté
2, 3, 4 Departamento de Física, Universidade Estadual Paulista, Bela Vista, 13506-900 Rio Claro, SP, Brazil Service de Physique de l’Etat Condensé, CEA, CNRS Université Paris-Saclay, CEA-Saclay, 91191 Gif-sur-Yvette, France Computational Science Research Center, Beijing 100094, China Sorbonne Université, CNRS, Laboratoire de Physique Théorique de la Matière Condensée, 75005 Paris, France (Dated: August 11, 2020)Using the main results of the Kuramoto theory of globally coupled phase oscillators combinedwith methods from probability and generalized function theory in a geometric analysis, we extendKuramoto’s results and obtain a mathematical description of the instantaneous frequency (phase-velocity) distribution. Our result is validated against numerical simulations, and we illustrate itin cases where the natural frequencies have normal and Beta distributions. In both cases, we varythe coupling strength and compare systematically the distribution of time-averaged frequencies (aknown result of Kuramoto theory) to that of instantaneous frequencies, focussing on their qualitativedifferences near the synchronized frequency and in their tails. For a class of natural frequency distri-butions with power-law tails, which includes the Cauchy-Lorentz distribution, we analyze rare eventsby means of an asymptotic formula obtained from a power series expansion of the instantaneousfrequency distribution.
I. INTRODUCTION
In large ensembles of interacting oscillatory units, or-der emerges when a group starts showing the same fre-quency. This frequency adjustment, called synchroniza-tion [1], is a crucial and ubiquitous phenomenon in manyareas of science, such as neurosciences [2], semiconduc-tor laser arrays [3], cardiac pacemaker cells [4], powergrids [5], cell metabolism [6], Josephson junction arrays[7], and chemical oscillators [8].Arthur Winfree was the first to present a theoreticaldescription of synchronization in a model of biological os-cillators [9]. Influenced by Winfree’s work, Yoshiki Ku-ramoto introduced his famous minimal model of coupledoscillators [10] and its mathematical analysis [11–14]. To-gether these seminal works occupy a proeminent place insynchronization research. Since Kuramoto’s early works,a large number of Kuramoto-like models appeared inlater studies discussing effects such as those resultingfrom phase shifts in the coupling function [15], noise [16],periodic external fields [17], higher modes [18], finite size[19], and complex coupling networks [20].The Kuramoto model consists in an infinitely large en-semble of oscillators coupled globally. The oscillators arereduced to their phase, they are characterized by theirindividual natural frequency, and their dynamics is first-order in time. The remarkable discovery of Kuramotois that when coupled strongly-enough, the oscillators canovercome their nominal frequency quenched disorder andsynchronize.Kuramoto showed that his model exhibits a transitionbetween an incoherent state, where instantaneous fre-quencies are completely desynchronized, and a partiallysynchronized state, in which some oscillators share thesame instantaneous frequency, and are thus phase locked.The theoretical framework developed by Kuramoto to an-alyze his model, hereafter "Kuramoto theory", comprisesa set of assumptions and analytical results describing sta- tionary collective states [21]. Kuramoto assumed thatthese states would be characterized by phase distribu-tions [30] with stationary profiles. These profiles mightbe uniform (incoherent state) or be a steadily rotatingtraveling wave profile (synchronized state).The most common characterization of synchronizationin the Kuramoto model is, however, not in terms of phasedistributions but simply in terms of a scalar order param-eter (which is zero in the incoherent state and takes finitevalues when oscillators synchronize). Here we pursue yetanother, finer, description in terms of the distribution ofinstantaneous frequencies . Although synchronization isa direct manifestation of the way instantaneous frequen-cies are distributed, this problem was not, to our knowl-edge, addressed so far. Not even by Kuramoto, who in-stead solved the problem of the distribution of “coupling-modified frequencies” [11, 12], i.e. instantaneous frequen-cies averaged over an infinitely long time.The main goal of this paper is to extend Kuramoto the-ory by presenting a detailed derivation of the instanta-neous frequency distribution without any time-averagingprocedure. Our work is essentially based on Kuramotoresults and provides a mathematical description of theinstantaneous frequency distribution in stationary states,whether incoherent or with synchronized oscillators.This paper is structured as follows. In Section II webriefly discuss important aspects of Kuramoto theory, in-cluding the Kuramoto model, the order parameter, phasedistributions of synchronized and desynchronized oscilla-tors, and the distribution of time-averaged frequencies.In Section III, results from Kuramoto theory, togetherwith methods from probability and generalized functiontheory, are used to develop a geometric analysis thatsolves the problem of the instantaneous frequency dis-tribution. In Section IV, we illustrate the properties ofthis distribution in two cases: in the first one, we considerthe classic case of an unbounded normal distribution ofnatural frequencies; in the second, we adopt a symmet- a r X i v : . [ n li n . C D ] A ug ric Beta distribution of natural frequencies, defined ona bounded interval. In Section V, we obtain a powerseries expansion of the instantaneous frequency distribu-tion and an asymptotic formula, which allow us to studyrare events (large instantaneous frequency occurrences)in cases where the natural frequency distribution haspower-law tails, e.g., Cauchy-Lorentz distributions. Ourconclusions and open problems left for further investiga-tion are presented in Sec. VI. In the Appendix, the for-mula of the instantaneous frequency distribution is com-pared to instantaneous frequency histograms obtainedfrom numerical simulations of the Kuramoto model. II. KURAMOTO THEORY
In this section we present some results from the theo-retical framework developed by Kuramoto to analyze hiscoupled oscillator model. For details about how these re-sults can be obtained we refer the reader to Refs. [11, 21].The Kuramoto model consists of an ensemble of N all-to-all coupled oscillators with randomly distributednatural frequencies ω i ( i = 1 , , . . . , N ) whose phases θ i evolve according to: ˙ θ i = ω i + KN N (cid:88) j =1 sin( θ j − θ i ) , (1)where K is the coupling strength. In Kuramoto theory, N is assumed to be a infinitely large number and the naturalfrequencies ω i are randomly distributed according to agiven probability density function g ( ω ) .Collective states in the Kuramoto model are usuallyanalyzed by using measures which quantify the level ofsynchronization. One such quantity, proposed by Ku-ramoto, is R = lim N → + ∞ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N N (cid:88) j =1 exp ( iθ j ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (2)which we call here order parameter. If the oscillator stateis represented by exp ( iθ ) , R is the magnitude of a com-plex number representing the mean state of oscillators.A fundamental surmise in Kuramoto’s theory is that sta-tionary collective states, when they exist, are character-ized by time-independent values of the order parameterreached after sufficiently long time. This is equivalentto the assumption that the distribution of phases has awell-defined, steady profile. In an ordered state, the syn-chronized oscillators adopt a common frequency denoted Ω hereafter.In order to simplify notations, from here on we use thedefinitions a = KR and ˜ χ = χ − Ω a , (3)where χ is a generic quantity. Kuramoto’ analysis considers the use of a frame rotat-ing with angular velocity equal to Ω . The dynamics ofan oscillator of natural frequency ω can then be written ˙ ψ = ω − Ω − a sin ψ, (4)where ˙ ψ = ˙ θ − Ω and ψ = θ − Ω t are descriptions in therotating frame of the oscillator’s instantaneous frequency ˙ θ and phase θ .As Kuramoto did in his early works, here we distin-guish two groups of oscillators: synchronized and desyn-chronized oscillators. Synchronized oscillators are thosewhose natural frequencies satisfy | ˜ ω | ≤ , which meansthat (4) has a stable fixed point given by ψ ∗ ( ω ) = arcsin (˜ ω ) . (5)Desynchronized oscillators have natural frequencies suchthat | ˜ ω | > . In this case, Eq. (4) has no fixed point:phases show oscillatory dynamics, which look like relax-ation oscillations: they change rapidly at some stage oftheir period, then evolve much more slowly for the restof the cycle (see Fig. 1(b) for some examples).According to Eq. (4), if R = 0 , instantaneous and nat-ural frequencies are the same. For heterogeneous naturalfrequencies, R = 0 characterizes the incoherent state,where oscillators are out of synchrony. Assuming thatthe natural frequency distribution, g , is symmetric andunimodal, as Kuramoto does in his early works, then Ω coincides with the center of symmetry of g . When os-cillators synchronize, the order parameter has a positivevalue R = aK , obtained by finding the value of a whichsolves K + π ˆ − π dψg (Ω + a sin ψ ) cos ψ = 1 . (6)Eq. (6) can be used if the assumption of a symmetricand unimodal g holds. For more general profiles, theorder parameter has to be computed using a more general(and more difficult to solve) equation whose solution isdefined in terms of both a and Ω (see Ref. [15]). Inthis paper, our analytical results do not depend on how a is computed or on any specific assumption about g .However, we will here consider Eq.(6) in the numericalexamples, since it is simpler and more widely known [31].If g is symmetric and unimodal, the critical couplingstrength, i.e. the value of K marking the transition be-tween the desynchronized and the partially synchronizedstates, is given by K c = 2 πg (Ω) , (7)Again, a more general expression has to be used for moregeneral forms of g [15].Let p ( ψ, ω ) denote the joint probability density in-volving the oscillator’s phase in the rotating frame andthe oscillator’s natural frequency. Applying Bayes’ rule, p ( ψ, ω ) can be expressed p ( ψ, ω ) = p ( ψ | ω ) g ( ω ) , (8)where p ( ψ | ω ) is the conditional phase density for a givennatural frequency ω .A detailed discussion about how to obtain the condi-tional density p ( ψ | ω ) can be found in Ref. [21]. Here weonly show p ( ψ | ω ) in its final possible forms. For | ˜ ω | ≤ , p ( ψ | ω ) = δ [ ψ − ψ ∗ ( ω )] , (9)where ψ ∗ ( ω ) is the stable fixed point of Eq. (4), given byEq. (5). For | ˜ ω | > , Eq. (4) has no fixed point, and thedensity p ( ψ | ω ) can be be written as p ( ψ | ω ) = ω − Ω2 π ˙ ψ (cid:114) − ω , (10)where ˙ ψ is given by Eq. (4).From formulas equivalent to Eqs. (9) and (10), Ku-ramoto obtained the phase distribution, given by n ( ψ ) = n S ( ψ ) + n D ( ψ ) , (11)where n S ( ψ ) , the phase distribution of synchronized os-cillators, is n S ( ψ ) = (cid:40) g (Ω + a sin ψ ) a cos ψ, | ψ | ≤ π , | ψ | > π (12)and n D ( ψ ) denotes the phase distribution of desynchro-nized oscillators, given by n D ( ψ ) = 12 π ˆ | x | >a xg (Ω + x ) x − a sin ψ (cid:114) − (cid:16) ax (cid:17) dx. (13)For R > , the phase distribution profile shape is fixedand travels a distance of Ω t during a time interval t inthe non-rotating frame (where phases are described by θ ) . This is the scenario with synchronized oscillators: aphase distribution as a steadily traveling wave. In theincoherent state, R = 0 , whence, according to Eqs. (12)and (13), phases are uniformly distributed, viz., n ( ψ ) = π .Kuramoto also obtained the distribution of “coupling-modified frequencies” [11, 12], which are the instanta-neous frequencies averaged over an infinitely long time,as thoroughly discussed in Ref.[21]. The infinite-time av-erage ω of an oscillator’s instantaneous frequency ˙ θ , canbe defined as ω = lim T → + ∞ T T ˆ ˙ θ ( t ) dt. (14) Kuramoto showed that the coupling-modified frequenciesare distributed accordingly to G ( ν ) = δ ( ν − Ω) S ( K ) + G D ( ν ) , (15)where S ( K ) = Ω+ a ˆ Ω − a g ( ω ) dω, (16) G D ( ν ) = | ˜ ν |√ ν g (cid:18) Ω + ν − Ω | ˜ ν | (cid:112) ν (cid:19) (17)for ν (cid:54) = Ω , and G D (Ω) = 0 .Equation (15) states that: (i) G ( ν ) exhibits asingularity at the synchronization frequency Ω ; (ii) G ( ν ) goes linearly to zero for ν near Ω ; and (iii) lim (cid:15) → + ´ Ω+ (cid:15) Ω − (cid:15) G ( ν ) dν = S ( K ) , i.e. the probability ofan oscillator having a time-averaged frequency arbitrar-ily near Ω is given by S ( K ) , which represents the fractionof synchronized oscillators.Figure 1 summarizes these findings, together with somenumerical illustration of the object of central interesthere, the distribution of instantaneous frequencies. InFigs. 1(a) and (b) we show time series of the instan-taneous frequencies of eight oscillators selected from atotal ensemble of N = 5 × oscillators with theirnatural frequencies distributed according to a normal(Gaussian) distribution centered at Ω = 0 . For this case K c (cid:39) . . In Fig. 1(a), we set K = 0 . < K c in thedesynchronized regime ( a = 0 ), and all oscillators quicklykeep their natural frequency. In Fig. 1(b), K = 1 . > K c ,in the synchronized regime ( a > ): the 4 oscillators withtheir natural frequency | ˜ ω | < synchronize to Ω = 0 ,while the others stay desynchronized and their instan-taneous frequencies exhibit relaxation-oscillation-like dy-namics.In Figs. 1(c) and (d) we use the same K values as in(a) and (b) and show: i) g the Gaussian distribution ofnatural frequencies; ii) G D the continuous part of thedistribution of time-averaged instantaneous frequenciesgiven by Eq. (17); and iii) numerically-determined nor-malized histograms of instantaneous frequencies (see Ap-pendix for details). As expected, all these distributionscoincide in the subcritical case shown in Fig. 1(c). For K = 1 . > K c (Fig. 1(d)), the histogram of instanta-neous frequencies shows a peak located at Ω = 0 whichcorrespond to the synchronized oscillators. Note thatthe instantaneous frequencies (shown in the histogram)are distributed in a qualitatively different way from theirtime-averaged counterparts. Remarkable differences arethe accumulation of desynchronized oscillators’ instanta-neous frequencies near zero and the fact that the tailsof the instantaneous frequency distribution are “fatter"than those of the time-averaged and nominal frequencies.Coming back to the periodic dynamics of individualdesynchronized oscillators in the synchronized regime (a) K = 0 . (b) K = 1 . g ( ) G D histogram (c) K = 0 . g ( ) G D histogram (d) K = 1 . \ T K=0.8K=1.8 (e) (f)
FIG. 1. (a) Time-series of instantaneous frequencies for
K < K c . (b) Same as (a), but for K > K c . (c) g and G D comparedto normalized histograms of instantaneous frequencies. (d) Same as (c), but for K > K c . (e) Oscillation frequency, i.e. numberof phase-cycles per time unit ( T − ), defined by Eq.(18). (f) Time-averaged instantaneous frequency ( ω ), defined by Eq. (19).Numerical results were obtained using a normal distribution of natural frequencies and a number of oscillators N = 5 × .See Appendix for numerical details. (Fig. 1(b)), we see from Eq. (4) that the amplitude ofoscillations of ˙ θ is a , and the minimum and maximumvalues reached are ˙ θ min = ω − a and ˙ θ max = ω + a , so thatthe middle value is the natural frequency of the oscillator.As discussed in Ref. [21], a desynchronized oscillator’sphase makes a complete π turn during a time-interval T ( ω ) , which is given by T ( ω ) = 2 π (cid:113) ( ω − Ω) − a . (18)The time-averaged instantaneous frequency, defined by(14), can be given in terms of T as ω = (cid:40) Ω + πT , ω > Ω + a Ω − πT , ω < Ω − a. (19)If a = 0 , we can write (18) and (19) as T = | ω − Ω | π , and ω = ω . If a > , (18) is the same as T = a π √ (cid:101) ω − , and(19) can be written as ω = Ω ± a √ (cid:101) ω − for ± (cid:101) ω > .These formulas allow us to infer the following proper-ties of the instantaneous frequency of non-synchronizedoscillators: i) as (cid:101) ω → ± ± , we have T → , which meansslow oscillations, or oscillations with large-time periods; ii) if (cid:101) ω → + , then ω → Ω + and ˙ θ min = ω − a → Ω + ,i.e. both the time-averaged and minimum instantaneousfrequencies have values close to and greater than Ω , thefrequency of synchronized oscillators (Similar propertiesare of course valid if (cid:101) ω → − ); iii) as | (cid:101) ω | → ∞ , thenwe have T → ∞ and ω → ω , viz. fast oscillations withtime-averages becoming close to their oscillation centers.Fig. 1(b) illustrates all these properties. The long peri-ods of time spent near Ω by the instantaneous frequencyof oscillators with | ˜ ω | > ∼ explains why the histogram inFig. 1(d) exhibits frequent occurrences near zero.Figs. 1(e) and (f) show plots of T and ω as functionsof ω created using Eqs. (18) and (19). The synchronizedfrequency Ω has been set to zero and we adopted thesame values of K as used in Figs.1(a)-(d). For K = 0 . ( a = 0 ), the graphs are shown in blue. The curves inorange correspond to the case where K = 1 . ( a (cid:39) ).In summary, many of the properties shown in Fig. 1are straight consequences of Kuramoto theory. Thoseregarding instantaneous frequencies can be qualitativelyexplained by it. This is the case of the accumulation ofinstantaneous frequencies near the synchronization fre-quency (Fig. 1(d)). We now proceed to the core of thiswork, which is the calculation of the full analytical ex-pression of the distribution of instantaneous frequencies. III. DISTRIBUTION OF INSTANTANEOUSFREQUENCIES
Our goal in this section is to obtain, based on the re-sults discussed in Section II, the distribution of instanta-neous frequencies in the Kuramoto model. This distribu-tion is a probability density function G ( ν ), which meansthat G ( ν ) dν is the probability of an oscillator showingits fixed frame instantaneous frequency , ˙ θ, in the interval [ ν, ν + dν ) .By using the random variable transformation theorem[24], we have G ( ν ) = ˆ + ∞−∞ ˆ + π − π δ [ ν − ˙ θ ( ω, ψ )] p ( ψ, ω ) dψdω, (20)where δ denotes the delta function, ˙ θ ( ω, ψ ) is a randomvariable transformation, given by ˙ θ ( ω, ψ ) = ω − a sin ψ, (21)and p ( ψ, ω ) is the joint probability density involving thephase in the rotating frame and the natural frequency.Equation (21) comes from Eq.(4), which decribes instan-taneous frequencies in the rotating frame.The delta function in Eq. (20), δ [ ν − ˙ θ ( ω, ψ )] , is con-centrated in a curve embedded in the two-dimensionalspace defined by the integration variables ψ and ω . Inorder to calculate the double integral in Eq. (20), we usea method proposed by Seeley [25, 26], which generalizesthe usual concept of one-dimensional delta functions todelta functions concentrated in manifolds with an arbi-trary number of dimensions.Let δ ( P ) denote a delta function concentrated in a n − -dimensional manifold M embedded in a n -dimensionalspace V . The manifold M is defined by P ( x ) = 0 , where P ( x ) is a function at x = ( x , . . . , x n ) , i.e. a point in V with coordinates x , . . . , x n . The delta function δ ( P ) can be defined by δ ( P ) = lim c → + Θ( P + c ) − Θ( P ) c , (22)where Θ( . ) is a Heaviside step function such that Θ( P ) = 1 for P ≥ , and Θ( P ) = 0 for P < .Consider a function ϕ ( x ) defined in V . From Eq. (22),we have ˆ V δ ( P ) ϕ ( x ) d n x = lim c → + ˆ − c ≤ P < ϕ ( x ) c d n x , (23)where d n x = dx . . . . .dx n . Let γ be the distance between the manifolds M P =0 and M P = − c , defined by the equations P = 0 and P = − c , re-spectively. Then, for P ( x ) = 0 , we have P [ x + γ ˆ n ( x )] = − c , i.e., if x ∈ M P =0 , then x + γ ˆ n ( x ) ∈ M P = − c ,where ˆ n ( x ) is the unit vector normal to M P =0 at x . As c → + , a first-order expansion of P [ x + γ ˆ n ( x )] resultsin P ( x ) − γ ˆ n ∇ P = − c , from which we obtain γ = c |∇ P | , (24)since P ( x ) = 0 and ˆ n ( x ) = ∇ P |∇ P | .Changing the infinitesimal volume element d n x by γdS , where γ is given by Eq. (24) and dS is an infinites-imal surface element of M P =0 , Eq. (23) can be rewrittenas ˆ V δ ( P ) ϕ ( x ) d n x = ˆ M P =0 ϕ ( x ) |∇ P | dS, (25)which means that the volume integral in the right-side ofEq. (23) can be changed by a surface integral on M P =0 .Suppose that M P =0 is a curve, and V, a two-dimensional space. In this particular case, Eq. (25) canbe written in the form ˆ V δ [ P ( x )] ϕ ( x ) d x = ˆ M P =0 ϕ ( x ) |∇ P | dl, (26)where x = ( x , x ) , d x = dx dx , and dl is an infinites-imal line element of M P =0 . Let P ( x , x ) be defined by P ( x , x ) = f ( x ) − x , where f ( x ) is a continuous func-tion, and the range of x is the the interval [ a, b ) . Then,the curve M P ( x ,x )=0 is the graph of x = f ( x ) with a ≤ x < b . Suppose that C M is a curve correspond-ing to a part of M P ( x ,x )=0 . This curve can be definedas a subset of M P ( x ,x )=0 by C M = { ( x , x ) ∈ V | x = f ( x ) and ψ a ≤ x ≤ ψ b < b } . An integral along C M , analogous to the right-hand side of Eq. (26), can bewritten as ˆ C ϕ ( x , x ) |∇ P ( x , x ) | dl = ψ b ˆ ψ a ϕ [ x , f ( x )] dx . (27)We are now able to compute the right-hand side ofEq. (20). Let us first put Eq. (20) in the form G ( ν ) = ˆ V δ [ P ν ( ψ, ω )] p ( ψ, ω ) dψdω, (28)where the integration manifold, V , is an infinite-lengthcylinder V = [ − π , + π ) × ( −∞ , + ∞ ) , P ν ( ψ, ω ) is givenby P ν ( ψ, ω ) = F ν ( ψ ) − ω, (29)with F ν ( ψ ) = a sin ψ + ν, (30)and − π ≤ ψ < + π .Then, M P ν =0 is a closed curve in V defined by M P ν =0 = { ( ψ, ω ) ∈ V | ω = F ν ( ψ ) } . (31) /2 0 + /2 + +3 /2 a + a M P =0 FIG. 2. Graph of F ν , Eq. (30), representing the curve M P ν =0 ,Eq. (31). The minimum is (cid:0) − π , ν − a (cid:1) , and the point ofmaximum is (cid:0) + π , ν + a (cid:1) . This curve is represented by the graph of F ν , shown inFig. 2. The position of the curve M P ν =0 in the integrationmanifold V is determined by the value of ν , which is theargument of G . And the height of M P ν =0 , as shown inFig. 2, is a . The curve shifts by varying ν and stretchesas the product a increases.Using the relation (26), we obtain from Eq. (28) theformula G ( ν ) = ˆ M Pν =0 p ( ψ, ω ) |∇ P ν ( ψ, ω ) | dl. (32)We can calculate the line integral in Eq. (32) by use of ageometric analysis based on dividing V in two disjointsregions, V S = (cid:2) − π , + π (cid:1) × [ − a + Ω , Ω + a ] and V D = V − V D . A sketch of both regions and M P ν =0 in differentlocations is shown in Figs.3(a)-(e). Depending on thelocation and height of M P ν =0 , this curve is completelyinside V D (Figs. 3 (a) and (e)) and can also be partly orentirely in V S (Figs. 3(b),(c) and (d)). We denote theparts in V D by M DP ν =0 , which are the blue curves, andthose in V S by M SP ν =0 , represented by the yellow curves.Using Eq. (8), we have that Eq. (32) is the same as G ( ν ) = G S ( ν ) + G D ( ν ) , (33)where G S ( ν ) = ˆ M SPν =0 p ( ψ | ω ) g ( ω ) |∇ P ν ( ψ, ω ) | dl, (34)and G D ( ν ) = ˆ M DPν =0 p ( ψ | ω ) g ( ω ) |∇ P ν ( ψ, ω ) | dl. (35)If M P ν =0 has no part inside V s , then the curve M SP ν =0 does not exist and G S ( ν ) = 0 . Similarly, if M DP ν =0 is anempty set, then we also have G D ( ν ) = 0 . /2 0 + /2 + +3 /2 a + aa + aV S V D V D a) < 2 M DP = 0 /2 0 + /2 + +3 /2 a + aa + aV S V D V D b) 2 < < 0 M SP = 0 M DP = 0 /2 0 + /2 + +3 /2 a + aV S V D V D c) = 0 M SP = 0 /2 0 + /2 + +3 /2 a + aa + aV S V D V D d) 0 < < +2 M SP = 0 M DP = 0 /2 0 + /2 + +3 /2 a + aa + aV S V D V D e) +2 < M SP = 0 M DP = 0 FIG. 3. The curve M P ν =0 in different locations depending onthe value of ν . The part of M P ν =0 immersed in V D is M DP ν =0 ,and the intersection between M P ν =0 and V S is M SP ν =0 . − < ˜ ν ≤ < ˜ ν < +2 ψ S ( ν ) - arcsin (1 + ˜ ν ) − π ψ S ( ν ) + π arcsin (1 − ˜ ν ) TABLE I. Integration limits of Eq. (37), given by ψ S ( ν ) and ψ S ( ν ) . For | ˜ ν | ≥ , ψ S ( ν ) and ψ S ( ν ) are not defined, and G S ( ν ) = 0 . We consider first the case in which M SP ν =0 exists. Apoint ( ψ, ω ) in M SP ν =0 satifies the conditions: i) ω = F ν ( ψ ) ; ii) − a + Ω ≤ ω ≤ Ω + a , as can also be seenin Figs. 3(b),(c) and (d). Condition ii) means that thedensity p ( ψ | ω ) in Eq. (34) is defined by Eq. (9). Then G S ( ν ) = ˆ M SPν =0 δ [ ψ − ψ ∗ ( ω )] g ( ω ) |∇ P ν ( ψ, ω ) | dl, (36)where ψ ∗ ( ω ) is given by Eq. (5). Integration along M SP ν =0 can be done in three steps: the first one is in-tegration along M S L P ν =0 , i.e. the subset of M SP ν =0 whoseprojection in the ψ -axis is contained in the closed interval (cid:2) − π , + π (cid:3) ; the second step is integration along M S R P ν =0 ,which is the subset of M SP ν =0 whose projection in the ψ -axis is contained in the open interval (cid:0) − π , + π (cid:1) ; the laststep consists in summing the results of both integrations.For a point ( ψ, ω ) in M S R P ν =0 , we have + π < ψ < + π and − π ≤ ψ ∗ ( ω ) ≤ + π . Then, ψ − ψ ∗ ( ω ) > and δ [ ψ − ψ ∗ ( ω )] = 0 , which means that the integral along M S R P ν =0 is zero, and the integral along M SP ν =0 in Eq. (36)reduces to the integral along M S L P ν =0 .If M S L P ν =0 is a non-empty set, the projection of M S L P ν =0 in the ψ -axis can be represented by the interval (cid:2) ψ S ( ν ) , ψ S ( ν ) (cid:3) , and we can use (27) and (36) to obtain G S ( ν ) = ˆ ψ S ( ν ) ψ S ( ν ) δ { ψ − ψ ∗ [ F ν ( ψ )] } g [ F ν ( ψ )] dψ, (37)where F ν is defined by Eq. (30), and the integration lim-its, ψ S ( ν ) and ψ S ( ν ) , are given in Table I.For | ˜ ν | ≥ , ψ S ( ν ) and ψ S ( ν ) are not defined, and G S ( ν ) = 0 . This case is illustrated by Figs. 3(a) and(e), which show M P ν =0 entirely outside V S . The case − < ˜ ν ≤ is related to configurations with M P ν =0 partly or completely inside V S , such as those depicted inFigs. 3(b)-(c). When < ˜ ν < +2 , there is still partialembedding of M P ν =0 in V s (See Fig. 3(d)).From the definitions of ψ ∗ and F ν (see Eqs. (5) and(30)), the delta function in the integrand of Eq. (37) isgiven by δ { ψ − ψ ∗ [ F ν ( ψ )] } = δ [ ψ − arcsin (sin ψ + ˜ ν )] , (38)which is singular when ψ = arcsin (sin ψ + ˜ ν ) . (39) By applying the sine function on both sides of Eq. (39),we see that the singularity occurs if and only if ν = Ω independently of the value taken by ψ . So Eq. (38) onlymakes sense in an integral where ν is the integration vari-able and will remain as a delta function despite of inte-gration in Eq. (37).From Eq. (38), δ { ψ − ψ ∗ [ F ν ( ψ )] } in its Gaussian formreads δ { ψ − ψ ∗ [ F ν ( ψ )] } = L [ ψ − arcsin (sin ψ + ˜ ν )] , (40)where L ( x ) = lim (cid:15) → + (cid:15) √ π exp (cid:20) − (cid:16) x(cid:15) (cid:17) (cid:21) . (41)This limit can be analyzed in an arbitrarily small openneighborhood N of radius (cid:15) centered in the singularitypoint, defined by ν = Ω . As expected for delta functionsin non-singular points, if v is not in N , i.e. | ν − Ω | ≥ (cid:15) ,the limit in (40) is zero, since the Gaussian is O ( (cid:15) m − ) as (cid:15) → + for any positive integer m . If v is inside N , i.e. | ν − Ω | < (cid:15) , then arcsin (sin ψ + ˜ ν ) = ψ + ˜ ν cos ψ + O (cid:0) ˜ ν (cid:1) , (42)for (cid:15) −→ + . By redefining (cid:15) as (cid:15) cos ψ and substituting(42) in (40), we obtain δ { ψ − ψ ∗ [ F ν ( ψ )] } = cos ψL (˜ ν ) , (43)where L (˜ ν ) is the Gaussian representation of δ (˜ ν ) . Sub-stituting (43) in (37), we have G S ( ν ) = δ (˜ ν ) ˆ ψ Sb ( ν ) ψ Sa ( ν ) g ( a sin ψ + ν ) cos ψ dψ, (44)which is the same as G S ( ν ) = aδ ( ν − Ω) ˆ + π − π g ( a sin ψ + Ω) cos ψ dψ. (45)Changing the integration variable ψ to ω = a sin ψ + Ω ,Eq. (45) results in G S ( ν ) = δ ( ν − Ω) S ( K ) , (46)where S ( K ) , as mentioned in Sec.II, is the fraction ofsynchronized oscillators, defined by Eq. (16). It is worthmentioning that G S ( ν ) is identical to the singular termin the distribution of time-averaged frequencies, as shownby Eq. (15).A similar geometric analysis can be used to calculate G D ( ν ) from Eq. (35), where integration is now performedalong the curve M DP ν =0 . As mentioned before, this curvecorresponds to the part of M P ν =0 in V D . A point ( ψ, ω ) in M DP ν =0 satifies the conditions: i) ω = F ν ( ψ ) ; ii) ω < Ω − a or ω > Ω + a (see orange curves in Figs. 3(a),(b),(d),and (e)). From condition i), (30), and (4) , we have that ˜ ν ≤ − − < ˜ ν < < ˜ ν < +2 +2 ≤ ˜ νψ D ( ν ) − π − π arcsin (1 − ˜ ν ) − π ψ D ( ν ) + π - arcsin (1 + ˜ ν ) + π + π TABLE II. Integration limits of Eq. (49), given by ψ D ( ν ) and ψ D ( ν ) . The limits are not defined for ν = Ω . (49) does notapply, and G D ( ν ) = 0 . ˙ ψ = ν − Ω . From condition ii), p ( ψ | ω ) is defined by (10),and Eq. (35) can then be rewritten as G D ( ν ) = 12 π ˜ ν ˆ M DPν =0 ˜ ωg ( ω ) |∇ P ν ( ψ, ω ) | (cid:114) − ω dl. (47)For ν = Ω , M DP ν =0 is an empty set, as shown in Fig. 3(c),and G D ( ν ) = 0 .Let U ( x ) be defined as U ( x ) = (cid:112) x − . (48)From Eq.(27) and condition i), (47) reads G D ( ν ) = 1 π | ˜ ν | ˆ ψ D ( ν ) ψ D ( ν ) g [ F ν ( ψ )] U (cid:104) ˜ F ν ( ψ ) (cid:105) dψ, (49)where ψ D ( ν ) and ψ D ( ν ) are given in Table II.As G S and G D are given by Eqs. (46) and (51), wecan now return to Eq. (33) and write G in its final form: G ( ν ) = δ ( ν − Ω) S ( K ) + G D ( ν ) , (50)where G D (Ω) = 0 , and G D ( ν ) = 1 π | ˜ ν | ˆ ψ + ν ψ − ν h ν ( ψ ) dψ, (51)if ν (cid:54) = Ω . In Eq. (51), the quantities h ν ( ψ ) , χ ν , ψ − ν , and ψ + ν are defined by h ν ( ψ ) = g ( a sin ψ + ν ) U (sin ψ + ˜ ν ) (52)and ψ ± ν = arcsin {− (˜ ν ±
2) Θ [ − ˜ ν (˜ ν ± ± } , (53)where Θ [ . ] denotes the Heaviside step function. Notethat ψ − ν and ψ + ν correspond to the previously defined ψ D ( ν ) and ψ D ( ν ) .The singular term in (50) means that lim (cid:15) → + Ω+ (cid:15) ˆ Ω − (cid:15) G ( ν ) dν = S ( K ) , (54) i.e. the probability of an oscillator with instantaneousfrequency in an infinitesimally small neighborhood of Ω is given by the fraction of synchronized oscillators.Our procedure to obtain Eq. (50) does not depend inany symmetry assumption related to g . But let us nowassume the situation where g (Ω + x ) = g (Ω − x ) for anypositive number x . This implies that G has the sameproperty, viz. if g (Ω + x ) = g (Ω − x ) , then G (Ω + x ) = G (Ω − x ) . To show this, it is sufficient showing that G D is also symmetric. For x > , we have G D (Ω ± x ) = aπx ˆ ψ +Ω ± x ψ − Ω ± x h Ω ± x ( ψ ) dψ, (55)and ψ − Ω ± x = − ψ +Ω ∓ x . In the formula for G D (Ω + x ) (seeEq. (55)), we can introduce the following changes: first,we change ψ − Ω+ x by − ψ +Ω − x and ψ +Ω+ x by − ψ − Ω − x ; second,we redefine ψ as − ψ . Then, G D (Ω + x ) = aπx ˆ ψ +Ω − x ψ − Ω − x h Ω+ x ( − ψ ) dψ (cid:48) , (56)Since g [Ω − ( a sin ψ − x )] = g ( a sin ψ + Ω − x ) (fromthe symmetry assumption of g ) and U (cid:0) − sin ψ + xa (cid:1) = U (cid:0) sin ψ − xa (cid:1) , we have h Ω+ x ( − ψ ) = h Ω − x ( ψ ) . Then,from Eq. (56), G D (Ω + x ) = G D (Ω − x ) , which provesour initial statement.According to Eq. (50), G consists of a delta peakand a distribution of instantaneous frequencies for non-synchronized oscillators ( G D ). G D is zero at Ω , wherethe delta peak is located. G has a similar form: the samedelta peak at Ω , and a distribution of time-averaged in-stantaneous frequencies for non-synchronized oscillators( G D ). G D is also zero at Ω . Since a synchronized oscil-lator’s instantaneous frequency goes to Ω as time goes toinfinity, the same happens to its average in time. Thisexplains why the fraction of synchronized oscillators, de-fined by S ( K ) , appears as a factor in the delta peaks ofboth G and G . IV. APPLICATION TO GAUSSIAN AND BETADISTRIBUTIONS
In this section we illustrate our analytical results aboutthe instantaneous frequency distributions on two distri-butions of natural frequencies, the normal (Gaussian) dis-tribution and the Beta distribution. Both are unimodal,but one has unbounded support, whereas the other liveson a finite interval (Beta distribution).
A. General features
The Gaussian natural frequency distribution consid-ered is g ( ω ) = 1 √ π exp (cid:18) − ω (cid:19) , (57) g ( ) G D , K = 1.596 G D , K = 1.61 G D , K = 1.64 G D , K = 1.67 G D , K = 1.8 (a) g ( ) G D , K = 0.42442 G D , K = 0.43 G D , K = 0.435 G D , K = 0.45 G D , K = 0.5 (b) FIG. 4. (a) Instantaneous frequency distribution ( G ) for natural frequencies normally distributed. The thin vertical lines showdifferent locations of Ω − and Ω + , defined by Ω ± = Ω ± a . Fractions of synchronized oscillators ( S ( K ) ): S (1 . . ; S (1 .
61) = 0 . S (1 .
64) = 0 . ; S (1 .
67) = 0 . ; S (1 .
8) = 0 . . (b) Graphs of G assuming a Beta(2,2) distribution of naturalfrequencies. S ( K ) : S (0 . . ; S (0 .
43) = 0 . , S (0 . . ; S (0 .
45) = 0 . ; S (0 .
5) = 0 . . for which Ω = 0 , and, according to Eq. (7), K c (cid:39) . .The Beta distribution considered reads g ( ω ) = (cid:40) ω (1 − ω ) ω ∈ [0 , ω / ∈ [0 , . (58)which is usually called Beta (2 , . All members of thefamily of Beta distributions have the support interval [0 , . So in this example oscillators have no naturalfrequencies outside the interval [0 , . Given the sym-metric shape, the synchronization frequency is Ω = 0 . .The critical coupling strength, given by Eq. (7), is K c (cid:39) . .In Fig. 4, we show these two distributions (thick bluecurves), but also the distribution of instantaneous fre-quencies for different values of the coupling strengthabove the critical coupling K c . (For subcritical couplingvalues, the instantaneous frequencies are the natural fre-quencies.) For K > K c , G ( ν ) = G S ( ν ) + G D ( ν ) , where G S ( ν ) and G D ( ν ) are defined by Eqs. (46) and (51),respectively. Since G S ( ν ) is a Dirac delta term, we rep-resent it by a black vertical line located in Ω .The thin colored curves show G D for different valuesof K . This continuous part of G obeys the symmetry of g , the distribution of natural frequencies. The tails of G D are fatter than those of g . In particular, G D extendsbeyond the interval of support of g in the Beta case. Inthe central region, G D ( ν ) < g ( ν ) . For K > ∼ K c , G D isvery close to g , but for a sharp drop near Ω , the synchro-nized frequency (orange curves). This drop, however,does not extend to zero: G D ( ν ) tends to a finite valuewhen ν → . Increasing K , G D develops a more com-plicated structure: the central region decreases, the tails grow, and some special values of ν appear, marked bythin gray vertical lines on the figure. They indicate thelocations of Ω − and Ω + , defined by Ω ± = Ω ± a . Thequantities Ω − , Ω , and Ω + are endpoints of intervals whichdefine the integration limits of G D (see Table II). Since G D is a piecewise function, the graph of G D consistsof four sub-graphs associated to the intervals ν ≤ Ω − , Ω − < ν < Ω , Ω < ν < Ω + , and Ω + ≤ ν .In Figures 5 (normal distribution) and 6 (Beta distri-bution), we compare the distributions of instantaneous,time-averaged and natural frequencies, again for differentvalues of K .Except for subcritical values of K or near the transition(panels (a) and (b) of each figure), the graphs of G D and G D are quite different from each other as ν → Ω . In par-ticular while G D ( ν ) show a big dip to zero for K > K c , G D ( ν ) approaches non-zero values near the synchroniza-tion frequency.For the Beta distribution (Fig. 6) the tails of G D and G D reach zero for large enough values of | ν | . As K in-creases, the support interval shrinks for G D , while it ex-pands for G D . B. A focus on tails
The tails of G describe rare events, viz. large instan-taneous frequencies with small occurrence probabilities.Here we analyze these events for the Gaussian and Betacases examined above.Figure 7 shows the decimal logarithms of G D (panels(a) and (c)) and G D (panels (b) and (d)) for the normaland Beta distributions studied above, using the same set0 g ( ) G D ( ) G D ( ) (a) K < K c (cid:39) . delta term g ( ) G D ( ) G D ( ) (b) K = 1 . delta term g ( ) G D ( ) G D ( ) (c) K = 1 . delta term g ( ) G D ( ) G D ( ) (d) K = 1 . delta term g ( ) G D ( ) G D ( ) (e) K = 1 . delta term g ( ) G D ( ) G D ( ) (f) K = 1 . FIG. 5. Comparison between the instantaneous frequency distribution, G ( ν ) , and the time-averaged frequency distribution, G ( ν ) . A normal distribution of natural frequencies, g ( ν ) , is assumed. For K < K c , g ( ν ) = G ( ν ) = G ( ν ) , and there is no deltaterm, which means that S ( K ) = 0 . of K values as in previous figures.For clarity, in the Gaussian case, the distributions areplotted as functions of ν , so that Gaussian tails ap-pear as straight lines. For large values of | ν | , the tails of G D stay above the tails of g , and the difference between log [ G D ( ν )] and log [ g ( ν )] increases with K . Neverthe-less, all G D distributions keep the same asymptotic tailas g , rescaled by a K -dependent factor (Fig. 7(a)). Thetails of G D are also Gaussian, and asymptotically iden-tical but below those of g , rescaled by a K -dependentfactor that decreases with K (Fig. 7(c)).In the case of the Beta distribution, both G D and G D have a bounded support, and they behave in aqualitatively-similar manner to g near the limit valuesof their support intervals. The tails of G D extends be-yond the support interval of g , and all the more so as K increases (Fig. 7(b)), while the tails of G D show theopposite tendency (Fig. 7(d)).Further information about the above results is pre-sented in Figure 8: In panels (a) and (b), we show thedecimal logarithms of the ratios G D /g and G D /g , whichshows clearly that G D “goes away” from g as K increases:in the central region G D becomes smaller and smallerthan g ; in the tails, the difference grows. In comparison,the behavior of G D is much more “gentle”. In the Betacase, we show how the bounds of the support interval varies with K (Fig 8(c)). For the average frequencies( G D ), the support shrinks almost linearly with K , whileit grows more slowly than linearly for G D . Finally, inFig. 8(d) we show A , the area of G D beyond the supportinterval of g . This quantity measures the overall like-lihood to observe instantaneous frequencies beyond thepossible nominal frequencies. Interestingly, A first growswith K , then decreases, in spite of the monotonous in-crease of the support of G D . (In Fig. 7(b), one can un-derstand that this comes from the increasingly trimodalnature of G D .) We also plot A , the area of g outside thesupport interval of G D , which indicates the overall weightof natural frequencies unobservable as time-averaged fre-quencies. It increases monotonously with K . V. RARE EVENTS AND POWER-LAW TAILS
In this section we consider natural frequency distri-butions with power-law tails and develop a power seriesexpansion of Eq. (51) in order to deepen our understand-ing of rare events. By rare events we mean occurrences oflarge instantaneous frequency values such that ν (cid:28) Ω − or ν (cid:29) Ω + . We assume that natural frequency distri-butions are smooth and have unimodal and symmetricprofiles centered at Ω = 0 .1 g ( ) G D ( ) G D ( ) (a) K < K c (cid:39) . delta term g ( ) G D ( ) G D ( ) (b) K = 0 . delta term g ( ) G D ( ) G D ( ) (c) K = 0 . delta term g ( ) G D ( ) G D ( ) (d) K = 0 . delta term g ( ) G D ( ) G D ( ) (e) K = 0 . delta term g ( ) G D ( ) G D ( ) (f) K = 0 . FIG. 6. Comparison between G and G . g is a Beta (2 , distribution. (a) For K < K c , g ( ν ) = G ( ν ) = G ( ν ) , and there isno delta term, which means that S ( K ) = 0 . (b-f) As K increases, the bounding interval of G decreases its width, while thebounding interval of G ( ν ) becomes larger. According to Table II, the instantaneous frequency dis-tribution for
Ω = 0 and | ν | > a can be written G ( ν ) = 1 π | ν | ˆ + π − π h ( ν + a sin ψ ) dψ. (59)with h ( x ) = u + ( x ) g ( x ) u − ( x ) where u ± ( x ) = √ x ± a . (60)Expanding h ( ν + a sin ψ ) as a Taylor series, we have h ( ν + a sin ψ ) = h ( ν ) + ∞ (cid:88) m =1 a m m ! h ( m ) ( ν ) sin m ψ, (61)where h ( m ) denotes the m th-order derivative of h . Sub-stituting Eq. (61) in Eq. (59), we obtain G ( ν ) = h ( ν ) | ν | + 1 π | ν | ∞ (cid:88) m =1 a m m ! h ( m ) ( ν ) ˆ + π − π sin m ψ dψ . (62)For any integer n > , ´ + π − π sin n ψ dψ = 2 ´ + π sin n ψ dψ and ´ + π − π sin n − ψ dψ = 0 . So only even order termsare present in Eq. (62). According to Eq. (3.621-3) in Ref.[27], ´ + π sin n ψ dψ = π (2 n − n )!! , whence G ( ν ) = h ( ν ) | ν | + 1 | ν | ∞ (cid:88) n =1 a n (2 n )! (2 n − n )!! h (2 n ) ( ν ) . (63)The Leibniz derivative rule allows us to write h (2 n ) as h (2 n ) ( ν ) = (cid:88) k + k + k =2 n (2 n )! k ! k ! k ! × u ( k )+ ( ν ) g ( k ) ( ν ) u ( k ) − ( ν ) , (64)where summation is taken over all partitions ( k , k , k ) of n into non-negative integers, g ( k ) and u ( k ) ± denote the k th-order derivatives of g . The latter are givenby u ( k ) ± ( ν ) = p k (cid:18) (cid:19) ( ν ± a ) − k u ± ( ν ) , (65)where p k ( q ) = k − (cid:89) l =0 ( q − l ) . (66)2 log [ g ( )]log[ G D ( )], K = 1.596log[ G D ( )], K = 1.61log[ G D ( )], K = 1.64log[ G D ( )], K = 1.67log[ G D ( )], K = 1.8 (a) normal. log [ g ( )] log [ G D ( )], K = 0.42442 log [ G D ( )], K = 0.43 log [ G D ( )], K = 0.435 log [ G D ( )], K = 0.45 log [ G D ( )], K = 0.5 (b) beta. log [ g ( )]log[ G D ( )], K = 1.596log[ G D ( )], K = 1.61log[ G D ( )], K = 1.64log[ G D ( )], K = 1.67log[ G D ( )], K = 1.8 (c) normal. log [ g ( )] log [ G D ( )], K = 0.42442 log [ G D ( )], K = 0.43 log [ G D ( )], K = 0.435 log [ G D ( )], K = 0.45 log [ G D ( )], K = 0.5 (d) beta. FIG. 7. Decimal logarithms of G D ( log [ G D ( ν )] ) and G D ( ν ) ( log (cid:2) G D ( ν ) (cid:3) ) assuming that g is a normal and a Beta (2 , PDF.
From Eqs. (64) and (65), it follows that h (2 n ) ( ν ) = u + ( ν ) u − ( ν ) (cid:88) k + k + k =2 n p k ( ) (2 n )! p k ( ) k ! k ! k ! × g ( k ) ( ν )( ν + a ) k ( ν − a ) k . (67)We can now use Eqs. (67), and (60) in Eq. (63) to obtainthe ratio G ( ν ) /g ( ν ) , which is given by Gg ( ν ) = (cid:114) − (cid:16) aν (cid:17) [1 + Λ( ν )] , (68)where Λ( ν ) = 1 g ( ν ) ∞ (cid:88) n =1 c n ( ν ) (cid:16) aν (cid:17) n (69) and c n ( ν ) = (2 n − n )!! (cid:88) k + k + k =2 n p k ( ) p k ( ) k ! k ! k ! × ν k g ( k ) ( ν ) (cid:0) aν (cid:1) k (cid:0) − aν (cid:1) k . (70)(If g is centered at a non-zero synchronization frequency Ω , more general formulas than Eqs. (68), (69), and (70)can be obtained by changing in them the terms aν by a ( ν − Ω) .)As an application of the result given by (68), let usnow consider a class of natural frequency distributions ofthe form g ( ν ) ∼ Cν − µ ( | ν | −→ ∞ ) , (71)3
12 6 0 6 120.60.00.61.21.8 l o g [ G D () / g ()] K=1.61K=1.64K=1.67K=1.8 (a) normal
12 6 0 6 120.60.00.61.21.8 l o g [ G D () / g ()] K = 1.61K = 1.64K = 1.67K = 1.8 (b) normal K ++ (c) beta K AA (d) beta FIG. 8. (a)Decimal logarithms of the ration G ( ν ) g ( ν ) : as ν increases, G ( ν ) g ( ν ) diverges. (b)Decimal logarithms of G ( ν ) g ( ν ) remains constantas ν increases. (c)Positive endpoints of the support intervals of G and G , denoted by ν + and ν + . (d) Area of G outside thesupport interval of g ( A ), and area of g outside the support interval of G ( A ). where µ is a positive integer, and C a real constant.Its derivatives read g ( k ) ( ν ) ∼ p k ( − µ ) ν − k g ( ν ) ( | ν | −→ ∞ ) , (72)where p k is given by Eq. (66), and k = 0 , ..., n . Since (cid:0) ± aν (cid:1) k ∼ as (cid:12)(cid:12) aν (cid:12)(cid:12) −→ , it follows from (70) and (72)that c n ( ν ) ∼ c n g ( ν ) (cid:16)(cid:12)(cid:12)(cid:12) aν (cid:12)(cid:12)(cid:12) −→ (cid:17) , (73)where the constant coefficient c n is defined by c n = (2 n − n )!! (cid:88) k + k + k =2 n p k (cid:0) (cid:1) p k ( − µ ) p k (cid:0) (cid:1) k ! k ! k ! . (74) Therefore, Λ( ν ) ∼| aν | → ∞ (cid:88) n =1 c n (cid:16) aν (cid:17) n , (75)and Gg ( ν ) ∼| aν | → (cid:113) − (cid:0) aν (cid:1) (cid:110) c (cid:0) aν (cid:1) + O (cid:104)(cid:0) aν (cid:1) (cid:105)(cid:111) . (76)This means: assuming that g ( ν ) has power-law tails ofform (71), G ( ν ) approaches asymptotically g ( ν ) for largeinstantaneous frequencies (compared to a ) or small orderparameter values (compared to (cid:12)(cid:12) Kν (cid:12)(cid:12) ).To illustrate this point, we consider the family of nat-4 g ( ) g ( ) g ( ) g ( ) FIG. 9. Graphs of g µ , defined by (77), for µ = 1 , , , . ural frequency distributions g µ ( ν ) = µπ (1 + ν µ ) sin (cid:18) π µ (cid:19) , (77)where µ is a positive integer. Formula (77) generalizesthe standard Cauchy-Lorentz distribution, which corre-sponds to the particular case g . Graphs of g µ ( ν ) areshown in Fig. 9 for µ = 1 , , , . By increasing µ , thetails of g µ ( ν ) gets thinner, and high natural frequencieshave lower occurrence probabilities. In Figs. 10(a)-(d),we show graphs of G D and G for different values of K considering the cases g and g .Figures 11(a) and (c) show decimal logarithms of theratio G D /g , where G D ( ν ) is computed using Eq. (17)with g ( ν ) = g µ ( ν ) . Like in the Gaussian and Beta ex-amples of Sec. IV, the graphs of G D are more and morebelow the graph of g as increasing K increases. Yet, for | ν | → ∞ , G D ( ν ) and g ( ν ) show the same asymptoticbehavior.In Figs. 11(b) and (d), we show decimal logarithms ofthe ratio G/g for the same K values used in Figs. 11(a)and (c). All graphs show that log[ Gg ( ν )] → as (cid:12)(cid:12) aν (cid:12)(cid:12) → , which is in agreement with formula (76). So G ( ν ) converges to g ( ν ) as ν increases, albeit this convergenceis restrained by increasing K .Another somewhat counterintuitive effect is related tothe tail thickness of g µ . Compared to the other distri-butions g µ> , g decays more slowly as ν increases, and log[ Gg ( ν )] decays more easily to zero. When µ increases , g µ ’s tails become thinner, and convergence of log[ Gg µ ( ν )] to zero requires larger values of ν . A simple explana-tion to this tail thickness effect is related to the criti-cal order parameter, which is defined by K ( µ ) c = πg µ (0) , if g = g µ (see Eq. (7)). Since the normalization con-dition remains valid for any µ , thinner tails result inhigher g (0) . If g µ (0) > . . . > g (0) > g (0) , then K ( µ ) c < . . . < K (2) c < K (1) c . So, for K fixed, the dif- ference K − K µc decreases with decreasing µ , R (and a )diminishes, and G resembles more closely g µ . When R is small, G ( ν ) converges more easily g µ ( ν ) as (cid:12)(cid:12) aν (cid:12)(cid:12) → .This is shown by Eq. (76) and Figs. 11(b) and (d). VI. CONCLUSION
Based on Kuramoto theory, we have obtained an an-alytical formulation of the instantaneous frequency dis-tribution in the Kuramoto model. Numerical data showexcellent agreement with our formula, provided they areobtained on very large collections of oscillators studiedin their steady state (see Appendix), i.e. in the limitswhere our results are expected to be valid. Access tothe distribution of instantaneous frequencies G extendsKuramoto theory, which was limited heretofore to theknowledge of G , the distribution of time-averaged, or“coupling-modified”, frequencies.Distributions G and G are functionals of the naturalfrequency distribution g . Irrespective of g , the synchro-nization scenario keeps the same basic features: beyondthe critical coupling strength value K c = πg (Ω) , a subsetof oscillators synchronize, so that both G and G com-prise a delta peak at the frequency Ω . This delta peakis identical for both G and G and represents the fractionof synchronized oscillators, which grows monotonouslywith K (at least for the g distributions considered here,see e.g. Fig. 4). As soon as K > K c , the continuous partof both G and G departs from g . Whereas G remains be-low g everywhere, and displays a characteristic dip nearthe synchronized frequency Ω , the continuous part of G has tails that pass over g and shows a maximum at Ω for K large enough. This last fact reflects the relaxation-oscillation-like dynamics of oscillators with natural fre-quency outside but close to the synchronization band.The distribution of instantaneous frequencies G typ-ically displays a rather complicated structure. It is infact trimodal for K large enough, even as the natural fre-quency distributions considered here are unimodal. Nev-ertheless, from the 3 qualitatively-different examples of g studied here, we have shown that G , like G , displays thesame tails as g : a normal g yields Gaussian tails for G and G which are just rescaled versions of those of g . Fora Beta g with a bounded support interval, G and G bothhave bounded supports, respectively wider and narrowerthan that of g .Due to the difficulty in extracting explicit expressionsfrom our results, most of the information about the dis-tribution of instantaneous frequencies presented here hasbeen obtained by numerical analysis of our formula. Yet,when dealing with rare events in the example of powerlawtailed distribution of natural frequencies, a power-seriesanalysis of the distribution of instantaneous frequencieshas allowed us to obtain an asymptotic expansion in fre-quency. This has confirmed that g , G , and G are asymp-totically equivalent in the limit of large frequencies.Beyond their intrinsic interest for a deeper understand-5 G D () g ( )K = 2.001K = 2.01K = 2.1K = 2.5 (a) µ = 1 G D () g ( )K = 2.001K = 2.01K = 2.1K = 2.5 (b) µ = 1 G D () g ( )K = 1.3334K = 1.3339K = 1.34K = 1.37 (c) µ = 3 G D () g ( )K = 1.3334K = 1.3339K = 1.34K = 1.37 (d) µ = 3 FIG. 10. G D and G for g = g µ , defined by (77). ing of synchronization, our results are useful when itcomes to choosing a numerical scheme and resolution tosimulate coupled oscillators: indeed a faithful simulationmust account properly for the largest instantaneous fre-quencies displayed by the system. As seen and quantifiedhere, these are larger than the largest natural frequencypresent, which implies, e.g., to choose higher-order inte-gration schemes and/or smaller timesteps than naivelysuggested by the natural frequencies at play.The approach followed here can easily be extendedto non-symmetric and/or non-unimodal distributions ofnatural frequencies. We also believe that important vari-ants of the Kuramoto model, such as the Kuramoto-Sakaguchi model [15] are amenable to the same type ofanalysis as developed here. More generally, we hope thatthis work opens new perspectives on synchronization phe-nomena beyond the usual order-parameter analysis. ACKNOWLEDGMENTS
This work was made possible through financial sup-port from Brazilian research agency FAPESP (grantn. 2019/12930-9 ). JDF warmly thanks Prof. JoaoPeres for valuable discussions. EDL thanks support fromBrazilian agencies CNPq (301318/2019-0) and FAPESP(2019/14038-6).
Appendix: Numerical checks
In this section we validate the formula of G D againstnumerical results from simulations of the Kuramotomodel. Our check is restricted to the Gaussian and Betaexamples discussed of Sec. IV. Simulations were per-formed with the numerical library ODEPACK [28]. The6 l o g [ G D () / g ()] K = 3.0K = 4.0K = 5.0K = 6.0 (a) µ = 1 l o g ( G D () / g ()) K = 3.0K = 4.0K = 5.0K = 6.0 (b) µ = 1 l o g [ G D () / g ()] K = 3.0K = 4.0K = 5.0K = 6.0 (c) µ = 3 l o g ( G D () / g ()) K = 3.0K = 4.0K = 5.0K = 6.0 (d) µ = 3 FIG. 11. Decimal logarithms of the ratios G D g ( ν ) and Gg ( ν ) . g ( ν ) = g µ ( ν ) , defined by Eq. (77). ODEPACK’s solver used in all simulations was
LSODA ,an hybrid implementation of Adams and BDF methods[29].
1. Gaussian case
In Figs. 12(a)-(f), we compare graphs of G D , the sameshown in Fig. 4, to numerically-obtained normalized his-tograms of the instantaneous frequencies.The histograms were obtained from numerical simu-lations of the Kuramoto model with N = 5 × . Inall simulations, the Kuramoto system of equations is nu-merically integrated from a initial time t = 0 to a finaltime t f = 5 × . Each histogram is created from theset of instantaneous frequencies { ˙ θ i ( t f ) } Ni =1 . Simulations are performed considering random samples of natural fre-quencies and initial phases. Initial phases are uniformlysampled: a sample { θ i ( t ) } Ni =1 is generated according toa uniform distribution in the interval between and π .Figure 13 shows the typical evolution of the numeri-cal order parameter. When the order parameter exhibitssmall fluctuations after a sufficiently long time, the corre-sponding histograms are in good agreement with the an-alytical curves. However, stronger order parameter fluc-tuations are observed near the transition ( K = 1 . ),and the histogram shown in Fig. 12(b), obtained withthe same value of K , does not fit properly the curve.Figures 14(a) and 14(b) show, in logarithmic scale, thegraph of G D and histograms of instantaneous frequenciesobtained from simulation data. The coupling strengthhas the same value used in Fig. 12(d), K = 1 . . Thehistograms are created with N = 5 × (Fig. 14(a))7 G D ( )histogram (a) K = 0 . G D ( )histogram (b) K = 1 . G D ( )histogram (c) K = 1 . G D ( )histogram (d) K = 1 . G D ( )histogram (e) K = 1 . G D ( )histogram (f) K = 1 . FIG. 12. Comparison between normalized histograms of instantaneous frequencies (in blue), obtained from numerical simula-tions of the Kuramoto model, and the curve of G D (in red). In all simulation, N = 5 × . time R K = 0.8K = 1.596K = 1.61K = 1.64K = 1.67K = 1.8
FIG. 13. Time-dependence of the numerical order parameter.Except near the transition, where K = 1 . , low fluctuationsare observed after a sufficiently long time. Created with thesame simulation near the transition, histogram of Fig. 12(b)does not fit properly the theoretical curve. and N = 5 × (Fig. 14(b)). In both of them, the syn-chronization peak is clearly visible. Large instantaneousfrequency occurrences (rare events) are more difficult toobserve. However, by increasing the number of oscilla- tors, rare events are more common, and the tails of G D fit better the histograms.
2. Beta case
In Figs. 15(a)-(f), we compare instantaneous fre-quency histograms to graphs of G D with g defined asthe Beta (2 , distribution, given by Eq. (58). Simula-tions were performed with N = 1 × oscillators, andintegration time was × . Histograms were createdby using the instantaneous frequency set { ˙ θ i ( t f ) } Ni =1 .Corresponding time series of the numerical order pa-rameter are shown in Fig. 16. Similarly to the case ofnormally-distributed natural frequencies, stronger order-parameter fluctuations are observed near the transition( K = 0 . ), away from which our analytical resultdescribes better the histograms.Behavior near the tails is shown in Fig. 17, wherewe compare again G D to instantaneous frequency his-tograms. The vertical axis has logarithmic scale for val-ues greater than − and linear scale between and − . We use the same coupling strength as in Fig. 15(d).As expected: i) rare events cannot be easily observed inthe histograms; ii) by increasing the number of oscilla-tors from N = 1 × to N = 2 × , these events aremore common, and a better fit is attainable in the tails.8 G D histogram (a) N = 5 × ; K = 1 . G D histogram (b) N = 5 × ; K = 1 . FIG. 14. Comparison (in logarithmic scale) between the graphof G D and histograms of instantaneous frequencies. The tailsof G D fit better the histograms when network size increasesfrom N = 5 × to N = 5 × . G D ( )histogram (a) K = 0 . G D ( )histogram (b) K = 0 . G D ( )histogram (c) K = 0 . G D ( )histogram (d) K = 0 . G D ( )histogram (e) K = 0 . G D ( )histogram (f) K = 0 . FIG. 15. Comparison of G D graphs to instantaneous frequency histograms. g is a Beta (2 , distribution, and the histogramswere obtained from simulations of the Kuramoto model with N = 1 × . The graphs provide better fits than near thetransition ( K = 0 . ). time R K = 0.2K = 0.42442K = 0.43K = 0.435K = 0.45K = 0.5
FIG. 16. Time-dependence of the numerical order parameter:low fluctuations after long time, except near the synchroniza-tion transition. G D frequency histogram (a) N = 10 ; K = 0 . G D frequency histogram (b) N = 2 × ; K = 0 . FIG. 17. Graph of G D and instantaneous frequency his-tograms. For the vertical axis, we use logarithmic scale forvalues greater than − and linear scale those between and − . By increasing the number of oscillators, a better fit isattainable in the tails.1