Phase reduction beyond the first order: the case of the mean-field complex Ginzburg-Landau equation
aa r X i v : . [ n li n . AO ] J u l Phase reduction beyond the first order: The case of the mean-field complex Ginzburg-Landauequation
Iv´an Le´on and Diego Paz´o
Instituto de F´ısica de Cantabria (IFCA), CSIC-Universidad de Cantabria, 39005 Santander, Spain (Dated: July 23, 2019)Phase reduction is a powerful technique that makes possible to describe the dynamics of a weakly perturbedlimit-cycle oscillator in terms of its phase. For ensembles of oscillators, a classical example of phase reductionis the derivation of the Kuramoto model from the mean-field complex Ginzburg-Landau equation (MF-CGLE).Still, the Kuramoto model is a first-order phase approximation that displays either full synchronization or in-coherence, but none of the nontrivial dynamics of the MF-CGLE. This fact calls for an expansion beyond thefirst order in the coupling constant. We develop an isochron-based scheme to obtain the second-order phase ap-proximation, which reproduces the weak coupling dynamics of the MF-CGLE. The practicality of our methodis evidenced by extending the calculation up to third order. Each new term of the power series expansioncontributes with additional higher-order multi-body (i.e. non-pairwise) interactions. This points to intricatemulti-body phase interactions as the source of pure collective chaos in the MF-CGLE at moderate coupling.
I. INTRODUCTION
Networks of nonlinear elements with oscillatory behavior(‘oscillators’) are found in a variety of disciplines, such asneuroscience or engineering [1–4]. It is an empirical factthat some phenomena arising in these systems can be under-stood in terms of interacting phase oscillators. This frame-work has proven to be useful modeling and engineering exper-imental setups composed of many rhythmic elements, operat-ing in a wide range of spatio-temporal scales, and interactingthrough very different physical processes. We may cite smallmotors—cell phone vibrators—interacting through an elasticplate [5], networks of (electro-)chemical oscillators [6, 7], ar-rays of Josephson junctions [8, 9] and globally coupled elec-trical self-oscillators [10, 11], or nanoelectromechanical os-cillators in a ring [12].Applying a phase reduction method [1, 13–15] is the rig-orous way of describing a weakly perturbed oscillator solelyin terms of its phase (the other degrees of freedom becomeenslaved). However, obtaining analytically the approximate‘phase-only model’ for a specific system is not an easy task.Moreover, phase reduction becomes inaccurate unless the dis-turbances are not sufficiently weak. While, according to com-mon wisdom phase reduction of oscillator ensembles yieldspairwise interacting phase oscillators [13], multi-body (i.e.non-pairwise) interactions may also be relevant in some con-texts. Apart from the idea of invoking hypothetical three-bodyinteracting limit-cycle oscillators [16], multi-body phase in-teractions naturally arise if the coupling is nonlinear [17], seealso [18]. Instead, for linear pairwise coupling, three-bodyinteractions are a distinctive element of second-order phaseapproximations, as recently highlighted in [12]. Recognizingthe ubiquity of multi-body interactions may also be importantfor reconstructing phase interactions from data [19].Much of our knowledge on nonlinear dynamics relies onminimal models that capture the essential mechanisms behindcomplex phenomena. For oscillatory dynamics, the conven-tional test bed is the normal form of the Hopf bifurcationabove criticality: the so-called Stuart-Landau oscillator. Con-cerning geometry, global coupling is a fruitful simplifying as- sumption [1, 13, 20]. These two ingredients are combined ina standard model of collective dynamics: the fully connectednetwork of Stuart-Landau oscillators, or the mean-field ver-sion of the complex Ginzburg-Landau equation (MF-CGLE)[21–36]. This system is particularly interesting for chaos the-ory since it exhibits both microscopic (extensive) and macro-scopic (collective) chaos, either combined or independently,depending on parameters [21–24, 26, 30, 33, 36]. Phasereduction of the MF-CGLE yields the Kuramoto-Sakaguchimodel [13, 22, 37], a first-order approximation that behavesin a pathological way (unless heterogeneities are present): itonly displays full synchrony or incoherence. Therefore, purecollective chaos and other phase dynamics of the MF-CGLEremain to be analytically described in terms of a phase model.Such a phase reduction should provide additional insights intothe nature of collective chaos (playing an analogous role to theKuramoto-Sivashinsky equation of phase turbulence).The aim of this paper is two-fold: we introduce a phase-reduction method, and we investigate the phase model ob-tained from the MF-CGLE. The paper is organized as fol-lows. In Sec. II we reexamine the phase dynamics of theMF-CGLE and the connection with its first-order phase re-duction (the Kuramoto model). In section III, we present oursystematic phase-reduction procedure, based on the direct useof isochrons, which delivers a well-controlled power expan-sion in the coupling strength parameter. Section IV is devotedto investigating the weak coupling limit of the MF-CGLE bymeans of the the second-order phase reduction, which unfoldsthe degeneracies of the Kuramoto model; we address the casesof a large ensemble of oscillators, as well as an small oneof four oscillators. Section V presents the third-order con-tribution to the phase reduction of the MF-CGLE. Finally, inSec. VI we discuss the implications of our work and someoutlooks.
II. MEAN-FIELD COMPLEX GINZBURG-LANDAUEQUATION
The MF-CGLE consists of N diffusively coupled Stuart-Landau oscillators governed by N coupled (complex-valued) - - - Re ( A j ) I m ( A j ) a ) - - - Re ( A j ) I m ( A j ) b ) - - - Re ( A j ) I m ( A j ) c ) - - - Re ( A j ) I m ( A j ) d ) FIG. 1. Snapshot of the positions A j for a population of N = 200 oscillators with c = 3 . The corresponding mean field ¯ A is markedby a red cross, and a thin solid line is the trajectory of ¯ A ( t ) for aninterval of 50 t.u. (a) NUIS state with Q ≈ . ( c = − . , ǫ =0 . ), (b) Quasiperiodic partial synchrony ( c = − , ǫ = 0 . ).(c) Pure collective chaos ( c = − , ǫ = 0 . ). (d) Collective andmicroscopic chaos ( c = − , ǫ = 0 . ) for N = 500 . ordinary differential equations: ˙ A j = A j − (1 + ic ) | A j | A j + ǫ (1 + ic )( ¯ A − A j ) . (1)Here, A j = r j e iϕ j is a complex variable (index j runs from1 to N ), and the mean field is ¯ A = N − P Nk =1 A k . Apartfrom the population size N , there are three free parameters inEq. (1): ǫ , c , and c . Parameter ǫ , controlling the couplingstrength, is positive in order to preserve the analogy with the(spatially extended) Ginzburg-Landau equation. Parameter c introduces a cross-coupling between real and imaginary partsof the A j ’s. This non-dissipative coupling, so-called ‘reac-tive’ [4], generically appears from center manifold reduction[13]. Finally, ‘nonisochronicity’ (or ‘shear’) parameter c inEq. (1) determines the dependence of the angular velocity ofone oscillator on its radial coordinate. There are two impor-tant symmetries in system (1): invariance under a global phaseshift A j → A j e iφ , and full permutation symmetry stemmingfrom the mean-field coupling. A. Phenomenology
For many parameter values, the global attractor of Eq. (1) iseither full synchronization (FS) A j = ¯ A = e − ic t or one inco-herent state with vanishing mean field ¯ A = 0 . In the latter casethe oscillators rotate freely, A j = √ − ǫ exp { i [ − c + ǫ ( c − c )] t + φ j } . Among all the states compatible with ¯ A = 0 , themost prominent one is the uniform incoherent state (UIS) in which the φ j are uniformly distributed in the thermodynamiclimit (for a finite ensemble, the φ j are evenly spaced, deserv-ing the name of splay state or ponies on a merry-go-roundstate). A continuum of nonuniform incoherent states (NUISs)coexists with UIS, but usually arbitrarily weak noise spreadsthe phases and UIS is eventually attained. Nonetheless, forcertain parameters values, such as those in Fig. 1(a), the UISis unstable and one NUIS sets in spontaneously [21, 25].In addition to FS, UIS, and NUIS, system (1) exhibits a richrepertoire of collective states including clustering [22, 27–29, 33, 35], diffusion-induced inhomogeneity (or chimera)[28, 29, 32], quasiperiodic partial synchronization (QPS)[22, 36], as well as collective and microscopic chaos [21–24, 26, 30, 33, 36]. In a QPS state, see e.g. Fig. 1(b), the meanfield ¯ A rotates uniformly, while the individual oscillators be-have quasiperiodically (since each oscillators ‘feels’ the peri-odic driving of the mean field). Remarkably, increasing cou-pling QPS may undergo a couple of secondary Hopf bifurca-tions resulting in a state of pure collective chaos [24, 36]. Withthis term we refer to a state in which the mean field behaveschaotically, while individual oscillators behave in seeminglychaotic-like fashion (neighboring oscillators remain close forever due to the absence of microscopic chaos). A shared fea-ture of NUIS, QPS and pure collective chaos [22, 24, 36] isthat the relative positions of the oscillators on top of a closedcurve are preserved, see Figs. 1(a), 1(b) and 1(c). This factsuggests that a description in terms of oscillators’ phases aloneis possible. In contrast to Fig. 1(c), Fig. 1(d) shows a chaoticregime in which phase description breaks down, as it involvesmicroscopic degrees of freedom and no phase ordering exists.Hence, our ultimate goal is to find a phase-reduced modelof Eq. (1) that captures as much as possible of the phase-describable states (NUIS, QPS, modulated QPS, pure collec-tive chaos,etc.). B. Basic phase diagrams
Before presenting our results it is convenient to review pre-vious results on the MF-CGLE. For fixed c and c values, letus denote by ǫ s and ǫ , the ǫ values of marginal linear stabilityfor FS and UIS. Closed formulas for ǫ s and ǫ are [21, 22]: ǫ s = − c c )1 + c , (2) ǫ (2 ǫ − c + 4( ǫ − ǫ − c c − ǫ ( ǫ − c + (3 ǫ − = 0 . (3)These formulas are also valid for finite ensembles, assuming ǫ refers to the splay state. To visualize the stability bound-aries in Eqs. (2) and (3), it is convenient to fix either c or c .Following [22] we choose to fix c , and display the loci of ǫ s and ǫ in the parameter plane ( c , ǫ ) . In the phase diagramsin Figs. 2(a) 2(b) and 2(c) we selected c = 3 , c = 2 and c = 1 , respectively. This choice is motivated by the fact thatmost previous works on the MF-CGLE adopt either c = 2 or c = 3 . One key observation is that, as ǫ s and ǫ approachzero, the boundaries converge to the condition c c = 0 , FIG. 2. Partial phase diagram of the MF-CGLE for c = 3 (a), (b), and (c). In each panel, the region with stable UIS is depicted in yellow,and the region with color gradation corresponds to stable NUIS, with a color gradient that indicates the actual Q value (see text); it becomesdarker as Q → . Stable FS is indicated by a blue hatched region. The stability boundaries of FS, UIS and ( Q = 1 )-NUIS are depicted byblue, black and red lines, respectively; following Eqs. (2), (3), and (4) (setting Q = 1 ). In panels (a) and (b), there is green-hatched regionwhere other phase-describable states like the ones shown in Figs. 1(b) and 1(c) are stable. which is the well known Benjamin-Feir-Newell criterion forthe stability of uniform oscillations in the complex Ginzburg-Landau equation in arbitrary dimension [4, 13, 38, 39]. Thereis a critical value c = √ . . . . at which the bound-aries ǫ s and ǫ become tangent at ǫ = 0 . Accordingly, for c = 1 , see Fig. 2(c), there is a region of bistability betweenUIS and synchrony, in contrast e.g. to Fig. 2(b).The stability of a NUIS depends exclusively on the meanfield Q = | N − P j exp(2 iϕ j ) | . The coupling constant ǫ Q atwhich one particular NUIS becomes unstable was obtained inRef. [25]: (cid:2) ǫ Q (2 ǫ Q − c + 4( ǫ Q − ǫ Q − c c − ǫ Q ( ǫ Q − c +(3 ǫ Q − (cid:3) [(2 − ǫ Q ) + ǫ Q c ] (4) = Q ǫ Q (1 − ǫ Q )(3 ǫ Q − ( c + 1)( c + 1) . This formula is the generalization of (3) with the importantqualitative information that the size of the stability region in-creases as Q grows, reaching its maximum for Q = 1 . At Q = 1 the NUIS collapses into a two-cluster state with equallypopulated groups. The value of Q is still far from breakingthe degeneracy of a NUIS, provided Q = 1 , since the valuesof all ‘higher-order’ mean fields f n = | N − P j exp( niϕ j ) | ( n > ) are free. Nevertheless, the conclusion based on nu-merical simulation is that any small amount of noise causes f n to converge to zero, and Q to take the smallest value amongall non-unstable (i.e. neutrally stable) NUISs. Therefore, it isassumed hereafter that the term NUIS is constrained to f n = 0 ( n > ).Figures 2 (a) and (b) include a green-hatched region, ad-jacent to the UIS region at moderate ǫ values, where otherphase-describable states are stable. These are QPS, modu-lated QPS and pure-collective chaos [24, 36]. We determinedthe boundary through simulations with N = 200 oscillators,but the result is insensitive if a larger N value is used. C. First-order phase reduction: Kuramoto-Sakaguchi model
At the lowest order, applying the classical averaging tech-nique [4, 13, 37] to Eq. (1) yields the Kuramoto-Sakaguchimodel [40]. In this model, each oscillator is described by aphase θ j , and it is coupled to the other ones by pairwise in-teractions of the form sin( θ i − θ j + α ) . In agreement withthe mean-field character of the system, oscillators are cou-pled through the Kuramoto order parameter Z ≡ R e i Ψ = N − P Nk =1 e iθ k , such that the ordinary differential equationsgoverning the dynamics are: ˙ θ j = Ω + ǫη R sin(Ψ − θ j + α ) , (5)with constants Ω ≡ − c + ǫ ( c − c ) , η ≡ p (1 + c )(1 + c ) ,and phase lag α = arg[1 + c c + ( c − c ) i ] . (6)Equation (5) is the disorder-free version of the paradig-matic Kuramoto-Sakaguchi model [13, 40] and related mod-els [41]. The dynamics of Eq. (5) is determined by the signof c c (Benjamin-Feir-Newell criterion): full synchrony—corresponding to R = 1 — is stable for c c > , andunstable for c c < . In the latter case, among infinitelymany oscillator densities with R = 0 , there is a convergenceto the UIS under an arbitrarily weak noise [22].As discussed above, the MF-CGLE has much richer dy-namics than its first-order phase reduction (5), even arbitrarilyclose to the ǫ = 0 limit. Therefore, it is mandatory to ex-tend the phase reduction to order O ( ǫ ) if we wish to avoiddegeneracies in the phase approximation. This is what we donext. III. SYSTEMATIC PHASE REDUCTION
In spite of the relevance of Eq. (1) no phase reduction be-yond the first order is currently available. Finding higher orderterms in the phase reduction is necessary to unfold the singu-larity at ( c , ǫ ) = ( − /c , , see Fig. 2. This path of investi-gation should allow us to discern which are the true behaviorsof the MF-CGLE in the small coupling limit | ǫ | ≪ . More-over, it might serve to shed light on the mechanisms behindcomplex dynamics found (so far) for moderate ǫ values.An isochron-based phase reduction approach is developedhere. It allowed us to obtain the phase reduction of the MF-CGLE up to order ǫ . In this section we give the details of ourphase reduction calculation. We anticipate that the results atsecond and third order in ǫ correspond to Eqs. (15) and (29)below. A. Isochrons
The concept of isochron [42, 43] is the cornerstone ofphase reduction methods [1, 13]. Isochrons foliate the attrac-tion basin of a stable limit cycle, each intersecting it at onepoint. The phase of that point is attributed to all points of theisochron, motivated by their convergence as time goes to in-finity (the so-called ‘asymptotic phase’ [44]). For the Stuart-Landau oscillator, polar coordinates ( r, ϕ ) relate to the phase θ according to [4, 13]: θ ( r, ϕ ) = ϕ − c ln r. (7)As mentioned above, on the limit cycle ( r = 1 ), θ = ϕ .The term “nonisochronicity” or “shear” for parameter c be-comes clear in light of Eq. (7), since c controls how much theisochrons deviate from radial lines. B. Isochron-based phase reduction
We continue the analysis writing Eq. (1) in polar coordi-nates: ˙ r j = r j (1 − ǫ − r j )+ ǫN N X k =1 r k (cid:20) cos( ϕ k − ϕ j ) − c sin( ϕ k − ϕ j ) (cid:21) , (8) ˙ ϕ j = − c r j − ǫc + ǫN r j N X k =1 r k (cid:20) c cos( ϕ k − ϕ j ) + sin( ϕ k − ϕ j ) (cid:21) . (9)After the change of variables ( r j , ϕ j ) → ( r j , θ j ) throughEq. (7), we get: ˙ r j = f ( r j ) + ǫg j ( r , θ ) , (10a) ˙ θ j = ǫh j ( r , θ ) . (10b)Here, we have also implemented the transformation θ j → θ j − c t (by moving to a rotating frame with angular ve-locity − c ). In this way, the time derivatives of the phasesin (10b) are proportional to ǫ , while the r j are fast variablesthat become enslaved to the dynamics of θ j . In Eq. (10) f ( r ) = r (1 − r ) , and functions g j and h j depend on thevectors r = ( r , r , . . . , r N ) T and θ = ( θ , θ , . . . , θ N ) T asfollows, g j ( r , θ ) = − r j + 1 N N X k =1 (cid:26) r k (cid:20) cos (cid:0) θ k − θ j + c ln r k r j (cid:1) − c sin (cid:0) θ k − θ j + c ln r k r j (cid:1)(cid:21)(cid:27) , (11a) h j ( r , θ ) = c − c + 1 N r j N X k =1 (cid:26) r k (cid:20) ( c − c ) cos (cid:0) θ k − θ j + c ln r k r j (cid:1) + (1 + c c ) sin (cid:0) θ k − θ j + c ln r k r j (cid:1)(cid:21)(cid:27) . (11b)The separation of time scales in Eq. (10) suggests using clas-sical perturbation techniques like averaging, adiabatic ap-proximation, or two-timing. However, the perturbation ap-proach described next proved to be both conceptually sim-ple and much less convoluted, permitting us to obtain thephase reduction up to cubic order in ǫ . Based on the em-pirical observation that, at small ǫ values, the oscillators fallon a closed curve and preserve their phase ordering, we as-sume that the radii are completely determined by the phases r j = r j ( θ , θ , . . . , θ N ) . We also postulate an expansion inpowers of ǫ for the radii: r = r (0) j + ǫr (1) j + ǫ r (2) j + · · · ; orin vector notation r = r (0) + ǫ r (1) + ǫ r (2) + · · · . Equation (10b) for θ j becomes: ˙ θ j = ǫh j ( r (0) , θ ) + ǫ (cid:16) ∇ r h j ( r (0) , θ ) (cid:17) · r (1) + ǫ h(cid:16) ∇ r h j ( r (0) , θ ) (cid:17) · r (2) + ( Mrr ) j i + · · · , (12)where ∇ r ≡ ( ∂ r , ∂ r , . . . , ∂ r N ) and ( Mrr ) j ≡ P k,l ∂ r k ∂ r l h j ( r (0) , θ ) r (1) k r (1) l . Now, the explicit depen-dence on the radii in (12) must be removed. This is accom-plished equating both sides of (10a) at the same order. Theorder O ( ǫ ) yields r (0) j = 1 , and (12) becomes (at the lowestorder) the Kuramoto-Sakaguchi model (5). At order ǫ : ˙ r (1) j = f ′ ( r (0) j ) r (1) j + g j ( r (0) , θ ) . (13)As r j depends exclusively on the phases, we can apply thechain rule: ˙ r j = ( ∇ θ r j ) · ˙ θ . At order ǫ , the time derivativevanishes: ˙ r (1) j = ( ∇ θ r (0) j ) · h = 0 . Hence Eq. (13) yields the result r (1) j = − g j ( r (0) , θ ) f ′ ( r (0) j ) = g j ( r (0) , θ )2 , (14) which can be inserted in (12) to obtain the ǫ contribution.Through elementary manipulations the second-order phase re-duction of Eq. (1) can be condensed into this expression: ˙ θ j = Ω + ǫη R sin(Ψ − θ j + α ) + ǫ η ( R Q sin(Φ − Ψ − θ j ) − X m =1 ( − R ) m sin[ m (Ψ − θ j ) + β ] ) + O ( ǫ ) . (15)The O ( ǫ ) term depends on Z as well as on the sec-ond Kuramoto-Daido order parameter [45] Z ≡ Q e i Φ = N − P Nk =1 e iθ k . To enhance the clarity of Eq. (15), wefound it convenient to define a phase lag β = arg(1 − c + 2 c i ) , (16)which turns out to be independent of c . The other constants inEq. (15) are the same as in Eq. (5); as the change to a rotatingframe has been reversed, the O ( ǫ ) term inside Ω is − c (asbefore). IV. SECOND-ORDER PHASE REDUCTION:THREE-BODY INTERACTIONS
In this section we study in detail the phase model obtainedfrom the second-order phase reduction of the MF-CGLE,i.e. the system of phase oscillators governed by Eq. (15). Ofthe three O ( ǫ ) contributions to Eq. (15), the first element ofthe sum ( m = 1 ) entails a parameter shift to the O ( ǫ ) inter-action, and it is therefore irrelevant in qualitative terms. Theother two terms in Eq. (15) correspond to three-body (i.e. non-pairwise) interactions: R sin[2(Ψ − θ j )+ β ] = 1 N X k,l sin( θ k + θ l − θ j + β ) (17) R Q sin(Φ − Ψ − θ j ) = 1 N X k,l sin(2 θ k − θ l − θ j ) (18)The price of working only with the phases is that two-bodyinteractions of the original MF-CGLE (1) become multi-bodyinteractions, as higher orders of ǫ are considered. In compari-son to Eq. (1) our phase model can be much more efficientlyanalyzed, both analytically and numerically. We devote the re-mainder of this section to analyze the phase model in Eq. (15).We note that, as expected, the model is invariant under globalphase shift θ j → θ j + φ . For the sake of making the presenta-tion simpler we assume constant Ω = 0 , since this can alwaysbe achieved by going to a rotating frame θ j → θ j + Ω t . - - - - c (cid:1) a ) - - - - - c (cid:1) b ) FIG. 3. Comparison between the bifurcation lines of FS, UIS andNUIS ( Q = 1 ) for the MF-CGLE (solid lines) and for the second-order phase reduction (dashed lines). Line colors are the same as inFig. 2. Panel (a) c = 3 , and panel (b) c = 1 . A. Full synchronization
The stability boundary of FS ( θ j = Ψ = Φ / ) is eas-ily calculated. In particular, for infinite N it is almost im-mediate: we simply assume one oscillator is infinitesimallyperturbed, say the first one, θ = Ψ + δθ . The evolu-tion of the perturbation obeys the linear equation dδθ /dt = ǫη (cid:2) cos α + ǫη (1 − cos β ) (cid:3) δθ . At threshold ( dδθ /dt = 0 )the coupling satisfies: ǫ s = − c c ) c (1 + c ) (19)where we have written cos α and cos β in terms of c and c .For illustration, the curve defined by (19) is represented by ablue dotted line in Figs. 3(a) and (b) for c = 3 and c = 1 ,respectively. Equation (19) is asymptotically exact as ǫ s → ,and deviates progressively from the FS boundary of the MF-CGLE (represented by a solid line) as ǫ s increases. B. Incoherent states
We adopt the thermodynamic limit and define a density ρ such that ρ ( θ, t ) dθ is the fraction of oscillators with phasesbetween θ and θ + dθ . Now θ ∈ [0 , π ) is a cyclic variable,i.e. ρ ( θ + 2 π, t ) = ρ ( θ, t ) , and we impose the normaliza-tion condition R π ρ ( θ, t ) dθ = 1 . The oscillator density ρ obeys the continuity equation because of the conservation ofthe number of oscillators: ∂ t ρ ( θ, t ) + ∂ θ [ v ( θ ) ρ ( θ, t )] = 0 . (20) Here v = ˙ θ is the ρ -dependent velocity of an oscillator withphase θ . We define the Fourier modes of ρ : ρ ( θ, t ) = 12 π ∞ X n = −∞ ρ n e inθ , (21)with ρ = 1 and ρ n = ρ ∗− n . The mean fields Z n reduce to Z n = Z π ρ ( θ, t ) e inθ dθ = ρ − n . Inserting the Fourier expansion (21) into the continuity equa-tion (20) allows us to rewrite our model in Fourier space: ˙ ρ n = n ǫη n e − iα ρ ρ n − − e iα ρ ∗ ρ n +1 + ǫη (cid:2) e − iβ ρ ( ρ n − − ρ ρ n − ) − e iβ ρ ∗ ( ρ n +1 − ρ ∗ ρ n +2 ) − ρ ∗ ρ ρ n +1 + ρ ρ ∗ ρ n − (cid:3)o (22)
1. Uniform incoherent state
The stability boundary of the UIS ( ρ ( θ ) = (2 π ) − ⇔ ρ n =0 = 0 ) is obtained linearizing the previous equation. Itis easy to notice that only the first mode may destabilize. Wehave for | ρ | ≪ : ddt δρ = ǫη (cid:20) e − iα + ǫηe − iβ (cid:21) δρ . (23)Neglecting the trivial marginal case ǫ = 0 , the stability bound-ary satisfies cos α +(1 / ǫ η cos β = 0 . Or in terms of c and c : ǫ = 4(1 + c c )( c − c ) . (24)In Figs. 3(a) and 3(b), we can contrast this formula with theexact one for the MF-CGLE, Eq. (3), for two c values.
2. Nonuniform incoherent states
According to (22), in an incoherent state ( ρ = 0 ) higher-order modes are at rest: ˙ ρ n = 0 ( n > ). The linearization of(22) around ρ = 0 and ρ n = 0 ( | n | ≥ ) is (schematically)as follows: ddt δρ δρ ∗ δρ δρ ∗ ... = • • · · ·• • · · ·• • · · ·• • · · · ... ... ... ... . . . δρ δρ ∗ δρ δρ ∗ ... , (25)where the • symbols denote nonzero elements. Clearly, thestructure of this equation yields an infinite set of vanishingeigenvalues plus two eigenvalues coming from the first two rows. The equation for δρ is hence, the only relevant one.The linear terms in δρ yield: ˙ δρ = ǫ η nh e − iα + ǫη (cid:0) e − iβ − | ρ | (cid:1)i δρ − h e iα + ǫη e iβ − i ρ δρ ∗ o . (26)All higher-order modes, save ρ , are absent in the equation.As ˙ ρ = 0 , we can choose the coordinate axes such that ρ = Q ∈ R . After some calculations we find that NUIS with aspecific Q value is marginally stable at: ǫ Q = 4(1 + c c )( c − c ) + η Q . (27)As occurs in the MF-CGLE the larger Q is, the larger is thestability region of the NUIS. Our empirical observation is that,for given c and c , if ǫ is set at a certain ǫ = ǫ Q ∗ the nu-merical integration of the system (either oscillators or Fouriermodes), under a very weak noise, always converges to a NUISwith ρ n ≥ = 0 ; and, | ρ | = Q ∗ . In other words, the sys-tem adopts the minimum value of | ρ | among all allowed byEq. (27).The state Q = 1 ( R = 0 ) —the last NUIS to destabilize—is singular, not only because it is just a two-cluster state withtwo equally populated groups, but also because in contrast tothe other NUIS, the instability is not oscillatory. Eq. (26) takesthe form ˙ δρ ∝ a δρ − a ∗ δρ ∗ what yields an additional zeroeigenvalue corresponding to the direction Im( δρ ) = 0 . C. Validity and accuracy
From our previous results, we conclude that the phase re-duction (15) is free of degeneracies. The boundaries of FS,UIS and NUIS with different Q values do not overlap. As adouble-check of the correctness of our analysis, we verifiedthat the boundaries (19), (24), and (27) obtained through thephase reduction are tangent to the equivalent boundaries of theMF-CGLE, Eqs. (2), (3) and (4), at ǫ = 0 . In Fig. 3 we depicttogether the boundaries of FS, UIS, and ( Q = 1) -NUIS of theMF-CGLE (solid lines) and phase reduction to second order(dashed lines) for two values of the nonisochronicity: c = 3 and c = 1 . These plots permit us to identify the range of ǫ inwhich the second-order approximation is accurate. For c = 1 the approximate bifurcation lines are accurate up to ǫ ≈ . ,while this range is certainly smaller for c = 3 .For general c , c values, the prefactor ( ǫη ) n appearing forfirst ( n = 1 ) and second ( n = 2 ) orders suggests to extrapolatethe relative smallness of ǫ to other c , and c values. Thus, ifin Fig. 3(b) accuracy is good up to ǫη ≈ . η , and η ≈ , wepropose ǫη < . (28)as a conservative range of validity of the second-order approx-imation. Nevertheless, Eq. (28) must be regarded with somecaution, since the third-order contribution to the phase reduc-tion expansion is not exactly proportional to ( ǫη ) , see Sec. V. D. Quasiperiodic partial synchronization
The phenomenon of QPS was originally reported in theMF-CGLE [24] as a state emerging from the destabilization ofthe UIS, see Fig. 1(b), though its finding is usually attributedto a model of phase oscillators [46]. As mentioned above, ina QPS state the mean-field rotates uniformly, but individualoscillators behave quasiperiodically. Each oscillator passesperiodically through a bottleneck located at the phase arg( ¯ A ) .The onset of QPS looks like a Hopf bifurcation undergone bythe UIS, but this is not the case because of the infinitely manyneutral directions pointing to nearby NUISs. It is also impor-tant to emphasize that stable QPS does not settles any time thatthe UIS becomes unstable. As can be appreciated in Fig. 2(a)and (b), QPS is only observed at moderate ǫ values when en-tering inside the green hatched region. Otherwise, what weobserve in the MF-CGLE is that the QPS state born at the in-stability of the UIS is a saddle. For parameter values with un-stable UIS and FS —outside the green-hatched regions— ini-tial conditions close to the UIS approach QPS for long time,eventually converging to one NUIS. If any small amount ofnoise is present, the NUIS with the lowest Q among the non-unstable ones is selected. The same behavior is displayed bythe second-order phase reduction, Eq. (15); see Fig. 4. Thelogarithmic scaling of the residence times near QPS indicatesa heteroclinic connection between UIS and QPS. The ampli-tude of the saddle QPS depends on the particular parametervalues. The state of QPS progressively grows as we moveaway from the UIS stability boundary, finally colliding withFS ( | ρ n | → ) at the point where FS becomes stable.All in all, these results confirm the correctness of our expan-sion, but at the same time prove the limitations of the second-order reduction, since the QPS attractor—found at moderate ǫ values, see Fig. 2(a) and (b)—is not reproduced. Time ( ) R T ( ) - - - R FIG. 4. Evolution of R ( t ) for N = 1000 phase oscillators governedby (15) initiated near the UIS state. For the selected parameter values( c = 3 , c = − . , ǫ = 0 . ) UIS is unstable but there are neu-trally stable NUISs. After a transient in the neighborhood of QPS( R ( t ) ≈ const . ), the system approaches a particular NUIS ( R = 0 ).From left to right the initial conditions are random perturbations ofthe UIS with R = 4 . × { − , − , − , − } . The ori-gin of times was shifted in all data sets to make the initial rise of R coincident. The inset shows the QPS transient time as a function of R . Note the logarithmic divergence of transient time T ∼ ln R (consistent with heteroclinicity). E. Clustering
Clustering is a much studied phenomenon in oscillator en-sembles [47]. In a clustered state there are several groups ofoscillators, each group formed by oscillators sharing the samephase. This kind of states are always possible in a mean-fieldmodel, so the relevant question is the stability. Indeed, theMF-CGLE is known to exhibit stable cluster at certain pa-rameter ranges [21, 22, 27–29, 33, 35], specifically for strongcoupling ( ǫ ≈ ).Are there stable clustered states at small coupling? Ourphase model allows us to address this question in an analyticalway. Nonetheless, the general problem is intractable and wedecided to restrict our study to states with two point clusters,where a fraction p of the population is in the A-state θ A , andthe remaining (1 − p ) fraction is in the B-state θ B = θ A .We now summarize the results; the corresponding calculationscan be found in Appendix I.As an illustrative example, Fig. 5 depicts the combinationsof phase difference ∆ = θ A − θ B and imbalance p corre-sponding to actual cluster solutions for three different c val-ues with fixed values of c and ǫ . Each panel is a typical situa-tion in a specific region of parameter space. At the FS thresh-old, between panels (a) and (b), there is an infinity of two-cluster solutions colliding with FS ( ∆ = 0 ). In consequencethere is a reconnection of the two-cluster solutions. In Fig. 5 - (cid:1) - (cid:1) / (cid:1) / (cid:1) (cid:2) p a ) - (cid:1) - (cid:1) / (cid:1) / (cid:1) (cid:2) p b ) - (cid:1) - (cid:1) / (cid:1) / (cid:1) (cid:2) p c ) FIG. 5. Two-cluster solutions of Eq. (15). Each panel representsthe fraction p of oscillators in one cluster as a function of thephase lag between clusters ∆ = θ A − θ B . We fix c = 3 and ǫ = 0 . and select three values of c in each panel. The arrowindicates the direction of increasing c . Solid (dashed) lines indi-cate stable (unstable) locking of the clusters. (a) Unstable FS region, c = − . , − . , − . ; (b) stable FS and not unstable ( Q = 1 )-NUIS region, c = − . , − . , − . ; (c) stable FS and unstable( Q = 1 )-NUIS region, c = − . , − . , . . solid lines represent stable locking of the clusters. However,these solutions are fragile against disintegration of the largestcluster. Our conclusion after an extensive exploration of pa-rameter space is that stable two-cluster states are not stableat small coupling. To be more precise, what we observe inour second-order phase reduction, Eq. (15), is that stable clus-tering is hardly found, and if so, it always requires moderatecoupling strengths, violating (28). And indeed, we could notreplicate clustering in the MF-CGLE for the parameter valuespredicted by Eq. (15).The stability analysis of the two-cluster solutions also con-firmed that slow switching [48] —a stable heteroclinic con-nection between two configurations of ∆ = θ A − θ B — is notpossible. F. Finite population, N = 4 This work focuses on the behavior of the MF-CGLE in thethermodynamic limit ( N → ∞ ). But the phase reduction (15)is valid for an arbitrary population size. In this section weconstruct a bifurcation diagram for N = 4 oscillators, onesize previously considered in the MF-CGLE context [35, 49].Here, this choice is motivated by the fact that in globally cou-pled systems this is the smallest size with a continuum ofstates with R = 0 [50], equivalent to the NUISs for N = ∞ . In analogy with its thermodynamic limit, the finite- N theKuramoto-Sakaguchi model has an exceptional transition be-tween FS and the splay state at c c = 0 . This degen-eracy can be broken down, for instance, adding higher-orderharmonics to the (pairwise) interactions [51]. In our case, de-generacy is broken down by the three-body interactions of thesecond-order in the phase-reduction expansion.Working with a small number of oscillators has the advan-tage that we can track all the stationary solutions, in partic-ular the clustered solutions. As there are orderings forthe oscillators phases, and phase ordering is preserved by thedynamics because of the mean-field interactions, we choosethe oscillators’ labels such that θ j ≤ θ j +1 . (We assumehere θ j ∈ [0 , π ) to avoid artificial degeneracies.) The setof phases { θ , θ , θ , θ } , may take several invariant config-urations. Apart from the trivial FS state { a, a, a, a } , thereexists a continuum of “NUIS-like” Z -symmetric states with { a, a + b, a + π, a + b + π } , where b ∈ (0 , π/ ∪ ( π/ , π ) .In the limit b → π/ the NUIS becomes the splay state (theanalogous of UIS). In addition, in the limits b → and b → π the NUIS collapses into a 2-cluster state with opposite phases.Apart from this one, other 2-cluster solutions are possible.Namely, for some parameter values there exist two symmetry-related 2-2 configurations { a, a, b, b } ( b = a + π ). Addi-tionally, one 3-1 cluster exists: designated as { a, a, a, b } or { b, a, a, a } . Three-cluster solutions, like { a, a, b, c } , do notexist in our phase reduction, in contrast to the MF-CGLE forstrong coupling [35]. Concerning the N = 4 analogous ofQPS, it is a periodic orbit, in which, due to the finiteness ofthe population, R and Q fluctuate around their average values.We use R , Q , and c to plot the bifurcation diagram inFig. 6. These coordinates have the drawback of collapsingmultiple equivalent states to a single point, hiding symmetries(e.g. pitchfork bifurcations). However, our choice intends toease the comparison with the previous section, and with thesame aim states are labeled borrowing the infinite- N terminol-ogy; namely, we use the labels UIS, NUIS, and QPS instead ofsplay state, Z -symmetric state, and limit cycle, respectively.Due to permutation symmetry FS destabilizes at point T inFig. 6, as three eigenvalues go through zero simultaneously.This comprises an equivariant transcritical bifurcation withthe - cluster, as well as a pitchfork bifurcation involving a - cluster. Moreover, at point T, QPS collapses into a hetero-clinic cycle. This coincidence of bifurcations is a known sce-nario in systems with full permutation symmetry [50]. Con-cerning UIS, it undergoes an oscillatory instability at point U,but this is not a standard Hopf bifurcation because of the neu-tral direction along the NUIS manifold. QPS is a saddle, andnot a stable limit cycle as it might have been naively expected.In Fig. 6 we took c = 3 , and the QPS branch connects pointsT and U in a simple way. In contrast to Fig. 6, for c = 1 FSand UIS coexist, and points U and T switch their relative posi-tions. In that case the QPS branch is completely reversed (notshown), and the QPS solution is fully unstable. Consistently,we found a range of c values in between, < c < , where(depending on ǫ ) the QPS branch develops a fold.In summary, the bifurcation diagram for N = 4 appearsto capture the global picture of the transition from UIS to FS. R - - - - c Q UIS NUIS QPSFS 2 - - FIG. 6. Bifurcation diagram for Eq. (15) with N = 4 oscillatorsand c = 3 , ǫ = 0 . . Solid (dashed) lines represent stable (unsta-ble) solutions. In the case of UIS and NUIS the solution depictedmust be understood as the one observed under arbitrarily weak noise(there is continuum of neutral solutions with Q larger than the solu-tion depicted). The saddle QPS orbit was continued by means of aNewton-Raphson algorithm, and the values of R and Q assigned inthe diagram correspond to their time averages. Considering more oscillators will increase the number of clus-ter solutions, see [27], but no essential new features.
V. THIRD-ORDER PHASE REDUCTION: FOUR-BODYINTERACTIONS
Our reason to deal with the third-order term now is to illus-trate the practicality of the phase reduction method, and get aglimpse of the power series expansion at higher orders. Eval-uating the cubic term in Eq. (12) yields the O ( ǫ ) correctionto Eq. (15): ǫ c (cid:26) C R sin(Ψ − θ j + γ ) + C R sin [2(Ψ − θ j ) + γ ] + C R Q sin(Φ − Ψ − θ j + γ ) + C R Q sin(Ψ − θ j + γ )+ C R sin(Ψ − θ j + γ ) + C R Q sin(Φ − θ j + γ ) + C R sin [3(Ψ − θ j ) + γ ] + C R P sin(Ξ − − θ j + γ )+ C R Q sin(Φ −
2Ψ + γ ) + DR (cid:27) . (29)This expression depends on the third Kuramoto-Daido or-der parameter Z ≡ P e i Ξ = N − P j e i θ j . The dependenceof constants { C j , γ j } j =1 ,..., and D on c and c is tabulatedin Appendix II. The structure of Eq. (29) deserves some wordshere. The terms proportional to C j with indices j = 1 , , arehigher-order corrections to Eq. (15), tantamount to a shift inparameter values. Four-body interactions appear in five dif-ferent forms, corresponding to indices j = 4 , . . . , . For illus-tration we expand a couple of these four-body contributions: R sin(Ψ − θ j ) = 1 N X k,l,n sin( θ k + θ l − θ n − θ j ) ,R P sin(Ξ − − θ j ) = 1 N X k,l,n sin(3 θ k − θ l − θ n − θ j ) . There are several qualitative features in Eq. (29) that deserveto be pointed out:1. The overall O ( ǫ ) contribution is not proportional to η —though some terms indeed are— in contrast to O ( ǫ ) and O ( ǫ ) , which are proportional to η and η , respec-tively.2. From Eqs. (15) and (29) we can expect that truncationof the power series to order ǫ n yields up to ( n + 1) -body interactions, but not higher-order non-pairwisecouplings. We can also expect that only Kuramoto-Daido order parameters Z k with k ≤ n appear.3. The last two terms in Eq. (29) are somewhat unex-pected, (nonetheless see [17]), since they depend onthe mean fields Z and Z , but not on θ j itself. Theyare hence irrelevant concerning synchronization bound-aries.4. As occurs with the O ( ǫ ) term, FS and (N)UIS statesare consistent with the MF-CGLE dynamics: (i) allterms in (29) are proportional to R ensuring that thecontribution to the oscillators’ frequencies vanishes in0one incoherent state; (ii) in the FS state, the contribu-tion also vanishes, as expected since the frequency ofFS in the MF-CGLE varies linearly with ǫ . Accord-ingly, it holds that D + P j C j sin γ j = 0 , cf. AppendixII.Unfortunately, there is not a recognizable pattern in the newterms appearing in the power series expansion, so it is notpossible to extrapolate to higher orders in ǫ .From Eqs.(15) and (29) we can derive the stability bound-ary of FS, NUIS (for UIS just set Q = 0 ) obtaining: c c ) + ǫ s c (1 + c ) + ǫ s c c (1 + c ) = 0 , (30) c c ) + ǫ Q (1 + c ) (cid:20) (1 − c ) − Q (1 + c ) (cid:21) + ǫ Q c ) (cid:20) (2 − c − c c + c c ) − Q (1 + c )( − c c ) (cid:21) = 0 . (31)In Fig. 7 we depict (a) ǫ , and (b) ǫ s from the previous ex-pressions and compare them with the result of the MF-CGLE,and with the second-order approximation. The slopes and thecurvatures of the bifurcation lines of the third-order phase re-duction agree with those of the MF-CGLE at ǫ = 0 . VI. DISCUSSIONA. Alternative phase reduction(s)
Our phase reduction is a genuine power series in the smallparameter ǫ . Another strategy to analyze (1) is to absorb the ǫA j term prior to the phase reduction. Specifically, setting t ′ = (1 − ǫ ) t , and κ = ǫ − ǫ , (32)we get dB j dt ′ = B j − (1 + ic ) | B j | B j + κ (1 + ic ) ¯ B, (33)where B j = A j exp( iǫc t ) / √ − ǫ . Applying our phase-reduction method to (33) we obtain an alternative phase re-duction in powers of κ (the result is not qualitatively differ-ent).Is it worth transforming (1) into (33)? In other words, is thephase reduction of (33) up to order κ n , superior to that of (1)at order ǫ n ? Certainly, phase reductions at order ǫ n and κ n arenot equivalent since κ = ǫ + ǫ + ǫ + · · · . Any truncationat order κ n involves all powers of ǫ . The relative accuracy ofthe phase reductions of (33) and (1) at the same order can beassessed comparing the bifurcation loci. Instead of applyingphase reduction to Eq. (33), the quickest strategy is to assumethe existence of an exact phase reduction involving all orders - / - - - c (cid:1) a ) - / - - - c (cid:1) b ) FIG. 7. Stability boundaries of (a) UIS and (b) FS obtained exactlyand from phase approximations, for c = 3 . The solid line corre-spond to the exact boundary of the MF-CGLE (1), while dotted anddashed lines correspond to second- and third-order phase approxi-mations, respectively. Blue lines are obtained from (15) and (29).Orange lines are the results if prior to phase reduction the MF-CGLEis transformed into (33), performing an isochron-base phase reduc-tion in powers of κ . in ǫ such that the exact critical value ǫ ∗ (the asterisk denotesan arbitrary state: UIS, FS, ...) satisfies: ∞ X n =1 a n ( c , c ) ǫ n ∗ + 1 + c c = 0 . (34)The coefficients a n depend on the specific instability we areconsidering.Phase reduction of (1) up to order n results in a truncationof (34) to order n − . For instance, the second-order phasereduction of (1) yields the linear relation [recall Eqs. (19) or(24)]: a ǫ ∗ + 1 + c c = 0 . (35)At the same order, the phase reduction of (33) results in ananalogous expression a ′ κ ∗ + 1 + c c = 0 . (36)Given that κ = ǫ + O ( ǫ ) , consistency with (34) de-termines a ′ = a . Thus the bifurcation locus estimatedfrom the phase reduction of (33) satisfies (in coordinate ǫ ) a ǫ ∗ / (1 − ǫ ∗ ) + (1 + c c ) = 0 , which is slightly differentfrom (35). Analogous reasoning permits us to obtain the bi-furcation lines for the third-order phase reduction of (33) fromEqs. (30) and (31).A comparison of the bifurcations lines of UIS and FS is dis-played in Figs. 7(a) and 7(b) for c = 3 . We see that the trans-formation of (1) into (33) allows us to obtain a phase modelthat captures better the stability boundary of UIS, but not ofFS. It is easy to understand why. Each strategy captures betterthe dynamics in which the quantities multiplying the couplingconstant are small. Thus, Eq. (1) is already a good startingpoint for states close to FS ( A j ≈ ¯ A ), while (33) works bet-ter close to incoherence ( ¯ A ≈ ¯ B ≈ ). Finally, note that inaddition to (1) and (33), there exists a continuum of alterna-tive, intermediate formulations, in which ǫA j is only partlyabsorbed by a coordinate transformation.1 B. Possible extensions of this work
The phase-reduction procedure presented in this work canbe easily implemented in other geometries, different from thefully connected network. In a networked architecture, phasereduction at first order in ǫ couples phases with the nearest-neighbor phases. At order ǫ , second nearest neighbors comealso into play [19], and progressively more distant nodes par-ticipate in the phase dynamics at higher orders. Also, the caseof non-locally coupled Stuart-Landau oscillators [52] is ana-lyzable with the phase reduction presented here. Concerningthe original complex Ginzburg-Landau equation, a partial dif-ferential equation of reaction-diffusion type, our phase reduc-tion procedure is very simple and efficient obtaining the coef-ficients of the second-order terms: ∇ θ , ( ∇ θ ) , ( ∇ θ ) ∇ θ , ∇ θ ∇ θ [13].Concerning the oscillator dynamics, the phase reductioncarried out here can be easily applied to planar oscillators withpolar symmetry ( λ − ω systems). In the latter case, analo-gously to (7), the isochrons satisfy θ = ϕ + χ ( r ) [1]. Even iffunction χ does not have a closed form, it is still possible toobtain the phase model using the derivatives of the isochronson the limit cycle. C. Relationship with other phase-reduction approaches
In this subsection, we comment on the progress of ourphase-reduction approach with respect to previous works,even if only directly applicable to λ − ω systems.An alternative way of obtaining the second-order phase re-duction of the MF-CGLE, Eq. (15), is applying the systematicaveraging formulation in Chap. 4 of the book by Kuramoto[13]. This calculation is, however, much more lengthy thanthe one presented in Sec. III. Not surprisingly, obtaining theorder ǫ with the averaging approach [13] is a totally imprac-tical task, while we succeeded with our method (with the as-sistance of symbolic software); see Eq. (29).Equation (15) can also be obtained assuming small varia-tions of the radii, i.e. setting ˙ r j = 0 . This procedure wasfollowed in [12, 53], with the difference being that there thesmall quantities are deviations from the reference limit cycle.Here, we pursued a bona-fide power expansion in terms of thecoupling constant ǫ , and the result differs from the one ob-tained following [12, 53]. In passing, we mention that insteadof assuming r j = r j ( θ , θ , . . . , θ N ) , once Eqs. (10) are de-rived, the two-timing approximation, such that the θ j dependonly on a slow time τ = ǫt , can also be applied.In contrast to our work, Refs. [17, 54] apply first-orderphase reduction obtaining multi-body phase interactions. Thereason is that those works invoke amplitude equations for anensemble close (but not asymptotically close) to a Hopf bi-furcation. The amplitude equation, which can be seen as ageneralization of Eq. (1), turns out to contain nonlinear in-teractions. The nonlinear coupling among the A j ’s leads tomulti-body interactions in the first-order phase reduction. Ap-plying second-order phase reduction, as described here, to theamplitude equations in [17] or [54] may be interesting. D. Towards a minimal phase model of pure collective chaos
Pure collective chaos has been found in several phase mod-els with heterogeneity [60] or delay [55]. Collective chaos inthe MF-CGLE, see Fig. 1(c), calls for a phase description interms of identical phase oscillators (without delays). The factthat we have not found evidence of collective chaos in our nu-merical simulations of the second-order phase reduction (15)—nor in the third-order one— can be reasonably attributed toa too restrictive truncation of the power expansion. We believethat a higher-order truncation will capture better the behaviorof the system at larger ǫ values, and eventually, will exhibitcollective chaos.As pairwise interactions through higher harmonics, like Q sin(Φ − θ j ) = N − P k sin[2( θ k − θ j )] , do not showup in the phase reduction of the MF-CGLE [56], multi-bodyphase interactions appear to be the most promising ingredi-ent to model collective chaos. In small ensembles of identicalphase oscillators, higher harmonics as well as multi-body in-teractions promote chaos alike, see [57] and [58], respectively.However, so far, collective chaos remains elusive in popula-tions of higher-order pairwise interacting identical phase os-cillators [59]. We believe multi-body interactions could be thekey element of collective chaos, instead.In the MF-CGLE with parameter values close to those inFig. 1(c), we found chaos with a population sizes as small as N = 6 . Does this say something about the order of the multi-body interactions needed in the phase reduction? Is this chaosconnected to collective chaos in the thermodynamic limit, asin [60]? E. Conclusions
Multi-body interactions are an unavoidable consequence ofphase reduction, but save for a few works [12, 16, 58, 61, 62],the role of multi-body phase interactions shaping exotic dy-namics remains largely unexplored. In the weak-couplingregime of the MF-CGLE, multi-body phase interactions areessential to describe all states apart from FS and UIS.In summary, in this work we achieve second- and third-order phase reductions of the MF-CGLE. In our view, higher-order phase reductions promise to be crucial for our under-standing of collective chaos and other exotic phenomena [12].Moreover, analytic higher-order phase reductions may alsoserve as test beds for numerical phase reductions recently im-plemented [63]. For these reasons, we regard phase reductionbeyond the first order as an exciting battleground of nonlineardynamics.
ACKNOWLEDGMENTS
We acknowledge support by MINECO (Spain) underProject No. FIS2016-74957-P. IL acknowledges support byUniversidad de Cantabria and Government of Cantabria un-der the Concepci´on Arenal programme.2
APPENDIX I: CLUSTERING
Our model (15) in a more convenient form (recall that inthe rotating frame
Ω = 0 ) reads: ˙ θ i = 1 N X j Γ( θ j − θ i )+ 1 N X j,k [ G ( θ j + θ k − θ i ) + G (2 θ j − θ k − θ i )] , (37)with Γ( x ) = ǫ (cid:20) ( c − c ) cos x + (1 + c c ) sin x (cid:21) − G ( x ) , (38) G ( x ) = − ǫ (1 + c ) (cid:20) c x + (1 − c )4 sin x (cid:21) , (39) G ( x ) = ǫ (1 + c )(1 + c )4 sin x. (40) Note that G is an odd function.Let us write first the evolution equation for cluster- A phase θ A , defining the phase difference ∆ = θ A − θ B : ˙ θ A = [ p Γ(0) + (1 − p )Γ( − ∆)] + (cid:2) p G (0) + 2 p (1 − p ) G ( − ∆) + (1 − p ) G ( − (cid:3) + (cid:2) − p (1 − p ) G (2∆) + ( − p + 3 p − G (∆) (cid:3) ; (41)the equivalent equation for the B -cluster is obtained with the substitution ∆ → − ∆ and p → (1 − p ) . The evolution of ∆( t ) obeys ˙∆ = (2 p − − p )Γ( − ∆) − p Γ(∆) + (cid:2) − p (1 − p ) G (2∆) + ( − p + 4 p − G (∆) (cid:3) + (cid:8) (2 p − G (0) + 2 p (1 − p ) [ G ( − ∆) − G (∆)] + (1 − p ) G ( − − p G (2∆) (cid:9) . (42)Setting ˙∆ = 0 we obtain a quadratic equation in p that canbe solved explicitly. We depict p (∆) in Fig. 5 for selectedparameter values. Note the symmetry of the curves becauseof the invariance under (∆ , p ) ↔ ( − ∆ , − p ) . There are ∆ values for which p is out of the range (0 , , indicating notwo-cluster states with those particular ∆ values exist. Con-versely, different values of ∆ may be consistent with the same p value, indicating the coexistence of multiple two-cluster so-lutions with the same sizes.
1. Stability
First of all, note that, one zero eigenvalue is always presentdue to the global phase shift invariance of the model, θ j → θ j + const . , and we ignore it hereafter. For the analysisthat follows it is simpler to assume the thermodynamic limit(eigenvalues are the same for finite N , but the calculation ismore convoluted.) As already known from previous studies[64], perturbations on a two-cluster solution can be decom-posed in three orthogonal modes. Two of them are the dis-integration of each respective cluster, and the third one is theunlocking of the two clusters. We denote λ A , λ B and λ L thecorresponding eigenvalues. For the stability of the A -cluster,we need to evaluate if one oscillator in the neighborhood ofthis cluster decays to it or departs (i.e., “evaporates”). The eigenvalue λ A is simply obtained linearizing around the state.The result is: λ A = − p Γ ′ (0) − (1 − p )Γ ′ ( − ∆) − p G ′ (0) − p (1 − p ) G ′ ( − ∆) − − p ) G ′ ( − − p G ′ (0) − (1 − p ) G ′ (∆) − p (1 − p ) G ′ (2∆) . (43)The eigenvalue λ B is obtained from λ A after the substitution p → (1 − p ) and ∆ → − ∆ , and viceversa. Finally, the lockingbetween the clusters is controlled by the eigenvalue obtainedlinearizing (42): λ L = − (1 − p )Γ ′ ( − ∆) − p Γ ′ (∆) − − p ) G ′ ( − − p (1 − p ) (cid:2) G ′ ( − ∆) + G ′ (∆) (cid:3) − p G ′ (2∆) − (4 p − p + 1) G ′ (∆) − p (1 − p ) G ′ (2∆) . (44)Stability requires λ A , λ B , λ L < . For small ǫ we summa-rized our findings in the main text, distinguishing three dif-ferent regions corresponding to the three panels of Fig. 5. Assaid in the main text we found stable clusters in the first re-gion (FS unstable), e.g. ǫ = 0 . , c = − , c = 2 . However,for these parameters the condition (28) does not hold, and infact the cluster solution destabilized when we implemented itin the MF-CGLE.3
2. No slow switching
With unstable two-cluster states, the system might still ex-hibit one nontrivial phenomenon called slow switching [48].In this phenomenon, the clusters switch between two differ-ent ∆ values with identical p value. The explanation for thisbehavior is a stable heteroclinic connection between the pairof two-cluster states that causes the system to switch for everbetween them with increasing residence times [48, 64]. Inpractice [59], switching either terminates in one of the unsta-ble two-cluster states (due to round-off errors), or it achievesa constant periodic switching (due to small noise). Accordingto [64], slow switching requires the coexistence of three two-cluster states ∆ ′ , ∆ ′′ , and ∆ ′′′ with an identical p value, suchthat < ∆ ′ < ∆ ′′′ < ∆ ′′ < π , and λ L < for ∆ ′ and ∆ ′′ , while λ L > for ∆ ′′′ . As may be seen in Fig. 5, findingparameter values with three solutions for ∆ at the same p isaready difficult —only for the green line in Fig. 5(a) do such p values exist. In addition, the condition for the eigenvaluesis even more stringent: e.g., in Fig. 5(a) the three points sharethe stability of λ A and λ B making the heteroclinic connectionimpossible. APPENDIX II: CONSTANTS C j , γ j AND D IN EQ. (29)
The dependences of constants { C j , γ j } j =1 ,..., on c and c is tabulated as follows: C j = q A j + B j , (45) γ j = arg( A j + iB j ) , (46) where A = 2( c c − c − c c + 2) A = − (3 c c − c − c c + 5) A = − c + 1)(2 c c − A = 2( c + 1)( c c + 1) A = 2( − c c + c + 2 c c A = 3( c + 1)( c c − A = − ( − c c + 9 c + 15 c c − A = (1 + c ) (5 c c − A = (1 + c )(1 + c c ) , and B = 2(4 c + c − c c ) B = ( c + 9 c c − c − c ) B = 2 c ( c + 1) B = 2( c + 1)( c − c ) B = ( c + 5 c c − c − c ) B = − c + 1)( c + c ) B = (cid:0) − c − c c + 9 c + 5 c (cid:1) B = (cid:0) c + 1 (cid:1) ( c + 5 c ) B = (1 + c )( c − c ) . Additionally, D = ( c + 1)( c − c ) (47) [1] A. T. Winfree, The Geometry of Biological Time (Springer, NewYork, 1980).[2] F. C. Hoppensteadt and E. M. Izhikevich,
Weakly connectedneural networks. (Spinger Verlag, N.Y., 1997).[3] S. H. Strogatz,
Sync: The emerging science of spontaneous or-der. (Hyperion Press, New York, 2003).[4] A. S. Pikovsky, M. G. Rosenblum, and J. Kurths,
Synchroniza-tion, a Universal Concept in Nonlinear Sciences (CambridgeUniversity Press, Cambridge, 2001).[5] D. Mertens and R. Weaver, “Synchronization and stimulatedemission in an array of mechanical phase oscillators on a reso-nant support,” Phys. Rev. E , 046221 (2011); “Individual andcollective behavior of vibrating motors interacting through aresonant plate,” Complexity , 45–53 (2011).[6] I. Z. Kiss, C. G. Rusin, H. Kori, and J. L. Hudson, “Engineeringcomplex dynamical structures: Sequential patterns and desyn-chronization,” Science , 1886–1889 (2007).[7] J. F. Totz, R. Snari, D. Yengi, M. R. Tinsley, H. Engel, andK. Showalter, “Phase-lag synchronization in networks of cou-pled chemical oscillators,” Phys. Rev. E , 022819 (2015).[8] K. Wiesenfeld and J. W. Swift, “Averaged equa-tions for Josephson junction series arrays,”Phys. Rev. E , 1020–1025 (1995).[9] K. Wiesenfeld, P. Colet, and S. H. Strogatz, “Synchro- nization transitions in a disordered Josephson series array,”Phys. Rev. Lett. , 404–407 (1996).[10] A. A. Temirbayev, Z. Zh. Zhanabaev, S. B. Tarasov,V. I. Ponomarenko, and M. Rosenblum, “Experimentson oscillator ensembles with global nonlinear coupling,”Phys. Rev. E , 015204(R) (2012).[11] L. Q. English, Z. Zeng, and D. Mertens, “Experimen-tal study of synchronization of coupled electrical self-oscillators and comparison to the Sakaguchi-Kuramoto model,”Phys. Rev. E , 052912 (2015).[12] M. H. Matheny, J. Emenheiser, W. Fon, A. Chapman, A. Sa-lova, M. Rohden, J. Li, M. Hudoba de Badyn, M. P´osfai,L. Duenas-Osorio, M. Mesbahi, J. P. Crutchfield, M. C.Cross, R. M. D’Souza, and M. L. Roukes, “Exotic statesin a simple network of nanoelectromechanical oscillators,”Science , eaav7932 (2019).[13] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, Berlin, 1984).[14] H. Nakao, “Phase reduction approach to synchronisation ofnonlinear oscillators,” Contemp. Phys. , 188–214 (2016).[15] B. Monga, D. Wilson, T. Matchen, and J. Moehlis, “Phase re-duction and phase-based optimal control for biological systems:a tutorial,” Biol. Cybern. , 11–46 (2019).[16] T. Tanaka and T. Aoyagi, “Multistable attractors in a net- work of phase oscillators with three-body interactions,”Phys. Rev. Lett. , 224101 (2011).[17] P. Ashwin and A. Rodrigues, “Hopf normal form with sn sym-metry and reduction to systems of nonlinearly coupled phaseoscillators,” Physica D , 14 – 24 (2016).[18] M. Rosenblum and A. Pikovsky, “Self-organized quasiperiod-icity in oscillator ensembles with global nonlinear coupling,”Phys. Rev. Lett. , 064101 (2007).[19] B. Kralemann, A. Pikovsky, and M. Rosenblum, “Re-constructing phase dynamics of oscillator networks,”Chaos , 025104 (2011); “Reconstructing effective phaseconnectivity of oscillator networks from observations,”New J. Phys. , 085013 (2014).[20] A. Pikovsky and M. Rosenblum, “Dynamics of globally cou-pled oscillators: Progress and perspectives,” Chaos , 097616(2015).[21] V. Hakim and W. J. Rappel, “Dynamics of the globally coupledcomplex Ginzburg-Landau equation.” Phys. Rev. A , R7347–R7350 (1992).[22] N. Nakagawa and Y. Kuramoto, “Collective chaos in a popu-lation of globally coupled oscillators,” Prog. Theor. Phys. ,313–323 (1993).[23] N. Nakagawa and Y. Kuramoto, “From collective oscillations tocollective chaos in a globally coupled oscillator system,” Phys-ica D , 74–80 (1994).[24] N. Nakagawa and Y. Kuramoto, “Anomalous lyapunov spec-trum in globally coupled oscillators,” Physica D , 307–316(1995).[25] M.-L. Chabanol, V. Hakim, and W.-J. Rappel, “Collectivechaos and noise in the globally coupled complex Ginzburg-Landau equation,” Physica D , 273 – 293 (1997).[26] M. Banaji and P. Glendinning, “Towards a quasi-periodicmean field theory for globally coupled oscillators,”Phys. Lett. A , 297 – 302 (1999).[27] M. Banaji, “Clustering in globally coupled oscillators,”Dynam. Syst. , 263–285 (2002).[28] H. Daido and K. Nakanishi, “Diffusion-induced inhomogeneityin globally coupled oscillators: Swing-by mechanism,” Phys.Rev. Lett. , 054101 (2006).[29] H. Daido and K. Nakanishi, “Aging and clustering in globallycoupled oscillators,” Phys. Rev. E , 056206 (2007).[30] K. A. Takeuchi, F. Ginelli, and H. Chat´e, “Lyapunov analy-sis captures the collective dynamics of large chaotic systems,”Phys. Rev. Lett. , 154103 (2009).[31] K. A. Takeuchi and H. Chat´e, “Collective lyapunov modes,”J. Phys. A: Math. Theor. , 254007 (2013).[32] G. C. Sethia and A. Sen, “Chimera states: The existence criteriarevisited,” Phys. Rev. Lett. , 144101 (2014).[33] Wai Lim Ku, Michelle Girvan, and Edward Ott, “Dynami-cal transitions in large systems of mean field-coupled Landau-Stuart oscillators: Extensive chaos and cluster states,” Chaos , 123122 (2015).[34] M. Rosenblum and A. Pikovsky, “Two types ofquasiperiodic partial synchrony in oscillator ensembles,”Phys. Rev. E , 012919 (2015).[35] F. P. Kemeth, S. W. Haugland, and K. Krischer, “Cluster singu-larity: The unfolding of clustering behavior in globally coupledStuart-Landau oscillators,” Chaos , 023107 (2019).[36] P. Clusella and A. Politi, “Between phase and amplitude oscil-lators,” Phys. Rev. E , 062201 (2019).[37] Y. Kuramoto, “Self-entrainment of a population of coupled non-linear oscillators,” in International Symposium on Mathemati-cal Problems in Theoretical Physics , Lecture Notes in Physics,Vol. 39, edited by Huzihiro Araki (Springer, Berlin, 1975) pp. 420–422.[38] D. Walgraef,
Spatio-Temporal Pattern Formation (Springer-Verlag, New York, 1997).[39] I. S. Aranson and L. Kramer, “The world of the complexGinzburg-Landau equation,” Rev. Mod. Phys. , 99 (2002).[40] H. Sakaguchi and Y. Kuramoto, “A soluble active rotatormodel showing phase transitions via mutual entrainment,”Prog. Theor. Phys. , 576–581 (1986).[41] E. Montbri´o and D. Paz´o, “Collective synchronization inthe presence of reactive coupling and shear diversity,”Phys. Rev. E , 046206 (2011); “Shear diversity prevents col-lective synchronization,” Phys. Rev. Lett. , 254101 (2011);D. Paz´o and E. Montbri´o, “The Kuramoto model with dis-tributed shear,” EPL (Europhys. Lett.) , 60007 (2011).[42] A. T. Winfree, “Patterns of phase compromise in biological cy-cles,” J. Math. Biol. , 73–93 (1974).[43] J. Guckenheimer, “Isochrons and phaseless sets,”J. Math. Biol. , 259–273 (1975).[44] E. A. Coddington and N. Levinson, Theory of Ordinary Differ-ential Equations (McGraw-Hill, New York, 1955) Chap. 13.[45] H. Daido, “Onset of cooperative entrainment in limit-cycle os-cillators with uniform all-to-all interactions: bifurcation of theorder function,” Physica D , 24–66 (1996).[46] C. van Vreeswijk, “Partial synchronization in populations ofpulse-coupled oscillators,” Phys. Rev. E , 5522–5537 (1996).[47] D. Golomb, D. Hansel, B. Shraiman, and H. Som-polinsky, “Clustering in globally coupled phase oscillators,”Phys. Rev. A , 3516–3530 (1992).[48] D. Hansel, G. Mato, and C. Meunier, “Clusteringand slow switching in globally coupled phase oscillators,”Phys. Rev. E , 3470–3477 (1993).[49] F. P. Kemeth, S. W. Haugland, and K. Krischer, “Symmetriesof chimera states,” Phys. Rev. Lett. , 214101 (2018).[50] P. Ashwin and J. W. Swift, “The dynamics ofn weakly coupledidentical oscillators,” J. Nonlin. Sci. , 69–108 (1992).[51] P. Ashwin, O. Burylko, and Y. Maistrenko, “Bifurcation to het-eroclinic cycles and sensitivity in three and four coupled phaseoscillators,” Physica D , 454 – 466 (2008).[52] Y. Kuramoto, “Phase- and center-manifold reduc-tions for large populations of coupled oscillatorswith application to non-locally coupled systems,”Int. J. Bifurcat. Chaos , 789–805 (1997).[53] A. Pikovsky and P. Rosenau, “Phase compactons,”Physica D , 56 – 69 (2006).[54] H. Kori, Y. Kuramoto, S. Jain, I. Z. Kiss, andJ. L. Hudson, “Clustering in globally coupled oscilla-tors near a hopf bifurcation: Theory and experiments,”Phys. Rev. E , 062906 (2014).[55] D. Paz´o and E. Montbri´o, “From quasiperiodic partial synchro-nization to collective chaos in populations of inhibitory neu-rons with delay,” Phys. Rev. Lett. , 238101 (2016); F. De-valle, E. Montbri´o, and D. Paz´o, “Dynamics of alarge system of spiking neurons with synaptic delay,”Phys. Rev. E , 042214 (2018).[56] Otherwise, the frequency of (N)UIS would depend nonlinearlyon ǫ , in disagreement with the MF-CGLE.[57] C. Bick, M. Timme, D. Paulikat, D. Rathlev, andP. Ashwin, “Chaos in symmetric phase oscillator networks,”Phys. Rev. Lett. , 244101 (2011).[58] C. Bick, P. Ashwin, and A. Rodrigues, “Chaos in genericallycoupled phase oscillator networks with nonpairwise interac-tions,” Chaos , 094814 (2016).[59] P. Clusella, A. Politi, and M. Rosenblum, “A min-imal model of self-consistent partial synchrony,” New J. Phys. , 093037 (2016).[60] C. Bick, M. J. Panaggio, and E. A. Martens, “Chaos in Ku-ramoto oscillator networks,” Chaos , 071102 (2018).[61] C. Bick, “Heteroclinic switching between chimeras,”Phys. Rev. E , 050201(R) (2018).[62] P. S. Skardal and A. Arenas, “Abrupt desynchronization and ex-tensive multistability in globally coupled oscillator simplices,” Phys. Rev. Lett. , 248301 (2019).[63] M. Rosenblum and A. Pikovsky, “Numerical phasereduction beyond the first order approximation,”Chaos , 011105 (2019).[64] H. Kori and Y. Kuramoto, “Slow switching in globally coupledoscillators: robustness and occurrence through delayed cou-pling,” Phys. Rev. E63