A minimal model of self-consistent partial synchrony

Many physical systems can be represented as networks of oscillators, different examples ranging from the mammalian brain [1], to power grids [2], out-of-equilibrium chemical reactions [3], spin-torque nanoscale oscillators [4–6], gene-controlled clocks in bacteria [7], and so on. A large number of books, chapters, and reviews devoted to the topic testify to the importance of this subject [8–17].

A general theory of oscillatory ensembles has not yet been developed. Indeed, such a theory requires taking into account many different features, such as the structure of the single units and their heterogeneity, as well as the topology and properties of the connections. Even in the simple context of globally coupled identical phase oscillators, it is not generally known what a kind of stationary regimes are to be expected. Roughly speaking, they can be classified by referring to the distribution of phases in the limit of a large number of oscillators (the so-called thermodynamic limit). If the mutual interaction leads to phase attraction (at least below a certain distance), the distribution of phases converges to a set of Dirac δ's which correspond to different clusters. Full synchrony is the extreme case, where all oscillators converge to the same trajectory, i.e. to a single cluster. In the presence of a mutual repulsion, a smooth phase distribution is observed instead. The splay state (or, equivalently, the asynchronous regime) is a prototypical example, characterized by a flat distribution of the phases and absence of a collective mode. Finally, one can encounter chimeras, where a big cluster coexists with a group of non-synchronized units [18]. Interestingly, in this state, the frequency of the cluster elements differs from the frequencies of asynchronous ones.

Most of the efforts have been devoted to the study of clustered [19] and chimera states [20, 21] and much less to the identification and analysis of regimes characterized by a smooth but non-uniform distribution of phases. Such regimes, typically characterized by a periodic collective evolution, are herein referred to as self-consistent partial synchrony (SCPS). The simplest form of SCPS is a 'rigid' rotation of the distribution, i.e. a regime where the instantaneous frequency of the oscillators coincides with that of the collective mode. Such a regime can emerge if the coupling strength vanishes for some finite value of the order parameter [22–24]. In a less trivial form of SCPS (of primary interest here) the (average) frequency of the single units and that of the mean field differ from each other. Moreover, the two frequencies are generally mutually incommensurate, i.e. no locking phenomena are observed, when a control parameter is continuously varied. Thus, the microscopic dynamics is quasiperiodic. Examples of such dynamics are: integrate-and-fire (IF) neurons interacting through finite-width pulses (the so-called α-functions) [25]; nonlinearly coupled Stuart–Landau systems [26]; a Kuramoto-like model obtained via phase reduction from the above mentioned Stuart–Landau ensemble. A variant of quasiperiodic partially synchronous dynamics has been detected in models beyond phase approximation, i.e. in globally coupled Hindmarsh–Rose neurons and Stuart–Landau oscillators; here the macroscopic and microscopic frequencies are equal only on average, but the motion of oscillators is additionally modulated by a generally incommensurate frequency [26, 27].

In the weak-coupling limit oscillators are effectively described by a Kuramoto–Daido phase model [28–31] with a suitable coupling function $G({\rm{\Delta }}\phi )$, where ${\rm{\Delta }}\phi$ is the phase difference between any two interacting oscillators. In the standard, widely used, Kuramoto–Sakaguchi model [9, 32, 33] the coupling function G is assumed to be perfectly sinusoidal. In the last years it has become increasingly clear that this is quite a special case: e.g., multiple clusters in this setup are not possible [17, 34]. A much richer dynamics, including formation of clusters [19] and of heteroclinic cycles (HCs) [35], is observed as soon as just one additional harmonic is added to G.

In this paper we further illustrate the richness of the Kuramoto–Daido model, by showing that SCPS spontaneously emerges in a minimal extension of the Kuramoto setup to a biharmonic model, where G is composed of just two harmonics. To further explore the ubiquity of SCPS, we study its emergence in an ensemble of linearly coupled Rayleigh oscillators. They are two-dimensional limit-cycle oscillators; performing numerically the reduction to a Kuramoto–Daido phase model, we reconstruct the coupling function which turns out to contain a few harmonics. The Kuramoto–Daido setup is shown to reproduce the dynamics of the original system.

The simplicity of the biharmonic model allows for a detailed analysis of SCPS, which can be seen as a stationary solution of a continuity equation in a suitably rotating frame. Accordingly, the phase-distribution can be accurately determined from the emergence of SCPS out of the splay state—through a Hopf bifurcation—to its collapse onto full synchrony as in [24, 36]. The additional stability analysis confirms the numerical evidence of the onset of an instability and allows identifying the unstable direction.

As SCPS in the biharmonic model coexists with two-cluster states, we revisit their stability properties, to understand under which conditions trajectories converge towards an HC [35, 37]. We find that more harmonics are needed to ensure that two-cluster states and the corresponding HCs are both unstable. Moreover, we find that HCs can be viewed at as a kind of quasiperiodic partial synchrony.

The paper is organized as follows. In section 2 the model is briefly introduced and the conditions for the stability of the fully synchronous and asynchronous regimes are recalled. The corresponding phase-diagram is thereby presented for a fixed amplitude of the second harmonic. In section 3, we analyze the occurrence of SCPS, show how it can be treated and finally develop the formalism needed to perform the stability analysis. Section 4 is devoted to a discussion of two-cluster states and HC analyzed in [35, 37]. Here, after briefly recalling some known properties, we present a general analysis of the stability properties of two-cluster states, in the perspective of shedding light on the general conditions under which such states can be effectively unstable (this is, for instance the case of the LIF model proposed in [25]). The theoretical predictions are then extensively tested in section 5. There we find that the mean-field frequency is not necessarily smaller than the frequency of the splay state as in the LIF model. Additionally we discover that SCPS can lose stability (through a Hopf bifurcation) and thereby lead to a collapse onto HCs. Actually, in order to avoid spurious effects due to finite computer accuracy, a minimal heterogeneity is added which makes the oscillators slightly different from one another. As a result, we find a bistable regime, where SCPS coexists with stable HC. In section 6, we discuss SCPS in two other setups that can be effectively described by a suitable Kuramoto–Daido model: a set of LIF neurons, and an ensemble of Rayleigh oscillators. A Kuramoto–Daido description of the former model was already presented in [38] where the validity of the approximation was investigated. Here we analyze two-cluster states verifying their instability. As for the Rayleigh oscillators, we reduce their description to a Kuramoto–Daido phase model and show that, at variance with the other setups formerly considered, here instability of the splay state is due to harmonics higher than the second one. The main results and the still open problems are summarized in the last section.

