Synchronization clusters emerge as the result of a global coupling among classical phase oscillators

When large ensembles of phase oscillators interact globally, and when bimodal frequency distributions are chosen for the natural frequencies of the oscillators themselves, Bellerophon states are generically observed at intermediate values of the coupling strength. These are multi-clustered states emerging in symmetric pairs. Oscillators belonging to a given cluster are not locked in their instantaneous phases or frequencies, rather they display the same long-time average frequency (a sort of effective global frequency). Moreover, Bellerophon states feature quantized traits, in that such average frequencies are all odd multiples (±(2n − 1), n = 1, 2...) of a fundamental frequency Ω1. We identify and investigate (analytically and numerically) several typical bifurcation paths to synchronization, including first-order and second-order-like. Linear stability analysis allows to successfully solve the critical transition point for synchronization. Our results highlight that the spontaneous setting of higher order forms of coherence could be achieved in classical Kuramoto model.


Introduction
Synchronization phenomena are ubiquitous in physics, chemistry, biology, engineering, and human society. In particular, synchronization in networked oscillators has attracted great attention in the last few decades, due to its potential in applications [1]. Recently, a special phase coherence (called the Bellerophon state) has been characterized in globally coupled nonidentical phase oscillators. Such a state is a quantized, time-dependent, clustered state which emerges close to a first-order-like transition to synchronization [2][3][4][5]. There are two main backgrounds of this work: the first is explosive synchronization, i.e. an abrupt, first-order-like, transition to a coherent (synchronous) state of globally coupled phase oscillators, and the second is the Bellerophon states which were so far found only in generalized Kuramoto models [3][4][5].
As for explosive synchronization, after the seminal works with Kuramoto model on scale-free networks [6,7], it was argued that a sufficient condition be the setting of proper correlations between the natural frequencies of the phase oscillators and the networkʼs node degrees. Later, Zhang et al proposed a frequencyweighted Kuramoto model, which can exhibit first-order-like synchronization transition on generic networkʼs topologies for typical frequency distributions [8][9][10][11][12]. Furthermore, it was revealed that the mechanism at the basis of such abrupt synchronization transition is similar to that underlying explosive percolation [2], where the formation of a giant component is controlled by a suppressive rule [10]. Finally, first-order-like synchronization transition was also found in adaptive and multi-layer networks [11]. As for the Bellerophon states, they have been observed so far in generalized Kuramoto models: either the frequency-weighted Kuramoto model [3,4] or the Kuramoto model with conformists and contrarians [5].
In our paper, we focus on the classical Kuramoto model, and we consider a bimodal distribution for the selection of the natural frequencies of the phase oscillators. In these conditions, we first point out that the regime Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. called standing wave by Crawford [13] and identified in [14] is actually not a standing wave but is, in fact, a Bellerophon state. Then, we carefully investigate five regimes on the parameter plane, and identify several typical bifurcation paths to synchronization, including both first-order and second-order-like. Moreover, the robustness of our findings is demonstrated by the fact that a similar scenario emerges for different bimodal frequency distributions. Finally, we apply linear stability analysis (to the entire system, and not to a reduced lower-dimensional one, as several other studies have done in the past by using the Ott-Antonsen method), and successfully solve the critical transition point for synchronization. All our theoretical predictions are well consistent with numerical simulations.

