We study the effect of structured higher-order interactions on the collective behavior of coupled phase oscillators. By combining a hypergraph generative model with dimensionality reduction techniques, we obtain a reduced system of differential equations for the system’s order parameters. We illustrate our framework with the example of a hypergraph with hyperedges of sizes 2 (links) and 3 (triangles). For this case, we obtain a set of two coupled nonlinear algebraic equations for the order parameters. For strong values of coupling via triangles, the system exhibits bistability and explosive synchronization transitions. We find conditions that lead to bistability in terms of hypergraph properties and validate our predictions with numerical simulations. Our results provide a general framework to study the synchronization of phase oscillators in hypergraphs, and they can be extended to hypergraphs with hyperedges of arbitrary sizes, dynamic-structural correlations, and other features.

Synchronization of networks of coupled oscillators is one of the most iconic problems in complex systems, with applications in biology,^{1,2} physics,^{3} and engineering.^{4–6} Usually, coupling between oscillators is assumed to be mediated by pair interactions. Recently, motivated by applications in physics^{7,8} and biology,^{1,9} there has been much interest in studying the effect of higher-order interactions, i.e., simultaneous interactions between multiple oscillators, on synchronization patterns.^{10–15} In this paper, we study the synchronization of coupled phase oscillators on complex hypergraphs. We use a hypergraph generative model and develop a mean-field analysis using dimensionality reduction techniques to obtain low-dimensional descriptions of synchronization in terms of hypergraph structural parameters. We find conditions on the hypergraph that result in bistability. Our results provide a general and flexible framework to study synchronization on hypergraphs.

## I. INTRODUCTION

Synchronization processes are present in many applications.^{16,17} Some common examples include collections of flashing fireflies,^{18,19} crickets chirping in unison,^{16} neuronal networks,^{20} cortical brain rhythms,^{1} and power grid dynamics.^{21,22} A paradigmatic model for synchronization is the Kuramoto model of phase oscillators,^{23,24} in which synchronization is mediated by pairwise interactions between oscillators. The Kuramoto model on complex networks has many applications and is one of the central models in complex science.^{25–27} Recently, with motivation from fundamental principles^{7,8} and applications to neuroscience,^{1,9} there has been much attention devoted to synchronization in networks with higher-order interactions, i.e., simultaneous interactions between multiple nodes. Higher-order interactions in coupled phase oscillator systems result in interesting phenomena like abrupt switching between incoherent and synchronized states, hysteresis, and bistability.^{10,12} So far, most of the analytical results have been obtained for the all-to-all coupling case, and there is no clear way to predict the effect of complex interaction structure on these phenomena. In this paper, we study the synchronization of phase oscillators on complex hypergraphs, i.e., networks with higher-order interactions and non-trivial connectivity. To do so, we restrict our attention to a specific but flexible hypergraph generative model that allows us to generate and study hypergraphs with tunable characteristics. By using this generative model in combination with the Ott–Antonsen ansatz,^{28} we are able to obtain low-dimensional descriptions of the system’s order parameters in terms of the hypergraph’s structural properties. We illustrate our approach with two examples of a hypergraph with interactions of sizes 2 (links) and 3 (triangles): a random hypergraph and a hypergraph constructed in such a way that the numbers of links and triangles at each node are correlated. We derive analytical conditions on the properties of these hypergraphs that result in synchronization, incoherence, or bistable behavior and validate our results with numerical simulations.

This paper is organized as follows. In Sec. II, we present our hypergraph generative model and the Kuramoto model on hypergraphs. In Sec. III, we use the Ott–Antonsen ansatz and a mean-field approximation to obtain low-dimensional descriptions for the local and global order parameters. In Sec. IV, we demonstrate our framework on two example hypergraphs. In Sec. V, we discuss our results and their limitations.

## II. MODEL

In this section, we introduce the hypergraph generative model and the Kuramoto model on hypergraphs.

### A. Hypergraph model