We hereby consider the Kuramoto–Daido type model of identical all-to-all coupled phase oscillators. Performing a transformation to the co-rotating coordinate frame we set the frequency to zero, so that the model reads

$$\begin{eqnarray}&&{\dot{\phi }}_{i}=\displaystyle \frac{1}{N}\displaystyle \sum _{j}G({\phi }_{j}-{\phi }_{i}),\end{eqnarray} \tag{ 1 }$$

where the coupling constant has been eliminated by performing a suitable rescaling of time. In most of the paper we consider the biharmonic coupling function

$$\begin{eqnarray}&&G(\phi )=\sin (\phi +{\gamma }_{1})+a\sin (2\phi +{\gamma }_{2}).\end{eqnarray} \tag{ 2 }$$

A standard way to classify the configurations of an ensemble of oscillators is via the set of complex order parameters

$$\begin{eqnarray}&&{Z}_{m}={R}_{m}{{\rm{e}}}^{{\rm{i}}{\beta }_{m}}=\displaystyle \frac{1}{N}\displaystyle \sum _{j}{{\rm{e}}}^{{\rm{i}}{m}{\phi }_{j}},\end{eqnarray} \tag{ 3 }$$

where R₁ is the famous Kuramoto order parameter [9, 32], which is equal to 1 in the case of full synchrony and is equal to 0 in splay states.

By making use of this definition, equation (1) can be rewritten as

$$\begin{eqnarray}&&{\dot{\phi }}_{k}={R}_{1}\sin ({\beta }_{1}-{\phi }_{k}+{\gamma }_{1})+{{aR}}_{2}\sin ({\beta }_{2}-2{\phi }_{k}+{\gamma }_{2}).\end{eqnarray} \tag{ 4 }$$

This model was already considered in [35, 37] with an emphasis on analysis of clustered states (see the next section) and has recently attracted a growing interest, especially in the presence of a distribution of oscillator frequencies [39, 40].

For a = 0 equations (1) and (2) reduce to the famous Kuramoto–Sakaguchi model. In this case it is known that the fully synchronous solution ${\phi }_{k}=\phi$ is stable if and only if $| {\gamma }_{1}| \lt \pi /2$, while the stability condition of the splay state is exactly opposite: asynchrony is stable for $\pi /2\lt | {\gamma }_{1}| \lt 3\pi /2$ and unstable otherwise. An SCPS solution can arise only at the border between stability and instability of the two regimes.

This pointwise region can be made structurally stable by assuming that ${\gamma }_{1}$ is a function of the order parameter R₁ and of the coupling strength $\varepsilon$ [24, 36] (such phase model can be obtained in the process of phase reduction of a system of nonlinearly coupled Stuart–Landau oscillators [26]). The resulting regime was called self-organized quasiperiodic dynamics.

In this paper we show that adding a second harmonic is yet a simpler way to generate SCPS. In the presence of a non-zero a, the stability of the fully synchronous state is determined by the condition

$$\begin{eqnarray}&&G^{\prime} (0)=-\cos {\gamma }_{1}-2a\cos {\gamma }_{2}\lt 0.\end{eqnarray} \tag{ 5 }$$

The marginal stability line is composed of two sinusoidal curves ${{\rm{\Gamma }}}_{+}$, shown in the stability diagram in figure 1. Altogether, the synchronous state is stable in the two yellow and green regions of the diagram.

Zoom In Zoom Out Reset image size

On the other hand, the stability of the asynchronous state can be assessed by linearizing the continuity equation for the probability density of the phases. As shown in [38], the evolution is diagonal in Fourier space and the leading eivengalue ${\delta }_{1}$ is that one characterizing the stability of the first Fourier mode of the perturbation. In the present setup

$$\begin{eqnarray}&&{\delta }_{1}=\pi [\cos {\gamma }_{1}+{\rm{i}}\sin {\gamma }_{1}],\end{eqnarray} \tag{ 6 }$$

so that ${\delta }_{1}$ depends only on ${\gamma }_{1}$. As a result, the splay state is stable in a rectangular region where $\cos {\gamma }_{1}\lt 0$ (see the cyan and green areas in figure 1). In the green areas, both the splay and the synchronous states are simultaneously stable, while in the white areas both of them are unstable.

Altogether, the diagram in figure 1 has a reflectional symmetry with respect to both the horizontal and vertical semi axes. Physically, there is only one symmetry: if $\phi \to -\phi$, the same dynamics is found upon mapping ${\gamma }_{1}\to 2\pi -{\gamma }_{1}$ and ${\gamma }_{2}\to 2\pi -{\gamma }_{2}$. The additional symmetry is, therefore, only accidental.

In the thermodynamic limit one can investigate the various regimes by studying the evolution of the probability density $P(\phi ,t)$ of oscillators with phase ϕ at time t. It satisfies the continuity equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial P}{\partial t}=-\displaystyle \frac{\partial }{\partial \phi }\left[\left(\int {\rm{d}}\psi G(\psi -\phi )P(\psi ,t)\right)P(\phi ,t)\right].\end{eqnarray} \tag{ 7 }$$

The splay state corresponds to a flat and constant density $P=1/(2\pi )$. As already illustrated in [38], SCPS typically manifests itself as a rotating non-uniform density, which emerges past a Hopf bifurcation. It is therefore convenient to perform a change of variables, introducing the angle $\theta =\phi -{\rm{\Omega }}t$ and $Q(\theta ,t)=P(\phi ,t)$. The corresponding evolution equation takes the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial Q}{\partial t}=\displaystyle \frac{\partial }{\partial \theta }\left[\left({\rm{\Omega }}-\int {\rm{d}}\psi G(\psi -\theta )Q(\psi ,t)\right)Q(\theta ,t)\right].\end{eqnarray} \tag{ 8 }$$

For a suitably chosen Ω there exists a stationary time-independent solution ${Q}_{0}(\theta )$, which satisfies the equation

$$\begin{eqnarray}&&\left[{\rm{\Omega }}-\int {\rm{d}}\psi G(\psi -\theta ){Q}_{0}(\psi )\right]{Q}_{0}(\theta )=\eta ,\end{eqnarray} \tag{ 9 }$$

where η is the probability flux. Notice that $2\pi \eta$ can be interpreted as the average microscopic frequency of the single oscillators in the moving frame. ${Q}_{0}(\psi )$ can be expanded in Fourier modes, (the harmonics coincide with the generalized order parameters defined in equation (3) for a finite ensemble of oscillators)

