Stochastic Kuramoto oscillators with inertia and higher-order interactions

Synchronization is a fundamental phenomenon observed across physics, biology, chemistry, and engineering; from the rhythmic flashing of fireflies, the coordinated firing of neurons in the brain, and the stability of power grids [1]. The Kuramoto model provides a pivotal analytical tool for studying synchronization [2, 3], and demonstrating how a system of interacting oscillators with diverse natural frequencies begins moving in unison as the interaction strength varies. This model is particularly valued for its simplicity and analytical tractability, which not only make theoretical analyses feasible but also enhance its utility in practical applications such as power system [4] and biological systems [5]. Additionally, research has explored how the adaptation function and phase lag parameter affect the synchronized state [6, 7, 8, 9]. Introducing noise to the Kuramoto model incorporates stochasticity, which reflects the intrinsic fluctuations found in real-world systems. Thus, this approach enables more accurate and realistic simulations, effectively mimicking real-world environments [10]. Campa and Gupta considered the case of heterogeneous noise and analyzed the impact of noise strength on dynamical evolution of coupled Kuramoto oscillators [12]. There exist few other studies investigating the impact of noise strength on synchronization profile in networks of phase oscillators [13, 10], and on the model used to study stability of large human connectome graph [11].

Furthermore, incorporation of inertia in the coupled Kuramoto system provides an application for modeling power grids [14, 15, 16]. An inclusion of inertia term in the Kuramoto model has been shown to lead first-order phase transitions characterized by abrupt changes in system dynamics in response to a small change in coupling strength [17]. Managing external perturbations or noise is crucial to prevent systemic failures in many complex systems [18]. Acebrón et. al. studied the influence of noise on critical transition points at which systems dynamics exhibit significant changes [19, 20]. Gupta et. al. outlined a detailed phase space diagram for coupled Kuramoto oscillators with pairwise interactions [21, 22]. Additionally, Cao et. al. described the effects of noise on cluster explosive synchronization in second-order Kuramoto oscillators on networks [23].

All these results were limited to Kuramoto model (with or without inertia) with pairwise interactions. Recent studies have emphasized importance of higher-order interactions in modeling real-world complex systems [24, 25, 26]. Incorporating higher-order interactions into the Kuramoto model has been shown to lead abrupt (de)synchronization transition [27, 28] and tiered synchronization [29, 30] in contrast to second-order transitions observed in pairwise interactions. In absence of noise, the second-order Kuramoto model with higher-order interactions revealed the presence of prolonged hysteresis [31]. Furthermore, studies on the second-order Kuramoto model with higher-order interactions, incorporating phase lag and coupling strengths ranging from negative to positive, have demonstrated the emergence of synchronization and frequency chimera states [32].

This Letter considers the second-order Kuramoto model with higher-order interactions and Gaussian white noise. The study investigates noise-induced transitions in the model. Using Fokker-Plank equation, we first derive distribution function for frequency and phases for identical oscillators, and observe that the frequency order parameter is modulated by inertia and remains unaffected by coupling strength of both pairwise and higher-order interactions. Second, we show that an increase in noise strength shifts the forward (backward) critical point associated with an abrupt jump from an incoherent to coherent state (vice versa) toward higher coupling values. By employing a perturbation analysis, we obtain the expression for critical coupling strength in the forward direction which shows a dependence on inertia and noise strength. Additionally, in an overdamped system, we analytically predict all (un)stable states using the Ott-Antonsen approach. Also, in presence of higher-order interactions, note a shift from first-order to second-order phase transitions as noise strength increases.

We consider a stochastic Kuramoto model with $2-$simplex interactions and inertia. The equation of motion of $N$ globally coupled oscillators is given as,

$$\begin{split}m\ddot{\theta_{i}}=&-\dot{\theta_{i}}+{\Omega_{i}}+\frac{K_{1}}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i})\\ &+\frac{K_{2}}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}\sin(2\theta_{j}-\theta_{k}-\theta_{i})+\xi_{i}(t),\end{split}$$

(1)

where $\theta_{i}$ and $\Omega_{i}$ indicate the phase and intrinsic frequency, respectively, of $i^{th}$ Kuramoto oscillator. The parameters $K_{1}$ and $K_{2}$ denote the coupling strength for pairwise and $2-$simplex interactions, respectively, and $m$ is the inertia term. The noise term $\xi_{i}$ is defined as white Gaussian noise, with $\langle\xi_{i}(t)\rangle=0$ and $\langle\xi_{i}(t)\xi_{j}(s)\rangle=2D\delta_{ij}\delta(t-s)$, where $D$ represents the noise strength. The order parameter is defined as:

$${r_{p}e^{\iota\psi_{p}}=\frac{1}{N}\sum_{j=1}^{N}e^{\iota p\theta_{j}}},$$

(2)

Here, for $p=1$, $0\leq r_{1}\leq 1$ quantifies the magnitude of global synchronization, and $\psi_{1}$ is the mean phase of all oscillators. $r_{1}=0$ implies a state in which all oscillators move incoherently around the circle, while $r_{1}=1$ indicates global synchronization. Moreover, the frequency order parameter is defined as follows,

$${se^{\iota\phi}=\frac{1}{N}\sum_{j=1}^{N}e^{\iota\omega_{j}}},$$

(3)

where, $0\leq s\leq 1$ represents the magnitude of frequency synchronization. We express Eq. 1 in the mean-field form by using the order parameter (Eq. 2) such that each oscillator interacts with the mean phase of all the other oscillators,

$$\begin{split}m\ddot{\theta_{i}}=-\dot{\theta_{i}}+&{\Omega_{i}}+{K_{1}r_{1}}\sin(\psi_{1}-\theta_{i})\\ &+{K_{2}r_{2}r_{1}}\sin(\psi_{2}-\psi_{1}-\theta_{i})+\xi_{i}(t).\end{split}$$

(4)

Further, expressing Eq. 4 as a system of simultaneous first-order differential equations we get,

	$\displaystyle\dot{\theta_{i}}$	$\displaystyle=\omega_{i},$
	$\displaystyle\dot{\omega_{i}}$	$\displaystyle=\frac{1}{m}[-\omega_{i}+\Omega_{i}+K_{1}r_{1}\sin(\psi_{1}-\theta_{i})$
		$\displaystyle\quad+K_{2}r_{2}r_{1}\sin(\psi_{2}-\psi_{1}-\theta_{i})+\xi_{i}(t)].$		(5)

We employ the Fokker-Planck equation (FPE) to derive the probability density function $\rho(\theta,\omega,\Omega,t)$ for the oscillators,

$$\begin{split}\frac{\partial\rho}{\partial t}=\frac{D}{m^{2}}\frac{\partial^{2}\rho}{\partial\omega^{2}}-\frac{1}{m}\frac{\partial}{\partial\omega}[(-\omega+\Omega+K_{1}r_{1}\sin(\psi_{1}-\theta)\\ K_{2}r_{2}r_{1}\sin(\psi_{2}-\psi_{1}-\theta))\rho]-\omega\frac{\partial\rho}{\partial\theta}.\end{split}$$

(6)

In the stationary state $\frac{\partial\rho}{\partial t}=0$, let us first consider that all oscillators are identical $g(\Omega)=\delta(\Omega)$. Further, the density distribution is $2\pi$ periodic in $\theta$ and decays as $\omega\rightarrow\pm\infty$. Following the normalization condition $\int_{0}^{2\pi}\int_{-\infty}^{\infty}\rho(\theta,\omega,\Omega,0)d\omega d\theta=1$. We look for the solution of the density function in the form of $\rho(\theta,\omega)=\eta(\omega)\chi(\theta)$, ultimately obtaining the expression for the distribution function of frequencies and phases. Moreover, simulation of Eq. Model: shows that the frequency order parameter is unaffected by $K_{1}$ and $K_{2}$ (Figs. 1). Hence, solving Eq. 6 in the following manner,

$$\frac{D}{m^{2}}\frac{d^{2}\eta}{d\omega^{2}}+\frac{\omega}{m}\frac{d\eta}{d\omega}+\frac{\eta}{m}=0.$$

(7)

$$\begin{split}-\omega\eta\frac{d\chi}{d\theta}-\frac{1}{m}&\Bigl{[}K_{1}r_{1}\sin(\psi_{1}-\theta)+K_{2}r_{2}r_{1}\\ &\sin(\psi_{2}-\psi_{1}-\theta)\Bigr{]}\chi\frac{d\eta}{d\omega}=0.\end{split}$$

(8)

By solving these differential equations, and applying boundary and normalization conditions, the distribution functions of frequencies and phases can be represented as follows,

$$\eta(\omega)=\sqrt{\frac{m}{2\pi D}}e^{\frac{-m\omega^{2}}{2D}}.$$

(9)