A hypergraph is a pair of nodes and hyperedges $(V,E)$, where $V$ is the set of nodes labeled $n=1,2\u2026,N$, and the set of hyperedges $E$ is a set of subsets of $V$. The $m$th order degree of a node $n$ is given by $kn(m)$, which gives the number of hyperedges with size $m$ that node $n$ is a part of. The *hyperdegree* of node $n$ is given by $kn={kn(1),kn(2),\u2026,kn(M)}$, where $M$ is the largest hyperedge size. For simplicity, we refer to hyperedges of sizes $2$ and $3$ as *links* and *triangles*, respectively. We denote by $N(k)$ the number of nodes with hyperdegree $k$, and define the hyperdegree distribution as $P(k)=N(k)/N$.

We will consider synchronization on a class of hypergraphs produced by the following generative model. For a given set of nodes $n=1,2,\u2026,N$ and a specified vector of target hyperdegrees $[k1,k2,\u2026,kN]$, the hyperedge ${i1,i2,\u2026,im}$ is created with probability $a(m)(ki1,ki2,\u2026,kim)$. By counting the expected number of hyperedges of size $m$ in two different ways, one finds that the functions $a(m)$ should be normalized such that

This model is a natural extension of latent feature models^{29} to hypergraphs and allows us to generate hypergraphs with heterogeneous and correlated hyperdegree distributions.^{30} The model can be easily extended to the case where hyperedges connect preferentially nodes with certain attribute variables such as nodal community index, oscillator frequency, or other dynamical parameters. On the other hand, the generative model is not able to capture features beyond the preference for hyperedges to connect certain types of nodes. An important class of hypergraphs that is not captured by this generative model is that of simplicial complexes, where triangles only connect triads of nodes that form a clique with pairwise connections (for simplicial complex generative models, see for example, Refs. 31–34).

Our subsequent results will apply to the “expected” network generated from this generative model. Such an approach is similar to the analysis of network processes based on the configuration model (e.g., Ref. 35) or the *annealed network approximation*.^{36–38} A similar approach has been successfully applied to the Kuramoto model on pairwise networks.^{39} The limitations of this approach are discussed in Sec. V.

### B. Higher-order Kuramoto model

The Kuramoto model of phase synchronization can be generalized to account for higher-order interactions in different ways. In Refs. 13–15, the synchronization of phases defined on the faces of a simplicial complex is studied. Here, following Refs. 10 and 12, we will instead consider synchronization mediated by the simultaneous, nonlinear interaction of all the phases belonging to the edges of a hypergraph. In this context, the Kuramoto model for the phases $\theta n$ of nodes $n=1,2,\u2026,N$ on a hypergraph $H$ can be generalized to

where the coupling term sums over all edges $e\u2208E$ containing node $n$, $Ke$ is the coupling to edge $e$, $\theta \u2192e$ is the vector of phases of oscillators in edge $e$ with $\theta n$ placed in the last component, and $P$ is a permutation of the remaining components. Adding over all permutations $P$ ensures symmetric coupling from all the other nodes in the hyperedge. The integer-valued vector $ce$ determines how the phases are combined inside the sine function and satisfies $ceT1=0$, where $1$ is a vector of ones. In the case of the pairwise all-to-all Kuramoto model, for example, $Ke=K/N$ and $ce=[1,\u22121]T$.

Here, we will study the case where there are hyperedges only of sizes 2 (links) and 3 (triangles), $Ke=K2$, $ce=[1,\u22121]T$ for links, and $Ke=K3$, $ce=[2,\u22121,\u22121]T$ for triangles. With these choices, Eq. (2) can be rewritten as

where we assume that the hypergraph is described by symmetric tensors with entries $Anm$ and $Bnjm$, where $Anm=1$ ($0$) if nodes $n,m$ are connected (not connected) by a link, and $Bnjm=1(Bnjm=0)$ if nodes $n,j,m$ are connected (not connected) by a triangle. However, the techniques that we present can be applied to the general case (2) as long as the hypergraph is generated (or can be approximated) with a generative model like the one discussed in Sec. II A.

