Maxime Lucas
email: maxime.lucas.1@univ-amu.fr

Giulia Cencetti

Federico Battiston
email: battistonf@ceu.edu

A multi-order 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 influences 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 significantly affect the emerging dynamics. However, these early studies have mostly considered interactions up to 4 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 higherorder interactions of all possible orders are considered, for any complex topology described by arbitrary hypergraphs, and for general coupling functions. To this scope, we introduce a multi-order 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 multi-order 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 multi-order Laplacian allows us to obtain a complete analytical characterization of the stability of synchrony in arbitrary higher-order networks, paving the way towards a general treatment of dynamical processes beyond pairwise interactions.

I. INTRODUCTION

The emergence of order in populations of interacting oscillators -a phenomenon known as synchronizationis ubiquitous in natural and man-made systems [START_REF] Pikovsky | Synchronization: A Universal Concept in Nonlinear Sciences[END_REF][START_REF] Strogatz | Sync: The emerging science of spontaneous order[END_REF]. 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 has been applied to a wide range of areas, including neuroscience [START_REF] Cumin | [END_REF], circadian rhythms [4], or the cardio-vascular system [5,6]. In particular, much attention has been devoted to unveil the relationship between the structure of the network of interactions and the emerging collective behavior [7,8]. As outcome of these investigations, noticeable examples include the discovery of abrupt synchronization induced by degree-frequency 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 relationship and influences between them at the level of pairs. However, many real-world systems are better mod-eled by including higher-order interactions, i.e. interactions between more than two nodes at a time. A typical example is that of human collaborations, which often occur at the level of groups. In this cases, a traditional network representation is misleading, as it would associate the same structure -a triangle -both to the case of a triplet of people collaborating on a single task, and to the case of three individuals collaborating as three distinct pairs in different projects. This indistinguishability can be solved by making use of higher-order network representations, such as hypergraphs [START_REF] Berge | Hypergraphs: combinatorics of finite sets[END_REF] and simplicial complexes [START_REF] Aleksandrov | [END_REF]. A stream of research has recently focused on correctly characterizing the structure of systems with higher-order interactions [13][14][15][16][17][18][19]. Interestingly, considering this additional level of complexity sometimes leads to changes in the emerging dynamics of a complex system, including social contagion [20,21], activity driven models [22], diffusion [23,24], random walk [START_REF] Zhou | Advances in neural information processing systems[END_REF][START_REF] Carletti | [END_REF] and evolutionary games [27].

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 [28,29], and from a network inference point of view [30][31][32]. From a theoretical point of view, higher-order interactions were considered in the context of global nonlinear coupling [30] and multi-population resonance [33], were shown to arise from phase reduc-tion beyond the first approximation [29,[34][35][36][37], and can facilitate chaotic behavior [34] and other exotic dynamical regimes [29]. 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 [38][39][40][41], showing the promotion of cluster synchronization [38][39][40] and explosive transitions [40,41]. 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 [42], where a Kuramoto phase oscillator is associated to 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 scope, 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 generalizes 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 multi-order 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 star-clique model, a toy model specifically generated to highlight some simple spectral properties of the multi-order 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 higher-order interactions in a population of coupled oscillators, and 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 θi = ω + γ 1 K (1) N j=1 A ij sin(θ j -θ i ) + γ 2 2! K (2) N j,k=1 B ijk sin(θ j + θ k -2θ i ) + γ 3 3! K (3) N j,k,l=1 C ijkl sin(θ j + θ k + θ l -3θ i ) + ... + γ D D! K (D) N j1,...,j D =1 M ij1...j D sin D m=1 θ jm -D θ i , (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, and the adjacency tensors M determine the topology. Just like A ij = 1 if there is a pairwise interaction (i, j) but 0 otherwise, B ijk = 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 which 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 3-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 The sinusoidal coupling functions are chosen as 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 the invariant. Even when restricting ourselves to 2π-periodic functions that vanish when oscillators are identical, there exists more choices at larger orders. At order 2, for example, the only other choice is sin(2θ j -θ k -θ i ), used e.g. in [41]. All these other choices are discussed in App. B and do not affect the generality of our framework. For the sake of the clarity of our presentation, however, we do not include them earlier.

A few studies have investigated systems similar to (1) analytically, but only in all-to-all schemes [38][39][40][41]43]. Insights on complex topologies, for instance on the onset of synchronization [41], were so far limited to numerical simulations, and a formal description of the process is still an open problem. In Sec. III, we introduce an analytical framework which allows us to overcome current limitations and investigate arbitrarily complex topologies. Such general patterns are formalized as hy- Here, the building blocks of the higher-order interactions consist in edges (1-simplices), triangles (2-simplices) and tetrahedra (3simplices). 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.

Simplicial complex of oscillators

pergraphs, 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 [START_REF] Pikovsky | Synchronization: A Universal Concept in Nonlinear Sciences[END_REF].

For the introduced model (1), 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 [START_REF] Pikovsky | Synchronization: A Universal Concept in Nonlinear Sciences[END_REF]. 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

δ ψi = + γ 1 K (1) N j=1 A ij (δψ j -δψ i ) + γ 2 2! K (2) N j,k=1 B ijk (δψ j + δψ k -2δψ i ) + γ 3 3! K (3) N j,k,l=1 C ijkl (δψ j + δψ k + δψ l -3δψ i ) + ... + γ D D! K (D) N j1,...,j D =1 M ij1...j D D m=1 δψ jm -D δψ i .
(2) 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 higherorder interactions up to any order d by extending the traditional Laplacian formalism to systems with any type of higher-order interactions.

III. MULTI-ORDER 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 [START_REF] Chung | Expanding graphs[END_REF][START_REF] Lu | International Workshop on Algorithms and Models for the Web-Graph[END_REF], to those more complicated associated to simplicial complexes [START_REF] Muhammad | Proc. of 17th International Symposium on Mathematical Theory of Net-works and Systems[END_REF][START_REF][END_REF][START_REF] Parzanchevski | [END_REF] and Hodge Laplacians [24,[START_REF] Lim | Proc. Sympos. Appl. Math[END_REF], to mention a few. Let us notice that these Laplacians describe the hierarchy among building blocks of the topology, and different orders are associated to 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 to edge adjacency; the second order Laplacian has a different size again, being based on the existing triangles, and so on.

Here, we propose a multi-order Laplacian: a new operator which 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 multi-order Laplacian.

A. Laplacian

We introduce a generalized Laplacian of order d

L (d) ij = dK (d) i δ ij -A (d) ij . (3)
with the Kronecker delta δ ij , and where we have defined at order d, the degree K (d) i , i.e. the number of distinct d-simplices node i is part of, and the adjacency matrix A (d) ij , i.e. the number of distinct d-simplices the pair of nodes (i, j) is part of,

K (d) i = 1 d! N j1,...,j D =1 M ij1...j D , (4)
A (d) ij = 1 (d -1)! N j2,...,j D =1 M ij1...j D . (5)
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 non-negative, as we shall see in the next section.

With those quantities, each term of the linearized equation (2) can be rewritten as

δ ψi = - γ d K (d) N j=1 L (d) ij δψ j . (6)
as shown in detail in App. A. We now have all the ingredients to treat oscillators at each order of interaction and build a multi-order Laplacian.

Multi-order interactions. We go back to our original system [START_REF] Pikovsky | Synchronization: A Universal Concept in Nonlinear Sciences[END_REF], 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

δ ψi = - N j=1 L (mul) ij δψ j , (7)
where we have defined

L (mul) ij = D d=1 γ d K (d) L (d) ij , (8)
the multi-order Laplacian L

(mul) ij
as a weighted sum of the Laplacian matrices of order d. The weight given to each order is proportional to γ d , and normalised 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 5simplices. 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) ij 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) ij of order d. And second, how to rewrite the full system, including all higher-order interactions, with a multi-order 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 multi-order 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 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 the Appendix A.

