Phase reduction beyond the first order: the case of the mean-field complex Ginzburg-Landau equation

Phase reduction is a powerful technique that makes possible describe the dynamics of a weakly perturbed limit-cycle oscillator in terms of its phase. For ensembles of oscillators, a classical example of phase reduction is 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 incoherence, but none of the nontrivial dynamics of the MF-CGLE. This fact calls for an expansion beyond the first order in the coupling constant. We develop an isochron-based scheme to obtain the second-order phase approximation, which reproduces the weak coupling dynamics of the MF-CGLE. The practicality of our method is evidenced by extending the calculation up to third order. Each new term of the power series expansion contributes with additional higher-order multi-body (i.e.non-pairwise) interactions. This points to intricate multi-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 as neuroscience or engineering [1][2][3][4].It is an empirical fact that some phenomena arising in these systems can be understood in terms of interacting phase oscillators.This framework has proven to be useful modeling and engineering experimental setups composed of many rhythmic elements, operating in a wide range of spatio-temporal scales, and interacting through very different physical processes.We may cite small motors-cell phone vibrators-interacting through an elastic plate [5], networks of (electro-)chemical oscillators [6,7], arrays of Josephson junctions [8,9] and globally coupled electrical self-oscillators [10,11], or nanoelectromechanical oscillators in a ring [12].
Applying a phase reduction method [1,[13][14][15] is the rigorous way of describing a weakly perturbed oscillator solely in terms of its phase (the other degrees of freedom become enslaved).However, obtaining analytically the approximate 'phase-only model' for a specific system is not an easy task.Moreover, phase reduction becomes inaccurate unless the disturbances are not sufficiently weak.While, according to common wisdom phase reduction of oscillator ensembles yields pairwise interacting phase oscillators [13], multi-body (i.e.non-pairwise) interactions may also be relevant in some contexts.Apart from the idea of invoking hypothetical three-body interacting limit-cycle oscillators [16], multi-body phase interactions naturally arise if the coupling is nonlinear [17], see also [18].Instead, for linear pairwise coupling, three-body interactions are a distinctive element of second-order phase approximations, as recently highlighted in [12].Recognizing the ubiquity of multi-body interactions may also be important for reconstructing phase interactions from data [19].
Much of our knowledge on nonlinear dynamics relies on minimal models that capture the essential mechanisms behind complex phenomena.For oscillatory dynamics, the conventional test bed is the normal form of the Hopf bifurcation above criticality: the so-called Stuart-Landau oscillator.Concerning geometry, global coupling is a fruitful simplifying as-sumption [1,13,20].These two ingredients are combined in a standard model of collective dynamics: the fully connected network of Stuart-Landau oscillators, or the mean-field version of the complex Ginzburg-Landau equation (MF-CGLE) [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36].This system is particularly interesting for chaos theory since it exhibits both microscopic (extensive) and macroscopic (collective) chaos, either combined or independently, depending on parameters [21-24, 26, 30, 33, 36].Phase reduction of the MF-CGLE yields the Kuramoto-Sakaguchi model [13,22,37], a first-order approximation that behaves in a pathological way (unless heterogeneities are present): it only displays full synchrony or incoherence.Therefore, pure collective chaos and other phase dynamics of the MF-CGLE remain to be analytically described in terms of a phase model.Such a phase reduction should provide additional insights into the nature of collective chaos (playing an analogous role to the Kuramoto-Sivashinsky equation of phase turbulence).
The aim of this paper is two-fold: we introduce a phasereduction method, and we investigate the phase model obtained from the MF-CGLE.The paper is organized as follows.In Sec.II we reexamine the phase dynamics of the MF-CGLE and the connection with its first-order phase reduction (the Kuramoto model).In section III, we present our systematic phase-reduction procedure, based on the direct use of isochrons, which delivers a well-controlled power expansion in the coupling strength parameter.Section IV is devoted to investigating the weak coupling limit of the MF-CGLE by means of the the second-order phase reduction, which unfolds the degeneracies of the Kuramoto model; we address the cases of a large ensemble of oscillators, as well as an small one of four oscillators.Section V presents the third-order contribution to the phase reduction of the MF-CGLE.Finally, in Sec.VI we discuss the implications of our work and some outlooks.ordinary differential equations: Here, A j = r j e iϕj is a complex variable (index j runs from 1 to N ), and the mean field is Ā = N −1 N k=1 A k .Apart from the population size N , there are three free parameters in Eq. ( 1): ǫ, c 1 , and c 2 .Parameter ǫ, controlling the coupling strength, is positive in order to preserve the analogy with the (spatially extended) Ginzburg-Landau equation.Parameter c 1 introduces a cross-coupling between real and imaginary parts of the A j 's.This non-dissipative coupling, so-called 'reactive' [4], generically appears from center manifold reduction [13].Finally, 'nonisochronicity' (or 'shear') parameter c 2 in Eq. ( 1) determines the dependence of the angular velocity of one oscillator on its radial coordinate.There are two important symmetries in system (1): invariance under a global phase shift A j → A j e iφ , and full permutation symmetry stemming from the mean-field coupling.

A. Phenomenology
For many parameter values, the global attractor of Eq. ( 1) is either full synchronization (FS) A j = Ā = e −ic2t or one incoherent state with vanishing mean field Ā = 0.In the latter case the oscillators rotate freely, ]t + φ j }.Among all the states compatible with Ā = 0, the most prominent one is the uniform incoherent state (UIS) in which the φ j are uniformly distributed in the thermodynamic limit (for a finite ensemble, the φ j are evenly spaced, deserving the name of splay state or ponies on a merry-go-round state).A continuum of nonuniform incoherent states (NUISs) coexists with UIS, but usually arbitrarily weak noise spreads the phases and UIS is eventually attained.Nonetheless, for certain parameters values, such as those in Fig. 1(a), the UIS is unstable and one NUIS sets in spontaneously [21,25].
In addition to FS, UIS, and NUIS, system (1) exhibits a rich repertoire 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 mean field Ā rotates uniformly, while the individual oscillators behave quasiperiodically (since each oscillators 'feels' the periodic driving of the mean field).Remarkably, increasing coupling QPS may undergo a couple of secondary Hopf bifurcations resulting in a state of pure collective chaos [24,36].With this term we refer to a state in which the mean field behaves chaotically, while individual oscillators behave in seemingly chaotic-like fashion (neighboring oscillators remain close for ever due to the absence of microscopic chaos).A shared feature of NUIS, QPS and pure collective chaos [22,24,36] is that the relative positions of the oscillators on top of a closed curve are preserved, see Figs. 1(a), 1(b) and 1(c).This fact suggests that a description in terms of oscillators' phases alone is possible.In contrast to Fig. 1(c), Fig. 1(d) shows a chaotic regime in which phase description breaks down, as it involves microscopic degrees of freedom and no phase ordering exists.Hence, our ultimate goal is to find a phase-reduced model of Eq. (1) that captures as much as possible of the phasedescribable states (NUIS, QPS, modulated QPS, pure collective chaos,etc.).