Higher-order interactions of the form of Eq. (3) can arise when a phase oscillator model is derived from an expansion beyond the first order of the complex Ginzburg–Landau equation (e.g., see Refs. 7 and 8). The diffusive-type coupling case where $ce=[1,1,\u22122]T$ for triangles has been studied for the all-to-all case in Ref. 12, so here we focus for simplicity on the form of the interactions in Eq. (3).

## III. DIMENSIONALITY REDUCTION

In this section, we use the Ott–Antonsen ansatz^{28} to derive a low-dimensional description of the dynamics and use it to find semi-analytical expressions for the order parameters. In order to accomplish this, we use a generalization of the ansatz in which oscillators are divided into subgroups of oscillators with the same hyperdegree, with oscillators in each subgroup assumed to be statistically equivalent.^{39–41} Furthermore, we neglect pair correlations among oscillators connected in triangles. These approximations are presented and discussed below. Using this procedure, we obtain a low-dimensional description in terms of the functions determining the probabilities of connection between the different subgroups, i.e., the functions $a(m)$. This low-dimensional description allows us to find conditions for synchronization and for the appearance of bistability of the synchronized and incoherent states.

Defining the local order parameters

we can rewrite Eq. (3) as

where we defined $Hn=K2Rn(1)+K3Rn(2)$.

Following Ref. 39, now we assume, based on the construction of the hypergraph from the generative model in Sec. II A, that nodes with the same hyperdegree $k$ are statistically equivalent, and make the identification

Moving to the continuum description in the limit as $N\u2192\u221e$, we define $f(\theta ,\omega ,k,t)$ to be the density of oscillators with phase $\theta $, natural frequency $\omega $, and hyperdegree $k$ at time $t$. Thus, we divide the population of oscillators into subpopulations characterized by their hyperdegree, which acts as a population parameter as in Refs. 39–41. In a mean-field approximation, the global order parameter $R(1)(k)$ can be written in terms of the connection probabilities $a(2)(k,k\u2032)$ and $a(3)(k,k\u2032,k\u2033)$ introduced in Sec. II A as

Similarly, the global order parameter $R(2)(k)$ can be written in terms of the joint density of two oscillators, $f2(\theta ,\omega ,\theta \u2032,\omega \u2032,k,k\u2032,t)$ as

where, to make further progress, we have neglected pair correlations and assumed that the joint density can be written as

We offer the following heuristic arguments to support this assumption: first, in the limits of total incoherence and total synchronization Eq. (9) is exact. Second, when each oscillator is connected to many others, the correlations between any specific pair of oscillators should be small. Thus, we anticipate that this approximation will be a good one either close to total synchrony or incoherence or for dense hypergraphs. This approximation is discussed further in Sec. V. For the regular Kuramoto model the effect of including pair (and higher) correlations has been studied in Ref. 42.

Due to the conservation of oscillators, the evolution of $f$ is governed by the continuity equation

To reduce the dimensionality of this system, we write $f$ as a Fourier series,

where c.c. denotes complex conjugate, and use the Ott–Antonsen ansatz^{28} $bn(\omega ,k,t)=(b(\omega ,k,t))n$. Substituting this ansatz in Eq. (10), one finds that the continuity equation is satisfied if $b(\omega ,k,t)$ satisfies the ODE

Assuming a Lorentzian distribution of frequencies, $g(\omega )=\Delta /(\pi [\Delta 2+(\omega \u2212\omega 0)2])$ and using contour integration to evaluate the integrals in Eqs. (13) and (14), we get

Inserting these in Eq. (12) and letting $\omega =\omega 0\u2212i\Delta $ and $b(k,t)=b(\omega 0\u2212i\Delta ,k,t)$, we get

Equation (17) provides a low-dimensional description of the dynamics in terms of the hypergraph generative functions $a(2)$ and $a(3)$. While the number of variables $b(k)$ might still be large, Eq. (17) allows us to study the bifurcations and fixed points of the system. For this, it is useful to define the global order parameters

which can be written in terms of $b(k,t)$ as

The factor of $2$ in the definition of $R(2)$ accounts for the fact that each triangle is counted twice in the calculation of $Rn(2)$; note that in the case of complete synchronization, $b=1$, the normalization (1) ensures $R(2)=1$.