$$\chi(\theta)=\frac{e^{K_{1}r_{1}\cos(\psi_{1}-\theta)+K_{2}r_{2}r_{1}\cos(\psi_{2}-\psi_{1}-\theta)}}{\int_{0}^{2\pi}e^{K_{1}r_{1}\cos(\psi_{1}-\theta)+K_{2}r_{2}r_{1}\cos(\psi_{2}-\psi_{1}-\theta)}d\theta}.$$

(10)

Fig. 1(d) illustrates the frequency distribution function derived analytically using Eq. 9 (solid lines), which match with the numerically obtained results from Eq. Model: (symbols).

In the incoherent state, phases are randomly distributed around the unit radius circle implying $r_{1}=r_{2}=0$. Using Eq. 6 to obtain the density function of the incoherent state given as,

$$\rho_{0}(\omega,\Omega)=\frac{1}{2\pi}\sqrt{\frac{m}{2\pi D}}e^{\frac{-m(\omega-\Omega)^{2}}{2D}}.$$

(11)

Giving a small perturbation ($\epsilon<<1$) to the incoherent state, we get

$$\rho(\theta,\omega,\Omega,t)=\rho_{0}(\omega,\Omega)+\epsilon\alpha(\theta,\omega,\Omega,t)+O(\epsilon^{2}).$$

(12)

The order parameter in the continuum limit can be written as,