B. Basic phase diagrams
Before presenting our results it is convenient to review previous results on the MF-CGLE.For fixed c 1 and c 2 values, let us denote by ǫ s and ǫ 0 , the ǫ values of marginal linear stability for FS and UIS.Closed formulas for ǫ s and ǫ 0 are [21,22]: These formulas are also valid for finite ensembles, assuming ǫ 0 refers to the splay state.To visualize the stability boundaries in Eqs. ( 2) and (3), it is convenient to fix either c 1 or c 2 .
Following [22] we choose to fix c 2 , and display the loci of ǫ s and ǫ 0 in the parameter plane (c 1 , ǫ).In the phase diagrams in Figs.2(a) 2(b) and 2(c) we selected c 2 = 3, c 2 = 2 and c 2 = 1, respectively.This choice is motivated by the fact that most previous works on the MF-CGLE adopt either c 2 = 2 or c 2 = 3.One key observation is that, as ǫ s and ǫ 0 approach zero, the boundaries converge to the condition 1 + c 1 c 2 = 0, which is the well known Benjamin-Feir-Newell criterion for the stability of uniform oscillations in the complex Ginzburg-Landau equation in arbitrary dimension [4,13,38,39].There is a critical value c 2 = √ 3 = 1.732 . . .at which the boundaries ǫ s and ǫ 0 become tangent at ǫ = 0. Accordingly, for c 2 = 1, see Fig. 2(c), there is a region of bistability between UIS and synchrony, in contrast e.g. to Fig. 2(b).
The stability of a NUIS depends exclusively on the mean field Q = |N −1 j exp(2iϕ j )|.The coupling constant ǫ Q at which one particular NUIS becomes unstable was obtained in Ref. [25]: This formula is the generalization of (3) with the important qualitative information that the size of the stability region increases as Q grows, reaching its maximum for Q = 1.At Q = 1 the NUIS collapses into a two-cluster state with equally populated groups.The value of Q is still far from breaking the degeneracy of a NUIS, provided Q = 1, since the values of all 'higher-order' mean fields f n = |N −1 j exp(niϕ j )| (n > 2) are free.Nevertheless, the conclusion based on numerical simulation is that any small amount of noise causes f n to converge to zero, and Q to take the smallest value among all non-unstable (i.e.neutrally stable) NUISs.Therefore, it is assumed hereafter that the term NUIS is constrained to Figures 2 (a) and (b) include a green-hatched region, adjacent to the UIS region at moderate ǫ values, where other phase-describable states are stable.These are QPS, modulated QPS and pure-collective chaos [24,36].We determined the 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 technique [4,13,37] to Eq. ( 1) yields the Kuramoto-Sakaguchi model [40].In this model, each oscillator is described by a phase θ j , and it is coupled to the other ones by pairwise interactions of the form sin(θ i − θ j + α).In agreement with the mean-field character of the system, oscillators are coupled through the Kuramoto order parameter Z 1 ≡ R e iΨ = N −1 N k=1 e iθ k , such that the ordinary differential equations governing the dynamics are: with constants ), and phase lag Equation ( 5) is the disorder-free version of the paradigmatic Kuramoto-Sakaguchi model [13,40] and related models [41].The dynamics of Eq. ( 5) is determined by the sign of 1 + c 1 c 2 (Benjamin-Feir-Newell criterion): full synchrony -corresponding to R = 1is stable for 1 + c 1 c 2 > 0, and unstable for 1 + c 1 c 2 < 0. In the latter case, among infinitely many oscillator densities with R = 0, there is a convergence to the UIS under an arbitrarily weak noise [22].
As discussed above, the MF-CGLE has much richer dynamics than its first-order phase reduction (5), even arbitrarily close to the ǫ = 0 limit.Therefore, it is mandatory to extend the phase reduction to order O(ǫ 2 ) if we wish to avoid degeneracies in the phase approximation.This is what we do next.

