Hysteretic transitions in the Kuramoto model with inertia

We report finite-size numerical investigations and mean-field analysis of a Kuramoto model with inertia for fully coupled and diluted systems. In particular, we examine, for a gaussian distribution of the frequencies, the transition from incoherence to coherence for increasingly large system size and inertia. For sufficiently large inertia the transition is hysteretic, and within the hysteretic region clusters of locked oscillators of various sizes and different levels of synchronization coexist. A modification of the mean-field theory developed by Tanaka, Lichtenberg, and Oishi [Physica D 100, 279 (1997)] allows us to derive the synchronization curve associated to each of these clusters. We have also investigated numerically the limits of existence of the coherent and of the incoherent solutions. The minimal coupling required to observe the coherent state is largely independent of the system size, and it saturates to a constant value already for moderately large inertia values. The incoherent state is observable up to a critical coupling whose value saturates for large inertia and for finite system sizes, while in the thermodinamic limit this critical value diverges proportionally to the mass. By increasing the inertia the transition becomes more complex, and the synchronization occurs via the emergence of clusters of whirling oscillators. The presence of these groups of coherently drifting oscillators induces oscillations in the order parameter. We have shown that the transition remains hysteretic even for randomly diluted networks up to a level of connectivity corresponding to a few links per oscillator. Finally, an application to the Italian high-voltage power grid is reported, which reveals the emergence of quasiperiodic oscillations in the order parameter due to the simultaneous presence of many competing whirling clusters.

Figure 1

Average order parameter $\overline{r}$ (black circles) and maximal locking frequency ${\Omega }_M$ (blue triangles) as a function of the coupling $K$ for two series of simulations performed following protocol (I) (filled symbols) and protocol (II) (empty symbols). The data refer to mass: $m=2$ (a) and $m=6$ (b). For $m=2$ ($m=6$) we set $\Delta K=0.2$ ($\Delta K=0.5$) and $K_M=10$ ($K_M=20$), in both cases $N=500$, $T_R=5000$, and $T_W=200$. The (magenta) diamonds indicate ${\Omega }_P^{\left(I\right)}=\frac{4}{\pi }\sqrt{\frac{K\overline{r}}{m}}$ for protocol (I), the (green) squares ${\Omega }_D^{\left(II\right)}=K\overline{r}$ for protocol (II), and the (black) dashed vertical line the coupling constant $K_{FC}^G$, whose expression is reported in Eq. (14).

Figure 2

Minimal ${\omega }_S$ giving rise to a state characterized by a finite level of synchronization (i.e., $\overline{r}>0$) as a function of the coupling constant $K$. The inset reports the minimal number $N_S$ of oscillators which should be initially locked in order to lead to the emergence of a coherent state, as a function of $K$. The vertical (green) dot-dashed line refers to the estimated $K_1^c$ and the (blue) dashed line indicates the estimated $K_2^c$. The data refer to simulations performed with protocol (S) for $N=16000$, $m=6$, with $T_W=2000$ and $T_R=20000$.

Figure 3

Average order parameter $\overline{r}$ vs the coupling constant $K$ for $m=2$ (a) and $m=6$ (b). Mean-field estimates: the dashed (solid) red curves refer to $r^I=r_L^I+r_D^I$ ($r^{II}=r_L^{II}+r_D^{II}$) as obtained by employing Eqs. (5) and (6) following protocol (I) (protocol (II)); the (green) dot-dashed curves are the solutions $r^0(K,{\Omega }_0)$ of Eq. (7) for different ${\Omega }_0$ values. The employed values from bottom to top are ${\Omega }_0=1.21$ and 1.71 (a) and ${\Omega }_0=0.79$, 1.09, 1.31, and 1.79 (b). Numerical simulations: (blue) filled circles have been obtained by following protocol (I) and then (II) starting from $K=0$ until $K_M=10$ ($K_M=20$) for mass $m=2$ ($m=6$) with steps $\Delta K=0.2$ ($\Delta K=0.5$); (orange) filled triangles refer to simulations performed by starting from a final configuration obtained during protocol (I) and by decreasing the coupling from such initial configurations. The insets display $N_L$ vs $K$ for the numerical simulations reported in the main figures; the value of $K_F^G$ [Eq. (14)] is also reported in the two cases. The numerical data refer to $N=500$, $T_R=5000$, and $T_W=200$.

