Complete synchronization in networks of Sakaguchi-Kuramoto oscillators with bi-harmonic coupling

Abstract

In heterogeneous networks of coupled oscillators, phase frustration typically prevents the emergence of complete synchronization in the Sakaguchi-Kuramoto (SK) model. In this study, we propose an analytical framework to overcome this barrier and induce complete synchronization in oscillators governed by phase-frustrated bi-harmonic coupling. We derive a general set of natural frequencies correlated with the network's degree heterogeneity, along with the parameters involved in the bi-harmonic coupling function that lead to complete synchronization $(r=1)$ in the presence of the harmonic coupling terms $(K_1, K_2 \neq 0)$. On top of that, we found hysteresis in the synchronization transition in the case of scale-free networks, indicating a first-order (discontinuous) phase transition, whereas Erdos--Renyi networks exhibit a second-order (continuous) synchronization transition. Furthermore, we use mean-field approximation to determine the critical coupling strength for the synchronization transition in the absence of first-harmonic coupling $(K_1=0)$. Here, the obtained optimal natural frequencies scale linearly with the node degree, and the critical coupling strength for the onset of synchronization is derived analytically from the self-consistent equations. In this specific regime, we observe distinct dynamical disparities: the second harmonic drives an explosive first order (or second order) transition for the second order parameter $(r_2)$, while the first order parameter $(r_1)$ remains suppressed in the forward direction but emerges during the backward transition. These findings remain robust with higher-order harmonic coupling schemes, as well as across a diverse range of synthetic and empirical networks, including scale--free, Erdos--Renyi, Zachary Karate Club and C.elegans neural network, demonstrating their general applicability.

Figures (6)

• Figure 1: Synchronization dynamics in a scale-free network ($N = 500$, $\gamma = 2.8$) under phase frustration $\alpha = \beta = 0.5$. Purple diamond and green square markers denote the synchronization levels of the order parameters ($r_1$ and $r_2$) obtained from normal and homogeneous natural frequency distributions, respectively. The blue solid and red dashed curves correspond to forward and backward adiabatic continuation of the coupling strength, computed using the proposed optimal frequency assignment. (a)--(b) Variation of order parameters with the second-harmonic coupling $K_2$ at fixed $K_1 = 0.3$. (c)--(d) Variation of order parameters with the first-harmonic coupling $K_1$ at fixed $K_2 = 0.1$. In both parameter sweeps, the system undergoes a discontinuous (first-order) synchronization transition characterized by pronounced hysteresis. Complete synchronization is achieved at the prescribed target points (cyan markers), namely $K_2^{(p)} = 0.4$ in panels (a)–(b) and $K_1^{(p)} = 0.5$ in panels (c)–(d).
• Figure 2: Order parameters $r_1$ and $r_2$ vs. coupling strength $K_2$. Numerical and analytical results for a scale-free network ($N=500$, $\gamma=2.8$, $\langle q \rangle \approx 9$, $\beta=0.5$). Blue and red lines denote simulation results via forward and backward integration, respectively. (a) Behavior of $r_1$ showing only a backward phase transition. (b) Bifurcation diagram of $r_2$ exhibiting a first-order transition and a distinct hysteresis area. Cyan dots at $K_2^{(p)}=4$ represent complete synchronization $(r_1,_2=1)$, provided $K_1=0$. Solid and dashed magenta curves indicate the stable and unstable analytical branches derived from Eqs. \ref{['eq:real']} and \ref{['eq:imaginary']}. The analytical threshold (green dashed line) $K^c_2 \approx 2.1846 < K^{(p)}_2=4$ from Eq. \ref{['Critical_value']} shows agreement with the numerical simulation.
• Figure 3: The effect of frequency deviation on system synchronization. The graph compares the simulated error $\eta$ (red symbols) against the theoretical scaling limit $\eta \sim \sigma^2$ (blue line) under varying noise strengths $\sigma$. The system exhibits two distinct regimes: a quadratic growth of error at low $\sigma$, followed by a saturation phase ($\eta \longrightarrow 1$) indicating total loss of synchronization at high $\sigma$. Results are shown for a scale-free network ($N = 500$, $\gamma = 2.8$); similar behavior is observed for other network types.
• Figure 4: Synchronization transitions in empirical networks. Variation of order parameters $r_1$ and $r_2$ as functions of the coupling strength $K_2$ for the Zachary Karate Club network ($N=34$) (a,b) and the C. elegans network ($N=131$) (c,d) with $\beta = 0.2$. Blue solid and red dashed curves denote forward and backward numerical simulations, respectively. In both networks, $r_1$ displays a backward transition, whereas $r_2$ exhibits a reversible, continuous transition with overlapping forward and backward paths. Magenta curves in (b) and (d) represent the analytical values of $r_2$ obtained from the self-consistency equations. Cyan dots at $(0.3,1)$ indicate the target synchronization points. The analytical critical couplings (green dashed lines), $K_2^{c} \approx 0.0708$ for the Karate Club and $K_2^{c} \approx 0.0476$ for C. elegans, closely match the numerical estimates.
• Figure 5: Third-order harmonic extension of the coupled oscillator model. Order parameters $r_1$, $r_2$, and $r_3$ as functions of the coupling strength for a scale-free network with degree exponent $\gamma = 2.8$ and size $N = 500$, where phase frustration $0.5$. Blue and red curves represent forward and backward simulations obtained using the proposed optimal frequency distribution. The maximum values (complete synchronization) of $r_1$, $r_2$, and $r_3$ appear at the target point $(4,1)$ (cyan dots) when $K_1 = 0$ and $K_2 = 0$. While $r_1$ and $r_2$ display only backward transitions, $r_3$ exhibits clear hysteresis, consistent with the behavior observed in the bi-harmonic case.