III. SYSTEMATIC PHASE REDUCTION
In spite of the relevance of Eq. ( 1) no phase reduction beyond the first order is currently available.Finding higher order terms in the phase reduction is necessary to unfold the singularity at (c 1 , ǫ) = (−1/c 2 , 0), see Fig. 2.This path of investigation should allow us to discern which are the true behaviors of the MF-CGLE in the small coupling limit |ǫ| ≪ 1.Moreover, it might serve to shed light on the mechanisms behind complex dynamics found (so far) for moderate ǫ values.
An isochron-based phase reduction approach is developed here.It allowed us to obtain the phase reduction of the MF-CGLE up to order ǫ 3 .In this section we give the details of our phase reduction calculation.We anticipate that the results at second and third order in ǫ correspond to Eqs. ( 15) and ( 29) below.

A. Isochrons
The concept of isochron [42,43] is the cornerstone of phase reduction methods [1,13].Isochrons foliate the attraction basin of a stable limit cycle, each intersecting it at one point.The phase of that point is attributed to all points of the isochron, motivated by their convergence as time goes to infinity (the so-called 'asymptotic phase' [44]).For the Stuart-Landau oscillator, polar coordinates (r, ϕ) relate to the phase θ according to [4,13]: As mentioned above, on the limit cycle (r = 1), θ = ϕ.
The term "nonisochronicity" or "shear" for parameter c 2 becomes clear in light of Eq. ( 7), since c 2 controls how much the isochrons deviate from radial lines.

B. Isochron-based phase reduction
We continue the analysis writing Eq. ( 1) in polar coordinates: After the change of variables (r j , ϕ j ) → (r j , θ j ) through Eq. ( 7), we get: Here, we have also implemented the transformation θ j → θ j − c 2 t (by moving to a rotating frame with angular velocity −c 2 ).In this way, the time derivatives of the phases in (10b) are proportional to ǫ, while the r j are fast variables that become enslaved to the dynamics of θ j .In Eq. ( 10) , and functions g j and h j depend on the vectors r = (r 1 , r 2 , . . ., r N ) T and θ = (θ 1 , θ 2 , . . ., θ N ) T as follows, The separation of time scales in Eq. ( 10) suggests using classical perturbation techniques like averaging, adiabatic approximation, or two-timing.However, the perturbation approach described next proved to be both conceptually simple and much less convoluted, permitting us to obtain the phase reduction up to cubic order in ǫ.Based on the empirical observation that, at small ǫ values, the oscillators fall on a closed curve and preserve their phase ordering, we assume that the radii are completely determined by the phases r j = r j (θ 1 , θ 2 , . . ., θ N ).We also postulate an expansion in powers of ǫ for the radii: r = r (0) l .Now, the explicit dependence on the radii in (12) must be removed.This is accomplished equating both sides of (10a) at the same order.The order O(ǫ 0 ) yields r (0) j = 1, and ( 12) becomes (at the lowest order) the Kuramoto-Sakaguchi model (5).At order ǫ: As r j depends exclusively on the phases, we can apply the chain rule: ṙj = (∇ θ r j ) • θ.At order ǫ, the time derivative vanishes: Hence Eq. ( 13) yields the result which can be inserted in (12) to obtain the ǫ 2 contribution.Through elementary manipulations the second-order phase reduction of Eq. ( 1) can be condensed into this expression: The O(ǫ 2 ) term depends on Z 1 as well as on the second Kuramoto-Daido order parameter [45] To enhance the clarity of Eq. ( 15), we found it convenient to define a phase lag which turns out to be independent of c 2 .The other constants in Eq. ( 15) are the same as in Eq. ( 5); as the change to a rotating frame has been reversed, the O(ǫ 0 ) term inside Ω is −c 2 (as before).