Figure 4

Maximal locking frequency ${\Omega }_M$ vs the coupling constant $K$. The initial state is denoted by the filled circle at $K_I=5$. The solid (red) curve indicates the frequency ${\Omega }_P^{\left(I\right)}$ and the dashed (green) curve the frequency ${\Omega }_D^{\left(II\right)}$. The numerical data refer to $N=500$, $T_R=5000$, $T_W=200$, $m=2$, and $\Delta K=0.1–0.05$.

Figure 5

Average order parameter $\overline{r}$ vs the coupling constant $K$ for various system sizes $N$: (a) $m=2$ and (b) $m=6$. The (red) solid and dashed curves are the theoretical estimates already reported in Fig. 3. The numerical data have been obtained by following protocol (I) and then protocol (II) from $K=0$ up to $K_M=10$ ($K_M=20$) for mass $m=2$ ($m=6$) with $\Delta K=0.2$ ($\Delta K=0.5$). The vertical dashed (maroon) lines refer to $K_1^{MF}$ reported in (8). Data have been obtained by averaging the order parameter over a time window $T_W=200$, after discarding a transient time $T_R\simeq 1000–80000$ depending on the system size; the larger $T_R$ have been employed for the larger $N$.

Figure 6

The upper symbols refer to the critical couplings $K_1^c$ (orange filled diamonds and red empty circles), while the lower ones to $K_2^c$ (cyan filled diamonds and blue empty circles) vs the system size $N$: (a) $m=0.8$, (b) $m=1$, (c) $m=2$, and (d) $m=6$. The filled symbols refer to estimates performed with protocol (S), while empty symbols, in panels (c) and (d), have been obtained with protocol (I) [protocol (II)] for $K_1^c$ ($K_2^c$). The dashed (black) lines in panels (a) and (b) are the mean-field values $K_1^{MF}$. This quantity is not reported in panels (c) and (d) for clarity reasons, due to its large value, namely, $K_1^{MF}=5.31$ for $m=2$ and $K_1^{MF}=13.27$ for $m=6$. For all panels the data have been derived by averaging in time over a window $T_W=2000$ and over 8 (5) different initial conditions for the protocol (S) [protocol (I) and (II)]. For each simulation an initial a transient time $T_R\simeq 20000$ ($T_R\simeq 1000–80000$) has been discarded for protocol (S) [protocol (I) and (II)].

Figure 7

The differences $K_1^{MF}-K_1^c\left(N\right)$ (filled orange diamonds) are reported vs the system size $N$ for $m=1$. The dashed (black) line is a power-law fit to the data: the difference vanishes, $5.27N^{-0.22}$ for $m=1$ (b). The data for $K_1^c$ are the same reported in panel (b) in Fig. 6.

Figure 8

Critical couplings $K_1^c$ (a) and $K_2^c$ (b) vs the mass $m$ for different system sizes $N$, namely, $N=1000$ (red diamond), 2000 (blue circles), 4000 (green triangles), 8000 (magenta squares), and 16 000 (black asterisks). The dashed (orange) line in (a) is the mean-field estimate $K_1^{MF}$, while the dashed (magenta) line in (b) is the value $K_m^{II}$ obtained by the TLO theory. The inset in panel (a) reports the critical rescaled couplings $\xi =(K_1^{MF}-K_1^c)/K_1^{MF}$ as a function of $m/N^{1/5}$. The estimates have been obtained with protocol (S), by averaging in time over a window $T_W=2000–4000$ and over eight different initial conditions. For each simulation an initial transient time $T_R\simeq 20000$ has been discarded.

Figure 9