$$\begin{eqnarray}&&{Q}_{0}(\psi )={\displaystyle \sum }_{m=-\infty }^{+\infty }{Z}_{m}{{\rm{e}}}^{-{\rm{i}}{m}\psi },\end{eqnarray} \tag{ 10 }$$

where ${Z}_{-m}={Z}_{n}^{* }$. The stationary solution can be obtained by solving equation (9)

$$\begin{eqnarray}&&{Q}_{0}(\theta )=\displaystyle \frac{\eta }{\pi {\rm{i}}\,[{Z}_{1}{{\rm{e}}}^{-{\rm{i}}(\theta -{\gamma }_{1})}+{{aZ}}_{2}{{\rm{e}}}^{-{\rm{i}}(2\theta -{\gamma }_{2})}-{\rm{c}}.{\rm{c}}.]+{\rm{\Omega }}}.\end{eqnarray} \tag{ 11 }$$

Since the phase of the solution is arbitrary, we are free to fix it by imposing that Z₁ is real. By considering that η can be determined by imposing a normalization on Q₀, the above equation contains four unknowns: ${\rm{\Omega }}$, Z₁, and Z₂ (the last variable is complex). They can be determined self-consistently by imposing

$$\begin{eqnarray}&&{Z}_{k}=\displaystyle \frac{1}{2\pi }\int {\rm{d}}\psi {{\rm{e}}}^{{\rm{i}}{k}\psi }{Q}_{0}(\psi )\end{eqnarray} \tag{ 12 }$$

for k = 1, 2. The solution can be found by searching for a fixed point in a four-dimensional space.

3.1. Stability analysis

Consider an infinitesimal perturbation $q(\theta ,t)$ of ${Q}_{0}(\theta )$ and linearize equation (8), making use of equation (9). One obtains

$$\begin{eqnarray}&&\displaystyle \frac{\partial q(\theta ,t)}{\partial t}=\displaystyle \frac{\partial }{\partial \theta }\left[\eta \displaystyle \frac{q(\theta ,t)}{{Q}_{0}(\theta )}-{Q}_{0}(\theta )\int {\rm{d}}\psi G(\psi -\theta )q(\psi ,t)\right].\end{eqnarray} \tag{ 13 }$$

By expanding the perturbation in Fourier series

$$\begin{eqnarray*}&&{\hat{q}}_{k}=\displaystyle \frac{1}{2\pi }\int {\rm{d}}\psi q(\psi ,t){{\rm{e}}}^{{\rm{i}}{k}\psi }\end{eqnarray*}$$

one can rewrite the integral in the previous equation as

$$\begin{eqnarray*}&&\int {\rm{d}}\psi G(\psi -\theta )q(\psi ,t)=-\pi {\rm{i}}({\hat{q}}_{1}{{\rm{e}}}^{-(\theta -{\gamma }_{1}){\rm{i}}}+{\hat{q}}_{2}a{{\rm{e}}}^{-(2\theta -{\gamma }_{2}){\rm{i}}}-{\rm{c}}.{\rm{c}}.)\equiv B(\theta )\end{eqnarray*}$$

so that

$$\begin{eqnarray}&&\displaystyle \frac{\partial q(\theta ,t)}{\partial t}=\displaystyle \frac{\partial }{\partial \theta }\left[\eta \displaystyle \frac{q(\theta ,t)}{{Q}_{0}(\theta )}-{Q}_{0}(\theta )B(\theta )\right].\end{eqnarray} \tag{ 14 }$$

This equation can be expressed as a Fourier series. With the help of equation (11), it is found

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}{\hat{q}}_{m}}{{\rm{d}}{t}} & = & m\{-{\rm{i}}{\rm{\Omega }}{\hat{q}}_{m}+\pi {{\rm{e}}}^{{\rm{i}}{\gamma }_{1}}[{Z}_{1}{\hat{q}}_{m-1}+{Z}_{m-1}{\hat{q}}_{1}]-\pi {{\rm{e}}}^{-{\rm{i}}{\gamma }_{1}}[{Z}_{-1}{\hat{q}}_{m+1}+{Z}_{m+1}{\hat{q}}_{-1}]\\ & & +a\pi {{\rm{e}}}^{{\rm{i}}{\gamma }_{2}}[{Z}_{2}{\hat{q}}_{m-2}+{Z}_{m-2}{\hat{q}}_{2}]-a\pi {{\rm{e}}}^{-{\rm{i}}{\gamma }_{2}}[{Z}_{-2}{\hat{q}}_{m+2}+{Z}_{m+2}{\hat{q}}_{-2}]\}.\end{eqnarray} \tag{ 15 }$$

The mth Fourier mode of the perturbation is coupled with the four nearest neighbor modes ($m-2$, $m-1$, $m+1$, and $m+2$) as well as with the first two modes; the latter coupling is mediated by the amplitude of higher components of the stationary solution ${\text{}}{Q}_{0}$. In other words, the corresponding matrix is sparse: it is pentadiagonal with two full rows. As each derivative is multiplied by m, ${\dot{\hat{q}}}_{0}=0$. This is a straightforward consequence of the conservation of the total probability; therefore ${\hat{q}}_{0}$ can be eliminated as it does not contribute to the eigenvalues. By further looking at the evolution equation for ${\hat{q}}_{2}$, we see that it involves ${\hat{q}}_{-1}$, so that the negative modes must be included as well. Since ${\hat{q}}_{-m}={\hat{q}}_{m}^{* }$, it is convenient to separate ${\hat{q}}_{m}$ into real and imaginary part (${\hat{q}}_{m}={u}_{m}+{{\rm{i}}{v}}_{m}$), so that we can exploit the relationships ${u}_{-m}={u}_{m}$, ${v}_{-m}=-{v}_{m}$ and thereby get rid of the negative m components. The relevant eigenvalues ${\mu }_{{\rm{R}}}+{\rm{i}}{\mu }_{I}$ of the resulting matrix can be then computed by considering a sufficiently large number of Fourier modes.

Clustered states represent another class of stationary solutions. In the biharmonic model such states have been already investigated in [35, 37]. Here below we summarize those results that are necessary to proceed ahead with our general considerations.

A two-cluster state consists of two families of oscillators with phases α and ψ, respectively. Both for the sake of simplicity and since they are the most widely observed, here, we mostly focus on the symmetric case of equal-size clusters. The phases of the two families follow the differential equations

$$\begin{eqnarray}&&\dot{\alpha }=(G(0)+G(\psi -\alpha ))/2,\quad \dot{\psi }=(G(\alpha -\psi )+G(0))/2.\end{eqnarray} \tag{ 16 }$$

The separation $\delta =\alpha -\psi$ satisfies the equation