IV. SECOND-ORDER PHASE REDUCTION: THREE-BODY INTERACTIONS
In this section we study in detail the phase model obtained from the second-order phase reduction of the MF-CGLE, i.e. the system of phase oscillators governed by Eq. (15).Of the three O(ǫ 2 ) contributions to Eq. ( 15), the first element of the sum (m = 1) entails a parameter shift to the O(ǫ) interaction, and it is therefore irrelevant in qualitative terms.The other two terms in Eq. ( 15) correspond to three-body (i.e.nonpairwise) interactions: The price of working only with the phases is that two-body interactions of the original MF-CGLE (1) become multi-body interactions, as higher orders of ǫ are considered.In comparison to Eq. (1) our phase model can be much more efficiently analyzed, both analytically and numerically.We devote the remainder of this section to analyze the phase model in Eq. ( 15).We note that, as expected, the model is invariant under global phase shift θ j → θ j + φ.For the sake of making the presentation simpler we assume constant Ω = 0, since this can always be achieved by going to a rotating frame θ j → θ j + Ωt.
For illustration, the curve defined by ( 19) is represented by a blue dotted line in Figs.3(a) and (b) for c 2 = 3 and c 1 = 1, respectively.Equation ( 19) is asymptotically exact as ǫ s → 0, 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 phases between θ and θ + dθ.Now θ ∈ [0, 2π) is a cyclic variable, i.e. ρ(θ + 2π, t) = ρ(θ, t), and we impose the normalization condition 2π 0 ρ(θ, t)dθ = 1.The oscillator density ρ obeys the continuity equation because of the conservation of the number of oscillators: Here v = θ is the ρ-dependent velocity of an oscillator with phase θ.We define the Fourier modes of ρ: with ρ 0 = 1 and ρ n = ρ * −n .The mean fields Z n reduce to Inserting the Fourier expansion (21) into the continuity equation (20) allows us to rewrite our model in Fourier space: 1. Uniform incoherent state The stability boundary of the UIS (ρ(θ) = (2π) −1 ⇔ ρ n =0 = 0) is obtained linearizing the previous equation.It is easy to notice that only the first mode may destabilize.We have for |ρ 1 | ≪ 1: Neglecting the trivial marginal case ǫ = 0, the stability boundary satisfies cos α+(1/4)ǫ 0 η cos β = 0. Or in terms of c 1 and c 2 : In Figs.3(a) and 3(b), we can contrast this formula with the exact one for the MF-CGLE, Eq. ( 3), for two c 2 values.

Nonuniform incoherent states
According to (22), in an incoherent state (ρ 1 = 0) higherorder modes are at rest: ρn = 0 (n > 2).The linearization of ( 22) around ρ 1 = 0 and ρ n = 0 (|n| ≥ 2) is (schematically) as follows: where the • symbols denote nonzero elements.Clearly, the structure of this equation yields an infinite set of vanishing eigenvalues plus two eigenvalues coming from the first two rows.The equation for δρ 1 is hence, the only relevant one.The linear terms in δρ 1 yield: All higher-order modes, save ρ 2 , are absent in the equation.As ρ2 = 0, we can choose the coordinate axes such that ρ 2 = Q ∈ R.After some calculations we find that NUIS with a specific Q value is marginally stable at: As occurs in the MF-CGLE the larger Q is, the larger is the stability region of the NUIS.Our empirical observation is that, for given c 1 and c 2 , if ǫ is set at a certain ǫ = ǫ Q * the numerical integration of the system (either oscillators or Fourier modes), under a very weak noise, always converges to a NUIS with ρ n≥3 = 0; and, |ρ 2 | = Q * .In other words, the system adopts the minimum value of |ρ 2 | among all allowed by Eq. ( 27).
The state Q = 1 (R = 0) -the last NUIS to destabilizeis singular, not only because it is just a two-cluster state with two equally populated groups, but also because in contrast to the other NUIS, the instability is not oscillatory.Eq. ( 26) takes the form δρ 1 ∝ a δρ 1 − a * δρ * 1 what yields an additional zero eigenvalue corresponding to the direction Im(δρ 1 ) = 0.

C. Validity and accuracy
From our previous results, we conclude that the phase reduction (15) is free of degeneracies.The boundaries of FS, UIS and NUIS with different Q values do not overlap.As a double-check of the correctness of our analysis, we verified that the boundaries ( 19), (24), and ( 27) obtained through the phase reduction are tangent to the equivalent boundaries of the MF-CGLE, Eqs. ( 2), ( 3) and ( 4), at ǫ = 0.In Fig. 3 we depict together the boundaries of FS, UIS, and (Q = 1)-NUIS of the MF-CGLE (solid lines) and phase reduction to second order (dashed lines) for two values of the nonisochronicity: c 2 = 3 and c 2 = 1.These plots permit us to identify the range of ǫ in which the second-order approximation is accurate.For c 2 = 1 the approximate bifurcation lines are accurate up to ǫ ≈ 0.05, while this range is certainly smaller for c 2 = 3.
For general c 1 , c 2 values, the prefactor (ǫη) n appearing for first (n = 1) and second (n = 2) orders suggests to extrapolate the relative smallness of ǫ to other c 1 , and c 2 values.Thus, if in Fig. 3(b) accuracy is good up to ǫη ≈ 0.05η, and η ≈ 2, we propose as a conservative range of validity of the second-order approximation.Nevertheless, Eq. ( 28) must be regarded with some caution, since the third-order contribution to the phase reduction expansion is not exactly proportional to (ǫη) 3 , see Sec.V.