(a) Average phase velocity $\langle {\dot{\theta }}_i\rangle$ of the oscillators vs their natural frequencies ${\Omega }_i$: (magenta) triangles refer to $K=2.5$, (green) diamond to $K=5$, (red) squares to $K=10$, and (black) circles to $K=15$. For each simulation an initial transient $T_R\simeq 5500$ has been discarded and the estimates have been obtained with protocol (I), by averaging in time over a window $T_W=5000$. (b) Order parameter $r\left(t\right)$ vs time for $m=6$ and $N=500$ and different coupling constants $K$: the (blue) solid curve corresponds to $K=1$, the (magenta) dot-dashed line to $K=2.5$, the (green) dashed line to $K=5$, the (red) dashed line to $K=10$, and the (black) solid line to $K=15$. The data have been obtained by employing protocol (I), and for each simulation an initial transient time $T_R\simeq 1500$ has been discarded and data are averaged over a time $T_W=5000$.

Figure 10

(a) Average phase velocity $\langle {\dot{\theta }}_i\rangle$ of the oscillators vs the corresponding natural frequency ${\Omega }_i$. The vertical dashed lines denote the three oscillators, $O_1$, $O_2$, and $O_3$, whose dynamical evolution is shown in panel (b). (b) The black curve represents the order parameter $r$ vs time; the other curves refer to the time evolution of the phase velocities $\dot{\theta }\left(t\right)$ of the three oscillators $O_1$ (red dot-dashed curve), $O_2$ (magenta solid line), and $O_3$ (dashed blue curve). For each simulation an initial transient time $T_R\simeq 1000$ has been discarded; the averages reported in (a) have been obtained over a time window $T_W=20000$. In both panels $K=5$, $m=6$, and $N=500$.

Figure 11

(a) Minima and maxima of the order parameter $r$ as a function of the coupling constant $K$. The (blue) dashed line refer to the theoretical estimate $r^I$, as obtained by employing Eqs. (5) and (6). (b) The number of oscillators in the drifting clusters $N_{DC}$ (filled black diamond) is reported vs the coupling $K$ together with the amplitude of the oscillations of the order parameter $r_{\max }-r_{\min }$ (empty red circles) rescaled by a factor 240. For each simulation an initial transient time $T_R\simeq 1500$ has been discarded. The estimates have been obtained with protocol (I), by averaging in time over a window $T_W=5000$, $m=6$, $N=500$.

Figure 12

(a) Average order parameter $\overline{r}$ vs the coupling constant $K$ for diluted networks for various $N_c$: 5 (filled black circles), 10 (red squares), 15 (green diamond), 25 (blue up triangles), 50 (orange left triangles), 125 (turquoise down triangles), 250 (right magenta triangles), 500 (violet crosses), 1000 (empty maroon circles), 2000 (black asterisks). (b) Critical constants $K_1^c$ and $K_2^c$ estimated for various values of the in-degree. The numerical data refer to $N=2000$; the upper inset refer to $N=1000$, the lower one to $N=500$. For all simulations $m=2$, $T_R=10000$, and $T_W=2000$; each series of simulations have been obtained by following protocol (I) and then (II) starting from $K=0$ until $K_M=20$ with steps $\Delta K=0.25$. The reported data have been obtained by averaging over 10–20 different series of simulations, each corresponding to a different realization of the random network and of the distribution of the frequencies $\left\{{\Omega }_i\right\}$. The error bars in panel (b) correspond to $\Delta K/2$.

Figure 13

Width of the hysteretic loop $W_h$, measured in correspondence of an order parameter value $\overline{r}=0.9$, as a function of the percentage of connected links $N_c/N$. The (green) circles refer to $N=500$, the (red) squares to $N=1000$, and the (black) diamonds to $N=2000$. The dashed lines refer to logarithmic fitting to the data in the range $0.01<N_c/N\le 1$. In the inset is graphically explained how $W_h$ has been estimated, starting from one of the curves reported in Fig. 12. The data refer to the same parameters and simulation protocols as in Fig. 12.

Figure 14

Average order parameter $\overline{r}$ vs the coupling constant $K$ for a diluted network with 70% of cut links. Mean-field estimates: the dashed (solid) red curves refer to $r^I=r_L^I+r_D^I$ ($r^{II}=r_L^{II}+r_D^{II}$) as obtained by employing Eqs. (5) and (6) following protocol I [protocol (II)]; the (green) dot-dashed curves are the solutions $r^0(K,{\Omega }_0)$ of Eq. (7) for different ${\Omega }_0$ values. The employed values from bottom to top are ${\Omega }_0=2.05$, 1.69, and 1.10. Numerical simulations: (blue) filled circles have been obtained by following protocol (I) and then (II) starting from $K=0$ until $K_M=20$ with steps $\Delta K=0.5$; (orange) filled triangles refer to simulations performed by starting from a final configuration obtained during protocol (I) and by decreasing the coupling from such initial configurations. The insets display $N_L$ vs $K$ for the numerical simulations reported in the main figures. The numerical data refer to $m=2$, $N=500$, $N_c=150$, $T_R=5000$, and $T_W=200$.