Pure 2-simplex interactions. We now rewrite the 3-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 3 phases as 2 identical terms of the difference of only 2 phases:

δ ψi = γ 2 2! K (2) N j,k=1 B ijk (δψ j + δψ k -2δψ i ), (9) = γ 2 K (2) N j,k=1 B ijk (δψ j -δψ i ), (10)
by using the symmetry B ijk = B ikj and using (δψ

j + δψ k -2δψ i ) = (δψ j -δψ i ) + (δψ k -δψ i ).
The phase difference in expression ( 10) 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 num- ber of distinct 2-simplices node i is part of, is K (2) i = 1 2! N j,k=1
B ijk where the factor 2! ensure 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) ij represent the number of 2-simplices shared by the pair (i, j), is A

(2) ij = N k=1 B ijk
which is a natural generalization of the usual pairwise adjacency matrix A. Indeed, just as A ij = 1 if i and j are part of a common 1-simplex interaction but 0 otherwise, A

(2) ij = 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 (10) as follows

δ ψi = γ 2 K (2) N j,k=1 B ijk (δψ j -δψ i ) (11) = γ 2 K (2)   N j=1 A (2) ij δψ j -δψ i 2!K (2) i   , (12)
= γ 2 K (2) N j=1 A (2) ij -2K (2) i δ ij δψ j , (13)
= - γ 2 K (2) N j=1 L (2) ij δψ j , (14)
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. ( 14), 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 which 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 showed in App. 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. [36,37], 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 scope, 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 Λ

