Tiered synchronization in coupled oscillator populations with interaction delays and higher-order interactions

In the interdisciplinary study of collective behavior and synchronization, the Kuramoto model is of pivotal importance due to its analytical tractability and versatility for modeling a wide range of behaviors.⁹ This versatility comes in part from its ability to incorporate a wide range of properties critical for different physical and biological systems. One such property is the presence of interaction delays,^10–12 which give rise to rich dynamics and multistability. Another property that has attracted more attention in recent years is the presence of higher-order interactions,^13–17 most notably motivated by applications of neuroscience^18–21 and physics.^22,23 In fact, the effect of higher-order interactions has already been investigated in both synchronization^24–36 and other kinds of collective behavior.^37–40 

Our understanding of the macroscopic dynamics of coupled oscillator populations with interactions delays and higher-order interactions has been further developed by applying the dimensionality reduction of Ott and Antonsen.^41,42 In the case of heterogeneous interaction delays, as the characteristic time delay and the mean natural frequency are increased, the onset of synchronization is likewise made larger, eventually transitioning from a supercritical bifurcation to a subcritical one as the curve of steady-state solutions folds over itself and a saddle-node bifurcation is born.^7,43–46 On the other hand, in systems with higher-order interactions, the higher-order interactions themselves do not alter the location of the onset of synchronization or the stability of the incoherent state but promote synchronization via nonlinear terms.⁸ Increasing the higher-order coupling strength eventually causes the curve of steady-state solutions to similarly fold over itself as a saddle-node bifurcation is born. Thus, in both cases, bistability emerges between the incoherent and synchronized states.

In this paper, we study populations of coupled phase oscillators with both interaction delays and higher-order interactions. The dynamics incorporate dyadic, triadic, and tetradic interactions, with heterogeneous time delays between each pair of oscillators. After applying the dimensionality reduction of Ott and Antonsen and analyzing the reduced macroscopic dynamics, we find that bistability remains a key feature, but the combination of interaction delays and higher-order interactions promotes an additional mechanism than that described above. Specifically, in addition to bistability between the incoherent state and the synchronized state owing to a subcritical pitchfork and a saddle node, we also observe a pair of saddle nodes that leaves the onset of synchronization supercritical, leading to tiered synchronization, i.e., bistability between a weakly synchronized state and a strongly synchronized state.

The remainder of this paper is organized as follows. In Sec. II, we present the governing equations and the dimensionality reduction using the Ott–Antonsen ansatz. In Sec. III, we present an analysis of the steady-state dynamics. In Sec. IV, we analyze the incoherent state. In Sec. V, we consider the special case of eliminating triadic coupling, leading to further analytical results and a sketch of an illustrative bifurcation diagram. In Sec. VI, we conclude with a discussion of our results.

In this work, we consider populations of coupled phase oscillators with both interaction delays and higher-order interactions. We consider heterogeneous time delays between oscillators⁷ with dyadic, triadic, and tetradic interactions,⁸ yielding

where $θi$ and $ωi$ are the phase and natural frequency of oscillator $i$⁠; $K1$⁠, $K2$⁠, and $K3$ are the respective 1-, 2-, and 3-simplex coupling strengths; and $τij$ is the interaction delay between oscillators $i$ and $j$ felt by oscillator $i$⁠. In general, we assume that natural frequencies and time delays are drawn from their respective distributions $g(ω)$ and $h(τ)$⁠. While the degree of synchronization is measured by the magnitude $r$ or the classical instantaneous Kuramoto order parameter, given by

we also define the so-called Daido order parameter,

as well as two different varieties of oscillator-specific time-delayed order parameters, given by

and

Using Eqs. (4) and (5), we may rewrite Eq. (1) as

where $∗$ denotes complex conjugate.

Next, we seek to derive a closed-form system governing the dynamics of the order parameter for which purpose, we consider the continuum limit of infinitely many oscillators, $N→∞$⁠. Note that in this limit, we may express the order parameter as the integral

and, similarly, the Daido order parameter as

