Synchronization in Networks of Heterogeneous Kuramoto-Sakaguchi Oscillators with Higher-order Interactions

TL;DR

The paper investigates synchronization in globally coupled Kuramoto–Sakaguchi oscillators subject to phase frustration, additive noise, and higher-order ($2$-simplex) interactions. It combines numerical simulations with Ott–Antonsen reduction to derive a low-dimensional amplitude equation and perform fixed-point and bifurcation analyses. The key finding is that higher-order interactions widen the bistable region by shifting the backward saddle-node while leaving the forward critical point fixed, with noise and phase lag modulating hysteresis and synchronization via Kramers escape dynamics. The results provide a unified framework for frustrated, noisy, higher-order oscillator networks with potential applications to biological, physical, and engineered systems.

Abstract

How do the combined effects of phase frustration, noise, and higher-order interactions govern synchronization in globally coupled heterogeneous Kuramoto oscillators? To address this question, we investigate a globally coupled network of Kuramoto-Sakaguchi oscillators that includes both pairwise (1-simplex) and higher-order (2-simplex) interactions, together with additive stochastic forcing. Systematic numerical simulations across a broad range of coupling strengths, phase-lag values, and noise intensities reveal that synchronization emerges through a nontrivial interplay among these parameters. In general, weak frustration combined with mutually reinforcing coupling promotes synchronization, whereas strong frustration favors coherence under repulsive coupling. Forward and backward parameter sweeps reveal the coexistence of synchronized and desynchronized states. The presence and width of this bistable region depend sensitively on phase frustration, noise intensity, and higher-order coupling strength, with higher-order interactions significantly widening the bistable interval. To explain these behaviors, we employ the Ott-Antonsen reduction to derive a low-dimensional amplitude equation that predicts the forward critical point in the thermodynamic limit, the backward saddle-node point, and the width of the bistable region. Higher order interactions widen this region by shifting the saddle-node point without affecting the forward critical point. Further analysis of Kramer's escape rate explains how noise destabilizes coexistence states and diminishes bistability. Overall, our results provide a unified theoretical and numerical framework for frustrated, noisy, higher-order oscillator networks, revealing that synchronization is strongly influenced by the combined action of phase frustration, stochasticity, and both pairwise and higher-order interactions.

Figures (6)

• Figure 1: A representative architecture of the globally coupled oscillator system where the nodes ($0-simplexes$) are represented by magenta color, the pairwise connections, that is, the edges($1-simplexes$), are represented by black lines, and higher-order interactions ($2-simplexes$, filled triangles) are represented by orange color. Panel (a) represents all-to-all coupling, whereas panel (b) represents a smaller version of simplices.
• Figure 2: Synchronization landscapes in the $(k_1, \beta)$ parameter space for a network of $N = 100$ frustrated stochastic Kuramoto oscillators incorporating higher-order (2-simplex) interactions. The degree of collective coherence is measured by the order parameter $R$. The color map encodes the value of $R$, with dark blue ($R = 0$) corresponding to a fully desynchronized state and red ($R = 1$) indicating complete synchronization. Each column represents a distinct noise strength, $D = 0$, $0.5$, and $1$, while each row corresponds to a fixed value of the 2-simplex coupling strength, $k_2 \in {-20, -10, 0, 10, 20}$. Specifically, panels $(a$–$c)$ depict $k_2 = -20$, $(d$–$f)$$k_2 = -10$, $(g$–$i)$$k_2 = 0$, $(j$–$l)$$k_2 = 10$, and $(m$–$o)$$k_2 = 20$. The figure highlights how varying higher-order coupling and noise strength jointly modulate the onset and extent of synchronization within the frustrated oscillator ensemble.
• Figure 3: Variation of the order parameter $R$ as a function of the pairwise coupling strength $k_{1}$ for forward sweep (magenta color) and backward sweep (cyan color) simulations, at a fixed higher-order interaction strength $k_{2}=10$ for a system of $N=100$ oscillators. Each column corresponds to a different noise strength: panels $(a-g)$ correspond to $D=0$, panels $(b-h)$ are for the case when $D=0.5$, and panels $(c-i)$ correspond to $D=1.0$; each row corresponds to a different phase-lag value: panels $(a-c)$ for $\beta=\pi/6$, panels $(d-f)$ for $\beta=\pi/4$, and panels $(g-i)$ for $\beta=\pi/3$. The natural frequencies $\omega_i$ are drawn from a Lorentzian distribution with $\omega_0=0$ and width $\Delta=0.5$. The initial condition is an equally populated bi-cluster state in which half of the oscillators have phase $\phi=0$ and the remaining half of the oscillators have phase $\phi=\pi$. The figure demonstrates that the width of the hysteresis loop decreases systematically with increasing phase lag $\beta$ and noise strength $D$.
• Figure 4: Variation of the order parameter $R$ as a function of the pairwise coupling strength $k_{1}$ for forward sweep (magenta color) and backward sweep (cyan color) simulations, at a fixed higher-order interaction strength $k_{2}=20$ for a system of $N=100$ oscillators. Each column corresponds to a different noise strength: panels $(a-g)$ correspond to $D=0$, panels $(b-h)$ are for the case when $D=0.5$, and panels $(c-i)$ correspond to $D=1.0$; each row corresponds to a different phase-lag value: panels $(a-c)$ for $\beta=\pi/6$, panels $(d-f)$ for $\beta=\pi/4$, and panels $(g-i)$ for $\beta=\pi/3$. The natural frequencies $\omega_i$ are drawn from a Lorentzian distribution with $\omega_0=0$ and a width $\Delta=0.5$. The initial condition is an equally populated bi-cluster state in which half of the oscillators have phase $\phi=0$, and the remaining half of the oscillators have phase $\phi=\pi$. The figure demonstrates that the width of the hysteresis loop decreases systematically with increasing phase lag $\beta$ and noise strength $D$.
• Figure 5: Comparison of theoretical and numerical order parameters $R$ as a function of the pairwise interaction strength $k_1$ for $N = 100$ coupled stochastic SK oscillators for a fixed higher-order (2-simplex) interaction strength $k_2=10$. Here, $\omega_0 =0$ and $\Delta=0.5$. The theoretical value of $R$ (blue solid line) is obtained by solving Eq. (\ref{['Eq. 12']}), while the numerical values (blue circles) are computed by integrating Eq. (\ref{['Eq. 1']}). The first, second, and third columns correspond to noise strengths $D = 0,\, 0.5$, and $1$, respectively, whereas the rows correspond to phase-lag values $\beta = \pi/6,\, \pi/4$, and $\pi/3$. Panels $(a$–$c)$ represent $\beta = \pi/6$, $(d$–$f)$ correspond to $\beta = \pi/4$, and $(g$–$i)$ to $\beta = \pi/3$. The results demonstrate good agreement between the theoretical predictions and numerical simulations across all parameter regimes.