Figure 15

Average order parameter $\overline{r}$ vs the parameter $K$ for the Italian high-voltage power grid network. The (blue) circles data have been obtained by following protocol (I) from $K=1$ up to $K_M=40$ with $\Delta K=1$. The other symbols refer to simulations performed following protocol (II) starting from different initial coupling $K_I$ down to $K=1$, namely, (orange) triangles $K_I=10$, (red) squares $K_I=9$, and (green) diamonds $K_I=7$. The dashed vertical (magenta) line indicates the value $K=9$. The reported data have been obtained by averaging the order parameter over a time window $T_W=5000$, after discarding an initial transient time $T_R\simeq 60000$. The numerical data refer to $\alpha =1/6$, $N=127$, $\langle N_c\rangle =2.865$.

Figure 16

Average phase velocity of each oscillator $\langle {\dot{\theta }}_i\rangle$ vs the oscillator index for different values of the coupling $K$. The oscillators have been reordered so that the first 93 are consumers and the last 34 sources. The data have been obtained by employing protocol (I), starting from zero coupling $K=0$ and with $\Delta K=1$. For each simulation an initial transient time $T_R\simeq 5000$ has been discarded and the average is taken over a window $T_W=5000$. The numerical data refer to the same parameters as in Fig. 15.

Figure 17

Average phase velocity of each oscillator $\langle {\dot{\theta }}_i\rangle$ vs the corresponding oscillator index, ordered following the geographical distribution from northern Italy to Sicily. The panels refer to different couplings. The colored clusters indicate Italian regions which remain connected for all the considered simulations: red symbols refer to Piedmont and Liguria, green symbols to Veneto and Friuli Venetia Giulia, blue symbols to Campania and Apulia, and magenta symbols to Sicily. The data have been obtained by employing protocol (II) starting from $K_I=12$ with $\Delta K=1$ down to $K=1$. For each simulation an initial transient time $T_R\simeq 50000$ has been discarded and the averages performed over a window $T_W=5000$. The numerical data refer to the same parameters as in Fig. 15.

Figure 18

Order parameter $r\left(t\right)$ vs time for the Italian high-voltage power grid network for different values of the parameter $K$. (a) The dotted (magenta) curve corresponds to $K=4$, the dashed (black) line to $K=7$, the solid (orange) line to $K=9$, the dot-dashed (cyan) line to $K=10$. The data have been obtained by employing protocol (I), and for each simulation an initial transient time $T_R\simeq 5000$ has been discarded. (b) The solid (magenta) curve corresponds to $K=1$, the dashed (black) line to $K=4$, the solid (orange) thick line to $K=6$, the dot-dashed (cyan) line to $K=7$. The data have been obtained by employing protocol (II), and for each simulation an initial transient time $T_R\simeq 50000$ has been discarded. The simulations refer to the same parameters employed in Fig. 15.

Figure 19

Italian high-voltage power grid network with Gaussian distributed powers $P_i$ with zero average and unitary standard deviation. (a) Order parameter $r\left(t\right)$ vs time for different values of the parameter $K$. The data have been obtained by employing protocol (I) starting from $K=0.5$ to $K=20$ with $\Delta K=0.5$: the dot-dashed (red) line corresponds to $K=2$, the dashed (green) line to $K=3$, and the solid (blue) curve to $K=4$. (b) Average phase velocity for each oscillator vs the corresponding oscillator index, ordered following the geographical distribution from northern Italy to Sicily for coupling $K=2$. The clusters are colored as in Fig. 17. For each simulation an initial transient of $T_R=100000$ has been discarded, and the reported data have been mediated over a time window $T_W=5000$ and refer to $\alpha =1/6$.