Bifurcations in the time-delayed Kuramoto model of coupled oscillators: Exact results
aa r X i v : . [ n li n . AO ] M a y Noname manuscript No. (will be inserted by the editor)
Bifurcations in the time-delayed Kuramoto model ofcoupled oscillators: Exact results
David M´etivier · Shamik Gupta
Received: date / Accepted: date
Abstract
In the context of the Kuramoto model of coupled oscillators with dis-tributed natural frequencies interacting through a time-delayed mean-field, wederive as a function of the delay exact results for the stability boundary betweenthe incoherent and the synchronized state and the nature in which the latter bifur-cates from the former at the critical point. Our results are based on an unstablemanifold expansion in the vicinity of the bifurcation, which we apply to boththe kinetic equation for the single-oscillator distribution function in the case of ageneric frequency distribution and the corresponding Ott-Antonsen(OA)-reduceddynamics in the special case of a Lorentzian distribution. Besides elucidating theeffects of delay on the nature of bifurcation, we show that the approach due to Ottand Antonsen, although an ansatz, gives an amplitude dynamics of the unstablemodes close to the bifurcation that remarkably coincides with the one derived fromthe kinetic equation. Further more, quite interestingly and remarkably, we showthat close to the bifurcation, the unstable manifold derived from the kinetic equa-tion has the same form as the OA manifold, implying thereby that the OA-ansatzform follows also as a result of the unstable manifold expansion. We illustrate ourresults by showing how delay can affect dramatically the bifurcation of a bimodaldistribution.
Keywords
Nonlinear dynamics and chaos · Synchronization · Coupled Oscilla-tors · Bifurcation analysis
David M´etivierCenter for Nonlinear Studies and Theoretical Division T-4 of Los Alamos National Laboratory,NM 87544, USATel.: +1 (505) 667-7052E-mail: [email protected] GuptaDepartment of Physics, Ramakrishna Mission Vivekananda University, Belur Math, Howrah711202, IndiaE-mail: [email protected] David M´etivier, Shamik Gupta
The Kuramoto model enjoys a unique status in the field of nonlinear dynamics [1, 2,3, 4, 5, 6, 7]. It provides arguably the minimal framework to model the phenomenonof spontaneous synchronization commonly observed in nature [8, 9], for example,among groups of fireflies flashing on and off in unison [10], in cardiac pacemakercells [11], in electrochemical [12] and electronic [13] oscillators, in Josephson junc-tion arrays [14] and in electrical power-grid networks [15], to name a few. Modelingthe individual units in a synchronizing system as nearly-identical limit-cycle os-cillators, the Kuramoto model arises on considering the individual limit-cycles tobe interacting weakly with one another, with the strength of coupling being thesame for every pair of oscillators [8, 7]. Specifically, the model comprises a set of N oscillators with distributed natural frequencies ω j ∈ [ −∞ , ∞ ]; j = 1 , , . . . , N .Denoting by θ j ∈ [0 , π ) the phase of the j -th oscillator, the dynamics of the modelis given by a set of N coupled equations of the form [1]d θ j ( t )d t = ω j + KN N X k =1 sin( θ k ( t ) − θ j ( t )) , (1)where K ≥ K by N ensures that theassociated term is well behaved in the thermodynamic limit N → ∞ , which isindeed the limit of interest to us in this paper and which also models the fact thatin this limit, the coupling between any pair of oscillators is weak, scaling as 1 /N .The frequencies { ω j } ≤ j ≤ N denote a set of quenched disordered random variablesdistributed according to a common distribution G ( ω ), with the latter obeying thenormalization R ∞−∞ d ω G ( ω ) = 1.The coherence among the phases of the oscillators in the Kuramoto model isconveniently measured by the complex order parameter r ( t ), defined as r ( t ) ≡ N N X j =1 e iθ j ( t ) , (2)where the quantity | r | measures the amount of coherence. A perfectly synchronizedstate corresponds to | r | = 1, while in an incoherent state, when the phases θ j areuniformly distributed in [0 , π ), one has | r | = 0. In terms of the order parameter,one may rewrite the equation of motion (1) asd θ j ( t )d t = ω j + K Im h r ( t ) e − i ( θ j ( t ) ) i , (3)which makes it evident that the phase of an oscillator at a given time evolves dueto the effect of the instantaneous mean field r built up in the system from theinteraction of the oscillator with all the other oscillators.The setting of the Kuramoto model is really the simplest one can conceive.It has indeed been very successful in addressing theoretically the emergence ofsynchrony in diverse dynamical settings involving coupled oscillators, showing inparticular in the limit N → ∞ and for a given distribution of natural frequenciesof the oscillators that the system exhibits a spontaneous transition from incoher-ence to synchrony as the coupling strength K is increased beyond a certain criticalvalue [1, 2]. Nevertheless, one does encounter situations where due to time delay ifurcations in the time-delayed Kuramoto model of coupled oscillators 3 in propagation of signals between the interacting units, the evolution of the dy-namical variables is determined not entirely by their instantaneous values but hasan essential dependence also on the values of the dynamical variables at an earliertime instant. Time delay is known to have important consequences, for example,for synchronization of biological clocks [16], in networks of digital phase-lockedloops [17], and in the case of information propagation through a neural networkin which delay is known to influence the temporal characteristics of oscillatorybehavior of neural circuits [18]. It has been argued that presence of even a smalltime delay may affect in general the global dynamics of ensembles of limit-cycleoscillators [19].To assess the effects of time delay, a pioneering work within the ambit of theKuramoto model was pursued by Yeung and Strogatz in Ref. [20], in which theyconsidered a generalization that allows for a time-delayed mean-field interactionbetween the oscillating units. Specifically, in terms of a time delay τ >
0, thephases evolve in time according to the dynamicsd θ j ( t )d t = ω j + KN N X k =1 sin [ θ k ( t − τ ) − θ j ( t ) − α ] . (4)In terms of r ( t ), the dynamics (4) takes the formd θ j ( t )d t = ω j + K Im h r ( t − τ ) e − i ( θ j ( t )+ α ) i , (5)which puts in evidence the presence of a delayed mean-field: the time evolutionof oscillator phases at time t is affected by the value of the mean-field r ( t − τ )measuring the macroscopic state of the system at a previous instant t − τ . Here, α ∈ ( − π/ , π/
2) is the so-called phase frustration parameter [21]. Note that settingthe phase frustration parameter and the time delay to zero, α = 0 , τ = 0, reducesthe dynamics (4) to that of the Kuramoto model (1).We may anticipate that introducing delay in the Kuramoto model may lead toa richer and more complex dynamical scenario, which would at the same time betheoretically more challenging to analyze. Indeed, Ref. [20] unraveled a range ofnew phenomena including bistability between synchronized and incoherent statesand unsteady solutions with time-dependent order parameters that do not occurin the original Kuramoto model (1). In the particular case of oscillators withidentical natural frequencies and with α = 0, Yeung and Strogatz were able toderive exact formulas for the stability boundaries between the incoherent andsynchronized states. For the general case of oscillators with distributed naturalfrequencies, still with α = 0, they adduced numerical results for the case of aLorentzian frequency distribution to suggest the occurrence of bifurcation of theincoherent state as a function of K , which could be either sub- or supercriticaldepending on the precise value of the delay parameter τ . A relevant work thatfollowed the work of Yeung and Strogatz is Ref. [22], which obtained the regionsof parameter space corresponding to synchronized and incoherent solutions, but forparticular frequency distributions. The complex dynamical scenario did not allowa straightforward analytical treatment to answer the following obvious questionsthat may be raised for generic frequency distributions : Can one predict analyticallyas a function of the time delay the nature of bifurcation of the incoherent statethat one would observe on changing the value of the coupling constant from low David M´etivier, Shamik Gupta to high values? What is the critical value of the coupling constant at which theincoherent state loses its stability? What are the effects of the phase frustrationparameter α ? Starting from the aforementioned developments and with an aim toanswer the questions just raised, in this work, we embark on a detailed analyticalcharacterization of bifurcation in model (4).In this work, we will consider G ( ω ) to be a distribution symmetric about itscenter given by ω = ω . The quantity ω coincides with the mean of G ( ω ) for casesin which the latter exists. Note that unlike (1), the dynamics (4) is not invariantunder the transformation θ j ( t ) → θ j ( t ) − ω t, ω j → ω j − ω ∀ j that tantamounts toviewing the dynamics in a frame rotating uniformly at frequency ω with respectto an inertial frame. From Eq. (4), it is evident that viewing the dynamics in sucha frame is equivalent to replacing α with α ′ ≡ α − ω τ . From now on, we willconsider such a choice of the reference frame and consequently the dynamicsd θ j ( t )d t = ω j + KN N X k =1 sin [ θ k ( t − τ ) − θ j ( t ) − α + ω τ ]= ω j + K Im h r ( t − τ ) e − i ( θ j ( t )+ α − ω τ ) i , (6)where the ω j ’s are to be regarded as distributed according to a distribution g ( ω )that is centered at zero: g ( ω ) ≡ G ( ω − ω ).The dynamics (6) in the limit N → ∞ is described by a kinetic equationfor the time evolution of the single-oscillator distribution function F ( θ, ω, t ) thatcounts the fraction of oscillators with natural frequency ω that have their phaseequal to θ at time t . The kinetic equation turns out to be an infinite-dimensionaldelay differential equation (DDE), which has the incoherent state as its stationarysolution existing for all values of K . In order to answer our queries raised above,namely, the bifurcation from the incoherent state that occurs as K is increasedto high values, as a first step, we specialize to a Lorentzian distribution for thefrequencies, and employ the celebrated Ott-Antonsen (OA) ansatz to deduce fromthe kinetic equation a simpler DDE satisfied by the order parameter r ( t ).The Ott-Antonsen (OA) approach offers a powerful exact method to study thedynamics of coupled oscillator ensembles [23, 24]. The approach allows to rewrite inthe thermodynamic limit the dynamics of coupled networks of phase oscillators interms of a few collective variables. In the context of the Kuramoto model (1) witha Lorentzian distribution of the oscillator frequencies, the ansatz studies the evo-lution in phase space by considering in the space D of all possible single-oscillatordistribution functions F ( θ, ω, t ) a particular class defined on and remaining con-fined to a manifold M in D under the time evolution of the phases. The OA ansatzobtains for this particular class of F ( θ, ω, t ) a single first-order ordinary differentialequation for the evolution of the synchronization order parameter r ( t ). The mostremarkable result of the approach, which explains its power and its widespreadapplicability in studying oscillator ensembles, is its ability to capture preciselyand quantitatively through this single equation all, and not just some, of the orderparameter attractors and bifurcations of the Kuramoto dynamics (which may beobtained either by performing numerical integration of the N coupled non-linearequations (1) for N ≫ r ( t ) in numerics, or, by simulating thekinetic system [3, 25, 26]), for a Lorentzian g ( ω ). The ansatz has since its proposi- ifurcations in the time-delayed Kuramoto model of coupled oscillators 5 tion been successfully applied to a variety of setups involving coupled oscillators.A few recent contributions are Refs. [27, 28, 29, 30, 31, 32, 33].Starting with the derived OA-ansatz-reduced dynamics for the order parame-ter, we perform both a linear and a nonlinear stability analysis of the incoherentstate, based on a general formalism to treat DDE [34, 35]. The linear analysis lo-cates the critical threshold K c above which a synchronized state bifurcates fromthe linearly-unstable incoherent state. On the other hand, the nonlinear analy-sis, developed in the spirit of Refs. [36, 37, 38, 39, 40, 41], treats the dynamical flowon the unstable manifold passing through K c , thereby obtaining the amplitudedynamics of the linearly unstable modes for K > K c . The form of amplitude dy-namics describing the flow in the regime of weak linear instability, namely, as K → K + c , allows to directly obtain the nature of bifurcation from the incoherentstate occurring as soon as K is increased beyond K c . In the second part of ouranalysis, we relax the choice of a Lorentzian g ( ω ) and perform the linear and non-linear stability analysis directly on the kinetic equation for F ( θ, ω, t ) to derive theamplitude dynamics of the linearly unstable modes close to the bifurcation of theincoherent state, that is, as K → K + c . Remarkably, we find that the amplitudeequation derived in the general case has the same form as that obtained fromthe OA-ansatz-reduced dynamics derived for Lorentzian g ( ω ), thus confirming thepower and general applicability of the OA ansatz. Viewing the analysis of the OA-ansatz-reduced dynamics vis-`a-vis that of the kinetic equation, we may commentthat the former is after all based on an ansatz that works only for Lorentzian g ( ω ),but allows to derive an amplitude dynamics that works for any K > K c . On theother hand, the analysis based on the kinetic equation is more general in the sensethat there is no ansatz and no specific choice of g ( ω ) involved, but a disadvan-tage is that the derived amplitude dynamics is valid only close to the bifurcation,namely, for K → K + c . Quite interestingly and remarkably, we find that close tothe bifurcation, the unstable manifold derived from the kinetic equation has thesame form as the OA manifold, so we may say that the OA-ansatz form may beobtained also as a result of the unstable manifold expansion.Furthermore, we demonstrate with our exact results a remarkable effect oftime delay: Considering a sum of two Lorentzians as a representative example ofa bimodal distribution, it is known that in the absence of delay, the bifurcation ofthe synchronized from the incoherent state is subcritical [42]. We show here thatpresence of even a small amount of delay suffices to completely change the natureof the bifurcation and make it supercritical !The paper is organized as follows. In Section 2, we discuss the characterizationof the dynamics (6) in the thermodynamic limit in terms of a kinetic equationfor the single-oscillator distribution function F ( θ, ω, t ). In Section 3, we apply forthe special case of a Lorentzian distribution for the natural frequencies the Ott-Antonsen ansatz to replace the kinetic equation with an equation for the order pa-rameter. The linear and nonlinear stability analysis of the incoherent state is thenpursued in detail in Section 5, based on a formalism for treating DDE summarizedin Section 4. A similar analysis for the case of general frequency distributions andapplied directly to the kinetic equation is taken up in Section 6, obtaining resultssimilar to those in Section 5 for Lorentzian g ( ω ). The paper ends with conclusionsand perspectives. David M´etivier, Shamik Gupta
In the limit N → ∞ , the state of the oscillator system (6) at time t may be char-acterized by the time-dependent single-oscillator distribution function F ( θ, ω, t ),defined such that F ( θ, ω, t )d θ gives the fraction of oscillators with frequency ω thathave their phase in [ θ, θ + d θ ] at time t . The function F ( θ, ω, t ) is 2 π -periodic in θ ,i.e., F ( θ + 2 π, ω, t ) = F ( θ, ω, t ), and is normalized as Z π d θ F ( θ, ω, t ) = g ( ω ) ∀ t. (7)Since the dynamics (6) conserves in time the number of oscillators with a givennatural frequency, the function F ( θ, ω, t ) evolves in time according to a kineticequation given by the continuity equation ∂F/∂t + ( ∂/∂θ )( F d θ/ d t ) = 0, which onusing Eq. (6) yields (for a rigorous derivation, see Ref. [43]) ∂F ( θ, ω, t ) ∂t + ω ∂F ( θ, ω, t ) ∂θ + K i ∂∂θ h(cid:16) r [ F ]( t − τ ) e − i ( θ + α − ω τ ) − r ∗ [ F ]( t − τ ) e i ( θ + α − ω τ ) (cid:17) F ( θ, ω, t ) i = 0 . (8)Here and in the following, ∗ stands for complex conjugation, and we have definedas functionals of F the quantity r [ F ]( t ) ≡ Z d θ d ω e iθ F ( θ, ω, t ) (9)as the thermodynamic limit of Eq. (2).Equation (8) is an infinite dimensional DDE. In the following section, we em-ploy the so-called Ott-Antonsen ansatz [23, 24] that allows to derive from Eq. (8)a DDE for the order parameter. The implementation of the Ott-Antonsen ansatz for the dynamics (6) has beendiscussed in Ref. [23], which we briefly recall here for use in later parts of thepaper. As is usual with OA ansatz implementation, we will make the specificchoice of a Lorentzian distribution for g ( ω ): g ( ω ) = ∆π ω + ∆ , (10)where ∆ > g ( ω ). Consider the func-tion F ( θ, ω, t ), which being 2 π -periodic in θ may be expanded in a Fourier seriesin θ : F ( θ, ω, t ) = g ( ω )2 π h ∞ X n =1 (cid:16) e F n ( ω, t ) e inθ + [ e F n ( ω, t )] ∗ e − inθ (cid:17)i , (11)where e F n ( ω, t ) is the n -th Fourier coefficient. Using R π d θ e inθ = 2 πδ n, , we checkthat the above expansion is consistent with Eq. (7). ifurcations in the time-delayed Kuramoto model of coupled oscillators 7 The OA ansatz considers in the expansion (11) a restricted class of Fouriercoefficients given by [23, 24] e F n ( ω, t ) = [ z ( ω, t )] n , (12)with z ( ω, t ) an arbitrary function with the restriction | z ( ω, t ) | < z ( ω, t ) may be analytically continued to the whole of thecomplex- ω plane, that it has no singularities in the lower-half complex- ω plane,and that | z ( ω, t ) | → ω ) → −∞ [23, 24]. Using Eqs. (11) and (12) in Eq. (9),one gets r ( t ) ≡ r [ F ]( t ) = Z ∞−∞ d ω g ( ω ) z ∗ ( ω, t ) . (13)On substituting Eqs. (11), (12), and (13) in Eq. (8) and on collecting and equatingthe coefficient of e inθ to zero, we get ∂z ( ω, t ) ∂t + iωz ( ω, t ) + K h e − i ( α − ω τ ) r ( t − τ ) z ( ω, t ) − e i ( α − ω τ ) r ∗ ( t − τ ) i = 0 . (14)Note that even if the initial condition does not respect the OA ansatz, it has beenshown that the OA manifold is attractive for the dynamics (8) with τ = 0, (see,for example, Ref. [44]).For the Lorentzian g ( ω ), Eq. (10), one may evaluate r ( t ) by using Eq. (13) toget r ( t ) = 12 iπ I C d ω z ∗ ( ω, t ) (cid:20) ω − i∆ − ω + i∆ (cid:21) = z ∗ ( − i∆, t ) , (15)where the contour C consists of the real- ω axis closed by a large semicircle in thelower-half complex- ω plane on which the integral in Eq. (15) gives zero contributionin view of | z ( ω, t ) | → ω ) → −∞ . The second equality in Eq. (15) is obtainedby applying the residue theorem to evaluate the complex integral over the contour C . Using Eqs. (14) and (15), we finally obtain the OA equation for the timeevolution of the synchronization order parameter as the DDE [23]d r ( t )d t + ∆r ( t ) + K h e − i ( α − ω τ ) r ∗ ( t − τ ) r ( t ) − e i ( α − ω τ ) r ( t − τ ) i = 0 , (16)whose solution requires as an initial condition the value of r ( t ) over an entireinterval of time t , namely, t ∈ [ − τ, τ = 0, Eq. (16) is a finite-dimensional ODE for r ( t ) that requires for its solution only the value of r ( t ) at t = 0 as an initial condition. In this case, it has been demonstrated that this singleequation contains all the bifurcations and attractors of the order parameter (2)as obtained through the evolution of the set of coupled differential equations (1)for a Lorentzian g ( ω ) and in the limit N → ∞ . In order to effectively describe thebifurcation of the dynamics (16) between an incoherent and a synchronized statein the case τ = 0, we need to put the system (16) into the so-called normal form,by reducing the dimensionality of the evolution (16) to get a simple ODE. Wetake up this program in Section 5 based on a general formalism for DDE that wesummarize in the next section. David M´etivier, Shamik Gupta
It is evident from the form of equations (8) and (16) that both may be cast in thegeneral form ∂∂t H ( t ) = M [ H ( t ) , H t ( − τ )] , (17)where the function H ( t ) is either the function F ( θ, ω, t ) or the function r ( t ) depend-ing on the case of interest. Here, we have introduced the notation H t ( ϕ ) ≡ H ( t + ϕ )for − τ ≤ ϕ ≤
0. To solve for H ( t ); t > H t ( ϕ ); − τ ≤ ϕ ≤
0. The time evolution of H t ( ϕ ) for − τ ≤ ϕ < ∂H t ( ϕ ) /∂t = lim δ → ( H t + δ ( ϕ ) − H t ( ϕ )) /δ =lim δ → ( H ( t + ϕ + δ ) − H ( t + ϕ )) /δ = d H ( t + ϕ ) / d ϕ = d H t ( ϕ ) / d ϕ . We may thenquite generally write ∂∂t H t ( ϕ ) = ( A H t )( ϕ ); − τ ≤ ϕ ≤ , (18)where we have ( A H t )( ϕ ) = dd ϕ H t ( ϕ ); − τ ≤ ϕ < , M [ H t ( ϕ ) , H t − τ ( ϕ )]; ϕ = 0 . (19)Around a stationary state H st of Eq. (18), we define the perturbation h t ( ϕ ) as H t ( ϕ ) ≡ H st + h t ( ϕ ). The time evolution of h t ( ϕ ) has the same form as Eq. (18),in which we now split the operator A into a part D that is linear in h and a part F that is nonlinear, to write( A h t )( ϕ ) = ( D h t + F [ h t ])( ϕ )= dd ϕ h t ( ϕ ) L h t ( ϕ ) + ( − τ ≤ ϕ < , N [ h t ]; ϕ = 0 . (20)We will find it convenient to decompose the linear operator L into a part L thatdoes not contain any term involving the delay and a part R containing the delayterms: L h t ( ϕ ) = L h t (0) + R h t ( − τ ).The adjoint of the linear operator D may be defined through the definition ofthe scalar product [34]( h t , h t ) τ ≡ ( h t (0) , h t (0)) + Z − τ d ξ ( h t ( ξ + τ ) , R h t ( ξ )) , (21)where the scalar product in the absence of any delay ( · , · ) is either( s, r ) = s ∗ r, with h (0) = s, h (0) = r (22)for the finite dimensional case, or, the standard L ( T × R ) scalar product( h, f ) = Z T × R h ∗ ( θ, ω ) f ( θ, ω ) d ω d θ, with h (0) = h ( θ, ω ) , h (0) = f ( θ, ω ) , (23)for the kinetic case (infinite dimensional). The adjoint operator D † , which satisfies( h t ( ϕ ) , D h t ( ϕ )) τ = ( D † h t ( ϕ ) , h t ( ϕ )) τ , (24) ifurcations in the time-delayed Kuramoto model of coupled oscillators 9 is obtained as( D † h t )( ϑ ) = − dd ϑ h t ( ϑ ); 0 < ϑ ≤ τ, L † h t ( ϑ ) = L † h t (0) + R † h t ( τ ); ϑ = 0 . (25)Note that h ( ϕ ) and h ( ϑ ) belong to different functional spaces, for example, ϕ ∈ [ − τ,
0] and ϑ ∈ [0 , τ ]. We now apply the formalism of the DDE discussed in the preceding section toanalyze Eq. (16), with the aim to reduce this infinite-dimensional equation to afinite-dimensional equation. The incoherent stationary state r st = 0 is always asolution of Eq. (16). It is of interest to study its stability as K is varied, whichmay be done by studying the behavior of perturbations about r st = 0, where theperturbations r t ( ϕ ) satisfydd t r t ( ϕ ) = ( D r t + F [ r t ])( ϕ ); − τ ≤ ϕ ≤ , (26)with ( D r t )( ϕ ) = dd ϕ r t ( ϕ ); − τ ≤ ϕ < , L r t ( ϕ ) = L r t (0) + R r t ( − τ ); ϕ = 0 , (27)and ( F [ r t ])( ϕ ) = ( − τ ≤ ϕ < , N [ r t ]; ϕ = 0 . (28)The adjoint of the linear operator D is given by( D † s t )( ϑ ) = − dd ϑ s t ( ϑ ); 0 < ϑ ≤ τ, L † s t ( ϑ ) = L † s t (0) + R † s t ( τ ); ϑ = 0 . (29)In these equations, we haveL r = − ∆r, (30)R r = K e i ( α − ω τ ) r, (31)L † r = − ∆r, (32)R † r = K e − i ( α − ω τ ) r, (33) N [ r t ] = − K e − i ( α − ω τ ) r ∗ t ( − τ ) r t (0) . (34)Small perturbations r t ( ϕ ) may be expressed as a linear combination of theeigenfunctions of the linear operator D . To this end, let us then solve the eigen-function equation ( D p )( ϕ ) = λp ( ϕ ) (35) for − τ ≤ ϕ <
0; we get p ( ϕ ) = p (0) e λϕ . The equation ( D p )( ϕ ) = λp ( ϕ ) for ϕ = 0gives λp (0) = − ∆p (0) + ( K/ e i ( α − ω τ ) p (0) e − λτ . With p (0) = 0, we thus get thedispersion relation Λ ( λ ) = λ + ∆ − K e − λτ + i ( α − ω τ ) = 0 . (36)The solution of the above equation gives the discrete eigenvalues λ l (with l ∈ Z )in terms of the Lambert-W function W l , as λ l = − ∆τ + W l (cid:18) Kτ e iα + ∆τ − iω τ (cid:19) τ . (37)Without loss of generality, we may take p (0) = 1. We thus conclude that p ( ϕ ) = e λϕ is an eigenfunction of the linear operator D for − τ ≤ ϕ ≤ λ provided that λ satisfies Λ ( λ ) = 0. In other words, a discrete set of eigenvaluescorrespond to D for all values of ϕ . Perturbations r t ( ϕ ) may be expressed as alinear combination of the corresponding eigenfunctions. It then follows that thestationary solution r st = 0 will be linearly stable under the dynamics (26) so longas all the eigenvalues λ have a real part that is negative. Vanishing of the real partof the eigenvalue with the smallest real part then signals criticality above which r st = 0 is no longer a linearly-stable stationary solution of Eq. (26). Denoting by λ i ; λ i ∈ R the imaginary part of the eigenvalue with the smallest real part, wethus have at criticality the following equations obtained from Eq. (36): K c α − ( ω + λ i ) τ ) = ∆, (38) K c α − ( ω + λ i ) τ ) = λ i . We want to study the behavior of r t ( ϕ ) as K → K + c , the goal being to uncoverthe weakly nonlinear dynamics occurring beyond the exponential growth takingplace due to the instability as K → K + c . To this end, we want to study thebehavior of r t ( ϕ ) on the unstable manifold, which by definition is tangential to theunstable eigenspace spanned by the eigenfunctions p ( ϕ ) at the equilibrium point( K = K c , λ = iλ i ). This manifold may be shown to be an attractor of the dynamicsfor the type of DDE under consideration [35, 45, 46] and is therefore of interest tostudy. To proceed, we need the eigenfunctions of the adjoint operator D † , whichwill be useful in discussing the unstable manifold expansion. It is easily checkedthat D † has the eigenfunction q ( ϑ ) = q (0) e − λ ∗ ϑ associated with the eigenvalue λ ∗ satisfying Λ ∗ ( λ ∗ ) = 0, that is, we get the same dispersion relation as for D . Wemay choose q (0) such that ( q ( ϕ ) , p ( ϕ )) τ = 1. Using Eq. (21), and noting that inthe present case, ( q (0) , p (0)) = q ∗ (0) p (0), we get q ∗ (0) + Z − τ d ξ q ∗ (0) e − λ ( ξ + τ ) K e i ( α − ω ) τ e λξ = 1 , (39)yielding q ∗ (0) = 11 + τ K e − ( λ + iω ) τ + iα = 1 Λ ′ ( λ ) (40) ifurcations in the time-delayed Kuramoto model of coupled oscillators 11 where we chose without loss of generality p (0) = 1.The unstable manifold expansion of r t ( ϕ ) for K > K c reads r t ( ϕ ) = A ( t ) p ( ϕ ) + w [ A ]( ϕ ) , (41)where w [ A ]( ϕ ), which is at least quadratic in A (in fact, one can prove that itis cubic in A in the present case), denotes the component of r t ( ϕ ) transverse tothe unstable eigenspace, so that ( q ( ϕ ) , w ( ϕ )) τ = 0. On using the latter equation,together with ( q ( ϕ ) , p ( ϕ )) τ = 1 in Eq. (41), we get A ( t ) = ( q ( ϕ ) , r t ( ϕ )) τ . The timeevolution of A ( t ) is then obtained as˙ A = ( q ( ϕ ) , ˙ r t ( ϕ )) τ = ( q ( ϕ ) , ( D r t + F [ r t ])( ϕ )) τ = ( q ( ϕ ) , A ( t ) λp ( ϕ ) + D w ( ϕ ) + F [ r t ]( ϕ )) τ = λA + q ∗ (0) N [ r t ] , (42)where the dot denotes derivative with respect to time. Here, in arriving at the sec-ond and the third equality, we have used Eqs. (26) and (41), while in obtaining thelast equality, we have used in the third step ( q ( ϕ ) , D w ( ϕ )) τ = ( D † q ( ϕ ) , w ( ϕ )) τ = λ ∗ ( q ( ϕ ) , w ( ϕ )) τ = 0. Since we can prove that w is O ( | A | A ), while we see that N [ r t ] is of order three in r t , the leading-order contribution to the nonlinearterm on the right hand side of Eq. (42) is obtained as N [ r t ] = N [ A ( t ) p ( ϕ )] + O ( | A | A ) = − ( K/ e − i ( α − ω τ ) p ∗ ( − τ ) p (0) | A ( t ) | A ( t )] + O ( | A | A ). Using this re-sult and Eq. (40) in Eq. (42), we get eventually the so-called normal form for thetime evolution of A as˙ A = λA − K e ( iω − λ ∗ ) τ − iα τ K e − ( λ + iω ) τ + iα | A | A + O ( | A | A ) . (43)The above is the desired finite-dimensional ordinary differential equation corre-sponding to the infinite-dimensional equation (26), which allows to decide thebifurcation behavior of r t ( ϕ ) as K → K + c . The relevant parameter to study thetype of bifurcation is given by the sign of the second term on the right hand side.Denoting this term by c , we then need to study the sign of the real part of c asthe real part of λ approaches zero, so that λ = iλ i is purely imaginary:Re( c ) = − K c e i ( ω + λ i ) τ − iα τ K c e − i ( ω + λ i ) τ + iα . (44)Note that a similar approach was pursued in [47]. It is of interest to consider the formalism of the DDE and the unstable mani-fold expansion applied in the previous section to the OA-ansatz-reduced dynam-ics (16), and to apply it directly to the kinetic system, Eq. (8), so as to reduce thisinfinite-dimensional equation to a finite-dimensional one. The advantage of suchan application stems from the fact that while the OA-ansatz-reduced-dynamics is obtained for a Lorentzian g ( ω ), the kinetic system is valid for any generic g ( ω ). Adisadvantage is that the computations for the kinetic system are only formal inthe sense that standard theorems proving the attracting property of the unstablemanifold do not apply straightforwardly for the type of kinetic equation underconsideration because of the existence of the continuous spectrum (see below). yattractiveness of a manifold is meant that almost all trajectories during the courseof the dynamical evolution come at long times arbitrarily close to the manifold sothat their eventual evolution coincides with evolution of trajectories lying on themanifold itself. The reader is referred to Refs. [48, 49] for a mathematical proof ofthe statement in the non-delayed model, that, is, for τ = 0.Let us give at the outset an outline of the current rather technical section.We start off with rewriting of the kinetic equation (8), of which the incoherentstate F st represents a stationary solution, in the form of a DDE for perturbations f t ( ϕ ) around F st . The DDE involves a linear evolution operator D and a nonlinearone. In the next step, we obtain the eigenvalues and the eigenvectors of D andof the corresponding adjoint operator D † . As is usual with linear stability analy-sis [50], the knowledge of the eigenvalues allows to locate the critical value K c ofthe coupling K above which the incoherent state F st becomes linearly unstable anda synchronized stationary state bifurcates from it. Following this, we perform for K > K c an unstable manifold expansion of F st along the two unstable eigenvectors(representing respectively the eigenvectors of D and D † ) and the unstable manifoldpassing through K c . By invoking a convenient Fourier expansion of the relevantquantities and working at K slightly greater than K c , we obtain the amplitudedynamics describing the evolution of perturbations f t ( ϕ ) in the regime of weaklinear instability, K → K + c . The form of the amplitude dynamics allows to directlyobtain the nature of bifurcation occurring as soon as K is increased beyond K c :The amplitude dynamics has a leading linear term, followed by a nonlinear (cubic)term, and according to elements of bifurcation theory well known in the literature[50], it is the sign of this cubic term that dictates the precise nature of the bifur-cation (positive and negative leading respectively to subcritical and supercriticalbifurcation). We show explicitly how the sign of the cubic term varies as a functionof the delay τ , for two separate cases, namely, that of a unimodal and a bimodalLorentzian frequency distribution. A remarkable result of this section is that closeto the bifurcation, the unstable manifold that we derive in this section based onthe kinetic equation (8) has the same form as the Ott-Antonsen manifold discussedin Section 3. We now proceed to a detailed derivation of our results.Similar to the preceding section, we rewrite Eq. (8) in the form of Eq. (18).We consider perturbations f t ( ϕ ) around the incoherent stationary state of (8),namely, F st = g ( ω ) / (2 π ), so that F t ( ϕ ) = F st ( ϕ ) + f t ( ϕ ). From Eq. (7), it followsthat R π d θ f t ( ϕ ) = 0. Perturbations f t will evolve according todd t f t ( ϕ ) = ( D f t + F [ f t ])( ϕ ); − τ ≤ ϕ ≤ , (45)with ( D f t )( ϕ ) = dd ϕ f t ( ϕ ); − τ ≤ ϕ < , L f t ( ϕ ) = L f t (0) + R f t ( − τ ); ϕ = 0 , (46) ifurcations in the time-delayed Kuramoto model of coupled oscillators 13 and ( F [ f t ])( ϕ ) = ( − τ ≤ ϕ < , N [ f t ]; ϕ = 0 . (47)The adjoint of the linear operator D is given by( D † h t )( ϑ ) = − dd ϑ h t ( ϑ ); 0 < ϑ ≤ τ, L † h t ( ϑ ) = L † h t (0) + R † h t ( τ ); ϑ = 0 . (48)Here, we haveL f = − ω ∂∂θ f, (49)R f = K g ( ω )2 π (cid:16) r [ f ] e − i ( θ + α − ω τ ) + r ∗ [ f ] e i ( θ + α − ω τ ) (cid:17) , (50)L † h = ω ∂∂θ h, (51)R † h = K π (cid:16) r [ g ( ω ) h ] e − i ( θ − α + ω τ ) + r ∗ [ g ( ω ) h ] e i ( θ − α + ω τ ) (cid:17) , (52) N [ f t ] = − K i ∂∂θ (cid:16)(cid:16) r [ f t ]( − τ ) e − i ( θ + α − ω τ ) − r ∗ [ f t ]( − τ ) e i ( θ + α − ω τ ) (cid:17) f t (0) (cid:17) . (53)To study the linear stability of F st , similar to what was done in the precedingsection, we first solve the eigenfunction equation( D P )( ϕ ) = λP ( ϕ ) (54)for − τ ≤ ϕ <
0; we get P ( ϕ ) = Ψe λϕ for arbitrary Ψ . Since we will in the followingexpand f t ( ϕ ) in terms of P ( ϕ ), we would need to choose Ψ as Ψ ( θ, ω ), where 2 π -periodicity of f t implies that so should be Ψ ( θ, ω ). Consequently, we may expand Ψ ( θ, ω ) in a Fourier series in θ , as Ψ ( θ, ω ) = (2 π ) − P ∞ k = −∞ ψ k ( ω ) e ikθ , so that P ( ϕ ) = (2 π ) − P ∞ k = −∞ ψ k ( ω ) e ikθ e λϕ . Using the equation ( D P )( ϕ ) = λP ( ϕ ) for ϕ = 0 and k = ± P ( ϕ ), it may be easily seen withthe condition r [ Ψ ] = r ∗ [ Ψ ] = 1 that p ( ϕ ) = ψ ( ω ) e iθ + λϕ and p ∗ ( ϕ ) give twoindependent eigenfunctions of D with eigenvalues λ and λ ∗ , respectively, wherethe latter satisfy Λ ( λ ) = Λ ∗ ( λ ∗ ) = 0, and ψ ( ω ) = K e − λτ + i ( α − ω τ ) g ( ω ) λ + iω , (55) Λ ( λ ) = 1 − K e − λτ + i ( α − ω τ ) Z d ω g ( ω ) λ + iω . (56)For k = ±
1, one has a continuous spectrum sitting on the imaginary axis; thisfeature is characteristic of kinetic equations of the type of Eq. (8) [51, 36, 52]. For
K > K c , when the incoherent stationary state is linearly unstable, the unstableeigenspace is spanned by the eigenfunctions p ( ϕ ) and p ∗ ( ϕ ).The stationary solution F st = g ( ω ) / (2 π ) will be neutrally stable, thanks tothe continuous spectrum sitting on the imaginary axis that generates a dynamicssimilar to that of Landau damping [53, 20], when there are no eigenvalues λ with a positive real part. Vanishing of the real part of the eigenvalue with the smallest realpart then signals criticality above which F st becomes a linearly unstable stationarysolution of Eq. (26). Denoting by λ i ; λ i ∈ R the imaginary part of the eigenvaluewith the smallest real part, we thus have at criticality the following equationsobtained from Eq. (56):cos ( α − ( ω + λ i ) τ ) = g ( − λ i ) πK c , (57)tan ( α − ( ω + λ i ) τ ) = PV R g ( ω ) ω + λ i πg ( − λ i ) , where PV stands for principal value. One then has to have cos ( α − ( ω + λ i ) τ ) > D † are given by q ( ϑ ) = e ψ ( ω ) e iθ − λ ∗ ϑ and q ∗ ( ϑ ) with eigenvalues λ ∗ and λ , respectively, where we fix e ψ ( ω ) by requiringthat ( q ( ϕ ) , p ( ϕ )) τ = 1 = R d θ d ω q ∗ (0) p (0)+ R − τ d ξ R d θ d ω q ∗ ( ξ + τ ) R p ( ξ ). We thusget e ψ ( ω ) = (( Λ ′ ( λ )) ∗ ( λ ∗ − iω )2 π ) − , and hence, q ( ϑ ) = 1( Λ ′ ( λ )) ∗ λ ∗ − iω e iθ − λ ∗ ϑ . (58)Similar to Eq. (41), we now decompose perturbations f t ( ϕ ) along the twounstable eigenvectors and the unstable manifold, as f t ( ϕ ) = A ( t ) p ( ϕ ) + A ∗ ( t ) p ∗ ( ϕ ) + w [ A, A ∗ ]( ϕ ) , (59)with the relations ( q ( ϕ ) , p ( ϕ )) τ = 1 , ( q ( ϕ ) , p ∗ ( ϕ )) τ = 0 , ( q ( ϕ ) , w ( ϕ )) τ = 0 yielding A ( t ) = ( q ( ϕ ) , f t ( ϕ )) τ . We require w ( ϕ ) to be at least quadratic in A .Let us define the following Fourier expansion needed for further analysis: f t = 12 π ∞ X k = −∞ ( f t ) k e ikθ , (60) N [ f t ] = 12 π ∞ X k = −∞ N k [ f t ] e ikθ , (61) w = 12 π ∞ X k = −∞ w k e ikθ . (62)Using Eq. (53), we then get N k [ f t ] = − kK (cid:16) r [ f t ]( − τ ) e − i ( α − ω τ ) ( f t ) k +1 (0) − r ∗ [ f t ]( − τ ) e i ( α − ω τ ) ( f t ) k − (0) (cid:17) . (63)Note that we have ( f t ) = 0, so that Eq. (59) gives w = 0; this feature is a majordifference with respect to a similar kinetic equation, the Vlasov equation [39]. Bysymmetry on the unstable manifold, we have [39, 46] w = O ( | A | A ) and for k > ifurcations in the time-delayed Kuramoto model of coupled oscillators 15 w k = O ( A k ). From Eq. (59), we get ( f t ) = Ap + w and for k = ±
1, ( f t ) k = w k .The amplitude A may be related to the order parameter close to the bifurcationby r = A ∗ + O ( | A | A ∗ ). Using Eq. (59), we obtain via the projection ( q, (45)) τ and(45) − (( q, (45)) τ p + c . c . ) the time evolution of A ( t ) and w as˙ A = λA + ( q, F [ f t ]) τ = λA + Z d ω e ψ ∗ N [ f t ] , (64)˙ w = D w + F [ f t ] − (cid:18) p Z d ω e ψ ∗ N [ f t ] + c . c . (cid:19) , (65)where c . c . stands for complex conjugate. Here, in arriving at the last equation, wehave used ( q, F [ f t ]) τ = R d ω d θ q (0) ∗ N [ f t ] = R d ω d θ e ψ ∗ ( ω ) e − iθ N [ f t ]. From (63),we get at first order N [ f t ] = − A | A | K r [ p ∗ ]( − τ ) e − i ( α − ω τ ) w , (0) + O ( A | A | ) , (66)where we have denoted the leading order of the second harmonic w = A w , + O ( A | A | ).Using the second harmonic of Eq. (65) and ˙( A ) = 2 A ˙ A = 2 A λ + O ( | A | A )gives for ϕ = 0, w , = w , (0) e ϕτ . The equation for ϕ = 0 with Eq. (63) for k = 2gives2 A λw , (0) = A (cid:16) − iωw , (0) + Kr ∗ [ p ]( − τ ) e i ( α − ω τ ) p (0) (cid:17) + O ( | A | A ) . (67)We thus get w , = (cid:18) K (cid:19) g ( ω )( λ + iω ) e i ( α − ω τ ) e − λτ e ϕτ . (68)Plugging Eqs. (68), (66), and (64) in Eq. (64), we obtain the desired normal formfor the time evolution of A ( t ):˙ A = λA + c | A | A + O ( | A | A ) , (69) where the cubic coefficient c is given by c = − K e − λτ − λ ∗ τ + i ( α − ω τ ) Λ ′ ( λ ) Z d ω g ( ω )( λ + iω ) = − K e − λ r τ R d ω g ( ω )( λ + iω ) R d ω g ( ω )( λ + iω ) + τ R d ω g ( ω )( λ + iω )= − K e − λ r τ R d ω g ( ω )( λ + iω ) R d ω g ( ω )( λ + iω ) + τ K e λτ − i ( α − ω τ ) (70) K → K + c −→ − K c (cid:18) πg ′′ ( λ i ) − i PV Z d ω g ′′ ( ω ) λ i + ω (cid:19) × " PV Z d ω g ′ ( ω ) λ i + ω − τK c cos( α − ( ω + λ i ) τ ) − i (cid:18) πg ′ ( λ i ) − τK c sin( α − ( ω + λ i ) τ ) (cid:19) − . (71)Here, we have used Eq. (56), together with the property that g ( ω ) = g ( − ω ), toobtain the last equality. The sign of Re( c ) given by the last equation gives thenature of the bifurcation as K → K + c . Note that contrary to similar unstable mani-fold analysis [38, 39, 41] c is not diverging as λ → + + λ i , which validates formallythe asymptotic analysis. Equations (57) and (71) suggest that at bifurcation, theeffects of changing τ at a fixed α are the same as those from changing τ at afixed α keeping α − ω τ constant. Interestingly, the sign of Re( c ) predicted byour analysis, which determines the nature of bifurcation, shows oscillations withdelay (see Fig. 1) that agree qualitatively with what is observed in many oscillatorsystems with delayed coupling or control [54, 55, 56].For a fixed ω and by varying τ , one may plot the sign of c by computing atcriticality K = K c ( τ ) and λ c ( τ ) = 0 + + iλ i ( τ ). The result for α = 0 is shown inFig. 1 for the case of the Lorentzian frequency distribution, Eq. (10), and also forthe case of a sum of two Lorentzians given by g ( ω ) = ∆π (cid:18) ω − ω c ) + ∆ + 1( ω + ω c ) + ∆ (cid:19) , (72)where ∆ is the width parameter of each Lorentzian and ± ω c their center frequencies( g has two separated maxima for ω c > ∆/ √ c )( τ ), a very small delay (here τ = 0 .
1) cansuppress the subcritical bifurcation present with a bimodal distribution in theabsence of delay and turn it into a supercritical bifurcation. In effect, a subcriticalbifurcation means that in the bifurcation regime, a small change in the value ofthe coupling K leads to a large change in the value of the order parameter, thatis, an incoherent state becomes with a small change of K a synchronized state,and vice versa. On the contrary, a supercritical bifurcation implies that a smallchange in K leads to only a small change in the order parameter, so that tuning ifurcations in the time-delayed Kuramoto model of coupled oscillators 17 τ K c Unimodal4 Re( c ) Unimodal τ K c Bimodal4 Re( c ) Bimodal Fig. 1
The figure shows for α = 0 the stability region of the incoherent state for the Lorentzianfrequency distribution (10) (Unimodal) with ∆ = 0 . ω = 3 and also for the case of sumof two Lorentzians (Bimodal), Eq. (72) with ∆ = 0 . , ω = 3 and ω c = 0 .
09. For
K > K c ,the incoherent state is unstable. The sign of 4 Re( c ) determines the super/sub-critical natureof the bifurcation from the incoherent state as K → K + c : a positive (respectively, a negative)sign implies a subcritical (respectively, a supercritical) bifurcation. The vertical dotted lines at τ = 1 and τ = 2 are where the bifurcation simulations were performed in Fig. 4 of Ref. [20] forthe case of the Lorentzian g ( ω ), Eq. (10), with the same parameters as us; the positive/negativesign of Re( c ) is fully consistent with the observed sub/super-critical bifurcations. It is knownthat for τ = 0 and for the distribution (72) with ∆/ √ < ω c < ∆ , one has subcriticalbifurcation [36,42]; we here see (the inset is a zoom for small delay for the bimodal case)that on introducing even a small amount of delay (here τ & . of K leads to a continuous change of the incoherent into a synchronized state, andvice versa. In the former case of a subcritical bifurcation, it is well known that anadiabatic tuning of K , in which K is tuned in time at a rate slow enough thatthe system is very close to the stationary state at every instant, and concomitantmonitoring of r leads to a hysteresis behavior of r as a function of K [50].While Fig. 1 represents the results for Lorentzians, it is evident from the natureof the expression (71) that it would be quite a daunting task to make just on thebasis of this expression general remarks on the nature of bifurcation for generalfrequency distributions, and every distribution has to be investigated on a case-by-case basis.It may be noted that for some special values of delay τ = τ n satisfying ω τ n =2 nπ for n ∈ Z , presence of delay in the cosine term of Eq. (6) has no effect. In thiscase, if the eigenvalue triggering the instability of the incoherent state is real (forthe distribution (72), it corresponds to | ω c | < ∆ ), it will be of multiplicity two,so that our derived two-dimensional unstable manifold is still valid. However, ifthere is a pair of complex eigenvalues (for the distribution (72), it corresponds to | ω c | > ∆ ), each one will have a multiplicity of two, and consequently, one shouldconsider instead a four-dimensional unstable manifold, as done for τ = τ = 0 inRef. [37]. For τ = τ n , there is a pair of complex eigenvalue of multiplicity one, sothat our two-dimensional unstable manifold expansion holds good. Kr ∞ τ = 0 τ = 0 . K c (0 . K c (0) Fig. 2
The figure shows the variation of the order parameter r as a function of the couplingconstant K , showing in particular the bifurcation behavior implied by Fig. 1 for the bimodalLorentzian distribution, Eq. (72), with ∆ = 0 . , ω = 3 , ω c = 0 .
09, and for two values of τ : τ = 0 and τ = 0 .
1. The data are obtained via numerical integration of the dynamics (4) fornumber of oscillators N = 64384 and with timestep δt = 10 − . For each value K , we run asimulation for a total time t = 2600 and compute r ∞ as the averaged of | r | ( t ) for t > K as the initial state of the runfor the next value of K . We first increase K , as in K → K + δK , with δK = 0 . . . K → K − δK . By this procedure, we clearly differentiate the subcritical bifurcation expectedfor τ = 0 (hysteresis behavior) from the supercritical one expected at τ = 0 . K c ( τ ) predicted in Fig. 1. For the Lorentzian distribution, (10), one may check that the normal form ob-tained from the OA-reduced-dynamics, Eq. (43), and the kinetic equation, Eq. (69),are the same.For generic g ( ω ), we may decompose F t = F st + f t as in Eq. (11), with Fouriercoefficients ( F t ) n / (2 π ) = g ( ω ) / (2 π )( α t ) n , where the ( α ) n ’s are the Fourier coeffi-cients on the unstable manifold. Using Eq. (7), we get ( α t ) = 1. From Eqs. (55)and (59), we get ( α t ) = A K e − λτ + i ( α − ω τ )+ λϕ λ + iω + O ( A | A | ) . (73)Notice that expression (73) explicitly satisfies the assumption made by the OAansatz that α ( ω, t ) → ω ) → −∞ . We cannot prove within the currentframework the validity of the OA assumption that | α ( ω, t ) | <
1. Similarly, Eqs.(68) and (59) give( α t ) = (cid:18) A K e − λτ + i ( α − ω τ )+ λϕ λ + iω (cid:19) + O ( A | A | ) . (74)As in Eq. (67), we may write equations for the k > w k, = A k w k, + O ( A k | A | ), obtaining w k, = w k, (0) e kϕτ , (75) k ( λ + iω ) w k, (0) = kK r ∗ [ p ]( − τ ) e i ( α − ω τ ) w k − , (0) + O ( | A | ) , (76) ifurcations in the time-delayed Kuramoto model of coupled oscillators 19 where we have used ˙( A ) k = kλA k + O ( A k | A | ), and the fact that the dominantcontribution in Eq. (63) always involves for k > f k − term and not f k +1 . Byinduction, we deduce for k ≥ α t ) k = ( α t ) k + O ( A k | A | ) , (77)and by taking the complex conjugate of the last equation, we obtain the cor-responding equation for k <
0. These equations clearly show that the unstablemanifold has exactly the same form as the OA manifold close to the bifurcation.Note that this feature is valid both in the absence and presence of delay. It isworthwhile to mention that the OA ansatz fails on adding a second harmonic(that is, with the form of interaction ∼ K sin θ + J sin(2 θ )) to Eq. (1), and even therelation (77) obtained with Eq. (76) is also not valid in this case. In fact in thiscase (without delay) the unstable manifold has been shown to be singular [38, 40].Nonetheless, the unstable manifold reduction still provided precious informationson the bifurcation. The singularities denote a profound change in the nature of theproblem. Studying how Eq. (76) is modified by the addition of a second harmonicor noise could be a starting point for investigations into generalizations of the OAansatz. In this work, we analyzed in detail the consequences of a time delay in the in-teraction in the case of the Kuramoto model of globally-coupled oscillators withdistributed natural frequencies, for generic choice of the frequency distribution g ( ω ). We derived as a function of the delay exact results for the stability boundarybetween the incoherent and the synchronized state and the nature in which thelatter bifurcates from the former at the critical point. Our results are obtained intwo independent ways: one, by considering the kinetic equation for the time evo-lution of the single-oscillator distribution, and two, by considering for the specificchoice of a Lorentzian distribution a reduced equation for the order parameter de-rived from the kinetic equation by invoking the celebrated Ott-Antonsen ansatz.In either case, the incoherent state, in which the oscillators are completely unsyn-chronized, is a stationary solution for all values of the coupling constant K betweenthe oscillators, but which is linearly stable only below a critical value K c of thecoupling. To examine how a stable synchronized state bifurcates from the incoher-ent state as the coupling crosses the value K c , we employed an unstable manifoldexpansion of perturbations about the incoherent state in the vicinity of the bifur-cation, which we applied both to the kinetic equation and to the correspondingOtt-Antonsen-reduced dynamics. We found that the nature of the bifurcation isdetermined by the sign of the coefficient of the cubic term in the equation describ-ing the amplitude dynamics of the unstable modes in the regime of weak linearinstability, namely, as K → K + c . Remarkably, we found that the amplitude equa-tion derived from the kinetic equation has the same form as that obtained fromthe OA-ansatz-reduced dynamics for the particular case of a Lorentzian g ( ω ), thusconfirming the power and general applicability of the OA ansatz. Moreover, quiteinterestingly, we found that close to the bifurcation, the unstable manifold hasthe same form as that of the OA manifold. This may have important bearings ontheir inter-relationship to be unravelled in future. As an explicit physical effect of the presence of delay, we demonstrated with our exact results that for a sum oftwo Lorentzians as a representative example of a bimodal frequency distribution,while absence of delay leads to a bifurcation of the synchronized from the incoher-ent state that is subcritical, even a small amount of delay changes completely thenature of the bifurcation and makes it supercritical. Acknowledgements
This work was initiated while DM was affiliated to Laboratory J.A.Dieudonn´e, Universit´e Cˆote d’Azur, Nice, France and was finalized with DM being affiliatedto Los Alamos National Laboratory (LANL). DM gratefully acknowledges the support of theU.S. Department of Energy through the LANL/LDRD Program and the Center for Non LinearStudies, LANL. The paper was written up during the visit of DM and SG to the InternationalCentre for Theoretical Physics – South American Institute for Fundamental Research, S˜aoPaulo, Brazil in May 2018 and during SG’s extended stay at the Universidade Federal de S˜aoCarlos and the Centro de Pesquisa em ´Optica e Fotˆonica, S`ao Carlos, Brazil during June 2018.The authors thank these institutions for warm hospitality and financial support.
References
1. Y. Kuramoto,
Chemical oscillations, waves, and turbulence (Springer-Verlag, Berlin, 1984).2. S. H. Strogatz,
From Kuramoto to Crawford: Exploring the onset of synchronization inpopulations of coupled oscillators , Physica D , 1 (2000).3. J. A. Acebr´on, L. L. Bonilla, C. J. P´erez Vicente, F. Ritort and R. Spigler,
The Kuramotomodel: A simple paradigm for synchronization phenomena , Rev. Mod. Phys. , 137 (2005).4. S. Gupta, A. Campa and S. Ruffo, Kuramoto model of synchronization: Equilibrium andnon-equilibrium aspects , J. Stat. Mech.: Theory Exp. R08001 (2014).5. F. A. Rodrigues, T. K. DM. Peron, P. Ji and J. Kurths,
The Kuramoto model in complexnetworks , Phys. Rep. , 1 (2016).6. S. Gherardini, S. Gupta and S. Ruffo,
Spontaneous synchronization and nonequilibriumstatistical mechanics of coupled phase oscillators , Contemporary Physics , 229 (2018).7. S. Gupta, A. Campa and S. Ruffo, Statistical physics of synchronization (Springer-Verlag,Berlin, 2018).8. A. Pikovsky, M. Rosenblum and J. Kurths,
Synchronization: a Universal concept in non-linear sciences (Cambridge University Press, Cambridge, 2001).9. S. H. Strogatz,
Sync: the emerging science of spontaneous order (Hyperion, New York,2003).10. J. Buck,
Synchronous rhythmic flashing of fireflies. II. , Q. Rev. Biol. , 265 (1988).11. C. S. Peskin, Mathematical aspects of heart physiology (Courant Institute of MathematicalSciences, New York, 1975).12. I. Kiss, Y. Zhai and J. Hudson,
Emerging coherence in a population of chemical oscillators ,Science , 1676 (2002).13. A. A. Temirbayev, Z. Zh. Zhanabaev, S. B. Tarasov, V. I. Ponomarenko and M. Rosen-blum,
Experiments on oscillator ensembles with global nonlinear coupling , Phys. Rev. E ,015204(R) (2012).14. S. P. Benz and C. J. Burroughs, Coherent emission from twodimensional Josephson junc-tion arrays , Appl. Phys. Lett. , 2162 (1991).15. M. Rohden, A. Sorge, M. Timme and D. Witthaut, Self-Organized synchronization indecentralized power grids , Phys. Rev. Lett. , 064101 (2012).16. L. Herrgen, S. Ares, L. G. Morelli, C. Schr¨oter, F. J¨ulicher and A. C. Oates,
Intercellularcoupling regulates the period of the segmentation clock , Curr. Biol. , 1244 (2010).17. L. Wetzel, D. J. J¨org, A. Pollakis, W. Rave, G. Fettweis and F. J¨ulicher, Self-organizedsynchronization of digital phase-locked loops with delayed coupling in theory and experiment .PLoS ONE , e0171590 (2017).18. F.-C. Blondeau and G. Chauvet, Stable, oscillatory, and chaotic regimes in the dynamicsof small neural networks with delay , Neural Netw. , 735 (1992).19. E. Niebur, H. G. Schuster and D. M. Kammen, Collective frequencies and metastabilityin networks of limit-cycle oscillators with time delay , Phys. Rev. Lett. , 2753 (1991).20. M. K. S. Yeung and S. H. Strogatz, Time delay in the Kuramoto Model of coupled oscil-lators , Phys. Rev. Lett. , 648 (1999).ifurcations in the time-delayed Kuramoto model of coupled oscillators 2121. H. Sakaguchi and Y. Kuramoto, A soluble active rotator model showing phase transitionsvia mutual entrainment , Prog. Theor. Phys. , 576 (1986).22. E. Montbri´o, D. Paz´o and J. Schmidt, Time delay in the Kuramoto model with bimodalfrequency distribution , Phys. Rev. E , 056201 (2006).23. E. Ott and T. M. Antonsen, Low dimensional behavior of large systems of globally coupledoscillators , Chaos , 037113 (2008).24. E. Ott and T. M. Antonsen, Long time evolution of phase oscillator systems , Chaos ,023117 (2009).25. N. J. Balmforth and R. Sassi, A shocking display of synchrony , Physica D, , 21 (2000).26. J. A. Carrillo, Y. P. Choi and L. Pareschi,
Structure preserving schemes for the continuumKuramoto model: phase transitions , arXiv preprint arXiv:1803.03886 (2018).27. M. Wolfrum, S. V. Gurevich and O. E. Omelchenko,
Turbulence in the Ott-Antonsenequation for arrays of coupled phase oscillators , Nonlinearity , 257 (2016).28. D. Paz´o and E. Montbri´o, From quasiperiodic partial synchronization to collective chaosin populations of inhibitory neurons with delay , Phys. Rev. Lett. , 23 (2016).29. E. A. Martens, C. Bick and M. J. Panaggio,
Chimera states in two populations withheterogeneous phase-lag , Chaos , 094819 (2016).30. C. R. Laing, Traveling waves in arrays of delay-coupled phase oscillators , Chaos ,094802 (2016).31. E. Ott and T. M. AntonsenJr., Frequency and phase synchronization in large groups:Low dimensional description of synchronized clapping, firefly flashing, and cricket chirping ,Chaos , 051101 (2017).32. D. S. Goldobin, A. V. Pimenova, M. Rosenblum and A. Pikovsky, Competing influence ofcommon noise and desynchronizing coupling on synchronization in the Kuramoto-Sakaguchiensemble , EPJ ST , 1921 (2017).33. X. Zhang, A. Pikovsky and Z. Liu,
Dynamics of oscillators globally coupled via two meanfields , Sci. Rep. , 2104 (2017).34. J. K. Hale, Linear functional differential equations with constant coefficients , Contribu-tions to Differential Equations , 291 (1963).35. J. K. Hale and S. M. V. Lunel, Introduction to functional differential equations (Springer-Verlag, New York, 1993).36. J. D. Crawford,
Amplitude expansions for instabilities in populations of globally-coupledoscillators , J. Stat. Phys. , 1047 (1994).37. J. D. Crawford, Universal trapping scaling on the unstable manifold for a collisionlesselectrostatic mode , Phys. Rev. Lett. , 656 (1994).38. J. D. Crawford, Scaling and singularities in the entrainment of globally coupled oscillators ,Phys. Rev. Lett. , 4341 (1995).39. J. D. Crawford, Amplitude equations for electrostatic waves: universal singular behaviorin the limit of weak instability , Physics of Plasmas , 97 (1995).40. J. D. Crawford and K. T. R. Davies, Synchronization of globally coupled phase oscillators:singularities and scaling for general couplings , Physica D , 1 (1999).41. J. Barr´e and D. M´etivier,
Bifurcations and singularities for coupled oscillators with inertiaand frustration , Phys. Rev. Lett. , 214102 (2016).42. E. A. Martens, E. Barreto, S. H. Strogatz, E. Ott, P. So and T. M. Antonsen
Exact resultsfor the Kuramoto model with a bimodal frequency distribution , Phys. Rev. E , 026204(2009).43. T. D. Frank, Kramers-Moyal expansion for stochastic differential equations with singleand multiple delays: Applications to financial physics and neurophysics , Phys. Lett. A ,552 (2007).44. Dietert H and Fernandez B.
The mathematics of asymptotic stability in the Kuramotomodel , Proceedings of the Royal Society A. (2018) 12;474(2220):20180467.45. J. Murdock,
Normal forms and unfoldings for local dynamical systems (Springer-Verlag,New York, 2006).46. S. Guo and J. Wu,
Bifurcation theory of functional differential equations (Springer-Verlag,New York, 2013).47. B. Niu and Y. Guo,
Bifurcation analysis on the globally coupled Kuramoto oscillators withdistributed time delays , Physica D , 23 (2014).48. H. Chiba,
A proof of the Kuramoto conjecture for a bifurcation structure of the infinite-dimensional Kuramoto model , Ergodic Theory and Dynamical Systems , 762 (2013).49. H. Dietert, Stability and bifurcation for the Kuramoto model , Journal de Math´ematiquesPures et Appliqu´ees , 451 (2016).2 David M´etivier, Shamik Gupta50. S. H. Strogatz,
Nonlinear Dynamics and Chaos: with Applications to Physics, Biology,Chemistry, and Engineering (Westview Press, Boulder, 2014).51. J. D. Crawford and P. D. Hislop,
Application of the method of spectral deformation to theVlasov-poisson system , Annals of Physics , 265 (1989).52. S. H. Strogatz and R. E. Mirollo,
Stability of incoherence in a population of coupledoscillators , J. Stat. Phys. , 613 (1991).53. S. H. Strogatz, R. E. Mirollo and P. C. Matthews, Coupled nonlinear oscillators below thesynchronization threshold: relaxation by generalized Landau damping , Phys. Rev. Lett. ,2730 (1992).54. A. Y. T. Leung, H. X. Yang and P. Zhu, Periodic bifurcation of Duffing-van der Pol os-cillators having fractional derivatives and time delay , Communications in Nonlinear Scienceand Numerical Simulation , 1142 (2014).55. X. Xu, H. Y. Hu and H. L. Wang, Stability, bifurcation and chaos of a delayed oscillatorwith negative damping and delayed feedback control , Nonlinear Dyn , 117 (2007).56. Chol-Ung Choe, Ryong-Son Kim, Hyok Jang, P. H¨ovel and E. Sch¨oll, Delayed-feedbackcontrol: arbitrary and distributed delay-time and noninvasive control of synchrony in net-works with heterogeneous delays , Int. J. Dynam. Control2