Alejandro Carballosa

Alberto P Muñuzuri

Stefano Boccaletti

Alessandro Torcini

Simona Olmi

Cluster States and π-Transition in the Kuramoto Model with Higher Order Interactions

Cluster states and π -transition in the Kuramoto model with higher order interactions

I. INTRODUCTION

The Kuramoto model [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF] represents the paradigmatic framework in which to investigate the synchronization phenomenon occurring when a large heterogeneous population of oscillators lock their phases at unison. The Kuramoto model has been successfully employed to recover the main synchronization features observable in many contexts, ranging from biological scenarios [START_REF] Winfree | [END_REF][3][4][5] to human behaviours [6]. Due to its success, the model has been thoroughly studied in many different variations and set-ups [7,8]. A definitive breakthrough in the studies dedicated to Kuramoto model is represented by the exact reduction methodology developed by Ott and Antonsen in 2008 [9]. The Ott-Antonsen (OA) approach allows us to rewrite, for Lorentzian distributed frequencies, the mean-field evolution of the network in terms of a complex macroscopic field, the so-called Kuramoto order parameter. This methodology has generated a renewal of interest in the field, allowing for low-dimensional investigation of phase oscillator networks.

The Kuramoto model has been firstly analyzed in a globally coupled configuration with an unimodal distribution of the natural frequencies. In this context, by varying the coupling strength, one observes a transition * corresponding author: simona.olmi@fi.isc.cnr.it from an asynchronous (AS) to a partially synchronized (PS) regime, where a finite fraction of oscillators are locked and rotate coherently with a common frequency. For a sufficiently strong coupling the system can eventually become fully synchronized (FS) with all oscillators phase locked [7]. Plenty of attention has been also devoted to the case where the frequency distribution is bimodal, i.e characterized by two symmetric peaks [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF]10]. In particular, for bimodal Lorentzian distributions, the phase diagram can be obtained analytically in the thermodynamic limit by employing the OA approach [11]. In particular in [11] it has been shown that, besides the previous mentioned regimes, there is also room for multistability and a limit-cycle solution of the order parameter, termed standing wave (SW) and characterized by two counter-rotating clusters with the same angular velocity. Travelling waves (TWs) are also expected in this set up whenever the bimodal distribution becomes asymmetric [11].The presence of hysteretic transitions [12] and other complex, time dependent states [13] has been also shown to be a consequence of introducing a bimodal distribution of frequencies.

Another generalization of the original Kuramoto model consisted in the modification of the coupling function by including either higher order harmonics [START_REF] Winfree | The geometry of biological time[END_REF][START_REF] Daido | [END_REF] or, more recently, by extending the pair-wise interaction to many body terms [START_REF] Battiston | Physics Reports Networks beyond pairwise interactions[END_REF]17]. For each extra harmonics appearing in the coupling term, the interaction of the single phase oscillator with all the others is mediated by an extra macroscopic field, termed Kuramoto-Daido order parameter [START_REF] Daido | [END_REF]18]. If the interaction is pair-wise, the Kuramoto-Daido order parameters appear linearly in the phase evolution equation, while many body interactions give rise to the emergence of non-linear combinations of these order parameters [START_REF] Battiston | Physics Reports Networks beyond pairwise interactions[END_REF]17,[19][20][21].

Higher harmonics, as well as many body terms, emerge naturally by considering the dynamics of N coupled dynamical systems in the proximity of a super-critical Hopf bifurcation [22] (or of a mean field Complex Ginzburg-Landau equation [23]), whenever such dynamics can be reduced to that of coupled phase oscillators, i.e. for sufficiently small coupling terms. In particular, at the lowest orders besides the usual Kuramoto-Sakaguchi model, a bi-harmonic coupling function as well as three body interactions terms always emerge. The latter ones can be either symmetric or asymmetric with respect to the reference oscillator, thus implying that its phase appears either as a higher order or a simple harmonic interaction term, respectively [START_REF] Battiston | Physics Reports Networks beyond pairwise interactions[END_REF]17]. The many body asymmetric interactions have been largely examined, because, in this case, the OA approach can be fruitfully employed to reduce the mean-field dynamics to the evolution of few macroscopic indicators, thus allowing for the derivation of analytic results [24].

The inclusion of higher harmonics in the coupling function gives rise to the emergence of a large number of coexisting stable or multi-stable phase-locked states at different locations on the unitary circle, termed cluster states [START_REF] Winfree | The geometry of biological time[END_REF][18][19][20]. A similar phenomenology has been also reported when considering only three body symmetric interactions; in such a case, the mean-field interaction is quadratic in the modulus of the Kuramoto order parameter R 1 and bi-harmonic in the single oscillator phase [19,21,25]. Therefore, bi-clusters in phase opposition are expected in this case. It should be noted that, if the bi-clusters are equally populated, the level of synchronization in the network, as measured by R 1 , will be exactly zero. As a matter of fact the bi-clusters always coexist with the asynchronous regime, which remains stable in the thermodynamic limit [19]. Finite size fluctuations may be responsible for the destabilization of the asynchronous regime and the emergence of asymmetrically populated clusters with an associated finite level of synchronization [21]. By considering as initial states the asymmetrically populated bi-clusters, one observes that they loose stability via an abrupt desynchronization towards the asynchronous regime, where the system remains forever [25]. The origin of this desynchronization transition has been recently related to a collective mode instability [26]. It is worth mentioning that the bi-clusters (CSs), here analyzed, are different from the previously mentioned TW solutions, since they are characterized by coexisting phase-locked clusters, as for the TWs, but where all oscillators are frequency locked.

Furthermore, in [24,27], it has been shown that, by considering both pairwise and higher order interactions in the system, the latter ones can stabilize synchronized states even for repulsive pair-wise coupling, a situation prohibited in the classical Kuramoto model. Furthermore, it was recently shown that, even with both pairwise and triadic coupling being repulsive, the two-cluster state can exist when the triadic strength is sufficiently large [28]. In social terms, this means that some level of agreement can exist in the network despite all the individuals being contrarians to the mainstream.

In this work we aim at studying the combined effect of introducing a bimodal frequency distribution in the Kuramoto model with the presence of higher order interactions, that can be either attractive or repulsive. In particular, by adapting the approach developed in [20] for a bi-harmonic coupling function to our set-up, we have performed a mean-field analysis based on the identification of self-consistent stationary solutions for the Kuramoto order parameter. The mean-field results, joined with direct numerical simulations, allow us to identify all the possible dynamical regimes and to characterize the transitions separating them, thus being able to obtain a quite detailed two dimensional phase diagram in terms of the pair-wise and higher order coupling terms.

This paper is structured as follows. In Section II we describe the model and the employed coherent indicators, as well as the simulation protocols. Section III is devoted to report the mean-field formulation and to the investigation and comprehension of the possible synchronization and de-synchronization transitions occurring for unimodal and bimodal Gaussian distributions of the natural frequencies, for attractive and repulsive higher-order interactions. The bimodal Kuramoto model with higher order interactions is further investigated in Section IV via numerical simulations. In particular, in this Section we report all the emergent dynamical regimes, as well as the phase diagrams obtained for fully synchronized and random initial states. The sub-section IV C is particularly devoted to the detailed characterization of the π-transition associated to the clustered regimes. A summary of the obtained results and a brief discussion are reported in Section V. A simplified version of the studied model depending on only one coupling parameter is studied and deeply analyzed in the Appendix.

II. METHODS

A. Models

In the original Kuramoto model [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF], the ensemble of N oscillators is coupled all-to-all and the evolution equation for a generic phase oscillator is given by

dθ i (t) dt = ω i + λ N N j sin (θ j -θ i ) i = 1, . . . , n (1)
with ω i being the natural frequency of the i-th oscillator and λ the coupling strength.

For globally coupled phase oscillators the evolution equation can be generally written in terms of complex mean-field variables, termed Kuramoto-Daido order parameters, which take into account the interaction with the rest of the network, specifically :

Z m = R m e jΨm = 1 N N i=1 e jmθi . (2)
For m = 1 one obtains the so-called Kuramoto order parameter, whose modulus R 1 measures the level of synchronization among the phase oscillators: R 1 ∝ 1/ √ N for an asynchronous state, while R 1 is finite (one) in a partially (fully) synchronized situation.

Since the coupling interaction in Eq. ( 1) is limited to the first harmonic and pair-wise, this can be rewritten simply in terms of Z 1 appearing linearly in the equation, as follows

dθ i (t) dt = ω i + R 1 sin(Ψ 1 -θ i ) i = 1, . . . , n . (3)
In the absence of interactions, it is expected that the phase of each oscillator evolves according to its natural frequency ω i . Consequently, it is natural to study how the system behaves as a function of the chosen frequency distribution g(ω). The most widely studied is a symmetric, unimodal distribution, peaked at some value ω 0 [7]. In this case, for increasing values of λ, the system passes from being asynchronous (AS), to partially synchronized (PS), and finally, for large enough values of λ, to a fully synchronized (FS) regime.

In this paper we consider an extension of the classical Kuramoto model by adding a term for many body interactions; specifically we will consider three body (or triadic) symmetric interactions. Therefore, the triadic interaction term will depend simultaneously on the phase differences between the reference oscillator i and all the other possible couples j and k, i.e. on (θ i -θ j ) and (θ i -θ k ). For both pair-wise and triadic interactions we will consider an all-to-all configuration, where all possible pairs and triangles are formed. Moreover we consider different weights for each interaction, being λ 1 the pairwise interaction strength and λ 2 the triadic interaction strength, namely:

dθ i (t) dt = ω i + λ 1 N N j sin (θ j -θ i ) + λ 2 2N 2 N j=1 N k=1 sin (θ j + θ k -2θ i ) . ( 4
)
This model can be rewritten in terms of the mean-field variable Z 1 as follows

dθ i (t) dt = ω i + λ 1 R 1 sin(Ψ 1 -θ i ) + λ 2 2 R 2 1 sin(2(Ψ 1 -θ i )) i = 1, . . . , N ; (5)
where, due to the triadic symmetric interaction, a quadratic term R 2 1 is now present, as well as a biharmonic interaction [19,21]. It should be noted that for asymmetric triadic interactions, as the one studied in [24], one would observe a single harmonic term for the reference oscillator phase joined to a nonlinear mean-field term of the type R 1 R 2 , therefore a quite different dynamics should be expected in such a case.

The natural frequencies {ω i } are usually assumed to be distributed as a bimodal Gaussian, which can be written as the superposition of two Gaussians with centers ±ω 0 and standard deviation ∆ as follows:

g B (ω) = 1 2 √ 2π∆ e -(ω+ω 0 ) 2 2∆ 2 + e -(ω-ω 0 ) 2 2∆ 2 . ( 6
)
In the following, the simulations will be performed by randomly generating {ω i } from the above distribution with the additional constraint that 50% of the oscillators are associated to the negatively (positively) peaked distribution, thus ensuring equally distributed oscillators. Moreover, we will also report some results for unimodal Gaussian distributed natural frequencies,

g U (ω) = 1 √ 2π∆ e -(ω) 2 2∆ 2
where the distribution is centered in zero and has standard deviation ∆.

B. Coherence indicators

In order to measure the level of synchronization of the oscillators, we will employ the modulus of the Kuramoto order parameter R 1 , as well as the modulus of the Kuramoto-Daido parameter R 2 that is particularly suited to reveal the presence of two phase clusters, as expected for the present model. However, while the R 2 parameter is indeed able to clearly discriminate two phase clusters in phase opposition, whenever the two clusters are located on the circle at an intermediate angle γ between 0 and π, this indicator shows a non monotonic dependence on the angle, as shown in Fig. 1. It is therefore desiderable to introduce a new order parameter monotonically varying with the angular separation of the two clusters, particularly suited to characterize bi-clusters :

R B = R (+) 1 + R (-) 1 2 -R 1 ; (7)
with

R (+) 1 (R (-)
1 ) being the modulus of the Kuramoto order parameter limited to the oscillators with positive (negative) natural frequencies. This indicator will be exactly zero whenever the oscillators are fully synchronized forming a single cluster, then it will vary smoothly with the angle γ between the two clusters, reaching the value one when γ = π, analogously to the indicator R 2 , as illustrated in Fig. 1 FIG. 1. Order parameters (b) against the inter-cluster angle γ (a). R1 is the Kuramoto order parameter, R2 is the order two Kuramoto-Daido parameter and RB is the clustering order parameter defined in eq. 7.

C. Simulation protocols

As the system shows many complex, multistable behaviors, we have examined the system evolution by considering many different initial conditions. These initial conditions are characterized by the different percentage η of initially phase-locked oscillators. Referring to η, we define the following two simulation protocols: protocol (I) refers to simulations started with η = 0, where all the oscillators have random initial phases; protocol (II) refers to an initial state where all the oscillators phases are set to the same value, therefore η = 1. Furthermore, protocols (I A) and (II A) will denote simulations started from η = 0.0 and η = 1.0 respectively, whenever one of the coupling parameters is varied quasi-adiabatically. In particular, once initialized the phases, the simulations are performed at a constant value of the control parameter, let us say λ 1 , for a certain duration T s , after discarding a transient T t . Then by considering the last configuration of the previous simulation as the new initial condition, the control parameter is increased (decreased) by a small amount ∆λ 1 and the simulation is repeated for the same transient and duration. This operation is repeated until the desired maximal (minimal) value of λ 1 is reached. In the next sections, the effect of 0 < η < 1 on the dynamical states will be also examined.

The network evolution equations (4) have been integrated by employing a 4th order Runge-Kutta scheme either with fixed time step dt = 0.01 -0.001 or via an adaptive Runge-Kutta-Fehlberg method. Regarding the terminology, in the following sections we will refer to the synchronizing (de-synchronizing) transition as the forward (backward) transition, where we mean that the transition occurs by increasing (decreasing) the coupling. We will also refer to the instantaneous angular velocity of each oscillator θi and to its time average by ⟨ θi ⟩.

III. MEAN-FIELD ANALYSIS

In this Section we will analyze the model in terms of a self-consistent mean-field approach specifically inspired by the analysis performed in [20] for a bi-harmonic coupling function, but essentially analogous to that employed by Kuramoto in his seminal work [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF]. In particular, by assuming that the system is partially synchronized, i.e. that R 1 is finite and the clustered oscillators are uniformly rotating with a common angular velocity ω c , the mean-field equation ( 5) can be rewritten as follows in terms of the rescaled phase ϕ = (θ

-Ψ 1 ) ∈ [-π : π], dϕ dt = Ω -λ 1 R 1 sin(ϕ) - λ 2 2 R 2 1 sin(2ϕ); ( 8
)
where Ψ 1 = ω c t is the global phase and ω = Ω + ω c is the natural frequency of the phase oscillator. The above equation ( 8) can be interpreted as the dynamical equation for an overdamped oscillator moving in a potential landscape given by

V (ϕ) = -Ωϕ -λ 1 R 1 cos(ϕ) - λ 2 R 2 1 4 cos(2ϕ) . (9)
For even frequency distributions and frequency locked clusters we can assume that Ω = 0; in such a case the potential is symmetrical to a change of the phase sign, i.e. V (ϕ) = V (-ϕ). By following [20], we can affirm that the potential always reveals two coexisting minima whenever

|u| = λ 1 λ 2 R 1 < 1 ; (10)
otherwise it displays a single minimum, analogously to the usual Kuramoto model. We will mostly focus on the region where the two minima coexist, by separately analysing the cases for positive and negative λ 2 -values. For |u| < 1 the potential presents three extrema located at ϕ = 0 ; ϕ = ±π and ϕ = ± φ = ± arccos(-u) , (11) since the phases ±π coincide.

As shown in Fig. 2 (a-b), for λ 2 > 0, the potential displays two minima of different heights in ϕ = 0 and ϕ = ±π, and two maxima in ± φ. Therefore, depending on the initial conditions and on the protocol used for the simulations, one or two clusters may emerge in the system. In particular we notice that, for increasingly larger positive (negative) values of λ 1 , the minimum corresponding to the cluster in ϕ = 0 (ϕ = ±π) becomes deeper and therefore more stable, as evident when looking at V (ϕ) (black solid lines) reported in Fig. 2 (a) and (b). Furthermore, for λ 1 ≤ -λ 2 R 1 (λ 1 ≥ λ 2 R 1 ) the minimum in ϕ = 0 (ϕ = ±π) disappears and the potential shows an unique minimum. For a negative λ 2 -value the potential presents two maxima at ϕ = 0 and ±π and two symmetric minima of same height at ± φ, as shown in Fig. 2 (c) and (d) (black solid lines).

The locked phases correspond to stationary solutions of (8) given by

Ω = λ 1 R 1 sin(ϕ) + λ 2 2 R 2 1 sin(2ϕ) = Y (ϕ, R 1 ) . ( 12
)
These solutions will be stable whenever ∂Y /∂ϕ = V " (ϕ) > 0, i.e. when the potential is concave in proximity of some minimum. The curves Ω = Y (ϕ, R 1 ) are shown in Fig. 2: solid (dashed) red lines denote the stable (unstable) regions. The function Ω = Y (ϕ, R 1 ) presents extrema in correspondence of the angles

±α = ± arccos -u + √ u 2 + 8 4 (13)
and

±β = ± arccos -u - √ u 2 + 8 4 . ( 14
)
We have denoted the absolute value of the natural frequencies at the extrema as Ω 1 = |Ω(±α)| and Ω 2 = |Ω(±β)| and the minimal frequency among the two as Ω min .

Looking at the curves Ω = Y (ϕ, R 1 ) shown in Fig. 2 for |u| < 1, it is evident that the system will always display two stable branches, which can be identified as follows:

• For λ 2 > 0 : the first branch corresponds to phases (frequencies) in the interval

Φ 1 ∈ [-α : α] ([-Ω 1 : Ω 1 ]
) and the second one to phases (frequencies) in the interval

Φ 2 ∈ [-π : -β] ∪ [β : π] ([-Ω 2 : Ω 2 ]);
• For λ 2 < 0 : the first branch corresponds to phases (frequencies) in the interval

Φ 1 ∈ [-α : -β] (-Ω 1 : Ω 2 ]
) and the second one to phases (frequencies) in the interval

Φ 2 ∈ [β : α] ([-Ω 2 : Ω 1 ]).
Moreover it is worth noticing that there is a region of coexistence between the two stable branches for -Ω min ≤ Ω ≤ Ω min , where oscillators with the same frequency Ω can be locked to 2 different phases (see red solid curves in Fig. 2). The self-consistent equation for the order parameter in the thermodynamic limit is analogous to the one for the Kuramoto model [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF] and it can be written as