$$\begin{split}r_{1}e^{\iota\psi_{1}}&=\int_{0}^{2\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{\iota\phi}\rho(\phi,\omega,\Omega,t)g(\Omega)d\Omega d\omega d\phi.\\ r_{2}e^{\iota\psi_{2}}&=\int_{0}^{2\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{2\iota\phi}\rho(\phi,\omega,\Omega,t)g(\Omega)d\Omega d\omega d\phi.\end{split}$$

(13)

Further, incorporating Eqs.( 12 and  13) into Eq. 6, and equating like terms in $\epsilon$ we obtain

$$\frac{\partial\alpha}{\partial t}+\omega\frac{\partial\alpha}{\partial\theta}-\frac{D}{m^{2}}\frac{\partial^{2}\alpha}{\partial\omega^{2}}-\frac{1}{m}\frac{\partial}{\partial\omega}\left((\omega-\Omega)\alpha\right)=-\frac{K_{1}}{m}\frac{\partial\rho_{0}}{\partial\omega}\int_{0}^{2\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\sin{(\phi-\theta)}\alpha(\phi,\omega,\Omega,t)g(\Omega)d\Omega d\omega d\phi.$$

(14)

This equation is exactly the same as in Ref. [20] (Eq. 9), and therefore the subsequent derivations are also the same as in Ref. [20]. For Lorentzian frequency distribution $g(\Omega)=\frac{\Delta}{\pi[(\Omega-\Omega_{0})^{2}+\Delta^{2}]}$ with mean $\Omega_{0}=0$ and standard deviation $\Delta$, we get $K_{c}=2\Delta(m\Delta+1)+\frac{2(2+3m\Delta)}{2+m\Delta}D+O(D^{2})$ for the forward transition point where oscillators jump from an incoherent to a synchronized state, applicable when noise strength is small ($D\ll 1$). For a negligible inertia $m$, $K_{c}=2(D+\Delta)$ is similar to the case of the Kuramoto Model with Gaussian white noise [33]. We find that $K_{c}$ is independent of $K_{2}$.

Numerical simulations are performed using the Euler method with a time step $dt=0.05$, for Eq. 1 by converting it to mean-field equation and then simultaneous first-order differential equations Eq. Model:. Results are obtained for $N=20,000$ and $r_{1}$ averaged over $8\times 10^{4}$ iterations after removing an initial transient by changing $K_{1}$ in the forward and backward directions.

Fig. 2 depicts the phase transition for low $m$ value indicating that critical transition point $K_{c}$ shifts towards higher positive coupling value with an increase in $D$. Fig. 2(a) represents the second-order phase transition in absence of higher-order interactions ($K_{2}=0$). We consider the Lorentzian frequency distribution with standard deviation $\Delta=0.5$, so $K_{c}=2D+1$ correctly determines the transition point, for ($D<1$). Furthermore, Fig. 2(b) for $K_{2}=5$ demonstrates that phase transition changes from first-order to second-order as we increase $D$, despite the presence of higher-order interactions.

For $m=1$, Fig. 3 (a, b) depicts the first-order phase transition. With an increase in $D$, $K_{c}$ shifts towards higher value. Also, the backward transition point $K_{1b}$ at which the oscillators jump from synchronized to incoherent states shifts towards higher value. Further, Fig. 3 (b) depicts that $K_{c}$ is not affected by the presence of $2-$simplex interactions term, as yielded by the perturbation analysis of the incoherent state that expression of $K_{c}$ comes out to be independent of $K_{2}$. Since $K_{2}$ impacts $K_{1b}$, thereby leading to an increase in the hysteresis width compared to the system with only pairwise interactions.

Fig. 4 analyzes forward $K_{c}$ and backward $K_{1b}$ transition points as a function of noise strength $D$, and examines the impact of changes in $m$, $K_{2}$, $\Delta$. Fig. 4(a) depicts that an increase in inertia $m$ shifts $K_{c}$ towards higher values, while $K_{1b}$ remains unaffected thereby increasing the hysteresis width. Fig. 4(b) demonstrates that increasing $K_{2}$ leads to a shift in $K_{1b}$ towards lower values with $K_{c}$ remaining the same. This also results in an increased hysteresis width. Fig. 4(c) illustrates that an increase in $\Delta$ value shifts both $K_{c}$ and $K_{1b}$ towards higher values. As reflected by Fig. 4(d) hysteresis width $W=K_{c}-K_{1b}$ decreases as $D$ increases.

In overdamped systems, where the inertia $m$ is negligible, the Ott-Antonsen approach [34] facilitates the dimensional reduction for intrinsic frequency drawn from the Lorentzian distribution as:

$$\dot{r_{1}}=-(D+\Delta)r_{1}+\frac{K_{1}(r_{1}-r_{1}^{3})}{2}+\frac{K_{2}(r_{1}^{3}-r_{1}^{5})}{2}.$$

(15)

Eq. 15 is derived by equating the coefficient of $e^{\iota\theta}$ in the context of pairwise interactions only (as discussed in [35]). We extend the analysis to include $2$-simplex interactions. The forward critical transition point $K_{c}=2(D+\Delta)$, is a threshold where oscillators jump from an incoherent to a synchronized state, and for $K_{1}>K_{c}$, $r_{1}=0$ solution becomes unstable. By equating the coefficients of $e^{\iota n\theta}$, where $n\in\mathbb{N}$, we get $\dot{r_{1}}=-(nD+\Delta)r_{1}+\frac{K_{1}(r_{1}-r_{1}^{3})}{2}+\frac{K_{2}(r_{1}^{3}-r_{1}^{5})}{2}$. As we can see that for the noisy case $n$ arrives explicitly in $\dot{r_{1}}$ expression, hence Ott-Antonsen approach does not result in the dimensional reduction of Eq. 1 for overdamped system. However, by substituting $n=1$ in $\dot{r_{1}}$, the analytical predictions closely matches with those achieved through numerical simulations (Fig. 5).

Fig. 5 (a) $K_{2}=0$ shows second-order phase transition, and (b) $K_{2}=5$ manifests first-order phase transition. Increasing $D$ results in decrease in the hysteresis width, finally yielding second-order phase transition.

This study investigates noise-induced phase transitions in the Kuramoto model with inertia and $2$-simplex interactions. First, using Fokker-Plank equation, we determine the distribution function of frequency and phase for identical oscillators, demonstrating that the frequency order parameter while depends on inertia $m$, remains unaffected by coupling strengths $K_{1}$ and $K_{2}$. Our analysis further reveals that noise strength is a critical factor governing both forward and backward transition points. Using the perturbation analysis for Lorentzian frequency distribution, we show that the forward transition point is independent of the $2$-simplex interaction term and rather depends on inertia and noise strength. By employing the Ott-Antonsen approach we derive approximate analytical predictions for the overdamped system. A significant finding is that as noise strength increases in the overdamped system, the phase transition shifts from the first-order to the second-order even in the presence of higher-order interactions. These results underscore the interplay between noise, inertia, and higher-order interactions in synchronized systems.

We hope that understanding these dynamical behaviours may be useful for studying stability of power grid systems. By incorporating noise and higher-order interactions into the Kuramoto model, we can more accurately predict critical transition points which are essential for preventing systemic failures. This improved modeling approach can aid in designing of more robust power grid systems, capable of maintaining stability under perturbations and noise conditions. Furthermore, few direct possible extensions of this work are to have more realistic noise models [36] and incorporation of multiplicative noise [37].

SJ and PR acknowledge Govt of India SERB Power grant SPF/2021/000136 and PMRF grant PMRF/2023/2103358, respectively.