$$\begin{eqnarray}&&\dot{\delta }=(G(-\delta )-G(\delta ))/2=-{G}_{{\rm{A}}}(\delta ),\end{eqnarray} \tag{ 17 }$$

where G_A is the anti-symmetric component of G. Two-cluster solutions are identified by the zeros of G_A; $\delta =0$ is always a solution which corresponds to a single cluster (vanishing distance between the two clusters).

In the biharmonic model, the symmetric and anti-symmetric component of the coupling function are

$$\begin{eqnarray}&&{G}_{{\rm{A}}}(\delta )=\sin \delta (\cos {\gamma }_{1}+2a\cos {\gamma }_{2}\cos \delta ),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{G}_{{\rm{S}}}(\delta )=\sin {\gamma }_{1}\cos \delta +a\sin {\gamma }_{2}\cos 2\delta .\end{eqnarray} \tag{ 19 }$$

There are various solutions of the equation ${G}_{{\rm{A}}}(\delta )=0$ besides ${\delta }_{0}=0$: ${\delta }_{0}=\pi$ corresponds to an antiphase two-cluster state. Two further solutions can be found by setting to zero the expression in parentheses in equation (18); however, these solutions represent only one physically meaningful state. Indeed, given a two-cluster state characterized by a separation ${\delta }_{0}$, the same state can be seen as characterized by a separation $2\pi -{\delta }_{0}$, if the two clusters are exchanged. These states exist only in the parameter region delimited by the curves where

$$\begin{eqnarray}&&\cos {\gamma }_{1}\pm 2a\cos {\gamma }_{2}=0.\end{eqnarray} \tag{ 20 }$$

This equation with a plus sign defines the curve ${{\rm{\Gamma }}}_{+}$ which coincides with the bifurcation line where the synchronous state loses stability, see figure 1. (In fact, $\cos {\delta }_{0}=1$ means ${\delta }_{0}=0$, i.e. the two-cluster solution bifurcates from the fully synchronous one.) The minus sign, instead, corresponds to the curve ${{\rm{\Gamma }}}_{-}$ where the two-cluster state becomes the antiphase one. As a result, non-trivial clustered solutions with $\delta \ne \pi$ exist only in the regions delimited by ${{\rm{\Gamma }}}_{+}$ and ${{\rm{\Gamma }}}_{-}$ (see the two pairs of eyelets e₁ and e₂ in figure 1).

The stability of a two-cluster state is determined by the value of the inter- and intra-cluster exponents. The inter-cluster exponent ${\lambda }_{{\rm{I}}}$ measures the stability against perturbation of the phase-separation between the two clusters. From equation (17) it follows

$$\begin{eqnarray}&&{\lambda }_{{\rm{I}}}=-{G}_{{\rm{A}}}^{\prime }({\delta }_{0}).\end{eqnarray} \tag{ 21 }$$

In the biharmonic model, ${G}_{{\rm{A}}}^{\prime }({\delta }_{0})=\cos {\gamma }_{1}\,\cos {\delta }_{0}+2a\,\,\cos {\gamma }_{2}\,\cos \,2{\delta }_{0}$. The intra-cluster exponents ${\lambda }_{E}^{\pm }$ measure the stability of the width of each cluster. From the equations of motion it is readily found that

$$\begin{eqnarray}&&{\lambda }_{E}^{\pm }=-[{G}_{{\rm{A}}}^{\prime }(0)+{G}_{{\rm{A}}}^{\prime }({\delta }_{0})\pm {G}_{{\rm{S}}}^{\prime }({\delta }_{0})].\end{eqnarray} \tag{ 22 }$$

The solution with the plus (minus) sign refers to the cluster that is lagging behind (leading) by ${\delta }_{0}$. The maximal eigenvalue is therefore

$$\begin{eqnarray}&&{\lambda }_{{\rm{M}}}=-[{G}_{{\rm{A}}}^{\prime }(0)+{G}_{{\rm{A}}}^{\prime }({\delta }_{0})-| {G}_{{\rm{S}}}^{\prime }({\delta }_{0})| ]/2.\end{eqnarray} \tag{ 23 }$$

The inter-cluster exponent of the antiphase solution, ${G}_{{\rm{A}}}^{\prime }(\pi )=-\cos {\gamma }_{1}+2a\cos {\gamma }_{2}$, is negative in between the upper and lower branch of ${{\rm{\Gamma }}}_{-}$ and vanishes along ${{\rm{\Gamma }}}_{-}$, confirming that the clustered solutions bifurcate out of the antiphase state. This region is almost complementary to the stability area of the synchronous state, since $G^{\prime} (\pi )G^{\prime} (0)\lt 0$ in the region where other clusters do not exist. The stability of the antiphase state to intra-cluster perturbations is controlled by (see equation (23))

$$\begin{eqnarray}&&{\lambda }_{{\rm{M}}}=-2a\cos {\gamma }_{2},\end{eqnarray} \tag{ 24 }$$

so that this state is stable within the region delimited by the purple curve in figure 1.

The non-trivial two-cluster state (${\delta }_{0}\ne \pi$) is unstable against inter-cluster perturbations within the eyelets e₂ (because of the one-dimensional nature of the equation for δ, it has opposite stability with respect to that of the synchronous solution), while it is stable within e₁. Its intra-cluster stability can be determined from equation (23). By taking into account that ${G}_{{\rm{A}}}({\delta }_{0})=0$, we find that

$$\begin{eqnarray}&&2{\lambda }_{{\rm{M}}}=-\cos {\gamma }_{1}(1-\cos {\delta }_{0})+| \sin {\delta }_{0}(\sin {\gamma }_{1}+4a\sin {\gamma }_{2}\cos {\delta }_{0})| .\end{eqnarray} \tag{ 25 }$$

Therefore, such a two-cluster state is unstable for $\pi /2\lt {\gamma }_{1}\lt 3\pi /2$. A detailed analysis for a = 0.2 reveals that the solution is unstable also within the eyelet e₁.

This is not yet the end of the story. The opposite sign of the two exponents ${\lambda }_{E}^{\pm }$ implies that, while the width of one cluster decreases, that of the other one diverges. However, as already discussed in [35], once the width of the 'exploding' cluster becomes of order 1, nonlinear effects (not captured by a linear stability analysis) induce a relative phase shift of the two clusters, so that the leading cluster becomes the lagging one: this implies a stability 'exchange'. As a consequence, the long term behavior can be assessed after averaging over the alternating periods of stability and instability. Symmetry reasons imply that the time duration of such two periods are equal to one another, so that the average exponent is, in general