ij (φ (d) α ) j = Λ (d) α (φ (d)
α ) i , with α = 1, . . . , N . Indeed, this eigenbasis can be used to project the perturbation vector

δψ i (t) = N α=1 c α exp(λ (d) α t) φ (d)
α , 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

λ (d) α = - γ d K (d) Λ (d) α . (15)
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: λ

(d) 1 ≥ λ (d) 2 ≥ • • • ≥ λ (d)
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 multi-order system (7) the stability of the synchronized equilibrium is determined by the interplay of all different orders, as encoded in the multi-order Laplacian. We can then analogously use the spectrum and eigenbasis of L (mul) , hence obtaining the N Lyapunov exponents

λ (mul) α = -Λ (mul) α . (16)
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 that determines the long-term behaviour the second Lyapunov exponent, λ

, i.e. the smallest non-zero 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
, the more stable the synchronized state.

IV. STABILITY IN ALL-TO-ALL HIGHER-ORDER NETWORKS

In this section we analyze a simple case -dubbed higher-order 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-to-all 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

θi = γ d R d 1 sin[d(Φ 1 -θ i )] = γ d Im[Z d 1 e -iθi ] (17)
which makes apparent the driving by the meanfield that is now nonlinear. Effectively, each oscillator is driven by the same meanfield 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. Take d = 2, then Eq. ( 17) has two stable fixed points with distance π, and that 2-cluster state has R 1 > 0 which keeps driving the system.

These considerations were used in [40] for the pure triplet case, to study cluster states and abrupt transition in the thermodynamic limit. Still in the pure triplet case, [39] shows that although the incoherent state is stable in the thermodynamic limit, finite-size effects can destabilize it and yield cluster states. In [41], the authors, combining interactions of orders up to 4, unveil the emergence of an abrupt transition to synchronization in the thermodynamic limit. Finally, in [43], the authors extended the Watanabe-Strogatz low-dimensional 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 ( 17) 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 primarily focus on higher-order interactions, it can be related via the nonlinear meanfield 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) ij = K (1) ij δ ij -A (1) ij .
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)- (5). In the all-to-all case, the degree reduces to K

(d) i = N -1 d
, in combinatorics notation. This, for order 1, yields the usual K

(1) i = N -1. The adjacency matrix reduces to

A (d) ij = N -2 d-1 (1 -δ ij )
, which yields the usual A

(1) ij = 1δ ij 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)

A (d) ij = [(N -2) • • • (N -d)/(d -1)!] A (1) ij , (18)
K (d) i = [(N -2) • • • (N -d)/d!] K (1)
i .

Now that we have these expressions, we have everything in hand to write an explicit formula for the multi-order 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 App. C.

