Quasi phase reduction of all-to-all strongly coupled $\lambda-\omega$ oscillators near incoherent states

The dynamics of an ensemble of $N$ weakly coupled limit-cycle oscillators can be captured by their $N$ phases using standard phase reduction techniques. However, it is a phenomenological fact that all-to-all strongly coupled limit-cycle oscillators may behave as"quasiphase oscillators", evidencing the need of novel reduction strategies. In response to this, we introduce here quasi phase reduction (QPR), a scheme suited for identical oscillators with polar symmetry ($\lambda-\omega$ systems). Applying QPR we achieve a reduction to $N+2$ degrees of freedom: $N$ phase oscillators interacting through one independent complex variable. This"quasi phase model"is asymptotically valid in the neighborhood of incoherent states, irrespective of the coupling strength. The effectiveness of QPR is illustrated in a particular case, an ensemble of Stuart-Landau oscillators, obtaining exact stability boundaries of uniform and nonuniform incoherent states for a variety of couplings. An extension of QPR beyond the neighborhood of incoherence, valid for moderate coupling, is also explored. Finally, a general QPR model with $N+2M$ degrees of freedom is obtained for coupling through the first $M$ harmonics.


I. INTRODUCTION
Dynamical reduction is a concept of paramount importance in nonlinear dynamics [1]. It may operate reducing the number of degrees of freedom or transforming the evolution equations into a canonical form. Classical reduction techniques are adiabatic elimination [2], center-manifold reduction [3] and phase reduction [4][5][6][7]. The latter has been crucial configuring our comprehension of oscillatory media and coupled self-sustained oscillators.
Large ensembles of coupled self-sustained oscillators are found in a variety of domains ranging from biology and technology to the social sciences, see e.g. [8][9][10][11] and references therein. It is well established that if the coupling among N limit-cycle oscillators is weak then phase reduction can be applied [4], and the dynamics becomes reliably described by N phase oscillators. This approach yields a minimal description of emergent phenomena in all-to-all coupled oscillators as, for instance, collective synchronization [12][13][14], quasiperiodic partial synchronization (QPS) [15,16] or nonuniform incoherent states (NUISs) [17]. If the coupling is strong, however, phase reduction is not applicable as evidenced by several forms of collective chaos, which manifestly elude phase reduction [18,19]. However, there are situations in which the oscillators are strongly coupled, but still, resemble phase oscillators as their ordering on top of a closed curve is preserved as time evolves. Apart from the uniform incoherent state (UIS) and NUISs, there are more complex phenomena such as QPS, modulated QPS or pure collective chaos in which identical oscillators behave as "quasiphase oscillators" on top of an unsteady closed curve [19][20][21]. Recent advances extending the standard phase reduction beyond the first order do not appear to be practical enough even to cover the moderate coupling regime [17,22]. Alternative methods based on phaseamplitude reduction or isostables fall short in the dimensionality reduction [6,[22][23][24].
In this paper we present quasi phase reduction (QPR), a method to capture-with the minimum number of degrees of freedom-the dynamics of all-to-all strongly coupled limitcycle oscillators near incoherent states. The method only covers λ − ω oscillators (two-dimensional systems with polar symmetry), but still, it is conceptually appealing since it yields a significant dimensionality reduction from 2N to N + 2 degrees of freedom. The reduced system consists of N phase oscillators (as in standard phase reduction) and one complexvalued variable mediating the interactions. Thereupon we can calculate analytically the stability boundary of incoherent states. Moreover, we explore an extension of QPR, keeping the N + 2 degrees of freedom, which correctly pinpoints a saddle QPS at moderate coupling in a specific model. Finally, a general QPR with N + 2M degrees of freedom is derived for coupling through the M th harmonic. Throughout the paper, the correctness of our approach is confirmed by numerical simulations with a popular λ − ω system called Stuart-Landau oscillator.
The paper is organized as follows. In Sec. II we introduce the λ − ω oscillator and the isochrons. Incoherent states in a particular system of globally coupled λ − ω oscillators are reviewed in Sec. III for illustrative purposes. Section IV presents QPR for a family of coupling functions. The results in Sec. IV are applied to Stuart-Landau oscillators in Sec. V. Sections VI and VII extend the results in Sec. IV beyond the lowest order, and to other coupling functions, respectively. The conclusions are summarized in Sec. VIII.
In this work we restrict to oscillators of the λ − ω type [9,25]. These are 2-dimensional systems with rotational symmetry, which admit the following representation of the evolution equations in polar coordinates: The overdot denotes time derivative as usual. Without lack of generality we assume the existence of a stable limit cycle at r = 1, i.e. λ(1) = 0. Moreover, the natural frequency of the oscillator is Ω = ω(1). The attraction rate to the limit cycle is given by the second Floquet exponent Λ = dλ dr r=1 < 0. Alternatively to Eq. (1), we may work with the complex variable A = re iφ , such that the dynamics of the λ − ω oscillator where function f satisfies f (Ae iα ) = e iα f (A). For simplicity, it is convenient to assume that f can be expressed as a series of the form where f n are complex coefficients. The existence of an attractive limit cycle of frequency Ω implies n n Re(f n ) = Λ and n f n = iΩ.

Isochrons
To account for the effect of perturbations, phase reduction approaches require extending the definition of the phase away from the limit cycle [4,5,7,12]. To do so, we seek a phase variable θ, such thatθ = Ω holds in the whole basin of attraction, not only on the limit cycle. The collection of points with the same phase or 'isochron' is defined based on the convergence to the same 'asymptotic phase' on the limit cycle. For λ − ω systems, polar symmetry yields a relation between the phase θ and the polar coordinates of the form [12] θ(r, φ) = φ − χ(r), with χ(1) = 0. The phase dynamics satisfiesθ =φ − dχ drṙ , and imposingθ = Ω, we solve the equation for χ(r): n Im(f n )r n − Ω n Re(f n )r n+1 dr (5) Depending on the specific oscillator type a closed analytical solution of χ(r) may exist or not. However, if deviations from the limit cycle are small it is enough to know the first coefficient of the Taylor expansion of χ around r = 1 . Differentiating Eq. (5) and evaluating the limit r → 1 by L'Hôpital's rule, we get: For later use, we condense this expression and a previous one into A. Stuart-Landau oscillators In this paper, we asses the validity of our theoretical findings with the Stuart-Landau oscillator: The limit cycle at |A| = 1 has a second Floquet exponent Λ = −2. Moreover, the isochrons are logarithmic spirals θ = φ − c 2 ln r, where c 2 is the so-called nonisochronicity parameter in Eq. (7). Therefore, χ 0 = c 2 for this system.

III. AN EXAMPLE: THE MEAN-FIELD COMPLEX GINZBURG-LANDAU EQUATION
Before presenting the QPR method, it is instructive to recall a well-studied model of globally coupled λ − ω oscillators in which incoherent states are observed. The mean-field complex Ginzburg-Landau equation (MF-CGLE) consist of N diffusively coupled Stuart-Landau oscillators [18,19]: here constants ǫ and c 1 determine the strength and the reactivity of the coupling, respectively; and A = 1 N N k=1 A k . As shown below, it is useful to absorb the local term ǫ(1+ic 1 )A j . Specifically, setting rescaling time (t → t 1−ǫ ) and going to a rotating frame with rescaled amplitude (A j → Aj √ 1−ǫ e −i(ǫc1+c2)t ), we get: At variance with the Stuart-Landau Eq. (7), the linear coefficient has nonzero imaginary part, as we have adopted a rotating frame such that Ω = 0. This is irrelevant but simplifies later calculations. In a broad region of parameter space the system (9) settles into an incoherent state (i.e. with zero mean field A = 0). This equality does not specifies the state of the system at it holds for a continuum of oscillator arrangements. Among them, the most prominent one is UIS corresponding to oscillators located over a circle with uniformly distributed phases (for a finite N , oscillators are evenly spaced deserving the name of splay state). For the remaining incoherent states there is a lack of uniformity in the oscillator distribution, and we use the acronym NUIS for them. (In a finite population this state is also called phase balanced configuration [27].) Figure 1(a) shows a partial phase diagram of Eq. (9) for a specific value of c 2 = 3. UIS is observed in the light (yellow) shaded region at the left of the black solid line [18,19]. In the other shaded region a NUIS settles spontaneously (UIS is unstable). The asterisk in the phase diagram indicates the parameter values for the snapshot of NUIS in Fig. 1(b). It is apparent in this figure that the oscillators are not evenly distributed, while A, represented by a red cross, settles at the origin.

IV. (N + 2)-QUASI PHASE REDUCTION
In this section we present our QPR method from 2N to N + 2 degrees of freedom for N strongly coupled λ − ω oscillators. The reduced system consists of N phases plus 2 global degrees of freedom, hence the name of QPR. The validity of the method requires a weak perturbation in the oscillators' motion, what holds in the neighborhood of the incoherent states (see details below), irrespective of the coupling strength.

A. Coupling
We will consider globally coupled oscillators: where κ is a positive coupling constant, and A denotes one or more mean fields of the set C A , A ⊆ C A . The set C A of mean fields is: Here |A| n A = 1 N N j=1 |A j | n A j and * stands for complex conjugation. Note that only the first harmonic in φ enters in the interaction. (The case with higher harmonics in the coupling is discussed in Sec. VII). Moreover, we demand the interaction function g in Eq. (10) to vanish when the mean fields in the argument vanish, i.e. g(A = 0) = 0. Among the possible couplings, the most preeminent one is diffusion, g(A) ∝ A, as in the MF-CGLE (9) introduced above [18,19]. Other examples with nonlinear coupling as g(A) ∝Ā + b|Ā| 2Ā and g(A) ∝Ā + b|A| 2 A have been considered in [28] and [29], respectively. It is also important to notice that we do not exclude symmetry breaking terms in the coupling such as Re(Ā), similarly to [30].

B. Preliminaries
The first step of the analysis is obtaining an equation for the dynamics of the phases. Equation (10) in polar coordinates writes,ṙ The phase dynamics is obtained through the change of variables in (4): here χ ′ (r j ) denotes the derivative of χ with respect r evaluated at r j . In order to reduce the dimensionality of the system we seek to remove the dependence on the radii r j . Using Eq. (12a) we write the evolution of a perturbation δr j off the limit cycle (r j = 1 + δr j ): It is obvious that the oscillators will be in the proximity of the limit cycle whenever κ|g(A)| ≪ −Λ. If this condition holds, the lowest order approximation is to set r j = 1 in (13) obtaining:θ This equation is not completely closed, as there are still dependences on the mean field(s) through g(A).

Small κ: Standard first-order phase reduction
To put our work in context, and for later comparison, we note that traditional first order phase reduction assumes κ ≪ −Λ, what automatically implies δr j ≃ 0 as noted above. Therefore the mean fields in g(A) can be approximated as where Z ≡ e iθ is the Kuramoto order parameter. Thus, at the lowest order, the coupling will only depend on Z. In this case we can make the replacement g(A) ≃Γ(Z) in Eq. (15) obtainingθ This system of N phase oscillators is the first-order phase reduction of (10). Unfortunately, this reduction works poorly if the coupling is not small; and even for asymptotically small coupling there are states of (10) not reproducible by Eq. (17).
Higher order terms proportional to κ 2 , κ 3 , etc. can be incorporated to (17) removing degeneracies and extending the validity of the phase model with N degrees of freedom [17]. However, if the coupling is strong, this procedure is either impractical (as the convergence rate of the series in powers of κ is not fast enough [17]) or plain wrong (if the expansion in κ is divergent). An alternative approach, suited for the strong coupling case, is presented next.
C. Small |g(A)|: Quasi phase reduction of Eq. (10) Incoherent states, the starting point of our analysis, are configurations of the oscillators compatible with A = 0. Accordingly, in an incoherent state, each oscillator evolves as if it experienced no input from the rest of the population. For a finite ensemble of oscillators, a typical incoherent state is the splay state, in which oscillators' phases are evenly spaced along the limit cycle [18,31]. The corresponding incoherent state in the thermodynamic limit (N → ∞) is the UIS, corresponding to a uniform phase density. In addition, if the oscillators' phases are non-uniformly distributed but the mean field still vanishes we designate the state as NUIS.
Strongly coupled oscillators may spontaneously settle into UIS or NUISs in wide regions of parameter space, see e.g. Fig. 1. Moreover, it is phenomenologically observed that there are also non-incoherent states in which strongly coupled oscillators behave as "quasiphase oscillators" [21], preserving their ordering on top of a closed curve that evolves in time. This occurs, in particular, in globally coupled Stuart-Landau oscillators when UIS loses its stability giving rise to a state called quasiperiodic partial synchrony, which after secondary instabilities yields pure collective chaos [19,21].
With the aim of describing the previous phenomena in a minimal way, we recall Eq. (15), since it suggests that some kind of perturbative approach in small g(A) is feasible, in analogy to small κ in standard phase reduction. As Eq. (15) is not closed due to g(A), we are tempted to consider g(A) as new variable. This is not the best choice as the evolution equation cannot be generally closed in terms of g(A). Instead, the complex variable B = A is the right choice, since, as shown below, any mean field |A| n A can be approximately expressed in terms of B and Z, and any g(A) in turn. Assuming the proximity of the oscillators to their fiducial limit cycles, r = 1 + δr, we expand φ = θ + χ 0 δr + O(δr 2 ). In this way the mean field B is Therefore, we can express the average δre iθ in terms of B and Z: and apply this equality to all the other mean fields, obtaining a linear dependence of |A| n A on B and Z: With the previous equation any g(A) can be approximated by a function of Z and B: Now the evolution of B is obtained averaging (10) over the whole population. Namely, The term f (A) is calculated using Eqs. (3), (6), and (18) : Finally, replacing (19) and (21) into (15) and (20) we obtain the (N + 2)-QPR of the globally coupled oscillator system in Eq. (10):θ This equation is the main result of this paper. Some important remarks follow.

Remarks on the
The QPR of (10) into (22) implies a drastic decrease in the number of degrees of freedom from 2N to N + 2: N phases plus a complex collective variable B. In contrast to standard phase reduction, there is an extra complex variable B. This is the key ingredient to make the strong coupling analyzable while preserving the population of phase oscillators. The theory is consistent since QPR (22) boils down to the standard phase reduction (17) in the κ → 0 limit. To see this, set Ω = 0 in (22) by going to a rotating frame (θ ′ j , B ′ ) = (θ j − Ωt, Be −iΩt ) if necessary, and note that B(t) → Z(t) as κ → 0 in Eq. (22b). In this way, Eq. (22a) reduces to (17) since Γ(Z, Z) =Γ(Z), cf. Eqs. (16) and (18).
Equation (22) can be regarded as a population of phase oscillators coupled through a sort of external medium B. Indeed, a similar model is obtained applying ordinary phase reduction (assuming weak coupling) to a model of 'dynamical quorum sensing' in which oscillators are coupled through a medium with intrinsic dynamics [32]. Here, in sharp contrast, there is no 'medium' in the original system (10), but our QPR endows the mean field of a virtual dynamical equation.
An important feature of Eq. (22) is that it is a quasiintegrable model that can be analyzed within the framework of the Watanabe-Strogatz theory [33,34]. Given a particular initial condition there are N − 3 constants of motion determining the fate of the system. This degeneracy of the model is not present in (10). Still, the system in Eq. (22) is useful at least because of two reasons: (i) we can use it to determine the stability (boundary) of incoherent states analytically, see next sections; and (ii) it is the starting point for higher-order QPR, see Sec. VI.

Stability of incoherent states
Equation (22) is the QPR of (10), irrespective of the number N of oscillators. In this section, we take the thermodynamic limit (N → ∞), and analyze the stability boundary of the incoherent states. The analysis requires defining a density ρ such that ρ(θ, t)dθ is the fraction of oscillators with phases between θ and θ + dθ at time t. Additionally, we impose the normalization condition 2π 0 ρ(θ, t)dθ = 1. The Kuramoto order parameter is now Z = 2π 0 ρ(θ, t)e iθ dθ. The oscillator density ρ obeys the continuity equation because of the conservation of the number of oscillators: This is a nonlinear equation since v =θ depends on ρ.
According to Eq. (18), all |A| n A are linear combinations of B and Z. Therefore all states with B = Z = 0 are incoherent states since Γ(0, 0) = 0. Obviously, there is an infinity of phase densities compatible with Z = 0, which rotate uniformly: ρ incoh (θ, t) = ρ incoh (θ − Ωt). Notably, it will be shown below that not all incoherent states become unstable simultaneously.
The analysis proceeds introducing the Fourier expansion of ρ: with coefficients ρ 0 = 1 and ρ −m = ρ * m . Inserting (24) into (23), and noting that Z = ρ 1 , we may rewrite our model (22) in Fourier space: In light of these equations it becomes apparent the existence of an infinite set of incoherent solutions characterized by ρ 1 = B = 0, and ρ m≥2 =ρ m e imΩt with arbitraryρ m≥2 . We distinguish between UIS, corresponding toρ m =0 = 0, and the remaining set of NUISs. The linear stability of (N)UIS is determined considering the evolution of infinitesimal perturbations of the form ρ m = (ρ m + δρ m )e imΩt and B = δBe iΩt . The linearization of Eq. (25) turns out to be: The right-hand sides of these equations only include perturbations in the subspace spanned by ρ 1 and B; note the shorthand notation δ = (δρ 1 , δρ * 1 , δB, δB * ) T , and the gradients ∇Γ defined in this subspace and evaluated at ρ 1 = B = 0. We have then an infinite set of vanishing eigenvalues corresponding to eigenvectors with δB = δρ 1 = 0. The naive expectation is that these neutral modes should decay to zero under arbitrarily weak noise, as observed in the UIS of the MF-CGLE [19]. Actually, this is not necessarily the case, as a specific example in the next section shows.
Hence, according to Eq. (26), the relevant infinitesimal instabilities develop in the subspace spanned by ρ 1 and B. We are led to analyze the 4 × 4 Jacobian matrix ruling the dynamics of δB and δρ 1 . In this Jacobian only the modeρ 2 andρ * 2 are present and it can be shown that the stability of all incoherent states can be classified by the value of |ρ 2 | = Q. That is, the stability boundary of an incoherent state depends exclusively on Q (other higher-order modes are irrelevant). This result was already proved in a particular case [18,35], but QPR shows that it is a general property of the coupling via the mean fields in C A .
Finally, we want to stress that the stability boundaries of (N)UIS obtained from (26) exactly match those of the original system (10). The QPR is exact in the limit g(A) → 0 that is precisely where the instability takes place.

V. QUASI PHASE REDUCTION FOR STUART-LANDAU OSCILLATORS
A. Linear coupling: mean-field complex Ginzburg-Landau equation A simple system to illustrate and test our previous findings is MF-CGLE presented in Sec. III. Written as in (9) the values of Λ = −2 and χ 0 = c 2 remain those indicated in Sec. II A and given that g(A) = (1 + ic 1 )A, it is straightforward to obtain Γ(Z, B) = (1+ic 1 )B. Hence, the quasi phase reduced model (22) becomes: (27) is similar to the Kuramoto-Sakaguchi model, but with the phase oscillators coupled through B instead of Z. Only in the limit κ → 0, B approaches Z and the standard first-order phase reduction is recovered [19].

Numerical results: Transient dynamics
To confirm the correctness of our approach we compare the transient behavior of the MF-CGLE, as written in Eq. (9), with its QPR (27). We track the evolution of the meanfield Z = e iθ for both systems near incoherent states, noting that for the MF-CGLE this average is e i(φ−c2 ln r) . In Fig. 2(ad) we initialized N = 50 oscillators randomly on the unit circle, i.e. near the UIS. The parameters used in Figs. 2(a,b) and 2(c,d) correspond to stable and unstable UIS, respectively.
The stability properties of UIS, decay/growth rate and oscillation frequency, are perfectly captured by the QPR equations (27). As expected, in Fig. 2(d) after a certain time interval, the mean field |Z| grows too much and the QPR equations become inaccurate (the MF-CGLE approaches a saddle quasiperiodic partial synchrony and eventually decays to a NUIS). In Figs. 2(e) and 2(f) we show that QPR also gives a good description of NUISs. With the same parameters that in Fig. 2(d), the oscillators were randomly set in the phase interval [0, π 2 ] ∪ [π, 3 π 2 ] of the unit circle. In this way B = Z ≃ 0 but Q = e i2θ ≃ 2/π. We can see in Fig. 2(f) that, as time evolves, Z decays to zero but Q converges to a nonzero constant value, because UIS is unstable but NUISs with large enough Q values are not.

Uniform incoherent state
A closed formula for the stability boundary of the UIS was already found in [18,19], so here we just wish to evidence how QPR permits to obtain it in a simple way. As mentioned above only the evolution of δρ 1 and δB must taken into account in Eq. (26). As ρ 2 = 0 in the UIS, we get: The characteristic equation is: The locus of the (oscillatory) instability is determined imposing λ = iΩ c . The critical coupling κ 0 satisfies: in agreement with [18,19] (be aware of the different parametrizations in each work).

Non-uniform incoherent state
As with the UIS, the stability boundary of each NUIS is determined through the evolution of δρ 1 and δB. Using Eq. (26) and inserting the specific value ofρ 2 = Q ≤ 1: The characteristic polynomial is a four degree polynomial with real coefficients: Although the zeros cannot be computed, the Routh-Hurwitz criterion [36] permits knowing if there is at least one root with nonnegative real part. For the fourth order polynomial P 4 (λ) plus two other inequalities that are always fulfilled. Equations (30) are precisely the exact Q-dependent NUIS stability boundaries of (9) [35].

The effect of arbitrarily weak noise
As already noted in Sec. IV C 2 a particular (N)UIS may be either unstable or neutrally stable, but not asymptotically stable. Thus, in the MF-CGLE, a continuum of neutrally stable incoherent states coexist in regions of parameter space. Hence, the question is the selective effect of arbitrarily weak noise. The color shading in the phase diagram of Fig. 1(a) has been made from Eqs. (28) and (30) on the assumption that the system adopts a phase density with Q = Q * , where Q * is the smallest Q value among all neutrally stable NUISs. In consistence a neutral UIS is attracting in the presence of weak noise, see [19]. As may be seen in Fig. 3, in the region where UIS is unstable the values of Q observed match almost perfectly with Q * , depicted by a black solid line.
Less intuitive is the behavior of the remaining modes, Z m (m > 2), that are irrelevant in the stability analysis. According to our Fig. 3 in the UIS region, all Z m go to zero, while in the NUIS region this is only the case for odd m. The even modes grow as Q increases. The last NUIS to destabilize is Q = 1, and corresponds to two equally populated point clusters in antiphase, i.e. a bi-delta phase density (Z 2 = Z 4 = Z 6 = · · · = e iξ ).

B. Nonlinear coupling, g(A) ∝ |A| n A
Recent papers by Krischer and coworkers [29,37] have studied chimera states in the MF-CGLE with an extra coupling term proportional to |A| 2 A. Here, instead of embarking on the exploration of the high-dimensional parameter space of that system, each coupling of the form |A| n A (n ∈ Z) is analyzed separately. Thus, the systems under consideration are:Ȧ (31) where parameter γ n is a real constant. In the particular case n = 0, Eq. (31) becomes the MF-CGLE (9), and γ 0 = c 1 accordingly.
Deriving the QPR of (31) requires calculating the function Γ(Z, B). Using (18) the result is straightforward: where the subscript n is used to indicate the dependence on the specific coupling considered. Finally, inserting Γ n into (22) we obtain the QPR of (31): and Z = Re iΨ . Prior to determining the exact stability boundaries of (N)UIS from Eq. (33), let us see what the standard first-order phase reduction predicts. For this, we take the limit κ → 0 + , observing that B collapses into Z, and Eq. (33a) becomes the Kuramoto-Sakaguchi model in Eq. (17). The crossover from perfect synchrony to incoherence, given by the Benjamin-Feir-Newell criterion (1 + c 1 c 2 = 0) in the diffusive case, generalizes to: by virtue of Eq. (16). Therefore, in a phase diagram of system (31) including axis κ, the stability boundaries of (N)UIS are expected to emanate from γ n = −c −1 2 at κ = 0. The exact stability boundaries of UIS and NUISs in the thermodynamic limit are obtained from (33) as explained in previous sections. The stability boundary of UIS is: We can see in the limit κ → 0 we recover (34). The stability boundary of a Q-dependent NUIS can be computed as was done in the linear coupling case; the interested reader can find its expression in the Appendix.
In Fig. 4 the stability boundaries lines of UIS and NUIS are depicted for five different values of n = −2, −1, . . . , 2. Taking n = 0 in panel (c) as the reference case, we see that augmenting n shrinks the region of incoherence. On the contrary, for n = −1 stable NUISs reach larger κ values, while the UIS region remains mostly unchanged. The boundaries for other negative n values are similar to those for n = −2 in Fig. 4(a). It is remarkable that for all n values there are regions in the parameter space where UIS is unstable, but certain NUISs are not. This means that, at least for certain initial conditions, the system spontaneously converges to a NUIS. According to FIG. 4. Partial phase diagrams for populations of Stuart-Landau oscillators with nonisochronicity parameter c2 = 3, coupled through g(A) = (1 + iγn)|A| n A with n = −2, 1, 0, 1, 2, in panels (a)-(e). The stability boundaries of UIS and the NUIS with Q = 1 are depicted by black and red lines, respectively. UIS is stable in the yellow region. In the shaded region the intensity of the red color indicates Q * (the smallest Q value among the non-unstable NUISs).
our numerical simulations, and as reasoned above, under weak noise the non-unstable NUIS with the smallest Q value is observed. In addition, save for n = −1, 0, there are also regions for small enough γ n where UIS is the last incoherent state to become unstable.
We want to remark that all the stability boundaries calculated are the exact results for (31). To our knowledge, only the case of the MF-CGLE had been solved so far. We believe using QPR (33) is the most effective method for computing these boundaries.

C. Other couplings
In this subsection we want to make some comments on the couplings where the (N + 2)-QPR scheme so far presented can be applied.
It is straightforward to consider a combination of nonlinear couplings such as: with complex σ n and µ n . In this case, Γ(Z, B) is simply a sum over terms like the bracketed part in the right-hand side of Eq. (32) and their complex conjugates.
Other nonlinear coupling considered in [28,38]: can be treated analogously to other couplings. Nonetheless, the stability boundaries of UIS and NUISs of this system are the same that for Eq. (9) because the nonlinear term does not contribute if |B| = |A| ≪ 1.
Our QPR approach does not exclude systems with nonlinear delayed feedback and/or couplings as h(A j )g(A) (provided h has polar symmetry h(A j e iφ ) = e iφ h(A j ) [39]), similar to those studied in [40].
Finally, the case of a scalar coupling like g(A) ∝ Re(A) ∝ A+Ā * , is particularly simple, as QPR may be further reduced to only one real-valued global variable, i.e. N + 1 degrees of freedom in total.

VI. EXPLORING THE NEXT ORDER OF QPR
As occurs with standard first-order phase reduction, extending the theory to the next order in the QPR scheme is not a trivial task. For QPR, a systematic expansion is even more troublesome as there is not a small coupling parameter (say, κ), but a small field g(A). This should be the goal of future works, but we think it may be instructive to pinpoint the difficult points, and examine the workable limit of small coupling.
The first step is to expand Eq. (13) to order δr j , what yields an augmented version of Eq. (15): where χ 1 = d 2 χ(r) dr 2 | r=1 . The value of δr j evolves in time as dictated by Eq. (14), which is coupled to θ j and to the mean field A. It is not obvious how to proceed next since A is not static.
A. Small κ Inspecting Eq. (14), we realize that if κ is small, then δr j (t) adjusts quickly to the current mean field: We can insert (37) With the new term, proportional to κ 2 , the Watanabe-Strogatz theory [33] cannot be applied. This is not a surprise, since the original model is not quasi-integrable. To proceed with the analysis, function g has to be written in terms of the the mean fields Z, B, and maybe others. To the lowest order we simply adopt the function Γ(Z, B) obtained above.

Mean-field complex Ginzburg-Landau equation
Let us see, how Eq. (38) applies to the particular case of the MF-CGLE, Eq. (9). As the unit oscillator is the Stuart-Landau oscillator, we insert χ 0 = −χ 1 = c 2 into Eq. (38). Moreover, we keep the evolution for B as before. This results in an extended QPR model: here β = arg(1 − c 2 1 + 2ic 1 ). Next, we test (39) against numerical simulations. We select constants c 1 , c 2 , and κ such that UIS and full synchrony are both unstable, but NUISs with Q above a certain value have not destabilized. As observed in [17], for small and moderate κ values there is a heteroclinic connection between UIS and a saddle QPS. Recall, that, in a QPS state the oscillator density rotates uniformly (as Z, Z 2 , etc, accordingly), but each individual oscillator exhibits quasiperiodic motion. For the numerical test in Fig. 5 we initialize N = 100 Stuart-Landau oscillators randomly in the unit circle for the full model (9), as well as the (N + 2)-QPR (27), and the extended (N + 2)-QPR (39) with identical initial phases and B = Z value. Two values of the coupling are selected κ = 0.2 and 0.5, in Figs. 5(a,b) and 5(c,d), respectively. The heteroclinic connection with the saddle QPS is captured by the extended model (39), in contrast to (27), which only reproduces the exponential instability of UIS. For both, the MF-CGLE (9) and Eq. (39), the final state is a NUIS. Unsurprisingly, the extended QPR (39) is more accurate for κ = 0.2 than for κ = 0.5, since we assumed a small κ in its derivation.
For κ values larger than those in Fig. 5 there is not a saddle QPS, but instead stable QPS branches off from UIS as it becomes unstable [20,21] . Remarkably, this also occurs for the extended model (39) at large enough κ (not shown).

VII. (N + 2M )-QUASI PHASE REDUCTION: HIGHER-ORDER HARMONICS
In this section we analyze QPR when we let the oscillators interact through higher-order harmonics: where A M are mean fields belonging to the set: These mean fields are the first M harmonics in φ. We show next that, provided the subset of (41) with m = M is not empty, the QPR of Eq. (40) possesses N + 2M degrees of freedom. Otherwise said, the largest harmonic of φ in the coupling determines the number of degrees of freedom of QPR.
where Z m = e imθ is the m-th Kuramoto-Daido order parameter. We can express δre imθ in terms of B m and Z m : Applying this equality to the averages |A| n A m with arbitrary n value yields: |A| n A m = r n+m e imφ ≃ Z m + (n + m + imχ 0 )δre imθ The evolution of the phases is, therefore, linked to the set of complex mean fields {B m } m=1,...,M , whose evolution equations remain to be determined. Recalling (40), we get: As a final note, we mention that is also possible to deal with a coupling function g |A| , where g is any function. In this case QPR can be accomplished using |A| and A as the collective variables (N + 3 degrees of freedom in total).

VIII. CONCLUSIONS
Phase reduction is a powerful technique that has deeply shaped our knowledge on the dynamics of oscillator ensembles. In spite of its enormous success, the description enabled by reduced phase models breaks down if the coupling between the oscillators is not weak. Recently, some works have extended standard phase reduction [17,22,41] in a perturbative fashion in the coupling constant. These approaches are however condemned to fail for strong coupling.
In this work we have given a new twist to the concept of phase reduction, introducing QPR for all-to-all strongly coupled λ − ω oscillators. This reduction procedure exploits the smallness of the collective oscillations near incoherent states, independently of the coupling strength. The reduced model has N + 2M degrees of freedom model corresponding to M dynamical complex variables mediating the interactions of N phase oscillators, akin to dynamical quorum-sensing models. We have studied in detail the case M = 1, corresponding to interactions via the first harmonic of the angle. Explicit stability boundaries for uniform and nonuniform incoherent states have been obtained for ensembles of Stuart-Landau oscillators.
Finally, an extension of QPR beyond the lowest order has been obtained for weak coupling. Nonetheless, a genuine well-controlled expansion to the next order remains to be developed. In parallel with this, some sort of generalization from global to more complex coupling topologies appears to be possible as well. All in all, we deem QPR as a promising path towards a comprehensive theory of collective phenomena in oscillator ensembles.