$$\begin{eqnarray}&&{\lambda }_{{\rm{a}}}=-[G{^{\prime} }_{{\rm{A}}}(0)+G{^{\prime} }_{{\rm{A}}}({\delta }_{0})]/2.\end{eqnarray} \tag{ 26 }$$

In the biharmonic model

$$\begin{eqnarray}&&{\lambda }_{{\rm{a}}}=-\cos {\gamma }_{1}(1-\cos {\delta }_{0})/2.\end{eqnarray} \tag{ 27 }$$

In the region of interest, ${\lambda }_{{\rm{a}}}\lt 0$ if $\cos {\gamma }_{1}\gt 0$, i.e the fluctuations of the cluster widths on average decrease, without ever collapsing onto it. This is nothing but an attracting HC. Because of the finite computer accuracy, these oscillations necessarily collapse on the otherwise unstable two-cluster state.

5.1. Microscopic analysis

We start by exploring the parameter region identified by the rectangle in figure 1, which includes the area where highly symmetric synchronous and asynchronous states and two-cluster states are all unstable. (Because of the above mentioned symmetry, the upper eyelet is characterized by an equivalent dynamics.)

The outcome of a direct integration of equations (1) and (2) (starting from an initial condition close to the splay state) is summarized in figure 2: the symbols identify the different asymptotic states, while the curves correspond to the marginal stability lines (determined theoretically, see figure 1). The simulation of equations (1) and (2) was performed for N = 1000; larger ensemble-size have been considered for several points without observing any essential difference. The transient time was $2\times {10}^{4}$ time units.

Zoom In Zoom Out Reset image size

In order to avoid the spurious formation of clusters in the HC states, we made the oscillators slightly heterogeneous. Namely, their frequencies (all equal to zero in equations (1) and (2)) have been taken as uniformly distributed in the interval $[-0.5\times {10}^{-12},0.5\times {10}^{-12}]$. This diversity, which is crucial for the detection of HCs, had no influence on the other dynamical states. It has been checked that variation of the inhomogeneity in the range ${10}^{-10}-{10}^{-14}$ has no essential effect on the parameters of the HC. For an automatic detection of the states, all oscillators characterized by phase differences smaller than 10⁻⁸ have been identified as belonging to the same cluster (this threshold was chosen by trial and error, after a visual inspection of the observed regimes). By monitoring the number of clusters, we have found that their number varies in time only in the case of HCs. This is due to the fact that the two cluster-widths greatly change over time, as illustrated in figure 3(b), where the phase of each oscillator is plotted versus time in a co-rotating frame. There we see that time intervals where the cluster amplitude is of order one alternate with relatively long periods where the width is extremely small, thus yielding spuriously detected clusters. The asymmetry between the behavior of the two clusters (one of them splits into two parts) follows from the non-perfectly equal size of the two clusters. The periodic variation of the cluster-width is reflected in periodic variation of the order parameter R₁ (figure 3(a)). Furthermore, the average frequency of the mean field is larger than that of oscillators, as can be appreciated from the plot of the mean-field phase (figure 3(b)); this difference emerges due to a permanent interchange of the leading and lagging clusters, discussed in the previous section. Thus, the HC state can be interpreted as a special form of quasiperiodic SCPS dynamics, with a smooth but non-stationary distribution. A movie showing the time evolution of HCs is included in the supplementary data (see supplementary movie 1).

Zoom In Zoom Out Reset image size

Furthermore, we noticed that some regimes are approached after a very long transient. This is particularly true close to the stability border of different states. The points that appear as SCPS for ${\gamma }_{2}\approx 4.75$, ${\gamma }_{1}\approx \pi /2$, converge to two-cluster states if the transient time is increased. Next, consider the thin 'belt' of SCPS states close to the theoretical curve, which denotes the border of stability of full synchrony. Computations show that upon increasing the transient time, some points that appear as SCPS states in figure 2 progressively converge towards an HC. Thus, the SCPS in the 'belt' domain is probably a very long transient regime. However, for the points in the main domain of SCPS, e.g. for ${\gamma }_{2}=\pi$, ${\gamma }_{1}=1.5$, the SCPS dynamics seems to be the asymptotic state, at least it survives for 10⁷ time units.

Further states have been occasionally detected (see the cyan symbol in figure 2). In particular three-cluster states are found for ${\gamma }_{2}=4.75$ and $1.58\leqslant {\gamma }_{1}\leqslant 1.66$, where the degree of synchronization may even oscillate in time (three-cluster states had been already observed in [35] for a slightly different amplitude of the second harmonic).

A detailed quantitative analysis has been performed along the line ${\gamma }_{2}=\pi$ in the parameter space. The results are shown in figure 4, where one can recognize three critical points: (i) ${\gamma }_{s}=\pi /2$ signals the loss of stability of the splay state; (ii) ${\gamma }_{f}=\arccos (2a)\approx 1.159$ signals the loss of stability of the fully synchronous state; (iii) ${\gamma }_{p}\approx 1.401$ signals the loss of stability of SCPS.

Zoom In Zoom Out Reset image size

Upon decreasing ${\gamma }_{1}$ from ${\gamma }_{s}$, SCPS is first born through a Hopf bifurcation from the splay state; R₁ and R₂ become strictly larger than zero and a finite mean-field frequency ${\rm{\Omega }}$, which corresponds to the frequency of the Hopf bifurcation, appears discontinuously. Simultaneously, the microscopic frequency ω deviates from zero because of the macroscopic modulation, without revealing any locking with ${\rm{\Omega }}$.

Interestingly, ${\rm{\Omega }}$ is positive and always larger than ω: this is at variance with the scenario observed in [38], where the opposite was found. However, we should also recall that an equivalent scenario is observed in the upper eyelet, once the transformation $\phi \to -\phi$ has been performed. In fact, such a change of variable would lead to the scenario observed in LIF neurons. Notice that for Kuramoto-like model of nonlinearly coupled oscillators [24, 36] both cases (${\rm{\Omega }}\gt \omega$ and ${\rm{\Omega }}\lt \omega$) are possible. Altogether, the existence of a finite difference between microscopic and macroscopic frequencies is a key signature of a quasiperiodic SCPS.

Upon further decreasing ${\gamma }_{1}$, we enter a region where the dynamics converges to HCs which are characterized by a sudden increase of the (average) R₁, while no discontinuity is exhibited by ${R}_{2}({\gamma }_{1})$. When, finally, ${\gamma }_{f}$ is approached, unsurprisingly R₁ and R₂ converge to 1, while both ${\rm{\Omega }}$ and ω converge to the frequency of the fully synchronous state. The nature of the bifurcation is not clear. We briefly comment in the next section, while discussing the stability of SCPS.