Pure d-simplex interactions. For the general case of order d, analogously injecting expressions (18) and (19) into definition (3) yields of Eq. ( 27) 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.

L (d) ij = dK (d) i δ ij -A (d) ij , (20)
= [(N -2) • • • (N -d)/(d -1)!] [K (1)
i δ ij -A (1) ij ], (21)
= [(N -2) • • • (N -d)/(d -1)!] L (1) ij , (22)
which can be rewritten

L (d) ij = d K (d) N -1 L (1) ij . (23)
Hence, we have shown that the Laplacian of order d is proportional to the usual pairwise Laplacian, for any order d. In addition, Eq. ( 23) indicates that L (d) ij is linearly growing with the order of the interactions d, and with K (d) , the number of d-simplex interactions 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 (23) serves as a limit case to understand the behaviour of L (d) ij 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.

Multi-order interactions. Oscillators in the higherorder network have multi-oscillator 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 multi-order Laplacian L (mul) , in Eq. ( 8), as we have shown in Sec. III. In the higher-order all-to-all setting, by injecting Eq. ( 23), this Laplacian reduces to

L (mul) ij = D d=1 γ d d N -1 L (1) ij . (24)
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) ij and L

(1) ij , 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 which determine the stability of the synchronized solution. We remind 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 Λ indicate instability and stability of the synchronized state, and confirms numerics in (a) and (b). The Lyapunov exponent vanishes when γ1 + 2γ2 = 0 (dashed black). (d) Analytical, non-zero, multi-order 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, non-zero, multi-order Lyapunov exponent as a function of D, for decaying coupling strengths γ d = γ d 1 , with γ1 = 0.6. The Lyapunov exponent converges to a negative value as higher-order interactions are taken into account.

order d: λ (d) 1 = 0 λ (d) 2,...,N = -γ d d N N -1 . ( 26
)
The second Lyapunov exponent is reported in Fig. 2(c) for different values of d as a function network size N . It appears clear that interactions of higher orders d stabilise 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 normalisation in system (1). This is illustrated with numerical simulations for the pure orders d = 1 and d = 2, respectively shown in Figs. 2(a) and (b). Indeed, trajectories converge faster with pure 2-simplex interactions than with pure pairwise ones.

Then, combining all orders, we obtain the multi-order Lyapunov exponents

λ (mul) 1 = 0 λ (mul) 2,...,N = - N N -1 D d=1 γ d d. ( 27
)
From this formula, we see that the more orders are taken into account, i.e. 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 relationship are associated to different interaction orders. Mixing attractive and repulsive couplings has been previously investigated, motivated by, e.g., analogies with inhibitory and excitatory connections between neurons [START_REF] Hong | [END_REF].

First, we restrict ourselves to only 1-simplex (2oscillator) and repulsive 2-simplex (3-oscillator). 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 : λ

