{'fraction_non_alphanumeric': 0.03379022945222674, 'fraction_numerical': 0.010870110960486018, 'mean_word_length': 4.523433166943982, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of fully coupled rotators with unimodal and bimodal frequency distribution

We analyze the synchronization transition of a globally coupled network of N phase oscillators with inertia (rotators) whose natural frequencies are unimodally or bimodally distributed. In the unimodal case, the system exhibits a discontinuous hysteretic transition from an incoherent to a partially synchronized (PS) state. For sufficiently large inertia, the system reveals the coexistence of a PS state and of a standing wave (SW) solution. In the bimodal case, the hysteretic synchronization transition involves several states. Namely, the system becomes coherent passing through traveling waves (TWs), SWs and finally arriving to a PS regime. The transition to the PS state from the SW occurs always at the same coupling, independently of the system size, while its value increases linearly with the inertia. On the other hand the critical coupling required to observe TWs and SWs increases with N suggesting that in the thermodynamic limit the transition from incoherence to PS will occur without any intermediate states. Finally a linear stability analysis reveals that the system is hysteretic not only at the level of macroscopic indicators, but also microscopically as verified by measuring the maximal Lyapunov exponent.


Introduction
The renowned Kuramoto model [14] for phase oscillators was generalized in 1997 by Tanaka, Lichtenberg and Oishi (TLO) [24,25] by including an additional inertial term. The TLO model revealed, at variance with the usual Kuramoto model, first order synchronization transitions even for unimodal distributions of the natural frequencies. TLO have been inspired in their extension by a work of Ermentrout published in 1991 [7]; in this paper Ermentrout has introduced a pulse coupled phase oscillator model with inertia to mimic the perfect synchrony achieved by a specific type of fireflies, the Pteroptix Malaccae (but also by certain species of crickets and humans). The peculiarity of these fireflies is that they are able to synchronize their flashing activity to some forcing frequency (even quite distinct from their own intrinsic flashing frequency) with an almost zero phase lag. This happens because they adapt their period of oscillation to that of the driving oscillator. After his introduction, the Kuramoto model with inertia has been employed to describe synchronization phenomena in crowd synchrony on Londons Millennium bridge [23], as well as in Huygens pendulum clocks [4]. Furthermore, phase oscillators with inertia (rotators) have recently found application in the study the self-synchronization in power and smart grids [6,8,16,20,21], as well as in the analysis of disordered arrays of underdamped Josephson junctions [26]. Cluster explosive synchronization has been reported for an adaptive network of Kuramoto oscillators with inertia, where the natural frequency of each oscillator is assumed to be proportional to the degree of the corresponding node [12]. Rotators arranged in two symmetrically coupled populations have recently revealed the emergence of intermittent chaotic chimeras [17], imperfect chimera states have been found in a ring with nonlocal coupling [11], and transient waves have been observed in regular lattices [13].
There is a wide literature devoted to coupled rotators with an unimodal frequency distribution, however only a really limited number of studies have been devoted to this model with a bimodal distribution, despite the subject being extremely relevant for the modelization of the power grids [8,18,20]. To our knowledge the synchronization transition in populations of globally coupled rotators with bimodal distribution has been previously analyzed only by Acebrón et al. in [1]. More specifically, the authors considered a model with white noise and a distribution composed by two δ-functions localized at ±Ω 0 . As suggested in [15], the presence of noise blurs the δ-functions in bell-shaped functions analogous to Gaussian distributions. Therefore one expects a similar phenomenology to the one observable for deterministic systems with bimodal Gaussian distributions of the frequencies. A multiscale analysis of the model, in the limit of sufficiently large Ω 0 , reveals the emergence from the incoherent state of stable standing wave solutions (SWs) and of unstable traveling wave solutions (TWs) via supercritical bifurcations, while partially synchronized stationary states (PSs) bifurcates subcritically from incoherence. However, the authors affirm that in the considered limit the bifurcation diagram coincides with that of the usual Kuramoto model without inertia [2].
In this article we analyze the synchronization transitions observable for unimodal and bimodal frequency distributions for a population of globally coupled rotators in a fully deterministic system. In particular, we will analyze the influence of inertia and of finite size effects on the synchronization transitions. Moreover, we will study the macroscopic and microscopic characteristics of the different regimes emerging during adiabatic increase and decrease of the coupling among the rotators. In particular, Section 1.2 will be devoted to the introduction of the model, of the indicators used to characterize the synchronization transition, and of the different protocols employed to perform adiabatic simulations. The Lyapunov linear stability analysis is introduced in SubSection 1.2.1. The results for unimodal distributions are reported in Sect. 1.3, with a particular emphasis to the TLO mean field theory and its extension to any generic state observable within the hysteretic region (SubSect. 1.3.2). The emergence of clusters of locked and whirling oscillators is described in details in SubSect. 1.3.3. The dynamics of the network for bimodal distributions is analyzed in Section 1.4. In particular, SubSect. 1.4.1 is devoted to two non overlapping distributions and SubSect. 1.4.2 to largely overlapping Gaussian distributions. Sect. 1.5 report the result of linear stability analysis for the considered distributions. Finally, in Sect. 1.6 the reported results are briefly summarized and discussed.

Model and Indicators
By following Refs. [24,25], we study the following version of the Kuramoto model with inertia for N fully coupled rotators : where θ i and Ω i are, respectively, the instantaneous phase and the natural frequency of the i-th oscillator, K is the coupling. In the following we will consider random natural frequencies Ω i Gaussian distributed according to: an unimodal distribution g(Ω) = 1 √ 2π e − Ω 2 2 with zero average and an unitary standard deviation or a bimodal symmetric distribution g(Ω) = , which is the overlap of two Gaussians with unitary standard deviation and with the peaks located at a distance 2Ω 0 .
To measure the level of coherence between the oscillators, we employ the complex order parameter [27] where r(t) ∈ [0 : 1] is the modulus and φ(t) the phase of the macroscopic indicator. An asynchronous state, in a finite network, is characterized by r ≃ 1 √ N , while for r ≡ 1 the oscillators are fully synchronized and intermediate r-values correspond to partial synchronization.
Another relevant indicator for the state of the rotator population is the number of locked oscillators N L , characterized by a vanishingly small average phase velocityω i ≡d θi dt , and the maximal locking frequency Ω M , which corresponds to the maximal natural frequency |Ω i | of the locked oscillators.
In general we will perform sequences of simulations by varying adiabatically the coupling parameter K with two different protocols. Namely, for the first protocol (I) the series of simulations is initialized for the decoupled system by considering random initial conditions for {θ i } and {ω i }. Afterwards the coupling is increased in steps ∆K until a maximal coupling K M is reached.
For each value of K, apart the very first one, the simulations is initialized by employing the last configuration of the previous simulation in the sequence. For the second protocol (II), starting from the final coupling K M achieved by employing the protocol (I), the coupling is reduced in steps ∆K until K = 0 is recovered. At each step the system is simulated for a transient time T R followed by a period T W during which the average value of the order parameter r and of the velocities {ω i }, as well as Ω M , are estimated.

Lyapunov Analysis
The stability of Eq. (1.1) can be analyzed by following the evolution of infinitesimal perturbations T = (δθ 1 , . . . , δθ N , δθ 1 , . . . , δθ N ) in the tangent space, whose dynamics is ruled by the linearization of Eq. (1.1) as follows: We will limit to estimate the maximal Lyapunov exponent λ M , by employing the method developed by Benettin et al. [3]. This amounts to follow the dynamical evolution of the orbit and of the tangent vector T for a time lapse T W by performing Gram-Schmidt ortho-normalization at fixed time intervals ∆t, after discarding an initial transient evolution T R .
Furthermore, the values of the components of the maximal Lyapunov vector T can give important information about the oscillators that are more actively contributing to the chaotic dynamics. It is useful to introduce the following squared amplitude component of the normalized vector for each rotator [9,17] The time averageξ i of this quantity gives a measure of the contribution of each oscillator to the chaotic dynamics.

Hysteretic Synchronization Transitions
In Fig. 1.1 the results for a sequence of simulations obtained by following protocol (I) and (II) are reported for a not too small inertia (namely, m = 2) and unimodal frequency distribution. For the first protocol the system remains incoherent up to a critical value K = K c 1 ≃ 2, wherer jumps to a finite value and then increases with K reachingr ≃ 1 for sufficiently large coupling. Starting from the last state and by reducing K one notices thatr assumes larger values than during protocol (I) and the system becomes incoherent at a smaller coupling, namely K c 2 < K c 1 . This is a clear indication of the hysteretic nature of the synchronization transition in this case.
For the chosen values of the inertia, we observe the creation of an unique cluster of N L locked oscillators withω i ≃ 0, for larger m the things will be more complex, as we will discuss in the following. The maximal locking frequency Ω M becomes finite for K > K c 1 and increases with K. The frequency Ω M attains a maximal value whenr ≃ 1, no more oscillators can be recruited in the large locked cluster. Once reached this value, even if K is reduced following the protocol (II), Ω M remains constant for a wide K interval. Then Ω M shows a rapid decrease towards zero by approaching K c 2 . This behavior will be explained in the following two sub-sections.

Mean field theory
In order to derive a mean field description of the dynamics of each single rotator, we can rewrite Eq. (1.1) by employing the order parameter definition (1.2) as follows this corresponds to the evolution equation for a damped driven pendulum. Eq. (1.5) admits, for sufficiently small forcing frequency Ω i , two fixed points: a stable node and a saddle. At larger frequencies Ω i > Ω P ≃ 4 π Kr m the saddle, via a homoclinic bifurcation, gives rise to a limit cycle. The stable limit cycle and a stable fixed point coexist until a saddle node bifurcation, taking place at Ω i = Ω D = Kr, leads to the disappearance of the fixed points and for Ω i > Ω D only the limit cycle persists. This scenario is correct for sufficiently large inertia; at small m one has a direct transition from a stable node to a periodic oscillating orbit at Ω i = Ω D = Kr [22]. Therefore for sufficiently large m there is a coexistence regime where, depending on the initial conditions, the single oscillator can rotate or stay quiet. The fixed point (limit cycle) solution corresponds to locked (drifting) rotators. The TLO theory [24,25] has explained the origin of the first order hysteretic transitions by considering two opposite initial states for the network: (I) the completely incoherent phase (r = 0) and (II) the completely synchronized one (r ≡ 1). In case (I) the oscillators are all initially drifting with finite velocities ω i ; by increasing K the oscillators with smaller natural frequencies |Ω i | < Ω P begin to lock (ω i = 0), while the other continue to drift. This is confirmed by the data reported in Fig. 1.1, where it is clear that the locking frequency Ω M is well approximated by Ω P . The process continues until all the oscillators are finally locked, leading to r = 1 and to a plateau in Ω M .
In the second case, initially all the oscillators are already locked, with an associated order parameter r ≡ 1. Therefore, the oscillators can start to drift only when the stable fixed point solution will disappear, leaving the system only with the limit cycle solution. This happens, by decreasing K, whenever |Ω i | ≥ Ω D = Kr. This is numerically verified, indeed, as shown in Fig. 1.1, where it is clear that the maximal locked frequency Ω M remains constant until, by decreasing K, it encounters the curve Ω D and then Ω M follows this latter curve down towards the asynchronous state. The case (II) corresponds to the situation observable for the usual Kuramoto model, where there is no bistability [14].
In both the considered cases there is a group of desynchronized oscillators and one of locked oscillators separated by a frequency, Ω P (Ω D ) in case (I) (case (II)). At variance with the usual Kuramoto model, both these groups contribute to the total level of synchronization, namely where r L (r D ) is the contribution of the locked (drifting) population. The contribution of the locked population is simply given by where θ P = sin −1 (Ω P /Kr) and θ D = sin −1 (Ω D /Kr) ≡ π/2. The contribution r D of the drifting rotators is negative and it has been estimated analytically by TLO by performing a perturbative expansion to the fourth order in 1/(mK) and 1/(mΩ). The obtained expression, valid for sufficiently large inertia, reads as By considering an initially desynchronized (fully synchronized) system and by increasing (decreasing) K one can get a theoretical approximation for the level of synchronization in the system by employing the mean-field expression (1.7), (1.8) and (1.6) for case I (II). In this way, two curves are obtained in the phase plane (K, r), namely r I (K) and r II (K). For a certain coupling K the system can attain all the possible levels of synchronization between r I (K) and r II (K).
Let us notice that the expression for r L and r D reported in Eqs. (1.7) and (1.8) are the same for case (I) and (II), only the integration extrema change in the two cases. These are defined by the frequency which discriminates locked from drifting neuron, that in case (I) is Ω P and in case (II) Ω D . It should be noticed that the value of these frequencies is a function of the order parameter r and of the coupling constant K, therefore one should solve implicit integrals to obtain r.
However, one could also fix the discriminating frequency to some arbitrary value Ω 0 and solve self-consistently the equations Eqs.
Therefore the part of the plane delimited by the curve r II (K), will be filled with the curves r 0 (K) obtained for different Ω 0 values (as shown in Fig. 1.2(a)). These solutions represent clusters of N L oscillators for which the maximal locking frequency and N L do not vary upon changing the coupling strength. In particular, for K > K c 2 these states can be observed in numerical simulations in the portion of the phase space delimited by the two curves r I (K) and r II (K) (see Fig. 1.2(b)).

Clusters of Locked and Whirling Oscillators
By observing the results reported in Fig. 1.2(a) for m = 6, it is evident that the numerical data obtained by following the procedure (II) are quite well reproduced from the mean field approximation r II (solid grey curve). This is not the case for the theoretical estimation r I (dashed grey curve), which does not reproduce the step-wise structure revealed for the data corresponding to protocol (I). This step-wise structure emerges only for sufficiently large inertia (as it is clear from Fig. 1.1 (b), where it is absent for m = 2); this is due to the break down of the independence of the whirling oscillators: namely, to the formation of clusters of drifting oscillators moving coherently at the same non zero velocity [24]. Oscillators join in small groups to the locked stationary cluster and not individually as it happens for smaller inertia; this is clearly revealed by the behavior of N L versus the coupling K as reported in Fig. 1.2(b).
Furthermore, once formed, these stationary locked clusters are particularly robust, as it can be appreciated by considering as initial condition a partially synchronized state obtained following protocol (I) for a certain coupling K S > K 1 . This state is characterized by a cluster of N L locked; if now we reduce the coupling K, the number of locked oscillators remain constant until we do not reach the descending curve obtained with protocol (II), see the black empty triangles in Fig. 1.2. On the other handr decreases slightly with K, this behavior is well reproduced by the mean field solutions of Eq. (1.9), namely r 0 (K, Ω 0 ) with Ω 0 = Ω P (K s , r I (K S )) = 4 π Ksr I m , these are shown in Fig. 1.2 (a) as black dot-dashed lines. As soon as, by decreasing K, the frequency Ω 0 becomes equal or smaller than Ω D , the order parameter has a rapid drop towards zero following the upper limit curve r II . These results indicate that hysteretic loops of any size are possible within the region delimited by the two curvesr obtained by following protocol I and II respectively, as shown in [18].
For sufficiently small m, the synchronization occurs starting from the incoherent state via the formation of an unique cluster of locked oscillators, and the size of this cluster increases with K, as evident from Fig. 1.3 (a) for m = 2. At the same time the value of r also increases with K and its evolution is characterized only by finite size fluctuations vanishing in the thermodynamic limit (see the inset of Fig. 1.3 (a)). As already mentioned, the situation is different for sufficiently large inertia, now the partially synchronized phase is characterized by the coexistence of the main cluster of locked oscillators with ω i ≃ 0, but also by the emergence of clusters composed by drifting oscillators with common finite velocities, see the data forω i reported in Fig. 1.3 (b) for m = 6. In particular, the clusters of whirling oscillators emerge always in couple and they are characterized by the same average velocity but opposite sign. These states are indicated as standing waves (SWs), therefore we have a SW coexisting with a partially synchronized stationary state (PS) (as shown in Fig. 1.3 (b). The effect of these extra clusters on the collective dynamics is to induce oscillations in the temporal evolution of the order parameter, as one can see from the inset of Fig. 1.3 (b). In presence of drifting clusters characterized by the same average velocity (in absolute value), as for m = 6 and K = 5 in Fig. 1.3 (b), r exhibits almost regular oscillations and the period of these oscillations corresponds to the one associated to the oscillators in the drifting cluster.

Bimodal Distribution
In this Section we consider a bimodal distribution, we will initially focus on two almost non overlapping Gaussians, namely we consider Ω 0 = 2, while Sub-Section 1.4.2 is devoted to overlapping Gaussians, examined for Ω 0 = 0.2.

Non Overlapping Gaussians
For Ω 0 = 2 and sufficiently small inertia (m = 1 and 2), we observe a very rich synchronization transition, as shown in Fig. 1.4 (a). In particular, by following protocol (I) we observe that the system leaves the incoherent state abruptly by exhibiting a jump to a finiter value at K T W ; above such value in the network emerges a single cluster of oscillators, drifting together with a finite velocity ≃ Ω 0 , this corresponds to Traveling Wave (TW) solution . By further increasing K a second finite jump of the order parameter at K SW denotes the passage to a Standing Wave (SW) solution, corresponding to two clusters of drifting oscillators with symmetric opposite velocities ≃ ±Ω 0 . A final jump at K P S leads the system to a Partially Synchronized (PS) phase, characterized by an unique cluster of locked rotators with zero average velocity. By increasing the coupling the PS state smoothly approaches the fully synchronized regime . Starting from this final state the return sequence of simulations, following protocol (II), displays a simpler phenomenology. The network stays in the PS regime, characterized by an order parameter larger than that measured during protocol (I) simulations, until K DS < K T W . For smaller coupling, the system leaves the PS state; however, depending on the realization of the natural frequencies and on the initial conditions, it can ends up or in a TW (most of the cases) or in a SW, or it can even reach directly the incoherent state (as shown in Fig. 1.4 (a) and Fig. 1.6 (a)). This first analysis clearly shows hysteretic effects and coexistence of macroscopic states with different level of synchronization for a wide range of couplings.
Let us now try to examine the observed transitions in terms of the maximal locking frequency, Ω M . This frequency is now defined in a different way with respect to the unimodal distribution, in the present case Ω M represents the maximal absolute value of the natural frequencies of the oscillators belonging to the main clusters present in the system, therefore in this estimation are considered both stationary and drifting clusters. As shown in Fig. 1.4 (b), Ω M increases with K for protocol (I) simulations. In particular Ω M shows a finite jump in correspondence of K = K T W , and then its evolution is reasonably well approximated by the curve Ω P + Ω 0 , where Ω P = 4 π Kr m . By approaching K P S the maximal frequency displays a constant plateau which extends beyond K P S , this indicates that the two symmetric drifting clusters merge at K = K P S giving rise to an unique locked cluster with zero average velocity, however no other oscillators join this cluster up to a larger coupling. Whenever this happens , Ω M starts again to increase, but this time it follows the curve Ω P = 4 π Kr m . Finally, forr ≃ 1 the maximal locking frequency attains a maximal value. Moreover, by reducing the coupling, following now protocol (II), Ω M remains stacked to such a value for a wide K interval. The fully synchronized cluster is difficult to break down due to the inertia effects. Finally, Ω M reveals a rapid decrease towards zero whenever it encounters the curve Ω D = Kr, initially it follows this curve, however as soon as the system desynchronizes towards a TW the decrease of Ω M is better described by the curve Ω D + Ω 0 . The observed behavior can be explained by the fact that for K < K P S (K < K (DS) ) for protocol (I) (protocol (II)), the network behaves as two independent sub-networks each characterized by an unimodal frequency distribution, one centered at Ω 0 and the other one at −Ω 0 . The extension of the analysis reported in Section 1.3.2 for an unimodal distribution not centered around zero simply amounts to shift the limiting curves Ω P and Ω D by Ω 0 . However, for sufficiently large coupling constant, once the system exhibits only one large cluster with zero velocity, the network behaves as a single entity and Ω M closely follows Ω P or Ω D as for a single unimodal distribution centered in zero.
Let us now describe the TW and SW states in more details with the help of the examples reported in Fig. 1.5 for inertia m = 2 and N = 2000. The TW is an asymmetrical cluster of whirling oscillators with a finite velocitȳ ω i ≃ Ω 0 , in particular in Fig. 1.5 (a)  in a range around Ω 0 , namely 0.67 ≤ Ω i ≤ 3.34. The effect of this cluster on the collective dynamics is to increase the average value of the order parameter without inducing any clear oscillating behavior in r(t). However, oscillators with positive natural frequencies are much more synchronized with respect to the ones with negative frequencies, as can be inferred by observing the order parameter r p (r n ) estimated only on the the sub-population of oscillators with positive (negative) natural frequencies and reported in Fig. 1.5 (b). We believe that the emergence of the TW state is related to the finite sampling of the distribution of the natural frequencies, which due to finite size effects can be non perfectly symmetric. The asymmetric cluster will emerge around +Ω 0 (−Ω 0 ) depending on the positive (negative) sign of the average natural frequency. The SWs are observable at larger coupling constants, this is characterized by two symmetrical clusters with opposite average velocities ≃ ±Ω 0 , as shown in Fig. 1.5 (c). The presence of these two clusters now induces clear periodic oscillations in the order parameter r(t), as observable in Fig. 1.5 (d). The period of the oscillations is related to |Ω 0 |, i.e. the average frequency of the clustered oscillators. However, at variance with the results reported in Fig. 1.3 (b) for the unimodal distribution the two symmetric clusters do not coexist with a cluster of locked oscillators with zero average velocity. By examining separately r p and r n reported in Fig. 1.3 (d), we notice that each sub-population is much more synchronized than the global one, in fact r p and r n have a higher average value than r with superimposed irregular oscillations. The data reported so far refer to a single system size, however finite size effects are quite relevant for this model, as shown in [18] for unimodal distributions. In Fig. 1.6 (a) we report the synchronization transition for several system sizes, namely 1, 000 ≤ N ≤ 50, 000 for a small inertia value (m = 1). We observe that K T W and K SW increase with the size N ; in particular, the incoherent state is observable on a wider coupling interval by increasing N (similarly to what reported in [18] for unimodal distributions). Finite size fluctuations induce transitions from the incoherent branch to the TW branch and from this to the SW branch. The fact that we do not observe transition back to the original states indicates that the energy barriers are higher from these sides. A quite astonishing result is the fact that the transition value K P S and K DS seem completely independent from N . The combination of these results seem to suggest that in the thermodynamic limit the incoherent state will loose stability at K P S and therefore the two branches corresponding to TW and SW will be no more visited, at least by following protocol (I).
By observing all the data reported in Fig. 1.6 (a) for various N and for protocol (I) and (II), it seems that there are clear indications that the two branches of solutions, corresponding to TW and SW, emerge via a supercritical bifurcation at the same coupling, namely K ≃ 3.8, while the transition to PS is clearly subcritical. These results confirm the analysis reported in [1] for a system with noise (in particular, see Fig. 17 in that paper). However, Acebrón et al. affirm that the SW is stable, while the TW is unstable. From our results, both branches seem to become inaccessible (in absence of noise) from the incoherent state, while at least a part of these branches appear to be reachable from the PS state by decreasing K below K DS following protocol (II). Another important difference with respect to the results reported in [1] is that the PS regime is clearly hysteretic revealing two coexisting branches of PS states visited by following protocol (I) or (II). As shown in Fig. 1.6 (b), for larger inertia (namely, m = 6), the transition from the incoherent state following protocol (I) occurs via the emergence of many small clusters leading finally to a SW state. In this case the critical value at which the incoherent state looses stability seems to saturate to a constant value already for N ≥ 2, 000. The value of K P S is also in this case insensible to the system size. For large inertia values, the TWs seem no more observable.
As a further aspect, we will report the numerical results of the dependence on the inertia of the critical coupling constant K P S , while the value of K DS ≃ 4.9 is independent not only by N , but also by the inertia. As shown in Fig.  1.7, K P S increases linearly with the inertia and this scaling is already valid for not too large inertia values. The linear scaling with the inertia is analogous to the scaling recently found within a theoretical mean-field analysis for the coupling K MF 1 , which delimits the range of linear stability of the asynchronous state [1,10]. In particular, the authors in [18] have shown for a Gaussian unimodal distribution of width σ that K MF 1 ≃ 2σ(0.64 + mσ), which shows a linear dependence on the inertia and a quadratic dependence on the variance of the frequency distribution.
In the final part of this sub-section we perform an analysis analogous to that reported in Sub-Sect. 1.3.3, in particular starting from states with a finite level of synchronization obtained by following protocol (I) we decrease the coupling and observe how these states evolve. In Fig. 1.8, we report the results of these simulations (shown as empty triangles) for two different inertia values, namely m = 1 and m = 10. Starting from PS states we observe that the cluster survives until the descending curve obtained with protocol (II) is encountered, analogously to the results reported in Fig. 1.3 for the unimodal distribution. Therefore any part of the hysteretic portion of the (K, r)-plane delimited by the PS curves obtained via protocol (I) or (II) is accessible . However, if one starts for m = 1 from a TW or a SW state, one observes only two curves (corresponding to the TW and SW branches previously discussed) which seem to end up at the same critical coupling which is smaller than K DS . Therefore it seems that there are no evidences of hysteresis for this small inertia for SW and TW solutions (as shown in Fig. 1.8 (a)). For large inertia values m = 10, since now, apart the SW solutions, there are solutions with many small clusters, the situations is much more complex. By starting from different values of K < K P S and by decreasing K, these curves seem all to end up at the same critical coupling smaller K P S , see Fig. 1.8 (b). These results suggest that for large inertia values is possible to observe a continuum of possible states even starting from states characterized by (many) drifting clusters and that these states coexist in a wide range of coupling.

Overlapping Gaussians
In this sub-section, we analyze a bimodal distribution, where the two Gaussians are largely overlapping, since Ω 0 = 0.2. In this case we expect to observe a phenomenology of the synchronization transition quite similar to the one seen for the unimodal case. In Fig. 1.9(a) is reported the average order parameterr versus the coupling constant K estimated by following protocol (I) and (II) for various inertia values and for a fixed system size, namely N = 10, 000. We observe that all the curves obtained for protocol (II) almost overlap irrespectively of the used inertia, while the protocol (I) curves reveal a strong dependence on m. In particular, the hysteretic region widens with m. For small inertia values, namely m = 1 and 2, there is a sudden transition from the asynchronous state to a PS state at K P S and neither traveling waves nor standing waves are observable: a single cluster at zero velocity emerges in correspondence of K P S and the order parameter never shows oscillating behavior in time.
For m = 6, it is possible to observe a scenario similar to the one reported in Fig. 1.3 (b), where not only a cluster at zero velocity is present, but also two symmetrical clusters at finite velocities emerge. In particular, following protocol (I) for K > 2.4 a small cluster of locked oscillators emerges; at larger coupling, namely K ≥ 3, two symmetrical clusters of whirling oscillators emerge and coexist with the zero velocity cluster. Finally, at K = 10.8 the PS regime arises, corresponding to a single large cluster of locked oscillators. Furthermore, in the range 3 ≤ K < 10.8 the order parameter reveals irregular oscillations. A more detailed analysis is needed to understand the origin of these oscillations as done in the next Section.
Finally, we examine the influence of the system size on the studied transitions: for m = 1 the results for the protocol (I) [protocol (II)] simulations are reported in Fig. 1.9(b) for sizes ranging from N = 1, 000 to N = 10, 000. It is immediately evident that the transition from synchronized state to the asynchronous state, following protocol (II), does not depend on N: for all considered sizes the transition happens in correspondence of K ≃ 2, analogously to what reported for unimodal distributions [18]. Starting from the incoherent regime and following protocol (I) the system reveals a jump to a finiter value for critical couplings increasing with N , quite similar once more to the results reported for unimodal distributions. We can conclude this sub-section by affirming that the phenomenology seen for bimodal, but largely overlapping, distributions should not differ much from the one observed for unimodal distributions.

Linear Stability Analysis
To better characterize the synchronization transitions and the stability of the observed states it is worth estimating the maximal Lyapunov exponent λ M following protocol (I) and (II) for an unimodal and a bimodal distributions. This quite time consuming analysis has been performed for a single inertia value m = 6 and a single system size N = 1, 000, the scaling of λ M with N will be discussed in the following for specific coupling constant values.
In general, we observe that once the system fully synchronizes, λ M vanishes; therefore for most of the simulations associated to protocol (II) corresponding to fully synchronized cluster down to the desynchronization transition, λ M is zero. This is not the case for protocol (I) simulations which reveal a positive λ M as soon asr is non zero. Thus indicating that not only the dynamics characterized in terms of the macroscopic order parameterr is hysteretic, but also at the level of the microscopic dynamics, investigated via λ M , the system has a clear hysteretic behavior.
The behavior of λ M with K exhibits chaotic dynamical states with windows of regularity for both unimodal and bimodal distribution with Ω = 2, as shown in Figs. 1.10 (a) and (b). As a general aspect, we observe the maximal level of chaoticity immediately after the transition from the incoherent state to partially coherence, where small clusters of synchronized oscillators and drifting oscillators coexist. The increase ofr is accompanied by a trend of λ M to decrease and finally to vanish forr → 1 An important aspect to understand is if this dynamics is weakly chaotic or not, in particular this amounts to verify if, in the thermodynamic limit, λ M will vanish or will remain finite. In order to test for this aspect, we have considered a configuration obtained by following protocol (I) for a specific coupling and analyzed λ M versus the system size for 200 ≤ N ≤ 32, 000. The results for unimodal distributions, as well as for bimodal ones with Ω 0 = 2 and Ω 0 = 0.2 are shown in Fig. 1.10 (e). It is clear for all the considered cases that the system remains chaotic for diverging system sizes.
As a final aspect we would like to understand which oscillators contribute more to the chaotic dynamics of the system; this can be understood by measuring the average squared amplitude of the components of the maximal Lyapunov vectorξ i (see the definition reported in Eq. 1.4). In particular, we consider the three cases analyzed in Fig. 1.10 (e) for N = 1, 000. The corresponding results are shown in Fig. 1.11. From panel (a) and (b) of the figure it is clear that for the unimodal distribution, as well as for the largely overlapping bimodal distributions, the chaotic activity is associated almost exclusively to the rotators which are outside the large clusters of locked oscil- lators withω i ≃ 0. Thus confirming recent results reported for two coupled populations of rotators with identical natural frequencies [17]. However, the situation for the bimodal distribution with Ω 0 = 2 is different; in particular, as shown in Fig. 1.11 (c), the network for this large value of the inertia and the considered coupling does not exhibit a cluster of locked oscillators withω i ≃ 0, but only drifting clusters. In this case the rotators outside and inside the clusters seem to contribute equally to the maximal Lyapunov vector, with the possible exclusion of a group of rotators with Ω i ≃ Ω 0 .

Conclusions
We have studied the synchronization transition for a globally coupled Kuramoto model with inertia for different frequency distributions. For the unimodal frequency distribution we have shown that clusters of locked oscillators of any size coexist within the hysteretic region. This region is delimited by two curves in the plane individuated by the coupling and the average value of the order parameter. Each curve corresponds to the synchronization (desynchronization) profile obtained starting from the fully desynchronized (synchronized) state. For sufficiently large inertia values, clusters composed by drifting oscillators with opposite velocities (standing wave state) emerge in addition to the locked oscillators clusters. The presence of clusters of whirling rotators induces oscillatory behavior in the order parameter.
For bimodal frequency distribution the scenario can become more complex since it is possible to play with an extra parameter: the distance between the peaks of the distributions. For simplicity we have analyzed only two cases: largely overlapping distributions (Ω 0 = 0.2), and almost not overlapping distributions (Ω 0 = 2). The phenomenology observed for Ω 0 = 0.2 resembles strongly that found for the unimodal distribution. The analysis of the non overlapping case reveals new interesting features. In particular, the transition from incoherence to coherence occurs via several states: namely, traveling waves, standing waves and finally partial synchronization. This scenario resembles that reported for the usual Kuramoto model for a bimodal distribution [5,15,19]. However, in our case the transition is always largely hysteretic, and for non overlapping distributions, traveling waves are clearly observable at variance, not only with the results for the Kuramoto model [5,15,19], but also with the theoretical phase diagram reported in [1] for oscillators with inertia. A peculiar aspect is that in the thermodynamic limit we expect a direct discontinuous jump from the incoherent to the coherent phase, without passing through any intermediate state. The critical coupling K P S required to pass from incoherence to partial synchronization is independent of the system size and grows linearly with inertia, while the partially synchronized state looses its stability at a smaller coupling K DS < K P S which is the same for any inertia value and system size.
Finally, by performing a linear stability analysis we have been able to show that the hysteretic behavior is not limited to macroscopic observables, as the level of synchronization, but it is revealed also by microscopic indicators as the maximal Lyapunov exponent. In particular, we expect that in a large interval of coupling values chaotic and non chaotic states will coexist.