Multiorder Laplacian for synchronization in higher-order networks

The emergence of synchronization in systems of coupled agents is a pivotal phenomenon in physics, biology, computer science, and neuroscience. Traditionally, interaction systems have been described as networks, where links encode information only on the pairwise inﬂuences among the nodes. Yet, in many systems, interactions among the units take place in larger groups. Recent work has shown that the presence of higher-order interactions between oscillators can signiﬁcantly affect the emerging dynamics. However, these early studies have mostly considered interactions up to four oscillators at time, and analytical treatments are limited to the all-to-all setting. Here, we propose a general framework that allows us to effectively study populations of oscillators where higher-order interactions of all possible orders are considered, for any complex topology described by arbitrary hypergraphs, and for general coupling functions. To this end, we introduce a multiorder Laplacian whose spectrum determines the stability of the synchronized solution. Our framework is validated on three structures of interactions of increasing complexity. First, we study a population with all-to-all interactions at all orders, for which we can derive in a full analytical manner the Lyapunov exponents of the system, and for which we investigate the effect of including attractive and repulsive interactions. Second, we apply the multiorder Laplacian framework to synchronization on a synthetic model with heterogeneous higher-order interactions. Finally, we compare the dynamics of coupled oscillators with higher-order and pairwise couplings only, for a real dataset describing the macaque brain connectome, highlighting the importance of faithfully representing the complexity of interactions in real-world systems. Taken together, our multiorder Laplacian allows us to obtain a complete analytical characterization of the stability of synchrony in arbitrary higher-order networks, paving the way toward a general treatment of dynamical processes beyond pairwise interactions.


