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

TL;DR

This paper develops a second-order phase reduction method for the mean-field complex Ginzburg-Landau equation, capturing complex collective dynamics like chaos that first-order models miss.

Contribution

It introduces an isochron-based scheme for second-order phase approximation, extending to third order, revealing multi-body interactions in oscillator ensembles.

Findings

Second-order approximation reproduces weak coupling dynamics of MF-CGLE.

Higher-order terms introduce multi-body interactions.

Multi-body interactions are linked to collective chaos.

Abstract

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.

Equations

$${\dot{A}}_j=A_j-(1+i{c}_2)|A_j{|}^2{A}_j+ϵ(1+i{c}_1)(\bar{A}-A_j).$$

$${ϵ}_s=-\frac{2(1+c_1{c}_2)}{1+c_1^2},$$

$${ϵ}_0(2{ϵ}_0-1)c_1^2+4({ϵ}_0-1)(2{ϵ}_0-1)c_1{c}_2$$

$$-{ϵ}_0({ϵ}_0-1)c_2^2+(3{ϵ}_0-2{)}^2=0.$$

$$[{ϵ}_Q(2{ϵ}_Q-1)c_1^2+4({ϵ}_Q-1)(2{ϵ}_Q-1)c_1{c}_2$$

$$-{ϵ}_Q({ϵ}_Q-1)c_2^2+(3{ϵ}_Q-2{)}^2][(2-3{ϵ}_Q{)}^2+{ϵ}_Q^2{c}_1^2]$$

$$=Q^2{ϵ}_Q(1-{ϵ}_Q)(3{ϵ}_Q-2{)}^2(c_1^2+1)(c_2^2+1).$$

$${\dot{\theta }}_j=\Omega +ϵ\eta R\sin (\Psi -{\theta }_j+\alpha ),$$

$$\alpha =\arg [1+c_1{c}_2+(c_1-c_2)i].$$

$$\theta (r,\phi )=\phi -c_2\ln r.$$

$$\dot{r_j}=r_j(1-ϵ-r_j^2)$$

$$\displaystyle+\frac{\epsilon}{N}\sum_{k=1}^{N}r_{k}\bigg{[}\cos(\varphi_{k}-\varphi_{j})-c_{1}\sin(\varphi_{k}-\varphi_{j})\bigg{]},$$

$$\dot{{\phi }_j}=-c_2{r}_j^2-ϵc_1$$

$$\displaystyle+\frac{\epsilon}{Nr_{j}}\sum_{k=1}^{N}r_{k}\bigg{[}c_{1}\cos(\varphi_{k}-\varphi_{j})+\sin(\varphi_{k}-\varphi_{j})\bigg{]}.$$

$${\dot{r}}_j=f(r_j)+ϵg_j(\mathbf{r},\mathbf{\theta }),$$

$${\dot{\theta }}_j=ϵh_j(\mathbf{r},\mathbf{\theta }).$$

$$g_j(\mathbf{r},\mathbf{\theta })$$

$$h_j(\mathbf{r},\mathbf{\theta })$$

$${\dot{\theta }}_j=ϵh_j({\mathbf{r}}^{(0)},\mathbf{\theta })+{ϵ}^2\left({\nabla }_r{h}_j({\mathbf{r}}^{(0)},\mathbf{\theta })\right)\cdot {\mathbf{r}}^{(1)}$$

$$+{ϵ}^3\left[\left({\nabla }_r{h}_j({\mathbf{r}}^{(0)},\mathbf{\theta })\right)\cdot {\mathbf{r}}^{(2)}+(\mathbf{Mrr}{)}_j\right]+\cdots ,$$

$${\dot{r}}_j^{(1)}=f^{\prime }(r_j^{(0)})r_j^{(1)}+g_j({\mathbf{r}}^{(0)},\mathbf{\theta }).$$

$${\dot{r}}_j^{(1)}=({\nabla }_{\theta }{r}_j^{(0)})\cdot \mathbf{h}=0.$$

$$r_j^{(1)}=-\frac{g_j({\mathbf{r}}^{(0)},\mathbf{\theta })}{f^{\prime }(r_j^{(0)})}=\frac{g_j({\mathbf{r}}^{(0)},\mathbf{\theta })}{2},$$

$${\dot{\theta }}_j=\Omega +ϵ\eta R\sin (\Psi -{\theta }_j+\alpha )+\frac{{ϵ}^2{\eta }^2}{4}\left\{RQ\sin (\Phi -\Psi -{\theta }_j)-\sum _{m=1}^2(-R{)}^m\sin [m(\Psi -{\theta }_j)+\beta ]\right\}+O({ϵ}^3).$$

$$\beta =\arg (1-c_1^2+2c_{1}i),$$

$$R^2\sin [2(\Psi -{\theta }_j)+\beta ]=\frac{1}{N^2}\sum _{k,l}\sin ({\theta }_k+{\theta }_l-2{\theta }_j+\beta )$$

$$RQ\sin (\Phi -\Psi -{\theta }_j)=\frac{1}{N^2}\sum _{k,l}\sin (2{\theta }_k-{\theta }_l-{\theta }_j)$$

