Predicting the Oscillatory Regimes of Global Synchrony Induced by Secondary Clusters

TL;DR

Inertia-bearing oscillator networks governed by the second-order Kuramoto model exhibit nonmonotonic, oscillatory global synchrony in the forward process due to secondary whirling clusters. The authors develop a self-consistent mean-field framework that resolves these secondary clusters, solving a reduced locked/drift problem and deriving quantitative criteria for their onset and termination. They identify an onset crossover mass $\tilde{m}^* \simeq 3.865$ and provide a self-consistent equation for the secondary-cluster order parameter $r_{(+)}$, enabling prediction of oscillatory regimes in the $(K,m)$ plane. This work offers practical guidance for stabilizing inertial synchronization systems, such as power grids, by predicting when oscillations occur and how to suppress them through parameter choices.

Abstract

Synchronization systems with effective inertia, such as power grid networks and coupled electromechanical oscillators, are commonly modeled by the second-order Kuramoto model. In the forward process, numerical simulations exhibit a staircase-like growth of global synchrony, reflecting temporal oscillations induced by secondary synchronized clusters of whirling oscillators. While this behavior has been observed previously, its governing conditions have not been quantitatively determined in terms of analytical criteria. Here, we develop a self-consistent theoretical framework that explicitly characterizes the secondary synchronized clusters. This analysis identifies an onset crossover mass $\tilde{m}^* \simeq 3.865$ for the emergence of secondary clusters and yields quantitative criteria for predicting both the crossover mass and the termination coupling strength at which they vanish. As a result, we determine the oscillatory regimes of coupling strengths over which global synchrony shows temporal oscillations, providing practical guidance for controlling and avoiding undesirable oscillatory behavior in inertial synchronization systems, such as power grids.

Figures (5)