where $f(ω,θ,t)$ is the density function that describes the fraction $f(ω,θ,t)dθdω$ of oscillators with phase and frequency, respectively, in $[θ,θ+dθ)$ and $[ω,ω+dω)$ at time $t$⁠. The time-delayed order parameters may then be written as

and

Importantly, Eqs. (9) and (10) reveal that in the continuum limit and at steady-state, the variation between the time-delayed order parameters vanish across different oscillators; therefore, we may drop the subscripts, i.e., $wi(1)=w(1)$ and $wi(2)=w(2)$ for all $i$⁠; that is, in the continuum limit, each time-delayed order parameter $wi(t)$ is defined by the same mean field via $z(t)$⁠. With this simplification of Eq. (6), we note that the density function must have a Fourier series that takes the form

where c.c. denotes the complex conjugate of the preceding term. Next, using the Ott–Antonsen ansatz,^41,42 which essentially posits solutions with geometrically decaying Fourier coefficients, i.e., $f^n(ω,t)=αn(ω,t)$⁠, all Fourier modes remarkably reduce to a single differential equation for the function $α$ given by

To connect the dynamics of $α$ to the order parameter, we consider the case of Lorentzian-distributed frequencies, letting

where $Δ$ and $ω0$ give the spread and mean of the natural frequencies. Specifically, for the proposed density function $f$⁠, Eq. (7) may first be integrated in the $θ$ direction, yielding

Equation (14) may further be evaluated using the Cauchy residue theorem, taking advantage of the simple pole of the frequency distribution $g(ω)$ at $ω=ω0−iΔ$⁠, resulting in

Thus, by evaluating Eq. (12) at $ω=ω0−iΔ$ and taking a complex conjugate, we obtain the following differential equation for $z$⁠:

To close the dynamics, we now seek differential equations for $w(1)$ and $w(2)$⁠. We begin by noting that Eqs. (9) and (10) may be rewritten as

and

where $z2(t)=∬f(ω,θ,t)e2iθ(t)dθdω$ is the Daido order parameter, which following the dimensionality reduction above is simply given by $z2(t)=z2(t)$⁠. By considering the special case of exponentially distributed time delays, namely, letting

so that the characteristic time delay between oscillators is given by $T$⁠, Eqs. (17) and (18) may be treated with the Laplace transform to obtain the following differential equations:

and

Thus, Eqs. (16), (20), and (21) constitute a closed system for the dynamics of the instantaneous and time-delayed order parameters.

To proceed with our analysis of the low dimensional dynamics given by Eqs. (16), (20), and (21), we begin by rewriting the dynamics in polar coordinates, yielding

We seek steady-state solutions where the global order parameters reach a fixed amplitude, $r˙=ρ˙(1)=ρ˙(2)=0$⁠, with a phase that processes at a constant rate, $ψ˙=ϕ˙(1)=Ω$ and $ϕ˙(2)=2Ω$⁠. (Note that $ϕ(2)$ processes at twice the velocity as $ψ$ and $ϕ(1)$ since $w(2)$ chases $z2$⁠, which processes with twice the angular velocity of $z$⁠.) Applying this to the time-delayed order parameters [Eqs. (24) and (25)], we get

and after using the trigonometric identity $cos2⁡x+sin2⁡x=1$⁠, we obtain

A similar treatment of Eqs. (26) and (27) yields

and subsequently,

Before returning to Eqs. (22) and (23), it will also be convenient to eliminate the trigonometric quantities in those equations. The terms $cos⁡(ϕ(1)−ψ)$ and $sin⁡(ϕ(1)−ψ)$ may be eliminated simply using Eqs. (28)–(30), yielding

The terms $cos⁡(ϕ(2)−ϕ(1)−ψ)$ and $sin⁡(ϕ(2)−ϕ(1)−ψ)$⁠, on the other hand, require some more care and the use of the trigonometric identities $cos⁡(x+y)=cos⁡(x)cos⁡(y)−sin⁡(x)sin⁡(y)$ and $sin⁡(x+y)=sin⁡(x)cos⁡(y)+cos⁡(x)sin⁡(y)$⁠. Rewriting the argument $ϕ(2)−ϕ(1)−ψ=(ψ−ϕ(1))+(ϕ(2)−2ψ)$ and using Eqs. (28)–(35) yields

