Collective dynamics of heterogeneously and nonlinearly coupled phase oscillators

Coupled oscillators have been used to study synchronization in a wide range of social, biological, and physical systems, including pedestrian-induced bridge resonances, coordinated lighting up of firefly swarms, and enhanced output peak intensity in synchronizing laser arrays. Here we advance this subject by studying a variant of the Kuramoto model, where the coupling between the phase oscillators is heterogeneous and nonlinear. In particular, the quenched disorder in the coupling strength and the intrinsic frequencies are correlated, and the coupling itself depends on the amplitude of the mean field of the system. We show that the interplay of these factors leads to a fascinatingly rich collective dynamics, including explosive synchronization transitions, hybrid transitions with hysteresis absence, abrupt irreversible desynchronization transitions, and tiered phase transitions with or without a vanishing onset. We develop an analytical treatment that enables us to determine the observed equilibrium states of the system, as well as to explore their asymptotic stability at various levels. Our research thus provides theoretical foundations for a number of self-organized phenomena that may be responsible for the emergence of collective rhythms in complex systems.

Figure 1

Phase diagrams of the time-averaged order parameter $\langle R\rangle$ vs the global coupling strength $K$, as obtained for different values of the exponent $\beta$. The panels show results for $\beta =1$ (left column), $\beta <1$ (middle column), and $\beta >1$ (right column). For (a1)–(a3), we have used $\beta =1$ and $g\left(\omega \right)=\frac{\mathrm{\Delta }}{2\pi }[\frac{1}{{(\omega -{\omega }_0)}^2+{\mathrm{\Delta }}^2}+\frac{1}{{(\omega +{\omega }_0)}^2+{\mathrm{\Delta }}^2}]$, $K_f=4/\sqrt{1+{({\omega }_0/\mathrm{\Delta })}^2}$, $K_c=2$, and $R_c=\sqrt{2}/2$. In (a1), ${\omega }_0=\mathrm{\Delta }=0.1$ and $K_f=2\sqrt{2}>K_c$, which leads to a first-order-like phase transition. In (a2), ${\omega }_0=\sqrt{3}/10$, $\mathrm{\Delta }=0.1$, and $K_f=2=K_c$, which, on the other hand, leads to a hybrid phase transition. In (a3), ${\omega }_0=\sqrt{15}/10$, $\mathrm{\Delta }=0.1$, and $K_f=1<K_c$, for which the system undergoes a tiered phase transition, where the oscillatory state arises in the region $K\in (1,2)$. For (b1)–(b3), we have used $\beta <1$, $g\left(\omega \right)=1/2$ with $\omega \in (-1,1)$, $K_f=0$, and $K_c={(\beta +1)}^{\frac{\beta +1}{2}}{\beta }^{-\frac{\beta }{2}}$. Under these conditions, the system undergoes a tiered phase transition with a vanishing onset from the oscillatory state to the two-cluster state. Parameter values are as follows: (b1) $\beta =0.2$, $K_c=1.31$, $R_c=0.41$; (b2) $\beta =0.5$, $K_c=1.61$, $R_c=0.58$; (b3) $\beta =0.8$, $K_c=1.86$, $R_c=0.67$. For (c1)–(c3), we have used $\beta >1,g\left(\omega \right)=1/\sqrt{2\pi }{e}^{-{\omega }^2/2}$, $K_f=+\infty$, and $K_c={(\beta +1)}^{\frac{\beta +1}{2}}{\beta }^{-\frac{\beta }{2}}$, leading to an irreversible abrupt desynchronization transition. Parameter values are as follows: (c1) $\beta =1.5$, $K_c=2.32$, $R_c=0.77$; (c2) $\beta =2$, $K_c=2.60$, $R_c=0.82$; (c3) $\beta =3$, $K_c=3.08$, $R_c=0.87$. For the simulations, we have used $N=100000$, $t=950$ with time step $\mathrm{\Delta }t=0.01$, and $T=50$. The red circles and the blue crosses represent the forward and backward continuations, respectively. The black lines depict theoretical predictions obtained with the mean-field approach, whereby solid and dashed lines depict stable and unstable branches, respectively.

Figure 2

Microscopic properties of the different oscillatory states. The top columns show the instantaneous phases in dependence on the natural frequencies, while the bottom columns show the effective frequencies $\langle {\omega }_i\rangle =1/T{\int }_t^{t+T}{\dot{\theta }}_i\left(t\right)dt$ in dependence on the natural frequencies. The insets show the time series of the order parameter. The following parameters have been used: in (a1) and (a2), $\beta =0.5$, $g\left(\omega \right)=1/2$ with $\omega \in (-1,1)$, $K=1.5$; in (b1) and (b2), $\beta =0.8$, $g\left(\omega \right)=1/2$ with $\omega \in (-1,1)$, $K=1.7$; in (c1) and (c2), $\beta =1.0$, $g\left(\omega \right)=\frac{\mathrm{\Delta }}{2\pi }[\frac{1}{{(\omega -{\omega }_0)}^2+{\mathrm{\Delta }}^2}+\frac{1}{{(\omega +{\omega }_0)}^2+{\mathrm{\Delta }}^2}]$ (i.e., $g\left(\omega \right)$ is a bimodal Lorentzian distribution with $\pm {\omega }_0$ being the locations of the two peaks of the distribution, and $\mathrm{\Delta }$ is the half width of the distribution) with ${\omega }_0=0.3,\mathrm{\Delta }=0.1$, $K=1.6$; and in (d1) and (d2), $\beta =1.0$, the bimodal Lorentzian distribution $g\left(\omega \right)=\frac{\mathrm{\Delta }}{2\pi }[\frac{1}{{(\omega -{\omega }_0)}^2+{\mathrm{\Delta }}^2}+\frac{1}{{(\omega +{\omega }_0)}^2+{\mathrm{\Delta }}^2}]$, with ${\omega }_0=0.3,\mathrm{\Delta }=0.1$, $K=1.8$. For the simulations, we have used the same values for $N$, $t$, $\mathrm{\Delta }t$, and $T$ as in Fig. 1.