(→2) 2 = -N N -1 (γ 1 +2γ 2 ) which is negative if γ 1 +2γ 2 > 0. The value of λ (→2) 2
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 3oscillator (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 [40] for phases oscillators with distributed frequencies. In addition, attractive pairwise interactions might be outplayed by repulsive 3-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 (a) an even and odd order d to be attractive and repulsive, respectively, by setting

γ 2n = -1 γ 2n+1 = +1. (28)
Now, the Lyapunov exponent is given by the alternating series

λ (mul) 2 = - N N -1 D d=1 d (-1) d+1 , (29)
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 reductions approaches developed for example in [32,34,35,37]. In these phase reductions 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 2 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

γ d = γ d 1 , (30)
which means that interactions between many oscillators are weak. In this case, the stability is given by the series

λ (mul) 2 = - N N -1 D d=1 d γ d 1 < 0. ( 31
)
which is always negative, and hence stability of the synchronized state is ensured. In addition, this is a geometric series which converges to -γ 1 /(γ 1 -1) 2 when D (and hence N ) tend 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 multi-order Laplacian framework to simplicial complexes with complex heterogeneous structures. We consider two cases: first, a toy model -the star-clique model-and second, a real network, i.e. a macaque brain dataset.

A. The star-clique model

As first example of complex topology with higher-order interactions, we consider a toy model that we call starclique [START_REF] Carletti | [END_REF]. We make use of this model in order to observe in a simple case the multi-order 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 work almost as two independent systems. We will see that, 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 1simplex 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 ij 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 Λ [START_REF] Pikovsky | Synchronization: A Universal Concept in Nonlinear Sciences[END_REF] α 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 on Fig. 4(b) for the Lyapunov exponents, i.e. scaled Lapla-cian eigenvalues (blue symbols).

Moving to higher-order interactions, d ≥ 2, the star and the clique actually are disconnected. In fact, in these pure d-simplex cases, there is no star, since it has no higher-order interactions. Hence the higher-order Laplacians only "see" the clique. This implies that spectrum of L

(2) ij 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 multi-order, reflect the structure of the two almost disconnected subgraphs, the star and the clique. Their spectrum is compared in Fig. 4(c). It is important to notice that λ 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 [START_REF] Rossi | AAAI (2015) we used the "macaque-rhesus-brain-2" dataset[END_REF], consisting of 91 nodes and 628 edges. We assume that each fully connected (d + 1)-node clique captures multi-order 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 Fig. 5, respectively (a), (b) and (c). We remark that a whole spectrum of models of coupled oscillators has been used to study the synchronization of neurons in the brain. 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 like integrate-andfire ones [START_REF] Politi | [END_REF].

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 only less than 1000 1-simplices, i.e. 2-oscillator interactions, and 10-simplices, i.e. 11-oscillator interactions, but 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 are very centralized: only few nodes 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 order 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 intra-order structure? To assess this, we compute the Laplacian (3) at each order d as well as the multi-order Laplacian (8), and 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 highers 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 of 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: a first step where higher orders cor-respond 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. multi-oscillator 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 multi-order 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, that are proportional to the Laplacian spectrum. We showed applications of the multi-order 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 turning 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. (23). As a consequence, the stability of the multi-order system can be linearly decomposed into the stability of each of the pure d-simplex systems. Indeed, in this special setting the multi-order 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, depend-ing on the highest order considered. Second, to link our work to higher-order phase reduction techniques, we considered attractive coupling strength γ d that decay 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 multi-order 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 do not influence the value of the first non-zero Lyapunov exponent. They do drastically change the other values of the spectrum, however, 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 multi-order 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 higher-order interactions, as compared to only pairwise. Specifically, higher-order interactions tend to stabilize synchronization, with more drastic change on the more negative part of the spectrum, due to the hub-like structure of the dataset. We showed in App. B that, for higher-order complex topologies, the choice of the coupling function affects the timescales to reach full synchrony. Yet, the precise effect of different coupling functions [53,54] 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 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 [55] to simplicial complexes [56]. Taken together, the multi-order Laplacian is a powerful tool that could be used widely not only to characterize populations of oscillators with higherorder interactions, but also for a wider general analytical treatment of dynamical processes beyond pairwise interactions.

  (b). Fig. (1)(c) shows a visualization of each pure order d.

FIG. 1 .

 1 FIG.1. Population of oscillators with higher-order interactions on simplicial complexes. (a) Example of simplicial complex: 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(3simplices). 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 .

 2 FIG. 2.Higher-order all-to-all: attractive coupling. Higher-order interactions increase the stability of synchronization. Synchronisation of N = 100 oscillators over time with (a) only 1-simplex (2-oscillator) interactions and (b) only 2-simplex (3oscillator) 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

2 FIG. 3 .

 23 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)). Times series are numerically integrated on a simplicial complex of N = 100 nodes. (c) Analytical, non-zero, multi-order 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