We note that Eqs. (30) and (33)–(37) may also be derived using Eqs. (17) and (18). We present this alternative derivation in the  Appendix.

We now have the ingredients necessary for returning to Eqs. (22) and (23). Specifically, using Eqs. (30) and (33)–(37), seeking the stationary state in Eqs. (22) and (23) yields

Equations (38) and (39) characterize the degree of synchronization via the amplitude $r$ of the order parameter and the angular velocity $Ω$ of the synchronized state depending on the coupling strengths $K1$⁠, $K2$⁠, and $K3$⁠, the characteristic time delay $T$⁠, and the width of the frequency distribution $Δ$⁠.

To verify the dimensionality reduction and steady-state bifurcation analysis, we now compare results from simulation to the solutions predicted by Eqs. (38) and (39). To overcome the numerical complexity of incorporating explicit time delays in simulations, we consider the dynamics rewritten as Eq. (6) with the time-delayed order parameter dynamics given by Eqs. (20) and (21). In particular, we consider the system dynamics as dyadic coupling $K1$ is adiabatically increased and then decreased for three combinations of higher-order coupling, $K2=0$ and $K3=0.7$⁠, $K2=K3=0.35$⁠, and $K2=0.7$ and $K3=0$⁠, and plot the simulation results using forward and backward triangles, respectively, in Figs. 1(a)–1(c). Simulations use $N=104$ oscillators and at each different value of $K1$ run through a transient of $2×105$ time steps with $Δt=4×10−3$ and then average $r$ over $2×105$ time steps. We then plot the analytical predictions given by Eqs. (38) and (39), solved numerically, using solid and dashed curves, indicating stability and instability, respectively. Other parameters are fixed at $T=0.4$⁠, $ω0=7.8$⁠, and $Δ=1$⁠. We note here that the small values of $K2$ and $K3$ compared to $K1$ are in line with phase reduction analyses that tend to characterize higher-order coupling as small in comparison with dyadic coupling.

The results plotted in Fig. 1 demonstrate a rich set of dynamics that come from the combination of time delays and higher-order interactions. In particular, for each choice of higher-order coupling used, the dynamics admit ranges of bistability. However, unlike the typical scenarios so far observed in coupled oscillator systems with time-delayed interactions and in coupled oscillator systems with higher-order interactions where bistability occurs between the incoherent state and a strongly synchronized state, here, there are regions of bistability between a weakly synchronized state and a strongly synchronized state. i.e., rather than bistability occurring directly from a subcriticality and a saddle node; here, supercriticality is maintained at the onset of synchronization and bistability emerges from the formation of a pair of saddle-node bifurcations. Thus, a tiered synchronization profile emerges, where a stable, weakly synchronized state, characterized by relatively small $r$⁠, exists beyond the supercritical Hopf bifurcation at the onset of synchronization and a stable, strongly synchronization state, characterized by larger values of $r$⁠, exists after the curve folds over onto itself twice through a pair of saddle-node bifurcations. Moreover, these two stable synchronized branches are connected via an unstable branch.

For comparison, we also plot in Fig. 2(a) similar results and parameters, but for the case of dyadic coupling only, i.e., with $K2=K3=0$⁠. We point out that the lack of higher-order interactions leads to a far less pronounced region of bistability. However, bistability does, in fact, exist for this choice of parameters, albeit in a very thin region of the coupling strength $K1$⁠. In Fig. 2(b), we present a zoomed-in view of this region. The presence of this bistability deserves a few remarks. First, to our knowledge, bistability of this nature (i.e., owing to a pair of saddle nodes after a supercritical pitchfork–Hopf bifurcation) has not been observed in systems with only interaction delays (i.e., without higher-order interactions). Second, these bistability regions appear to be so thin that they are virtually unobservable in direct simulations of oscillator systems, even for large enough systems where finite-fluctuations have been all but eliminated. Note that the simulation results in Fig. 2(a) show no trace of bistability—for this, we need the analytically predicted curve. Taking these two points together, it appears that the presence of higher-order interactions remains an important ingredient for bistability of this nature in oscillator systems. In Secs. IV and V, we, respectively, analyze the incoherent state and consider a special case that allows us to sketch the bifurcation diagram of the system.