The sudden jump observed for ${\gamma }_{1}={\gamma }_{p}$ suggests the possible existence of multistability. Accordingly, we have performed two additional series of simulations, by progressively increasing (decreasing) ${\gamma }_{1}$, and choosing the new initial condition as a slightly perturbed version of the final configuration for the previous value of ${\gamma }_{1}$. As shown in figure 5, a bistability region is indeed found where HCs and SCPS are simultaneously stable.

Zoom In Zoom Out Reset image size

5.2. Macroscopic analysis

For a more detailed characterization of SCPS, we have determined the corresponding probability distribution by solving equation (11), as discussed in section 3. This method allows to obtain the phase distribution even when it is unstable. As it can be seen in figure 4 (see the thin solid lines), the results are consistent with the direct numerical simulations wherever SCPS is stable. Moreover, it is found that SCPS exists also in the interval $[{\gamma }_{f},{\gamma }_{p}]$ and it reconnects to the fully synchronous state. As for the microscopic frequency ω, given by

$$\begin{eqnarray*}&&\omega ={\rm{\Omega }}+2\pi \eta ,\end{eqnarray*}$$

it also agrees with the numerical simulations.

The stability analysis carried out in section 3 allows determining ${\gamma }_{p}$ by solving the eigenvalue problem (15). A typical spectrum is plotted in figure 6. There we see that, with a few exceptions, the eigenvalues tend to align along the imaginary axis. Besides the zero associated to the conservation of the total probability, there exists a second zero eigenvalue which follows from the invariance under phase-translation of the probability density. In panel (b) we see that the real part of the eigenvalues decreases exponentially with their imaginary component. The plateau seen for large ${\mu }_{{\rm{I}}}$ is a consequence of the finite numerical accuracy of the simulations. Interestingly, finite-size effects are practically absent (at least in this parameter region): a larger number of Fourier modes leads to eigenvalues with a larger frequency (imaginary component) and an exponentially small real part.

Zoom In Zoom Out Reset image size

In practice, SCPS is marginally stable, as there is an infinite number of eigenvalues with a practically vanishing real part. The evolution of a uniform distribution initially confined to an interval of size $2\pi -{{\rm{\Delta }}}_{0}$ offers the chance to appreciate the role of the weakly attracting directions. As seen in figure 7, the gap size goes to zero, but it does so in an extremely slow way, namely as $1/\mathrm{ln}t$.

Zoom In Zoom Out Reset image size

In the past, the evidence of a similarly slow convergence was found for the splay state itself. In the context of pulse-coupled oscillators with an analytic velocity field a similar set of exponentially decreasing real parts had been observed [41, 42]. In the context of Kuramoto–Daido models, the strength of the real part is directly proportional to the amplitude of the Fourier component of the coupling function (see equation (C2) in [38]), so that, when the number of Fourier modes is finite, infinitely many strictly marginal directions are present (in the thermodynamic limit). This is reminiscent of the Watanabe–Strogatz theorem [43, 44] which implies that (for a strictly mono-harmonic coupling function) infinitely many directions are not only linearly marginally stable but actually correspond to conservation laws. Our results show that the existence of conservation laws breaks down already when two harmonics are considered. In fact, although exactly marginally stable directions are detected in the analysis of the splay state (here only two Fourier modes are present, so that only two directions can be strictly (un)stable), the same is no longer true for SCPS, as all modes have a weak but finite stability.

In figure 8 we plot the real and imaginary part of the most unstable eigenvalue versus ${\gamma }_{1}$. The data confirms that linear instability occurs below ${\gamma }_{p}$: the bifurcation is of Hopf type. SCPS is maximally unstable around ${\gamma }_{1}\approx 1.33$. By further decreasing ${\gamma }_{1}$, both the real and the imaginary parts decrease to zero, while approaching ${\gamma }_{f}$. From the dependence of the real part, it seems that the scaling behavior is quadratic in the distance from the critical point. The nature of the bifurcation is unclear: simulations close to ${\gamma }_{f}$ are not reliable. It is nevertheless instructive to notice that the weak instability is consistent with the observation of SCPS over extremely long time scales mentioned in the beginning of this section. The correct identification of the critical point is further confirmed in figure 2, where the triangles are the outcome of the stability analysis for three different choices of ${\gamma }_{2}$.

Zoom In Zoom Out Reset image size

So far, we have discussed the occurrence of SCPS in a simple Kuramoto–Daido setup where the coupling function is composed of just two Fourier harmonics. The Kuramoto–Daido representation generally holds in the weak coupling limit; however, the bi-harmonic coupling function may be too simple a model and it is therefore worth comparing with other setups. As recalled in the introduction, SCPS was first observed in pulse-coupled LIF neurons [25]. In [38], it has been shown that the system can be effectively described by a Kuramoto–Daido model with a coupling function containing several non-negligible harmonics. In the first subsection we show that such harmonics contribute to destabilizing the two-cluster states that are, in fact, never observed. In the second subsection we discuss an ensemble of two-dimensional Rayleigh oscillators each described by two variables, showing that SCPS can be generated in this case as well. A phase reduction works also in this latter case, where the higher harmonics play a different role: they are responsible for the destabilization of the splay state, giving rise to more structured cluster states.

6.1. LIF neurons

One of the important open questions in the study of ensembles of phase-oscillators is that of determining a priori whether a given coupling function G gives rise to either smooth distributions, or macroscopic clusters, or both. In this perspective it is instructive to explore the difference between the scenario seen in the biharmonic model and that observed in an ensemble of LIF neurons. In such a case the model can be reduced to a Kuramoto–Daido setup with

$$\begin{eqnarray}&&G(\varphi )={g}_{1}({g}_{2}-1+\varphi ){{\rm{e}}}^{\alpha (1-\varphi )/\nu }+{g}_{3}{{\rm{e}}}^{(1-\varphi )/\nu }-{g}_{4}\end{eqnarray} \tag{ 28 }$$

(see [38] for a precise definition of the various parameters—notice that here we are using a different notation—$\varphi$ instead of ϕ—since in the above equation the phase is normalized between 0 and 1). Correspondingly we change notations in this paragraph. Straightforward calculations reveal that two-cluster states exist also in the above model, the main difference being, however, that now (for $\alpha =6$), the average exponent ${\lambda }_{{\rm{a}}}\gt 0$, so that two-cluster states are effectively unstable, i.e. no any spurious cluster is appearing due to the finite precision of computations. This explains why they have never been seen in numerical simulations.