R 1 e jΨ1 = e jΨ1
dϕ dΩ ρ(ϕ|Ω) e jϕ g(Ω) [START_REF] Daido | [END_REF] where g(Ω) is the distribution of the natural frequencies and ρ(ϕ|Ω) the conditional probability distribution function characterizing the oscillator population. The oscillators can be divided in drifting or locked depending on their natural frequency. By following [21] we denote with Y max (Y min ) the maximum (minimum) of the function Y (ϕ, R 1 ) (12). According to this notation, the locked (drifting) oscillators will have

Y min ≤ Ω ≤ Y max (Ω > Y max and Ω < Y min ).
The probability distribution ρ d for the drifting oscillators is proportional to the inverse of their phase velocity [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF]:

ρ d (ϕ|Ω) = C |Ω -Y (ϕ, R 1 )| ; C = 1 π -π dϕ|Ω -Y | . ( 16
)
The distribution ρ l of the locked oscillators in the region of coexistence can be written as (17) where 0 ≤ σ ≤ 1 is the redistribution factor of the oscillators among the two branches. Outside the coexistence region, the distribution ρ l is simply given by

ρ l (ϕ|Ω) = σδ(ϕ -Φ 1 (Ω)) + (1 -σ)δ(ϕ -Φ 2 (Ω)) ;
ρ l (ϕ|Ω) = δ(ϕ -Φ x (Ω)) , x = 1, 2 , (18)
depending if the oscillator characterized by the phase ϕ is locked to the first or second branch. For λ 2 > 0, by assuming the frequency distribution to be even (g(Ω) = g(-Ω)) and by taking into account the symmetry of the integrand, the expression for the Kuramoto order parameter takes the form

R 1 = 2 α 0 dY (ϕ)S 1 (Y (ϕ))g(Y (ϕ)) cos(ϕ) + π β dY (ϕ)S 2 (Y (ϕ))g(Y (ϕ)) cos(ϕ) + U π -π dΩdϕ Cg(Ω) cos(ϕ) |Ω -Y (ϕ, R 1 )| , (19)
where S 1 = σ and S 2 = 1 -σ in the coexistence regions, S 1 = S 2 = 1 outside the coexistence region and

dY (ϕ) = [λ 1 R 1 cos(ϕ) + λ 2 R 2 1 cos(2ϕ)]dϕ . (20)
The integration region U = (-∞, Y min ) ∪ (Y max , +∞) that appears in the third integral in (19), is the frequency domain where the oscillators are drifting. For λ 2 < 0 the expression for R 1 becomes

R 1 = -α -β dY (ϕ)S 1 (Y (ϕ))g(Y (ϕ)) cos(ϕ) + α β dY (ϕ)S 2 (Y (ϕ))g(Y (ϕ)) cos(ϕ) + U π -π dΩdϕ Cg(Ω) cos(ϕ) |Ω -Y (ϕ, R 1 )| . (21)
In this case, due to the symmetry of the potential wells (see Fig. 2 (c-d)), the value of R 1 does not depend on the value of the redistribution factor σ, as we have verified.

Therefore we can set S 1 = S 2 = 1/2 inside the coexistence region and S 1 = S 2 = 1 outside, without any loss of generality. It should be noted that, for the parameter values considered in this paper, the contribution of the drifting oscillators to ( 19) and ( 21) is usually negligible. The self-consistent approach allows to find coherent solutions of the system in the thermodynamic limit corresponding to finite values of the Kuramoto order parameter, but not to determine their stability properties.

A. Unimodal Distribution of the Natural Frequencies

Let us first consider the case where natural frequencies are distributed according to an unimodal Gaussian distribution g U (ω) with ∆ = 0.25. In this case we have characterized the synchronization transition in terms of the order parameter R 1 , evaluated by performing simulations with protocols (I A) and (II A). In the following, if not stated otherwise, we avoid explicitly indicating the time average of R 1 (t) by introducing a new symbol and we intend for R 1 its time averaged value. Moreover, we fix the parameter λ 2 to a positive or negative value and we vary quasi-adiabatically λ 1 .

Positive λ2

The simulation results for λ 2 = 15 and N = 10000 are shown in Fig. 3 (a). For protocol (I A) the partially synchronized state emerges for λ (AS) 1 ≃ 0.4, corresponding to the destabilization of the asynchronous regime of the usual Kuramoto model [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF], i.e. λ (AS) 1 = 2 πg(0) . In the present case the higher order interactions are proportional to R 2 1 , as shown in ( 8), therefore they cannot affect the stability of the incoherent regime, which is usually analysed in proximity of the transition to coherence via an expansion of the model at the first order in R 1 [29]. For this reason we obtain a critical λ (AS) 1

value corresponding to the one for the usual Kuramoto model. Moreover the system becomes fully synchronized (R 1 = 1) for λ 1 ≃ 13.7 and the approach to the fully synchronized regime occurs via a series of plateaus where R 1 remains constant when varying λ 1 . The partially synchronized regime is always characterized by two phase clusters differently populated, as shown in Figs. 3 (b-e), which correspond to the two minima of the potential at ϕ = 0 and ϕ = ±π displayed in Fig. 2 (a). By increasing λ 1 the main cluster, corresponding to the minimum at ϕ = 0 of V (ϕ), becomes more and more populated at the expenses of the secondary cluster. Finally only this cluster remains. This can be clearly appreciated by measuring the redistribution factor σ directly from the simulations, where σ corresponds to the fraction of oscillators present in the minimum located at ϕ = 0 (see the inset of Fig. 3 (a)). In particular we observe that σ ≥ 0.5, meaning that the cluster at ϕ = 0 is always more populated than that at ϕ = ±π. For protocol (II A), the system is initialized with all the phases equal to θ i = 0. Therefore the oscillators are all located in the minimum ϕ = 0 of the potential and they remain fully synchronized until λ

(P S) 1 (σ = 1) ≡ λ (F S) 1
≃ -11.8, when they abruptly desynchronize. In particular, λ (P S) 1 (σ) indicates the critical value at which the coherent regime characterized by σN ((1 -σ)N ) oscillators in the minimum ϕ = 0 (ϕ = ±π) looses stability. A peculiar difference with the usual Kuramoto model is that, as soon as λ 1 > λ (AS) 1

, two phase clusters emerge, corresponding to the two minima in the potential shown in Fig. 2 (a), at an angular distance of π. This is confirmed by the snapshots of the phase oscillators presented in Figs. 3 (b-e), which always reveal two clusters in the partially synchronized regime located at a distance π in phase.

To get some more understanding we estimated the selfconsistent solutions (19) for different σ-values, thus obtaining the evolution of R 1 as a function of λ 1 , shown in Fig. 3 (a) for a few σ-values. These curves for σ > 0.5 have all a similar shape. The partially synchronized state emerge abruptly via a bifurcation, that we can identify as a saddle-node, at a value λ (P S) 1 (σ) < 0, with a finite order parameter value

R 1 (σ) = 2σ -1 ,
corresponding to a 2 cluster state with σN oscillators in ϕ = 0 and (1 -σ)N oscillators in ϕ = ±π. The critical value λ

(P S) 1 (σ) should be larger than λ 1 = -λ 2 R 1 (σ),
where the minimum at ϕ = 0 in V (ϕ) emerges. By increasing λ 1 , the value of the order parameter remains equal to R 1 (σ) within a certain interval, giving rise to a plateau and finally it approaches R 1 = 1 for λ 1 = λ 2 , corresponding to the value where the minimum in V (ϕ) at ϕ = ±π disappears and all the oscillators are fully synchronized at ϕ = 0. We have estimated a few selfconsistent curves corresponding to σ-values actually measured during the finite size quasi-adiabatic simulation. The plateaus observed during the simulation performed via protocol (I A) coincide with those obtained by the self-consistent approach for the corresponding σ-value, as clearly evident in Fig. 3 (a). Furthermore, the selfconsistent results for σ = 1 perfectly coincide with the simulations obtained with protocol (II A) and the desynchronization transition occurs at λ (F S) 1 = -11.8 as in the finite size simulations.

For the case σ = 0.5, corresponding to equally populated clusters, the self-consistent approach reveals that this state emerges at λ (P S) 1 (σ = 0.5) ≃ 0.25 with an order parameter value that is quite small (R 1 ≃ 0.04).

Moreover, according to the self-consistent approach for σ = 0.5, no plateau is observable (see the cyan line in Fig. 3 (a)). Similar results have been found for σ < 0.5 with curves R 1 = R 1 (λ 1 ) lying slightly below the one obtained for σ = 0.5.

To better clarify the situation, we expect the incoherent regime to become unstable at λ AS 1 = 2 πg(0) ≃ 0.4 in the thermodynamic limit. However there is a wide region of multistability in the interval λ 1 ∈ [λ 11.8; 15], where two cluster state solutions coexist, although populated with a different percentage σ of oscillators. Therefore we expect the synchronization transition to be generically hysteretic, and to observe a scenario somehow similar to that already observed for the Kuramoto model with inertia [30]. The hysteretic nature of the transition is evident by considering the simulation performed by employing protocol (II A).