Before analyzing the nonlinear effects that are present in the dynamics, we first focus our attention on the incoherent state described by $z=w(1)=w(2)=0$⁠. Note that our steady-state equation (38) implies that the incoherent state is always a solution, and therefore, the stability properties of the incoherent state determine the onset of synchronization. By eliminating the incoherent state from Eq. (38) and letting $r→0+$ in Eqs. (38) and (39), we obtain a simplified set of equations describing this critical point corresponding to the collision between the synchronized an incoherent branches in Fig. 1, given by

Equations (40) and (41) may be combined to find

which may be inserted back into Eq. (40) and solved for $K1$ to yield the critical dyadic coupling strength

We note that from Eq. (43), we see that this critical coupling strength depends monotonically on both the mean natural frequency $ω0$ and the characteristic time delay $T$⁠, but not the width of the frequency distribution, $Δ$⁠; i.e., increasing either $ω0$ or $T$ causes $K1c$ to increase, but increasing $Δ$ may increase or decrease $K1c$⁠.

Alternatively, it is easy to check that this critical coupling strength corresponds exactly to the first crossing of the eigenvalues of the complex-valued Jacobian of Eqs. (16), (20), and (21) for the incoherent state, given by

Specifically, the incoherent state is asymptotically stable (with all eigenvalues located in the left-half complex plane) for $K1<K1c$⁠, after which stability is lost at $K1=K1c$⁠. Given the rotating nature of solutions (described by the angular velocity $Ω$⁠), this is a Hopf bifurcation that may be either supercritical or subcritical depending on the nature of the steady-state solutions. (In Fig. 1, it is supercritical for all thee cases.) However, in the appropriate rotating reference frame, it may viewed as a pitchfork bifurcation. In Fig. 1, we delineate the stable and unstable portions of the incoherent branch in solid and dashed curves along $r=0$⁠. Importantly, we note that the onset of synchronization is unaffected by the triadic and tetradic coupling strengths, $K2$ and $K3$⁠, implying that in terms of the macroscopic dynamics, the higher-order interactions offer only nonlinear effects.

We now turn our attention to the special case where triadic coupling is eliminated, i.e., we set $K2=0$⁠, leaving only dyadic and tetradic coupling. As we will see, this special case allows for a modest simplification of the general case, thereby allowing for further analytical results and a fuller picture of the bifurcation diagram. We begin by noting that when $K2$ is set to zero, Eqs. (38) and (39) reduce to

which may be combined to solve for $Ω$⁠, yielding

Inserting Eq. (47) into Eq. (45) and solving for $K1$ then yields

While Eq. (48) describes $K1$ as a function of $r$ (rather than vice versa), it turns out to be extremely useful for exploring the system dynamics. First, we use it to plot multiple synchronization profiles in Fig. 3 for $K3=−0.5$⁠, $1$⁠, $2$⁠, $3$⁠, and $4$ (red to purple, right to left). Other parameters are the same as those used in Fig. 1 except for $K2=0$⁠. In fact, these example choices of $K3$ sweep out a collection of curves that illustrate a rich range of system dynamics. In particular, starting at $K3=−0.5$ (right-most curve), the system undergoes a single supercritical Hopf bifurcation delineating the incoherent a synchronized states. As $K3$ is then increased, the dynamics qualitatively change, as the top portion of the curve develops a fold over onto itself but then folds back. This folding corresponds to a pair of saddle node bifurcation and separates a single stable synchronized branch into two stable synchronized branches, corresponding to weak and strong synchronization. These branches are connected by another branch that is unstable. This is illustrated by the $K3=1$⁠, $2$⁠, and $3$ curves. Finally, for even larger $K3$⁠, the bottom saddle-node disappears as the weakly synchronized branch vanishes as the onset of synchronization becomes a subcritical Hopf bifurcation, as illustrated by the $K3=4$ curve.