It is natural to ask whether this holds true for arbitrarily large α, when the SCPS becomes increasingly close to full synchrony. We proceed by expanding G_A around 0 in the limit of large α. One first finds

$$\begin{eqnarray*}&&{G}_{{\rm{A}}}^{\prime }(0)=-\displaystyle \frac{\nu }{a+g\nu }({{\rm{e}}}^{\tau }-1),\end{eqnarray*}$$

which is negative and finite. Comparison with equation (17), implies that the synchronous solution is always unstable. As for the second derivative, the leading order

$$\begin{eqnarray*}&&{G}_{{\rm{A}}}^{\prime\prime }(0)={\alpha }^{3}\displaystyle \frac{{{\rm{e}}}^{\tau }-1}{2\nu (a+g\nu )}\end{eqnarray*}$$

is positive and increasingly large. As a result, on the basis of the first two polynomial terms of ${G}_{{\rm{A}}}(\delta )$, one finds that it vanishes also for

$$\begin{eqnarray}&&{\delta }_{0}=\displaystyle \frac{2{\nu }^{2}}{{\alpha }^{3}}.\end{eqnarray} \tag{ 29 }$$

This phase shift identifies a two-cluster state; its value decreases as the cubic power of the pulse-width (equal to $1/\alpha$). Under this quadratic approximation, the slope of G_A in 0 and ${\delta }_{0}$ are equal and opposite to one another, so that ${\lambda }_{T}=0$. This marginal value requires going one order beyond in the perturbation analysis. The computation of the third derivative shows that it is negative and of order ${\alpha }^{4}$. As a result its contribution to the derivative in ${\delta }_{0}$ is of order ${\alpha }^{4}{\delta }_{0}^{2}\approx 1/{\alpha }^{2}$. This makes the sum of the derivatives in 0 and ${\delta }_{0}$ slightly negative and proves that the two-cluster state is always effectively unstable.

If we recall the definition of G_A, one might be surprised to see that its expansion around 0 contains the second, even, power of φ. One should however remember that the original function G is continuous but not even of class ${{ \mathcal C }}^{(1)};$ so it is not unnatural to expect a discontinuity of some derivative in zero.

In smooth models, in the vicinity of the bifurcation where the synchronized solution loses its stability (e.g., for the biharmonic model) one can write

$$\begin{eqnarray}&&{G}_{{\rm{A}}}(\varphi )=-\varepsilon \varphi +A{\varphi }^{3}.\end{eqnarray} \tag{ 30 }$$

If $A\gt 0$ and $\varepsilon \gt 0$, then there exists a second zero ${\delta }_{0}=\sqrt{\varepsilon /A}$, which corresponds to the clustered solution. Therefore

$$\begin{eqnarray}&&{\lambda }_{T}=2\varepsilon -3A{\varphi }_{0}^{2}=-\varepsilon \lt 0,\end{eqnarray} \tag{ 31 }$$

so that smooth interaction functions necessarily lead to effectively stable two-cluster states (at least with respect to intracluster perturbations).

Even more, equation (27) yields that in the biharmonic model no effectively unstable cluster may exist: only HCs. For these clusters to exist, higher harmonics are needed.

6.2. Rayleigh oscillators

In order to provide further evidence on the ubiquity of SCPS, we present a numerical study of globally coupled identical Rayleigh oscillators⁴ . The equations are

$$\begin{eqnarray}&&{\ddot{x}}_{k}-\zeta (1-{\dot{x}}_{k}^{2}){\dot{x}}_{k}+{\omega }^{2}{x}_{k}=\varepsilon \mathrm{Re}[{{\rm{e}}}^{{\rm{i}}\gamma }(X+{\rm{i}}{Y})],\end{eqnarray} \tag{ 32 }$$

where $X={N}^{-1}{\sum }_{k}{x}_{k}$ and $Y={N}^{-1}{\sum }_{k}{\dot{x}}_{k}$ are two mean fields, while $\varepsilon$ is the coupling strength. Finally, the control parameter γ accounts for a phase shift of the coupling term: it determines whether the interaction is attractive or repulsive.

It is well-known that uncoupled units in equation (32) exhibit limit-cycle oscillations, while the nonlinearity parameter ζ determines the stability of the limit cycle. Below we consider $\zeta =5;$ for this parameter the transversal Lyapunov exponent is −7.358. Therefore, adiabatic elimination of the amplitude is rather meaningful.

An appropriate order parameter is

$$\begin{eqnarray}&&\rho ={\rm{rms}}(X)/{\rm{rms}}(x),\end{eqnarray} \tag{ 33 }$$

where 'rms' means root mean square of the time evolution. The splay state is characterized by a constant mean field and thereby $\rho =0$. In the fully synchronous state, for identical oscillators, the microscopic and macroscopic dynamics are equivalent to one another so that $\rho =1$.

The scenario resulting for N = 1000, $\varepsilon =0.05$, $\zeta =5$, and transient time $7.5\times {10}^{5}$ is reported in figure 9. It is reminiscent of that observed in the biharmonic model, with some differences. SCPS establishes itself in between the parameter range where full synchrony is stable (below $\gamma \approx -0.6$) and the region where ρ is negligible (above $\gamma \approx 0.05$). In this case, the mean field is slower than the individual oscillators, as in the upper eyelet of the biharmonic model, see figure 2. The regime characterized by a vanishing ρ is not asynchronous but a symmetric nine-cluster state (see also below). In fact, the splay state turns out to be unstable in the full range of γ values that we have explored. Finally, in the interval $[-0.6,-0.2]$ we see a slow convergence towards a two-cluster state (the large value of the order parameter ρ is due to the closeness between the two clusters).

Zoom In Zoom Out Reset image size

The dynamics of the coupled system equation (32) can be better understood by performing a phase reduction. This can be done by introducing a phase variable for each individual oscillator, making reference to the uncoupled limit (i.e. $\varepsilon =0$): ${\phi }_{k}=2\pi {t}_{k}/T$, where t_k is the time elapsed from the passage through a chosen origin x = 0, $\dot{x}\gt 0$, while T is the oscillation period.

In the weak coupling limit, the original Rayleigh oscillators (32) can be mapped onto a Winfree model

$$\begin{eqnarray}&&{\dot{\phi }}_{k}={\omega }_{0}+\displaystyle \frac{\varepsilon }{N}{\rm{\Gamma }}({\phi }_{k})\displaystyle \sum _{j}Z({\phi }_{j}),\end{eqnarray} \tag{ 34 }$$