D. Quasiperiodic partial synchronization
The phenomenon of QPS was originally reported in the MF-CGLE [24] as a state emerging from the destabilization of the UIS, see Fig. 1(b), though its finding is usually attributed to a model of phase oscillators [46].As mentioned above, in a QPS state the mean-field rotates uniformly, but individual oscillators behave quasiperiodically.Each oscillator passes periodically through a bottleneck located at the phase arg( Ā).The onset of QPS looks like a Hopf bifurcation undergone by the UIS, but this is not the case because of the infinitely many neutral directions pointing to nearby NUISs.It is also important to emphasize that stable QPS does not settles any time that the UIS becomes unstable.As can be appreciated in Fig. 2(a) and (b), QPS is only observed at moderate ǫ values when entering inside the green hatched region.Otherwise, what we observe in the MF-CGLE is that the QPS state born at the instability of the UIS is a saddle.For parameter values with unstable UIS and FS -outside the green-hatched regions-initial conditions close to the UIS approach QPS for long time, eventually converging to one NUIS.If any small amount of noise is present, the NUIS with the lowest Q among the nonunstable ones is selected.The same behavior is displayed by the second-order phase reduction, Eq. ( 15); see Fig. 4. The logarithmic scaling of the residence times near QPS indicates a heteroclinic connection between UIS and QPS.The amplitude of the saddle QPS depends on the particular parameter values.The state of QPS progressively grows as we move away from the UIS stability boundary, finally colliding with FS (|ρ n | → 1) at the point where FS becomes stable.
All in all, these results confirm the correctness of our expansion, but at the same time prove the limitations of the secondorder reduction, since the QPS attractor-found at moderate ǫ values, see Fig. 2(a) and (b)-is not reproduced.

E. Clustering
Clustering is a much studied phenomenon in oscillator ensembles [47].In a clustered state there are several groups of oscillators, each group formed by oscillators sharing the same phase.This kind of states are always possible in a mean-field model, so the relevant question is the stability.Indeed, the MF-CGLE is known to exhibit stable cluster at certain parameter ranges [21, 22, 27-29, 33, 35], specifically for strong coupling (ǫ ≈ 1).
Are there stable clustered states at small coupling?Our phase model allows us to address this question in an analytical way.Nonetheless, the general problem is intractable and we decided to restrict our study to states with two point clusters, where a fraction p of the population is in the A-state θ A , and the remaining (1 − p) fraction is in the B-state θ B = θ A .We now summarize the results; the corresponding calculations can be found in Appendix I.
As an illustrative example, Fig. 5 depicts the combinations of phase difference ∆ = θ A − θ B and imbalance p corresponding to actual cluster solutions for three different c 1 values with fixed values of c 2 and ǫ.Each panel is a typical situation in a specific region of parameter space.At the FS threshold, between panels (a) and (b), there is an infinity of twocluster solutions colliding with FS (∆ = 0).In consequence there is a reconnection of the two-cluster solutions.In Fig.  solid lines represent stable locking of the clusters.However, these solutions are fragile against disintegration of the largest cluster.Our conclusion after an extensive exploration of parameter space is that stable two-cluster states are not stable at small coupling.To be more what we observe in our second-order phase reduction, Eq. ( 15), is that stable clustering is hardly found, and if so, it always requires moderate coupling strengths, violating (28).And indeed, we could not replicate clustering in the MF-CGLE for the parameter values predicted by Eq. (15).
The stability analysis of the two-cluster solutions also confirmed that slow switching [48] -a stable heteroclinic connection between two configurations of ∆ = θ A − θ B -is not possible.