(F S) 1 ; λ 2 ] = [-

Negative λ2

In the case λ 2 = -100 and N = 10000 oscillators, by following protocol (I A), the asynchronous regime looses stability via a continuous transition in the proximity of λ

(AS) 1 = 2
πg(0) , as expected (see blue circles in Fig. 4 (a)). In the partially synchronized regime two equally populated clusters are always observable and the increase of λ 1 renders the two clusters more synchronized, i.e. the oscillators within each cluster tend to have more and more similar phases, as shown in Fig. 4 (b). This behaviour was expected from the shape of the potential V (ϕ), that in this case presents two minima of equal heights as observable in Fig. 2 (b).

By following protocol (II A) we observe that, when decreasing λ 1 , the Kuramoto order parameter R 1 de- creases smoothly with the control parameter and then, at λ (π) 1 ≃ 3, it presents an abrupt jump to a higher value (see maroon circles in Fig. 4 (a)). Here we indicate with λ (π) 1

the critical value at which a discontinous transition is observable when performing protocol (II A), where the system jumps from a lower value of R 1 to a higher one. Once performed the jump, by decreasing λ 1 , the value of R 1 continuously diminishes towards zero and the asynchronous regime is smoothly achieved for λ

(P S) 1 ≃ -2.1.
This discontinuity in R 1 has been previously reported in [13], but its origin is still unclear. An important aspect to notice is that, before the jump, the oscillators are organised in two equally populated clusters, characterized by either positive or negative natural frequencies while, after the jump, the two clusters still present mainly positive (negative) natural frequencies, but together with a small group of oscillators with opposite ω i sign (as visible in Fig. 4 (c)). According to protocol (II A), the oscillators are all initialized with the same phase θ i = 0, that now corresponds to a local maximum of the potential, see Fig. 2 (c-d). Therefore, depending on the sign of their natural frequencies, the oscillators split in 2 clusters characterized, each one, by the same sign of ω i . In this way, the oscillators migrate towards one of the two minima of V (ϕ) located in correspondence of the positive or the negative phases.

Also in the present case we can estimate the possible values attained by R 1 as a function of λ 1 in the mean field limit corresponding to protocol (I A), by solving the selfconsistent equation Eq. (21). As already mentioned, in this case, due to the symmetries of the potential V (ϕ), we do not observe any dependence on σ, therefore we set σ = 0.5 in (21). The results are shown as a cyan line in Fig. 4 (a): the self-consistent solution for R 1 reveals two successive saddle-node bifurcations at λ 1 ≃ 2.41 and λ 2 ≃ 1.38. However a good agreement with the numerical data is observable for λ 1 ≥ 2.0 only. This is due to the fact that the phase clusters emerge only above λ 1 ≥ 2.0 while, below such a value, the system do not reveal any clear cluster.

To estimate the results corresponding to protocol (II A) via the self-consistent approach, we solved Eq. ( 21) by integrating only on positive (negative) natural frequencies depending if the oscillators are in the cluster with positive (negative) phases. Once more the agreement with the simulations is quite good from large λ 1 values down to the discontinuous transition value λ (π) 1 . The self-consistent solution (magenta line in Fig. 4 (a)) reveals a sharp transition at λ (π) 1 ≃ 3.0, that we have termed π-transition for reasons that will become evident in the following. Since we have no information about the stability of the solutions themselves, we can argue that the solution is stable (unstable) for λ 1 > λ

(π) 1 (λ 1 < λ (π) 1 )
and that the 2 solutions annihilate in a saddle-node bifurcation at λ (π)

1 . The destabilization of the phase cluster states characterized by either positive or negative ω i can be understood by following the evolution in time of the order parameters R 1 and R B at λ 1 = 3.0, as shown in the inset of Fig. 4 (c). The value of R 1 (R B ) is initially constant then, abruptely, R 1 vanishes (R B reaches one) and, afterwards, the value of R 1 (R B ) becomes larger (smaller). This behaviour can be interpreted as follows: the two cluster states are initially separated by an angle γ < π, then they recede and γ reaches exactly π, corresponding to an anti-phase situation. The vanishing of R 1 implies a vanishing of the potential V (ϕ), therefore the 2 clusters destabilize and reform with a smaller angle γ between them. Moreover they reform with a different composition inside each cluster, which now includes oscillators with positive and negative natural frequencies. Therefore, the π-transition was at the origin of the discontinuity in R 1 reported in [13], and in the following, we will examine this transition in more detail.

The π-transition: dependence on η

The same analysis to the one reported in the previous sub-section for protocol (IIA), has been performed here by varying the different percentage η of initially phaselocked oscillators. In particular we have calculated the time average value of the order parameter R 1 by follow-ing protocol (II A), here adapted to take into account different initial states where only a percentage η of the oscillators' phases are set to the same value. The results of these simulations, shown in Fig. 5 (a), reveal that, irrespectively of the chosen η value, the order parameter R 1 decreases smoothly with the control parameter and then, almost at the same critical value, it presents an abrupt jump to a higher value. By decreasing η the jump occurs at smaller and smaller λ (π) -values while, at the same time, the entity of the jump in R 1 decreases. Finally the asynchronous regime is smoothly achieved at a λ (P S) 1 value that increases for decreasing η. From the numerical data, we can infer that, in the vanishing η limit, the curve will become continuous and coincident with the one obtained by following the protocol (I A), shown as a dashed cyan line in Fig. 5 (a). An interesting aspects is that, for η < 1, the two cluster states are always composed by oscillators with positive and negative {ω i } before and after the transition, as shown in Fig. 5 (b-e). The main difference that we observe in the 2 clusters before and after the π-transition is that, once they reform, the distribution of the {ω i } in each cluster does not cover any more a compact interval, but there are holes in the support of the distributions.

It should be remarked that all the discontinuous transitions observed for finite η are associated to the fact that the two clusters reach an anti-phase configuration, as we have verified by inspecting the time evolution of R 1 and R B at the transition value.

B. Bimodal Distribution of the Natural Frequencies

We will now examine the bimodal case, where the natural frequencies are distributed by following g B (ω) with ω 0 = 1 and ∆ = 0.25. Also in this case we will compare numerical simulation results for protocols (I A) and (II A) with mean-field results obtained through the selfconsistent approach.

Positive λ2

Due to the bimodal frequency distribution, the asynchronous regime is characterized by two groups of oscillators rotating in opposite directions. On one hand, in the unimodal framework, the oscillators synchronize around the common frequency Ω = 0, which corresponds to the most probable natural frequency value, and form a phase cluster that smoothly grows with λ 1 . On the other hand, here we pass from 2 groups of oscillators rotating with opposite frequencies around ±ω 0 to locked oscillators with Ω = 0 divided in two phase clusters. Therefore, we expect the partial synchronization to emerge via a discontinuous transition for a coupling λ (AS) 1 larger than in the unimodal case [12]. Indeed, as shown in Fig. 6 (a) for λ 2 = 15, the simulations for N = 10000 and protocol (I A) show a discontinuous transition at λ (AS) 1 ≃ 0.7 (blue circles in Fig. 6 (a)). The 2 cluster states formed after the transition find themselves at a distance π, as in the unimodal case. However they are now quite asymmetric, as evident from the value of σ, measured during the simulations and shown in the inset of Fig. 6 (a). By increasing λ 1 the main cluster at ϕ = 0 becomes more and more populated and, for λ 1 ≥ λ 2 = 15, all the oscillators reside in such a cluster (see Fig. 6 (b-e)). Therefore, once the partially synchronized regime is emerged, the scenario turns out to be quite similar to that observed for the unimodal frequency distribution with λ 2 > 0. Indeed, also in this case a small plateau is observable for σ = 0.876.

The self-consistent approach estimated by employing Eq. ( 19) for σ > 0.5 reveals, in complete analogy with the unimodal case, that a partially synchronized regime emerges via a saddle-node bifurcation at λ (P S) 1 (σ), characterized by a finite value of the order parameter R 1 (σ) = 2σ -1. In this regard we have reported a few of these curves in Fig. 6 (a). The solution shown as magenta diamonds, corresponds to σ = 0.768, which is the minimal σ value measured via the simulations in the PS regime. The magenta solution reproduces quite well the numerical simulations for λ 1 > 0.7 up to the complete synchronization, apart for the plateau, which is instead captured by the self-consistent solution for σ = 0.876 (orange symbols). Finally, the curve for σ = 1 (green symbols) perfectly reproduces the simulations obtained by following protocol (II A) (maroon circles) at the associated de-synchronization transition. At λ 1 ≃ -9 the evaporation of few oscillators towards the second minimum results in the reduction of R 1 , with a consequent abrupt system de-synchronization, as shown in Fig. 6 (e).

The self-consistent solution for σ = 0.5 (cyan sym-bols) shows the formation of a PS state via a saddlenode bifurcation for λ 1 > λ (P S) 1 (σ = 0.5) ≃ 0.6 with a finite value of the order parameter R 1 ≃ 0.36. The value λ (P S) 1 (σ = 0.5) can represent a reasonable approximation of the numerically observed transition. As in the unimodal case, we have a wide interval of λ 1 -values where we can have coexistence of incoherent and partially synchronized states. Moreover, also in this case, hysteretic synchronization transitions are observable.

From the analysis reported in [11], one expects the transition from the asynchronous to the partially synchronous regime to be mediated by the emergence of SWs. However, for the present choice of parameters, we do not observe stable SWs. Indeed for ω 0 = ±1 and ∆ = 0.25, SWs are observable only in the range λ 2 ∈ [-11.7; 2], as we will discuss in the following.

Negative λ2

For λ 2 = -100, we observe that, at variance with the unimodal case, the transition from the asynchronous to the partially synchronized state is discontinuous. In more details, the network simulations for protocol (I A) show that the transition is steep but occurring in a finite interval λ (AS) 1 ∈ [0.6 : 1.4], probably due to finite size and finite time effects. The self-consistent estimation of the order parameter R 1 is obtained by solving Eq. ( 21). Analogously to the unimodal case we observe that, due to the symmetric shape of the potential, the integral (21) does not depend on σ, therefore we fix σ = 0.5. From the numerical simulations the PS emerges at λ (P S) 1 ≃ 3.8. Despite this time the two cluster state is already observable immediately after the transition, the mean-field results and numerical simulations are in reasonable agreement only for λ 1 > 5, as shown in Fig. 7 (b). The analysis of the network simulations performed following protocol (II A) reveals that (as in the unimodal case) oscillators are in a two cluster state characterized by only positive (negative) natural frequencies (see Fig. 7 (c)). By estimating the self-consistent integral in (21), limited to either positive or negative frequencies, we obtain two branches of solutions that merge in a saddlenode bifurcation at λ (π) 1 ≃ 8.80 (magenta line in Fig. 7 (a)). From the direct simulations, by decreasing λ 1 , we observe a very good agreement with the mean-field results down to the saddle-node bifurcation. As shown in the inset of Fig. 7 (c), where the time evolutions of R 1 and R B are reported at λ 1 = 8.8, R 1 (R B ) decreases (increases) in time until it reaches the zero (one) value corresponding to two clusters in phase opposition (i.e. with an angular distance γ = π). Immediately after, the clusters rearrange at a different angle distance γ < π. Each cluster is, in this case, still populated by oscillators with all positive or all negative natural frequencies, but the system is now characterized by a higher level of synchronization R 1 than before the transition. Therefore we can affirm that this is another example of π-transition, analogous to the one reported for the unimodal distribution. By further decreasing λ 1 , the system becomes less and less synchronized until it becomes completely asynchronous at λ (P S) 1 ≃ -12, due to an abrupt transition. Unfortunately, with the employed self-consistent approach we cannot predict the new clusterized state appearing below λ (π) 1 for protocol (II A). We leave this analysis for future work. It is important to notice that, for larger λ 2 ≥ -35, the state formed after the saddle-node bifurcation is made by clusters containing mostly oscillators with positive (negative) natural frequencies plus a few oscillators with negative (positive) ω i , analogously to what observed for positive λ 2 .

IV. NUMERICAL INVESTIGATIONS

This Section will be devoted to the numerical analysis of the dynamical regimes observable for the model ( 4) by considering mostly natural frequencies bimodally distributed as g B (ω).

A. Dynamical regimes

In this sub-section we will characterize all the states identified in our model by direct numerical simulations, apart the fully synchronized regime (FS) that corresponds to a group of oscillators locked to the same phase and traveling with the same angular velocity, for which R 1 = 1. The other observed states are presented in Fig. 8 and can be classified as follows:

• Asynchronous State (AS): as shown in panel (a), this state presents uniformly distributed phases, while the oscillators are drifting, driven by their natural frequencies (e). The order parameters R 1 and R B are essentially zero apart finite size fluctuations (i).

• Standing Waves (SWs): these solutions correspond to two clusters of oscillators (b) moving with opposite angular velocity (f), thus inducing oscillations in the order parameters, as shown in panel (j).

• Bimodal Clusters (BCs): these clusters are made up of two oscillators groups located in correspondence to the peaks of the bimodal distribution. Even though each cluster is locked to a different phase value (c), they are rotating with the same angular velocity, that, for symmetrically populated clusters, is zero, as shown in panel (g). The value of R 1 is tipically finite and not zero, since the two clusters are generically not located in phase opposition, while R B has a finite value that depends on the phase shift γ between the clusters (k). As we have seen in the previous Section, these solutions emerge by performing simulations with protocol (II) for negative λ 2 . In addition to this, for not too negative λ 2 (e. g. λ 2 -30) and λ 1 < λ

(π)
1 , are observable clusters composed not only of oscillators associated to one peak of the distribution, but also a few associated to the other peak.

• Uniform Clusters (UCs): they are constituted by two groups of oscillators phase locked at two different phases, where the oscillators of each cluster can have both positive and negative natural frequencies (d). Furthermore the two groups are usually not equally populated and this can give rise to a finite common angular velocity (h). Finally R B is essentially zero, since in this case it is not possible to identify two clusters on the basis of their natural frequencies, while R 1 is usually small but not zero since the angle γ between the clusters is tipically not π (l). These solutions usually emerge by performing simulations with protocol (I), as seen in the previous Section.

The typical solutions that emerge due to the bimodal Gaussian distribution are the SWs, observable also in the usual Kuramoto model in the absence of many body interactions [11], while the UCs and BCs have been already observed for unimodal distributions in presence either of a second harmonic or a three body interaction with a sufficiently negative λ 2 .

B. Phase Diagrams for Synchronized and Random Initial States

We will now numerically investigate the phase diagram of the model ( 4) as a function of the coupling parameters λ 1 and λ 2 , for bimodal Gaussian distributions centered at values ±ω 0 with a fixed distribution width ∆ = 0.25.

In particular, we have compared the results for the unimodal case (corresponding to ω 0 = 0) with two bimodal cases (centered in ω 0 = ±1 and ω 0 = ±3 respectively), where the overlap of the two Gaussians is definitely negligible for the chosen ∆-value. Both these distributions satisfy the criterion reported in [11] to discriminate unimodal from bimodal distributions: i.e. ω 0 > ∆/ √ 3. In order to characterize the macroscopic dynamics of the network we consider the modulus R 1 of the Kuramoto order parameter averaged in time and we measure the level of synchronization in the system by starting with two different sets of initial conditions corresponding to protocol (I) and (II), respectively.

As a general result, from Fig. 9 it is quite evident that asynchronous states (ASs) are observable for sufficiently negative λ 1 , for both unimodal and bimodal distributions, while large positive λ 1 values favour the emergence of fully synchronized regimes (FS). Furthermore, performing simulations with protocol (I) or (II) leads to large regions of the phase diagram in the plane (λ 1 , λ 2 ) where coexisting solutions are observable.

For protocol (I) (upper row in Fig. 9) we observe that in general, for the unimodal frequency distribution, we have a transition from an AS regime towards a PS one via the emergence of UCs that occurs for λ

(AS) 1 = 2
πg(0) . From the mean-field analysis, for sufficiently positive λ 1 , a FS regime corresponding to all oscillators in the same minimum of the potential V (ϕ) at ϕ = 0 is expected to occur for λ

(F S) 2
= λ 1 : this prediction is well satisfied for the unimodal distribution (dashed white line in panel (a)). For bimodal distributions we observe that λ 8) for ω 0 = 1 (ω 0 = 3). Therefore, the transition towards the FS regime occurs when the potential V (ϕ) still displays 2 minima. For negative λ 2 the transition fo FS become steeper. Furthermore, for bimodal frequency distributions, the emergence of synchronized states is now mediated by SWs, observable for not too negative, slightly positive λ 2 and in regions whose size increases for increasing ω 0 values, as expected from [10]. Moreover, the UCs observable for sufficiently positive (negative) λ 2 are indeed different, since one has not equally populated (equally populated) UCs due to the shape of the potential V (ϕ), as discussed in detail in Section III.