I. INTRODUCTION
The emergence of order in populations of interacting oscillators-a phenomenon known as synchronization-is ubiquitous in natural and man-made systems [1,2].Typical examples of synchronization include the flashing of fireflies, or the clapping of an audience.In the last decades, synchronization has been the subject of intense research, and it has been applied to a wide range of areas, including neuroscience [3], circadian rhythms [4], or the cardiovascular system [5,6].In particular, much attention has been devoted to unveiling the relationship between the structure of the network of interactions and the emerging collective behavior [7,8].As an outcome of these investigations, noticeable examples include the discovery of abrupt synchronization induced by degreefrequency correlation [9], and cluster synchronization induced by structural symmetries [10].
Most interacting systems have so far been represented as networks, a collection of nodes and links describing relationships and influences between them at the level of pairs.However, many real-world systems are better modeled by including higher-order interactions, i.e., interactions between more than two nodes at a time [11].A typical example is that of human collaborations, which often occur at the level of groups.In this case, a traditional network representation is misleading, as it would associate the same structure-a triangle-both with the case of a triplet of people collaborating on a single task, and with the case of three individuals collaborating as three distinct pairs in different projects.This indistinguishability can be solved by making use of higherorder network representations, such as hypergraphs [12] and simplicial complexes [13].A stream of research has recently focused on correctly characterizing the structure of systems with higher-order interactions [14][15][16][17][18][19][20].Interestingly, considering this additional level of complexity sometimes leads to changes in the emerging dynamics of complex systems, including social contagions [21,22], activity-driven models [23], diffusion [24,25], random walks [26,27], and evolutionary games [28].
Higher-order interactions can influence the nature of the dynamics also for systems of coupled oscillators.A few studies have recently investigated their effect experimentally [29,30], and from a network inference point of view [31][32][33].From a theoretical point of view, higher-order interactions were considered in the context of global nonlinear coupling [31] and multipopulation resonance [34], they were shown to arise from phase reduction beyond the first approximation [30,[35][36][37][38], and they can facilitate chaotic behavior [35] and other exotic dynamical regimes [30].The Kuramoto model, where phase oscillators interact in pairs, is often invoked as the most simple way to describe the emergence of synchronization in a population of interacting nodes with local dynamics.Only a few works have so far considered higher-order generalizations [39][40][41][42], showing the promotion of cluster synchronization [39][40][41] and explosive transitions [41,42].Yet, for all these studies, analytical insights are limited to all-to-all coupling settings, disregarding the rich architecture of interactions of real-world systems.Interestingly, a different type of model was introduced in [43], where a Kuramoto phase oscillator is associated with each simplex.
In this work, we provide a full analytical treatment for synchronization in a population of coupled oscillators where arbitrary higher-order interactions of all orders are possible.
To this end, we study a generalization of the Kuramoto model of identical phase oscillators to general group interactions.These can be conveniently described by hypergraphs, in the most flexible mathematical representations of higher-order interactions, which also generalize other commonly used formalisms, such as simplicial complexes.For this model, we show that the stability of the fully synchronized state is determined by the eigenvalues of a newly defined multiorder Laplacian, which takes into account the higher-order complex topology of interactions.We validate this Laplacian framework on several toy models describing higher-order interactions of increasing complexity.First, we investigate all-to-all interactions at all orders, for which the eigenvalues can be derived fully analytically.We further characterize this system by considering three subcases: (i) attractive coupling only, (ii) interplay between attractive and repulsive orders, (iii) and decaying coupling strength (that we link to higher-order phase reduction studies).We confirm our analytical findings with numerical simulations.Second, we consider the starclique model, a toy model specifically generated to highlight some simple spectral properties of the multiorder Laplacian.Finally, we investigate the effect of higher-order interactions on synchronization on a real macaque brain dataset.Taken together, our work sheds new light on the effect of higherorder interactions in a population of coupled oscillators, and it unveils how new emergent phenomena can be captured analytically through the introduction of a suitable Laplacian framework, which naturally generalizes the traditional approach to networks beyond pairwise interactions.

II. GENERALIZED HIGHER-ORDER INTERACTIONS
We study the effect of higher-order interactions in a population of N identical phase oscillators.Specifically, we consider the dynamics of oscillators with the most general topology describing many-body interactions of any order d = 1, . . ., D, (1) a natural generalization of the Kuramoto model, where ω is the natural frequency of each oscillator, γ 1 , γ 2 , . . ., γ D are the coupling strengths at each order, the K (1) , . . ., K (D)  are the average degrees at order 1, . . ., D [explicitly defined in Eq. ( 4)], and the adjacency tensors M determine the topology.Just like A i j = 1 if there is a pairwise interaction (i, j) but 0 otherwise, B i jk = 1 if there is a triplet interaction (i, j, k) but 0 otherwise, and similarly for all orders.Note that the interactions are assumed undirected, i.e., the adjacency tensors are invariant under any permutation of their indices.These adjacency tensors encode the most general topology of higher-order interactions that can be formalized as hypergraphs or simplicial complexes, for example.The largest value that D can take is N − 1, which corresponds to N-oscillator interactions, the highest order possible.The general interaction scheme of (1) is illustrated with an example in Fig. 1.There, a three-oscillator interaction is represented by a 2-simplex, and any (d + 1)-oscillator interaction is represented by a d-simplex (also called a simplex of order d), as illustrated in Fig. 1(b).Figure 1(c) shows a visualization of each pure order d.
The sinusoidal coupling functions are chosen as a natural generalization of those used in the Kuramoto model: when describing the dynamics of oscillator i, they are symmetric with respect to i, meaning that any permutation of the other indices leaves them invariant.Even when restricting ourselves to 2πperiodic functions that vanish when oscillators are identical, there exist more choices at larger orders.At order 2, for example, the only other choice is sin(2θ j − θ k − θ i ), used, e.g., in [42].All these other choices are discussed in Appendix B and do not affect the generality of our framework.For the sake of clarity of our presentation, however, we do not include them earlier.Sinusoidal couplings naturally arise when applying phase reduction to more realistic nonlinear models, for the Kuramoto model at first order, but also at higher orders [35,38], e.g., for nanoelectromechanical oscillators [30].
A few studies have investigated systems similar to (1) analytically, but only in all-to-all schemes [39][40][41][42]44].Insights into complex topologies, for instance the onset of synchronization [42], have so far been limited to numerical simulations, and a formal description of the process is still an open problem.In Sec.III A, we introduce an analytical framework that allows us to overcome current limitations and investigate arbitrarily complex topologies.Such general 1. Population of oscillators with higher-order interactions on simplicial complexes.(a) Example of simplicial complex of oscillators: the red node has higher-order interactions of orders up to 3 (4-oscillator) with the other oscillators.(b) Here, the building blocks of the higher-order interactions consist in edges (1-simplices), triangles (2-simplices), and tetrahedra (3-simplices).A d-simplex represents a (d + 1)oscillator interaction.(c) The simplicial complex (a) can be decomposed into its pure d-simplex interactions: the red node has seven 1-simplex interactions (blue), four 2-simplex interactions (orange), and one 3-simplex interaction (green), but no 4-or more simplex interactions.In (a), the faces of the tetrahedra belong to 2-simplices (orange) and to a 3-simplex (green) so that they are depicted in gray.
patterns are formalized as hypergraphs, which are the mathematical structures that allow for the most general encoding of higher-order interactions.These include-but are not limited to-simplicial complexes, more constrained higher-order representations where a group interaction of order d requires all its lower order subgroup interactions, illustrated in Fig. 1 for model (1).
For the introduced model (1), the existence of the fully synchronized state, θ i = θ j for all i and j, is trivially guaranteed and implies the solution θ i (t ) = ωt.In this paper, we focus on the stability of that fully synchronized state.For convenience, we start by going to the rotating reference frame ψ i = θ i − ωt.This is equivalent to applying the transformation θ i → ψ i and ω → 0 to the original system (1).In this new reference frame, the synchronized solution is given by ψ i (t ) = 0 for all i = 1, . . ., N.
The linear stability of the synchronized state is determined by the dynamics of heterogeneous perturbations δψ i (t ), which satisfy the linearized dynamics ( These equations are straightforwardly obtained by linearizing system (1).It is worth noticing, however, that the same equations could also be derived by assuming a different (not necessarily sinusoidal) shape of the coupling functions in the original system, via the linearization mechanism.For networks with pairwise interactions only, the dynamics of those perturbations-and hence the stability of the system-is typically assessed by using the so-called Laplacian formalism.In the next section, we see how we can characterize the stability of system (1) with higher-order interactions up to any order d by extending the traditional Laplacian formalism to systems with any type of higher-order interactions.

III. MULTIORDER LAPLACIAN
Different generalizations of the Laplacian operator have been proposed so far in the literature to include higher orders of interactions: from the simplest versions for uniform hypergraphs [45,46], to those more complicated associated with simplicial complexes [47][48][49] and Hodge Laplacians [25,50], to mention a few.Let us notice that these Laplacians describe the hierarchy among building blocks of the topology, and different orders are associated with Laplacian matrices of different sizes, where the order zero is the traditional node point of view; the first order represents the edge perspective where the Laplacian size is equal to the number of pairwise connections, and each entry is associated with edge adjacency; the second-order Laplacian has a different size again, being based on the existing triangles, and so on.
Here, we propose a multiorder Laplacian: an operator that generalizes the typical pairwise Laplacian framework and allows us to analytically describe the effect of higher-order couplings on node oscillatory dynamics for any simplicial interactions.This is different from the previously defined Laplacians, where interactions between, e.g., triangles, are seen from the point of view of triangles and not from the point of view of the nodes in those triangles.
First, we show that each d-simplex interaction term in Eq. ( 2) can be written in terms of a generalized Laplacian of order d.Second, we display how the full system (2) can be written in terms of a multiorder Laplacian.

A. Laplacian
We introduce a generalized Laplacian of order d, with the Kronecker delta δ i j , and where we have defined at order d the degree K (d ) i , i.e., the number of distinct d-simplices that node i is part of, and the adjacency matrix A (d )  i j , i.e., the number of distinct d-simplices that the pair of nodes (i, j) is part of, ( Note that these definitions are natural generalizations of their pairwise counterparts to which they reduce when d = 1.This newly defined Laplacian can be shown to have the expected properties of a standard Laplacian matrix: it is symmetric, and its rows sum to zero.Moreover, its eigenvalues are all nonnegative, as we shall see in the next section.With those quantities, each term of the linearized equation ( 2) can be rewritten as as shown in detail in Appendix A. We now have all the ingredients to treat oscillators at each order of interaction and build a multiorder Laplacian.

Multiorder interactions
We go back to our original system (1), with interactions at orders d = 1, . . ., D combined.We know that the stability of the synchronized solution of system (1) is determined by system (2), which we now can write where we have defined the multiorder Laplacian L (mul)   i j as a weighted sum of the Laplacian matrices of order d.The weight given to each order is proportional to γ d , and normalized by the average degree of order d.Hence, by definition, the multiorder Laplacian gives an equal weight to each order, even if the network contains more, say, 2-simplices than 5-simplices.Notice that this normalization is not included in the definition of Laplacians of pure order d.This newly defined Laplacian reduces to the usual operator when D = 1, i.e., when only pairwise interactions are taken into account.Finally, we note that L (mul)   i j depends on D, whose maximum value is limited by the network size, since when D = N − 1, all possible orders are considered.
In this section, our generalized framework showed us two things.First, how to rewrite interactions at each order with a Laplacian matrix L (d )  i j of order d.And second, how to rewrite the full system, including all higher-order interactions, with a multiorder Laplacian matrix L (mul) .Additionally we showed how the latter matrix is just a weighted sum of the former matrices.
So far, we have an analytical expression for the Laplacian matrix of a given simplicial complex.It is the eigenvalues of this Laplacian that quantify the stability or instability of the synchronized state of the system.We note that, even though the multiorder Laplacian is the weighted sum of each Laplacian of order d, its eigenvalues cannot be obtained in general as a linear combination of the Laplacians at each order d, as the eigenvalue operator is nonlinear.In general, we need to numerically compute the eigenvalues of the Laplacian.While this is generally true, special cases exist where the eigenvalues can nonetheless be summed and the system characterized in a fully analytical manner.We present such a case in Sec.IV.
Before this, let us make a short didactic digression and step back to analyze the pure order 2. We show the details of the Laplacian derivation in this specific case and leave those for order 3 and d in Appendix A.

Pure 2-simplex interactions
We now rewrite the three-oscillator interaction term in Eq. ( 2), i.e., those interactions represented by filled orange triangles in Fig. 1, with our Laplacian framework.The first simple mathematical step represents an important point of the derivation, allowing us to write the 2-simplex interaction of three phases as two identical terms of the difference of only two phases: by using the symmetry B i jk = B ik j and using (δψ The phase difference in expression ( 9) is similar to the pairwise case, and this allows us to naturally generalize the Laplacian formalism to 2-simplex interactions.Indeed, the 2-degree K (2)   i of node i, i.e., the number of distinct 2-simplices that node i is part of, is , where the factor 2! ensures each 2-simplex is counted only once.For example, in Fig. 1(a), the red node has a degree of order 2, K (2)   i = 4.The adjacency matrix of order 2, whose entries A (2)  i j represent the number of 2-simplices shared by the pair (i, j), is A (2)  i j = N k=1 B i jk , which is a natural generalization of the usual pairwise adjacency matrix A. Indeed, just as A i j = 1 if i and j are part of a common 1-simplex interaction but 0 otherwise, A (2)  i j = n if i and j are part of n common (but distinct) 2-simplex interactions, but 0 otherwise.
With these definitions in hand, we can now perform the second important step of our procedure, which will take us straight to the Laplacian.We rewrite the 2-simplex interaction term (9) as follows: where we obtained the last line by defining the Laplacian of order 2 as in Eq. ( 3), as a natural generalization of the usual pairwise Laplacian.With Eq. ( 10), we have explicitly shown the link between the structure of the 2-simplex interactions and the dynamics of the oscillators by casting it into a Laplacian form.Notice that the two essential analytical steps that allowed us to write the Laplacian of order 2 can be straightforwardly generalized at each order.In Appendix A, we show the same procedure for order 3 and generic order d.
It is important to mention another relevant feature of our formalism.In Eq. ( 1), we considered the interactions at orders higher than 1 occurring by means of coupling functions designed to be natural generalizations of that at the first order.Specifically, we chose sinusoidal functions that are symmetric with respect to index i and vanish at synchronization.We note here that this two-step derivation above also allows us to treat any other such choice of coupling function, as shown in Appendix B. Indeed, by performing the same first step, one notices that terms that do not contain the phase i vanish, and the terms left yield a fraction c 0 /d of the Laplacian of order d, where c 0 is the integer coefficient that multiplies θ i in the coupling function.Physically, this means that these other choices, which arise naturally in phase reduction studies beyond the first approximation (see, e.g., Refs.[37,38]), show a slower convergence to synchronization.

B. Stability and Lyapunov exponents
We are finally able to study the stability of our system of oscillators.The synchronized state is stable if the perturbation δψ i on each node i converges to zero.
Let us first consider pure d-simplex interactions.To this end, we need to solve Eq. ( 6) to obtain the temporal evolution of the perturbation.To do so, we can make use of the Laplacian eigenvalues (d )  α and eigenvectors φ (d ) α defined by this eigenbasis can be used to project the perturbation vec- α , where the c α are real constants.By plugging this solution into system (6), we can decouple our system of N equations and obtain the N Lyapunov exponents of the synchronized state The Lyapunov exponents are a measure a stability: the system with pure d-simplex interactions is stable if all values λ (d ) α are negative, so that the perturbations δψ i tend to zero over time.By convention, the Lyapunov exponents are ordered: N .The Laplacian eigenvalues are non-negative by definition so that if all γ d > 0 (attractive coupling), the synchronized state is always stable up to a global phase shift: 0 = λ (d )  1 > λ (d ) 2 .In the multiorder system (7), the stability of the synchronized equilibrium is determined by the interplay of all different orders, as encoded in the multiorder Laplacian.We can then analogously use the spectrum and eigenbasis of L (mul) , hence obtaining the N Lyapunov exponents The Laplacian eigenvalues are non-negative by definition so that if all γ d > 0 (attractive coupling), the synchronized state is always stable.The Lyapunov exponent determines the longterm behavior of the second Lyapunov exponent, λ (mul) 2 , i.e., the smallest nonzero one.Its value determines the resilience of the system to perturbations, i.e., how fast the system comes back to the stable state after a perturbation.In particular, the more negative the λ (mul) 2 is, the more stable is the synchronized state.

IV. STABILITY IN ALL-TO-ALL HIGHER-ORDER NETWORKS
In this section, we analyze a simple case-dubbed higherorder all-to-all-which can be solved analytically.This setting is a generalization of the usual "all-to-all" (or global) coupling scheme in traditional networks.Indeed, in networks with pairwise interactions, all-to-all coupling indicates that every possible pairwise interaction takes place.Similarly, in a network with higher-order interactions, higher-order all-toall indicates that every possible (d + 1)-oscillator interaction occurs, for all orders d.This setting is the only one that has been studied analytically, with a main focus on cluster states.In this homogeneous case, each term of order d in the original system (1) can be written in terms of the order parameter amplitude R 1 and its phase 1 as, up to a normalization factor, which makes apparent the driving by the mean field that is now nonlinear.Effectively, each oscillator is driven by the same mean field with strength R d 1 and with a dth harmonic.Pure harmonics are known to yield stable cluster states, which can be checked via a self-consistency argument.If we take d = 2, then Eq. ( 13) has two stable fixed points with distance π , and that two-cluster state has R 1 > 0, which keeps driving the system.These considerations were used in [41] for the pure triplet case to study cluster states and abrupt transition in the thermodynamic limit.Still in the pure triplet case, Ref. [40] shows that although the incoherent state is stable in the thermodynamic limit, finite-size effects can destabilize it and yield cluster states.In [42], the authors, combining interactions of orders up to 4, unveiled the emergence of an abrupt transition to synchronization in the thermodynamic limit.Finally, in [44], the authors extended the Watanabe-Strogatz lowdimensional description to any pure higher harmonics l for the general system θi = ω + Im[He −ilθ i ], where H is a function of the generalized order parameters, which can take the form (13) and yield higher-order interactions of any order.With their framework, the authors tracked the basins of attraction of the clusters in cluster states.While this treatment does not focus primarily on higher-order interactions, it can be related via the nonlinear mean-field coupling.
Here, in contrast with these studies, we focus on the stability of full synchrony and provide a full analytical description of its spectrum.In addition, we investigate the effect considering a variable number of orders, mixing attractive and repulsive couplings, and decaying coupling strengths as observed in phase reduction studies.In particular, we show that, due to the absence of complex topology, the Lyapunov spectrum of the full system reduces to a linear combination of those at each pure order.We remark that this is not true in complex topologies, for which the aggregated spectrum is determined by a nonlinear combination of each order.

A. Higher-order Laplacians are proportional to the traditional pairwise Laplacian
We now show that at each order d, the higher-order allto-all Laplacian L (d ) is proportional to the usual pairwise all-to-all Laplacian L (1) , defined by L (1)  i j = K (1)  i j δ i j − A (1)  i j .Consequently, since the analytical spectrum of the latter is known, the analytical spectrum of L (d ) can also be obtained, and in turn, that of L (mul) .
First, we need to explicitly write down our degree and adjacency matrices defined in general in Eqs. ( 4) and (5).In the all-to-all case, the degree reduces to ), in combinatorics notation.This, for order 1, yields the usual K (1)   i = N − 1.The adjacency matrix reduces to A (d )  i j = ( N−2 d−1 )(1 − δ i j ), which yields the usual A (1)  i j = 1 − δ i j at order 1.Additionally, the following identities will prove useful in the next section to relate the quantities at order d to their traditional counterparts (at order 1), Now that we have these expressions, we have all we need to write an explicit formula for the multiorder Laplacian (3) in the next section.In the higher-order all-to-all setting, all nodes have the same degree of order d.Hence, in this section, we write K (d )   i and K (d ) interchangeably.More details are provided in Appendix C.

Pure d-simplex interactions
For the general case of order d, analogously injecting expressions ( 14) and (15) into definition (3) yields which can be rewritten Hence, we have shown that the Laplacian of order d is proportional to the usual pairwise Laplacian for any order d.In addition, Eq. ( 17) indicates that L (d ) i j is linearly growing with the order of the interactions d, and with K (d ) , which is the number of d-simplex interactions that each oscillator has.
We will discuss the implications that these dependencies have on the stability in more detail in Sec.IV B. For now, we only stress that the analytical formula (17) serves as a limit case to understand the behavior of L (d )  i j in more complex coupling schemes than the higher-order all-to-all scheme.Examples of those will be investigated in Sec.V, in which a full analytical derivation is not always possible.

Multiorder interactions
Oscillators in the higher-order network have multioscillator interactions with the other oscillators at all orders d = 1, . . ., D. In general, the stability of the synchronized solution is determined by the eigenvalues of the multiorder Laplacian L (mul) , in Eq. ( 8), as we have shown in Sec.III A. In the higherorder all-to-all setting, by injecting Eq. ( 17), this Laplacian reduces to Hence, the stability of the synchronized solution is only determined by the usual pairwise Laplacian, as well as by the strength γ d of each interaction of order d, and by the order d of those interactions.We give the full analytical spectrum of L (mul) for the higher-order all-to-all case in the next section.We will see that, as a consequence of the proportionality between L (d )  i j and L (1)  i j , the multiorder eigenvalue spectrum is just a linear combination of the spectra at each order.However, this is only true in the all-to-all scheme, which overlooks the complexity of real-world topologies.

B. Spectrum and stability
In this section, we give an analytical formula for the eigenvalues of the higher-order all-to-all Laplacian matrix, and consequently for the Lyapunov exponents that determine the stability of the synchronized solution.We remind the reader that the attractiveness or repulsiveness of interaction at a given order d is determined by the sign of the coupling strength γ d : positive and negative coupling strengths correspond to attractive and repulsive interactions, respectively.We consider different scenarios, where we tune at will the sign and intensity of the coupling strengths γ d .

Attracting couplings at all orders
We have shown in the previous section that the d-order Laplacians are intrinsically connected to the usual pairwise Laplacian L (1) .Therefore, this latter shapes the Laplacian spectrum at each order d, and consequently also at the multiorder.For higher-order all-to-all networks, the spectrum of L (1) is degenerate and given by from which we derive the Lyapunov exponents of each order d: The second Lyapunov exponent is reported in Fig. 2(c) for different values of d as a function of network size N.It appears clear that interactions of higher orders d stabilize the synchronized state more (more negative λ (d ) 2 as d increases).Interestingly, this is true despite each order being given an equal weight through the normalization in system (1).This is illustrated with numerical simulations for the pure orders d = 1 and 2, respectively, shown in Figs.2(a) and 2(b).Indeed, trajectories converge faster with pure 2-simplex interactions than with pure pairwise ones.
Then, combining all orders, we obtain the multiorder Lyapunov exponents From this formula, we see that the more orders are taken into account, i.e., the more D is increased, the more negative λ (mul) 2 is, as shown in Fig. 2(d).Physically, additional attractive higher-order interactions tend to stabilize the synchronized state, as expected.

Attractive even orders and repulsive odd orders
Here, we apply our analytical framework to investigate the interplay between attractive and repulsive interactions.For simplicity, we study the case where attractive and repulsive relationships are associated with different interaction orders.Mixing attractive and repulsive couplings has been previously investigated, motivated by, e.g., analogies with inhibitory and excitatory connections between neurons [51].
First, we restrict ourselves to only 1-simplex (2-oscillator) and repulsive 2-simplex (3-oscillator) interactions.Depending on the sign of γ 1 , 1-simplex interactions are attractive (γ 1 > 0) or repulsive (γ 1 < 0).The attractiveness or repulsiveness of 2-simplex interactions depends identically on the sign of γ 2 .Physically, attractive interactions favor synchronization by stabilizing the synchronized state.By contrast, repulsive interactions will favor incoherence by destabilizing the synchronized state.The result of this interplay between the interactions at various orders depends on the respective coupling strengths γ d : is shown in Fig. 3(c) for a range of positive and negative values of γ 1 and γ 2 .We see that synchronization can be stable even if the traditional pairwise (1-simplex) interactions are repulsive, as long as 3-oscillator (2-simplex) interactions are attractive enough to counterbalance them, as illustrated in Fig. 3(a).This highlights a potential benefit of considering higher-order interactions: indeed they might stabilize the systems, in cases when the purely pairwise system is unstable.This confirms a similar result obtained in [41] for phase oscillators with distributed frequencies.In addition, attractive pairwise interactions might be outplayed by repulsive three-oscillators ones.This can result in an unstable synchronized state, as illustrated in Fig. 3(b).
Second, we generalize to all orders up to D with alternating signs.Specifically, we consider all interactions of an even and odd order d to be attractive and repulsive, respectively, by setting 1 , with γ 1 = 0.6.The Lyapunov exponent converges to a negative value as higher-order interactions are taken into account.Now, the Lyapunov exponent is given by the alternating series which diverges as D → ∞ (which requires N → ∞).For increasing but finite values of D, however, the Lyapunov exponent alternates between positive and negative values, as illustrated in Fig. 3(d).In other words, if the highest order of interaction considered is odd, the synchronized state is stable.However, if it is even, the synchronized state is unstable.This is due to the fact that the contribution to λ (mul) 2 of each order d is proportional to d: adding one repulsive or attractive order outplays all lower-order interactions.Finally, we note that the factor N/(N − 1) → 1 in the limit of large networks N → ∞, so that it can be seen as a finite-size correction factor.

Weaker higher orders: Link to phase reduction
Here, we make a brief link between our formalism and the higher-order phase reduction approaches developed, for example, in [33,35,36,38].In these phase reduction studies, the authors obtain a phase model from an original network of nonlinear oscillators by performing a sophisticated perturbative expansion in a small parameter.This small parameter is usually linked to the original pairwise coupling strength γ 1 .The authors find that, at higher orders in the expansion, i.e., at higher powers of γ 1 , there appear terms including higher-order interactions, i.e., interactions between more than two phases.
Motivated by these studies, we consider a scenario of decaying coupling strengths.Specifically, we set 0 γ 1 < 1, and couplings at higher orders as powers of the pairwise coupling strength which means that interactions between many oscillators are weak.In this case, the stability is given by the series which is always negative, and hence the stability of the synchronized state is ensured.In addition, this is a geometric series that converges to −γ 1 /(γ 1 − 1) 2 when D (and hence N) tends to infinity.This convergence is illustrated in Fig. 3(e).
Even though the qualitative result, i.e., stability, was expected since interactions at all orders are attractive, such a quantitative result has more predictive power.In particular, we note that −γ 1 /(γ 1 − 1) 2 → −∞ as we go away from the domain of validity of the perturbative regime γ 1 → 1.

V. STABILITY IN HIGHER-ORDER NETWORKS WITH ARBITRARY TOPOLOGY
In this section, we apply our multiorder Laplacian framework to simplicial complexes with complex heterogeneous structures.We consider two cases: first, a toy model-the starclique model-and second, a real network, i.e., a macaque brain dataset.

A. The star-clique model
As a first example of complex topology with higher-order interactions, we consider a toy model that we call star-clique [27].We make use of this model in order to observe in a simple case the multiorder Laplacian properties and thus the stability conditions for synchronization.The star-clique, as its name indicates, is composed of two subnetworks: a star of size N s and a clique of size N c , which we treat analogously to a higher-order all-to-all network of size N c from the previous section.The subnetworks are connected only by one link from the center of the star to a single node of the clique, as illustrated in Fig. 4(a).By construction of this simplicial complex, the two subnetworks are almost disconnected, and they work almost as two independent systems.We will see, however, that at each pure order d 2 the subnetworks are actually disconnected.This allows us to use our analytical understanding from the previous section to comprehend the present case with complex topology.
We start by analyzing the structure of the 1-simplex interactions.The Lyapunov spectrum of the pure 1-simplex interactions is shown in blue in Fig. 4(b).The spectrum reflects the quasi-independence of the two subnetworks.Indeed, since the star and the clique are almost independent, the adjacency matrix A i j is a slightly perturbed matrix with two independent blocks, one for each subnetwork.The same holds for the Laplacian matrix, and hence its eigenvalues and the Lyapunov exponents.We know that the N s Laplacian eigenvalues (1)  α of a star network with pairwise interactions are given by {0, 1, . . ., 1, N s + 1}.We also know that the N c Laplacian eigenvalues (1)  α of a pairwise all-to-all network are given by {0, N c , . . ., N c }. Hence, the spectrum of N eigenvalues (1)   α of the 1-simplex Laplacian L (1) is the union of both spectra, only slightly perturbed.This can be seen in Fig. 4(b) for the Lyapunov exponents, i.e., scaled Laplacian eigenvalues (blue symbols).
Moving to higher-order interactions, d 2, the star and the clique actually are disconnected.In fact, in these pure dsimplex cases, there is no star, since it has no higher-order interactions.Hence the higher-order Laplacians only "see" the clique.This implies that the spectrum of L (2)  i j is the union of N s zeros and the spectrum of a higher-order all-to-all that we derived analytically in the previous section: one element zero and the others with value −γ 2 2N c /(N c − 1).Analogously at orders d higher than 2, the spectrum is given by N s + 1 zeros and N c − 1 elements −γ d dN c /(N c − 1).See For what concerns the spectrum of the total Laplacian, it is not exactly the sum of the spectrum of order d, as in the all-to-all case, but almost.Indeed, the only thing that distinguishes it from the previous case is that the pairwise Laplacian of the star-clique network is not composed of two disconnected blocks, but almost.Let us observe that if the pure orders d 2 can be simply described by the all-to-all higher order spectrum, the order d = 1 (traditional pairwise) and the multiorder reflect the structure of the two almost disconnected subgraphs, namely the star and the clique.Their spectrum is compared in Fig. 4(c).It is important to notice that λ (1)  2 and λ (mul) 2 are identical.The two spectra are, however, globally different, and in particular λ (1)  N and λ (mul) N have very different values, suggesting that in the repulsive case (all γ d < 0) the two Laplacians would yield dynamics with very different timescales.

B. Real network: Macaque brain
In this section, we demonstrate the use of our multiorder Laplacian framework on a real dataset.We use a macaque brain dataset publicly available at [52], consisting of 91 nodes and 628 edges.The nodes represent cortical areas, and the edges represent links among those, experimentally determined by retrograde tracing [53].We assume that each fully connected (d + 1)-node clique captures multiorder interactions that can be described as a d-simplex, and we study the dynamics of (1) on top of the resulting simplicial complex.The obtained simplicial complex as it appears at orders 1, 4, and 9 is depicted in Figs.5(a)-5(c), respectively.A whole spectrum of models of coupled oscillators has been used to study the synchronization of neurons in the brain.Here we do not intend this to be a realistic model of a functional brain, but rather we use this dataset as an empirical basis to test our theoretical toy model of synchronization in systems with higher-order interactions.While there are clear limitations in applying extremely simplified neuronal models to capture real brain dynamics, idealized phase models are still interesting, as they can sometimes capture well more realistic dynamics such as integrate-and-fire ones [54].
The full distribution of the number of d-simplices present in the simplicial complex is shown in Fig. 5(a).The highest order of interaction is 10.As expected, the distribution follows a bell curve: the simplicial complex contains fewer than 1000 1simplices, i.e., 2-oscillator interactions, and 10-simplices, i.e., 11-oscillator interactions, but it has over 6000 5-simplices.In addition to the distribution of d-simplices across orders d, the dynamics is affected by the structure of the d-simplex interactions at each order d.Specifically, we wonder how 1-, 2-, ..., and d-simplices are distributed among nodes.It turns out that at each order d, the interactions very centralized: only few have many d-simplex interactions, whereas the majority of the nodes have very few or no d-simplex interaction at all.This is shown for orders 4 and 9 in Figs.5(b) and 5(c), respectively.The higher the order d is, the more nodes have zero d-simplex interactions.
How is the stability of the synchronized state affected by this inter-and intraorder structure?To assess this, we compute the Laplacian (3) at each order d as well as the multiorder Laplacian (8), and we obtain the Lyapunov exponents from the eigenvalue spectra.First, we assess the stability of each pure case, at orders d from 1 to 4, by computing their respective Lyapunov spectra shown in Fig. 5(e).Orders higher than 4 are not shown for visualization purposes.These spectra reflect the structure we described above: the higher the order d, the more nodes that have no d-simplex interactions, i.e., that are disconnected from the main component in the pure d-simplex network.These disconnected nodes translate into eigenvalues of value 0 in the spectra.We observe also two other stages of the spectra: In the first step, higher orders correspond to lower Lyapunov exponents, which is similar to what happens for the star-clique model because of the clique subnetwork and may reflect, also in this example, the aggregated nature of the network at each order d.Then, we observe a strong descending behavior for the last eigenvalues, which confirms the stronger stability of the higher orders.
Finally, we compare the stability of the full system with combined higher-order interactions against what is obtained on the corresponding projecting graph, where nodes coupled at any order are linked with a pairwise edge.This is illustrated in Fig. 5(f).

VI. SUMMARY AND CONCLUSIONS
In summary, in this work we have studied the effects of higher-order interactions, i.e., multioscillator interactions, on the synchronization of identical phase oscillators on hypergraphs.We considered interactions up to any order in an extension of the Kuramoto model by introducing a multiorder Laplacian, which makes the system amenable to analytical treatment, and determines the stability of the fully synchronized state.In particular, we obtained a quantitative measure of that stability by computing the Lyapunov exponents of the state, which are proportional to the Laplacian spectrum.We showed applications of the multiorder Laplacian in settings of increasing complexity.We emphasize that such analytical treatment is especially important for higher-order systems, as numerical simulations become slower with each additional order of interactions, quickly becoming unfeasible.
As a first example, we considered a case that is fully tractable analytically: the higher-order all-to-all setting.In this setting, all nodes interact in all possible d-simplex interactions for any order d.In this case, we showed that the Laplacian of each order d, L (d ) , is proportional to the traditional pairwise Laplacian L (1) , and to the value of d itself; see Eq. (17).As a consequence, the stability of the multiorder system can be linearly decomposed into the stability of each of the pure d-simplex systems.Indeed, in this special setting the multiorder Lyapunov exponents are merely a weighted sum of the Lyapunov exponents of order d.We confirm our findings with numerical simulations.Finally, this fully tractable higher-order all-to-all setting also serves as a limit case to help us understand what happens in more complex topologies.
In this context, we investigated analytically two additional scenarios.First, we investigated the interplay between orders with attractive interactions (γ d > 0) and orders with repulsive interactions (γ d < 0).For example, we showed that their opposite effect on synchronization means that repulsive pairwise couplings can be countered effectively by higher-order attractive couplings.In general, we showed that when odd orders are attractive and even orders are repulsive, taking more higher orders into account can stabilize or destabilize the system, depending on the highest order considered.Second, to link our work to higher-order phase reduction techniques, we considered the attractive coupling strength γ d , which decays as the order d increases.We derived the convergence of the Lyapunov exponents as higher orders are considered.
Then, we considered two cases with more complex topologies.First, we applied the multiorder Laplacian framework to a toy model: the star-clique network.With this model, we illustrated spectral properties of the multiorder Laplacian, and of the Laplacian at each order.In that specific topology, the star subnetwork does not contain higher-order interactions and hence does not influence the value of the first nonzero Lyapunov exponent.Higher-order interactions in the clique subnetwork, however, do change the other values of the spectrum, as we showed.Second, we considered a real-world topology: a macaque brain network.In this setting, we show how the complex shape of the multiorder Laplacian spectrum can be understood from the structure of the simplicial complex.More importantly, our analysis confirms changes in all Lyapunov exponents due to the inclusion of higherorder interactions, as compared to only pairwise interactions.Specifically, higher-order interactions tend to stabilize synchronization, with more drastic change on the more negative part of the spectrum, due to the hublike structure of the dataset.We showed in Appendix B that, for higher-order complex topologies, the choice of the coupling function affects the timescales, causing them to reach full synchrony.Yet, the precise effect of different coupling functions [55,56] on other dynamical regimes on such systems is an open question.
In conclusion, in this paper we introduced a multiorder Laplacian to assess the stability of synchronization in populations of oscillators with higher-order interactions.Our framework has two main strengths: (i) it is a natural generalization of the traditional and well-known pairwise Laplacian framework, and (ii) it can be applied to arbitrary hypergraphs describing any structured group interactions up to any order d.This is in contrast with previous studies, where analytical insights were provided for the higher-order all-to-all coupling scheme only.Besides, other than its inherent ability to deal with complex topologies, the Laplacian formalism is valid for any N, allowing one to investigate finite-size effects.
Our framework promises to find wide applicability.Indeed, it has very recently been used to extend the ideas behind the master stability function [57] to simplicial complexes [58].Taken together, the multiorder Laplacian is a powerful tool that could be used widely not only to characterize populations of oscillators with higher-order interactions, but also for a wider general analytical treatment of dynamical processes beyond pairwise interactions.
We briefly show how to rewrite the 3-simplex interactions of Eq. ( 2) in a similar way to the 2-simplex case from Sec. III A. Although it is similar, the case of 3-simplices better prepares us for the next step: the general d-simplex case.The first important step is to reduce the expression with four phases into an expression of only two phases: by using the invariance of C i jkl under index permutations, and where the factor 2! = 3!/3 comes from the intermediate step where the expression is written as the sum of three identical terms.At order 3, definitions (3)-( 5) for the degree K (3)  i , the adjacency matrix A (3)  i j , and the Laplacian L (3)  i j read C i jkl , (A2) C i jkl , (A3) L (3)  i j = 3K (3)  i δ i j − A (3)  i j . (A4) We remind the reader that the degree of order 3, K (3)  i , of node i is the number of distinct 3-simplex interactions it is part of, and A (3)  i j is the number of shared 3-simplices including nodes i and j.With these definitions, we can now perform the second important step and rewrite (A1) as 2!A (3)  i j δψ j − δψ i 3!K (3) A (3)  i j − 3K (3)  i δ i j δψ j = − γ 3 L (3)  i j δψ j (A5) in a similar way to the usual case of order 1.This shows that the dynamics of the 3-simplex interactions, i.e., of four oscillators, can be rewritten in terms of a Laplacian of order 3. Indeed, once again such an operator fulfills the requirements to be a Laplacian, being symmetric, positive-semidefinite, and zero row-sum.

Order d
The rewrite of order d follows the same two important steps as the derivation above.This shows that, at any order d, the dynamics caused by interactions of order d (d + 1 oscillators) is determined by the matrix L (d ) .Such a matrix fulfills the properties necessary to be a Laplacian so that its eigenvalues are non-negative and include at least one zero.
the number of ways to pick d − 1 oscillators among the N − 1 oscillators left.This number is given by the combinatorics formula A (d ) i j = ( N−2 d−1 )(1 − δ i j ), where (1 − δ i j ) ensures that entries with i = j are equal to zero.
Second, we proceed similarly to write explicitly the degree of order d.In a usual all-to-all setting with only pairwise interactions, every oscillator has a pairwise interaction with all N − 1 oscillators left, i.e., K (1)   i = N − 1.The degree of order d is equal to the number of (d + 1)-oscillator interactions of which oscillator i is a part.Hence, it is given by the number of ways to pick d oscillators out of the N − 1 left, which can be written K (d ) i = ( N−1 d ) in combinatorics notation.Note that it scales with the size of the system as K (d ) FIG.1.Population of oscillators with higher-order interactions on simplicial complexes.(a) Example of simplicial complex of oscillators: the red node has higher-order interactions of orders up to 3 (4-oscillator) with the other oscillators.(b) Here, the building blocks of the higher-order interactions consist in edges (1-simplices), triangles (2-simplices), and tetrahedra (3-simplices).A d-simplex represents a (d + 1)oscillator interaction.(c) The simplicial complex (a) can be decomposed into its pure d-simplex interactions: the red node has seven 1-simplex interactions (blue), four 2-simplex interactions (orange), and one 3-simplex interaction (green), but no 4-or more simplex interactions.In (a), the faces of the tetrahedra belong to 2-simplices (orange) and to a 3-simplex (green) so that they are depicted in gray.

FIG. 2 .
FIG. 2. Higher-order all-to-all: attractive coupling.Higher-order interactions increase the stability of synchronization.Synchronization of N = 100 oscillators over time with (a) only 1-simplex (2-oscillator) interactions and (b) only 2-simplex (3-oscillator) interactions, obtained by numerical integration.The oscillators with only 2-simplex interactions synchronize faster than those with only 1-simplex ones.This is confirmed by (c): the analytical second Lyapunov exponent λ (d ) 2 of Eq. (21) as a function of network size N, at each order d = 1, . . ., 4. The Lyapunov exponent is more negative for larger orders d of interaction.(d) Lyapunov exponent of the full system, as a function of D, the largest order taken into account.The higher the order of interactions taken into account, the more negative the Lyapunov exponent is, i.e., the most stable the synchronized state is.Parameters are N = 100 and γ d = 1 for all d.

FIG. 3 .
FIG. 3. Higher-order all-to-all: (a)-(d) interplay of attractive and repulsive coupling orders, and (e) decaying coupling strength.Phases over time with (a) which synchronize even with repulsive 1-simplex interactions due to attractive 2-simplex interactions [γ 1 = −0.5, γ 2 + −0.5, black plus in (c)], but do not synchronize with attractive 1-simplex interactions because of repulsive 2-simplex interactions [γ 1 = +0.5, γ 2 = −0.5, black star in (c)].Time series are numerically integrated on a simplicial complex of N = 100 nodes.(c) Analytical, nonzero, multiorder Lyapunov exponent (with D = 2) for a range of 1-and 2-simplex coupling strengths (γ 1 , γ 2 ).Positive and negative coupling strengths correspond to attractive and repulsive coupling, respectively.Negative values (blue) and positive values (red) of λ (mul) 2 indicate instability and stability of the synchronized state, and confirm numerics in (a) and (b).The Lyapunov exponent vanishes when γ 1 + 2γ 2 = 0 (dashed black).(d) Analytical, nonzero, multiorder Lyapunov exponent as a function of the largest order taken into account, D, for attractive and repulsive coupling at even and odd orders, respectively.The synchronized state changes its stability as higher orders are taken into account.(c) Analytical, nonzero, multiorder Lyapunov exponent as a function of D for decaying coupling strengths γ d = γ d1 , with γ 1 = 0.6.The Lyapunov exponent converges to a negative value as higher-order interactions are taken into account.

FIG. 4 .
FIG.4.Synchronization in the star-clique simplicial complex (a) with N s = 6 and N c = 7.The highest-order of interaction is D = 6.In the clique, faces are part of simplices of five different orders, and hence they are depicted in gray.(b) Lyapunov spectra for each case of pure d-simplex interactions.(c) Lyapunov spectrum with higher-order interactions combined (black), γ d = 1 for all d, compared to the traditional pure 1-simplex case (blue).

Fig. 4 (
b) for the spectrum at different orders d.

FIG. 5 .
FIG. 5. Real macaque brain dataset: structure and dynamics.Panels (a), (b), and (c) depict the simplicial complex as it appears at orders d = 1, 4, and 9, respectively, and they report the associated distributions of degree K (d ) i across nodes.At each order, most nodes have zero or few interactions, but a few nodes have many.(d) Number of d-simplices in the simplicial complex for each d.The highest order of interaction in the network is D = 10.(e) Lyapunov spectra with only 1-simplex interactions, and with higher-order interactions (with d up to D = 4).(f) Comparison of the Lyapunov spectra: multiorder (black) and projected graph (blue), i.e., D = 1.The spectrum is affected by higher-order interactions.