where ${\rm{\Gamma }}(\phi )$ is the phase response curve (PRC), while $Z(\phi )$ is the forcing function. ${\rm{\Gamma }}(\phi )$ can be obtained by following a standard approach: it corresponds to the phase shift imposed by an infinitesimal kick when the phase of the oscillator is ϕ. The resulting PRC is reported in figure 10 (see the red dotted curve). The forcing function Z is instead obtained by expressing the coupling term due to the jth oscillator $\mathrm{Re}[{{\rm{e}}}^{{\rm{i}}\gamma }({x}_{j}+{\rm{i}}{\dot{x}}_{j})]$ as a function of ${\phi }_{j}$ (see the blue dotted curve in figure 10). Finally, following [12, 38], one can further map the Winfree model (34) onto a Kuramoto–Daido model, by computing the convolution of Γ and Z, i.e.

$$\begin{eqnarray}&&{G}_{{\rm{R}}}(\phi )=\int {\rm{\Gamma }}(\phi -\psi )Z(\psi ){\rm{d}}\psi .\end{eqnarray} \tag{ 35 }$$

The resulting coupling function corresponds to the black curve in figure 10. Fourier analysis shows that even modes are absent: this follows from the symmetry of the limit cycle: adding π to the phase results in changing the sign of both the PRC and forcing term.

Zoom In Zoom Out Reset image size

In order to test the validity of the Kuramoto–Daido reduction, it has been simulated for eight harmonics. The resulting scenario is in close agreement with that one exhibited by the original ensemble of Rayleigh oscillators, including the observation of nine-cluster states. Within the Kuramoto–Daido representation, one can easily perform a stability analysis of both the splay and synchronous state. As discussed in [38], the stability of the splay state is determined by the imaginary components of the Fourier modes of G_R. In table 1 we show the contribution of the first eight non-vanishing components for $\gamma =0.2$, where a zero order parameter is observed. In fact, the splay state turns out to be unstable because of the contribution of the 7th and higher harmonics. This is at variance with the biharmonic model, where the stability was determined by the first Fourier mode.

Table 1.  Eigenvalues associated to the stability of the splay state for γ = 0.2. The index refers to the Fourier mode which they are associated with.

Index	Real part	Imaginary part
1	−1.6 × 10−2	2.1 × 10−1
3	−5.6 × 10−2	1.5 × 10−2
5	−1.6 × 10−2	1.9 × 10−2
7	3.0 × 10−3	1.0 × 10−2
9	4.0 × 10−3	1.7 × 10−3
11	1.5 × 10−3	5.0 × 10−4
13	3.8 × 10−4	4.3 × 10−4
15	7.9 × 10−5	1.8 × 10−4

As for the fully synchronous state, its stability is determined by the sign of ${G}_{{\rm{R}}}^{\prime }(0)$. The bifurcation point where it changes stability is $\gamma \simeq -0.573$ in agreement with the numerical simulations of the original model (32).

For larger values of γ and below −0.2, the Kuramoto–Daido model possesses two-cluster states. Following the analysis outlined in section 4, we conclude that two-cluster states are effectively stable. In fact, for $\gamma =-0.4$, the inter-cluster exponent ${\lambda }_{{\rm{I}}}=-{G}_{{\rm{A}}}^{\prime }({\delta }_{0})=-0.0669$ is negative, while ${\lambda }_{E}^{+}=0.0869$, and ${\lambda }_{E}^{-}=-0.1146$. Accordingly the average exponent ${\lambda }_{{\rm{a}}}=-0.01385$ is negative, in spite of one of the two intra-cluster exponents being positive. Altogether, the scenario is similar to that observed in the biharmonic model.

In this paper we have shown that SCPS is a general phenomenon, arising in many setups of globally coupled oscillators. This regime naturally emerges if the system is close to the border between full synchrony and asynchrony. In fact, the minimal requirement for SCPS to arise is the presence of two harmonics in the coupling function $G(\phi )$. This condition is naturally fulfilled in weakly coupled oscillators away from the Hopf bifurcation. As an example we have indeed verified that SCPS arises also in an ensemble of Rayleigh oscillators, where an approximate coupling function containing eight Fourier modes suffices to quantitatively reproduce the dynamics of the original model. Altogether, SCPS is yet another dynamical regime that cannot be produced by the standard Kuramoto–Sakaguchi model. In fact, both SCPS with and without frequency difference can be obtained in a model, based on strictly sinusoidal coupling function, but it requires a dependence on the coupling strength and/or the phase shift of the sine-function on the order parameter, i.e. SCPS can be observed in case of nonlinear coupling only. This limitation disappears for our minimal Kuramoto–Daido model.

The mathematical structure of Kuramoto–Daido models allows for a semi-analytic treatment: it is not only possible to determine the probability distribution of the phases, but also to perform a linear stability analysis and thereby determine the parameter range, where SCPS can be effectively observed. In the biharmonic model and the Rayleigh oscillators, the loss of stability of SCPS drives the system towards a HC, which can itself be interpreted as a (more structured) form of SCPS: here, besides a difference between the microscopic and mean field frequencies, a pulsation of the amplitude is present. In LIF neurons, instead, SCPS is always stable, while two-cluster states are always unstable.

It would be interesting to discover whether and under which conditions other kinds of bifurcation can drive SCPS towards more complex forms of collective dynamics. This question is related to that of identifying the number of relevant collective variables. A fairly trivial answer can be given in the case of perfect clusters, as it boils down to studying low-dimensional networks composed of a few 'supernodes' (the clusters themselves). The question is much less trivial in the context of smooth distributions such as those associated to SCPS. Possibly, the first example of a complex behavior in a globally coupled partially synchronized system was given in [45], where the Morris–Lecar neuronal oscillators were analyzed numerically. However, this is a setup where phase-reduction is not globally possible. More recently, evidence of a chaotic collective behavior has been found in a population of quadratic IF neurons [46], where, however, the higher dynamical complexity is triggered by the presence of delayed interactions. So the question whether identical phase-oscillators can lead to collective chaos is still open.

Another open general problem is that of using the information encoded in the coupling function to predict whether SCPS and/or cluster states can be generated. In the case of a sinusoidal coupling, the situation is simplified by the fact that clusters are not possible: their existence is excluded by the Watanabe–Strogatz theory [43, 44], see also [17, 34]. Thus, when the splay and synchronous states are both unstable, SCPS is the only possible solution. In more general contexts cluster states sometimes coexist with SCPS, as well as with chimera-like solutions.

We wish to acknowledge A Pikovsky and M Zaks for useful discussions. This work has been financially supported by the EU project COSMOS (642563).