Moreover, while Eq. (48) does not characterize $r$ as a function of $K1$⁠, but rather vice versa, for purposes of sketching the bifurcation diagram of the system, this format is quite convenient since the saddle-node bifurcations described in Figs. 1 and 3 may be found and described using the derivative $∂K1∂r$⁠. Treating Eq. (43) analytically remains difficult, but the appropriate bifurcation conditions may be easily identified numerically by searching for local minima and maxima of $K1$ as a function of $r$⁠. In Fig. 4, we plot the bifurcation diagram over dyadic coupling $K1$ and tetradic coupling $K3$⁠. Bifurcation curves indicate Hopf and saddle-node bifurcations, labeled accordingly. In fact, the pair of saddle-node bifurcations is born at a codimension-two point, indicated by the lower black circle. The saddle-node curves split as $K3$ increases, and eventually, the higher saddle-node curve collides with the pitchfork bifurcation at another codimension-two point, above and below which the Hopf is subcritical and supercritical.

In addition to providing a representatively full picture of the macroscopic dynamics, the bifurcation diagram in Fig. 4 demonstrates the point made previously where bistability owing to a supercritical Hopf bifurcation followed by a pair of saddle nodes can be achieved in the absence of higher-order interactions but, in fact, is so subtle that it is difficult to observe in simulations and occurs for such a thin region of $K1$⁠. In particular, the codimension-two point at the formation of the two saddle nodes in Fig. 4 lies just under the $K3=0$ line, indicating that at $K3=0$⁠, i.e., in the absence of higher-order interactions, the pair of saddle nodes exists (for the parameters chosen). Again, we emphasize that to our knowledge, this kind of transition has not been observed in previous work, likely due to the thin parameter range of bistability (see also Fig. 1).

In this paper, we have studied the synchronization dynamics of populations of coupled phase oscillators with both interaction delays and higher-order interactions. After employing the dimensionality reduction of Ott and Antonsen,^41,42 we presented an analysis of the steady-state solutions. We showed that the combination of interactions delays and higher-order interactions promote a bistability that, unlike what has been observed previously when either interactions delays or higher-order interactions are present, leaves the onset of synchronization supercritical and creates two stable synchronized states via a pair of saddle nodes.

We also highlight two other pieces of interest that emerge from our analysis. First, as is evident from the contribution of triadic and tetradic coupling strengths in low dimensional equations, the presence of higher-order interactions offers only nonlinear effects to the macroscopic system dynamics. This is further emphasized by the fact that the onset of synchronization does not depend on the triadic or tetradic coupling strengths. Second, while a mechanism for bistability that gives rise to tiered synchronization discussed above is promoted by the combination of interactions delays and higher-order interactions, we find that the double saddle node can, in fact, be observed with only interaction delays. It appears that the regions of parameter space that give rise to such transitions are simply very small and have not been previously observed.

Last, we have considered here the case of heterogeneous oscillators and heterogeneous time delays, where natural frequencies and time delays are distributed via a Lorentzian distribution and an exponential distribution, respectively. While these choices were made in order to facilitate analytical treatment of the system, specifically aiding in the derivation of the low dimensional dynamics, it remains to be seen whether different choices of natural frequency distributions and/or interaction delays qualitatively change the macroscopic dynamics, a task left for future work.

P.S.S. acknowledges support from the NSF grant (No. MCB-2126177). C.X. acknowledges the National Natural Science Foundation of China (Grant No. 11905068) and the Scientific Research Funds of Huaqiao University (Grant No. ZQN-810).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

We begin our alternative derivation with Eq. (17). Using the polar form for $z$ and $w(1)$⁠, choosing $ψ=Ωt$⁠, and denoting $φ(1)=ψ−ϕ(1)$⁠, we have that

Taking an absolute value then recovers Eq. (30) in the main text. Moreover, taking real and imaginary parts, we recover Eqs. (34) and (35).

On the other hand, if we begin with Eq. (18), we may similarly write

Taking an absolute value then recovers Eq. (33) in the main text. Moreover, taking real and imaginary parts, we obtain

which can be used along with Eqs. (34) and (35) to yield Eqs. (36) and (37) in the main text.