$${ϵ}_s=\frac{-2(1+c_1{c}_2)}{c_1^2(1+c_2^2)}$$

$${\partial }_t\rho (\theta ,t)+{\partial }_{\theta }[v(\theta )\rho (\theta ,t)]=0.$$

$$\rho (\theta ,t)=\frac{1}{2\pi }\sum _{n=-\infty }^{\infty }{\rho }_n{e}^{in\theta },$$

$$Z_n={\int }_0^{2\pi }\rho (\theta ,t)e^{in\theta }d\theta ={\rho }_{-n}.$$

$${\dot{\rho }}_n=\frac{n}{2}ϵ\eta \left\{e^{-i\alpha }{\rho }_1{\rho }_{n-1}-e^{i\alpha }{\rho }_1^{\ast }{\rho }_{n+1}+\frac{ϵ\eta }{4}\left[e^{-i\beta }{\rho }_1({\rho }_{n-1}-{\rho }_1{\rho }_{n-2})-e^{i\beta }{\rho }_1^{\ast }({\rho }_{n+1}-{\rho }_1^{\ast }{\rho }_{n+2})-{\rho }_2^{\ast }{\rho }_1{\rho }_{n+1}+{\rho }_2{\rho }_1^{\ast }{\rho }_{n-1}\right]\right\}$$

$$\frac{d}{dt}\delta {\rho }_1=\frac{ϵ\eta }{2}\left[e^{-i\alpha }+\frac{ϵ\eta e^{-i\beta }}{4}\right]\delta {\rho }_1.$$

$${ϵ}_0=\frac{4(1+c_1{c}_2)}{(c_1^2-1)(1+c_2^2)}.$$

$$\frac{d}{dt}\begin{pmatrix}\delta\rho_{1}\\ \delta\rho_{1}^{*}\\ \delta\rho_{2}\\ \delta\rho_{2}^{*}\\ \vdots\end{pmatrix}=\left(\begin{array}[]{cc|ccc}\bullet&\bullet&0&0&\cdots\\ \bullet&\bullet&0&0&\cdots\\ \hline\cr\bullet&\bullet&0&0&\cdots\\ \bullet&\bullet&0&0&\cdots\\ \vdots&\vdots&\vdots&\vdots&\ddots\\ \end{array}\right)\begin{pmatrix}\delta\rho_{1}\\ \delta\rho_{1}^{*}\\ \delta\rho_{2}\\ \delta\rho_{2}^{*}\\ \vdots\end{pmatrix},$$

Full text

I Introduction

Networks of nonlinear elements with oscillatory behavior (‘oscillators’) are found in a variety of disciplines, such as neuroscience or engineering Winfree (1980); Hoppensteadt and Izhikevich (1997); Strogatz (2003); Pikovsky et al. (2001). 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 Mertens and Weaver (2011a); *mertens11b, networks of (electro-)chemical oscillators Kiss et al. (2007); Totz et al. (2015), arrays of Josephson junctions Wiesenfeld and Swift (1995); Wiesenfeld et al. (1996) and globally coupled electrical self-oscillators Temirbayev et al. (2012); English et al. (2015), or nanoelectromechanical oscillators in a ring Matheny et al. (2019).

Applying a phase reduction method Winfree (1980); Kuramoto (1984); Nakao (2016); Monga et al. (2019) 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 Kuramoto (1984), 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  Tanaka and Aoyagi (2011), multi-body phase interactions naturally arise if the coupling is nonlinear Ashwin and Rodrigues (2016), see also Rosenblum and Pikovsky (2007). Instead, for linear pairwise coupling, three-body interactions are a distinctive element of second-order phase approximations, as recently highlighted in Matheny et al. (2019). Recognizing the ubiquity of multi-body interactions may also be important for reconstructing phase interactions from data Kralemann et al. (2011); *kralemann14.

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 assumption Winfree (1980); Kuramoto (1984); Pikovsky and Rosenblum (2015). 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) Hakim and Rappel (1992); Nakagawa and Kuramoto (1993, 1994, 1995); Chabanol et al. (1997); Banaji and Glendinning (1999); Banaji (2002); Daido and Nakanishi (2006, 2007); Takeuchi et al. (2009); Takeuchi and Chaté (2013); Sethia and Sen (2014); Ku et al. (2015); Rosenblum and Pikovsky (2015); Kemeth et al. (2019); Clusella and Politi (2019). 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 Hakim and Rappel (1992); Nakagawa and Kuramoto (1993, 1994, 1995); Banaji and Glendinning (1999); Takeuchi et al. (2009); Ku et al. (2015); Clusella and Politi (2019). Phase reduction of the MF-CGLE yields the Kuramoto-Sakaguchi model Kuramoto (1975, 1984); Nakagawa and Kuramoto (1993), 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 phase-reduction 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.

II mean-field complex Ginzburg-Landau equation

The MF-CGLE consists of $N$ diffusively coupled Stuart-Landau oscillators governed by $N$ coupled (complex-valued) ordinary differential equations:

[TABLE]