(F S) 2 ≃ 1.1λ 1 + C, with C = 2 (C =
For protocol (II) (lower row in Fig. 9) we observe 2 main differences with respect to protocol (I) simulations. The first one is detectable for positive λ 2 : the FS regime extends to quite negative λ 1 values where it abruptly de-synchronizes for values (λ

(F S) 2 , λ (F S) 1
), in agreement with the mean-field results, as shown in the previous section. The numerical simulations show a linear dependence of the critical λ

(F S) 2
value for the destabilization of the FS state versus the corresponding λ

(F S) 1 , in particu- lar λ (F S) 2 = sλ (F S) 1 + C with s ≃ 1.1 -1.
2 increasing with the value of ω 0 (white dashed lines in Fig. 9).

The second difference concerns the emergence of BCs at sufficiently negative λ 2 . As already discussed in detail in Section III, these states emerge perfectly symmetric at large λ 1 and destabilize at some smaller λ (π) 1 value of the parameter. The π-transition lines are clearly visible in the panels (d,e,f), at sufficiently negative λ 2 and positive λ 1 , as a discontinous variation in R 1 leading from a lower to a higher value of the order parameter by decreasing λ 1 . This is due to the fact that the BCs find themselves in phase opposition at the π-transition (i.e. with an angular distance γ = π). Immediately after, the clusters rearrange at a different angle distance in a configuration characterized by a higher level of synchrony that can eventually include a few oscillators belonging to the other peak of the bimodal distribution.

C. Clustered regimes

To better characterize the two coexisting cluster states BCs and UCs, we have investigated the dependence of the average angular velocity W = 1 N N j=1 ⟨ θj ⟩ on the distance of the peaks from the center of the frequency distribution. The BCs emerge when starting with fully synchronized initial conditions (protocol (II)), as shown in Fig. 10 (a). In this case we observe that W becomes essentially zero whenever ω 0 ≥ 1, i.e. whenever the 2 peaks are sufficiently separated. On the other hand, for uniformly distributed initial phases (corresponding to η = 0 and protocol (I)), the system ends up in UCs and W exhibits large oscillations depending on the initial conditions. Moreover W remains finite even for large ω 0 ≃ 6 × ∆, see Fig. 10 (b). This behaviour finds its explanation in the fact that, for η = 1, for a sufficiently large separation between the peaks, one observes the emergence of two equally populated clusters with opposite natural frequencies while, for random initial conditions, the oscillators belong to one of the two clusters irrespectively of the sign of their natural frequency. This imbalance is responsible for a finite common angular velocity W for all the oscillators. However, in both cases, all the oscillators are frequency locked, i.e. they can move coherently all together around the circle or they can be at rest, but these two situations can be considered as identical when examined in a comoving reference frame moving with velocity W .

To further explore how the cluster states depend on the initial conditions, we perform several simulations with different initial degree of the phase-locking η, for definitely negative λ 2 . In this region, for λ 1 = 10, we have previously identified UCs for protocol (I) (η = 0) and BCs for protocol (II) (η = 1). As shown in Fig. 11 (a), by sweeping η from 0 to 1, the clustering order parameter R B (Eq. ( 7)) grows from zero to one. Moreover this grow does not depend particularly on λ 2 ∈ [-1000 : -100], thus indicating that we pass smoothly from a situation where the phase clusters are not associated to the frequency distribution of the ω i (as observable for the UCs, see panel (d) in Fig. 8) to a situation where each cluster state involves only oscillators with positive or negative ω i (panel (c) in Fig. 8). By increasing the level of the initial synchronization of the phase oscillators η, we favour the emergence of phase clusters consisting of oscillators that are in prevalence associated to positive or negative ω i .

Let us now consider negative λ 1 values. For negative λ 1 = -10 we have previously identified an AS (BC) regime for η = 0 (η = 1). By varying η we expect to find a transition between these two regimes, that indeed occurs for the considered negative λ 2 values, as show in Fig. 11 (b). The transition is abrupt and it occurs at some η c > 0: this value decreases noticeably when increasing As a final aspect, we wish to better characterize the emergence of BCs in terms of the clustering order param-eter R B and the angle γ between the phase clusters, introduced in the Sub-section II B. As already shown, the BCs emerge for protocol (II A) simulations, therefore we limit to these ones. In particular, we have redone the phase diagram reported in Fig. 9 (e) by employing the clustering order parameter R B , as shown in Fig. 12 (a). This panel reveals that, across the negative semi-plane of λ 2 , we have indeed a very strong level of cluster-synchronicity associated to BC states and that, despite the presence of the SN bifurcation line (still clearly visible), the BCs remain well defined before and after traversing such a bifurcation line. The white arrows drawn in the diagram tell the position and sense of the cuts shown in panels (b) and (c), where we display the de-synchronization transitions obtained via quasi-adibatic simulations by varying λ 1 and λ 2 , respectively. As a matter of fact, we started with clusterized initial conditions with η = 1 and varied the coupling parameters adiabatically, in contrast to the simulations shown in panel (a) where the initial conditions are newly generated for each pixel.

In panel (b) we report R 1 , R B and γ versus λ 1 for the case λ 2 = -30. This simulation corresponds to the protocol (II A) previously examined. We observe that, while R 1 decreases for decreasing λ 1 , the parameter R B increases, thus indicating that the unique cluster present at λ 1 = 20, corresponding to γ = 0, begins to split in two for decreasing λ 1 , as indicated by the growing value of γ. This grows is abruptly interrupted at λ (π) 1 when the two clusters reach the phase opposition, corresponding to γ = π. Immediately after, they rearrange in a configuration with a smaller angle γ, while the new BC state is characterized by a higher R 1 and a smaller R B . By further decreasing λ 1 the angle tends to increase as R B , while R 1 decreases. The BC state destroys abruptly at some negative λ 1 value and the system becomes asynchronous. This latter transition occurs also when one maintains λ 1 constant and increases λ 2 , as suggested by the vertical white arrow in panel (a). This transition is reported in more detail in panel (c). As λ 2 increases, the clusters begin to de-synchronize slowly, thus R 2 decreases monotonically, until a critical value is reached and both clusters disappear abuptly. An interesting remark is that the oscillators that are dropping out from the clusters are not random, but they are the ones with the natural frequencies furthest from the center of the distribution. This can be verified looking both at the snapshots of the phases and the time-averaged velocities for two different λ 2 values, corresponding to the points circled in pink and blue in panel (c), along the knee of the transition. Specifically, to the λ 2 values highlighted in pink and blue in panel (c), correspond the dynamical behaviors shown in panel (d, e), where the color code is maintained: the snapshots of the phases (the average velocities) are reported in panel (d) ((e)). From these figures it is clear that only the oscillators with natural frequencies in proximity of ω i = 0 are synchronized and their number decrease by approaching the transition towards the incoherent state.

V. SUMMARY AND OUTLOOK

In this work, we have considered a higher-order version of the Kuramoto model combining the usual pair-wise first harmonic interaction term with a symmetric three body interaction for unimodal and bimodal Gaussian distributions of the natural frequencies, with a particular emphasis on the latter distribution. We have characterized, via a self-consistent approach, the different synchronization (de-synchronization) transitions occurring by increasing (decreasing) the parameter λ 1 controlling the pair-wise interactions, starting from de-synchronized (fully synchronized) initial conditions. To better understand the origin of the transitions at the mean-field level we have analyzed the dynamics of a reference oscillator as that of an overdamped particle moving in a potential landscape V (ϕ) (11). Depending on the sign of the parameter λ 2 controlling the higher order interactions, the potential has striking different characteristics: for positive λ 2 the potential shows the existence of two minima with different depth (controlled by λ 1 ), corresponding to two clusters in anti-phase; for negative λ 2 we have two symmetric minima at a distance γ in phase, which is controlled by λ 1 .

For unimodal Gaussian distribution and positive λ 2 the scenario of the possible transitions is reported in Fig. 3. In particular, we observe that the incoherent state looses stability via a super-critical Hopf bifurcation by giving rise to two asymmetric clusters at the same critical coupling found for the usual Kuramoto model [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF]. The 2 clusters are usually not equally populated, tipically a fraction σ of the oscillators will be in one minimum and a fraction 1 -σ in the other one. In this context the FS regime, corresponding to all the oscillators located in the main minimum of the potential (σ = 1), remains stable even for very negative λ 1 values and it looses stability abruptly exhibiting a discontinous transition to the incoherent state at λ (F S) 1 . In the mean-field context, such transition is expected to occur at λ (F S) 1 > -λ 2 , whenever the minimum of V (ϕ), corresponding to ϕ = 0, disappears. In a wide range of parameters we observe the coexistence between the incoherent regime and the coherent states characterized by different values of σ. The synchronization transition is clearly hysteretic in this case.

As shown in Fig. 4, for negative λ 2 and unimodal Gaussian distribution, we observe that the incoherent regime looses stability via a smooth transition giving rise to two equally populated phase clusters composed by oscillators with positive and negative natural frequencies, termed Uniform Clusters (UCs). On the other hand, starting from a synchonized initial condition at large λ 1 , a two cluster state arises with equally populated clusters, each including oscillators with either positive or negative natural frequencies: the so-called Bimodal Clusters (BCs). The BCs are stable for λ 1 > λ (π) 1 , and disappear at λ (π) 1 via the so-called π-transition : a SN bifurcation, where the two clusters reach a perfect phase opposition. This induces a vanishing of the order parameter as well as of the mean-field potential V (ϕ), thus giving rise to the abrupt transition. Immediately after, a two cluster state is reformed but it cannot be considered anymore a BC, since now the oscillators of the two previous clusters are mixed up. In particular, in each cluster we have mostly oscillators with ω i of a specific sign, but there are also a few oscillators with ω i of opposite signs (see Fig. 4 (c)). This new state becomes unstable at a definitely negative λ 1 value. In this case, depending on the level of synchronization η of the initial conditions, we observe cluster states, which present π-transitions at smaller and smaller λ (π) 1 > 0 values for decreasing η, followed by a two cluster state disappearing at λ

(P S) 1 that approaches λ (AS) 1
for η → 0 (as shown in Fig. 5). Indeed for η → 0 the π-transition disappears and the continuous transition obtained starting from the incoherent state is recovered, suggesting that no hysteretic behaviour is present in this context.

Whenever we introduce bimodal Gaussian distributions we observe, with respect to the previous scenario, that, on one side, for sufficiently positive or negative λ 2 , the incoherent state looses stability at λ (AS) 1 via discontinuous transitions. On the other side, at intermediate values of λ 2 , the transition from the AS to the coherent regime is mediated by the emergence of SWs, observable in a region that enlarges with increasing distance between the two peaks of the bimodal distributions, as displayed in Fig. 9. On the contrary, the scenario describing the loss of stability of the FS states initialized with η = 1 is essentially identical to the one observed for the unimodal distribution, with the only difference that, for sufficiently negative λ 2 , the clusterized state that is reformed, after the π-transition involving the BC state, has indeed a higher level of coherence, but it can be still identified as a BC (see Fig. 7 (c)).

In this paper we have put special emphasis on the negative (repulsive) coupling since this situation is still far from being fully explored. This was previously analysed in [28] for unimodal Gaussian distributions, where the authors observed that, when both couplings become repulsive, synchrony still persists in the system in the form of BC states. Here we have shown that BCs emerge in a much wider region of the negative half-plane for bimodal Gaussian distributions and that the π-transition can be better characterized by introducing a new order parameter R B (Eq. ( 7)) that takes into account the level of synchronization within each single cluster. In particular, R B growths smoothly with the angle γ between the two clusters, while the Kuramoto-Daido order parameter R 2 , usually employed to characterize a two cluster state, does not show a monotonic dependence on γ, as shown in Fig. 1. By employing this new indicator we observe that the π-transition is characterized by a desynchronization of the system (in terms of R 1 ), which is due to the fact that the 2 clusters have reached an angular distance γ = π, thus being in anti-phase. After the transition, the clusters find themselves at an angular distance γ smaller than before the transition. The level of synchronization within each cluster can remain unchanged, what changes abruptly is the angular distance γ between the two clusters, and this explains the jump in R 1 at the π-transition shown in Fig. 12.

Another interesting result of our work, is that the negative higher-order interactions render the SWs observable in a wider range of parameters with respect to the usual Kuramoto model. We must note that, in [31], it has been carried out an analysis similar to the one here reported. In particular in [31] the authors consider four body asymmetric interactions with Lorentzian bimodal distributions of the natural frequencies, thus being able to obtain analytically the mean-field evolution via the OA approach. Their main result is that the higher order interactions essentially enlarge the bi-stable regions found in [11]. Nonetheless, our higher order term is quite different from their one; moreover we have analysed repulsive interactions, while in [31] they limited themselves to attractive ones. In the future it would be worth examining how the presence of additional peaks in the natural frequency distribution would affect the emergence of the SWs, as well as the scenario here reported.

Furthermore, another general aspect that we have highlighted by performimg quasi-adiabatic simulations with a simplified version of the original model (Eq. (A.1)), is that, for attractive higher order interactions, the synchronization transitions are hysteretic with a large region of coexistence, while this is not the case for repulsive higher order ones, as clearly shown in Fig. 14.

Recently, the analysis of an enlarged Kuramoto model encompassing a first harmonic pair-wise interaction plus symmetric and anti-symmetric three body interactions has revelead the emergence of collective chaos [32]. In particular, in [32], the authors have obtained mean-field results for unimodal Gaussian frequency distributions by employing a truncated Hermite-Fourier decomposition of the oscillator density. It would be worth performing an analysis to interpret the observed transition to collective chaos in terms of our mean-field potential landscape V (ϕ) and to compare the results of the self-consistent approach with those reported in [32].

Finally, it would be interesting to combine the usual pair-wise first harmonic interaction term with three body interactions also for the Kuramoto model with inertia, thus extending the results obtained in [30,[START_REF] Olmi | Control of Self-Organizing Nonlinear Systems, Understanding Complex Systems[END_REF], where it has been investigated in detail the hysteretic transition to synchronization for a network of Kuramoto oscillators with inertia, in presence of the usual pair-wise first harmonic interaction term, for unimodal and bimodal Gaussian distributions of the natural frequencies.

ACKNOWLEDGMENTS

We acknowledge useful interactions with Nina La Miciotta, Arkady Pikovsky and Diego Pazó. AT received financial support by the Labex MME-DII (Grant No ANR-11-LBX-0023-01), by the ANR Project ERMUNDY (Grant No ANR-18-CE37-0014), and by CY Generations (Grant No ANR-21-EXES-0008), all part of the French programme "Investissements d'Avenir". APM gratefully acknowledges financial support by the Spanish Ministerio de Economía y Competitividad and European Regional Development Fund under contract PID2022-138322OB-I00 AEI/FEDER, UE, and by Xunta de Galicia (together with AC) under Research Grant No. 2021-PG036. All the programs acknowledged by APM are co-funded by FEDER (UE).

Appendix: Simplified Version of the Model

In this Appendix, we present a simplified version of the Kuramoto model with higher order interactions and show that it recovers all the transitions and dynamical behaviors that we have reported so far. In particular, we notice that one of the parameters entering Eq. ( 4) is actually redundant and the formulation of the model can be simplified by rescaling the time as t ′ = t • |λ 2 | (the frequencies as ω ′ i = ω i /|λ 2 |) and by introducing an unique coupling parameter ξ = λ 1 /|λ 2 |. The new version of the model reads as :

dθ i (t) dt ′ = ω ′ i + ξ N N j=1 sin (θ j -θ i ) + sign(λ 2 ) 1 2N 2 N j=1 N k=1
sin (θ j + θ k -2θ i ) , (A.1) which corresponds to two distinct ordinary differential equations depending on the sign of λ 2 , one with attractive higher order interactions (plus sign) and another with repulsive one (minus sign).

A detailed description of the phase diagram of the system according to the actual parameters is shown in Figure 13, where the top panels correspond to attractive three-body interactions and the lower ones to repulsive In each panel the dynamical states found with the twoparametric model are tagged, similarly to what done in Figure 9. For the acronyms reported in the panels we refer to the states described in the main text. For each value of the symmetric peaks centers, we re-scaled the standard deviation ∆ accordingly, so that the relation ω ′ 0 /∆ ′ = 4 is fulfilled. This scaling choice was mainly motivated in order to compare the results of this Appendix with the ones presented in the previous sections.

In Figure 14 we show two cuts of the two dimensional phase diagrams for ω ′ 0 = 0.1 and 0.5. Again, top panels account for positive, attractive triadic interactions and lower panels for negative, repulsive ones. In (a), for ω ′ 0 = 0.1, we observe a direct abrupt synchronization transition analogous to the one shown in Fig. 6 (a). The new feature to highlight here is that, if the backward simulations (obtained by decreasing ξ) are initialized with some clustered state with finite R 1 value, this state is maintained until ξ becomes sufficiently negative to destroy such a state. The critical negative value of ξ is larger and larger for higher level of synchronization of the initial state, since the depth of the minimum of the potential where the synchronized state resides depends on R 1 , as previously shown. These results indicate that there is a large region of multistability where clusters of different level of synchronization can coexist. For larger ω ′ 0 = 0.5 the synchronization transition is preceeded by the emergence of SW, as expected (see panel (b)). The SWs are characterized by a large value of the standard deviation of R 1 , since R 1 is now oscillating in time. Also in the present case we have a clear hysteretic region.

For repulsive triadic interaction, the SW regime occupies a quite large portion of the phase diagram for sufficiently distant peaks, i.e. sufficiently large ω ′ 0 . This is evident both in panels (c) and (d). In particular for ω ′ 0 = 0.5, the partially synchronized regimes are always characterized by standing waves and we are not able to observe a FS regime. For negative triadic interactions no hysteresis is observable for sufficiently large ω ′ 0 . At lower ω ′ 0 we can find a coexisting region for the clustered regimes, UCs and BCs, as in the standard model. This is shown in Fig. 15 for ω ′ 0 = 0.01, where the order parameters R 1 and R B are reported as a function of ξ for quasi-adiabatic simulations. As displayed in Fig. 15 (a), we observe a smooth transition in R 1 by increasing ξ, that leads to an UC regime, as it is evident from the value of the clustering order parameter R B in panel (b). In the backward simulations, the BC is clearly present and it de-synchronizes via the π-transition, before reforming. As shown in panel (b), the order parameter R B tends to one approaching the π-transition where the BC desynchronizes. After the π-transition, the cluster states are reformed and both R 1 , R B smoothly de-synchronize for decreasing ξ-values.

FIG. 2 .

 2 FIG. 2. Potential V (ϕ) (11) and frequencies Ω(ϕ) (12) versus the phase ϕ for positive values of λ2 = 15 (a-b) and negative one λ2 = -100 (c-d). Black (red) lines refer to V (ϕ) (Ω(ϕ)), while dashed green (blue) vertical lines denote the angles ±α (±β). The red solid (dashed) curves identify the stable (unstable) branches associated to phases locked oscillators. Parameters are set in both panels to R1 = 0.5 and Ω = 0 for the potential, while panel (a) (panel (b)) refers to λ1 = 3 (λ1 = -3) and panel (c) (panel (d)) to λ1 = 10 (λ1 = -10).

FIG. 3 .

 3 FIG. 3. Unimodal Gaussian natural frequencies distribution and λ2 > 0. (a) R1 versus λ1 : the blue and maroon circles refer to simulations performed by following protocols (I A) and (II A), respectively. The value reported for R1 represents the average in time over a simulation time Ts. The diamonds of different colours denote the self consistent evaluation of the order parameter R1, for different σ values, obtained by employing (19). In the inset it is shown σ versus λ1 as estimated by following protocol (I A). In panels (b-e) are reported snapshots of the natural frequencies ωi versus their phases θi for different values of λ1 for simulations obtained by following protocol (I A). Analogous snapshots are shown in (f) for simulations done by following protocol (II A). When necessary, the snapshots are arbitrarly shifted along the x-axis to improve readability. The displayed results refer to λ2 = 15, ∆ = 0.25 for gU (ω), ∆λ1 = 0.2, Tt = 10, Ts = 200, integration step dt = 0.001 for the adiabatic simulations, and a network size N = 10000.

FIG. 4 .

 4 FIG. 4. Natural frequencies distributed according to an unimodal Gaussian distribution and λ2 < 0. (a) R1 versus λ1 : the blue and maroon circles refer to simulations performed by following protocols (I A) and (II A), respectively. The value reported for R1 represents the average in time over a simulation time Ts. The cyan (magenta) line refers to the self consistent evaluation of the order parameter R1 by employing (21) estimated in the whole range of frequencies (restricted to positive frequencies). In panels (b) and (c) are reported snapshots of the natural frequencies ωi versus their phases θi for different values of λ1 for simulations obtained by following protocols (I A) and (II A), respectively. In the inset of panel (c) we report the instantaneous values of R1(t) (black curve) and RB(t) (red curve) versus time for λ1 = 3 obtained by following protocol (II A). The displayed results refer to λ2 = -100, ∆ = 0.25 for gU (ω), ∆λ1 = 0.2 Tt = 10, Ts = 200, integration step dt = 0.001 for the adiabatic simulations, and a network size N = 10000.

FIG. 5 .

 5 FIG. 5. The π transition. (a) R1 versus λ1 for different initial η values (i.e. different percentage of initially synchronized oscillators): the black dots, blue and magenta triangles and orange stars refer to simulations performed by following protocol (II A). The value reported for R1 represents the average in time over a simulation time Ts. The cyan line refers to the simulations performed by following protocol (I A) and already shown in Fig. 4. In panels (b)-(e) are reported snapshots of the natural frequencies ωi versus their phases θi for different initial values of η, for simulations obtained by following protocols (II A). In all panels black dots correspond to λ1 = 2.2, red dots to λ1 = 6. The displayed results refer to λ2 = -100, ∆ = 0.25 for gU (ω), ∆λ1 = 0.2 Tt = 10, Ts = 200, integration step dt = 0.001 for the adiabatic simulations, and a network size N = 10000.

FIG. 6 .

 6 FIG. 6. Natural frequencies distributed according to a bimodal Gaussian distribution and λ2 > 0. (a) R1 versus λ1 : the blue and maroon circles refer to simulations performed by following protocols (I A) and (II A), respectively. The value reported for R1 represents the average in time over a simulation time Ts. The diamonds of different colours denote the self consistent evaluation of the order parameter R1, for different σ values, obtained by employing (19). In the inset it is shown σ versus λ1 as estimated by following protocol (I A). In panels (b-e) are reported snapshots of the natural frequencies ωi versus their phases θi for different values of λ1 for simulations obtained by following protocol (I A). Analogous snapshots are shown in (f) for simulations done by following protocol (II A). The displayed results refer to λ2 = 15, ω0 = 1 and ∆ = 0.25 for gB(ω), ∆λ1 = 0.2, Tt = 10, Ts = 200, integration step dt = 0.001 for the adiabatic simulations, and a network size N = 10000.

FIG. 7 .

 7 FIG. 7. Natural frequencies distributed according to a bimodal Gaussian distribution and λ2 < 0. (a) R1 versus λ1 : the blue and maroon circles refer to simulations performed by following protocols (I A) and (II A), respectively. The value reported for R1 represents the average in time over a simulation time Ts. The cyan (magenta) line refers to the self consistent evaluation of the order parameter R1 by employing (21) estimated in the whole range of frequencies (restriced to positive frequencies). In panels (b) and (c) are reported snapshots of the natural frequencies ωi versus their phases θi for different values of λ1 for simulations obtained by following protocols (I A) and (II A), respectively. In the inset of panel (c) we report the instantaneous values of R1(t) (black curve) and RB(t) (red curve) versus time for λ1 = 8.8 obtained by following protocol (II A). The displayed results refer to λ2 = -100, ω0 = 1 and ∆ = 0.25 for gB(ω), with ∆λ1 = 0.2, Tt = 10, Ts = 200, integration step dt = 0.001 for the adiabatic simulations, and a network size N = 10000.

FIG. 8 .

 8 FIG.8. Dynamical regimes: asynchronous state (AS) obtained for (λ1,λ2)=(-10,-10); standing wave (SW) for (λ1,λ2)=(3,-5); bimodal cluster (BC) for (λ1,λ2)=(10,-1000) and uniform cluster (UC) coexisting with the BC for the same parameters. Row-wise, instantaneous phases θi (a-d) and time-averaged phase velocities ⟨ θi⟩ (e-h) versus the natural frequencies ωi of the oscillators, together with the corresponding time evolution of the order parameters R1 and RB (i-l), for the considered dynamical state. UCs (BCs) are obtained by employing the simulation protocol (II) (protocol (I)). For all cases the bimodal distribution gB(ω) with ω0 = 1.0 and ∆ = 0.25 is considered and the system size is N = 10000.

FIG. 9 .

 9 FIG. 9. Modulus of the Kuramoto order parameter R1 versus the coupling parameters λ1 and λ2 for three different cases : unimodal Gaussian distribution of the natural frequencies with ω0 = 0 (a, b) and bimodal Gaussian distribution with ω0 = 1 (c, d) and ω0 = 3 (e, f). In all cases the standard deviation of the Gaussians is fixed to ∆ = 0.25. The upper (lower) row refers to results obtained by following the protocol (I) (protocol (II)). The reported abbreviations refer to different dynamical regimes : asynchronous state (AS); fully synchronous case (FS); standing waves (SW); bimodal cluster (BC) and uniform cluster (UC). The white dashed lines mark the onst of the FS regime. The size of the network has been fixed to N = 5000 and the parameter R1 has been averaged over a time Ts = 200 once discarding a transient time Tt = 100.

λ 2 .

 2 Therefore this parameter has a de-stabilising effect on the asyncronous regime, favouring the emergence of the BCs, since the depth of the two symmetric minima in the potential increases with |λ 2 | (see Fig.2 (b)).

FIG. 11 .

 11 FIG. 11. Clustering order parameter RB for increasing values of the initial degree of phase-locking η and for different values of the triadic coupling λ2 ∈ [-1000 : -100]. In (a), λ1 = 10. The extreme points correspond to the panels (c) and (d) of Figure 8. In (b), where λ1 = -10, we pass from an AS to a BC regime at different critical points ηc as a function of the λ2 value. Parameters are ω0 = 1.0, ∆ = 0.25 and N = 5000.

FIG. 12 .

 12 FIG. 12. Numerical characterization of the π-transition and of the BC state. In (a), the order parameter RB is a function of the coupling strengths λ1 and λ2, only for the simulation protocol (II A). The white arrows give the sense and direction to the de-synchronization transitions shown in panels (b) and (c), obtained by varying quasi-adiabatically the coupling parameters. The value reported for R1, RB, γ in panels (b) and (c) represent the average in time over a simulation time Ts. The angle γ is rescaled by a factor π to plot all variables on the same scale. For (b) λ2 = -30, while for (c) λ1 = -20. The pink and blue circles highlighted in panel (c) refer to (λ1, λ2)=(-20, -450) and (λ1, λ2)=(-20, -250) respectively. Using the same colorcode (pink and blue), these states are analyzed in panels (d,e) in terms of their instantaneous phases and averaged velocities. The snapshot presented in (d) is related to the instantaneous phases against their natural frequencies. In (e) the time-averaged velocity of each oscillator is reported against its natural frequency. For all panels, the bimodal parameters were fixed at ω0 = 1 and ∆ = 0.25.

FIG. 13 .

 13 FIG.13. Kuramoto order parameter R1 (a-d) and clustering order parameter RB (e-h) as a function of the coupling parameter ξ and of the peak position of the bimodal Gaussian distribution ω ′ 0 in rescaled units. Upper row refers to attractive triadic interactions, while the lower row to repulsive ones. These simulations were performed quasi-adiabatically as explained in the text, by either increasing ξ, as in panels (a,c,e,g), or decreasing ξ, as in panels (b,d,f,h). For each value of ω ′ 0 , we fixed ∆ ′ = 0.25ω ′ 0 . The total number of oscillators is N = 5000. For every simulation Ts = 200 -1000 and Tt = 100 -500 with modification of ξ with a step ∆ξ = 0.02 in the quasi-adiabatic simulations.

  interactions. In particular in such a figure are shown the order parameters R 1 (panels a-d) and R B (e-h) for quasi-adiabatic simulations, where we start at η = 0 and ξ = -1.0 and we increase quasi-adiabatically ξ up to ξ = 1.0. Then ξ is decreased quasi-adiabatically, starting from the state at ξ = 1.0 obtained from the previous simulations, down to ξ = -1.0. The parameter is always varied in steps of ∆ξ = 0.02. The forward synchronizing transitions are shown in panels (a,c,e,g) while the backward de-synchronizing ones are shown in panels (b,d,f,h). Note that R 1 and R B are measured from the same simulations and thus they are complementary.

FIG. 15 .

 15 FIG. 15. The π-transition. Kuramoto order parameter R1 (a) and clustering order parameter RB (b) versus ξ for repulsive triadic interactions and ω ′ 0 = 0.01. The simulation is run adiabatically starting at η = 0 with a total number of N = 5000 oscillators.