FIG. 4 .

 4 FIG.4. Synchronisation in the star-clique simplicial complex (a) with Ns = 6 and Nc = 7. The highest-order of interaction is D = 6. In the clique, faces are part of simplices of 5 different orders, and hence are depicted in grey. (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 traditional pure 1-simplex case (blue).

FIG. 5 .

 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 reporting the associtated 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: multi-order (black), and projected graph (blue), i.e. D = 1. The spectrum is affected by higher-order interactions.

Fig. 4 (

 4 b) for the spectrum at different orders d.

  identical. The two spectra are however globally different and in particular λ (1) N and λ (mul) N are very have very different values, suggesting that in the repulsive case (all γ d < 0) the two Laplacian would yields dynamics with very different timescales. B. Real network: macaque brain

Acknowledgments

The authors would like to thank Duccio Fanelli, Iacopo Iacopini, and Arkady Pikovsky for fruitful discussions, as well as Alain Barrat for feedback on the manuscript. M.L. acknowledges partial support from the "European Cooperation in Science & Technology" (COST): Action CA15109. F.B. acknowledges partial support from the ERC Synergy Grant 810115 (DYNASNET).

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 3simplices better prepares us for the next step: the general d-simplex case. The first important step is to reduce the expression with 4 phases into an expression of only 2 phases

C ijkl (δψ j -δψ i ), (A2) by using the invariance of C ijkl under index permutations, and where the factor 2! = 3!/3 comes from the intermediate step where the expression is written as the sum of 3 identical terms.

At order 3, definitions (3)-( 5) for the degree K

i , the adjacency matrix A

(3) ij , and the Laplacian L

We remind that the degree of order 3, K

i , of node i is the number of distinct 3-simplex interactions it is part of, and A

(3) ij 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 (A2) as

in a similar way to the usual case of order 1. This shows that the dynamics of the 3-simplex interactions, i.e. of 4 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 semi-definite and zero row-sum.

Order d

The rewrite of order d follows the same two important steps as the derivation above. First, we rewrite the term of d + 1 phases as d terms of 2 phases, and see that they are all equal, yielding

With definitions (3)-( 5), we can now perform the second step, which is to rewrite the this difference of two phases in the terms of the Laplacian of order d:

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.

Appendix B: Alternative higher-order coupling functions

As we mentioned in the main text, even restricting ourselves to sine coupling functions that vanish at synchronization, other choices are possible. Here, we show that our framework can readily be used for these other choices. We refer to the coupling functions in system (1), of which there is only one at each order.

Order 2

At order 2, the only other choice is the asymmetric function sin(2θ j -θ k -θ i ). As in the previous appendix, the first step is to separate the 3-phase term into three 2-phase terms:

The second term is as in the symmetric case of the main text, and hence yields 1 2 L

(2) ij . The first term is different, however, but we can use the symmetry between j and k in this term to show that it vanishes:

because B ijk = B ikj . Hence, the full term of order 2 with this function is half that with the symmetric function (see Eq. ( 14))

Remarkably, this means that, for pure triplets, this choice of coupling function will lead to a convergence that is 2 times slower that the symmetric one.

Order d

The same two steps can be used to treat any arbitrary d, by noticing that, in the first step, each term that includes the phase i yields 1/d Laplacian of order d, and all other terms vanish as shown for order 2 above. Hence, for oscillators coupled via a general function

with integer coefficient c j such that sin(0) = 0, the linearized dynamics is determined solely by c 0 as

Appendix C: Higher-order degree and adjacency matrices for all-to-all coupling

In this section, we explicitly write the higher-order adjacency matrix and connectivity defined in section IV A, for each order d, in the higher-order all-to-all setting.

First, we start with the adjacency matrices of order d. In a traditional all-to-all setting, the pairwise adjacency matrix has all entries equal to one, but 0 on the diagonal, which can be written A

(1) ij = 1 -δ ij . At order d, the generic entry (i, j) of the adjacency matrix A (d) ij is equal to the number of distinct d + 1-oscillator interactions including both oscillators i and j, as described above. Hence, the matrix is simply given by the number of ways to pick d -1 oscillators among the N -1 oscillators left. This number is given by the combinatorics formula A

, where (1 -δ ij ) 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 oscillators 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 that oscillator i is part of. Hence, it is given by the number of ways to pick d oscillators out of the N -1 left, which can be written K