Discordant synchronization patterns on directed networks of identical phase oscillators with attractive and repulsive couplings
DDiscordant synchronization patterns on directed networks of identical phase oscillators withattractive and repulsive couplings
Thomas Peron ∗ Instituto de Ciˆencias Matem´aticas e de Computac¸ ˜ao, Universidade de S˜ao Paulo, S˜ao Carlos 13566-590, S˜ao Paulo, Brazil
We study the collective dynamics of identical phase oscillators on globally coupled networks whose interac-tions are asymmetric and mediated by positive and negative couplings. We split the set of oscillators into twointerconnected subpopulations. In this setup, oscillators belonging to the same group interact via symmetriccouplings while the interaction between subpopulations occurs in an asymmetric fashion. By employing thedimensional reduction scheme of the Ott-Antonsen (OA) theory, we verify the existence of traveling-wave and π -states, in addition to the classical fully synchronized and incoherent states. Bistability between all collectivestates is reported. Analytical results are generally in excellent agreement with simulations; for some parametersand initial conditions, however, we numerically detect chimera-like states which are not captured by the OAtheory. I. INTRODUCTION
The Kuramoto model of coupled phase oscillators has be-come over the years a paradigmatic tool for the study of emer-gent synchronization phenomena in nonlinear sciences. Inits first formulation, globally coupled oscillators interact viathe sine of the differences of their phases; this interaction isweighted by a positive coupling strength, and by increasing itsmagnitude, the phases are gradually pulled towards a commonvalue, creating then a synchronization phase transition [1–3]. Initially conceived as a solvable extension of the modelproposed by Winfree [4], the Kuramoto model attracted greatattention thanks to its analytical tractability and its later dis-covered potential to describe synchronization phenomena ina diverse set of systems, such as in optomechanical cells [5],Josephson junctions [6], chemical oscillators [7], power net-works [8], and even the synchrony among violin players [9].For a long list of examples of the use of Kuramoto models inreal applications see the reviews in Refs. [1–3].Many variations of the original model by Kuramoto havebeen inspired by particular features found in different physi-cal systems [1]. One example is the seminal work carried outby Daido [10], who, inspired by spin-glass models, treated thecouplings between oscillators in a Kuramoto model as randomvariables which could be either positive or negative. Daido’sresults provided evidence for an analogous glass phase tran-sition in oscillatory systems; however, the precise conditionsfor the existence of those “oscillator glasses” have remainedunclear, and still some debate surrounds the problem [11–15].After the early works by Daido and others [10–14, 16], thediscussion on oscillator glasses was brought back to atten-tion by Hong and Strogatz, who in a series of papers [17–19]further exploited the role of negative and positive couplings.In the first coupling setting considered by them [17, 18], os-cillators were divided into two globally coupled subpopula-tions characterized by distinct coupling strengths. In a sce-nario resembling sociodynamical models, oscillators withinthe first subpopulation were modeled to have the tendency toalign with the mean-field (conformist oscillators), whereas the ∗ [email protected] second subpopulation was defined by oscillators that are re-pelled by the other units (contrarian oscillators). In the secondmodel [19], a fraction of the oscillators was considered to pro-vide positive coupling inputs to other nodes, while the remain-ing contributed with negative couplings; that is, in mathemat-ical terms, the coupling variable was placed inside the sum-mation term of the interaction function. Despite being verysimilar, the two coupling formulations have been shown toyield significantly different collective dynamics; in fact, onlythe model in Refs. [17, 18] was found to lead to different tran-sitions other than between incoherence and classical partiallysynchronized states.Although Hong and Strogatz did not bring new evidenceto support or discard the existence of oscillator-glasses, theirpapers motivated several other studies on discordant synchro-nization patterns – i.e., states characterized by the separationof the population of oscillators into partially synchronizedclusters – induced by the coexistence of attractive and repul-sive couplings (see, e.g., [15, 20–28]). Of particular interesthere is the work by Sonnenschein et al. [24], where the au-thors unified the coupling settings of Refs. [17–19] into a sin-gle model that also included the influence of stochastic fluc-tuations on the frequencies. More specifically, in Ref. [24],Kuramoto oscillators were set to interact concomitantly via acoupling K i , which was placed outside the summation termof the interaction function, regulating the neighboring influ-ence perceived by oscillator i ; and a coupling G i placed in-side the sum, endowing the oscillators with the ability to con-tribute differently to the mean-field. By employing the di-mension reduction framework offered by the Gaussian Ap-proximation [29], the authors showed that all states previouslyreported in Refs. [17–19] (namely, traveling waves and π -states, and conventional incoherent and synchronous states)persisted in the model with both types of couplings, but undernew routes outlining the transitions between different states.The relevance of mixing positive and negative couplingsin phase-oscillator models actually goes beyond the theoret-ical interest in oscillator glasses: It turns out that certaintypes of physical, biological and chemical systems can indeedbe described as oscillators coupled through attractive and re-pulsive interactions. Noteworthy real-world examples show-ing similar characteristics and phenomena to those described a r X i v : . [ n li n . AO ] J a n above include laser arrays [30, 31] and electrochemical os-cillators [32, 33]. Furthermore, the balance between phase-attraction and -repulsion has been recently shown to be a keyfactor in the regulation of circadian rhythms by pacemakerscells in the suprachiasmatic nucleus (SCN) [34].Motivated by the aforementioned contributions, here we in-vestigate the model in Ref. [24] of identical oscillators in theabsence of stochastic fluctuations acting on the frequencies.We divide the oscillators into two subpopulations asymmet-rically coupled, and employ the theory by Ott and Antonsen(OA) [35, 36] to obtain a reduced set of equations that de-scribes the evolution of the system. By studying the linearstability of the reduced system, we analytically derive sev-eral conditions that delineate the transitions between synchro-nized, incoherent, traveling-waves and π -states. Interestingly,we find the dynamics of the present model to be overall moreintricate than its stochastic version [24], with wider regions inthe parameter space exhibiting coexistence between differentsynchronization patterns. As we shall see, simulations withlarge populations of oscillators in general confirm with excel-lent agreement the predictions by the theory; for a small setof parameters, however, we report strong deviations from thedynamics yielded by the reduced system. II. MODEL
Following Sonnenschein et al. [24], we study here the sys-tem made up of N identical Kuramoto oscillators whose equa-tions are given by ˙ θ i = ω + K i N N (cid:88) j =1 G j sin( θ j − θ i ) , (1)where i = 1 , ..., N , and ω is the natural frequency. Noticethat, in contrast to Ref. [24], we do not consider identical os-cillators under the influence of stochastic fluctuations in thephase dynamics; instead, the only source of disorder is theone inflicted by the coupling strengths. We henceforth refer toparameters K i and G j as the in- and out-coupling strengths,respectively. As defined in Eq. (1), these couplings set theinteractions between the oscillators to be asymmetric (or di-rected): oscillator i contributes to the dynamics of neighbor-ing nodes with weight G i , while the input arising from otheroscillators is weighted by K i . The first model by Hong andStrogatz [17] is recovered when G i = 1 ∀ i in Eq. (1) (no out-coupling strengths), while the second model investigated bythe same authors [19] is obtained by symmetrizing the in-coupling strengths, i.e., K i = 1 ∀ i . Bifurcation conditionshave been calculated recently for a stochastic system with acoupling scheme similar to Eq. (1) [37]. Other similar formsof the coupling setting of Eq. (1) have also been addressed re-cently in Refs. [15, 38], considering a phase frustration termin the interaction function (Kuramoto-Sakaguchi model [39]),and in populations of asymmetrically coupled R¨ossler oscilla-tors [40]. III. DIMENSIONAL REDUCTION
In the continuum limit N → ∞ , we rewrite the originalEq. (1) by omitting the sub-indexes as ˙ θ = ω + KR sin(Θ − θ ) , (2)where R and Θ are the “weighted” order parameter and themean-field phase, respectively, defined by Re i Θ = (cid:90) (cid:90) Gr K,G e iψ K,G P ( K, G ) dKdG. (3) P ( K, G ) is the joint distribution of in- and out-couplingstrengths; variables r K,G e iφ K,G are the local order parametersthat quantify the synchrony within subpopulations: z K,G = r K,G e iψ K,G = (cid:90) ρ ( θ, t | K, G ) e iθ dθ, (4)where ρ ( θ, t | K, G ) is the probability density function of ob-serving an oscillator with phase θ at time t for a givencoupling pair ( K, G ) . Henceforth we adopt the notation ρ K,G ( θ, t ) ≡ ρ ( θ, t | K, G ) . The normalization condition (cid:82) π − π ρ K,G ( θ, t ) dθ = 1 leads to the following continuity equa-tion: ∂ρ K,G ∂t + ∂∂θ { ρ K,G [ ω + KR sin(Θ − θ )] } = 0 . (5)Next we expand the phase density ρ K,G ( θ, t ) in a Fourier se-ries and apply the ansatz by Ott and Antonsen [35, 36] to itscoefficients to get ρ K,G ( θ, t ) = 12 π (cid:32) N (cid:88) n =1 [ α K,G ( ω , t )] n e inθ + c . c (cid:33) , (6)where α K,G ( t ) ≡ α ( K, G, t ) , and c.c. stands for the complexconjugate. Substituting Eq. (6) into Eq. (5) yields ˙ α K,G + iω α K,G + K α K,G R − R ∗ ) = 0 , (7)Inserting Eq. (5) into Eq. (4), we have that the local orderparameters become z K,G = α ∗ K,G . Hence, for a general dis-tribution of coupling strengths P ( K, G ) , we get ˙ r K,G = K − r K,G ) (cid:104)(cid:104) G (cid:48) r K (cid:48) ,G (cid:48) cos( ψ K,G − ψ K (cid:48) ,G (cid:48) (cid:105)(cid:105) ˙ ψ K,G = ω − K r K,G + r − K,G ) (cid:104)(cid:104) G (cid:48) r K (cid:48) ,G (cid:48) sin( ψ K (cid:48) ,G (cid:48) − ψ K,G ) (cid:105)(cid:105) (8)where (cid:104)(cid:104) ... (cid:105)(cid:105) = (cid:82) (cid:82) P ( K (cid:48) , G (cid:48) ) ...dK (cid:48) dG (cid:48) .The boundaries of the asynchronous state ( R = 0 )can be obtained straightforwardly for arbitrary distributions P ( K, G ) . By considering small perturbations δr K,G aroundthe incoherent state r K,G = 0 , and setting ψ K,G = 0 , withoutloss of generality, we get from Eq. (8) ˙ δr K,G = K (cid:104)(cid:104) G (cid:48) δr K (cid:48) ,G (cid:48) (cid:105)(cid:105) . (9) FIG. 1. Schematic illustration of a finite system composed of twointertwined subpopulations. Oscillators belonging to subpopula-tion 1 (2) interact among themselves via couplings K G ( K G ) ,while oscillators from different subpopulations interact asymmetri-cally with effective couplings K G and K G , as depicted above. By multiplying the previous equation by G and averaging overthe distribution P ( K, G ) , we rewrite Eq. (9) in terms of a per-turbation to the global order parameter δR = (cid:104)(cid:104) G (cid:48) δr K (cid:48) ,G (cid:48) (cid:105)(cid:105) , ˙ δR = [ (cid:104)(cid:104) KG (cid:105)(cid:105) / δR , which leads to the critical condition (cid:104)(cid:104) KG (cid:105)(cid:105) = 0 . (10)Therefore, for (cid:104)(cid:104) KG (cid:105)(cid:105) < , the oscillators remain incoherent,while for (cid:104)(cid:104) KG (cid:105)(cid:105) > the incoherent state loses stability, anda partially synchronized state sets in. Notice that Eq. (10)is similar to the condition obtained in Ref. [24] for identicaloscillators subjected to Gaussian noise.Let us consider now the case in which the oscilla-tors are coarse-grained into n intertwined subpopulationswith joint distribution of couplings given by P ( K, G ) = n (cid:80) nq =1 δ [( K, G ) − ( K q , G q )] . Substituting the previous ex-pression for P ( K, G ) into Eqs. (8) yields ˙ r q = − K q n (1 − r q ) n (cid:88) p =1 G p r p cos( ψ p − ψ q )˙ ψ q = ω − K q n ( r q + r − q ) n (cid:88) p =1 G p r p sin( ψ p − ψ q ) . (11)In what follows we investigate a special case of the abovesystem, namely, the setup of two intertwined subpopula-tions [24]. In this case, we have n = 2 , and Eqs. (11) arereduced to ˙ r = K − r )[ r G + r G cos δ ] , ˙ r = K − r )[ r G + r G cos δ ] , ˙ δ = − sin δ (cid:2) ( r + r − ) r K G + ( r + r − ) r K G (cid:3) , (12)where we have defined the phase-lag δ = ψ − ψ . An illustra-tion of a finite network with two intertwined subpopulationscan be seen in Fig. 1. We measure the global synchronization with the classical Kuramoto order parameter as: r ( t ) e i Φ( t ) = 12 [ r ( t ) e iψ ( t ) + r ( t ) e iψ ( t ) ] (13)Note that r ( t ) e i Φ( t ) in the above equation is different fromthe “weighted” order parameter Re i Θ = [ r ( t ) G e iψ ( t ) + r ( t ) G e iψ ( t ) ] [Eq. (3)], which can be larger than one.Equations (12) are very similar to the set of equations foridentical oscillators under the influence stochastic fluctua-tions obtained via Gaussian approximation [24, 29]. Actually,the only difference between the reduced system obtained inRef. [24] and Eq. (12) is that the former exhibits terms r , instead of r , in the equations for ˙ r , ; and terms propor-tional to ( r − , + r , ) in place of ( r − , + r , ) in the equationfor ˙ δ . Notice also that Eqs. (12) could be obtained via theWatanabe-Strogatz theory [41, 42] under uniform distributionof constants of motion (Ott-Antonsen manifold) [18, 36].From Eqs. (12), we expect to observe the following sta-tionary states for the two subpopulation system (see the illus-tration in Fig. 2): (i) the classical incoherent state in which r = r , = 0 ; (ii) the perfectly synchronized state in which r , = 1 and δ = 0 (we denominate this state as “zero-lagsync” state); (iii) partially synchronized states characterizedby r , < and δ = 0 , and which we refer to as “blurredzero-lag sync” states; (iv) the so-called “ π -state” for which thesubpopulations are perfectly synchronized ( r , = 1 ), whileremaining diametrically opposed in the phase space ( δ = π ),yielding, hence, a vanishing global synchronization ( r = 0 );(v) “blurred“ π -states, in which at least one of the subpopu-lations is partially synchronized ( r , < ) and the peaks oftheir phase distributions are separated by δ = π ; and, finally,(vi) the traveling-wave (TW) state [17, 24, 43] in which thesubpopulations can be either partially or fully synchronized, < r , ≤ , while keeping a constant phase-lag separa-tion within < δ < π . The interesting feature of this stateis that, in contrast to standard formulations of the Kuramotomodel, the oscillators no longer rotate with a common fre-quency given by the frequency ω of the co-rotating frame–or ¯ ω = (cid:82) ωg ( ω ) dω in the case of non-identical oscillators,where g ( ω ) is a frequency distribution [17, 43]–; instead, theysettle on a stationary new rhythm whose magnitude will alsodepend on the coupling parameters. Deviations from the meanfrequency ¯ ω can be calculated either by the average (or mean-ensemble) frequency Ω = 1 N N (cid:88) j =1 (cid:104) ˙ θ j (cid:105) t , (14)where (cid:104) ... (cid:105) t denotes a long-time average; or by the lockingfrequency ˙Θ , defined in Eq. (3) [see also Eq. (18)]. TWstates typically appear in phase-oscillator systems when cer-tain symmetry patterns are broken in the model, such as by thepresence of a phase frustration in the sine coupling term [39],asymmetric coupling strength distributions [17, 18], or by nat-ural frequencies asymmetrically distributed [43]. FIG. 2. Long-time profile on the complex unit circle of the states observed for a finite system [Eq. (1)] with two intertwined subpopulationscoupled as depicted in Fig. 1: (a) Incoherent state, (b) Zero-lag sync, (c) π -state, (d) TW1, (e) blurred zero-lag sync, and (f) blurred π -state.Only in the TW1 state [panel (d)], the collective frequencies ˙Θ [Eq. (18)], ˙Φ [Eq. (13)] and Ω [Eq. (14)] are different from zero, and theoscillators travel across all possible phase values in the co-rotating frame defined by the natural frequency. The configuration of the TW2 stateis obtained by interchanging sub-indexes and colors in panel (d). IV. BIFURCATION ANALYSIS OF THE REDUCEDSYSTEM WITH TWO SUBPOPULATIONS
For convenience, we adopt the following parametrizationfor the couplings: K , = K ± ∆ K and G , = G ± ∆ G , (15)where K and G are the average in- and out-couplingstrengths, respectively; parameters ∆ K and ∆ G are definedas the corresponding coupling mismatches. In our calcula-tions we always consider positive mismatches ( ∆ K, ∆ G > ). Therefore, if | K | < ∆ K/ or | G | < ∆ G/ , half of thecouplings are positive (attractive) and half are negative (re-pulsive). If one of these conditions is satisfied, we say that theoscillators interact via mixed couplings .By setting ˙ δ = 0 , we uncover two possible fixed-point so-lutions for phase-lag δ : δ = mπ , m ∈ Z (16) r + r − ) r K G + ( r + r − ) r K G . (17) Equation (16) corresponds to the solution of partially synchro-nized states with no separation between populations (even m )and π -states (odd m ), while Eq. (17) gives the condition forthe existence of TW states. We can verify that TWs appearonly for < δ < π by rewriting the equations for ˙ ψ , to-gether with Eq. (17) as: lim t →∞ ˙ ψ , = lim t →∞ ˙Θ = ω − sin δ r + r − K G r . (18)Therefore, spontaneous drifts in the collective frequencies oc-cur only for intermediate values of the phase-lag δ ; otherwise,for δ = mπ , oscillators rotate with collective frequencies ˙ ψ , = ω , which here is set to ω = 0 .From the parametrization in Eq. (15), we realize the criticalconditions K , = 0 , or K = ± ∆ K . (19)When one of the above conditions holds, it follows thatone subpopulation is deprived of receiving inputs from othernodes (including from the same subpopulation), and its oscil-lators have instead only out-couplings towards nodes externalto their subpopulation. Similarly, if one of the out-couplingsvanishes, G , = 0 , the corresponding subpopulation ceasesto influence the dynamics of the rest of the network and startsacting only as link receiver. As we shall see, these conditionsplay an important role in the phase diagram of Eq. (12). In thesequel, we calculate the coupling ranges in which the statesdiscussed in the previous section appear. A. Incoherent state
The first critical condition of Eqs. (12) is given by the stabil-ity analysis of the incoherent state performed in the last sec-tion. For the case of two subpopulations with P ( K, G ) = δ [( K, G ) − ( K , G )] + δ [( K, G ) − ( K , G )] , Eq. (10)reads K G + K G = 0 , and by solving it in terms of theaverage in-coupling strength we have K = − ∆ K ∆ G G . (20)The above equation, therefore, delineates the boundary of theincoherent state. B. π -states For the π -state, the fixed point solutions read r , = 1 and δ = π . Linear stability reveals that the π -state is stable foraverage in-couplings given by − ∆ K < K < min (cid:26) ∆ K , ∆ K ∆ G G (cid:27) , for G > (cid:26) − ∆ K , ∆ K ∆ G G (cid:27) < K < ∆ K , for G < . (21)The other state characterized by subpopulations diametri-cally opposed is the blurred π -state, and its fixed points aredefined by r G = r G and δ = π . Linearizing the dy-namics about this state, we find that the corresponding Jaco-bian matrix has a single nonzero eigenvalue, λ = K G (1 − r ) / K G (1 − r ) / . Hence, and because couplings G and G must have the same sign so that r G = r G is aphysical solution, we have that blurred π -states appear when K < ∆ K ∆ G G , for | G | > ∆ G . (22)Thus, both π -states with fully synchronous subpopulations( r , = 1 ) and with partial synchronization ( r , < ) areyielded by systems described by Eqs. (12). Observe also that r G = r G defines a one-parameter family of fixed points.Furthermore, from Eq. (22) we notice that incoherent andblurred π -states may coexist in a large region of the param-eter space defined by coupling strengths G and K . C. Zero-lag sync and partially synchronized states
For the zero-lag sync state, we linearize Eqs. (12) around r , = 1 and δ = 0 , and seek a zero eigenvalue of the re-lated Jacobian matrix. By following this calculation, we findthat the zero-lag sync state emerges for average in-couplingstrengths given by K > max (cid:26) ∆ K , ∆ K ∆ G G (cid:27) , for G > K < min (cid:26) − ∆ K , ∆ K ∆ G G (cid:27) , for G < . (23)Next, by linearizing Eqs. (12) around an arbitrary partiallysynchronized solution ( r G = − r G and δ = 0 ), we findthat such states, which we have denominated blurred zero-lagsync states, must occur for parameters in the range K < ∆ K ∆ G G , for | G | < ∆ G . (24)Any state that satisfies r G = − r G and δ = 0 is a fixedpoint of Eqs. (12). From this we foresee that several par-tially synchronized states should coexist in the region delim-ited by Eq. (24). As in the case of blurred π -states, blurredzero-lag sync solutions yield a single nonzero Jacobian eigen-value, λ = K G (1 − r ) / K G (1 − r ) / ; hence,the latter states do not coexist with zero-lag sync nor with π -states, because in regions where r , = 1 we have λ = 0 ,and the blurred zero-lag sync states lose stability. Notice fur-ther that in order to r G = − r G be a physical solution,out-couplings G and G must have opposite signs; thus, par-tially synchronized states with δ = 0 are expected to appearonly in the presence of mixed-out coupling strengths, i.e., for | G | < ∆ G/ , as indicated in Eq. (24). D. Traveling waves
We now turn our attention to the stationary TW states. Nu-merical results show us that two possible TW states are mani-fested by the system (12). In the first state, which we refer toas “TW1”, the first subpopulation remains fully synchronized( r = 1 ), whereas the second one exhibits partial synchroniza-tion ( r < ). We label as “TW2” the opposite situation, i.e.,when r = 1 and r < . In both states we have < δ < π and Ω > . By setting r = 1 in Eqs. (12), we find the fol-lowing fixed point solutions for r and δ : r = (cid:114) − K G K G + 2 K G and cos δ = − G G r . (25)Since < r < , we have that the solution of the TW1 existsfor − K G < K G < . (26)By writing Eq. (26) in terms of the parametrization in Eq. (15)and considering ∆ K, ∆ G > , we find that the regions with (a) (b)(c) (d) FIG. 3. Bifurcation diagram of the reduced system [Eq. (12)] for (a,b) ∆ K = 8 and ∆ G = 2 ; and (c,d) K = 3 and G = 2 . Zero-lag synccorresponds to the perfectly synchronized state ( r , = 1 and δ = 0 ). “ π -state” refers to the state in which r , = 1 and phase-lag separation δ = π . Similarly, blurred zero-lag and blurred π -states denote the states with r , < along with δ = 0 and δ = π , respectively. Travelingwave states are characterized by r = 1 , r < (TW1), and r = 1 , r < (TW2). Both TW1 and TW2 exhibit Ω (cid:54) = 0 [Eq. (14)]. Solidlines are obtained from Eqs. (20)-(24) and Eqs. (27)-(29). Dashed line in panels (a) and (c) depict the condition K ∆ G − ∆ KG = 0 forwhich the couplings are symmetric, i.e., when the network connections are undirected. Panels (b) and (d) show zoomed-in regions of panels(a) and (c), respectively. TW1 are outlined by the following conditions: ∆ K ∆ G G
Let us now compare the results of the bifurcation analysisin the previous section with numerical simulations of the fi-nite original dynamics [Eq. (1)]. Figure 4 shows the evolutionof order parameters, phase-lag separation δ , and collective fre-quencies, Ω and ˙Θ , as a function of K for different choices of G . All the simulations are performed by integrating Eqs. (1)with the Heun’s method with a time step dt = 0 . and con-sidering total number of oscillators N = 10 (see the captionof Fig. 4 for more details). For G = 0 [Fig. 4 (a), (d), and(g)] we observe that the system transitions from TW2 to π -state, and then subsequently to TW1, as correctly predictedby the critical conditions depicted in the diagram of Fig. 3. InFig. 4(b), we see that at K = − the subpopulations abruptlysynchronize as they switch from blurred π -states to π -state. Asimilar discontinuous transition of the local order parameters r , was observed for similar parameter configurations in thestochastic version of the system (1) [24]. Abrupt transitionsare also seen in Fig. 4(c), but this time as a consequence ofthe transition from blurred zero-lag sync to TW2 state. InFig. 4(c) we also observe irregular points in the “B. zero-lagsync” region. Those points correspond to partially synchro-nized states and display such an irregular pattern because ofthe one-parameter family of solutions that exists in that re-gion; specifically, different initial conditions drive the systemto different stationary states that satisfy r G = − r G . Thesolid branches in “B. zero-lag sync” area correspond to TW2solutions, which are also stable for K (cid:46) − . in Fig. 4(c)[see also Fig. 3(b)], but are not obtained numerically with theinitial conditions used in Fig. 4. We shall return to this pointshortly.Notice in Fig. 4(g)-(i) that | Ω | ≤ | ˙Θ | . The reason for thisresides in the fact that Ω is a microscopic average of the in-stantaneous frequencies (cid:104) ˙ θ i (cid:105) , while ˙Θ [and equivalently ˙Φ inEq. (13)] quantifies how fast the center of the bulk formed by TW2 -state TW1 -stateBlurred -state TW1
Zero-lagsync
TW2 -state TW1
Zero-lagsyncB. zero-lagsync (a) (b) (c)(d) (e) (f)(g) (h) (i)
FIG. 4. Order parameters r , , r , phase-lag δ , and collective frequencies ˙Θ [Eq. (18)] and Ω [Eq. (14)] for (a,d,h) G = 0 , ∆ K = 8 , and ∆ G = 2 ; (b,e,h) G = 2 , ∆ K = 8 , and ∆ G = 2 ; and (c,f,i) G = 2 , ∆ K = 3 , and ∆ K = 10 and ∆ G = 2 . Dots are obtained bynumerically integrating the original system [Eq. (1)] using the Heun’s method with N = 10 oscillators. For each coupling value K , thequantities are averaged over t ∈ [500 , with a time-step dt = 0 . . In all panels, initial conditions θ i ( t = 0) ∀ i are randomly distributedaccording to a uniform distribution between [ − π, π ] . Solid lines correspond to the analytic solutions obtained in Sec. IV. entrained oscillators rotates. Therefore, oscillators that are notlocked with the mean-field contribute to the sum in Eq. (14)with (cid:104) ˙ θ i (cid:105) t ≈ , thus reducing the value of Ω in comparisonwith its upper bound ˙Θ . The latter frequency offers in thepresent case the slight advantage of being calculated directlyfrom the solutions in Eq. (18). Analogously, Ω can be es-timated analytically (not shown here) through the ensembleaverage Ω = (cid:82) π − π (cid:82) (cid:82) ˙ θρ ( θ, t | K, G ) dKdGdθ .To conclude this section, in Fig. 5 we compare the theoret-ical results with simulations considering coupling parametersover a G × K grid. Our goal with this approach is to inspectfor a larger set of parameters whether the stability analysisperformed in the previous section correctly predicts the stabil-ity regions shown in Fig. 3. For each coupling pair ( G , K ) ,Eqs. (1) are integrated numerically, and the global variablesare averaged over t ∈ [500 , with a time step dt = 0 . .As can be seen in Fig. 5, the boundaries of the collective statesare predicted very accurately by the theory. Seeking to verifythe bistable behavior of the model, in Fig. 6 we show sim-ulations results for a zoomed region of the space in Fig. 5considering different initial conditions: in Fig. 6(a), phases θ i are initiated with values distributed uniformly at random be-tween [ − π, π ] . Initial conditions were chosen differently forFig. 6(b); specifically, for each point in the grid G × K , thephases of populations 1 and 2 were chosen from Gaussian dis-tributions with standard deviation σ = 2 , and means θ and θ , which were taken uniformly at random between [ − π, π ] .By initiating the oscillators in this way, we observe in Fig. 6that the system converges to TW2 in the region where thisstate was predicted to coexist with partial synchronization. VI. ACCURACY OF THE OTT-ANTONSEN REDUCTION
Although we have observed a good agreement betweensimulations and the theory for the states previously discussed,there are also dynamical patterns which seem not to be cap-tured by the OA reduction. Figure 6 shows an example: therewe observe a set of points with ˙Θ = 0 in the TW2 area, i.e., aregion where one would expect stationary states with ˙Θ (cid:54) = 0 .By inspecting the temporal trajectories of the local order pa-rameters r , in Fig. 7, we see that such states do have a dif-ferent nature from traveling waves. The trajectories in Figs. 7actually resemble breathing chimera states [44] in which onesubpopulation remains fully locked, while the other exhibitsan oscillating synchrony. In the figure, we compare the timeevolution of the original model (1) with the numerical inte-gration of the reduced system [Eq. (12)] using the same ini-tial conditions. As it is seen, while the finite subpopulation1 shows oscillating synchrony, the solution provided by thetheory converges to a constant value. Although chimera stateshave been studied extensively with the OA reduction, Figs. 6and 7 suggest that such solutions in the present model mightlie outside the OA manifold. Interestingly, Refs. [17, 18, 24]did not report solutions akin to the ones shown in Fig. 7. VII. CONCLUSION
In this paper we have studied a variant of the Kuramotomodel in which identical oscillators are coupled via in-and out-coupling strengths, which in turn can have positive
FIG. 5. Comparison between simulations (colormaps) and theory (solid lines). (a) Total order parameter r ; (b) local order parameter r ofsubpopulation 1; (c) phase-lag δ measuring the separation between the two subpopulations; and (d) average frequency Ω [Eq. (14)]. For eachpair of coupling ( G , K ) , Eqs. (1) are evolved numerically with the Heun’s method considering N = 10 oscillators and with integrationtime step dt = 0 . . The long-time behavior of each parameter is quantified by averaging the trajectories over t ∈ [500 , . For all ( G , K ) , the initial phases θ i ( t = 0) are drawn uniformly at random over the interval [ − π, π ] . Coupling mismatch parameters: ∆ K = 8 and ∆ G = 2 . and negative values. Similarly to the setting considered inRef. [24], heterogeneity in the interactions was introducedby dividing the oscillators into two mutually coupled sub-populations (each one characterized by a distinct pair of cou-plings), so that connections within the same subpopulation re-main symmetric, while connections between subpopulationsare asymmetric. In the infinite size limit, we applied thetheory by Ott and Antonsen [35] to obtain a reduced set ofequations. With the reduced description of the original sys-tem, we performed a thorough bifurcation analysis wherebya rich dynamical behavior was revealed. We showed thatthe present system exhibits different types of π -states andtraveling-waves, along with classical incoherent and partiallysynchronized states. Though the transitions among thesestates bear some similarity to those uncovered for the modelwith Gaussian white noise in Ref. [24], we have found thatour model exhibits a more intricate long-term dynamics thanthat of observed for its noisy version. The reason for this con-clusion resides in the different types of π - and zero-lag sync states uncovered (which may consist of either perfectly or par-tially synchronized subpopulations), and in the observationof wider regions in the parameter space displaying bistabil-ity. These findings for the seemingly simpler system are inline with previous studies [17, 18] comparing the dynamicsof identical oscillators with that of non-identical ones. [Al-though the system in Eq. (1) and the one of Ref. [24] areboth models of identical oscillators, the inclusion of Gaussianwhite noise yields equivalent phenomenology–with respect tothe linearized dynamics–to the case of phase oscillators withnatural frequencies drawn from Lorentzian distributions; seethe discussion in Ref. [45].]Despite the excellent agreement between simulations andtheory, for a small set of parameter combinations we veri-fied dynamical states which turned out not to be reproducedby the reduced system. As discussed in Sec. VI, we veri-fied vanishing values for the temporal average of the mean-field frequency for couplings inscribed in a TW region. Byvisualizing the time-series of such unanticipated states, we0 FIG. 6. Comparison between simulations (colormaps) and theory(solid lines) for different sets of initial conditions: (a) θ i ( t = 0) randomly chosen from the uniform distribution [ − π, π ] (same as inFig. 5); (b) for each coupling pair ( G , K ) the initial phases ofsubpopulation 1 and 2 were drawn from Gaussian distributions withstandard deviation σ = 2 , and means θ and θ , respectively, whichwere chosen uniformly at random between [ − π, π ] . The “+” marksa point in the diagram for which the behavior observed in the simu-lation departs from the dynamics predicted by the theory. Figure 7shows the temporal evolution of the collective variables at the “+”point depicted in panel (b). Other parameters: N = 10 , ∆ K = 8 and ∆ G = 2 . In both panels the resolution of the grid is × couplings. Integration was performed with the Heun’s method usingtime step dt = 0 . . found that one local parameter evolved with an oscillatory dy-namics akin to breathing chimera states [44], in sharp con-trast to the evolution predicted by our calculations for thesame parameters and initial conditions. It is worth noting,nonetheless, that disagreements of this nature are somewhatexpected to occur: for the identical frequencies case there ex-ists a one-parameter family of invariant manifolds (of whichthe OA-manifold is a special solution) that are neutrally stablewith respect to perturbations in directions transverse to them-selves [36, 42, 46, 47]. Hence, there could be certain types of (a)(b)(c) FIG. 7. Temporal evolution of (a) local order parameters r , , (b)phase-lag δ , and (c) locking frequency ˙Θ [Eq. (18)]. Solid linesare obtained from simulations, while dashed lines correspond to theresults yielded by the numerical integration of the reduced system[Eq. (12)]. In panel (a) the solid and dashed lines of r overlap eachother at r = 1 . Average in- and out-coupling strengths are takenfrom the “+” point in Fig. 6(b), that is, ( G , K ) = ( − . , − . .Other parameters: N = 10 oscillators, ∆ K = 8 , ∆ G = 2 , and dt = 0 . . perturbations that may drive the system away from the man-ifold contemplated by the OA ansatz, thus generating unex-pected results such as the ones discussed in Sec. VI. Anotherdeviation from the theory was observed in the appearance ofzero-lag sync and π -states over a region in the parameter spaceinitially believed to manifest traveling-waves solely [see Fig. 4(c)]. Future works should further investigate the emergence ofchimeras and other states with respect to perturbations off theOA manifold for populations of identical oscillators coupledasymmetrically.As mentioned in the introduction, there are many systemswhose dynamics can be modeled by phase oscillators inter-acting via positive and negative couplings. Our results canthereby serve as a guide in the search for clustered statesin different contexts. For instance, as shown in Ref. [34],phase-attraction or phase-repulsion alone cannot account for1 FIG. 8. Comparison between simulations (colormaps) and theory (solid lines). (a) Total order parameter r ; (b) local order parameter r ofsubpopulation 1; (c) phase-lag δ measuring the separation between the two subpopulations; and (d) average frequency Ω [Eq. (14)]. For eachpair of coupling ( G , K ) , Eqs. (1) are integrated numerically with the Heun’s method considering N = 10 oscillators and with integrationtime step dt = 0 . . The long-time behavior of each parameter is quantified by averaging the trajectories over t ∈ [500 , . For allcoupling pairs ( G , K ) , Eqs. (1) were initiated with the exact same configuration for θ i ( t = 0) : phases in subpopulation 1 were randomlychosen from a Gaussian distribution with mean θ (cid:39) . and standard deviation σ (cid:39) . , yielding r ( t = 0) (cid:39) . ; phases ofsubpopulation 2 were also Gaussian distributed, but with mean θ (cid:39) . and standard deviation σ (cid:39) . , yielding r ( t = 0) (cid:39) . .Coupling mismatch parameters: ∆ K = 8 and ∆ G = 2 . G × K grid resolution: × couplings. the regulation of circadian rhythms; a phase model incorpo-rating mixed couplings linked asymmetrically, on the otherhand, does reproduce the outcome of experiments with neu-ronal networks of the suprachiasmatic nucleus (SCN). There-fore, system (1) with two intertwined subpopulations may bea suitable model to describe the synchronization between thedorsal and ventral subregions of the SCN [34].Finally, there are a number of potentially relevant exten-sions for the present model: given the non-trivial behavioruncovered here, it would be interesting, for instance, to inves-tigate more than two coupled populations, as well as to studythe effect of attractive and repulsive interactions on the col-lective dynamics of oscillators with higher-order harmonicsin the coupling function [48, 49]. One could also consider os-cillators coupled on structures that allow interactions beyondthe classical pairwise, such as hypergraphs [50] and simpli-cial complexes [51]. On the experimental domain, realizations of the dynamical transitions reported here may be obtainedin populations of chemical [52, 53] and optical arrays [54],which are experimental setups that have been shown to re-produce chimeras and other dynamical states found in phaseoscillator models. ACKNOWLEDGMENTS
The author thanks Bernard Sonnenschein, Chen ChrisGong, Deniz Eroglu, Paul Schultz, Vinicius Sciuti, and BrunoMessias for useful conversations. This research was fundedby FAPESP (Grant No. 2016/23827-6) and carried out usingthe computational resources of the Center for MathematicalSciences Applied to Industry (CeMEAI) funded by FAPESP(Grant No. 2013/07375-0).2
Appendix A: Supplementary diagram
In this section we recalculate the diagrams of Fig. 5 by ini-tiating the oscillators differently than in Sec. V. Specifically,in Fig. 8 we choose the phases of each subpopulation to bedistributed according to distinct Gaussian distributions whosepeaks are separated by a phase-lag (see the caption of Fig. 8for details). Comparing Fig. 5 with Fig. 8 we see that a blurred π -state region in the former is converted into a π -state area inthe latter [notice the regions with r = 0 in Fig. 5(b) and r = 1 in Fig. 8(b)]. In addition to the bistability betweenTWs and blurred zero-lag sync states confirmed in Sec. VI,another significant difference between Fig. 5 and Fig. 8 liesin the “Incoherence/Blurred π -state” areas: in Fig. 5 these re-gions exhibit small values for the order parameters (a conse-quence of choosing the initial conditions uniformly at randombetween [ − π, π ] ), along with δ = π ; in Fig. 8, on the otherhand, we observe higher values for r , thus confirming thatmultiple local synchronization levels are possible in those re-gions as revealed by the analysis in Sec. IV. [1] Juan A Acebr´on, Luis L Bonilla, Conrad J P´erez Vicente, F´elixRitort, and Renato Spigler, “The Kuramoto model: A simpleparadigm for synchronization phenomena,” Reviews of ModernPhysics , 137 (2005).[2] Francisco A Rodrigues, Thomas K DM Peron, Peng Ji, andJ¨urgen Kurths, “The Kuramoto model in complex networks,”Physics Reports , 1–98 (2016).[3] Alex Arenas, Albert D´ıaz-Guilera, Jurgen Kurths, YamirMoreno, and Changsong Zhou, “Synchronization in complexnetworks,” Physics Reports , 93–153 (2008).[4] Arthur T Winfree, “Biological rhythms and the behavior of pop-ulations of coupled oscillators,” Journal of Theoretical Biology , 15–42 (1967).[5] Georg Heinrich, Max Ludwig, Jiang Qian, Bj¨orn Kubala, andFlorian Marquardt, “Collective dynamics in optomechanical ar-rays,” Physical Review Letters , 043603 (2011).[6] Kurt Wiesenfeld, Pere Colet, and Steven H Strogatz, “Syn-chronization transitions in a disordered Josephson series array,”Physical Review Letters , 404 (1996).[7] Istv´an Z Kiss, Yumei Zhai, and John L Hudson, “Emergingcoherence in a population of chemical oscillators,” Science ,1676–1678 (2002).[8] Florian D¨orfler and Francesco Bullo, “Synchronization in com-plex networks of phase oscillators: A survey,” Automatica ,1539–1564 (2014).[9] Shir Shahal, Ateret Wurzberg, Inbar Sibony, Hamootal Duadi,Elad Shniderman, Daniel Weymouth, Nir Davidson, and MotiFridman, “Synchronization of complex human networks,” Na-ture Communications , 1–10 (2020).[10] Hiroaki Daido, “Quasientrainment and slow relaxation in a pop-ulation of oscillators with random and frustrated interactions,”Physical Review Letters , 1073 (1992).[11] LL Bonilla, CJ P´erez Vicente, and JM Rubi, “Glassy synchro-nization in a population of coupled oscillators,” Journal of Sta-tistical Physics , 921–937 (1993).[12] JC Stiller and G Radons, “Dynamics of nonlinear oscillatorswith random interactions,” Physical Review E , 1789 (1998).[13] Hiroaki Daido, “Algebraic relaxation of an order parameter inrandomly coupled limit-cycle oscillators,” Physical Review E , 2145 (2000).[14] JC Stiller and G Radons, “Self-averaging of an order parameterin randomly coupled limit-cycle oscillators,” Physical ReviewE , 2148 (2000).[15] Dima Iatsenko, Peter VE McClintock, and Aneta Stefanovska,“Glassy states and super-relaxation in populations of coupledphase oscillators,” Nature communications , 4118 (2014).[16] Dami´an H Zanette, “Synchronization and frustration in oscil-lator networks with attractive and repulsive interactions,” EPL (Europhysics Letters) , 190 (2005).[17] Hyunsuk Hong and Steven H Strogatz, “Kuramoto model ofcoupled oscillators with positive and negative coupling param-eters: an example of conformist and contrarian oscillators,”Physical Review Letters , 054102 (2011).[18] Hyunsuk Hong and Steven H Strogatz, “Conformists and con-trarians in a Kuramoto model with identical natural frequen-cies,” Physical Review E , 046202 (2011).[19] Hyunsuk Hong and Steven H Strogatz, “Mean-field behavior incoupled oscillators with attractive and repulsive interactions,”Physical Review E , 056210 (2012).[20] Ernest Montbri´o and Diego Paz´o, “Collective synchronizationin the presence of reactive coupling and shear diversity,” Phys-ical Review E , 046206 (2011).[21] Dmytro Iatsenko, Spase Petkoski, PVE McClintock, and A Ste-fanovska, “Stationary and traveling wave states of the Ku-ramoto model with an arbitrary distribution of frequenciesand coupling strengths,” Physical Review Letters , 064101(2013).[22] Hyunsuk Hong, “Periodic synchronization and chimera in con-formist and contrarian oscillators,” Physical Review E ,062924 (2014).[23] Isabel M Kloumann, Ian M Lizarraga, and Steven H Strogatz,“Phase diagram for the Kuramoto model with van hemmen in-teractions,” Physical Review E , 012904 (2014).[24] Bernard Sonnenschein, Thomas K DM Peron, Francisco A Ro-drigues, J¨urgen Kurths, and Lutz Schimansky-Geier, “Col-lective dynamics in two populations of noisy oscillators withasymmetric interactions,” Physical Review E , 062910(2015).[25] Bertrand Ottino-L¨offler and Steven H Strogatz, “Volcano tran-sition in a solvable model of frustrated oscillators,” PhysicalReview Letters , 264102 (2018).[26] Jinha Park and B Kahng, “Metastable state en route to traveling-wave synchronization state,” Physical Review E , 020203(2018).[27] Dustin Anderson, Ari Tenzer, Gilad Barlev, Michelle Girvan,Thomas M Antonsen, and Edward Ott, “Multiscale dynam-ics in communities of phase oscillators,” Chaos: An Interdisci-plinary Journal of Nonlinear Science , 013102 (2012).[28] Jinha Park and B Kahng, “Competing synchronization on ran-dom networks,” Journal of Statistical Mechanics: Theory andExperiment , 073407 (2020).[29] Bernard Sonnenschein and Lutz Schimansky-Geier, “Approxi-mate solution to the stochastic Kuramoto model,” Physical Re-view E , 052111 (2013).[30] Chene Tradonsky, Micha Nixon, Eitan Ronen, Vishwa Pal, Ro-nen Chriki, Asher A Friesem, and Nir Davidson, “Conversion of out-of-phase to in-phase order in coupled laser arrays withsecond harmonics,” Photonics Research , 77–81 (2015).[31] Vishwa Pal, Simon Mahler, Chene Tradonsky, Asher AFriesem, and Nir Davidson, “Rapid fair sampling of the x yspin hamiltonian with a laser simulator,” Physical Review Re-search , 033008 (2020).[32] Hiroshi Kori, Istv´an Z Kiss, Swati Jain, and John L Hudson,“Partial synchronization of relaxation oscillators with repulsivecoupling in autocatalytic integrate-and-fire model and electro-chemical experiments,” Chaos: An Interdisciplinary Journal ofNonlinear Science , 045111 (2018).[33] Michael Sebek and Istv´an Z Kiss, “Plasticity facilitates patternselection of networks of chemical oscillations,” Chaos: An In-terdisciplinary Journal of Nonlinear Science , 083117 (2019).[34] Jihwan Myung, Sungho Hong, Daniel DeWoskin, ErikDe Schutter, Daniel B Forger, and Toru Takumi, “Gaba-mediated repulsive coupling between circadian clock neuronsin the scn encodes seasonal time,” Proceedings of the NationalAcademy of Sciences , E3920–E3929 (2015).[35] Edward Ott and Thomas M Antonsen, “Low dimensional be-havior of large systems of globally coupled oscillators,” Chaos:An Interdisciplinary Journal of Nonlinear Science , 037113(2008).[36] Arkady Pikovsky and Michael Rosenblum, “Dynamics of glob-ally coupled oscillators: Progress and perspectives,” Chaos:An Interdisciplinary Journal of Nonlinear Science , 097616(2015).[37] JM Meylahn, “Two-community noisy Kuramoto model,” Non-linearity , 1847 (2020).[38] Vladimir Vlasov, Elbert EN Macau, and Arkady Pikovsky,“Synchronization of oscillators in a Kuramoto-type model withgeneric coupling,” Chaos: An Interdisciplinary Journal of Non-linear Science , 023120 (2014).[39] Hidetsugu Sakaguchi and Yoshiki Kuramoto, “A soluble activerotater model showing phase transitions via mutual entertain-ment,” Progress of Theoretical Physics , 576–581 (1986).[40] Thomas K DM Peron, J¨urgen Kurths, Francisco A Rodrigues,Lutz Schimansky-Geier, and Bernard Sonnenschein, “Travel-ing phase waves in asymmetric networks of noisy chaotic at-tractors,” Physical Review E , 042210 (2016).[41] Shinya Watanabe and Steven H Strogatz, “Constants of motionfor superconducting josephson arrays,” Physica D: NonlinearPhenomena , 197–253 (1994).[42] Arkady Pikovsky and Michael Rosenblum, “Partially integrabledynamics of hierarchical populations of coupled oscillators,” Physical Review Letters , 264103 (2008).[43] Spase Petkoski, Dmytro Iatsenko, Lasko Basnarkov, and AnetaStefanovska, “Mean-field and mean-ensemble frequencies of asystem of coupled oscillators,” Physical Review E , 032908(2013).[44] Daniel M Abrams, Rennie Mirollo, Steven H Strogatz, andDaniel A Wiley, “Solvable model for chimera states of coupledoscillators,” Physical Review Letters , 084103 (2008).[45] Bastian Pietras, Nicol´as Deschle, and Andreas Daffertshofer,“Equivalence of coupled networks and networks with multi-modal frequency distributions: Conditions for the bimodal andtrimodal case,” Physical Review E , 052211 (2016).[46] Erik A Martens, “Bistable chimera attractors on a triangu-lar network of oscillator populations,” Physical Review E ,016216 (2010).[47] Seth A Marvel, Renato E Mirollo, and Steven H Strogatz,“Identical phase oscillators with global sinusoidal couplingevolve by m¨obius group action,” Chaos: An InterdisciplinaryJournal of Nonlinear Science , 043104 (2009).[48] Chen Chris Gong and Arkady Pikovsky, “Low-dimensionaldynamics for higher-order harmonic, globally coupled phase-oscillator ensembles,” Physical Review E , 062210 (2019).[49] Per Sebastian Skardal and Alex Arenas, “Abrupt desynchro-nization and extensive multistability in globally coupled oscil-lator simplexes,” Physical Review Letters , 248301 (2019).[50] Guilherme Ferraz de Arruda, Giovanni Petri, and YamirMoreno, “Social contagion models on hypergraphs,” PhysicalReview Research , 023032 (2020).[51] Ana P Mill´an, Joaqu´ın J Torres, and Ginestra Bianconi, “Ex-plosive higher-order Kuramoto dynamics on simplicial com-plexes,” Physical Review Letters , 218301 (2020).[52] Jan Frederik Totz, Julian Rode, Mark R Tinsley, KennethShowalter, and Harald Engel, “Spiral wave chimera statesin large populations of coupled chemical oscillators,” NaturePhysics , 282–285 (2018).[53] Dumitru C˘alug˘aru, Jan Frederik Totz, Erik A Martens, andHarald Engel, “First-order synchronization transition in a largepopulation of strongly coupled relaxation oscillators,” Scienceadvances , eabb2637 (2020).[54] Aaron M Hagerstrom, Thomas E Murphy, Rajarshi Roy, PhilippH¨ovel, Iryna Omelchenko, and Eckehard Sch¨oll, “Experimen-tal observation of chimeras in coupled-map lattices,” NaturePhysics8