• Figure 1: Formation of secondary clusters at $K=6$ and $m=6 > \tilde{m}^*$ with the marginal crossover mass $\tilde{m}^*\simeq 3.865$ (see Fig. \ref{['fig:critical_m']}) for $N=5\,000$. (a) Evolution of the order parameter $R$. (b) Natural frequency $\Omega_i$ versus the temporal-averaged phase velocity $v$ for all oscillators, computed over a time window $\Delta t = 200$ at $t = 10\,000$. The black solid and blue dotted vertical lines indicate the cluster boundaries $\Omega_{(\pm),0}$ and $\Omega_{(\pm)}$, corresponding to Eq. (\ref{['eq:omega_boundaries']}). The red dashed vertical line denotes the theoretical boundary separating the main cluster $\mathbb{N}_0$ from the secondary clusters $\mathbb{N}_{(\pm)}$ (see Sec. \ref{['subsec:solution']}). (c) Temporal-averaged phase velocity $v$ for the oscillators belonging to each cluster. For a representative time point, the mean and standard deviation of $v$ within each cluster are shown, and the minimum and maximum values in each cluster are highlighted by colored solid lines. (d) Cluster synchrony $\rho_C$ for $\mathbb{N}_0$ and $\mathbb{N}_{(\pm)}$. The triangle, circle, and inverted triangle denote $\mathbb{N}_0$, $\mathbb{N}_{(+)}$, and $\mathbb{N}_{(-)}$, respectively.
• Figure 2: Interaction of oscillators for $m=6$, $K=6$, and $N=5\,000$ with a temporal averaging window $\Delta t = 200$ at three representative times: (a, b) the initial state, (c, d) the stage where $\rho_0 \approx 1$ but $\rho_{(\pm)} < 1$, and (e, f) the stage where both $\rho_0 \approx 1$ and $\rho_{(\pm)} \approx 1$. The left panels (a–c) show $\theta_i$ and $\dot{\theta}_i$ in a rotating frame, where each point represents an oscillator with angular coordinate $\theta_i$ and radial coordinate $\dot{\theta}_i$. The gray dashed circle indicates $\dot{\theta}_i = 0$; larger (smaller) radial distance corresponds to more positive (negative) phase velocity. Colors denote intrinsic frequencies $\Omega_i$. Symbols and colored borders mark the pre-defined clusters from Fig. \ref{['fig:r_v_evolve']}(b): green circles for $\mathbb{N}_0$, pink triangles for $\mathbb{N}_{(+)}$, and cyan inverted triangles for $\mathbb{N}_{(-)}$. The magnitude of the global synchrony $R$ is also shown with mean phase $\phi = 0$. The right panels (d–f) display the interaction matrix $P_{ij}$, sorted by intrinsic frequency, with colored boxes indicating the cluster boundaries.
• Figure 3: Comparison between the theoretical analysis and numerical simulation. We use $N=5\,000$ oscillators on the $K\text{--}m$ plane ($0 \le K, m \le 10$) with $\mathrm{d}K=\mathrm{d}m=0.1$, evaluated at $t = 10\,100$ with a temporal-averaging window $\Delta t = 50$. (a) Scatter plot comparing the global order parameter $R$ obtained from simulation and theory [Eq. (\ref{['eq:z_decompose']})] for the parameter regime where secondary clusters exist. The color indicates the mass, and the solid line denotes $y=x$. (b) Phase diagram of the main-cluster synchrony $r_0$ for $\mathbb{N}_0$ on the parameter space, where color indicates the synchrony level. The dashed and solid curves mark the theoretical critical boundaries obtained from the self-consistent equation [Eq. (\ref{['eq:z_lock_drift']})]: the minimum $K_{\rm min}$ required to sustain nonzero $R$ (dashed) and the minimum $K_c$ at which $R$ first becomes nonzero (solid) [Fig. \ref{['fig:forward']}(b)]. Inset: density plot comparing simulation and theory, with $y=x$ (solid line). (c) Density plot of the secondary synchrony $r_{(+)}$ from simulation and theory [Eq. (\ref{['eq:rplus_selfcon']})], with $y=x$ marked by the solid line. Inset: simulation density plot of $r_{(+)}$ and $r_{(-)}$. (d) Phase diagram of the secondary-cluster synchrony $r^{\mathrm{(lock)}}_{(+)}$ for $\mathbb{N}_{(+)}$ on the parameter space. The dashed and solid curves represent the predicted onset boundaries of secondary synchronization given by $r_{(+)} = |R^{\mathrm{(drift)}}/2|$ and $r_{(+)} \to 0^{+}$, respectively. (e) Density plot comparing $r^{\mathrm{(lock)}}_{(+)}$ between simulation and the theoretical expression [Eq. (\ref{['eq:rplus_lock']})], with the solid line showing $y=x$. (f) Temporal fluctuation $\sigma_t$ of the global synchrony $R$ from simulation. Inset: density plot of $\sigma_t$ versus $r_{(+)}$, with the theoretical line $y = \sqrt{2 - \left({\pi}/{8}\right)^2 - \left({\pi}/{4}\right)^4}\, x$ from Eq. (\ref{['eq:sigma_rp']}).
• Figure 4: The crossover mass $m^{*}$ at which the secondary clusters can emerge or vanish. (a) The quantity $r_{(+)} + {R^{\mathrm{(drift)}}}/{2}$ as a function of $K$ for various $m$ from $m = 3.8$ to $m = 6$. Secondary clusters can exist when $r_{(+)} + {R^{\mathrm{(drift)}}}/{2} \ge 0$. The marginal value for the onset of the secondary clusters is $\tilde{m}^{*} \simeq 3.865$. (b) The secondary synchrony $r_{(+)}$ as a function of $m$ at $K = 7$ with $t = 10\,100$, $\Delta t = 50$, and $\mathrm{d}m = 0.1$. The arrows indicate the crossover points: $m^{*} = 4.3(1)$ from simulation, $m^{*} = 4.174(1)$ from theory based on $r_{(+)} \to 0^{+}$, and $m^{*} = 4.419(1)$ from theory based on $r_{(+)} = |R^{\mathrm{(drift)}}/2|$ (our conjecture). The shaded region represents the error bar computed over 100 realizations.
• Figure 5: The analytic and simulation results of the forward process (increase in $K$, $\mathrm{d}K=0.01$) at $m=6$ and $N=5\,000$, evaluated at $t=2\,500$ with $\Delta t=50$. (a) The global synchrony $R$ versus the coupling strength $K$. The critical points $K_{\mathrm{min}}$ and $K_c$ are theoretically evaluated as the smallest $K$ for which $r_0>0$ and the value of $K$ at which $r_0 \to 0^{+}$, respectively. (b) The main-cluster synchrony $r_0$ versus $K$. (c) The secondary synchrony $r_{(+)}$ versus $K$. (d) The intrinsic frequency $\Omega_i$ versus $K$ for all oscillators. White denotes drifting oscillators in $\mathbb{N}_{\Omega \gtrless \Omega_{(\pm)}}$, gray denotes the main cluster $\mathbb{N}_0$, and blue and red denote the secondary clusters $\mathbb{N}_{(\pm)}$. The two dashed vertical lines indicate the theoretical boundaries $r_{(+)} = \left|R^{\mathrm{(drift)}}/2\right|$ and $r_{(+)} \to 0^{+}$, and the vertically shaded regions mark the plateau region of $R$. Inset: histogram of the terminal point $K_t$ of the oscillatory regime over 100 realizations.