Here, $A_{j}=r_{j}e^{i\varphi_{j}}$ is a complex variable (index $j$ runs from 1 to $N$), and the mean field is $\bar{A}=N^{-1}\sum_{k=1}^{N}A_{k}$. Apart from the population size $N$, there are three free parameters in Eq. (1): $\epsilon$, $c_{1}$, and $c_{2}$. Parameter $\epsilon$, 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’ Pikovsky et al. (2001), generically appears from center manifold reduction Kuramoto (1984). 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}\to A_{j}e^{i\phi}$, and full permutation symmetry stemming from the mean-field coupling.

II.1 Phenomenology

For many parameter values, the global attractor of Eq. (1) is either full synchronization (FS) $A_{j}=\bar{A}=e^{-ic_{2}t}$ or one incoherent state with vanishing mean field $\bar{A}=0$. In the latter case the oscillators rotate freely, $A_{j}=\sqrt{1-\epsilon}\exp\{i[-c_{2}+\epsilon(c_{2}-c_{1})]t+\phi_{j}\}$. Among all the states compatible with $\bar{A}=0$, the most prominent one is the uniform incoherent state (UIS) in which the $\phi_{j}$ are uniformly distributed in the thermodynamic limit (for a finite ensemble, the $\phi_{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 Hakim and Rappel (1992); Chabanol et al. (1997).

In addition to FS, UIS, and NUIS, system (1) exhibits a rich repertoire of collective states including clustering Nakagawa and Kuramoto (1993); Banaji (2002); Daido and Nakanishi (2006, 2007); Ku et al. (2015); Kemeth et al. (2019), diffusion-induced inhomogeneity (or chimera) Daido and Nakanishi (2006, 2007); Sethia and Sen (2014), quasiperiodic partial synchronization (QPS) Nakagawa and Kuramoto (1993); Clusella and Politi (2019), as well as collective and microscopic chaos Hakim and Rappel (1992); Nakagawa and Kuramoto (1993, 1994, 1995); Banaji and Glendinning (1999); Takeuchi et al. (2009); Ku et al. (2015); Clusella and Politi (2019). In a QPS state, see e.g. Fig. 1(b), the mean field $\bar{A}$ 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 Nakagawa and Kuramoto (1995); Clusella and Politi (2019). 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 Nakagawa and Kuramoto (1993, 1995); Clusella and Politi (2019) 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 phase-describable states (NUIS, QPS, modulated QPS, pure collective chaos,etc.).

II.2 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 $\epsilon_{s}$ and $\epsilon_{0}$, the $\epsilon$ values of marginal linear stability for FS and UIS. Closed formulas for $\epsilon_{s}$ and $\epsilon_{0}$ are Hakim and Rappel (1992); Nakagawa and Kuramoto (1993):

[TABLE]

These formulas are also valid for finite ensembles, assuming $\epsilon_{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 Nakagawa and Kuramoto (1993) we choose to fix $c_{2}$, and display the loci of $\epsilon_{s}$ and $\epsilon_{0}$ in the parameter plane $(c_{1},\epsilon)$. 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 $\epsilon_{s}$ and $\epsilon_{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 Kuramoto (1984); Walgraef (1997); Pikovsky et al. (2001); Aranson and Kramer (2002). There is a critical value $c_{2}=\sqrt{3}=1.732\ldots$ at which the boundaries $\epsilon_{s}$ and $\epsilon_{0}$ become tangent at $\epsilon=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}\sum_{j}\exp(2i\varphi_{j})|$. The coupling constant $\epsilon_{Q}$ at which one particular NUIS becomes unstable was obtained in Ref. Chabanol et al. (1997):

[TABLE]

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\neq 1$, since the values of all ‘higher-order’ mean fields $f_{n}=|N^{-1}\sum_{j}\exp(ni\varphi_{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 $f_{n}=0$ ($n>2$).

Figures 2 (a) and (b) include a green-hatched region, adjacent to the UIS region at moderate $\epsilon$ values, where other phase-describable states are stable. These are QPS, modulated QPS and pure-collective chaos Nakagawa and Kuramoto (1995); Clusella and Politi (2019). We determined the boundary through simulations with $N=200$ oscillators, but the result is insensitive if a larger $N$ value is used.

II.3 First-order phase reduction: Kuramoto-Sakaguchi model

At the lowest order, applying the classical averaging technique Kuramoto (1975, 1984); Pikovsky et al. (2001) to Eq. (1) yields the Kuramoto-Sakaguchi model Sakaguchi and Kuramoto (1986). In this model, each oscillator is described by a phase $\theta_{j}$, and it is coupled to the other ones by pairwise interactions of the form $\sin(\theta_{i}-\theta_{j}+\alpha)$. In agreement with the mean-field character of the system, oscillators are coupled through the Kuramoto order parameter $Z_{1}\equiv R\,e^{i\Psi}=N^{-1}\sum_{k=1}^{N}e^{i\theta_{k}}$, such that the ordinary differential equations governing the dynamics are:

[TABLE]

with constants $\Omega\equiv-c_{2}+\epsilon(c_{2}-c_{1})$, $\eta\equiv\sqrt{(1+c_{2}^{2})(1+c_{1}^{2})}$, and phase lag

[TABLE]

Equation (5) is the disorder-free version of the paradigmatic Kuramoto-Sakaguchi model Kuramoto (1984); Sakaguchi and Kuramoto (1986) and related models Montbrió and Pazó (2011a); *MP11; *PM11. 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=1$— is 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 Nakagawa and Kuramoto (1993).

As discussed above, the MF-CGLE has much richer dynamics than its first-order phase reduction (5), even arbitrarily close to the $\epsilon=0$ limit. Therefore, it is mandatory to extend the phase reduction to order $O(\epsilon^{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},\epsilon)=(-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 $|\epsilon|\ll 1$. Moreover, it might serve to shed light on the mechanisms behind complex dynamics found (so far) for moderate $\epsilon$ 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 $\epsilon^{3}$. In this section we give the details of our phase reduction calculation. We anticipate that the results at second and third order in $\epsilon$ correspond to Eqs. (15) and (29) below.

III.1 Isochrons

The concept of isochron Winfree (1974); Guckenheimer (1975) is the cornerstone of phase reduction methods Winfree (1980); Kuramoto (1984). 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’ Coddington and Levinson (1955)). For the Stuart-Landau oscillator, polar coordinates $(r,\varphi)$ relate to the phase $\theta$ according to Kuramoto (1984); Pikovsky et al. (2001):

[TABLE]

As mentioned above, on the limit cycle ($r=1$), $\theta=\varphi$. 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.

III.2 Isochron-based phase reduction

We continue the analysis writing Eq. (1) in polar coordinates:

[TABLE]

After the change of variables $(r_{j},\varphi_{j})\to(r_{j},\theta_{j})$ through Eq. (7), we get:

[TABLE]

Here, we have also implemented the transformation $\theta_{j}\to\theta_{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 $\epsilon$, while the $r_{j}$ are fast variables that become enslaved to the dynamics of $\theta_{j}$.

In Eq. (10) $f(r)=r(1-r^{2})$, and functions $g_{j}$ and $h_{j}$ depend on the vectors $\mathbf{r}=(r_{1},r_{2},\ldots,r_{N})^{T}$ and $\boldsymbol{\theta}=(\theta_{1},\theta_{2},\ldots,\theta_{N})^{T}$ as follows,

[TABLE]

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 $\epsilon$. Based on the empirical observation that, at small $\epsilon$ 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}(\theta_{1},\theta_{2},\ldots,\theta_{N})$. We also postulate an expansion in powers of $\epsilon$ for the radii: $r=r_{j}^{(0)}+\epsilon r_{j}^{(1)}+\epsilon^{2}r_{j}^{(2)}+\cdots$; or in vector notation $\mathbf{r}=\mathbf{r}^{(0)}+\epsilon\mathbf{r}^{(1)}+\epsilon^{2}\mathbf{r}^{(2)}+\cdots$. Equation (10b) for $\theta_{j}$ becomes:

[TABLE]

where $\boldsymbol{\nabla}_{r}\equiv(\partial_{r_{1}},\partial_{r_{2}},\ldots,\partial_{r_{N}})$ and $(\mathbf{Mrr})_{j}\equiv\tfrac{1}{2!}\sum_{k,l}\partial_{r_{k}}\partial_{r_{l}}h_{j}(\mathbf{r}^{(0)},\boldsymbol{\theta})r_{k}^{(1)}r_{l}^{(1)}$. 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(\epsilon^{0})$ yields $r_{j}^{(0)}=1$, and (12) becomes (at the lowest order) the Kuramoto-Sakaguchi model (5). At order $\epsilon$:

[TABLE]

As $r_{j}$ depends exclusively on the phases, we can apply the chain rule: $\dot{r}_{j}=(\boldsymbol{\nabla}_{\theta}r_{j})\cdot\dot{\boldsymbol{\theta}}$. At order $\epsilon$, the time derivative vanishes:

[TABLE]

Hence Eq. (13) yields the result

[TABLE]

which can be inserted in (12) to obtain the $\epsilon^{2}$ contribution. Through elementary manipulations the second-order phase reduction of Eq. (1) can be condensed into this expression:

[TABLE]

The $O(\epsilon^{2})$ term depends on $Z_{1}$ as well as on the second Kuramoto-Daido order parameter Daido (1996) $Z_{2}\equiv Q\,e^{i\Phi}=N^{-1}\sum_{k=1}^{N}e^{2i\theta_{k}}$. To enhance the clarity of Eq. (15), we found it convenient to define a phase lag

[TABLE]

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(\epsilon^{0})$ term inside $\Omega$ 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(\epsilon^{2})$ contributions to Eq. (15), the first element of the sum ($m=1$) entails a parameter shift to the $O(\epsilon)$ interaction, and it is therefore irrelevant in qualitative terms. The other two terms in Eq. (15) correspond to three-body (i.e. non-pairwise) interactions:

[TABLE]

[TABLE]

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 $\epsilon$ 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 $\theta_{j}\to\theta_{j}+\phi$. For the sake of making the presentation simpler we assume constant $\Omega=0$, since this can always be achieved by going to a rotating frame $\theta_{j}\to\theta_{j}+\Omega t$.

IV.1 Full synchronization

The stability boundary of FS ($\theta_{j}=\Psi=\Phi/2$) is easily calculated. In particular, for infinite $N$ it is almost immediate: we simply assume one oscillator is infinitesimally perturbed, say the first one, $\theta_{1}=\Psi+\delta\theta_{1}$. The evolution of the perturbation obeys the linear equation $d\delta\theta_{1}/dt=\epsilon\eta\left[\cos\alpha+\frac{\epsilon\eta}{4}(1-\cos\beta)\right]\delta\theta_{1}$. At threshold ($d\delta\theta_{1}/dt=0$) the coupling satisfies:

[TABLE]

where we have written $\cos\alpha$ and $\cos\beta$ in terms of $c_{1}$ and $c_{2}$. 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 $\epsilon_{s}\to 0$, and deviates progressively from the FS boundary of the MF-CGLE (represented by a solid line) as $\epsilon_{s}$ increases.

IV.2 Incoherent states

We adopt the thermodynamic limit and define a density $\rho$ such that $\rho(\theta,t)d\theta$ is the fraction of oscillators with phases between $\theta$ and $\theta+d\theta$. Now $\theta\in[0,2\pi)$ is a cyclic variable, i.e. $\rho(\theta+2\pi,t)=\rho(\theta,t)$, and we impose the normalization condition $\int_{0}^{2\pi}\rho(\theta,t)d\theta=1$. The oscillator density $\rho$ obeys the continuity equation because of the conservation of the number of oscillators:

[TABLE]

Here $v=\dot{\theta}$ is the $\rho$-dependent velocity of an oscillator with phase $\theta$. We define the Fourier modes of $\rho$:

[TABLE]

with $\rho_{0}=1$ and $\rho_{n}=\rho_{-n}^{*}$. The mean fields $Z_{n}$ reduce to

[TABLE]

Inserting the Fourier expansion (21) into the continuity equation (20) allows us to rewrite our model in Fourier space:

[TABLE]

IV.2.1 Uniform incoherent state

The stability boundary of the UIS ($\rho(\theta)=(2\pi)^{-1}\Leftrightarrow\rho_{n\neq 0}=0$) is obtained linearizing the previous equation. It is easy to notice that only the first mode may destabilize. We have for $|\rho_{1}|\ll 1$:

[TABLE]

Neglecting the trivial marginal case $\epsilon=0$, the stability boundary satisfies $\cos\alpha+(1/4)\epsilon_{0}\eta\cos\beta=0$. Or in terms of $c_{1}$ and $c_{2}$:

[TABLE]

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.

IV.2.2 Nonuniform incoherent states

According to (22), in an incoherent state ($\rho_{1}=0$) higher-order modes are at rest: $\dot{\rho}_{n}=0$ ($n>2$). The linearization of (22) around $\rho_{1}=0$ and $\rho_{n}\neq 0$ ($|n|\geq 2$) is (schematically) as follows:

[TABLE]

where the $\bullet$ 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 $\delta\rho_{1}$ is hence, the only relevant one. The linear terms in $\delta\rho_{1}$ yield:

[TABLE]

All higher-order modes, save $\rho_{2}$, are absent in the equation. As $\dot{\rho}_{2}=0$, we can choose the coordinate axes such that $\rho_{2}=Q\in\mathbb{R}$. After some calculations we find that NUIS with a specific $Q$ value is marginally stable at:

[TABLE]

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 $\epsilon$ is set at a certain $\epsilon=\epsilon_{Q^{*}}$ the numerical integration of the system (either oscillators or Fourier modes), under a very weak noise, always converges to a NUIS with $\rho_{n\geq 3}=0$; and, $|\rho_{2}|=Q^{*}$. In other words, the system adopts the minimum value of $|\rho_{2}|$ among all allowed by Eq. (27).

The state $Q=1$ ($R=0$) —the last NUIS to destabilize— is 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 $\dot{\delta\rho_{1}}\propto a\,\delta\rho_{1}-a^{*}\,\delta\rho_{1}^{*}$ what yields an additional zero eigenvalue corresponding to the direction $\mathrm{Im}(\delta\rho_{1})=0$.

IV.3 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 $\epsilon=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 $\epsilon$ in which the second-order approximation is accurate. For $c_{2}=1$ the approximate bifurcation lines are accurate up to $\epsilon\approx 0.05$, while this range is certainly smaller for $c_{2}=3$.

For general $c_{1}$, $c_{2}$ values, the prefactor $(\epsilon\eta)^{n}$ appearing for first ($n=1$) and second ($n=2$) orders suggests to extrapolate the relative smallness of $\epsilon$ to other $c_{1}$, and $c_{2}$ values. Thus, if in Fig. 3(b) accuracy is good up to $\epsilon\eta\approx 0.05\eta$, and $\eta\approx 2$, we propose

[TABLE]

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 $(\epsilon\eta)^{3}$, see Sec. V.

IV.4 Quasiperiodic partial synchronization

The phenomenon of QPS was originally reported in the MF-CGLE Nakagawa and Kuramoto (1995) 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 van Vreeswijk (1996). 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(\bar{A})$. 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 $\epsilon$ 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 non-unstable 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 ($|\rho_{n}|\to 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 second-order reduction, since the QPS attractor—found at moderate $\epsilon$ values, see Fig. 2(a) and (b)—is not reproduced.

IV.5 Clustering

Clustering is a much studied phenomenon in oscillator ensembles Golomb et al. (1992). 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 Hakim and Rappel (1992); Nakagawa and Kuramoto (1993); Banaji (2002); Daido and Nakanishi (2006, 2007); Ku et al. (2015); Kemeth et al. (2019), specifically for strong coupling ($\epsilon\approx 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 $\theta^{A}$, and the remaining $(1-p)$ fraction is in the B-state $\theta^{B}\neq\theta^{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 $\Delta=\theta^{A}-\theta^{B}$ and imbalance $p$ corresponding to actual cluster solutions for three different $c_{1}$ values with fixed values of $c_{2}$ and $\epsilon$. 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 two-cluster solutions colliding with FS ($\Delta=0$). In consequence there is a reconnection of the two-cluster solutions. In Fig. 5 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 precise, 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 Hansel et al. (1993) —a stable heteroclinic connection between two configurations of $\Delta=\theta^{A}-\theta^{B}$— is not possible.

IV.6 Finite population, $N=4$

This work focuses on the behavior of the MF-CGLE in the thermodynamic limit ($N\to\infty$). 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 Kemeth et al. (2018, 2019). 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$ Ashwin and Swift (1992), equivalent to the NUISs for $N=\infty$.

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 Ashwin et al. (2008). 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 $\theta_{j}\leq\theta_{j+1}$. (We assume here $\theta_{j}\in[0,2\pi)$ to avoid artificial degeneracies.) The set of phases $\{\theta_{1},\theta_{2},\theta_{3},\theta_{4}\}$, may take several invariant configurations. Apart from the trivial FS state $\{a,a,a,a\}$, there exists a continuum of “NUIS-like” $\mathbf{Z}_{2}$-symmetric states with $\{a,a+b,a+\pi,a+b+\pi\}$, where $b\in(0,\pi/2)\cup(\pi/2,\pi)$. In the limit $b\to\pi/2$ the NUIS becomes the splay state (the analogous of UIS). In addition, in the limits $b\to 0$ and $b\to\pi$ 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\neq a+\pi$). 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 Kemeth et al. (2019). 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, $\mathbf{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 Ashwin and Swift (1992). 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 $\epsilon$) 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. Considering more oscillators will increase the number of cluster solutions, see Banaji (2002), 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(\epsilon^{3})$ correction to Eq. (15):

[TABLE]

This expression depends on the third Kuramoto-Daido order parameter $Z_{3}\equiv Pe^{i\Xi}=N^{-1}\sum_{j}e^{i3\theta_{j}}$. The dependence of constants $\{C_{j},\gamma_{j}\}_{j=1,\ldots,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,\ldots,8$. For illustration we expand a couple of these four-body contributions:

[TABLE]

There are several qualitative features in Eq. (29) that deserve to be pointed out:

The overall $O(\epsilon^{3})$ contribution is not proportional to $\eta^{3}$ —though some terms indeed are— in contrast to $O(\epsilon)$ and $O(\epsilon^{2})$, which are proportional to $\eta$ and $\eta^{2}$, respectively. 2. 2.

From Eqs. (15) and (29) we can expect that truncation of the power series to order $\epsilon^{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\leq n$ appear. 3. 3.

The last two terms in Eq. (29) are somewhat unexpected, (nonetheless see Ashwin and Rodrigues (2016)), since they depend on the mean fields $Z_{1}$ and $Z_{2}$, but not on $\theta_{j}$ itself. They are hence irrelevant concerning synchronization boundaries. 4. 4.

As occurs with the $O(\epsilon^{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 $\epsilon$. Accordingly, it holds that $D+\sum_{j}C_{j}\sin\gamma_{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 $\epsilon$.

From Eqs.(15) and (29) we can derive the stability boundary of FS, NUIS (for UIS just set $Q=0$) obtaining:

[TABLE]

[TABLE]

In Fig. 7 we depict (a) $\epsilon_{0}$, and (b) $\epsilon_{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 $\epsilon=0$.

VI Discussion

VI.1 Alternative phase reduction(s)

Our phase reduction is a genuine power series in the small parameter $\epsilon$. Another strategy to analyze (1) is to absorb the $\epsilon A_{j}$ term prior to the phase reduction. Specifically, setting $t^{\prime}=(1-\epsilon)t$, and

[TABLE]

we get

[TABLE]

where $B_{j}=A_{j}\exp(i\epsilon c_{1}t)/\sqrt{1-\epsilon}$. Applying our phase-reduction method to (33) we obtain an alternative phase reduction in powers of $\kappa$ (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 $\kappa^{n}$, superior to that of (1) at order $\epsilon^{n}$? Certainly, phase reductions at order $\epsilon^{n}$ and $\kappa^{n}$ are not equivalent since $\kappa=\epsilon+\epsilon^{2}+\epsilon^{3}+\cdots$. Any truncation at order $\kappa^{n}$ involves all powers of $\epsilon$. 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 in $\epsilon$ such that the exact critical value $\epsilon_{*}$ (the asterisk denotes an arbitrary state: UIS, FS, …) satisfies:

[TABLE]

The coefficients $a_{n}$ depend on the specific instability we are considering.

Phase reduction of (1) up to order $n$ results in a truncation of (34) to order $n-1$. For instance, the second-order phase reduction of (1) yields the linear relation [recall Eqs. (19) or (24)]:

[TABLE]

At the same order, the phase reduction of (33) results in an analogous expression

[TABLE]

Given that $\kappa=\epsilon+O(\epsilon^{2})$, consistency with (34) determines $a^{\prime}_{1}=a_{1}$. Thus the bifurcation locus estimated from the phase reduction of (33) satisfies (in coordinate $\epsilon$) $a_{1}\epsilon_{*}/(1-\epsilon_{*})+(1+c_{1}c_{2})=0$, which is slightly different from (35). Analogous reasoning permits us to obtain the bifurcation lines for the third-order phase reduction of (33) from Eqs. (30) and (31).

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}\approx\bar{A}$), while (33) works better close to incoherence ($\bar{A}\approx\bar{B}\approx 0$). Finally, note that in addition to (1) and (33), there exists a continuum of alternative, intermediate formulations, in which $\epsilon A_{j}$ is only partly absorbed by a coordinate transformation.

VI.2 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 $\epsilon$ couples phases with the nearest-neighbor phases. At order $\epsilon^{2}$, second nearest neighbors come also into play Kralemann et al. (2011); *kralemann14, and progressively more distant nodes participate in the phase dynamics at higher orders. Also, the case of non-locally coupled Stuart-Landau oscillators Kuramoto (1997) 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: $\nabla^{4}\theta$, $(\nabla^{2}\theta)^{2}$, $(\nabla\theta)^{2}\nabla^{2}\theta$, $\nabla\theta\nabla^{3}\theta$  Kuramoto (1984).

Concerning the oscillator dynamics, the phase reduction carried out here can be easily applied to planar oscillators with polar symmetry ($\lambda-\omega$ systems). In the latter case, analogously to (7), the isochrons satisfy $\theta=\varphi+\chi(r)$ Winfree (1980). Even if function $\chi$ 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.

VI.3 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 $\lambda-\omega$ 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 Kuramoto (1984). This calculation is, however, much more lengthy than the one presented in Sec. III. Not surprisingly, obtaining the order $\epsilon^{3}$ with the averaging approach Kuramoto (1984) 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 $\dot{r}_{j}=0$. This procedure was followed in Pikovsky and Rosenau (2006); Matheny et al. (2019), 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 $\epsilon$, and the result differs from the one obtained following Pikovsky and Rosenau (2006); Matheny et al. (2019). In passing, we mention that instead of assuming $r_{j}=r_{j}(\theta_{1},\theta_{2},\ldots,\theta_{N})$, once Eqs. (10) are derived, the two-timing approximation, such that the $\theta_{j}$ depend only on a slow time $\tau=\epsilon t$, can also be applied.

In contrast to our work, Refs. Kori et al. (2014); Ashwin and Rodrigues (2016) 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 Ashwin and Rodrigues (2016) or Kori et al. (2014) may be interesting.

VI.4 Towards a minimal phase model of pure collective chaos

Pure collective chaos has been found in several phase models with heterogeneity Bick et al. (2018) or delay Pazó and Montbrió (2016); *devalle18. 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 $\epsilon$ values, and eventually, will exhibit collective chaos.

As pairwise interactions through higher harmonics, like $Q\sin(\Phi-2\theta_{j})=N^{-1}\sum_{k}\sin[2(\theta_{k}-\theta_{j})]$, do not show up in the phase reduction of the MF-CGLE 111Otherwise, the frequency of (N)UIS would depend nonlinearly on $\epsilon$, in disagreement with the MF-CGLE., 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 Bick et al. (2011) and Bick et al. (2016), respectively. However, so far, collective chaos remains elusive in populations of higher-order pairwise interacting identical phase oscillators Clusella et al. (2016). 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 multi-body interactions needed in the phase reduction? Is this chaos connected to collective chaos in the thermodynamic limit, as in Bick et al. (2018)?

VI.5 Conclusions

Multi-body interactions are an unavoidable consequence of phase reduction, but save for a few works Tanaka and Aoyagi (2011); Bick et al. (2016); Bick (2018); Matheny et al. (2019); Skardal and Arenas , 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 third-order phase reductions of the MF-CGLE. In our view, higher-order phase reductions promise to be crucial for our understanding of collective chaos and other exotic phenomena Matheny et al. (2019). Moreover, analytic higher-order phase reductions may also serve as test beds for numerical phase reductions recently implemented Rosenblum and Pikovsky (2019). For these reasons, we regard phase reduction beyond the first order as an exciting battleground of nonlinear dynamics.

Acknowledgements.

We acknowledge support by MINECO (Spain) under Project No. FIS2016-74957-P. IL acknowledges support by Universidad de Cantabria and Government of Cantabria under the Concepción Arenal programme.

Appendix I: Clustering

Our model (15) in a more convenient form (recall that in the rotating frame $\Omega=0$) reads:

[TABLE]

with

[TABLE]

Note that $G_{2}$ is an odd function.

Let us write first the evolution equation for cluster-$A$ phase $\theta^{A}$, defining the phase difference $\Delta=\theta^{A}-\theta^{B}$:

[TABLE]

the equivalent equation for the $B$-cluster is obtained with the substitution $\Delta\to-\Delta$ and $p\to(1-p)$. The evolution of $\Delta(t)$ obeys

[TABLE]

Setting $\dot{\Delta}=0$ we obtain a quadratic equation in $p$ that can be solved explicitly. We depict $p(\Delta)$ in Fig. 5 for selected parameter values. Note the symmetry of the curves because of the invariance under $(\Delta,p)\leftrightarrow(-\Delta,1-p)$. There are $\Delta$ values for which $p$ is out of the range $(0,1)$, indicating no two-cluster states with those particular $\Delta$ values exist. Conversely, different values of $\Delta$ may be consistent with the same $p$ value, indicating the coexistence of multiple two-cluster solutions with the same sizes.

VI.5.1 Stability

First of all, note that, one zero eigenvalue is always present due to the global phase shift invariance of the model, $\theta_{j}\to\theta_{j}+\mathrm{const.}$, and we ignore it hereafter. For the analysis that follows it is simpler to assume the thermodynamic limit (eigenvalues are the same for finite $N$, but the calculation is more convoluted.) As already known from previous studies Kori and Kuramoto (2001), perturbations on a two-cluster solution can be decomposed in three orthogonal modes. Two of them are the disintegration of each respective cluster, and the third one is the unlocking of the two clusters. We denote $\lambda_{A}$, $\lambda_{B}$ and $\lambda_{L}$ the corresponding eigenvalues. For the stability of the $A$-cluster, we need to evaluate if one oscillator in the neighborhood of this cluster decays to it or departs (i.e., “evaporates”). The eigenvalue $\lambda_{A}$ is simply obtained linearizing around the state. The result is:

[TABLE]

The eigenvalue $\lambda_{B}$ is obtained from $\lambda_{A}$ after the substitution $p\to(1-p)$ and $\Delta\to-\Delta$, and viceversa. Finally, the locking between the clusters is controlled by the eigenvalue obtained linearizing (Appendix I: Clustering):

[TABLE]

Stability requires $\lambda_{A},\lambda_{B},\lambda_{L}<0$. For small $\epsilon$ we summarized our findings in the main text, distinguishing three different regions corresponding to the three panels of Fig. 5. As said in the main text we found stable clusters in the first region (FS unstable), e.g. $\epsilon=0.1$, $c_{1}=-9$, $c_{2}=2$. However, for these parameters the condition (28) does not hold, and in fact the cluster solution destabilized when we implemented it in the MF-CGLE.

VI.5.2 No slow switching

With unstable two-cluster states, the system might still exhibit one nontrivial phenomenon called slow switching Hansel et al. (1993). In this phenomenon, the clusters switch between two different $\Delta$ values with identical $p$ value. The explanation for this behavior is a stable heteroclinic connection between the pair of two-cluster states that causes the system to switch for ever between them with increasing residence times Hansel et al. (1993); Kori and Kuramoto (2001). In practice Clusella et al. (2016), switching either terminates in one of the unstable two-cluster states (due to round-off errors), or it achieves a constant periodic switching (due to small noise). According to Kori and Kuramoto (2001), slow switching requires the coexistence of three two-cluster states $\Delta^{\prime}$, $\Delta^{\prime\prime}$, and $\Delta^{\prime\prime\prime}$ with an identical $p$ value, such that $0<\Delta^{\prime}<\Delta^{\prime\prime\prime}<\Delta^{\prime\prime}<2\pi$, and $\lambda_{L}<0$ for $\Delta^{\prime}$ and $\Delta^{\prime\prime}$, while $\lambda_{L}>0$ for $\Delta^{\prime\prime\prime}$. As may be seen in Fig. 5, finding parameter values with three solutions for $\Delta$ at the same $p$ is aready difficult —only for the green line in Fig. 5(a) do such $p$ values exist. In addition, the condition for the eigenvalues is even more stringent: e.g., in Fig. 5(a) the three points share the stability of $\lambda_{A}$ and $\lambda_{B}$ making the heteroclinic connection impossible.

Appendix II: Constants $C_{j}$, $\gamma_{j}$ and $D$ in Eq. (29)

The dependences of constants $\{C_{j},\gamma_{j}\}_{j=1,\dots,9}$ on $c_{1}$ and $c2$ is tabulated as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

Additionally,

[TABLE]