In the following, we will demonstrate the application of this formalism to selected examples.

## IV. EXAMPLES

In this section, we apply our theory to two examples: a random hypergraph analogous to an Erdös–Rényi network, and a hypergraph where the triangle and link degrees are correlated.

### A. Random hypergraph

We start by considering the hypergraph analog of an Erdös–Rényi network, i.e., a hypergraph where a link connects every pair of nodes with probability $p2$ and a triangle connects every triad of nodes with probability $p3$. Synchronization on this hypergraph was studied numerically in Ref. 10. In terms of the average numbers of links and triangled per node, $\u27e8k\u27e9$ and $\u27e8q\u27e9$, using Eq. (1), we obtain

where we assumed $N\u226b1$. Inserting these in Eq. (17), we find that all $b(k)$ satisfy the same equation,

where

Since all $b(k)$ approach the same attractors, we look for stationary rotating solutions of the form $b(k)=bei\Omega t$, $V1(t)=V1ei\Omega t$, $V2(t)=V2e2i\Omega t$. Each $b(k)$ is assumed to have the same complex phase dictated by the fourth and fifth terms of Eq. (23). After separating real and imaginary parts, we find that $\Omega =\u2212\omega 0$ and that $b$ satisfies $b=0$ or (note that $V1=b$, $V2=b2$)

Solving for $b$ and noting that, from Eq. (19), $|R(1)|=b$, we get

where $K^2=\u27e8k\u27e9K2$ and $K^3=2\u27e8q\u27e9K3$. This generalizes the all-to-all result [Eq. (5) in Ref. 10] to random hypergraphs by properly rescaling the dyadic and triadic coupling strengths. (The factor of $2$ in $K^3$ can be understood in the context of the normalization used in Ref. 10 by noting that in the all-to-all case the mean triangle degree is $\u27e8q\u27e9\u2248N2/2$.) Depending on the values of $K^2$ and $K^3$, Eq. (26) can have zero, one, or two real solutions. As noted in Ref. 10, for $K^3<2$, the system undergoes a supercritical pitchfork bifurcation from incoherence ($R(1)$ = 0) to synchronization ($|R(1)|$ ¿ 0) at $K^2=2$. For $K^3>2$, the system is incoherent for $K^2<22K^3\u2212K^3$. At $K^2=22K^3\u2212K^3$, there is a saddle-node bifurcation where a pair of stable and unstable synchronized solutions appear. At $K^2=2$, the unstable solution disappears in a subcritical pitchfork bifurcation. The phase diagram, mirroring that for the all-to-all case in Ref. 10, is shown in Fig. 1.

### B. Correlated links and triangles

Now, we move to an example where the structure of links is correlated with the structure of triangles. We assume that a prescribed degree sequence is given, ${k1,k2,\u2026,kN}$, and links are created between nodes according to the Chung–Lu model,^{32} so that

Following Ref. 30, we consider a model where the probability that a triangle connects nodes with degrees $k,k\u2032$, and $k\u2033$ is given by

The normalization is chosen using Eq. (1) so that $\u27e8k(2)\u27e9=\u27e8k(3)\u27e9=\u27e8k\u27e9$. By construction, the expected degrees $kn(2)$ and $kn(3)$ for a given node $n$ coincide, so we call this model *correlated links and triangles*.^{30} The model is illustrated in Fig. 2. Since a node is only characterized by a single degree $k$, in the rest of this section, we index all quantities by $k$ only, i.e., we write $bk$ instead of b(**k**). Using the forms for $a(2)$ and $a(3)$ above, Eq. (17) becomes

Defining

Eq. (29) can be rewritten as

where $\Delta $ comes from the Lorentzian distribution of frequencies ($g(\omega )=\Delta /(\pi [\Delta 2+(\omega \u2212\omega 0)2])$). As compared to Eqs. (13) and (14), the ODE on Eq. (32) is much simplified. We have reduced Eqs. (7) and (8) to a closed set of ODEs in terms of variables $bk(t)$ coupled to two global variables $U1$ and $U2$. Now, seeking a stationary rotating solution, we let $bk(t)=bkei\Omega t$, $U1(t)=U1ei\Omega t$ and $U2(t)=U2ei\Omega t$. Then Eq. (32) becomes