F. Finite population, N = 4
This work focuses on the behavior of the MF-CGLE in the thermodynamic limit (N → ∞).But the phase reduction (15) is valid for an arbitrary population size.In this section we construct a bifurcation diagram for N = 4 oscillators, one size previously considered in the MF-CGLE context [35,49].Here, this choice is motivated by the fact that in globally coupled systems this is the smallest size with a continuum of states with R = 0 [50], equivalent to the NUISs for N = ∞.
In analogy with its thermodynamic limit, the finite-N the Kuramoto-Sakaguchi model has an exceptional transition between FS and the splay state at 1 + c 1 c 2 = 0.This degeneracy can be broken down, for instance, adding higher-order harmonics to the (pairwise) interactions [51].In our case, degeneracy is broken down by the three-body interactions of the second-order in the phase-reduction expansion.
Working with a small number of oscillators has the advantage that we can track all the stationary solutions, in particular the clustered solutions.As there are 3! orderings for the oscillators phases, and phase ordering is preserved by the dynamics because of the mean-field interactions, we choose the oscillators' labels such that θ j ≤ θ j+1 .(We assume here θ j ∈ [0, 2π) to avoid artificial degeneracies.)The set of phases {θ 1 , θ 2 , θ 3 , θ 4 }, may take several invariant configurations.Apart from the trivial FS state {a, a, a, a}, there exists a continuum of "NUIS-like" Z 2 -symmetric states with {a, a + b, a + π, a + b + π}, where b ∈ (0, π/2) ∪ (π/2, π).In the limit b → π/2 the NUIS becomes the splay state (the analogous of UIS).In addition, in the limits b → 0 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 symmetryrelated 2-2 configurations {a, a, b, b} (b = a + π).Additionally, 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 not exist in our phase reduction, in contrast to the MF-CGLE for strong coupling [35].Concerning the N = 4 analogous of QPS, it is a periodic orbit, in which, due to the finiteness of the population, R and Q fluctuate around their average values.
We use R, Q, and c 1 to plot the bifurcation diagram in Fig. 6.These coordinates have the drawback of collapsing multiple equivalent states to a single point, hiding symmetries (e.g.pitchfork bifurcations).However, our choice intends to ease the comparison with the previous section, and with the same aim states are labeled borrowing the infinite-N terminology; namely, we use the labels UIS, NUIS, and QPS instead of splay state, Z 2 -symmetric state, and limit cycle, respectively.
Due to permutation symmetry FS destabilizes at point T in Fig. 6, as three eigenvalues go through zero simultaneously.This comprises an equivariant transcritical bifurcation with the 3-1 cluster, as well as a pitchfork bifurcation involving a 2-2 cluster.Moreover, at point T, QPS collapses into a heteroclinic cycle.This coincidence of bifurcations is a known scenario in systems with full permutation symmetry [50].Concerning UIS, it undergoes an oscillatory instability at point U, but this is not a standard Hopf bifurcation because of the neutral direction along the NUIS manifold.QPS is a saddle, and not a stable limit cycle as it might have been naively expected.In Fig. 6 we took c 2 = 3, and the QPS branch connects points T and U in a simple way.In contrast to Fig. 6, for c 2 = 1 FS and UIS coexist, and points U and T switch their relative positions.In that case the QPS branch is completely reversed (not shown), and the QPS solution is fully unstable.Consistently, we found a range of c 2 values in between, 1 < c 2 < 3, where (depending on ǫ) the QPS branch develops a fold.
In summary, the bifurcation diagram for N = 4 appears to capture the global picture of the transition from UIS to FS.  15) with N = 4 oscillators and c2 = 3, ǫ = 0.1.Solid (dashed) lines represent stable (unstable) solutions.In the case of UIS and NUIS the solution depicted must be understood as the one observed under arbitrarily weak noise (there is continuum of neutral solutions with Q larger than the solution depicted).The saddle QPS orbit was continued by means of a Newton-Raphson algorithm, and the values of R and Q assigned in the diagram correspond to their time averages.
Considering more oscillators will increase the number of cluster solutions, see [27], but no essential new features.

V. THIRD-ORDER PHASE REDUCTION: FOUR-BODY INTERACTIONS
Our reason to deal with the third-order term now is to illustrate the practicality of the phase reduction method, and get a glimpse of the power series expansion at higher orders.Evaluating the cubic term in Eq. ( 12) yields the O(ǫ 3 ) correction to Eq. ( 15): This expression depends on the third Kuramoto-Daido order parameter Z 3 ≡ P e iΞ = N −1 j e i3θj .The dependence of constants {C j , γ j } j=1,...,9 and D on c 1 and c 2 is tabulated in Appendix II.The structure of Eq. ( 29) deserves some words here.The terms proportional to C j with indices j = 1, 2, 3 are higher-order corrections to Eq. ( 15), tantamount to a shift in parameter values.Four-body interactions appear in five different forms, corresponding to indices j = 4, . . ., 8. For illustration we expand a couple of these four-body contributions: There are several qualitative features in Eq. ( 29) that deserve to be pointed out: 1.The overall O(ǫ 3 ) contribution is not proportional to η 3 -though some terms indeed are-in contrast to O(ǫ) and O(ǫ 2 ), which are proportional to η and η 2 , respectively.
2. From Eqs. ( 15) and ( 29) we can expect that truncation of the power series to order ǫ n yields up to (n + 1)body interactions, but not higher-order non-pairwise couplings.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 unexpected, (nonetheless see [17]), since they depend on the mean fields Z 1 and Z 2 , but not on θ j itself.They are hence irrelevant concerning synchronization boundaries.
4. As occurs with the O(ǫ 2 ) term, FS and (N)UIS states are consistent with the MF-CGLE dynamics: (i) all terms in (29) are proportional to R ensuring that the contribution to the oscillators' frequencies vanishes in one incoherent state; (ii) in the FS state, the contribution also vanishes, as expected since the frequency of FS in the MF-CGLE varies linearly with ǫ.Accordingly, it holds that D + j C j sin γ j = 0, cf.Appendix II.
Unfortunately, there is not a recognizable pattern in the new terms appearing in the power series expansion, so it is not possible to extrapolate to higher orders in ǫ.
From Eqs.( 15) and ( 29) we can derive the stability boundary of FS, NUIS (for UIS just set Q = 0) obtaining: In Fig. 7 we depict (a) ǫ 0 , and (b) ǫ s from the previous expressions and compare them with the result of the MF-CGLE, and with the second-order approximation.The slopes and the curvatures of the bifurcation lines of the third-order phase reduction agree with those of the MF-CGLE at ǫ = 0.