Model and theoretical analysis
Let us start by introducing the classical Kuramoto model [15], which describes the evolution of an ensemble of N globally coupled phase oscillators: where dot denotes a temporal derivative, and i=1, K, N. The instantaneous phases θ i (t) evolve as q w = ( ) t i i when the coupling strength κ is set to zero. The oscillators' natural frequencies ω i are taken from a distribution g (ω) which is usually considered symmetric and centered at zero, i.e. g(ω) = g(−ω) in the thermodynamic limit (  ¥ N ). To monitor the setting of coherence in the ensemble as κ increases, one can rely on the (complex valued) order parameter defined as where 0r1 is the modulus of the mean-field, and Ψ is the average phase. r=0 (r=1) corresponds to a fully incoherent (fully coherent) state, while intermediate values of r characterize partially coherent states where a portion of the ensemble is organized into one, or many synchronized clusters, and coexists with a sea of not entrained oscillators. For the case of a unimodal distribution g(ω), a continuous transition from incoherence to full coherence is found at the critical point k = p ( ) c g 2 0 . For bimodal distributions the results are somehow inconclusive, and point to the formation, at κ>κ c , of two symmetric counter-rotating clusters of synchronized oscillators, a state which was (inaccurately, as we will shortly see) named as standing wave.
For a direct comparison with existing literature, we start by considering the same family of symmetric bimodal distributions used in [14], i.e. the sum of two Lorentzian distributions is the Lorentzian distribution having 2Δ as the half width at half maximum of each peak, and ±ω 0 as center frequencies. Note that equation (3) includes also a unimodal distribution as the trivial case for w D > 3 0 , i.e. when the two peaks are close enough to each other.
Several detailed studies of equation (1) with a bimodal distribution were carried out [13,14,16,17], and the conclusion is that the model exhibits the route to synchronization depicted in figure 1. As reported in [14], three types of attractor are found: the incoherent (green area in figure 1(a)) and partially synchronized states (yellow area) corresponding to the trivial and nontrivial fixed points, and what was initially thought to be a standing wave (red, dark and light brown areas) corresponding to a limit-cycle solution. The transitions between these states are mediated by transcritical (TC), saddle-node (SN), Hopf (HB), and homoclinic (HC) bifurcations. However, the bifurcation and stability analysis were made using a reduction method developed by Ott and Antonsen [18] and, as the Authors themselve alert, some of the actual systemʼs behavior might be lost due to the fact that the reduced system only represents a special restricted class of all possible solutions of the original system. Therefore, it is desirable to conduct the analysis on the full system, and to avoid any simplification.
Let us therefore carry out a direct linear stability analysis of equation (1) in the thermodynamical limit (i.e. when  ¥ N ) for the case in which the distribution of natural frequencies is given by equation (3), by making use of a method similar to that of [9,19]. In the continuum limit, a density function ρ(θ, ω, t) can be introduced, such that ρ(θ, ω, t) dθ accounts for the fraction of oscillators of natural frequency ω whose phases are between θ and θ+dθ at time t. ρ satisfies the normalization condition for all ω and all t, and its evolution is governed by the continuity equation where υ(θ, ω, t) is the angular velocity. In its mean-field form, equation (1) can be rewritten as where it is clear that the different oscillators are only 'seeing' the mean-field quantities r and Ψ.
Let us now analyze the stability of the incoherent state ρ 0 (θ, ω, t)=1/(2π) (υ=ω) by linearizing the continuity equation in the limit of small coupling strengths. Following the methodology of [9], one obtains that the characteristic equation for the discrete eigenvalue λ (whose real part determines the stability of the incoherent state) is Equation (8) explicitly relates the coupling strength κ with the eigenvalue λ and, therefore, when the sign of the Re[λ] changes from negative to positive, the incoherent state loses its stability. Once the frequency distribution of equation (3) is inserted in equation (8), one obtains and In general, the eigenvalue λ is complex, i.e. λ=a + ib with , . In the following, we separately discuss the three possible cases: (i) a>0, (ii) a=0, and (iii) a<0.

(i) [a>0]
In this case, f 1 (ω) and f 2 (ω) have poles ω 1 =ω 0 −iΔ and ω 2 =−ω 0 −iΔ respectively, and thus the integral in equation (9) can be solved by conveniently choosing a contour line in the lower half complex plane. The result is  to determine the critical coupling strength for the synchronization transition, and one is presented with two possible situations: (a) Soft mode instability, for Δ<ω 0 (bimodal regime). The critical coupling strength is In this case, a pair of complex eigenvalues crosses the imaginary axis at the bifurcation point, which typically leads to limit-cycle oscillations. Equation (14) defines the half line reported (in solid magenta) in figure 1(b), and labeled as HB.
There is only one real eigenvalue crossing the origin along the real axis at the bifurcation point, and no periodic solutions exist after the synchronization transition. The critical coupling strength is which defines the semicircle reported (in blue and green) in figure 1(b) and labeled as TC (for transcritical).
(ii) [a=0, λ=ib ]. f 1 and f 2 have poles at ω 1 =ω 0 +iΔ and ω 2 =−ω 0 +iΔ (in the upper half complex plane), and ω 3 =−b (in the real axis). By choosing an appropriate contour line, integration of equation (9) gives: Note that κ is real as long as the rhs of the above equation is a purely imaginary number, which is only possible if both the lhs and rhs are equal to 0. This leads to b(b 2 +Δ 2 −ω 0 2 )=0 and k  ¥ f which is physically unreasonable. There are two eigenvalues in this circumstance. One is λ 2 =0, the other is l w = -D i 2 0 2 2 (with Δ<ω 0 , bimodal) that is purely imaginary.