The imaginary part of Eq. (33) gives $\Omega =\u2212\omega 0$, and the real part simplifies to

Solving for $bk$, we get

where we chose the solution that satisfies $bk\u21920$ when $K3=0$, $K2\u21920$. Inserting this expression into the definition of $U1$ and $U2$, we find the self-consistent equations

Using Eqs. (27) and (28) in Eqs. (18) and (19), we find that the order parameters $R(1)$ and $R(2)$ can be expressed in terms of $U1$ and $U2$ as

Note that setting $K3=0$ in Eq. (35), one recovers after some manipulation the degree-based mean-field approximation for the network Kuramoto model [i.e., Eq. (25) in Ref. 43 or Eq. (13) in Ref. 44].

With the self-consistent Eqs. (35)–(37), we now proceed to determine the nature of the bifurcation from incoherence ($U1,U2=0$) to synchronization ($U1,U2>0$) with a perturbative approach. Expanding (35) for small $U1$, $U2$ we obtain up to cubic order in $U1$ (note that $U2\u223cU12$)

Letting $U1\u21920+$ to find the onset of synchronization, the leading order terms give the critical coupling strength

Next, solving for $U1$, we find, after canceling the incoherent solution, that close to the transition $U1$ satisfies

where

Thus, a bifurcation occurs at $K2=K2c$, which is independent of $K3$ and equal to the critical constant for the network Kuramoto model in the mean-field approximation.^{43,44} The bifurcation is supercritical for $a>0$ and subcritical for $a<0$. Evaluating Eq. (44) at $K2=K2c$, we find that the transition is subcritical (and, therefore, explosive and with hysteretic behavior) for

For regular networks with $ki=k$, there is bistability for $K3>K2c/2=1/\u27e8k\u27e9$. For networks with a diverging fourth moment [such as networks with a power-law degree distribution with exponent $\gamma \u2208(4,5)$ in the limit $N\u2192\u221e$], $K3c$ diverges and there is no bistability.

Now, we validate our theoretical results with numerical simulations. First, we generate a sequence of $N=5000$ target degrees ${k1,k2,\u2026,kN}$ drawn randomly from a uniform distribution in ${30,31,\u2026,70}$. Then, we create links and triangles connecting nodes according to Eqs. (27) and (28) and generate a synthetic hypergraph. To each oscillator, we assign a frequency drawn from a Lorentzian distribution with $\Delta =1$ and $\omega 0=0$ by setting $\omega n=tan\u2061(\pi (2n\u2212N\u22121)(N+1))$. For this hypergraph, we have $K2c\u22480.038$ and $K3c\u22480.021$.

In Fig. 3, we show the steady-state value of $|R(1)|$ as $K2$ is adiabatically increased and then decreased (blue circles and pink crosses, respectively) for $K3=0.02<K3c$ [Fig. 3(a)], and $K3=0.05>K3c$ [Fig. 3(b)]. For each $K2$, Eq. (3) was solved numerically using Heun’s method with a time step $\Delta t=0.002$ for 100 time units, and the value of $|R(1)|$ was averaged for the last 4 time units. For $K3=0.02<K3c$ [Fig. 3(a)], the transition to synchronization is continuous and there is no hysteresis. On the other hand, for $K3=0.05>K3c$ [Fig. 3(b)], the transition is explosive, and there is a hysteresis loop as $K2$ is increased and then decreased (indicated with arrows). In general, the numerical solution of Eq. (3) agrees well with the numerical solution of the self-consistent Eqs. (36) and (37), shown as black lines, except for $K2\u2248K2c$. The dashed black line corresponds to an unstable solution of Eqs. (36) and (37). The observation that higher-order interactions promote bistability and hysteresis are consistent with findings in Refs. 10 and 12, where higher-order interactions occur via all-to-all simplicial complexes. We note the discrepancy between the order parameters predicted by the mean-field theory and those calculated numerically. This could be a result of the finite network size we used for numerical simulations or our neglect of pair correlations in Eq. (9).

