A minimal model of self-consistent partial synchrony

We show that self-consistent partial synchrony in globally coupled oscillatory ensembles is a general phenomenon. We analyze in detail appearance and stability properties of this state in possibly the simplest setup of a biharmonic Kuramoto–Daido phase model as well as demonstrate the effect in limit-cycle relaxational Rayleigh oscillators. Such a regime extends the notion of splay state from a uniform distribution of phases to an oscillating one. Suitable collective observables such as the Kuramoto order parameter allow detecting the presence of an inhomogeneous distribution. The characteristic and most peculiar property of self-consistent partial synchrony is the difference between the frequency of single units and that of the macroscopic field.

. As a result, the caption of figure 1 is to be modified, stating that the splay state is stable only in the central (cyan) square region ( 2  3 2 1,2   p g p ). The entire analysis performed in the paper remains unaffected, since it refers to the eyelet e 1 .
(C) The equation at the end of page 9 has to be changed to 2 n pn = W -.
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 socalled 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][23][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 Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 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][29][30][31] with a suitable coupling function f D G ( ), where f D 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.

The model
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 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 A standard way to classify the configurations of an ensemble of oscillators is via the set of complex order parameters where R 1 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 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)  This pointwise region can be made structurally stable by assuming that g 1 is a function of the order parameter R 1 and of the coupling strength e [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 The marginal stability line is composed of two sinusoidal curves G + , shown in the stability diagram in figure 1. Altogether, the synchronous state is stable in the two yellow and green regions of the diagram. (cyan and green areas); the fully synchronous solution is stable above the upper and below the lower solid curves (yellow and green areas); two-cluster (anti-phase) states are stable inside the area delimited by the purple curve (i.e. below g p = 2 2 and above g p = 3 2 2 ). G + and Gidentify two pairs of sinusoidal curves where the synchronous state and the antiphase two-cluster states solution lose stability, respectively. e 1 and e 2 denote two pairs of eyelets (bounded by G + and G -) where non-trivial dynamics between synchrony and asynchrony can be expected. The box identifies the parameter region numerically investigated in this paper, see figure 2.
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 d 1 is that one characterizing the stability of the first Fourier mode of the perturbation. In the present setup d p g g = + cos i sin , 6 so that d 1 depends only on g 1 . As a result, the splay state is stable in a rectangular region where g < cos 0 1 (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 f f  -, the same dynamics is found upon mapping g p g  -2 1 1 and g p g  -2 2 2 . The additional symmetry is, therefore, only accidental.

Self-consistent partial synchronization
In the thermodynamic limit one can investigate the various regimes by studying the evolution of the probability density f P t , ( )of oscillators with phase f at time t. It satisfies the continuity equation The splay state corresponds to a flat and constant density p = P 1 2 ( ). 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 q f = -Wt and q f For a suitably chosen Ω there exists a stationary time-independent solution q Q 0 ( ), which satisfies the equation where η is the probability flux. Notice that ph 2 can be interpreted as the average microscopic frequency of the single oscillators in the moving frame. y Q 0 ( ) can be expanded in Fourier modes, (the harmonics coincide with the generalized order parameters defined in equation (3) for a finite ensemble of oscillators) The stationary solution can be obtained by solving equation Since the phase of the solution is arbitrary, we are free to fix it by imposing that Z 1 is real. By considering that η can be determined by imposing a normalization on Q 0 , the above equation contains four unknowns: W, Z 1 , and Z 2 (the last variable is complex). They can be determined self-consistently by imposing The solution can be found by searching for a fixed point in a four-dimensional space.

By expanding the perturbation in Fourier series
one can rewrite the integral in the previous equation as ò y y q y This equation can be expressed as a Fourier series. With the help of equation (11), it is found 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 Q 0 . In other words, the corresponding matrix is sparse: it is pentadiagonal with two full rows. As each derivative is multiplied by m, = q 0 0 . This is a straightforward consequence of the conservation of the total probability; therefore q 0 can be eliminated as it does not contribute to the eigenvalues. By further looking at the evolution equation for q 2 , we see that it involvesq 1 , so that the negative modes must be included as well. Since

Clusters and HCs
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 equalsize clusters. The phases of the two families follow the differential equations a y a y a y = The separation d a y =satisfies the equation where G A is the anti-symmetric component of G. Two-cluster solutions are identified by the zeros of G A ; d = 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 d d g g d = + G a sin cos 2 cos cos , 18 There are various solutions of the equation corresponds to an antiphase twocluster 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 d 0 , the same state can be seen as characterized by a separation p d -2 0 , if the two clusters are exchanged. These states exist only in the parameter region delimited by the curves where g g  = a cos 2 cos 0. , i.e. the two-cluster solution bifurcates from the fully synchronous one.) The minus sign, instead, corresponds to the curve Gwhere the twocluster state becomes the antiphase one. As a result, non-trivial clustered solutions with d p ¹ exist only in the regions delimited by G + and G -(see the two pairs of eyelets e 1 and e 2 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 l I measures the stability against perturbation of the phase-separation between the two clusters. From equation (17) it follows In the biharmonic model, The intra-cluster exponents l  E measure the stability of the width of each cluster. From the equations of motion it is readily found that The solution with the plus (minus) sign refers to the cluster that is lagging behind (leading) by d 0 . The maximal eigenvalue is therefore The inter-cluster exponent of the antiphase solution, , is negative in between the upper and lower branch of Gand vanishes along G -, 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 in the region where other clusters do not exist. The stability of the antiphase state to intracluster perturbations is controlled by (see equation (23)) so that this state is stable within the region delimited by the purple curve in figure 1.
The non-trivial two-cluster state (d p ¹ 0 ) is unstable against inter-cluster perturbations within the eyelets e 2 (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 1 . Its intra-cluster stability can be determined from equation (23). By taking into account that cos 1 cos sin sin 4 sin cos . 25 Therefore, such a two-cluster state is unstable for p g p < < 2 3 2 1 . A detailed analysis for a = 0.2 reveals that the solution is unstable also within the eyelet e 1 . This is not yet the end of the story. The opposite sign of the two exponents l  E 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 In the biharmonic model , 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.

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 was2 10 4 time units.
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 10 , 0.5 10 12 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 −8 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 1 (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).
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 g » 4.   Evolution of the phases of all oscillators (black) and of the mean field (red)-after subtracting the average growth. For rather long epochs oscillators are grouped into two clusters with almost identical phases. However, these states are unstable and alternately each group widens until its width becomes of order one, and then shrinks again. Notice, that there is no transfer of elements between the two groups, but the mean field frequency is larger than that of oscillators. Indeed, we see that in the coordinate frame co-rotating with the frequency w, the phase b 1 drifts away. Parameters are g p = 2 , g = 1.35 1 , N=1000 (the dynamics is preserved for the ensemble size as large as N = 5000). The initial conditions are a perturbed splay state. Size of the two groups is close but not equal to N 2. 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 g p = 2 , g = 1.5 1 , the SCPS dynamics seems to be the asymptotic state, at least it survives for 10 7 time units.
Further states have been occasionally detected (see the cyan symbol in figure 2). In particular three-cluster states are found for g = 4.75 2 and   g 1.58 1.66 1 , 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 g p = Upon decreasing g 1 from g s , SCPS is first born through a Hopf bifurcation from the splay state; R 1 and R 2 become strictly larger than zero and a finite mean-field frequency W, 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 W.
Interestingly, W 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 f f 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 ( w W > and w W < ) 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 g 1 , we enter a region where the dynamics converges to HCs which are characterized by a sudden increase of the (average) R 1 , while no discontinuity is exhibited by g R 2 1 ( ). When, finally, g f is approached, unsurprisingly R 1 and R 2 converge to 1, while both W 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. suggests the possible existence of multistability. Accordingly, we have performed two additional series of simulations, by progressively increasing (decreasing) g 1 , and choosing the new initial condition as a slightly perturbed version of the final configuration for the previous value of g 1 . As shown in figure 5, a bistability region is indeed found where HCs and SCPS are simultaneously stable.

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 g g , f p [ ] and it reconnects to the fully synchronous state. As for the microscopic frequency ω, given by it also agrees with the numerical simulations.
The stability analysis carried out in section 3 allows determining g 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  and g = 1.45 1 . Circles and crosses correspond to a truncation after 100 and 500 modes, respectively. 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 m 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.
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 p -D 2 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 t 1 ln . 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 g 1 . The data confirms that linear instability occurs below g p : the bifurcation is of Hopf type. SCPS is maximally unstable around g » 1.33 1 . By further decreasing g 1 , both the real and the imaginary parts decrease to zero, while approaching g 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 g 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 g 2 .

Other setups
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 pulsecoupled 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.

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 [38] for a precise definition of the various parameters-notice that here we are using a different notation -j instead of f-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 a = 6), the average exponent l > 0 a , 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 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 a n n is positive and increasingly large. As a result, on the basis of the first two polynomial terms of d G A ( ), one finds that it vanishes also for d n a = 2 . . This makes the sum of the derivatives in 0 and d 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 j. One should however remember that the original function G is continuous but not even of class  ; 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 j ej j = -+ G A . 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.

Rayleigh oscillators
In order to provide further evidence on the ubiquity of SCPS, we present a numerical study of globally coupled identical Rayleigh oscillators 4 . The equations are z w e -- 1˙a re two mean fields, while e 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 z = 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 r = X x rms rms , 33 ( ) ( ) ( ) where 'rms' means root mean square of the time evolution. The splay state is characterized by a constant mean field and thereby r = 0. In the fully synchronous state, for identical oscillators, the microscopic and macroscopic dynamics are equivalent to one another so that r = 1.
The scenario resulting for N=1000, e = 0.05, z = 5, and transient time7.5 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 g » -0.6) and the region where ρ is negligible (above g » 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 twocluster state (the large value of the order parameter ρ is due to the closeness between the two clusters).
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. e = 0): , where t k is the time elapsed from the passage through a chosen origin x=0, > x 0 , while T is the oscillation period. In the weak coupling limit, the original Rayleigh oscillators (32) can be mapped onto a Winfree model where f G( ) is the phase response curve (PRC), while f Z ( ) is the forcing function. f G( ) 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 f. 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 as a function of f 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.
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. 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 g = 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.
As for the fully synchronous state, its stability is determined by the sign of ¢ G 0 R ( ). The bifurcation point where it changes stability is g - 0.573 in agreement with the numerical simulations of the original model (32).

Conclusions
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 f G ( ). 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 chimeralike solutions.