Numerical results
To have a full confirmation of the bifurcations and transitions predicted from the linear stability analysis of equation (8) and reported in figure 1(b), we performed extensive numerical simulations 7 of the full Kuramoto model (equation (1)) along the lines Δ/ω 0 =0.36, 0.8, 0.92, 1.0, 1.16, and 1.96 that traverse the regions I, II, III, Line EA, IV, and V, marked in figure 1(b).
Each panel forming figure 2 corresponds to one of these cases, and shows the order parameter r for the forward (backward) synchronization transition as the coupling strength κ increases (decreases). In particular, figure 2(a) depicts a case in Region I (Δ<0.55ω 0 ) where the system bifurcates from the incoherent (I) state (green area) to a Bellerophon state (Bs) and then to a partially synchronous (PS) state (yellow area). Both transitions are continuous, and there is no hysteresis loop. The two points in the transition curve marked as A and B correspond to Bs (whose detailed description will be provided below). In Region II (0.55ω 0 <Δ<0.85ω 0 , reported in figure 2(b)), the transition  I Bs is continuous, while the transitions Bs  PS  Bs are first-order-like with hysteresis.
The scenario reported in figure 2(c) (describing the emerging dynamics in Region III, i.e. for 0.85ω 0 <Δ<ω 0 ) is almost identical to that reported in figure 2(b), except for the backward transition at which the system experiences the passage from PS to I. The hysteresis area shows bistability, and there are actually two separate regimes where the PS state coexists with the I and the Bellerophon states. As the frequency distribution becomes closer to a unimodal one (see the insets at the left hand side of the panels), Bellerophon states are no longer present and the system bifurcates directly from I to PS and from PS to I following a first-order -like transition with hysteresis (panels d and e belonging to Line EA and Region IV, respectively) or without hysteresis (panel f, region V). . This is an important clarification to previous works' astray conclusions that such a region of the phase diagram sustains standing waves characterized by two counter-rotating groups of oscillators whose frequencies are locked plus a group of desynchronized oscillators. The Bellerophon states, instead, are multiclustered states emerging in symmetric pairs as the coupling increases, such that oscillators in each pair of clusters rotate with the same average velocity but in opposite direction, and whose phases and instantaneous frequencies are not locked [3]. Figure 3 shows a complete characterization of the Bellerophon states marked as points A and B in figure 2(a). For each oscillator i in the ensemble, the figure reports the instantaneous phases θ i (panels a1 and b1) and frequencies q i (panels a2, b2), and the long-time averaged frequency q á ñ i (panels a3, b3) of each oscillator in the ensemble, as a function of the natural frequency ω i . One can clearly see a symmetric multicluster structure (2 clusters labeled as '1' and '−1' for point A, and 8 visible clusters ±1,±3,±5,±7 for point B) portrayed as plateaus in panels a3 and b3, coexisting with some oscillators whose behavior is incoherent.
Remarkably, such a collective organization implies the emergence of quantized traits: the plateaus (in the staircases of the time averaged frequency) are, indeed, quantized in the sense that the q á ñ i are odd multiples (±(2n−1), n=1, 2...) of the lowest clusterʼs frequency, Ω 1 . Figure 3(a4) shows the real and imaginary values of the order parameter (equation (2)) for the two populations of oscillators with positive and negative frequencies splitting in two (red and green) oval orbits, reflecting the oscillatory motion of the total order  figure 2(a). Snapshots of the instantaneous phase θ i (a1 and b1), the instantaneous speed q i (a2 and b2), and the average speed q á ñ i (a3 and b3) versus the natural frequencies ω i of the oscillators. (a4) Local order parameters of the two counter-rotating clusters (red oval for positive frequencies and green oval for negative), and global order parameter in blue. Insets: time evolution of the modulus r and the phase Ψ of the global order parameter. (b4) Enlargement of the average speeds of the coherent clusters in panel b3. They correspond to odd-numbered multiples of the principle frequency Ω 1 . Insets: local order parameters in the complex plane for clusters 1, 3, 5, and 7. contrarians), the conclusion is that the Bellerophon state is in fact a generic organization of globally coupled phase oscillators occurring at intermediate values of the coupling strength, not limited to specific dynamical model nor to special arrangements in the frequency distributions. Even more importantly, Bellerophon states are multi-clustered states where higher order of coherence are set among the nonidentical oscillators (only longterm average frequencies are locked in each cluster of oscillators).
Furthermore, our work shows that explosive synchronization is inherent in classical Kuramoto model with bimodal frequency distributions. Its occurrence does not require any special scheme, such as a dynamicstopology correlation [6], a frequency weighted coupling [8], or an adaptive coupling [11]. In our model, as Δ/ω 0 increases, the frequency distribution gradually changes from a bimodal to a unimodal one. As shown in figure 2, explosive synchronization does not occur for purely unimodal distributions, nor it occurs when the two peaks of the bimodal distributions are too separated. It only occurs, interestingly, in a limited range of the bimodal case where the two peaks are not too much separated (as it can be seen in figure 2). The present results raise an important question: what is the necessary condition for first-order synchronization transitions? In other words, what is the mechanism underlying such phenomena?
We emphasize that Bellerophon states are generic states of partial (or weak) coherence occurring in globally coupled nonidentical oscillators, when frequencies are widely distributed. Typically, such states occur in the regime where the control parameter is at an intermediate value. In such case, on the one hand the coupling is not strong enough to completely entrain all oscillators (synchronized state), and on the other hand it is large enough to achieve certain correlations among them. As we know, fully or strong synchronization sometimes turns out to be harmful in real situations, such as the spiral wave in human heart and strong synchronization of EEG in brain during epileptic seizure. Therefore, the further investigation on Bellerophon states will certainly enhance our understanding of collective behaviors towards real-world circumstances.