To further illustrate the bistable nature of the system, in Fig. 4, we plot $|R(1)(t)|$ vs $t$ for fixed $K3=0.05$ and $K2=0.005$ (a), $K2=0.03$ (b), and $K2=0.07$ (c), corresponding to the incoherent, bistable, and synchronized regimes, respectively. For each value of $K2$, we use five different initial conditions with $|R(1)(0)|\u22480,0.2,0.4,0.6,0.8$. The solid (dashed) red lines indicate the stable (unstable) solutions of the steady-state self-consistent Eqs. (36) and (37). The values of $|R(1)(t)|$ approach the values predicted by the mean-field theory, including both stable values in the bistable regime [Fig. 4(b)].

## V. DISCUSSION

In this paper, we explored the synchronization of phase oscillators on hypergraphs with heterogeneous structures, generalizing the results in Refs. 10 and 12 to more complex scenarios. The mean-field approximation allowed us to predict the onset of synchronization, explosive transitions between synchronized and incoherent states, and their bistability as a function of system parameters. In the absence of hyperedges of size larger than 2, we recover a smooth transition between incoherent and synchronized states as found in the standard network Kuramoto model.^{25,43} Sufficiently strong higher-order interactions lead to an abrupt transition and bistability of incoherent and synchronized states (see Ref. 45 for a broader perspective of this issue). For a hypergraph with correlated links and triangles, we showed that the onset of synchronization and the onset of bistability depend on the moments of the degree distribution. For the hypergraph model we considered, higher-order interactions only affect the onset of bistability, but not the onset of synchronization (however, see additional discussion on this point below). We have also verified that similar results hold true for networks with power law and bimodal degree distributions.

The main limitations of our study are the requirement for hypergraphs to be produced by the generative model of Sec. II, the use of the mean-field approximation, and the use of approximation (9). [Here, we refer to the approximation that all nodes with the same hyperdegree are statistically equivalent as the mean-field approximation, rather than neglecting pair correlations in Eq. (9)]. The generative model we used assumes that the presence of a hyperedge connecting a group of nodes depends only on a set of pre-determined quantities of these nodes, which might not capture the generative mechanisms behind some real-world or model hypergraphs. For example, a simplicial complex model where triangles only join triples of nodes that are already forming a clique (as assumed in some studies^{10,46}) is not included in the class of models that the generative model in Sec. II covers. In such a model, correlations between the states of nodes belonging to the same triangle could be non-negligible, and, thus, approximation (9) could break down. In that case, the techniques introduced in Ref. 42 could be needed to account for pair correlations. For example, for the SIS model on a simplicial complex, Ref. 47 finds that the epidemic threshold is only predicted correctly when accounting for pair correlations. In addition, Ref. 48 recently noted that synchronization properties in the strongly synchronized regime differ between simplicial complexes and random hypergraphs. Exploring the limitations and possible extensions of our method for simplicial complexes is an interesting problem left for future work.

Despite the limitations discussed above, our framework constitutes a flexible method to study the synchronization of phase oscillators on complex hypergraphs. While we demonstrated our framework in a particular case [hypergraphs constructed following Eqs. (27) and (28)], we emphasize that the techniques presented here allow for the study of a much larger class of systems. Examples include hypergraphs with independently chosen link and triangle degree distributions, correlations between hyperedge degrees and frequencies, and varying degrees of correlations between link and triangle degrees. The techniques presented here open a way to understand the effects that a large class of structural properties of hypergraph connectivity can have on the synchronization of coupled oscillators.

## ACKNOWLEDGMENTS

J.G.R. and S.A. acknowledge useful discussions with Nicholas Landry. J.G.R. and P.S.S. acknowledge a useful discussion with Christian Bick. S.A. was partially supported by NSF (Grant No. DMS-2205967).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Sabina Adhikari:** Data curation (equal); Formal analysis (equal); Investigation (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Juan G. Restrepo:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Per Sebastian Skardal:** Conceptualization (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.