VI. DISCUSSION
A. Alternative phase reduction(s) Our phase reduction is a genuine power series in the small parameter ǫ.Another strategy to analyze (1) is to absorb the ǫA j term prior to the phase reduction.Specifically, setting t ′ = (1 − ǫ)t, and we get where B j = A j exp(iǫc 1 t)/ √ 1 − ǫ.Applying our phasereduction method to (33) we obtain an alternative phase reduction in powers of κ (the result is not qualitatively different).
Is it worth transforming (1) into (33)?In other words, is the phase reduction of (33) up to order κ n , superior to that of (1) at order ǫ n ?Certainly, phase reductions at order ǫ n and κ n are not equivalent since κ = ǫ + ǫ 2 + ǫ 3 + • • • .Any truncation at order κ n involves all powers of ǫ.The relative accuracy of the phase reductions of ( 33) and (1) at the same order can be assessed comparing the bifurcation loci.Instead of applying phase reduction to Eq. ( 33), the quickest strategy is to assume the existence of an exact phase reduction involving all orders  1), while dotted and dashed lines correspond to second-and third-order phase approximations, respectively.Blue lines are obtained from ( 15) and (29).Orange lines are the results if prior to phase reduction the MF-CGLE is transformed into (33), performing an isochron-base phase reduction in powers of κ.
in ǫ such that the exact critical value ǫ * (the asterisk denotes an arbitrary state: UIS, FS, ...) satisfies: The coefficients a n depend on the specific instability we are considering.
A comparison of the bifurcations lines of UIS and FS is displayed in Figs.7(a) and 7(b) for c 2 = 3.We see that the transformation of (1) into (33) allows us to obtain a phase model that captures better the stability boundary of UIS, but not of FS.It is easy to understand why.Each strategy captures better the dynamics in which the quantities multiplying the coupling constant are small.Thus, Eq. ( 1) is already a good starting point for states close to FS (A j ≈ Ā), while (33) works better close to incoherence ( Ā ≈ B ≈ 0).Finally, note that in addition to (1) and ( 33), there exists a continuum of alternative, intermediate formulations, in which ǫA j is only partly absorbed by a coordinate transformation.

B. Possible extensions of this work
The phase-reduction procedure presented in this work can be easily implemented in other geometries, different from the fully connected network.In a networked architecture, phase reduction at first order in ǫ couples phases with the nearestneighbor phases.At order ǫ 2 , second nearest neighbors come also into play [19], and progressively more distant nodes participate in the phase dynamics at higher orders.Also, the case of non-locally coupled Stuart-Landau oscillators [52] is analyzable with the phase reduction presented here.Concerning the original complex Ginzburg-Landau equation, a partial differential equation of reaction-diffusion type, our phase reduction procedure is very simple and efficient obtaining the coefficients of the second-order terms: Concerning the oscillator dynamics, the phase reduction carried out here can be easily applied to planar oscillators with polar symmetry (λ − ω systems).In the latter case, analogously to (7), the isochrons satisfy θ = ϕ + χ(r) [1].Even if function χ does not have a closed form, it is still possible to obtain the phase model using the derivatives of the isochrons on the limit cycle.

