A minimal model of 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 to distribution of phases to an oscillating one. Suitable collective observables such as the Kuramoto order parameter allow detecting the presence of a inhomogeneous distribution. The characteristic and most peculiar property of partial synchrony is the difference between the frequencies of single units and that of the macroscopic field.

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 phenomena are observed, when a control parameter is continuously varied. Thus, the microscopic dynamics is quasiperiodic. Examples of such dynamics are: integrateand-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 [27,26].
In the weak-coupling limit oscillators are effectively described by a Kuramoto-Daido phase model [28,29,30,31] with a suitable coupling function G(∆φ), where ∆φ is the phase difference between any two interacting oscillators. In the standard, widely used, Kuramoto-Sakaguchi model [32,9,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 [34,17]. A much richer dynamics, including formation of clusters [19] and of heteroclinic cycles [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, 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 [36,24]. The additional stability analysis confirms the numerical evidence of the onsed 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 heteroclinic cycle (HC) [35,37]. We find that more hamornics are needed to ensure that two-cluster states and the corresponding HCs are both unstable. Moreover, we find that heteroclinic cycles can be viewed at as a kind of quasiperiodic partial synchrony.
The paper is organized as follows. In Sec. II 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 Sec. III, we analyze the occurrence of SCPS, show how it can be treated and finally develop the formalism needed to perform the stability analysis. Sec. IV 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 twocluster 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 Sec. V. 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 Sec. VI, 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 Ref. [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 readṡ where the coupling constant has been eliminated by performing a suitable rescaling of time.
In most of the paper we consider the 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 [32,9], 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, Eq. (1) can be rewritten aṡ 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 Eqs. (2,1) reduce to the famous Kuramoto-Sakaguchi model. In this case it is known that the fully synchronous solution φ k = φ is stable if and only if |γ 1 | < π/2, while the stability condition of the splay state is exactly opposite: asynchrony is stable for π/2 < |γ 1 | < 3π/2 and unstable otherwise. A partially synchronous 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 γ 1 is a function of the order parameter R 1 and of the coupling strength ε [36,24] (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 Γ + , shown in the stability diagram in Fig. 1. Altogether, the synchronous state is stable in the two yellow and green regions of the diagram. On the other hand, the stability of the asynchronous state, which we analyze in the thermodynamic limit N → ∞, is still determined only by the amplitude of the first mode of the perturbation (see [38]) so that it depends only on γ 1 . As a result, the splay state is stable in a rectangular region where cos γ 1 < 0 (see the cyan and green areas in Fig. 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 Fig. 1 has a reflectional symmetry with respect to both the horizontal and vertical semi axes. Physically, there is only one symmetry: if φ → −φ, the same dynamics is found upon mapping γ 1 → 2π − γ 1 and γ 2 → 2π − γ 2 . The additional symmetry is, therefore, only accidental. The splay state is stable in the rectangular region π/2 < γ 1 < 3π/2 (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 γ 2 = π/2 and above γ 2 = 3π/2). Γ + and Γ − identify 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 Γ + and Γ − ) where nontrivial dynamics between synchrony and asynchrony can be expected. The box identifies the parameter region numerically investigated in this paper, see Fig. 2.

Partial synchronization
In the thermodynamic limit one can investigate the various regimes by studying the evolution of the probability density P (φ, t) of oscillators with phase φ at time t. It satisfies the continuity equation The splay state corresponds to a flat and constant density P = 1/(2π). As already illustrated in Ref. [38], SCPS typically manifests itself as a rotating nonuniform density, which emerges past a Hopf bifurcation. It is therefore convenient to perform a change of variables, introducing the angle θ = φ − Ωt and Q(θ, t) = P (φ, t). The corresponding evolution equation takes the form For a suitably chosen Ω there exists a stationary time-independent solution Q 0 (θ), which satisfies the equation where η is the probability flux. Notice that 2πη can be interpreted as the average microscopic frequency of the single oscillators in the moving frame. Upon expanding Q 0 (ψ) in Fourier modes, (the harmonics coincide with the generalized order parameters defined in Eq. (3) for a finitesize ensemble of oscillators), we can solve Eq. (9), obtaining 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: Ω, Z 1 , and Z 2 (the last variable is complex). They can be determined self-consistently by imposing for k = 1, 2. The solution can be found by searching for a fixed point in a four-dimensional space.

Stability analysis
Consider an infinitesimal perturbation q(θ, t) of Q 0 (θ) and linearise Eq. (8), making use of Eq. (9). One obtains By expanding the perturbation in Fourier series, one can rewrite the integral in the previous equation as This equation can be expressed as a Fourier series. With the help of Eq. (11), it is found The m-th Fourier mode of the perturbation is coupled with the four nearest neighbour 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; thereforeq 0 can be eliminated as it does not contribute to the eigenvalues. By further looking at the evolution equation forq 2 , we see that it involvesq −1 , so that the negative modes must be included as well. Sinceq −m =q * m , it is convenient to separateq m into real and imaginary part (q m = u m + iv 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 µ R + iµ I of the resulting matrix can be then computed by considering a sufficiently large number of Fourier modes.

Clusters and heteroclinic cycles
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, The separation δ = α − ψ satisfies the equatioṅ where G A is the anti-symmetric component of G. Two-cluster solutions are identified by the zeros of G A ; δ = 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 There are various solutions of the equation G A (δ) = 0 besides δ 0 = 0: δ 0 = π corresponds to an antiphase two-cluster state. Two further solutions can be found by setting to zero the expression in parentheses in Eq. (18); however, these solutions represent only one physically meaningful state. Indeed, given a two-cluster state characterized by a separation δ 0 , the same state can be seen as characterized by a separation 2π − δ 0 , if the two clusters are exchanged. These states exist only in the parameter region delimited by the curves where This equation with a plus sign defines the curve Γ + which coincides with the bifurcation line where the synchronous state loses stability, see Fig. 1. (In fact, cos δ 0 = 1 means δ 0 = 0, i.e. the two-cluster solution bifurcates from the fully synchronous one.) The minus sign, instead, corresponds to the curve Γ − where the two-cluster state becomes the antiphase one. As a result, nontrivial clustered solutions with δ = π exist only in the regions delimited by Γ + and Γ − (see the two pairs of eyelets e 1 and e 2 in Fig. 1). The stability of a two-cluster state is determined by the value of the inter-and intracluster exponents. The inter-cluster exponent λ I measures the stability against perturbation of the phase-separation between the two clusters. From Eq. (17) it follows In the biharmonic model, G A (δ 0 ) = cos γ 1 cos δ 0 + 2a cos γ 2 cos 2δ 0 . The intra-cluster exponents λ ± 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 δ 0 . The maximal eigenvalue is therefore, The inter-cluster exponent of the antiphase solution, G A (π) = − cos γ 1 + 2a cos γ 2 , is negative in between the upper and lower branch of Γ − and vanishes along Γ − , 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 (π)G (0) < 0 in the region where other clusters do not exist. The stability of the antiphase state to intra-cluster perturbations is controlled by (see Eq. (23)), so that this state is stable within the region delimited by the purple curve in Fig. 1.
The nontrivial two-cluster state (δ 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 Eq. (23). By taking into account that G A (δ 0 ) = 0, we find that Therefore, such a two-cluster state is unstable for π/2 < γ 1 < 3π/2. 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 λ ± 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 In the region of interest, λ a < 0 if cos γ 1 > 0, i.e the fluctuations of the cluster widths on average decrease, without ever collapsing onto it. This is nothing but an attracting heteroclinic cycle (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 Fig. 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 Eqs. (1,2) (starting from an initial condition close to the splay state) is summarized in Fig. 2: the symbols identify the different asymptotic states, while the curves correspond to the marginal stability lines (determined theoretically, see Fig. 1). The simulation of Eqs. (1,2) was performed for N = 1000; larger ensemblesize have been considered for several points without observing any essential difference. The transient time was 2 · 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 Eqs. (1,2)) have been taken as uniformly distributed in the interval [−0.5 · 10 −12 , 0.5 · 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 heteroclinic cycle. 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 Fig. 3b, 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 (Fig. 3a). 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 (Fig. 3b); 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 nonstationary distribution. A movie showing the time evolution of heteroclinic cycles 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 γ 2 ≈ 4.75, γ 1 ≈ π/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 -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 ω, the phase β 1 drifts away. Parameters are γ 2 = π, γ 1 = 1.35, 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. A movie showing the time evolution of this regime is included in the supplementary data (see supplementary movie 1). the transient time, some points that appear as SCPS states in Fig. 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 γ 2 = π, γ 1 = 1.5, the partially synchronous 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 Fig. 2). In particular three-cluster states are found for γ 2 = 4.75 and 1.58 ≤ γ 1 ≤ 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 γ 2 = π in the parameter space. The results are shown in Fig. 4, where one can recognize three critical points: (i) γ s = π/2 signals the loss of stability of the splay state; (ii) γ f = arccos(2a) ≈ 1.159 signals the loss of stability of the fully synchronous state; (iii) γ p ≈ 1.401 signals the loss of stability of SCPS.
Upon decreasing γ 1 from γ 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 Ω, 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 Ω.
Interestingly, Ω 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 φ → −φ 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 [36,24] both cases (Ω > ω and Ω < ω) are possible. Altogether, the existence of a finite difference between microsocpic and macroscopic frequencies is a key signature of a quasiperiodic SCPS.  Upon further decreasing γ 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 R 2 (γ 1 ). When, finally, γ f is approached, unsurprisingly R 1 and R 2 converge to 1, while both Ω 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 γ 1 = γ p suggests the possible existence of multistability. Accordingly, we have performed two additional series of simulations, by progressively increasing (decreasing) γ 1 , and choosing the new initial condition as a slightly perturbed version of the final configuration for the previous value of γ 1 . As shown in Fig. 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 Eq. (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 Fig. 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 [γ 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 Sec. III allows determining γ p by solving the eigenvalue problem (15). A typical spectrum is plotted in Fig. 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 µ 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 2π − ∆ 0 offers the chance to appreciate the role of the weakly attracting directions. As seen in Fig. 7, the gap size goes to zero, but it does so in an extremely slow way, namely as a/ ln t.
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 Eq. (C2) in [38]), so that, when the number  Figure 6. Eigenvalues of SCPS for γ 2 = π and γ 1 = 1.45. Circles and crosses correspond to a truncation after 100 and 500 modes, respectively. 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 Fig. 8 we plot the real and imaginary part of the most unstable eigenvalue versus γ 1 . The data confirms that linear instability occurs below γ p : the bifurcation is of Hopf type. SCPS is maximally unstable around γ 1 ≈ 1.33. By further decreasing γ 1 , both the real and the imaginary parts decrease to zero, while approaching γ 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 γ f are not reliable. It is nevertheless instructive to notice that the weak instability is consistent with the observation of partial synchrony over extremely long time scales mentioned in the beginning of this section. The correct identification of the critical point is further confirmed in Fig. 2, where the triangles are the outcome of the stability analysis for three different choices of γ 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 pulse-coupled LIF neurons [25]. In Ref. [38], it has been shown that the system can be effectively described by a Kuramoto-Daido model with a coupling function containing several nonnegligible 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 twodimensional Rayleigh oscillators esch 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 (see Ref. [38] for a precise definition of the various parameters -notice that here we are using a different notation -ϕ 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 α = 6), the average exponent λ a > 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 which is negative and finite. Comparison with Eq. (17), implies that the synchronous solution is always unstable. As for the second derivative, the leading order G A (0) = α 3 e τ − 1 2ν(a + gν) is positive and increasingly large. As a result, on the basis of the first two polynomial terms of G A (δ), one finds that it vanishes also for This phase shift identifies a two-cluster state; its value decreases as the cubic power of the pulse-width (equal to 1/α). Under this quadratic approximation, the slope of G A in 0 and δ 0 are equal and opposite to one another, so that λ 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 α 4 . As a result its contribution to the derivative in δ 0 is of order α 4 δ 2 0 ≈ 1/α 2 . This makes the sum of the derivatives in 0 and δ 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 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 G A (ϕ) = −εϕ + Aϕ 3 .
If A > 0 and ε > 0, then there exists a second zero δ 0 = ε/A, which corresponds to the clustered solution. Therefore so that smooth interaction functions necessarily lead to effectively stable two-cluster states (at least with respect to intracluster perturbations). Even more, Eq. (27) yields that in the biharmonic model no effectively unstable cluster may exist: only heteroclinic cycles. 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 ‡. The equations arë where X = N −1 k x k and Y = N −1 kẋ k are two mean fields, while ε 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 Eqs. (32) exhibit limit-cycle oscillations, while the nonlinearity parameter ζ determines the stability of the limit cycle. Below we consider ζ = 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 where "rms" means root-mean-square of the time evolution. The splay state is characterized by a constant mean field and thereby ρ = 0. In the fully synchronous state, for identical oscillators, the microscopic and macroscopic dynamics are equivalent to one another so that ρ = 1. The scenario resulting for N = 1000, ε = 0.05, ζ = 5, and transient time 7.5 · 10 5 is reported in Fig. 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 γ ≈ −0.6) and the region where ρ is negligible (above γ ≈ 0.05). In this case, the mean field is slower than the individual oscillators, as in the upper eyelet of the biharmonic model, cf. Fig. 2. The regime characterized by a vanishing ρ is not asynchronous but a symmetric 9-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).
The dynamics of the coupled system Eq. (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. ε = 0): φ k = 2πt k /T , where t k is the time elapsed from the passage through a chosen origin x = 0,ẋ > 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 Γ(φ) is the phase response curve (PRC), while Z(φ) is the forcing function. Γ(φ) 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 Fig. 10 (see the red dotted curve). The forcing function Z is instead obtained by expressing the coupling term due to the jth oscillator Re e iγ (x j + iẋ j ) as a function of φ j (see the blue dotted curve in Fig. 10). Finally, following Refs. [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 Fig. 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 9-cluster states. (c) Figure 9. Self consistent partial synchrony is observed in model (32)  Within the Kuramoto-Daido representation, one can easily perform a stability analysis of both the splay and synchronous state. As discussed in Ref. [38], the stability of the splay state is determined by the imaginary components of the Fourier modes of G R . In Tab. 1 we show the contribution of the first eight nonvanishing components for γ = 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 R (0). The bifurcation point where it changes stability is γ −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 Sec. 4, we conclude that two-cluster states are effectively stable. In fact, for γ = −0.4, the inter-cluster exponent λ I = −G A (δ 0 ) = −0.0669 is negative, while λ + E = 0.0869, and λ − E = −0.1146. Accordingly the average exponent λ 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.

Conclusions
In this paper we have shown that self-consistent partial synchrony (SCPS) is a general phenomenon, arising in many setups of globally coupled oscillators. This regime naturally  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. 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(φ). 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 partial synchrony 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 heteroclynic cycle, 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 integrate-and-fire neurons [46], where, however, the higher dynamical complexity is triggered by the presence of delayed interactions. So the question whether identical phaseoscillators 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 [34,17]. 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.