C. Relationship with other phase-reduction approaches
In this subsection, we comment on the progress of our phase-reduction approach with respect to previous works, even if only directly applicable to λ − ω systems.
An alternative way of obtaining the second-order phase reduction of the MF-CGLE, Eq. ( 15), is applying the systematic averaging formulation in Chap. 4 of the book by Kuramoto [13].This calculation is, however, much more lengthy than the one presented in Sec.III.Not surprisingly, obtaining the order ǫ 3 with the averaging approach [13] is a totally impractical task, while we succeeded with our method (with the assistance of symbolic software); see Eq. (29).
Equation (15) can also be obtained assuming small variations of the radii, i.e. setting ṙj = 0.This procedure was followed in [12,53], with the difference being that there the small quantities are deviations from the reference limit cycle.Here, we pursued a bona-fide power expansion in terms of the coupling constant ǫ, and the result differs from the one obtained following [12,53].In passing, we mention that instead of assuming r j = r j (θ 1 , θ 2 , . . ., θ N ), once Eqs. ( 10) are derived, the two-timing approximation, such that the θ j depend only on a slow time τ = ǫt, can also be applied.
In contrast to our work, Refs.[17,54] apply first-order phase reduction obtaining multi-body phase interactions.The reason is that those works invoke amplitude equations for an ensemble close (but not asymptotically close) to a Hopf bifurcation.The amplitude equation, which can be seen as a generalization of Eq. ( 1), turns out to contain nonlinear interactions.The nonlinear coupling among the A j 's leads to multi-body interactions in the first-order phase reduction.Applying second-order phase reduction, as described here, to the amplitude 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 models with heterogeneity [60] or delay [55].Collective chaos in the MF-CGLE, see Fig. 1(c), calls for a phase description in terms of identical phase oscillators (without delays).The fact that we have not found evidence of collective chaos in our numerical simulations of the second-order phase reduction (15) -nor in the third-order one-can be reasonably attributed to a too restrictive truncation of the power expansion.We believe that a higher-order truncation will capture better the behavior of the system at larger ǫ values, and eventually, will exhibit collective chaos.
As pairwise interactions through higher harmonics, like Q sin(Φ − 2θ j ) = N −1 k sin[2(θ k − θ j )], do not show up in the phase reduction of the MF-CGLE [56], multi-body phase interactions appear to be the most promising ingredient to model collective chaos.In small ensembles of identical phase oscillators, higher harmonics as well as multi-body interactions promote chaos alike, see [57] and [58], respectively.However, so far, collective chaos remains elusive in populations of higher-order pairwise interacting identical phase oscillators [59].We believe multi-body interactions could be the key element of collective chaos, instead.
In the MF-CGLE with parameter values close to those in Fig. 1(c), we found chaos with a population sizes as small as N = 6.Does this say something about the order of the multibody interactions needed in the phase reduction?Is this chaos connected to collective chaos in the thermodynamic limit, as in [60]?

E. Conclusions
Multi-body interactions are an unavoidable consequence of phase reduction, but save for a few works [12,16,58,61,62], the role of multi-body phase interactions shaping exotic dynamics remains largely unexplored.In the weak-coupling regime of the MF-CGLE, multi-body phase interactions are essential to describe all states apart from FS and UIS.
In summary, in this work we achieve second-and thirdorder phase reductions of the MF-CGLE.In our view, higherorder phase reductions promise to be crucial for our understanding of collective chaos and other exotic phenomena [12].Moreover, analytic higher-order phase reductions may also serve as test beds for numerical phase reductions recently implemented [63].For these reasons, we regard phase reduction beyond the first order as an exciting battleground of nonlinear dynamics.

FIG. 2 .
FIG. 2. Partial phase diagram of the MF-CGLE for c2 = 3 (a), 2 (b), and 1 (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 becomes darker as Q → 1. Stable FS is indicated by a blue hatched region.The stability boundaries of FS, UIS and (Q = 1)-NUIS are depicted by blue, black and red lines, respectively; following Eqs.(2), (3), and (4) (setting Q = 1).In panels (a) and (b), there is green-hatched region where other phase-describable states like the ones shown in Figs.1(b) and 1(c) are stable.

FIG. 4 .
FIG.4.Evolution of R(t) for N = 1000 phase oscillators governed by(15) initiated near the UIS state.For the selected parameter values (c2 = 3, c1 = −0.38,ǫ = 0.1) UIS is unstable but there are neutrally 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 of the UIS with R0 = 4.3 × {10 −7 , 10 −9 , 10 −11 , 10 −13 }.The origin 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 R0.Note the logarithmic divergence of transient time T ∼ ln R0 (consistent with heteroclinicity).

FIG. 6 .
FIG.6.Bifurcation diagram for Eq.(15) with N = 4 oscillators and c2 = 3, ǫ = 0.1.Solid (dashed) lines represent stable (unstable) solutions.In the case of UIS and NUIS the solution depicted must be understood as the one observed under arbitrarily weak noise (there is continuum of neutral solutions with Q larger than the solution depicted).The saddle QPS orbit was continued by means of a Newton-Raphson algorithm, and the values of R and Q assigned in the diagram correspond to their time averages.

FIG. 7 .
FIG. 7. Stability boundaries of (a) UIS and (b) FS obtained exactly and from phase approximations, for c2 = 3.The solid line correspond to the exact boundary of the MF-CGLE (1), while dotted and dashed lines correspond to second-and third-order phase approximations, respectively.Blue lines are obtained from (15) and(29).Orange lines are the results if prior to phase reduction the MF-CGLE is transformed into(33), performing an isochron-base phase reduction in powers of κ.